韋志龍 蔣 勤

(河海大學港口海岸與近海工程學院,南京 210024)

引言

二相流是自然界中常見的流體運動形式.其中,水氣二相流與水利工程、海岸工程和海洋工程等實際工程問題密切相關(guān).水氣二相流因具有界面變形復雜、密度差大和紊動強烈等特點,需要較高的界面追蹤和控制方程對流項的數(shù)值求解精度,是計算流體力學的難點和熱點問題.相關(guān)研究結(jié)果表明,當流體速度低于0.3 倍的聲速時,流體運動可視為不可壓縮流進行處理[1].因此,在對開敞水域的自由表面流工程問題進行數(shù)值模擬分析時,為簡化模型,一般情況下可以忽略空氣的可壓縮性,將其視為不可壓縮流體.據(jù)此,本研究的水氣二相流的模型亦采用水和空氣的不可壓縮假定.

如上所述,水氣二相流中界面兩側(cè)的介質(zhì)密度相差較大,而密度的微小變動將會引起界面兩側(cè)壓強的劇烈變化;此外,水氣二相流還具有紊動劇烈的特點;這些均要求所選取的對流項的求解方法在具有較高的精度的同時,能夠很好地處理變量的突變問題.加權(quán)基本無震蕩(weighted essentially nonoscillatory,WENO)格式[2]根據(jù)變量變化曲線的光滑程度選取加權(quán)值.對于2r?1 階WENO 格式,在曲線光滑段具有2r?1 階精度,在間斷點處退化為基本無震蕩(essentially non-oscillatory,ENO)格式[3],并保持基本無震蕩特性,被廣泛用于實際工程中流體運動問題,特別是空氣動力學問題的數(shù)值模擬[4-6].相對而言,WENO 格式在水氣二相流中的應(yīng)用成果較少.鑒于其高精度、自適應(yīng)、標準化等特點,本文選取WENO 格式離散對流項,同時采用三階總變差遞減(total variation diminishing,TVD)龍格庫塔(Runge-Kutta,RK)法[7]對時間項進行離散,求解流體運動的動量方程.

另一方面,對二相流的界面處理以VOF[8](volume of fluid)法應(yīng)用最為廣泛[9-11],VOF 模型理論上具有質(zhì)量守恒、物理意義明確等特點,可分為幾何重構(gòu)法和代數(shù)法[12].在多維情況下,前者的界面幾何重構(gòu)過程過于復雜且計算效率低下,相比而言代數(shù)法對界面追蹤的計算過程更為簡單.Xiao 等[13]提出的雙曲正切函數(shù)界面捕捉法(tangent of hyperbola for interface capturing,THINC)為典型的代數(shù)法.該方法采用分段雙曲正切函數(shù)重構(gòu)界面,將流體各物理量沿界面法向方向的變化視為連續(xù)函數(shù),通過模型中的界面形態(tài)參數(shù)控制不同介質(zhì)界面的厚度,進而確定界面網(wǎng)格的體積分數(shù)通量,不僅大大簡化了界面追蹤的計算流程,而且能夠保證流體質(zhì)量守恒的計算精度.此外,Yokoi[14]提出了將加權(quán)線性界面計算(weighted line interface calculation,WLIC)和THINC 法相結(jié)合的改進方法,以取代之前將一維THINC 法直接用于多維界面運動問題求解的作法,有效地提高了THINC 法對多維界面追蹤的計算精度.THINC/WLIC法的計算過程簡單,適用于二維和三維問題的流體界面追蹤.除此之外,還有一些對THINC 法的改進方法[12,15-17],但大多包含多維雙曲正切函數(shù)的重構(gòu),涉及的參數(shù)過多,程序?qū)崿F(xiàn)較困難,且計算量亦相對較大.

迄今,尚未見到國內(nèi)外關(guān)于將WENO 格式與THINC 法進行耦合的相關(guān)模型研究成果.考慮到WENO 格式的高精度、自適應(yīng)等特點適用于二相流的對流項的計算,THINC/WLIC 法在自由表面追蹤時既簡單易行又具有較高的求解精度的特點,適合于復雜界面的追蹤計算,本文嘗試將兩者進行耦合,建立WENO-THINC/WLIC 耦合水氣二相流數(shù)值模型,并通過對典型流體力學問題的模擬分析與模型驗證,探索適用于大多數(shù)工程領(lǐng)域中水氣二相流運動模擬的高精度數(shù)值計算方法.

1 控制方程

本模型中水氣兩相均視為不可壓縮,流體運動的控制方程為NS 方程,其二維形式可表示為

式中,i,j=1,2;ui為速度分量;fi為質(zhì)量力,在本模型中為重力;τi j為總應(yīng)力,對于牛頓流體有τi j=?pδi j+2μSij?2μδijSkk/3,其中p為壓強,Sij=(?ui/?xj+?uj/?xi)/2,δi j為Kronecker 算子.

基于VOF 方法追蹤水氣交界面,水和空氣的體積分數(shù)αm滿足方程

式中,m=1,2;αm滿足α1+α2=1.

流體的密度和黏性系數(shù)通過下式計算求得

式中,λ 為密度ρ 或黏性系數(shù)μ,λ1和λ2分別為對應(yīng)的水體或氣體的特性參數(shù).

2 數(shù)值方法

2.1 分步計算法

采用分步計算法[18]對控制方程進行離散求解,計算過程分為3 個階段:對流項、非對流項(I)、非對流項(II).從時間刻tn到tn+1的具體求解過程如下.

(1)對流項

方程右側(cè)壓強梯度?pn+1/?xi使用中心差分法離散.

以上為單個時間步的計算流程,單步計算結(jié)束后返回步驟(1),根據(jù)設(shè)定的CFL 條件更新?t,如此循環(huán)可以得到設(shè)定時段的計算結(jié)果.

2.2 WENO 格式

式(5)的一般形式為

式中,? 為變量u或v.

現(xiàn)以? 在x方向的離散格式為例說明WENO 格式的計算過程

2.3 時間項離散

本文采用三階TVD RK 格式對時間項進行離散.具體過程如下

式中,? 為基本變量,本模型中為ui;L 為式(1) 和式(2) 中對流項和等號右側(cè)對應(yīng)的數(shù)值解;角標(1)和角標(2)為當前時間步tn和下一時間步tn+1之間的子時間步.

2.4 THINC/WLIC 算法

將水氣二相流中水的體積分數(shù)輸運方程式(3)改寫為守恒格式并忽略角標可得

式中,? 為水的體積分數(shù).

其一維形式為

在WLIC 方法中,多維界面重構(gòu)是根據(jù)界面法向量分別對水平界面和豎向界面賦權(quán)重來實現(xiàn)的.在主方向使用一維THINC 法計算體積分數(shù)通量,在次主方向使用簡單線性界面法(simple line interface calculation,SLIC) 計算,從而實現(xiàn)計算量和精度之間的平衡,以較小的計算量獲取較高的精度.如圖1 所示,設(shè)在網(wǎng)格(i,j)內(nèi),體積分數(shù)為αi,j.根據(jù)WLIC 理論,該網(wǎng)格內(nèi)界面重構(gòu)表達式為

式中,n為界面法向量,χ 為界面特征表達式,χx和χy分別為豎向界面和水平界面的特征表達式,ωx和ωy為對應(yīng)的權(quán)重.

圖1 界面加權(quán)分解Fig.1 Weighted decomposition of the interface

以上各量需滿足如下條件

ωx和ωy由下式確定

式中,nx,i,j和nx,i,j分別是界面法向量n的x方向和y方向分量,具體計算過程見文獻[14].

待權(quán)重確定后,通過網(wǎng)格邊壁(i+1/2,j)的通量的計算公式如下

其中,式(33) 采用一維THINC 法計算,式(34) 采用SLIC 計算.網(wǎng)格邊壁(i,j+1/2)的通量計算過程類似.

3 數(shù)值模型驗證

3.1 界面處理

采用本研究建立的數(shù)值模型對Zalesak’s disk 旋轉(zhuǎn)[21]和shearing flow 問題[22]進行模擬分析,以驗證本數(shù)值模擬方法對大變形界面的追蹤能力.其中,驗證案例的計算域均為[0,1]×[0,1],均采用100×100,200×200,400×400 和800×800 四種密度的結(jié)構(gòu)化網(wǎng)格進行模擬計算,并設(shè)置CFL=0.25.

(1)Zalesak’s disk 旋轉(zhuǎn)問題

Zalesak’s disk 旋轉(zhuǎn)問題常被用來檢驗界面處理方法的適用性和模擬精度.如圖2 和圖3 中黑色實線所示,Zalesak’s disk 是一個帶有缺口的圓形,在外加速度場(u,v)=(y?0.5,0.5?x)的作用下繞點(0.5,0.5)勻速旋轉(zhuǎn),在t=2 時回到原處.Zalesak’s disk 在計算域上的初值定義為

圖2 和圖 3 分別給出對應(yīng)于 200×200 和400×400 兩種網(wǎng)格的模擬結(jié)果.圖中a,b,c,d 分別為t=0.5,t=1.0,t=1.5,t=2.0 時刻的模擬結(jié)果,黑色實線為理論值.可見,在光滑界面處數(shù)值模擬結(jié)果與理論值基本吻合,幾乎沒有形變,在界面轉(zhuǎn)折突變處隨著界面運動而逐漸光滑.在網(wǎng)格加密后,這種界面形變得到了很好的抑制.由此可見,模擬得到的Zalesak’s disk 的形狀在旋轉(zhuǎn)過程中與理論值基本一致,表明本數(shù)值模型具有很好的界面追蹤能力.為定量分析界面模擬精度,定義誤差計算公式為

圖2 Zalesak’s disk 數(shù)值模擬結(jié)果(網(wǎng)格密度:200×200)Fig.2 Simulation results of Zalesak’s disk with 200×200 mesh

圖3 Zalesak’s disk 數(shù)值模擬結(jié)果(網(wǎng)格密度:400×400)Fig.3 Simulation results of Zalesak’s disk with 400×400 mesh

式中,N為網(wǎng)格數(shù)量,H(x)dx為數(shù)值模擬得到網(wǎng)格i內(nèi)的體積分數(shù),為網(wǎng)格i內(nèi)的體積分數(shù)理論參考值.由于網(wǎng)格尺寸倍減,可采用下式計算界面模擬精度

式中,Eq表示qth對應(yīng)網(wǎng)格的模擬誤差.不同網(wǎng)格尺度下t=2.0 時刻的本模型模擬誤差見表1.由表可知,采用THINC/WLIC 法本模型對Zalesak’s disk 旋轉(zhuǎn)問題的界面計算精度達到了一階精度.

表1 Zalesak’s disk 數(shù)值模擬誤差Table 1 Numerical errors in the Zalesak’s disk test

(2)Shearing flow 問題

Shearing flow 問題是驗證大變形情況下數(shù)值方法追蹤界面能力的另一典型算例.通過模擬初值為圓心在(0.5,0.75)、半徑為0.15 的圓,其內(nèi)部體積分數(shù)為1,外部為0 條件下,在隨時間變化的外加速度場驅(qū)動下圓盤的形狀變化來檢驗數(shù)值模型的界面計算精度.其中,外加速度場用如下的流函數(shù)表示

分析流函數(shù)的時間變化可知,在0

圖4 Shearing vortex 數(shù)值模擬(網(wǎng)格密度:200×200)Fig.4 Simulation of shearing vortex with 200×200 mesh

圖5 Shearing vortex 數(shù)值模擬(網(wǎng)格密度:400×400)Fig.5 Simulation of shearing vortex with 400×400 mesh

表2 Shearing flow 數(shù)值模擬誤差Table 2 Numerical errors in the shearing flow test

此外,對界面追蹤模擬過程中計算域內(nèi)總質(zhì)量的變化進行了計算分析.設(shè)定M=∑αi,M0為初始時刻計算域內(nèi)總質(zhì)量.在不同網(wǎng)格尺度條件下,本模型得到的質(zhì)量變化比|M?M0|/M0均保持在10?13以下,可以認為THINC/WLIC 界面追蹤法基本上可以保持質(zhì)量守恒.

3.2 線性液艙晃蕩問題

根據(jù)線性勢流理論可以得到線性液艙晃蕩問題的解析解,該解可用來檢驗數(shù)值計算模型的收斂精度.如圖6 所示,假定初始時刻水體和氣體都處于靜止狀態(tài),艙壁為不可滲可滑移邊界.在t>0 后艙體以a的加速度水平向左運動,相當于艙內(nèi)流體在受到豎直向下的重力g還受到水平向右、大小為a的質(zhì)量力的作用.計算中,取幾何尺寸和物理參數(shù)與Li等[23]的模擬案例一致.具體地,艙體寬度L=1.00 m,初始水深h1=1.00 m,初始氣體高度h2=1.25 m,重力g=9.81 m/s2,水平加速度a=0.01g,水體密度ρ1=1000 kg/m3,氣體密度ρ2=1 kg/m3.

根據(jù)線性勢流理論,水面高程 η(x,t) 的解析解[24]為

其中

圖6 線性液艙晃蕩問題示意圖Fig.6 Setting of the linear sloshing problem

分別選取64 × 144,128 × 288 和256 × 576 三種均勻網(wǎng)格,并設(shè)置CFL=0.1 進行模擬.需要指出的是,為了與理論解進行比較,在數(shù)值模擬時忽略了控制方程中黏性項的影響.圖7 給出了本模型得到的艙體左右邊壁處的水面高程(α=0.5) 時間變化的計算結(jié)果與解析解的比較.其中,64 × 144 解析度對應(yīng)的網(wǎng)格尺寸約為0.016 m × 0.016 m,大于水面高程振幅,此時模型模擬結(jié)果已較準確地再現(xiàn)了水面高程的變化趨勢,說明THINC/WLIC 法能較好地模擬微小的水面變形;256×576 解析度網(wǎng)格對應(yīng)的網(wǎng)格尺寸約為0.004 m × 0.004 m,由圖7 可見,模型模擬結(jié)果與解析解基本吻合.

圖7 艙體左右邊壁處水面高程變化Fig.7 Temporal profile of the interface elevations on the left and right walls

為分析模型的計算精度,增加512×1152 網(wǎng)格算例并將t=4 s 時刻的密度場與理論值進行對比分析.

定義

表3 線性液艙晃蕩問題數(shù)值計算誤差Table 3 Numerical errors in the linear sloshing test

本模型采用有限差分法對控制方程進行求解,而有限差分法在數(shù)值守恒性方面偏弱,為確定本模型在質(zhì)量、動量和能量守恒方面的表現(xiàn),進行了以下的分析.定義計算域總質(zhì)量M=∑ρi?x?y,M0為初始時刻計算域內(nèi)總質(zhì)量,定義質(zhì)量損失為(|M?M0|)/M0.計算域內(nèi)總能量為

圖8 為對上述4 個特征變量數(shù)值守恒性的評估結(jié)果.其中,黑色線條代表質(zhì)量損失,綠色線條代表能量損失,藍色線條代表水平方向動量損失,紅色線條代表豎直方向動量損失.

由圖8 可知,數(shù)值質(zhì)量損失保持在10?10量級以下,可認為本案例中數(shù)值計算模型滿足質(zhì)量守恒;數(shù)值能量損失在10?4量級以下,且隨著網(wǎng)格加密能量損失減小.由此可見,本耦合模型在流體的質(zhì)量守恒方面具有良好表現(xiàn),在流體的動量和能量守恒方面與采用有限體積法得到模擬結(jié)果相比略有不足,但采用較高精度的計算網(wǎng)格會提高模擬結(jié)果的能量和動量守恒性.

圖8 質(zhì)量、能量及動量的模擬誤差Fig.8 Simulation errors on the conservation of mass,energy and momentum of fluid

3.3 潰壩問題

為驗證本研究提出的WENO-THINC/WLIC水氣耦合數(shù)值模型的模擬精度,針對干床面潰壩問題進行了數(shù)值模擬,并與Lobovsk′y 等[25]的物理模型試驗結(jié)果進行了比較.如圖9 所示,數(shù)值模擬的計算條件與Lobovsk′y 等的物理模型試驗相同.水槽長為1.61 m,寬0.60 m.在距離水槽左邊壁600 mm 處設(shè)置一豎直擋板,擋板左側(cè)為水體,右側(cè)為空氣.水體的高度包括高水位H=600 mm 和低水位H=300 mm兩種工況.物理模型試驗中,擋板上接定滑輪,連接一重物,在重物的作用下?lián)醢灞怀殡x,左側(cè)水體涌向右側(cè),形成潰壩流.參照Lobovsk′y 等的模型試驗,數(shù)值模擬中設(shè)置了與物模試驗相同的物理量監(jiān)測點,包括右側(cè)邊壁上的兩個壓力監(jiān)測點和水槽中部4 個水位監(jiān)測點.其中,壓力監(jiān)測點P1 距水槽底部30 mm,P2 距水槽底部80 mm,水位監(jiān)測線具體位置見圖9.

圖9 潰壩模擬試驗示意圖(單位:mm)Fig.9 Setting of the dam-breaking simulation test(unit:mm)

首先進行了潰壩過程中水流運動學特征的模擬結(jié)果與試驗結(jié)果的對比,包括潰壩水流形態(tài)、前鋒位移和沿程的水位變化等.其中,相關(guān)數(shù)值模擬案例均采用?x=?y=0.01 m 的結(jié)構(gòu)化網(wǎng)格并設(shè)置CFL=0.1.圖10 為不同時刻潰壩水體運動形態(tài)的模擬結(jié)果與試驗結(jié)果.由圖可見,在低水位H=300 mm情況下,數(shù)值模擬的水體運動形態(tài)與物理模型試驗結(jié)果基本一致,潰壩水體前鋒撞擊右側(cè)邊壁后形成的水舌及空氣空腔的形狀和位置與物模試驗結(jié)果也基本相同.

圖10 不同時刻潰壩水體運動形態(tài)對比Fig.10 Comparison of dam-breaking water shapes at different times

圖11 為潰壩水體前鋒位移隨時間變化曲線.將模擬得到的潰壩前鋒位移與物理模型試驗結(jié)果[25-27]進行定量對比可知,在潰壩水體前鋒沖擊右側(cè)邊壁前,數(shù)值模擬結(jié)果與試驗結(jié)果吻合良好,說明本模型對潰壩前期水體形狀的模擬精度較高.

圖11 潰壩前鋒線位移計算結(jié)果與試驗結(jié)果比較Fig.11 Comparison of calculated and measured displacements of the wave front

與潰壩水體前鋒位移相比,H1 ～H4 四個監(jiān)測點處的水位變化能給出更詳細的潰壩水體沿程變化過程.考察從無量綱時間t?=0 到t?=10 內(nèi)各監(jiān)測點水面變化,可以看到潰壩波撞擊邊壁后形成的水舌和二次波的變化特征.圖12 給出了H=300 mm 條件下模擬得到的4 個水位監(jiān)測點處水位變化與前人物理模型試驗研究結(jié)果[25]及其他數(shù)值模擬結(jié)果[28-29]的比較.由圖可見,本模型得到的潰壩前鋒首次到達各水位監(jiān)測點的時間與物模試驗一致;同時,二次波到達的時間和水位上升趨勢與試驗結(jié)果基本吻合.本模型和其他數(shù)值模型結(jié)果在潰壩波撞擊邊壁前均與物理模型試驗結(jié)果較為吻合;不同點在于,本模型對潰壩水體沖擊右側(cè)壁面以后的后期水體運動模擬結(jié)果相比其他數(shù)值模型更為接近物理模型試驗結(jié)果.Peltonen 等[28]使用的OpenFOAM 內(nèi)置的標準VOF求解器interFOAM 所得結(jié)果在H3 測點后期出現(xiàn)水位嚴重偏大現(xiàn)象、H4 測點后期出現(xiàn)水位偏小現(xiàn)象,其改進的GFMFOAM (ghost fluid method FOAM) 方法在H2 測點后期出現(xiàn)水位嚴重偏小現(xiàn)象;Nguyen 和Yark[29]改進的VOF 法在H2 測點處存在非物理的觸頂現(xiàn)象、H3 測點中后期出現(xiàn)水位明顯偏大的問題.綜合來看,潰壩問題在水體前鋒撞擊右側(cè)邊壁后,水氣界面形變復雜,水氣紊動更為劇烈,此階段數(shù)值模擬較為困難.相比而言,本模型對潰壩水體整個運動過程的模擬結(jié)果與物理模型試驗結(jié)果更為吻合.可見,本數(shù)值計算模型較好地描述了潰壩水體的運動學特征,能夠較為準確地再現(xiàn)潰壩水體這種大變形自由表面流的運動過程.

圖12 H1 ～H4 處水位變化Fig.12 Water level elevations at the four locations H1 ～H4

此外,潰壩水體對下游床面或建筑物的沖擊力是潰壩問題研究中的另一個重要物理參數(shù).通過考察高水位H=600 mm 條件下下游壓強的計算結(jié)果,對本模型再現(xiàn)潰壩水體動力學特征的能力進行分析.其中,設(shè)定計算網(wǎng)格分別為?x=?y=10 mm和?x=?y=5 mm 兩種情況,對應(yīng)的網(wǎng)格數(shù)量為60×161 和120×322,并取CFL=0.1,將模擬得到的P1 和P2 兩點的壓強與試驗結(jié)果[25]和他人數(shù)值結(jié)果[30]進行對比.如圖13 和圖14 所示,當潰壩波未到達右側(cè)邊壁時,監(jiān)測點處的壓強為大氣壓,相對壓強為0;當潰壩前鋒撞擊邊壁時,壓強急劇增大,之后緩慢變小.沿右邊壁,沖擊壓強峰值隨高度上升而下降,即P1 點沖擊壓強峰值應(yīng)大于P2 點峰值.本模型得到結(jié)果與物理模型試驗相比,沖擊壓強峰值出現(xiàn)時間相同,峰值大小相當,下降趨勢也大致相同.考慮到試驗數(shù)據(jù)為典型壓力變化,存在一定的波動空間,可認為本模型數(shù)值模擬壓強結(jié)果與物理模型試驗結(jié)果較為一致.與Sun 等[30]使用MPS (moving particle semiimplicit)法的結(jié)果相比,本模型結(jié)果的壓力波動較小.此外,網(wǎng)格尺寸對壓強影響較小,但采用低密度網(wǎng)格的模擬結(jié)果出現(xiàn)少許壓力局部突變,網(wǎng)格加密后這種現(xiàn)象減少.總體而言,使用網(wǎng)格尺寸?x=?y=0.01 m得到的計算結(jié)果已與試驗結(jié)果吻合較好,也就是說本模型通過較少的網(wǎng)格數(shù)即可得到較精確的解.

圖13 P1 處壓強歷時曲線Fig.13 Pressure history at P1

圖14 P2 處壓強歷時曲線Fig.14 Pressure history at P2

4 結(jié)論

本研究通過采用五階WENO 格式求解對流項,三階TVD RK 進行時間離散,THINC/WLIC 法追蹤界面,建立了WENO-THINC/WLIC 耦合水氣二相流數(shù)值計算模型.首先,選取Zalesak’s disk 旋轉(zhuǎn)和shearing flow 問題,驗證THINC/WLIC 法對外加速度場下大變形水氣界面的追蹤能力;其次,以線性液艙晃蕩問題為例分析本研究提出的WENO-THINC/WLIC 耦合模型的計算誤差以及對質(zhì)量、動量和能量的數(shù)值守恒性;最后,以潰壩水流運動問題為例檢驗耦合模型模擬不可壓縮水氣二相流大變形運動問題的能力.主要結(jié)論如下:

(1) THINC/WLIC 法可有效捕捉復雜界面變形,模擬收斂精度可達一階并保持物質(zhì)質(zhì)量守恒.

(2) 所建立的WENO-THINC/WLIC 耦合模型對線性液艙晃蕩問題的模擬結(jié)果與解析解吻合度較高,對微小形變的捕捉較為準確;通過誤差分析可知模型模擬精度為一階,數(shù)值分析結(jié)果表明本模型保持質(zhì)量守恒和能量守恒.

(3)耦合模型對潰壩水流運動的計算結(jié)果與物理模型試驗結(jié)果吻合良好,并優(yōu)于interFOAM、VOF 以及MPS 等數(shù)值模型的模擬結(jié)果,更夠準確地捕捉了潰壩過程中潰壩水體的運動學特征和動力學特征,證明本模型對不可壓縮水氣二相流問題的數(shù)值模擬具有較好的適用性.

將改進的WENO 格式[23,31-36]和THINC 法代入本研究提出的WENO-THINC 耦合模型預期會得到更高精度的模擬結(jié)果.