氣液兩相流問題的高精度數(shù)值模擬

覃柏英， 林賢坤， 榮吉利， 馮志偉

(1. 廣西科技大學(xué) 理學(xué)院，柳州　545006; 2. 廣西科技大學(xué) 汽車與交通學(xué)院，柳州　545006;3. 北京理工大學(xué) 宇航學(xué)院，北京　100081)

氣液兩相流問題的高精度數(shù)值模擬

覃柏英1， 林賢坤2， 榮吉利3， 馮志偉3

(1. 廣西科技大學(xué) 理學(xué)院，柳州545006; 2. 廣西科技大學(xué) 汽車與交通學(xué)院，柳州545006;3. 北京理工大學(xué) 宇航學(xué)院，北京100081)

摘要：采用加權(quán)三階ENN格式離散Euler方程、LS運輸方程及其重新初始化方程的空間導(dǎo)數(shù)，耦合三階Runge-Kutta法離散它們的時間導(dǎo)數(shù)，LS法追蹤運動界面和MGFM 法定義界面邊界條件，探討了氣液兩相流問題的高精度數(shù)值模擬。通過一維大密度大壓力比實驗的數(shù)值模擬，以及二維水下爆炸過程中沖擊波的產(chǎn)生、傳播、反射、透射和水面隆起等演化過程的數(shù)值模擬，驗證了該方法的可行性和穩(wěn)定性，可實現(xiàn)氣液兩相流的高精度數(shù)值模擬及其運動界面的高分辨率追蹤。

關(guān)鍵詞：氣液兩相流；ENN格式；WENN格式；Level Set法；MGFM 法；界面追蹤

因氣液兩種介質(zhì)可壓縮性、密度、壓力等存在巨大差異，氣液兩相流的數(shù)值模擬一直是可壓縮多介質(zhì)流數(shù)值模擬的難點和熱點，其研究重點是運動界面的精確追蹤和界面邊界條件的穩(wěn)健構(gòu)造。對于運動界面的追蹤，目前主要有波前追蹤(Front Tracking)[1]和界面捕捉(Front Capturing)[2]兩大類方法。前者難于處理發(fā)生拓?fù)渥兓慕缑?，且拓展到高維時算法非常復(fù)雜，后者計算簡單，格式統(tǒng)一，易于處理界面的拓?fù)渥兓?。而對于界面邊界條件的構(gòu)造，因界面兩側(cè)流體具有完全不同的特性，密度、速度和壓力在界面處都有可能出現(xiàn)大梯度變化，傳統(tǒng)方法難以實現(xiàn)穩(wěn)健的構(gòu)造。

作為界面捕捉法其中一種的Level Set(LS)法[3]，可處理復(fù)雜和發(fā)生拓?fù)渥冃蔚倪\動界面，且其特征隱含于LS函數(shù)，可較精確描述運動界面。基于LS法，F(xiàn)edkiw等[4]引入Ghost Fluid Method(GFM)，通過界面邊界條件定義，將多介質(zhì)流動問題轉(zhuǎn)化為單介質(zhì)流問題進行求解。但GFM未考慮到界面兩側(cè)流體的性質(zhì)及相互作用，計算強激波與界面作用問題時出現(xiàn)數(shù)值振蕩。劉鐵鋼等[5]通過界面處Riemann 問題的構(gòu)造和求解，由其解定義界面邊界條件，得到引入了較小守恒誤差的Modified Ghost Fluid Method(MGFM)，從而能準(zhǔn)確計算界面處流體狀態(tài)出現(xiàn)大梯度變化的問題。

對可壓縮多介質(zhì)流的數(shù)值模擬，目前已發(fā)展了WENO、NND格式和RKDG有限元法等[6-11]，其中無振蕩無自由參數(shù)耗散的NND格式，具有物理概念清楚，差分格式簡單，易于編程和計算量小等優(yōu)點。為了提高精度，在NND基礎(chǔ)上，張涵信等[12]推導(dǎo)出精度達(dá)三階的ENN格式(ENN3)。

為了進一步提高ENN3的精度，采用加權(quán)方法，可獲得加權(quán)ENN3格式(WENN3)，本文通過實例對其精度進行分析。同時，采用WENN3離散氣液兩相流的Euler方程、LS方程及其重新初始化方程的空間導(dǎo)數(shù)，三階Runge-Kutta法離散時間導(dǎo)數(shù)，耦合LS法追蹤運動界面和MGFM法定義界面邊界條件，對氣液兩相流的高精度數(shù)值模擬和運動界面的高分辨率追蹤進行探討和研究。

1基本方程

1.1Euler方程組

二維可壓縮流體力學(xué)Euler方程組為

(1)

式中，U=[ρ，ρu，ρv,E]T，F(xiàn)=[ρu，ρu2+p，ρuv,u(E+p)]T，G=[ρv，ρuv,ρu2+p，v(E+p)]T。ρ,u,v,E,p分別為密度、x和y方向速度、能量和壓力。

采用Stiffened狀態(tài)方程描述氣液兩種流體，則

p=(γ-1)ρe-γπ

(2)

式中,e為內(nèi)能，π為具有黏性張量量綱的常數(shù)，γ對氣體為比熱比，液體為擬合常數(shù)。不同流體γ和π有確定值。理想氣體γ=1.4,π=0，水γ=7.15,π=3 309。

1.2LS方程及其重新初始化方程

在計算區(qū)域Ω中，氣液分別在區(qū)域Ω1和Ω2，氣液交界面為∑。設(shè)d為t時刻Ω中的點x到∑的距離，則可定義LS函數(shù)φ(x,t)為[13]：當(dāng)x∈Ω1時，φ(x,t)=-d；當(dāng)x∈∑時，φ(x,t)=0；當(dāng)x∈Ω5時，φ(x,t)=d。設(shè)u,v是流體x和y方向速度，則φ(x,t)滿足

(3)

因數(shù)值方法的耗散效應(yīng)，經(jīng)幾個時間步長的迭代后，φ(x,t)不再滿足是x點到界面的符號距離，需重新初始化，可求解如下重新初始化方程到穩(wěn)定解來實現(xiàn)：

(4)

1.3二維轉(zhuǎn)化為一維問題求解

對于二維Euler方程組(1)和方程(3)，可采時間分裂法求解，即將其分裂成兩個一維問題分別求解：

(5)

(6)

由Un，利用式(5)的第一式可求出Un+1/2；再由Un+1/2，利用第二式可求出Un+1。由φn，利用式(6)的第一式可求出φn+1/2；再由φn+1/2，利用第二式可求出φn+1。

對于二維重新初始化方程(4)，將按下式進行求解：

(7)

式中,φy的求解類似于φx，方程(7)的求解類似于一維。

2WENN3有限差分法

2.1算法理論

已知函數(shù)f(x)在各網(wǎng)格中心xi的函數(shù)值f(xi)，設(shè)Δfi+k+1/2=fi+k+1-fi+k，k=-2,-1,0,1，函數(shù)f(x)在右邊界xi+1/2的左狀態(tài)值由ENN3計算[15]：

(8a)

(8b)

其中式(8a)和(8b)中等號右邊的第三部分，都選擇絕對值最小的項。從而

(9)

為了提高精度，對ENN3引入加權(quán)思想，構(gòu)造加權(quán)三階ENN有限差分法(WENN3)，即設(shè)[16]

(10)

(11)

(12)

式中，為了避免分母為零，取ε=10-6。dr為線性權(quán)，βr為光滑因子：

2.2精度分析

(14)

表1　ENN3和WENN3的誤差與精度

由表1可知，ENN3的精度只有三階，而WENN3的精度可達(dá)五階。因此，相比ENN3，WENN3的精度得到了提高，更適合氣液兩相流的高精度數(shù)值模擬及其運動界面的高分辨率追蹤。

3改進GFM法(MGFM)

(15)

圖1　采用MGFM 對氣體各網(wǎng)格的密度、速度和壓力賦值Fig. 1 The grid unit and method to compute normal vector

4數(shù)值模擬

對于一維問題的Euler方程組(5)、LS方程(6)及其重新初始化方程(4)，其空間導(dǎo)數(shù)可采用ENN3或WENN3求解，時間導(dǎo)數(shù)可采用三階Runge-Kutta法提高精度?，F(xiàn)給出求解空間導(dǎo)數(shù)的具體思路與方法。

4.1Euler方程組的求解

(16)

對?F+/?x,?F-/?x進行空間離散，有

(17)

(18)

4.2LS方程的求解

對方程(6)采用凍結(jié)系數(shù)法，并對?φ/?x空間離散，有

(19)

4.3LS值的重新初始化

將一維重新初始化方程(4)寫成如下形式：

(20)

(21)

(22)

5數(shù)值實驗

5.1一維算例

算例1的密度比為1 000，算例2的壓力比為8 000，算例3的密度比和壓力比都為1 000，分別為大密度比大壓力比的氣液兩相流問題。取網(wǎng)格數(shù)為N=500，計算時間分別為T=0,1,0.015,0.001 9,0.02，條件數(shù)為CFL=0.1，3個算例的密度、速度和壓力分布如圖2～圖4所示。

圖2　算例1的密度、速度和壓力分布圖Fig.2 The distribution of density, velocity and pressure of first example

圖3　算例2的密度、速度和壓力分布圖Fig.3 The distribution of density, velocity and pressure of second example

圖4　算例3的密度、速度和壓力分布圖Fig.4 The distribution of density, velocity and pressure of third example

由圖2～圖4可見沖擊波、稀疏波和接觸間斷。相對于ENN3，WENN3的密度、速度和壓力的數(shù)值解與精確解的吻合度很高，在接觸間斷左右側(cè)未出現(xiàn)非物理振蕩，且接觸間斷的耗散寬度和運動界面寬度被控制在4個網(wǎng)格內(nèi)。對于ENN3，算例4和6速度的數(shù)值解與其精確解不吻合，同時，算例6的壓力也不吻合。

5.2二維算例

算例4來自文獻[17]，算例5和6來自文獻[6]。算例4為大壓力比氣液兩相流問題，氣液壓力比近9 000。算例5和6分別是深水和近水面爆炸問題，數(shù)值模擬水中炸藥爆炸時產(chǎn)生的爆轟產(chǎn)物，即高壓氣體與周圍水的相互作用。兩個算例都存在大壓力比，高壓氣體與水壓力比為15 000。同時，算例6還存在大密度比，空氣與水密度比達(dá)720多。取網(wǎng)格數(shù)為N=300×300，條件數(shù)為CFL=0.25，由初始界面開始，三個算例運動界面變化如圖5所示，其中界面間相隔50個時間步長。不同時刻的密度、壓力等值線分別如圖6～圖9所示。

圖5　三個算例的運動界面變化圖Fig. 5 The moving interface changes of three examples

(a) T=0.017 985 0　　　　　(b) T=0.037 994 1　　　　　　(c) T=0.057 989 3　　　　　　(d) T=0.077 500 0圖7　算例4不同時刻的壓力分布圖Fig.7 The pressure distribution in different time of four example

(a) T=0.002 360 7　　　　　(b) T=0.004 820 2　　　　　　(c) T=0.007 508 1　　　　　　(d) T=0.009 787 4圖8　算例5不同時刻的密度分布圖Fig.8 The density distribution in different time of fifth example

(a) T=0.002 360 7　　　　　(b) T=0.004 820 2　　　　　　(c) T=0.007 508 1　　　　　　(d) T=0.009 787 4圖9　算例5不同時刻的壓力分布圖Fig.9 The pressure distribution in different time of fifth example

從圖6～圖9可看出，高壓氣體向外膨脹過程中壓縮其周圍液體，產(chǎn)生向外傳播的沖擊波，并將大部分能量轉(zhuǎn)化為沖擊波能量，同時在氣液交界面處產(chǎn)生向內(nèi)傳播的稀疏波。隨著高壓氣體能量的轉(zhuǎn)化，其能量迅速衰減，沖擊波接收的能量也迅速減少，當(dāng)氣液壓力接近時，沖擊波不再接收氣體能量，開始脫離氣體獨立向外傳播。氣液交界面在高壓氣體作用下繼續(xù)向外膨脹過程中，由于沖擊波向外傳播速度快，致使兩者間距離逐漸增加。從算例4結(jié)果可見，在近相同時刻T=0.057 989 3，本文結(jié)果與文獻[17]結(jié)果相符合。同時，由算例5結(jié)果可見，在近相同時刻T=0.002 360 7，T=0.004 820 2和T=0.009 787 4，本文結(jié)果與文獻[6]結(jié)果相符合。結(jié)果在接觸間斷處仍然存在振蕩，可能是MGFM在定義界面邊界條件引入的誤差導(dǎo)致，但相比文獻結(jié)果，振蕩明顯減小。因此，本文方法較準(zhǔn)確模擬了深水的水下爆炸過程，及其沖擊波的產(chǎn)生和傳播。

從密度和壓力分布圖10和11可看出，對于近水面爆炸問題的算例6，在T=0.003 680 7前，高壓氣體也是向外膨脹過程中壓縮其周圍水，產(chǎn)生向外傳播的沖擊波，同時在氣液交界面處產(chǎn)生向內(nèi)傳播的稀疏波。之后如T=0.004 874 4時，沖擊波與空水交界面發(fā)生相互作用，導(dǎo)致空水交界面開始向上運動，同時在空水交界面處產(chǎn)生反射波和透射波。之后如、T=0.008 654 4時，空水交界面處產(chǎn)生的反射波向下傳播到氣液交界面處，與氣液交界面發(fā)生相互作用，產(chǎn)生反射波和透射波，透射波傳入氣體中，反射波向上運動，與氣液交界面作用，使其加速向上運動；之后如T=0.013 118 3時，因沖擊波在空水交界面和氣液交界面間不斷反射，使得氣液交界面和空水交界面間的距離逐漸減小，并且水面不斷向上隆起。從圖10可見，在近相同時刻，本文結(jié)果與文獻[6]結(jié)果相符合。計算結(jié)果在接觸間斷處也仍然存在振蕩，這也可能是MGFM定義界面邊界條件引入的誤差導(dǎo)致，但相比文獻[6]的結(jié)果，振蕩也明顯減小。因此，本文方法較準(zhǔn)確模擬了近水面爆炸的過程，及其中沖擊波的產(chǎn)生、傳播、反射、透射和水面隆起等演化過程。

(a) T=0.003 680 7　　　　　(b) T=0.004 847 4　　　　　　(c) T=0.008 654 4　　　　　　(d) T=0.013 118 3圖10　算例6不同時刻的密度分布圖Fig.10 The density distribution in different time of sixth example

(a) T=0.003 680 7　　　　　(b) T=0.004 847 4　　　　　　(c) T=0.008 654 4　　　　　　(d) T=0.013 118 3圖11　算例6不同時刻的壓力分布圖Fig.11 The pressure distribution in different time of sixth example

6結(jié)論

本文對ENN3采用加權(quán)方法，獲得WENN3，通過實例對兩者精度進行分析，明確ENN3的精度只有三階，而WENN3的精度可達(dá)到五階。采用WENN3離散Euler方程、LS運輸方程及重新初始化方程的空間導(dǎo)數(shù)，耦合三階Runge-Kutta法離散其時間導(dǎo)數(shù)，LS法追蹤運動界面和MGFM 法定義界面邊界條件，通過大密度大壓力比的數(shù)值實驗，探討了氣液兩相流問題的高精度數(shù)值模擬。算例表明，采用MGFM法定義界面邊界條件，能較準(zhǔn)確計算界面處流體狀態(tài)出現(xiàn)的大梯度變化，而采用LS法，可較精確描述運動界面的變化，實現(xiàn)運動界面的高分辨率追蹤。同時，WENN3保持ENN3優(yōu)點的同時精度也得到提高，能夠很好抑制界面處的非物理振蕩，實現(xiàn)氣液兩相流的高精度數(shù)值模擬。數(shù)值模擬結(jié)果驗證了本文方法的可行性和穩(wěn)定性，較真實模擬了二維水下爆炸過程中沖擊波的產(chǎn)生、傳播、反射、透射和水面隆起等演化過程，可實現(xiàn)氣液兩相流的高精度數(shù)值模擬及其運動界面的高分辨率追蹤，較好抑制了運動界面處的非物理振蕩。

參 考 文 獻

[1] Marshall G. A front tracking method for one-dimensional moving boundary problems[J].SIAM Journal on Scientific Computing,1986,7(1):252-263.

[2] Osher S,Sethian J A. Fronts propagating with curvature dependent speed:Algorithms based on Hamilton-Jacobi formulations. Journal of Computational Physics,1988,79(1):12-49.

[3] Koren B,Lewis M R,van Brummelen E H,et al. Riemann-problem and level-set approaches for homentropic two-fluid flow computations[J]. Journal of Computational Physics,2002,181:654-674.

[4] Fedkiw R P,Aslam T,Merriman B,et al. A non-oscillatory Eulerian approach to interfaces in multimaterial flows(the ghost gluid method)[J]. Journal of Computational Physics,1999,152(2):457-492.

[5] Liu T G,Khoo B C,Yeo K S. Ghost fluid method for strong shock impacting on material interface[J]. Journal of Computational Physics,2003,190(2)：651-681.

[6] 趙海濤,王成,寧建國. 可壓縮多介質(zhì)問題的高精度數(shù)值模擬[J]. 高壓物理學(xué)報,2013(2):261-267.

ZHAO Hai-tao,WANG Cheng,NING Jian-guo. High resolution of compressible multi-medium flow[J].Chinese Journal of High Pressure Physics,2013(2):261-267.

[7] Gryngarten L D,Menon S. A generalized approach for sub-and super-critical flows using the Local Discontinuous Galerkin method[J]. Computer Methods in Applied Mechanics and Engineering,2013,253(1):169-185.

[8] Zhu Jun,Liu Tie-gang,Qiu Jian-xian,et al. RKDG methods with WENO limiters for unsteady cavitating flow[J].Computers and Fluids,2012,57(30):52-65.

[9] Wang Chun-wu,Shu Chi-wang. An interface treating technique for compressible multi-medium flow with Runge-Kutta discontinuous Galerkin method[J]. Journal of Computational Physics,2010,229(23):8823-8843.

[10] Zhu Jun,Qiu Jian-xian,Liu Tie-gang,et al. RKDG methods with WENO type limiters and conservative interfacial procedure for one-dimensional compressible multi-medium flow simulations[J]. Applied Numerical Mathematics,2011,61:554-580.

[11] 魏海鵬,符松. 不同多相流模型在航行體出水流場數(shù)值模擬中的應(yīng)用[J]. 振動與沖擊,2015,34(4):48-52.

WEI Hai-peng,FU Song. Multiphase models for flow field numerical simulation of a vehicle rising from water[J]. Journal of Vibration and Shock,2015,34(4):48-52.

[12] 張涵信,賀國宏,張雷. 高精度差分求解氣動方程的幾個問題[J]. 空氣動力學(xué)學(xué)報,1993(4):347-356.

ZHANG Han-xin,HE Guo-hong,ZHANG Lei. Some important for high order accurate difference solving gas dynamics equations[J]. Acta Aerodynamica Sinica,1993(4):347-356.

[13] Shen Yi-qing,Wang Ru-quan. A fifth-order accurate weighted ENN difference scheme[J]. Communications in Nonlinear Science and Numerical Simulation,1997(2):111-116.

[14] 陳榮三,蔚喜軍. 二維多介質(zhì)可壓縮流的RKDG有限元方法[J]. 計算物理,2006(6):699-705.

CHEN Rong-san,YU Xi-jun. An RKDG finite element method for Two-dimensional compressible multi-media fluids[J]. Chinese Journal of Computational Physics,2006(6):699-705.

High precision numerical simulation for gas-liquid two-phase flow

QIN Bo-ying1, LIN Xian-kun2, RONG Ji-li3, FENG Zhi-wei3

(1. College of Science, Guangxi University of Science and Technology, Liuzhou 545006, China;2. College of Automobile and Transportation Engineering, Guangxi University of Science and Technology, Liuzhou 545006, China;3. School of Aerospace Engineering, Beijing Institute of Technology, Beijing 100081, China)

Abstract:Coupled with the third-order Runge-Kutta, LS and modified ghost fluid method MGFM, the weighted ENN3 scheme (WENN3) was proposed to study the numerical simulation for gas-liquid two-phase flow. WENN3 was applied to discrete the spatial derivative of Euler equation, level set (LS) transport equation and its reinitialization equation. The third-order Runge-Kutta was applied to discrete their time derivative. LS method was applied to track the moving interface and MGFM method was applied to define interface boundary conditions. Through one-dimensional simulations with high density and high pressure ratio, and two-dimensional underwater explosion simulations, the evolution process of shock wave generation, propagation, reflection, transmission and water surface uplift were analyzed, the numerical simulation for gas-liquid two-phase flow were performed, and the correctness and reliability of the proposed method were verified. The results of numerical simulations showed that this method has high resolution, high precision, and many advantages in suimulating gas-liquid two-phase flow with high density, high pressure ratio and strong discontinuity, and it can realize high resolution tracking of moving interface.

Key words:gas-liquid two-phase flow; ENN3; WENN3; level set (LS) method; MGFM; moving interface tracking

基金項目：國家自然科學(xué)基金項目(51209042;1272057)

收稿日期：2014-12-16修改稿收到日期：2015-04-27

通信作者林賢坤 男，博士，教授，1976年生

中圖分類號:O354;O368;O382

文獻標(biāo)志碼：A

DOI:10.13465/j.cnki.jvs.2016.09.013

第一作者 覃柏英 女，碩士，講師，1979年生

E-mail:Linxk0209@163.com