積鼎科技依托 VirtualFlow 數(shù)值計算仿真平臺，運(yùn)用先進(jìn)數(shù)值模擬方法，可開展氣液兩相流多相動力學(xué)及其特性的數(shù)值模擬研究。在眾多工程應(yīng)用中，高密度和高粘度比的多相流發(fā)揮著關(guān)鍵作用，然而，對這類流動進(jìn)行穩(wěn)健且準(zhǔn)確模擬的方法仍有待完善。本文針對工業(yè)界廣泛應(yīng)用的兩款多相流數(shù)值計算（CMFD）軟件 —— 國外商軟與積鼎科技的 VirtualFlow，在狹小管道兩相流仿真預(yù)報方面展開全面對比。兩款軟件核心差異在于采用不同的常用兩相流模擬方法：國外商軟采用 VOF 方法，VirtualFlow 則采用 Level Set 方法。對比發(fā)現(xiàn)，兩款軟件預(yù)報的流動拓?fù)浣Y(jié)構(gòu)存在顯著不同。對于氣泡流，VirtualFlow 的預(yù)報結(jié)果能夠捕捉到氣泡內(nèi)部的再循環(huán)流動，而國外商軟的 VOF 方法計算結(jié)果中未觀察到明顯的再循環(huán)流動。在段塞流場景下，國外商軟與 VirtualFlow 關(guān)于段塞形成及頻率的結(jié)果偏差顯著。其中，VirtualFlow 的計算結(jié)果呈現(xiàn)出周期性的段塞形成，與前人（Chen 2002）的實驗結(jié)果基本一致；而國外商軟基于 VOF 方法的計算結(jié)果未能捕捉到周期性的段塞形成過程。

 一、引言

近年來，工業(yè)界始終在同時推動微流動工程應(yīng)用部件的性能發(fā)展和小型化發(fā)展，特別在芯片實驗室、生物MEMS和微冷卻電子設(shè)備等領(lǐng)域。在這些部件的微流動通道中，會發(fā)生傳熱和傳質(zhì)過程，可以通過使用多相流來增加傳熱和傳質(zhì)的過程。更進(jìn)一步地深入探索兩相流機(jī)理特性，如界面拓?fù)浣Y(jié)構(gòu)和壓降等方面，可以進(jìn)一步提高微流動工程應(yīng)用部件性能的重要控制參數(shù)的合理性。

近年來，CMFD技術(shù)在工業(yè)界和學(xué)界的應(yīng)用越來越多，CMFD方法中的自由表面跟蹤方法如今已成為商業(yè)軟件中不可或缺的重要組成部分。我們要認(rèn)識到，各種數(shù)值算法和數(shù)值模型目前在這一領(lǐng)域還無法提供物理上100%完全可靠的結(jié)果。我們還需要進(jìn)一步致力于客觀描述并盡量消除數(shù)值仿真在結(jié)果上的不足支出以盡可能獲得可靠的模擬結(jié)果。積鼎科技正在基于目前已相對確定的流動機(jī)理對數(shù)值仿真開發(fā)策略做更宏大的展望，始終致力于建立完善更穩(wěn)健和更準(zhǔn)確的數(shù)值仿真方法并推廣，目標(biāo)超越目前微流體部件設(shè)計的經(jīng)驗主義范式。

本文采用兩個目前較為流行的CMFD軟件：國外商軟與VirtualFlow研究了微型通道中氣液兩相流的特性。這兩款軟件具有不同的自由界面追蹤（Interface Tracking, IT）能力：Fluent采用Volume of Fluid (VOF)方法，而VirtualFlow采用Level Set (LS)方法。VOF和LS是應(yīng)用最廣泛的兩種自由界面追蹤(IT)方法，目前已成功廣泛應(yīng)用于多相流動的數(shù)值模擬領(lǐng)域。Taha等采用國外商軟中的VOF方法研究了上升的Taylor氣泡，并報告了氣相中的強(qiáng)再循環(huán)模式。Akbar、Ghiaasiaan以及Quian、Lawal也用國外商軟中的VOF方法研究了水平通道中的兩相流。他們的計算工況中，氣泡的流動方向與周圍液體的運(yùn)動平行，氣泡與周圍液體處于并流狀態(tài)，數(shù)值結(jié)果顯示氣泡內(nèi)未顯示再循環(huán)流場。Lakehal團(tuán)隊以及Fukagata分別采用VirtualFlow軟件以及其他搭載Level Set方法的CFD軟件進(jìn)行了類似的數(shù)值仿真研究，在他們的計算結(jié)果中都明確發(fā)現(xiàn)了由于氣泡周圍的流動而顯著增強(qiáng)的傳熱。

VOF方法與LS方法作為目前最為流行的兩種多相流自由表面捕捉方法，二者應(yīng)當(dāng)在充分收斂的計算條件下得到流動規(guī)律與仿真質(zhì)量相對一致的結(jié)果。然而，在文獻(xiàn)中，我們發(fā)現(xiàn)對于分別采用VOF和LS方法的水平管道中的兩相流動的數(shù)值仿真，得到的仿真結(jié)果之間似乎存在明顯的差異。

我們仔細(xì)研究了國外商軟和VirtualFlow針對微通道兩相流的數(shù)值仿真結(jié)果，旨在回答以下問題：二者的計算結(jié)果是否再現(xiàn)了隨著氣體體積分?jǐn)?shù)的增加而發(fā)生的預(yù)期流態(tài)變化？建立的模型中的哪些差異可以解釋數(shù)值結(jié)果中觀察到的偏差？導(dǎo)致入口氣體射流破裂的界面不穩(wěn)定性的產(chǎn)生對自由表面追蹤方法數(shù)值計算的敏感度如何？并最終通過詳細(xì)的比較闡述國外商軟中的VOF方法用于預(yù)測段塞流型的數(shù)值方法的不足之處。核心在于其無法預(yù)測不斷增長的界面不穩(wěn)定性的形成以及注入氣體的最終段塞形成。我們建議或許可以采取緩和策略來改進(jìn)針對這類問題的建模。

 二、仿真框架

2.1 本文采用的兩款CMFD軟件的說明

本文的數(shù)值模擬使用了兩種不同的CMFD軟件，分別為ANSYS公司的國外商軟與積鼎科技的VirtualFlow。國外商軟采用VOF方法，而VirtualFlow采用LS方法。VOF和Level Set作為被廣泛使用的兩種自由表面追蹤方法，其各自有繼承了一些有據(jù)可查的、積極的和消極的自身的優(yōu)劣之處，如圖1所示。本文的重點(diǎn)不是通過開發(fā)修正這些劣勢，而是在評估VOF和LS在數(shù)值計算過程中的表現(xiàn)時需要合理地重視并盡量規(guī)避這些弊端。國外商軟與VirtualFlow搭載的都是基于有限體積法求解多流體Navier-Stokes方程的壓力基求解器。當(dāng)然，使用的離散化方案和時間積分方法是不同的。表1列出了模擬中應(yīng)用的方法，VOF（c）中界面的重建是用二階CISAM格式進(jìn)行求解的，而LS函數(shù)（φ）是用QUICK線性迎風(fēng)格式進(jìn)行求解的。

圖1:VOF和LS方法各自的優(yōu)劣之處

  

表1：離散化方案和應(yīng)用的時間積分方法

Scheme description	國外商軟	VirtualFlow
Convection scheme	MUSCL (3rd order)	QUICK(3rd order)
Time integration	Euler (1st order)	Runge–Kutta (3rd order)
P-V coupling	PISO	SIMPLE
P interpolation	PRESTO	Standard
Colour function	CISAM	QUICK
Reinitialization	–	WENO (3rd order)

2.2問題描述

本文針對前人的實驗結(jié)果的基礎(chǔ)上進(jìn)行了二維軸對稱的數(shù)值模擬。將空氣（ρ_g=1.22kg/m³，μ_g=1.78×10^?5kg/ms）和水（ρ_l=998kg/m³，μ_l=0.001kg/ms）注入一根共流水管中。水管截面的外環(huán)圈向管內(nèi)供水，管軸線上的小同心管向管內(nèi)供氣。假設(shè)在t=0時，域中沒有完全發(fā)展的速度剖面，并且可以忽略入口邊緣效應(yīng)。根據(jù)Bretherton推論中的極限，忽略重力效應(yīng)，Bo<0.842（邦德數(shù)定義為重力和表面張力間的比率）。管道的軸向長度為30D，因此完全可以獲得完全發(fā)展的兩相流動。

為了避免國外商軟中VOF模型出現(xiàn)的數(shù)值干燥（管道壁上液膜的數(shù)值破裂，氣體與固體壁面直接接觸），需要細(xì)化近壁網(wǎng)格。使用VirtualFlow的模擬因為沒有出現(xiàn)數(shù)值干涸，所以可以在等距網(wǎng)格間距的情況下進(jìn)行。速度和壓力的收斂標(biāo)準(zhǔn)均設(shè)置為10^?6。邊界條件保持與Lakehal等人詳細(xì)描述的實驗相同，初始條件如圖2所示。

圖2：模擬的初始條件和邊界條件。

表2：入口空隙率α_i，氣體來流速度U_G和液體來流速度U_L

三、計算技術(shù)評估

3.1界面拓?fù)浣Y(jié)構(gòu)和速度剖面

本文研究了兩種不同的兩相流拓?fù)浣Y(jié)構(gòu)，氣泡流和段塞流。氣泡流計算結(jié)果在圖3中展示出來了，根據(jù)直觀印象，兩款軟件的計算結(jié)果似乎生成了相同的流拓?fù)浣Y(jié)構(gòu)(見圖3(a)-(b))，但是氣泡的大小和破裂脫落頻率可以直觀觀察到明顯的差異。國外商軟的VOF方法的模擬結(jié)果相對于VirtualFlow的LS顯示在較低的脫落頻率下產(chǎn)生較大的氣泡。進(jìn)一步觀察氣泡內(nèi)和周圍流場的流動表明，在VirtualFlow的模擬結(jié)果中，氣泡內(nèi)存在再循環(huán)流場。該現(xiàn)象究竟屬于物理解還是數(shù)值上的非物理解還有待通過與實驗結(jié)果的比較來證明。

（a）國外商軟計算結(jié)果在T=0.0076s時的密度等值線

（b）VirtualFlow計算結(jié)果在T=0.0236s時的密度等值線

（c）國外商軟計算得到的T=0.0076s時的壓力云圖

（d）VirtualFlow計算得到的T=0.0236s時的壓力云圖

圖3：國外商軟和VirtualFlow結(jié)果的比較氣泡流，三分之一的管道被抽出，0?10D。

圖3的c圖與d圖中黑色線為氣體積分?jǐn)?shù)c=0.5或LS函數(shù)φ=0的等值線，視作氣液兩相交界線，黑色箭頭線為流場流線。

對于段塞流情況（圖4），國外商軟無法提供具有能夠形成周期性形成段塞的物理規(guī)律的數(shù)值解決方案。段塞的產(chǎn)生似乎是由數(shù)值擾動引起的，可能是由于表面張力處理產(chǎn)生的雜散流動。段塞流和殘差較高的結(jié)果（10^-4，泡狀流，見圖6）與自入口射入的氣體射流的破裂過程具有相同的性質(zhì)。氣體射流在入口處形成細(xì)長的氣體頸部，該氣體頸部折斷并縮回到段塞中。另一方面，VirtualFlow的計算結(jié)果中產(chǎn)生了具有典型特征的周期性段塞流模式。氣液界面處的波紋是微通道中段塞流的典型特征，由氣泡前后曲率的差異導(dǎo)致的表面張力的不平衡造成。這些特性可以在圖中清楚地觀察到。這些結(jié)果與Chen等人實驗結(jié)果呈現(xiàn)出的段塞流狀態(tài)以及Fukagata等和Lakehal等的數(shù)值模擬一致。

（a）國外商軟在T=7.83×10^-3s時的密度等值線（破裂后縮成段塞）

（b）VirtualFlow在T=0.0155s時的密度等值線

（c）Chen等關(guān)于段塞流型的實驗結(jié)果

圖4：對比國外商軟和VirtualFlow關(guān)于段塞流計算結(jié)果(展示1/3管道:5?15D)

3.2 VirtualFlow中的高階時間積分效應(yīng)

為了研究時間積分算法精度格式對計算結(jié)果可能造成的影響，我們使用VirtualFlow進(jìn)行了兩個附加的低階精度格式算法模擬，即Euler一階和Runge–Kutta二階積分算法。由于非定常模擬是在0.1的CFL數(shù)限制下運(yùn)行的，因此沒有觀察到時間積分方案的強(qiáng)烈影響（圖5）。

（a）歐拉一階

（b）龍格-庫塔二階

（c）龍格–庫塔三階

圖5：對比VirtualFlow采用不同的時間積分精度格式得到的速度矢量和流向速度等值線

3.3收斂準(zhǔn)則的影響

利用國外商軟對較高殘差(10^-4)的氣泡流進(jìn)行數(shù)值仿真得到的結(jié)果（圖6）顯示，國外商軟對界面不穩(wěn)定性的預(yù)測對殘差具有很強(qiáng)的敏感性，殘差將直接影響氣體射流破裂的預(yù)測結(jié)果。計算結(jié)果顯示氣體射流將在入口處形成一個細(xì)長的氣體頸，形態(tài)類似從水龍頭滴下的液體，隨后氣泡將以較低的頻率破裂成蠕動的單個氣泡，其拓?fù)浣Y(jié)構(gòu)類似于段塞流。這里給出的結(jié)果與圖3的a、b圖是相同的，都屬于非物理的數(shù)值現(xiàn)象。

圖6：國外商軟計算得到的氣液分布(T=0.0076s)證明殘差對氣泡流拓?fù)浣Y(jié)構(gòu)具有重要影響

四、對擴(kuò)散界面單元特性的處理

分析造成上述差異的部分原因可能是兩款軟件對擴(kuò)散界面單元特性的處理方式不同。液相和氣相之間缺乏剪切應(yīng)力的傳輸可能會在流場中引入誤差。在界面跟蹤方法中，基于界面單元中顏色函數(shù)(c_i，j)的分布，在單元(i, j)中計算材料特性(μ, ρ)。VirtualFlow和國外商軟都使用算術(shù)平均值。

諧波平均在數(shù)學(xué)上更適用于剪切應(yīng)力傳遞描述。

其中剪切應(yīng)力S_ij給出為，

實際上，當(dāng)界面運(yùn)動與外部流動方向垂直時，直接用算術(shù)平均值進(jìn)行處理就是最佳的。但我們發(fā)現(xiàn)，剪切傳遞的缺乏不能僅僅依靠平均的算法來全部解釋。Ho等人用相場法模擬了水平管道中的氣泡流動。他們的研究結(jié)果表明，在界面上應(yīng)用流體物理性質(zhì)的算術(shù)平均值來處理無法捕捉氣相中存在的強(qiáng)烈的再循環(huán)流場對界面的影響。

對比國外商軟的VOF方法，VirtualFlow中的LS方法對界面性質(zhì)卻有不同的處理。在VirtualFlow中，使用兩個相鄰細(xì)胞中心值之間的諧波平均值（基于顏色函數(shù)）來計算單元的表面粘度。由于國外商軟作為商軟不對外開放其底層代碼，所以目前尚不清楚國外商軟中是否使用了相同的方法來計算細(xì)胞表面粘度。在VirtualFlow中，界面寬度是恒定的，最小為一個單元的尺寸。由于顏色函數(shù)的平滑變化，粘度在兩個相位值之間平滑變化。但國外商軟的VOF方法會導(dǎo)致粘度的急劇變化并且會產(chǎn)生更厚的界面（至少兩個單元）。因此，界面單元進(jìn)行數(shù)值計算的流體物理性質(zhì)不同，這將導(dǎo)致界面的厚度會影響相之間的剪切應(yīng)力的傳遞。

我們進(jìn)一步對泡狀流的工況提取了氣體入口的截面流動情況，以描述兩款軟件計算得到的流場和界面剪切應(yīng)力的差異（圖7）。通過比較表明，使用國外商軟可以獲得更平滑的速度場變化。除了粘度的處理外，表面張力項的處理在小規(guī)模流動中也很重要。此外，針對密度的一致處理對于大密度比的情況變得尤其非常重要，尤其是對于壓力校正方程。零發(fā)散條件在VirtualFlow中更容易獲得（每個時間步長的相對收斂條件為10^-6）。國外商軟中的殘差將更難控制，尤其是對于段塞流的情況，即使我們不斷調(diào)整松弛因子也很難將數(shù)值殘差調(diào)整到比較理想的結(jié)果。

（a）兩款軟件計算得到T=7.8×10^-4時入口附近軸向速度對比，黑色線為自由界面

（b）兩款軟件計算結(jié)果在T=7.8×10^-4時入口速度在界面處有不同的速度梯度

圖7：兩款軟件在偏轉(zhuǎn)流場和速度剖面方面計算結(jié)果的對比

 五、結(jié)論

本文對國外商軟與 VirtualFlow 這兩款 CMFD 軟件在氣泡流與段塞流流型的兩相流模擬方面展開系統(tǒng)比較。數(shù)值驗證結(jié)果顯示，國外商軟的 VOF 方法在預(yù)測段塞流準(zhǔn)確度上存在不足，而 VirtualFlow 整體表現(xiàn)更優(yōu)。依據(jù) Chen 等人的實驗結(jié)果，VirtualFlow 的 LS 方法能更準(zhǔn)確地捕捉流動的非穩(wěn)態(tài)模式，在定性上可同時準(zhǔn)確呈現(xiàn)氣泡流和段塞流的流動物理狀態(tài)。

在數(shù)值計算成本方面，國外商軟與 VirtualFlow 在規(guī)則網(wǎng)格上的計算時間基本相當(dāng)。但國外商軟需更高的近壁網(wǎng)格分辨率以避免數(shù)值干涸問題，這意味著在實際工程計算中，國外商軟會產(chǎn)生更高的計算成本。

Nourgraliev 等人的研究指出，當(dāng)自由表面流動由剪切力主導(dǎo)時，幾乎所有擴(kuò)散界面模型都難以準(zhǔn)確預(yù)測所有界面不穩(wěn)定模式的增長過程。本文通過實際算例表明，能否捕捉到剪切力對主要界面不穩(wěn)定性增長的影響至關(guān)重要，因為這直接決定了氣泡和段塞的形成。分析結(jié)果顯示，國外商軟無法生成段塞流的原因之一，可能是其難以確保在每個時間步都達(dá)到高精度的零發(fā)散條件。

眾所周知，準(zhǔn)確預(yù)測多相流型對于精確估計傳熱傳質(zhì)過程意義重大。盡管當(dāng)前針對微尺度傳熱的 CMFD 軟件仍無法提供 100% 完全可靠的結(jié)果，但值得慶幸的是，可借助 PIV 等實驗技術(shù)，滿足學(xué)界及工業(yè)界對物理不穩(wěn)定流型實驗數(shù)據(jù)的大量需求，進(jìn)而進(jìn)一步驗證不同 CMFD 軟件預(yù)測結(jié)果的可靠程度差異。