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

李泠泉，楊愛明

(復旦大學 航空航天系，上海 200433)

?

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

李泠泉，楊愛明

(復旦大學 航空航天系，上海 200433)

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

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

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

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

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

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

1控制方程

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

(1)

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

(2)

即

φt+V·φ=0.

(3)

上式中，V為流體速度.

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

·V=0.

(4)

在上式中引入壓強關(guān)于偽時間的導數(shù)項，可以得到

(5)

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

(6)

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

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

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

(7)

其中：

2數(shù)值方法

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

(8)

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

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

其中：

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

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

(9)

殘值R的定義為

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

2.1HLLC近似Riemann解

下面將詳述如何利用HLLC近似Riemann解來給出單元界面上的無黏通量.對一個法向為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頁)所示的三波模型進行近似.相比于HLL近似Riemann解，HLLC近似Riemann解包含了接觸間斷，適合于模擬自由面流動問題.

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

(11)

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

(12)

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

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

(13)

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

(14)

其中

2.2高階WENO重構(gòu)

單元交界面兩側(cè)高斯點處的守恒變量WL和WR通過高階WENO格式重構(gòu)得到.重構(gòu)過程是針對原始變量進行的，得到原始變量的重構(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個單元界面包圍，這4個界面分別記為(m+1/2，n)，(m-1/2，n)，(m，n+1/2)和(m，n-1/2).

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

(16)

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

(17)

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

2.3雙時間推進方法

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

(18)

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

(19)

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

(20)

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

3算例與分析

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

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

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

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

本算例模擬一液泡在部分盛水容器中的上升過程，驗證本文方法處理多相流的能力.為了方便起見，把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，不考慮表面張力作用.計算幾何模型如圖4所示，計算域大小為3L×3.5L，初始水深2.5L，空氣厚度為1L，直徑為D0=2R=L的圓形油泡浸沒在水中，圓心到容器底部距離為1.5L.計算網(wǎng)格為60×70的笛卡爾網(wǎng)格.初始速度場為0，泡內(nèi)壓強均取液泡中心處流體靜壓強：

p(x，y)=(ρoilR+0.5ρwaterL+ρairL)g.

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

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

4結(jié)論

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

參考文獻：

[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方法求解具有自由面的流動問題 [J].力學學報，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.

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

收稿日期：2016-01-10

基金項目：國家自然科學基金(11172070)

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

中圖分類號：O 359.1

文獻標志碼：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