鄒 赫， 王企鯤,2*， 薛壯壯

(1.上海理工大學(xué) 能源與動力工程學(xué)院， 上海 200093；2.上海理工大學(xué) 上海市動力工程多相流動與傳熱重點實驗室， 上海 200093)

依據(jù)非牛頓流體是否具有彈性，可將其劃分為純黏性流體與黏彈性流體[1-2]。黏彈性流體特征：同時具有黏性和彈性2種性質(zhì)，如血液[3]、唾液[4]和PVP聚合物溶液[5]等。大量研究表明黏彈性流體具有許多不同于牛頓流體特殊物理現(xiàn)象[6-10]，如：湍流減阻效應(yīng)、剪切稀變效應(yīng)、擠出脹大效應(yīng)、Weissenberg效應(yīng)和Kaye效應(yīng)，在工程實際中具有廣泛的應(yīng)用前景[11-12]。

目前國內(nèi)外學(xué)者對黏彈性流體的研究主要集中于理論與實驗2方面。然而，描述黏彈性流體的數(shù)學(xué)模型十分復(fù)雜，理論研究僅限于簡單流動，結(jié)果不具有普適性；實驗研究結(jié)果直觀，但很難獲得流場中的重要內(nèi)特性參數(shù)，如流場變形率、彈性應(yīng)力等；因而限制了對黏彈性流體流動機理的深入研究。

隨著計算機技術(shù)的迅速發(fā)展，數(shù)值模擬被廣泛應(yīng)用于黏彈性流體流動機理的研究。相較于牛頓流體，黏彈性流體本構(gòu)方程因其形式復(fù)雜且非線性強，因此，對其開展數(shù)值研究的難度遠高于牛頓流體。上述問題，引起眾多學(xué)者的關(guān)注，并相繼展開研究：Alves等[13]和Hulsen等[14]分別應(yīng)用有限體積法和有限元法對黏彈性流體圓柱繞流問題進行數(shù)值研究，指出此類方法無法在局部區(qū)域獲得收斂，表明需采用更為先進的數(shù)值計算方法或更高的網(wǎng)格精細度；Xiong等[15]通過直接數(shù)值模擬(Direct Numerical Simulation，DNS)計算固定翼型(NACA0012)在黏彈性流體中的流動狀態(tài)，指出由于黏彈性流體本構(gòu)方程的雙曲特性會導(dǎo)致計算量和計算所需時間增大。

上述方法均需要編寫專門的程序代碼，針對性較強，且計算復(fù)雜性及難度相對較高，嚴重限制數(shù)值模擬在彈性流體研究領(lǐng)域中的推廣應(yīng)用。目前FLUENT自帶的二次開發(fā)模塊UDF對黏彈性流動進行數(shù)值模擬的研究對編程方法原理的論述較少、缺乏可靠性驗證，在數(shù)值計算穩(wěn)定性與收斂性方面的缺陷并沒有提出相應(yīng)的有效解決方案[16]。

因此，課題組從原理性角度出發(fā)，提出新型黏彈性流體數(shù)值模擬方法。將FLUENT核心求解器中牛頓流體動量方程缺少的彈性應(yīng)力項以源項形式編寫進UDF程序并載入FLUENT，實現(xiàn)對黏彈性流動的數(shù)值模擬。為驗證該方法的可靠性，基于Oldroyd-B本構(gòu)模型的驅(qū)動方腔流與FENE-P本構(gòu)模型的圓柱繞流2個實例作為驗證，并提出將流場計算中非定常問題定常化。

1 黏彈性流體控制方程

黏彈性流體其基本方程主要由連續(xù)性方程、動量方程及本構(gòu)方程組成：

連續(xù)方程：

?ui/?xi=0。

(1)

動量方程：

(2)

對于均勻的高聚物溶液所形成的黏彈性流體，其本構(gòu)方程可統(tǒng)一寫為[17]：

(3)

(4)

式中：μp為聚合物溶液中溶質(zhì)的動力黏度；λ為黏彈性流體的松弛時間；Cij稱為黏彈性流體變形率張量；f(r)為Perterlin函數(shù)；α為已知參數(shù)。

根據(jù)f(r)和α的不同取值，形成黏彈性流體本構(gòu)方程[18]，如表1所示。

表1 黏彈性流體模型分類Table 1 Classification of viscoelastic fluid models

2 基于UDF二次開發(fā)黏彈性流體數(shù)值模擬方法

2.1 UDF 代碼編譯

為便于描述，課題組以基于Oldroyd-B本構(gòu)模型的黏彈性流體為例，來論述UDF程序的具體實施過程。該方法具有普適性，可被推廣至具有其它本構(gòu)關(guān)系的黏彈性流體。

根據(jù)式(1)～(4)，令其中α=0，f(r)=1，并在式(4)等號右側(cè)加入能起到穩(wěn)定計算作用的人工黏性項，便可得到基于Oldroyd-B本構(gòu)模型的黏彈性流體運動控制方程：

(5)

(6)

對于動量方程式(5)，在FLUENT軟件中所計算的牛頓流體動量方程的基本形式為：

(7)

式中Fi為方程的源項。

對于牛頓流體而言，利用源項可額外定義流體的各種體積力(如重力、慣性力等)。將式(7)與式(5)對比，令

(8)

如此便可將黏彈性流體中彈性應(yīng)力項以源項形式添加入FLUENT軟件核心求解器所使用的動量方程之中，這是隨后UDF代碼的編寫的基礎(chǔ)。對于黏彈性流體本構(gòu)方程式(6)，需要將其轉(zhuǎn)化為FLUENT軟件中標準的用戶自定義標量輸運方程形式后，才能為軟件自動求解。

在FLUENT軟件的UDF編程系統(tǒng)中，式(9)為用戶自定義標量輸運方程的一般形式，該式從左邊到右邊分別為非定常項(瞬態(tài)項)、對流項、耗散項和源項：

(9)

式中φ代表用戶自定義標量(或通用變量)。

為迎合式(7)的基本形式，將式(6)2邊同乘以流體密度ρ，即得：

(10)

將式(9)與Oldroyd-B本構(gòu)模型式(10)的各項對比，令：

φ=Cij；

(11)

г=ρκ；

(12)

(13)

這樣便將黏彈性流體Oldroyd-B模型本構(gòu)方程轉(zhuǎn)化為FLUENT軟件標準用戶自定義標量輸運方程的形式，以便于隨后的UDF代碼的編制。

對于式(11)中求解變量Cij，可直接用FLUENT軟件中自帶的自定義標量數(shù)組函數(shù)C_UDSI(c,t,0:2)來定義。以二維問題為例，Cij具有3個獨立變量：Cxx,Cyy,Cxy，可依次對應(yīng)于標量數(shù)組函數(shù)C_UDSI以下分量：C_UDSI(c,t,0)，C_UDSI(c,t,1)，C_UDSI(c,t,2)。對于式(12)中的耗散系數(shù)г，可作為物性參數(shù)，直接在FLUENT的材料面板中輸入，無需編寫UDF程序。

計算的流程如圖1所示，在計算開始之前，需要對用戶通過FLUENT輸入的值或默認值使用UDF定義邊界條件，對所計算的區(qū)域以及用戶自定義標量場進行初始化。再進行循環(huán)計算，同時求解動量方程(5)。通過連續(xù)性方程(1)對速度進行更新，最后計算本構(gòu)方程(6)。在完成一次計算后，對收斂性進行判斷，直到收斂為止，完成計算。

圖1 數(shù)值計算流程圖Figure 1 Flow chart of numerical solution

2.1 可靠性驗證

2.1.1 驗證算例Ⅰ——Oldroyd 本構(gòu)模型的驅(qū)動方腔流

為了驗證本研究所編寫的UDF代碼的正確性，課題組以文獻[20]和文獻[16]中黏彈性流體頂蓋驅(qū)動方腔流動進行同參數(shù)的數(shù)值模擬對比。經(jīng)網(wǎng)格無關(guān)性驗證后網(wǎng)格分布及邊界條件設(shè)置如圖2所示。方腔的高度H=10 mm，長高比為1∶1。方腔內(nèi)部充滿著一種黏彈性流體，方腔頂板以恒定速度U沿x軸正向移動，其余3個壁面均保持靜止。

圖2 幾何模型、網(wǎng)格分布以及邊界條件的設(shè)置Figure 2 Setting of geometric model, grid distribution and boundary conditions

選用二維Oldroyd-B模型作為黏彈性流體的本構(gòu)模型，壁面采用無滑移條件。計算采用雙精度，壓力與速度的耦合采用“SIMPLEC”算法，采用標準格式離散壓力方程，動量方程和用戶自定義標量方程通過采用“QUICK”格式來離散，校驗工況雷諾數(shù)為Re為15，120和333，魏森貝格數(shù)Wi=0.05，溶劑黏性比β=0.3。

計算結(jié)果如圖3所示。本研究計算結(jié)果與文獻中結(jié)果十分吻合。

圖3 計算數(shù)據(jù)對比Figure 3 Validation of numerical results

2.2.2 驗證算例Ⅱ——FENE-P本構(gòu)模型的圓柱繞流

利用類似地方法，課題組編制了基于FENE-P模型的黏彈性流動計算程序，并與文獻[21]中的黏彈性流體圓柱繞流結(jié)果相對比，以驗證程序的正確性。文獻[16]中，圓柱繞流的雷諾數(shù)Re為100，威森貝格數(shù)Wi為0～10，分子最大拉伸無量綱長度L為10和100，聚合物溶液的溶劑黏性比β為0.9。

圖4 計算數(shù)據(jù)對比Figure 4 Validation of numerical results

3 黏彈性流體數(shù)值計算方法研究

課題組選用三維黏彈性流體球形顆粒的力學(xué)特性作為示例算例，進一步說明本研究所提出黏彈性流體數(shù)值計算程序是可靠的。

3.1 模型簡介

研究的對象為懸浮于液相中的單個剛性球狀顆粒，利用相對運動模型進行建模，其模型和邊界條件如圖5(a)所示。

由于來流為黏彈性流體，經(jīng)網(wǎng)格無關(guān)性驗證，表明有關(guān)黏彈性流體流動的數(shù)值研究對網(wǎng)格的精度要求相比于牛頓流體更高。因此，課題組采用106的計算網(wǎng)格數(shù)量來進行數(shù)值模擬，其網(wǎng)格示意圖如圖5(b)所示，流體為基于FENE-P本構(gòu)模型的黏彈性流體。

圖5 三維模型示意圖Figure 3 3D model diagram

3.2 定常問題非定?；惴?/h3>

本研究計算的工況為為較低Re的層流、定常流場，在商用FLUENT軟件中，往往會選擇穩(wěn)態(tài)方法進行流場計算，控制方程式(5)和(6)中，不再含有對時間的偏導(dǎo)項，控制方程形式得到了簡化，整個計算量也相應(yīng)大幅下降。然而實際計算經(jīng)驗表明：利用一般的穩(wěn)態(tài)算法計算這一類黏彈性流動問題，數(shù)值收斂性較差，并且需要在計算過程中根據(jù)實際的計算斂散性情況，不斷地調(diào)節(jié)亞松弛因子，計算才可能達到收斂。因此，若要達到較高計算精度，總的迭代次數(shù)是非常多的，需要很長的計算時間。

課題組的計算實踐表明采用“定常問題的非定?；狈椒梢杂行У亟鉀Q上述問題。所謂“定常問題的非定?；狈椒ǎ侵讣词故嵌ǔＡ鲃訂栴}，仍按照非定常的方程，按時間推進的方法進行計算，只是在每個物理時間步上僅作極少量次數(shù)的內(nèi)迭代。這種處理方法使每個時間步上所獲得的數(shù)值解并非是收斂的物理解，只有在計算物理時間步數(shù)足夠多時，最終的計算結(jié)果才收斂至定常流動的物理解。

圖6所示為基于上述驗證算例Ⅰ和Ⅱ所需迭代步數(shù)，F(xiàn)LUENT中殘差的精度標準均設(shè)置為10-09，Δt為步長，而其數(shù)值計算至收斂的時間相比于穩(wěn)態(tài)計算至少縮短一半。

圖6 總迭代次數(shù)對比圖Figure 6 Comparison chart of total iterations

此種方法對于三維算例同樣適用。但是由于時間步長的取值對于數(shù)值計算的迭代次數(shù)有著極大的影響，若時間步長取得太小，會大大增加計算周期；若時間步長取得太大，又會引起計算發(fā)散，從而無法得到收斂解。因此，同網(wǎng)格無關(guān)性的原理一樣，選擇一個對計算結(jié)果影響較小而又可以滿足穩(wěn)定性要求的合適時間步長也是非常有必要的。

表2中d=2,Re=10,Wi=1時為本研究三維數(shù)值計算的時間步長驗證結(jié)果。不難發(fā)現(xiàn)：當數(shù)值計算可以正常收斂時，采用不同的時間步長并不影響最終的數(shù)值結(jié)果。隨著分子最大拉伸長度L的增大，其合適時間步長會越?。浑S著溶劑黏性比β的降低，其合適時間步長也會減小。綜合計算效率和成本，時間步長0.000 20 s被應(yīng)用于本研究中。在一些較難收斂的工況下，時間步長0.000 01 s也被使用。

表2 計算結(jié)果Table 2 Calculation results

4 結(jié)語

課題組提出了一種新型基于FLUENT中二次開發(fā)模塊UDF的黏彈性流體數(shù)值模擬方法，以O(shè)ldroyd-B模型的二維方腔驅(qū)動流和FENE-P模型的二維圓柱繞流驗證方法的可靠性，并以三維黏彈性流體FENE-P模型的圓球的力學(xué)特性為例，提出加速和穩(wěn)定收斂的方案。結(jié)論如下：

1) 課題組所提出的新型計算方法所編譯的黏彈性流體代碼描述其流場的內(nèi)特性參數(shù)與現(xiàn)有相關(guān)文獻研究結(jié)果吻合度較好，具有較高的準確性。

2) 對于較低雷諾數(shù)定常情況下的黏彈性流體流場計算，可采用定常問題非定常化的處理方案選用合適的時間步長，可有效加快收斂速度，并節(jié)省計算成本。