李泠泉，楊愛明

(復(fù)旦大學(xué) 航空航天系，上海 200433)

?

基于HLLC近似Riemann解的自由面流動(dòng)數(shù)值模擬

李泠泉，楊愛明

(復(fù)旦大學(xué) 航空航天系，上海 200433)

摘要：對(duì)于自由面流動(dòng)問題，發(fā)展了一種新的數(shù)值模擬方法.在該方法中，自由面的捕捉通過求解Level Set方程實(shí)現(xiàn).該方程與流體力學(xué)控制方程耦合形成雙曲系統(tǒng)方程，并采用人工壓縮方法進(jìn)行求解.為了提高計(jì)算的分辨率，發(fā)展了一種有限體積框架下的高階WENO格式，并構(gòu)造了基于HLLC近似Riemann算子的通量函數(shù).對(duì)若干典型自由面流動(dòng)問題進(jìn)行了數(shù)值模擬，計(jì)算結(jié)果表明了該方法的有效性.

關(guān)鍵詞：自由面流動(dòng)； Level Set方程； WENO格式； HLLC近似Riemann算子； 人工壓縮方法

很多自然科學(xué)與工程技術(shù)問題都涉及自由面流動(dòng)問題，比如貼近海面飛行的飛行器與自由面相互作用的問題、船舶的水動(dòng)力學(xué)問題等.含自由面的流動(dòng)一般比較復(fù)雜，自由面的拓?fù)浣Y(jié)構(gòu)也會(huì)經(jīng)常發(fā)生變化，尤其是涉及到復(fù)雜外形物體與自由面的相互作用時(shí)，這類問題會(huì)變得更加復(fù)雜.因?yàn)楹茈y得到這類問題的解析解，所以數(shù)值模擬成為求解含自由面流動(dòng)問題的重要手段，自由面流動(dòng)的數(shù)值模擬也日益成為計(jì)算流體力學(xué)領(lǐng)域的熱門研究問題.

本文采用LS(Level Set)方法對(duì)自由面流動(dòng)問題進(jìn)行數(shù)值分析.該方法最早由Osher等[1]提出來，后被廣泛應(yīng)用于各種自由面問題的數(shù)值模擬.

傳統(tǒng)的LS方法將LS方程與流體力學(xué)方程分開求解，即在得到流場信息后再推進(jìn)LS方程.這種求解方法不能實(shí)現(xiàn)流場與自由面的同步演化，自由面的演化過程存在一定的滯后現(xiàn)象[2-6].另外，廣泛應(yīng)用于自由面流動(dòng)分析的數(shù)值方法大多只有2階精度[2，4]，在某些分辨率要求較高的問題中，2階精度的數(shù)值方法對(duì)某些細(xì)微的流動(dòng)結(jié)構(gòu)捕捉不到.

基于此，本文將不可壓流體力學(xué)方程與LS方程耦合求解，通過引入偽時(shí)間導(dǎo)數(shù)項(xiàng)，得到關(guān)于偽時(shí)間的雙曲系統(tǒng)方程，并采用人工壓縮方法進(jìn)行求解.為了提高數(shù)值方法的分辨率，本文在有限體積框架下發(fā)展了一種基于HLLC(Harten-Lax-van Leer-Contact)近似Riemann解的高階WENO格式.

1控制方程

在LS方法中，用一個(gè)函數(shù)Ф(x，t)的零值面來描述自由面.該方法的優(yōu)點(diǎn)是可以精確地捕捉到自由面拓?fù)浣Y(jié)構(gòu)的變化.以二相流為例，假設(shè)物理區(qū)域Ω1中為流體A，Ω2中為流體B，Γ(t)為兩種流體的交界面.假設(shè)t=0，用d(x，Γ(0))表示計(jì)算點(diǎn)x到自由面Γ(0)的距離，定義LS函數(shù)為計(jì)算點(diǎn)到自由面的符號(hào)距離函數(shù)：

(1)

該函數(shù)要能夠在任意時(shí)間描述自由面，即需滿足

(2)

即

φt+V·φ=0.

(3)

上式中，V為流體速度.

對(duì)于不可壓縮流體，連續(xù)性方程為

·V=0.

(4)

在上式中引入壓強(qiáng)關(guān)于偽時(shí)間的導(dǎo)數(shù)項(xiàng)，可以得到

(5)

其中：β為人工壓縮系數(shù)；p為壓強(qiáng)；τ為偽時(shí)間.如果壓強(qiáng)在偽時(shí)間上收斂，那么連續(xù)性方程也成立.類似地，對(duì)于2維問題，在動(dòng)量方程中引入動(dòng)量關(guān)于偽時(shí)間的導(dǎo)數(shù)項(xiàng)，可以得到

(6)

上式中：g為重力加速度；P=p/β；u和v分別為速度V在x和y方向上的分量，黏性應(yīng)力項(xiàng)為

本文中μ為分子黏性系數(shù)，暫不考慮湍流流動(dòng).

將(3)、(5)和(6)式耦合，可以得到如下的系統(tǒng)控制方程：

(7)

其中：

2數(shù)值方法

對(duì)任意控制體將(7)式進(jìn)行積分，可以得到積分形式的控制方程：

(8)

上式中： Ω代表控制體單元，?Ω代表其單元邊界，S代表面積矢量，無黏通量H和黏性通量Hv為

H=Ei+Fj，Hv=Evi+Fvj.

其中：

i=(1，0)，j=(0，1).

用Wi，j和Gi，j表示單元上的體積平均，可以得到半離散形式的控制方程：

(9)

殘值R的定義為

上式中，無黏通量的計(jì)算采用了高階高斯積分，αk和Hk為高斯積分點(diǎn)處的權(quán)重和函數(shù)值，N為高斯點(diǎn)個(gè)數(shù)，本文中N=2.黏性通量的離散采用中心平均格式.

2.1HLLC近似Riemann解

下面將詳述如何利用HLLC近似Riemann解來給出單元界面上的無黏通量.對(duì)一個(gè)法向?yàn)閚=(nx，ny)的控制面，無黏通量為

其中q=unx+vny.該通量的Jacobian矩陣為

該Jacobian矩陣的特征值為

λ1=q-c，λ2=q，λ3=q，λ4=q+c，

(10)

其中

HLLC近似Riemann解是Toro等[7]在HLL近似Riemann解[8]的基礎(chǔ)上發(fā)展出來的一種新的近似求解Riemann問題的方法.在該方法中，Riemann問題的解用圖1(見第324頁)所示的三波模型進(jìn)行近似.相比于HLL近似Riemann解，HLLC近似Riemann解包含了接觸間斷，適合于模擬自由面流動(dòng)問題.

下面給出方程(7)所對(duì)應(yīng)的HLLC近似Riemann解.WL和WR代表單元界面左右兩側(cè)的重構(gòu)守恒變量.對(duì)于接觸間斷有

(11)

對(duì)于左右兩邊的非線性波，利用Rankine-Hugoniot條件有

(12)

上式中左右波速的估計(jì)采用如下方法：

上式中，上標(biāo)～代表相應(yīng)物理量的Roe平均值.將(12)式展開成分量形式，并利用接觸間斷關(guān)系式(11)，可以求出接觸間斷左右狀態(tài)的值：

(13)

其中K=L，R.最后可以得到如下的HLLC通量：

(14)

其中

2.2高階WENO重構(gòu)

單元交界面兩側(cè)高斯點(diǎn)處的守恒變量WL和WR通過高階WENO格式重構(gòu)得到.重構(gòu)過程是針對(duì)原始變量進(jìn)行的，得到原始變量的重構(gòu)值之后再轉(zhuǎn)換為守恒變量值.原始變量的定義為

下面給出重構(gòu)的具體過程.

如圖2所示，若已知原始變量中任一分量q在單元(m，n)上的平均值

(15)

其中m=i-2，i-1，i，i+1，i+2；n=j-2，j-1，j，j+1，j+2.單元(m，n)被4個(gè)單元界面包圍，這4個(gè)界面分別記為(m+1/2，n)，(m-1/2，n)，(m，n+1/2)和(m，n-1/2).

基于該單元平均值，首先在ξ方向上做重構(gòu)得到η方向上的面平均值，在圖2中用粗實(shí)線表示出單元界面(i+1/2，n)：

(16)

基于面平均值，在η方向上再進(jìn)行重構(gòu)，得到高斯積分點(diǎn)處左右兩側(cè)的原始變量值，在圖2中以實(shí)心圓點(diǎn)表示：

(17)

上述重構(gòu)過程均采用5階WENO重構(gòu)，具體方法可以參考文獻(xiàn)[9-11].

2.3雙時(shí)間推進(jìn)方法

時(shí)間推進(jìn)采用雙時(shí)間方法.首先忽略掉方程(9)中的偽時(shí)間導(dǎo)數(shù)項(xiàng)，并將物理時(shí)間導(dǎo)數(shù)采用二階向后差分近似：

(18)

上式中，上標(biāo)n代表物理時(shí)間站位，Wn和Wn-1為已知量，Wn+1為未知量.將Wn+1重新記為W*，并引入未知量的偽時(shí)間導(dǎo)數(shù)項(xiàng)，則有

(19)

(19)式在偽時(shí)間方向上的定常解W*即為物理時(shí)間上的時(shí)間精確解Wn+1.偽時(shí)間導(dǎo)數(shù)采用1階向后差分近似：

(20)

上式中，上標(biāo)m代表偽時(shí)間站位.當(dāng)m→∞時(shí)，(W*)m+1→Wn+1.方程(20)的求解采用高效的隱式LU-SSOR(Lower-Upper-Symmetric Successive Overrelaxation)方法，具體過程可以參考文獻(xiàn)[12].

3算例與分析

3.1Rayleigh-Taylor不穩(wěn)定問題

本算例計(jì)算一個(gè)空氣/氦氣的Rayleigh-Taylor不穩(wěn)定性問題，用來驗(yàn)證本文方法捕捉自由面形狀變化的能力.計(jì)算域?yàn)?.01m×0.04m，初始時(shí)交界面為波幅為5cm的余弦波.空氣與氦氣的密度分別為1.225kg/m3和0.1694kg/m3，黏性系數(shù)為1.77625×10-5kg·m-1·s-1和1.941×10-5kg·m-1·s-1.初始時(shí)流體速度均為0.計(jì)算網(wǎng)格大小為64×256.本算例不考慮表面張力的影響.

圖3(見第326頁)給出了6個(gè)時(shí)間(t=0.000，0.047，0.066，0.088，0.118，0.250s)自由面的本文計(jì)算結(jié)果，并和Puckett等[13]的計(jì)算結(jié)果作了對(duì)比，雖然本文采用的網(wǎng)格比文獻(xiàn)[13]中采用的網(wǎng)格要稀疏，但兩者得到的自由面形狀的發(fā)展趨勢基本一致.在重力場的作用下，密度大的空氣向下運(yùn)動(dòng)，氦氣則向上運(yùn)動(dòng)，自由面的擾動(dòng)幅度逐漸增大；在t=0.088s時(shí)，自由面的下部形成類似蘑菇云的形狀；在t=0.118s時(shí)，蘑菇云形狀被拉長，此時(shí)流場的重心要比初始時(shí)的重心低，重力勢能被逐漸釋放出來，流場的速度逐漸變大，重力勢能轉(zhuǎn)變?yōu)榱鲌龅膭?dòng)能.用本文的方法預(yù)測出在t=0.250s時(shí)空氣已經(jīng)達(dá)到容器的底部，自由面拓?fù)浣Y(jié)構(gòu)發(fā)生變化，出現(xiàn)了翻轉(zhuǎn)和破碎.

3.2三相流中的液泡上升問題

本算例模擬一液泡在部分盛水容器中的上升過程，驗(yàn)證本文方法處理多相流的能力.為了方便起見，把3種模型流體介質(zhì)記為水、油和空氣.3種流體的密度比取為100∶50∶1，黏性系數(shù)比取為100∶50∶1.水的密度取為1.0×103kg/m3，黏性系數(shù)取為1.7921×10-3kg·m-1·s-1.雷諾數(shù)定義為

本文中取Re=200，不考慮表面張力作用.計(jì)算幾何模型如圖4所示，計(jì)算域大小為3L×3.5L，初始水深2.5L，空氣厚度為1L，直徑為D0=2R=L的圓形油泡浸沒在水中，圓心到容器底部距離為1.5L.計(jì)算網(wǎng)格為60×70的笛卡爾網(wǎng)格.初始速度場為0，泡內(nèi)壓強(qiáng)均取液泡中心處流體靜壓強(qiáng)：

p(x，y)=(ρoilR+0.5ρwaterL+ρa(bǔ)irL)g.

一個(gè)Level Set函數(shù)Ф只能標(biāo)記2種流體，而本算例包含3種流體和2個(gè)交界面，需要再引入一個(gè)Level Set函數(shù)ψ.定義Ф為到水/油交界面的符號(hào)距離函數(shù)，ψ為到水/氣交界面的符號(hào)距離函數(shù)，Ф和ψ的符號(hào)設(shè)定如圖4所示.

圖5給出了本文計(jì)算得到的自由面形狀與流場，每幅圖所對(duì)應(yīng)的無量綱時(shí)間t分別為0.5、1.0、1.5、2.0、2.5、3.0、3.5、4.0，從圖中可以看出，在浮力作用下油泡逐漸上升.在t=1.5時(shí)，油泡底部變平；在t=2.0時(shí)，油泡變?yōu)槟I形；隨著油泡進(jìn)一步接近水/空氣界面，油泡對(duì)該自由面形成擾動(dòng)，可以看到水/空氣界面發(fā)生明顯變化.本文計(jì)算結(jié)果與Pan等[14]的計(jì)算結(jié)果相比，兩者的自由面位置基本一致.雖然本文采用的網(wǎng)格比文獻(xiàn)[14]采用的網(wǎng)格要稀疏，但是在t=3.0時(shí)間后，本文結(jié)果捕捉到了更細(xì)節(jié)的自由面結(jié)構(gòu).

4結(jié)論

本文發(fā)展了一種采用Level Set方程與流體力學(xué)控制方程相耦合的求解自由面流動(dòng)的數(shù)值方法.利用人工壓縮方法，將可壓縮流動(dòng)的計(jì)算方法引入到不可壓自由面流動(dòng)的數(shù)值模擬.為了提高自由面流動(dòng)計(jì)算的分辨率，本文提出了一種有效的WENO -HLLC格式，構(gòu)造了基于HLLC近似Riemann解的通量計(jì)算格式.對(duì)2個(gè)典型自由面流動(dòng)問題進(jìn)行數(shù)值模擬，計(jì)算結(jié)果與相關(guān)參考文獻(xiàn)的計(jì)算結(jié)果吻合很好，說明了本文算法的有效性.

參考文獻(xiàn)：

[1]OSHER S， SETHIAN J A. Fronts propagating with curvature-dependent speed： Algorithms based on Hamilton-Jacobi formulations [J].JournalofComputationalPhysics， 1988，79(1)： 12-49.

[2]PRICE W G， CHEN Y G. A simulation of free surface waves for incompressible two-phase flows using a curvilinear level set formulation [J].InternationalJournalforNumericalMethodsinFluids， 2006，51(3)： 305-330.

[3]KURIOKA S， DOWLING R. Numerical simulation of free surface flows with the level set method using an extremely high-order accuracy WENO advection scheme [J].InternationalJournalofComputationalFluidDynamics， 2009，23(3)： 233-243.

[4]SUSSMAN M， SMEREKA P， OSHER S. A level set approach for computing solutions to incompressible two-phase flow [J].JournalofComputationalPhysics， 1994，114(1)： 146-159.

[5]GROOSS J， HESTHAVEN J S. A level set discontinuous Galerkin method for free surface flows [J].ComputerMethodsinAppliedMechanicsandEngineering， 2006，195(25)： 3406-3429.

[6]袁德奎，陶建華.用Level-Set方法求解具有自由面的流動(dòng)問題 [J].力學(xué)學(xué)報(bào)，2000，32(3)： 264-271.

[7]TORO E F， SPRUCE M， SPEARES W. Restoration of the contact surface in the HLL-Riemann solver [J].ShockWaves， 1994，4(1)： 25-34.

[8]HARTEN A， LAX P D， LEER B. On upstream differencing and Godunov-type schemes for hyperbolic conservation laws [J].SIAMReview， 1983，25(1)： 35-61.

[9]SHU C W. High order ENO and WENO schemes for computational fluid dynamics[M]//High-order methods for computational physics. Berlin， Heidelberg： Springer， 1999： 439-582.

[10]BALSARA D S， SHU C W. Monotonicity preserving weighted essentially non-oscillatory schemes with increasingly high order of accuracy [J].JournalofComputationalPhysics， 2000，160(2)： 405-452.

[11]TITAREV V A， TORO E F. Finite-volume WENO schemes for three-dimensional conservation laws [J].JournalofComputationalPhysics， 2004，201(1)： 238-260.

[12]YANG A， CHEN S， YANG L，etal. An upwind finite volume method for incompressible inviscid free surface flows [J].Computers&Fluids， 2014，101： 170-182.

[13]PUCKETT E G， ALMGREN A S， BELL J B，etal. A high-order projection method for tracking fluid interfaces in variable density incompressible flows [J].JournalofComputationalPhysics， 1997，130(2)： 269-282.

[14]PAN D， CHANG C H. The capturing of free surfaces in incompressible multi-fluid flows [J].InternationalJournalforNumericalMethodsinFluids， 2000，33(2)： 203-222.

文章編號(hào)：0427-7104(2016)03-0321-08

收稿日期：2016-01-10

基金項(xiàng)目：國家自然科學(xué)基金(11172070)

作者簡介：李泠泉(1989—)，女，碩士研究生；楊愛明(1971—)，男，博士，副教授，通訊聯(lián)系人，E-mail： amyang@fudan.edu.cn.

中圖分類號(hào)：O 359.1

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

Numerical Simulation of Free Surface Flow Based on HLLC Riemann Solver

LI Lingquan， YANG Aiming

(DepartmentofAeronauticsandAstronautics，F(xiàn)udanUniversity，Shanghai200433，China)

Abstract：A new numerical method of simulating free surface flow problems is developed. In this method， Level Set function is used to capture free surface. Hyperbolic equations are formed by coupling the Level Set function and the hydrodynamic equations. Artificial compressibility method is used to solve the equations. In order to improve the resolution of calculation， a high order WENO scheme based on finite volume method is developed. The flux function based on HLLC Riemann solver is constructed. Several classic free surface flow problems are simulated to prove that the present method is effective.

Keywords：free surface flow； Level Set； WENO scheme； HLLC Riemann solver； artificial compressibility