VOF在兩相流交界面捕捉的應用

王健生，鐘易成

(南京航空航天大學 能源與動力學院，江蘇 南京 210016)

0 引言

相界面的捕捉涉及到流體力學、航空航天等眾多領域。運動界面的模擬主要分為兩種類型：界面追蹤和界面捕捉。界面捕捉方法的代表有VOF方法和LevelSet方法。VOF方法是處理界面問題最常用的歐拉方法之一，它通過捕捉網格單元中流體體積的變化來處理任意自由面，由HIRT和NICHOLS提出。PARK I R[1]等提出了用于求解VOF的RHRIC方法，提供了比HRIC方法更高的精度。本文結合不可壓縮有黏RANS方程，采用MHRIC[2]對流通量計算方法求解VOF方程。MHRIC使用的ULTIMATE QUICKEST格式比RHRIC方法精度更高。方程以壓力和速度等作為獨立的原始變量，邊界條件容易處理，方便程序編寫和調試，計算時間短、存儲量小，非常適合研究多項流交界面的運動變化。

1 RANS方程

1.1 主控方程

在低速空氣動力學問題以及水動力相關的液體流動時，流體均可看作不可壓縮的。所以在氣液兩相流的數值模擬中，只考慮不可壓縮的情況。在VOF模型中，認為流體為低速不可壓縮流體，故采用不可壓縮N-S方程算法，求解的變量從[ρuvwT]T變成了[puvwT]T。

不可壓縮多相流體連續方程、動量方程和能量方程分別為：

(1)

對流體控制方程進行無量綱化處理，并通過建立壓強和密度的虛擬關系，在連續方程中引入壓力的虛擬時間導數項，將橢圓型的連續方程變成雙曲型方程，通過偽時間推進求解得到速度散度場，從而使質量和動量都能守恒。N-S方程則變為如下形式：

(2)

式中：τ為虛擬時間；β為虛擬壓縮參數，一般取5～15；Re為雷諾數；Pr為普朗特數。可將上列微分形式的N-S方程寫成積分形式，并通過高斯公式(式3)將體積積分轉換為面積分，得到如式(4)形式的積分形式ALE控制方程。

(3)

(4)

其中：

由于上式中存在黏性相，需額外引入必要的湍流方程。引入多少個附加的湍流量，就要同時求解多少個附加的微分方程，從而使雷諾平均的N-S方程封閉。

1)主控方程離散

對積分形式的控制方程，空間上運用有限體積法在任意網格控制單元上離散，物理時間方向上采用k階后向差分離散，每個物理時間步運用偽時間推進求解，偽時間用一階后向差分，可得到控制體i處離散后的方程為：

(5)

(6)

其中：j為控制體i的鄰居控制體，ij表示控制體i和控制體j的交接面；nface為當前控制體邊界面的數量，如圖1所示。離散方程中包括兩個時間層的推進，一個是物理時間層t、n，是真實物理意義層面的時間；另一個是偽時間層τ、m，其作用在于將每個物理時間步當作一個偽時間層面的定常問題來推進求解。因此，對非定常問題的求解，在每個物理時間步內都嵌套一個偽時間步的迭代循環；偽時間步迭代求解收斂后的結果，即認為是第n個物理時間層的非定常流場，隨即進入下一個物理時間步的求解循環。

圖1 控制體示意圖

數列{φn}用來控制方程離散的物理時間精度，表1列出了幾種取值方式。數列{φn}的取值統一了幾種不同時間精度非定常問題和定常問題的算法。當所有φn為0時，流場計算為定常模式，非零數取到φn-3時，該方法可達到非定常的三階時間精度。

表1 向后差分格式的系數定義

2)對流通量

離散方程中，控制體交界面的對流通量計算為Roe迎風格式：

(7)

如果控制體交界面左右兩側的流場值QL、QR取為左、右控制體中心值，即QL=Qi，QR=Qj，則對流通量的計算為一階空間精度；若要達到二階空間精度，則QL、QR需要由左右控制體梯度重構得到。

(8)

3)黏性通量

對于黏性通量的計算，控制體交界面處的速度、溫度的梯度通過將左右控制體的梯度平均并沿rij修正得到。

(9)

4)隱式格式迭代求解

若殘差項RES由未知的n時間層的值來計算，則求解為隱式格式。隱式求解中，需要對方程式(6)進行線性化處理，即對未知的n時間層的RES泰勒展開：

(10)

忽略高階項后，離散方程可以線性化處理為

(11)

通過迭代求解該方程，即可更新得到下一時間步的流場值。

1.2 湍流模型方程

根據求解的附加微分方程的數目，一般可將渦黏性模式湍流模型劃分為3類：零方程模式或半方程模型、一方程模式和兩方程模式。一方程湍流模型有Baldwin-Barth(B-B)[3]、Spalart-Allmaras(S-A)[4]模型。兩方程湍流模型，如k-ε模型、k-ω模型等。本文使用的湍流模型為S-A模型和k-ω模型。方程均經過無量綱化，寫成積分形式的ALE格式。限于篇幅，不進行詳細敘述。

2 VOF算法

2.1 算法理論

VOF方法的基本原理是通過研究網格單元中流體和網格體積比函數F來確定自由面，追蹤流體的變化，而非追蹤自由液面上質點的運動。VOF方法可以處理自由面重入等強非線性現象，所需計算時間短、存儲量少。VOF方法提出了相界面構造的基本思想，成為運動相界面數值模擬方法的一個全新開端。與MAC方法相比，VOF方法節省了大量的計算機存儲空間，對計算機硬件的要求較低，應用于三維問題計算時優勢會更加明顯。VOF方法的計算相對簡單，相界面銳利程度高，可以表示復雜相界面的結構和變化；在描述復雜相界面和處理三維相界面的融合與破碎問題時，VOF方法遠比其他相界面追蹤方法優越。

假設有N種流體，定義一標量函數——流體體積函數αq，表示計算單元i內某一組分q流體體積Vq占總體積Vi的比率：

(12)

式中：αq=0表示計算單元內沒有流體q；αq=1表示計算單元內只有流體q；0<αq<1表示該計算單元為多種流體混合區域，其中存在交界面。具體如圖2所示。

圖2 αq示意圖

混合流體的特性——密度、黏性、比熱如式(13)所示。

(13)

2.2 VOF輸運方程

對于體積分數αq，滿足以下輸運方程：

(14)

其中：u為流體速度矢量；Sq為組分q的質量源項(默認為0)。其積分形式為

(15)

由于網格均使用靜網格，所以體積對時間的偏導數為0，即

(16)

則輸運方程的積分形式可以表示成

(17)

其中：F(αq)=αqθ；θ=u·n。

與主控方程形式保持一致，增加虛擬時間項，則輸運方程的積分形式為

(18)

2.3 對流通量計算

目前使用較多的對流通量計算主要分為幾何重構格式和空間代數離散格式。本文采用MHRIC格式，該格式是一種復合NVD(normalized variable diagram)格式，其中包含了迎風與順風差分的非線性組合。與其他面通量求解格式相比，例如隱式一階/二階迎風、QUICK[5]格式，MHRIC格式可以得到更精確的結果，同時其計算量遠遠小于幾何重構格式。實現方法如下：

(19)

如圖3所示，其中D為順風單元，C為中心單元，U為迎風單元。

圖3 上游單元的網格結構

對于非結構網格而言，如圖4所示。迎風單元的值可以通過以下方式虛擬構造：反向延長線段CD至U，使得長度UC=CD。U的體積函數值可以用C、D處的值及C處的梯度外插得到：

圖4 MHRIC定義

(20)

(21)

(22)

(23)

(24)

最后可以得到面f上的體積函數：

(25)

通過隱式格式離散后得到的方程為大型稀疏系數矩陣的線性方程組，該方程與章節1.1中RANS方程的求解方式相同，在代碼中體現為調用的函數相同，僅函數中的參數不同。通過迭代求解，本文使用LU-SGS[7]方法求解。得到下一時間步的體積分數αq之后，可以得到混合流體的密度、黏性、比熱，進而求解RANS方程。

2.4 VOF算法隱式求解流程

如圖5所示，VOF算法的求解流程與湍流模型類似。不同的是通量計算格式是采用MHRIC格式。VOF界面捕捉模型求解完成后可獲得各控制體的體積分數，通過體積分數判斷交界面，并更新交界面的密度、黏性和比熱；且交界面影響主控方程的表面張力和體積力計算，因此還需更新表面張力和體積力。

圖5 VOF算法隱式求解過程

3 算例驗證

以經典的NACA0012、油箱晃蕩為數值模型，針對基于虛擬壓縮法的不可壓縮RANS算法和VOF算法進行算例驗證，計算結果與Fluent商業軟件計算結果對比，從收斂性、準確度等方面對本文編寫的數值代碼進行全面的考察。

3.1 二維NACA0012

二維Naca0012翼型較簡單，目的在于考察不可壓縮RANS算法的有效性和收斂性。二維Naca0012模型生成網格如圖6所示，網格數為20 364。

圖6 二維Naca0012翼型網格

NACA0012采用k-w湍流模型，同時考慮流動過程中熱量交換。入口采用壓力遠場邊界，來流速度為102 m/s，來流溫度300 K，兩側的出口邊界類型選擇壓力出口邊界，且出口靜壓為101 000 Pa。除進出口邊界層外，二維Naca0012翼型的上下表面均為壁面邊界且壁面類型為無滑移表面。根據上述數值模擬邊界條件，分別在Fluent和自研軟件CFD2PHASE開展數值模擬求解，當達到穩定后，對兩種軟件的求解結果進行對比如下。計算過程中，CFD2PHASE的殘差分布隨迭代步數的變化情況如圖7所示(本刊為黑白印刷，如有疑問請咨詢作者)。此時在自研軟件中迭代至30 000步以后仍在繼續降低且殘差值進一步降低至10×10-6以下。

圖7 CFD2PHASE的殘差分布示意圖

壓力分布對比情況如圖8所示，兩個軟件求解出的壓力大小及分布情況基本一致，均在機翼前段和尾端達到最大值，在中間的上下表面壓力達到最小。

圖8 壓力分布示意圖

根據上述兩個軟件求解出的不同方向上的速度分布如圖9所示，此時求解出的速度大小、分布基本一致，相差較小。

圖9 速度分布示意圖

提取中間截面位置上壓力分布數值進行對比，所提取的位置在計算區域中的分布如圖10所示。提取出的數據結果對比如圖11所示。進一步對比可知，二者之間相差較小，主要差值出現在壓力較低區域位置。求解結果基本與商業軟件一致。

圖10 壓力數據提取示意圖

圖11 壓力對比結果示意圖

3.2 油箱晃蕩

為了驗證本文提出的VOF算法的效果，首先驗證初始界面構造的準確性與精細性，如圖12所示，其結果均達到了令人較為滿意的程度。對簡化后的油箱在周期速度的作用下進行晃蕩模擬，油箱晃蕩模型設定網格最大尺寸為2 mm，最終生成網格數量為678 906。圖13為網格示意圖。該算例在重力的情況下，選擇層流模型，主相選擇空氣，次相選擇水。

圖12 初始界面

圖13 油箱網格圖

油箱內部充滿空氣，液面高度為400 mm。對計算結果后處理時，在油室中創建一個截面，截面與運動方向平行。得到此截面上密度變化，即氣液交界面隨時間變化情況，如圖14所示。

圖14 油室中氣液交界面隨時間變化情況

圖14是每一秒的氣液交界面圖。模型簡諧振動周期約為1.84 s，振幅隨著時間慢慢增大，可以看出，在1～5 s內，由于運動剛開始振幅較小，液面晃蕩不明顯；5 s過后液面晃蕩較大。為了和Fluent中計算結果對比，同樣在距底面400 m處即初始自由液面處取點1，如圖15所示，Fluent測壓點壓力變化曲線如圖16所示，自研軟件CFD2phase結果如圖17所示。

圖15 測壓點位置

圖16 Fluent中結果

圖17 點1壓力變化圖

對比圖16和圖17發現，兩個計算結果顯示，CFD2phase模擬的點1壓力變化規律基本符合Fluent中壓力結果，由于設置的參考壓力為0，即油箱液面上部空氣區域壓力為0。最低壓力為0的時候，液面由于搖晃低于點1，得到壓力值為0。隨著油箱內部晃蕩越激烈，點1在受到液流沖擊的壓力最大值也越來越大。但本文結果中在每個周期有壓力突變點，壓力突然升高，是因為本文瞬態計算時間步長為0.05 s，步長較長，不能準確反映壓力突變點前后的曲線走勢。點1位置在初始液面上，所以點1初始壓力值在0左右。

4 結語

本文描述了虛擬壓縮法求解不可壓N-S方程的離散求解過程，通過引入壓力的偽時間導數項，將連續方程和動量方程進行聯立，同步推進直接求解壓力場和速度場。求解時只需給出速度的邊界條件以及壓力和速度的初始值，從而避免了大量的求解壓力計算時間，為數值求解創造有利的條件。利用NACA0012壓力遠場算例驗證了二維湍流及遠場邊界模型的可行性。在此基礎上，利用VOF算法，采用MHRIC方式計算面對流通量來捕捉兩相流的交界面，并采用周期速度作用下的油箱內油液面晃蕩模型來驗證本文VOF算法的準確性，所得結果基本正確，體現出本文算法計算時間短，所需內存少的優點。在后期工作中，可將VOF與LevelSet算法耦合形成CLSVOF[8]算法，可進一步提高兩相流交界面捕捉的精度。