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

(1.上海理工大學(xué) 能源與動(dòng)力工程學(xué)院， 上海 200093；2.上海理工大學(xué) 上海市動(dòng)力工程多相流動(dòng)與傳熱重點(diǎn)實(shí)驗(yàn)室， 上海 200093)

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

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

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

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

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

1 黏彈性流體控制方程

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

連續(xù)方程：

?ui/?xi=0。

(1)

動(dòng)量方程：

(2)

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

(3)

(4)

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

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

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

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

2.1 UDF 代碼編譯

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

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

(5)

(6)

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

(7)

式中Fi為方程的源項(xiàng)。

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

(8)

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

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

(9)

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

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

(10)

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

φ=Cij；

(11)

г=ρκ；

(12)

(13)

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

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

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

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

2.1 可靠性驗(yàn)證

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

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

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

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

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

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

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

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

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

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

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

3.1 模型簡(jiǎn)介

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

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

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

3.2 定常問(wèn)題非定?；惴?/h3>

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

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

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

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

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

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

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

4 結(jié)語(yǔ)

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

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

2) 對(duì)于較低雷諾數(shù)定常情況下的黏彈性流體流場(chǎng)計(jì)算，可采用定常問(wèn)題非定?；奶幚矸桨高x用合適的時(shí)間步長(zhǎng)，可有效加快收斂速度，并節(jié)省計(jì)算成本。