基于WENO-THINC/WLIC 模型的水氣二相流數(shù)值模擬1)

韋志龍 蔣 勤

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

引言

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

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

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

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

1 控制方程

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

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

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

式中,m=1,2;αm滿(mǎn)足α1+α2=1.

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

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

2 數(shù)值方法

2.1 分步計(jì)算法

采用分步計(jì)算法[18]對(duì)控制方程進(jìn)行離散求解,計(jì)算過(guò)程分為3 個(gè)階段:對(duì)流項(xiàng)、非對(duì)流項(xiàng)(I)、非對(duì)流項(xiàng)(II).從時(shí)間刻tn到tn+1的具體求解過(guò)程如下.

(1)對(duì)流項(xiàng)

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

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

2.2 WENO 格式

式(5)的一般形式為

式中,? 為變量u或v.

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

2.3 時(shí)間項(xiàng)離散

本文采用三階TVD RK 格式對(duì)時(shí)間項(xiàng)進(jìn)行離散.具體過(guò)程如下

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

2.4 THINC/WLIC 算法

將水氣二相流中水的體積分?jǐn)?shù)輸運(yùn)方程式(3)改寫(xiě)為守恒格式并忽略角標(biāo)可得

式中,? 為水的體積分?jǐn)?shù).

其一維形式為

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

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

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

以上各量需滿(mǎn)足如下條件

ωx和ωy由下式確定

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

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

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

3 數(shù)值模型驗(yàn)證

3.1 界面處理

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

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

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

圖2 和圖 3 分別給出對(duì)應(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 時(shí)刻的模擬結(jié)果,黑色實(shí)線(xiàn)為理論值.可見(jiàn),在光滑界面處數(shù)值模擬結(jié)果與理論值基本吻合,幾乎沒(méi)有形變,在界面轉(zhuǎn)折突變處隨著界面運(yùn)動(dòng)而逐漸光滑.在網(wǎng)格加密后,這種界面形變得到了很好的抑制.由此可見(jiàn),模擬得到的Zalesak’s disk 的形狀在旋轉(zhuǎn)過(guò)程中與理論值基本一致,表明本數(shù)值模型具有很好的界面追蹤能力.為定量分析界面模擬精度,定義誤差計(jì)算公式為

圖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)的體積分?jǐn)?shù),為網(wǎng)格i內(nèi)的體積分?jǐn)?shù)理論參考值.由于網(wǎng)格尺寸倍減,可采用下式計(jì)算界面模擬精度

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

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

(2)Shearing flow 問(wèn)題

Shearing flow 問(wèn)題是驗(yàn)證大變形情況下數(shù)值方法追蹤界面能力的另一典型算例.通過(guò)模擬初值為圓心在(0.5,0.75)、半徑為0.15 的圓,其內(nèi)部體積分?jǐn)?shù)為1,外部為0 條件下,在隨時(shí)間變化的外加速度場(chǎng)驅(qū)動(dòng)下圓盤(pán)的形狀變化來(lái)檢驗(yàn)數(shù)值模型的界面計(jì)算精度.其中,外加速度場(chǎng)用如下的流函數(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

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

3.2 線(xiàn)性液艙晃蕩問(wèn)題

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

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

其中

圖6 線(xiàn)性液艙晃蕩問(wèn)題示意圖Fig.6 Setting of the linear sloshing problem

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

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

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

定義

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

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

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

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

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

3.3 潰壩問(wèn)題

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

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

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

圖10 不同時(shí)刻潰壩水體運(yùn)動(dòng)形態(tài)對(duì)比Fig.10 Comparison of dam-breaking water shapes at different times

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

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

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

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

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

圖13 P1 處壓強(qiáng)歷時(shí)曲線(xiàn)Fig.13 Pressure history at P1

圖14 P2 處壓強(qiáng)歷時(shí)曲線(xiàn)Fig.14 Pressure history at P2

4 結(jié)論

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

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

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

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

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