劉 娜，陳藝冰

（1.北京應(yīng)用物理與計算數(shù)學(xué)研究所，北京100088；2.中國工程物理研究院高性能數(shù)值模擬軟件中心，北京100088；3.北京應(yīng)用物理與計算數(shù)學(xué)研究所計算物理實驗室，北京100088）

多介質(zhì)流體力學(xué)計算的譜體積方法＊

劉 娜1，2，3，陳藝冰1，3

（1.北京應(yīng)用物理與計算數(shù)學(xué)研究所，北京100088；2.中國工程物理研究院高性能數(shù)值模擬軟件中心，北京100088；3.北京應(yīng)用物理與計算數(shù)學(xué)研究所計算物理實驗室，北京100088）

針對高維及多物理耦合計算耗費大等困難，設(shè)計適合多介質(zhì)流動模擬的模板緊致、易于并行、高階精度、計算耗費小的譜體積方法。該方法是求解雙曲型守恒率譜體積方法的直接推廣，針對多介質(zhì)流動物質(zhì)界面捕捉的困難，利用擬守恒格式的思想避免物質(zhì)界面處的非物理振蕩。數(shù)值模擬結(jié)果表明，本方法具有高階精度、高分辨率，且節(jié)約計算量，并且可以有效避免物質(zhì)界面處非物理振蕩。

流體力學(xué)；譜體積方法；高精度；剛性氣體；多介質(zhì)；模板緊致

多介質(zhì)流體力學(xué)的模型可以用擴(kuò)展歐拉方程組描述。早期，使用守恒格式直接求解，但在接觸間斷（物質(zhì)界面）附近，會產(chǎn)生非物理的壓力和速度振蕩。隨后，很多學(xué)者對此進(jìn)行分析，并提出各種補(bǔ)救措施［1－2］，其中R.Abgrall等提出的擬守恒型格式，由于其簡單、健壯得到了廣泛應(yīng)用。此后，K.M.Shyue［3－4］、Y.B.Chen等［5］進(jìn)一步發(fā)展了這個方法。這些工作通常采用MUSCL限制器，精度不超過二階。

當(dāng)前，實際工程問題中對精細(xì)化設(shè)計提出越來越高的要求，不少更高階（高于2階）精度的格式被應(yīng)用于求解多介質(zhì)問題。其中，E.Johnsen等［6］將WENO格式推廣于求解擴(kuò)展歐拉方程組，取得成功。他的算法極大提高了對物質(zhì)界面的分辨率，但WENO格式進(jìn)行重構(gòu)時，需要利用目標(biāo)單元及其周圍單元的信息，隨著格式階數(shù)的增大，需要不斷擴(kuò)展模板，在高維情形復(fù)雜度將急劇增加。J.Zhu等［7］采用DG格式來求解這個問題。與WENO相比，DG格式僅采用局部模板即可獲得高精度。實際上，DG格式是近年來廣受關(guān)注的高精度緊致格式的典型代表之一。除了DG方法之外，譜體積（SV）格式［8－9］也是一類成功的高精度緊致格式。它是譜方法和有限體積方法的結(jié)合體，其主要思想是根據(jù)細(xì)分的控制體的單元均值，重構(gòu)出譜體積元中的近似解。本質(zhì)上，DG方法和SV方法都是采用目標(biāo)單元內(nèi)部信息來獲得網(wǎng)格單元內(nèi)解的逼近，因此模板都較為緊致，且在并行計算時不需要通信周圍單元的信息，是一種高并行度的格式。與DG方法相比，SV方法的一個優(yōu)點在于節(jié)省了表面積分和體積分的計算，并且根據(jù)Z.J.Wang［8］的分析，SV方法對CFL數(shù)的限制也更低。

本文中，針對高維及多物理耦合計算耗費大等困難，設(shè)計適合多介質(zhì)流動模擬的譜體積方法。

1 控制方程

非定常可壓縮流體可以用歐拉方程來描述，以一維為例可以記為：W／t＋F（W）／x＝0，其中W＝（ρ，ρu，ρE）T，F(xiàn)（W）＝（ρu，ρu2＋p，（ρE＋p）u）T。剛性氣體狀態(tài)方程為p＝（γ－1）ρe－γp，其中γ為比熱比，p為參考壓力。由于流體流動時狀態(tài)量γ隨著流動保持不變，因此滿足［10］：γ／t＋uγ／x＝0。定義質(zhì)量比分Y＝ρ1／（ρ1＋ρ2），建立γ與Y的關(guān)系：則質(zhì)量比分同樣滿足：

式中：W＝（ρ，ρu，ρE，ρY）T，F(xiàn)（W）＝（ρu，ρu2＋p，（ρE＋p）u，ρuY）T。

2 守恒譜體積方法

直接將譜體積方法推廣應(yīng)用于多介質(zhì)Euler方程組（2）。下面給出格式設(shè)計的基本框架：將計算區(qū)域Ω劃分為N個子區(qū)間｛Si＝［xi－1／2，xi＋1／2］Ω，i∈Z｝，每個Si稱為譜體積元，hi＝xi＋1／2－xi－1／2表示其區(qū)間長度。假設(shè)要設(shè)計一個k階精度格式，則將每個譜體積元Si進(jìn)一步分化成k個控制體，即Ci，j＝（xi，j－1／2，xi，j＋1／2），j＝1，…，k。在每個控制體上定義控制體單元平均值tn）dx，j＝1，…，k。將控制方程（2）在Ci，j上積分，可以得到如下半離散格式：

要完成k階譜體積格式，還需要3個步驟。

步驟1：重構(gòu)。利用給定的tn時刻的控制體單元平均值進(jìn)行重構(gòu)。在每個譜體積單元Si內(nèi)部，由k個控制體的單元平均值，重構(gòu)出在譜體積元Si內(nèi)的一個至多k－1次多項式Wni（x）。這樣得到的重構(gòu)多項式可能會引起數(shù)值振蕩，為了避免高階重構(gòu)帶來的數(shù)值振蕩，對重構(gòu)多項式利用基于原始變量的修正TVB minmod斜率限制器的思想進(jìn)行限制。詳細(xì)過程可以參見文獻(xiàn)［8］。經(jīng)過限制，在整個計算區(qū)域的重構(gòu)函數(shù)具有如下性質(zhì)：（1）在譜體積單元內(nèi)部不需要限制的控制體界面處，重構(gòu)多項式函數(shù)是連續(xù)的；（2）在譜體積單元界面處或譜體積單元內(nèi)部限制器起作用的控制體單元界面處，重構(gòu)多項式函數(shù)是間斷的。

步驟2：計算數(shù)值通量。根據(jù)重構(gòu)多項式的性質(zhì)，在譜體積單元內(nèi)部狀態(tài)連續(xù)的控制體單元界面處，直接使用連續(xù)的數(shù)值通量，而在譜體積單元界面處以及譜體積單元內(nèi)部經(jīng)過限制后的控制體單元界面處，由于兩側(cè)的狀態(tài)值不同，使用間斷數(shù)值通量。常用的數(shù)值通量計算方法，如Lax－Friedrichs、HLLC、Roe、BGK通量等都可以使用，本文中選擇最簡便的Lax－Friedrichs數(shù)值通量。

步驟3：時間方向離散。時間方向可以采用4階Runge－Kutta方法離散。

3 穩(wěn)定性分析

根據(jù)第2節(jié)譜體積方法在多介質(zhì)守恒擴(kuò)展歐拉方程組上的直接推廣計算多介質(zhì)問題，會產(chǎn)生物質(zhì)界面處壓力、速度的非物理振蕩。在第5.2節(jié)將顯示一維多介質(zhì)運(yùn)動界面問題的計算結(jié)果。R.Abgrall等［1－2］指出，這些非物理振蕩是由于對組分方程的守恒離散造成的。下面以Lax－Friedrichs通量為例進(jìn)行分析，為了簡便起見，時間方向考慮Euler方法離散，則質(zhì)量、動量、能量方程的離散可以表示為：

若n時刻速度u＝u0和壓力p＝p0為常數(shù)，則根據(jù)（5）結(jié)合（4），可得：，即在n＋1時刻速度可以保持常數(shù)：進(jìn)一步，將狀態(tài)方程代入（6），則n＋1時刻滿足

4 擬守恒譜體積方法

根據(jù)第3節(jié)的穩(wěn)定性分析，對第2節(jié)中的數(shù)值格式做一些修正。將質(zhì)量比分方程（1）寫成如下形式Y(jié)／t＋f／x－Yg／x＝0，其中f＝uY，g＝u。質(zhì)量比分方程的半離散格式修正為：

在物質(zhì)界面處，速度u保持常數(shù)，則離散格式（8）滿足壓力無振蕩的必要條件（7），且離散格式與質(zhì)量比分方程（1）相容。

5 數(shù)值模擬

數(shù)值模擬用以檢驗格式的精度和性能，計算中一維算例CFL數(shù)取1.0，二維取0.5。

5.1 精度測試

在區(qū)域［0，2］上有兩種介質(zhì)組成的正弦波：

周期邊界條件。表1中給出t＝1.0時各階格式精度測試結(jié)果，結(jié)果顯示2～5階格式都可以達(dá)到相應(yīng)收斂率。

表1 格式的數(shù)值精度Table 1 Numerical accuracy of present schemes

5.2 一維多介質(zhì)運(yùn)動界面問題

在區(qū)域［0，1］上有兩種介質(zhì)流體，以相同的速度勻速向右移動，界面勻速向右移動，初值條件為：

圖1顯示了采用守恒譜體積格式100個譜體積單元在t＝0.1時的壓力和速度，可以看出，在物質(zhì)界面處出現(xiàn)速度與壓力的數(shù)值振蕩，與第3節(jié)中的分析結(jié)果一致。為方便看清圖像，圖中只顯示了奇數(shù)階結(jié)果。圖2顯示了格式修正后的計算結(jié)果，結(jié)果顯示，此時在物質(zhì)界面處不會產(chǎn)生速度和壓力的數(shù)值振蕩，物質(zhì)界面附近密度和比熱比的局部放大圖顯示高階格式具有更高的分辨率。

圖1 守恒譜體積格式的一維多介質(zhì)運(yùn)動界面問題Fig.1 One－dimensional moving interface problem for conservative spectral volume scheme

圖2 擬守恒譜體積格式的一維多介質(zhì)運(yùn)動界面問題Fig.2 One－dimensional moving interface problem for quasi－conservative spectral volume scheme

5.3 高壓力比氣液激波管問題

在區(qū)域［－0.2，1］上求解初值問題：

圖3顯式了t＝0.000 2時刻使用720個譜體積單元1、3、5階擬守恒譜體積格式計算的結(jié)果，結(jié)果顯示，此時在物質(zhì)界面處不會產(chǎn)生速度和壓力的數(shù)值振蕩。圖4顯示了密度及其在物質(zhì)界面和激波附近的局部放大圖，結(jié)果表明高階格式具有更高分辨率。

圖3 擬守恒譜體積格式的高壓力比氣液激波管問題Fig.3 High pressure ratio gas－liquid shock tube problem for quasi－conservative spectral volume scheme

圖4 擬守恒譜體積格式的高壓力比氣液激波管問題Fig.4 High pressure ratio gas－liquid shock tube problem for quasi－conservative spectral volume scheme

5.4 三點問題

初值條件為：

初始空氣（輕介質(zhì)）位于水（重介質(zhì)）上方，兩種介質(zhì)流體保持靜止，靜止高密度高壓氣體位于其右側(cè)分別與低壓空氣和水產(chǎn)生作用，分別產(chǎn)生向左運(yùn)動的稀疏波、向右移動的激波，以及推動物質(zhì)界面向右移動。由于右側(cè)上方的空氣較輕，因此右側(cè)上方的激波和物質(zhì)界面向右移動的速度更快，于是在3種狀態(tài)的流體的交匯點產(chǎn)生漩渦。形成漩渦后，上層空氣與下層水的交界面發(fā)生彎曲，并與激波發(fā)生作用，形成馬赫桿。圖5顯示了t＝1.0時使用600×200個網(wǎng)格2～5階格式計算的結(jié)果，可以看到，各階格式都能計算出此時基本的波結(jié)構(gòu)，并且隨著格式精度的增加能夠捕捉到更精細(xì)的結(jié)構(gòu)。

圖5 擬守恒譜體積格式的三點問題密度等值線圖Fig.5 Density contours of triple point problem for quasi－conservative spectral volume scheme

6 結(jié) 論

以多介質(zhì)擬守恒格式為基礎(chǔ)，通過高精度譜體積重構(gòu)，構(gòu)造了一類求解多介質(zhì)問題的高階精度擬守恒譜體積格式，有效提高了求解多介質(zhì)問題的精度，并且繼承譜體積格式模板緊致、易于并行等優(yōu)點。數(shù)值結(jié)果表明，新格式求解多介質(zhì)問題具有高階精度，且模板緊致，并且不會在物質(zhì)界面處產(chǎn)生非物理振蕩。

參考文獻(xiàn)：

［1］Abgrall R.How to prevent pressure oscillations in multicomponent flow calculation：A quasi conservative approach［J］.Journal of Computational Physics，1996，125（1）：150－160.

［2］Abgrall R，Karni S.Computations of compressible multifluids［J］.Journal of Computational Physics，2001，169（2）：594－623.

［3］Shyue K M.An effcient shock－capturing algorithm for compressible multicomponent problems［J］.Journal of Computational Physics，1998，142（1）：208－242.

［4］Shyue K M.A fluid－mixture type algorithm for compressible multicomponent flow with Mie－Gruneisen equation of state［J］.Journal of Computational Physics，2001，171（2）：678－707.

［5］Chen Y B，Jiang S.A non－oscillatory kinetic scheme for multicomponent flows with the equation of state for a stiffned gas［J］.Journal of Computational Mathematics，2011，29（6）：661－683.

［6］Johnsen E，Colonius T.Implementation of WENO schemes in compressible multicomponent flow problems［J］.Journal of Computational Physics，2006，219（2）：715－732.

［7］Zhu J，Qiu J X，Liu T G，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（4）：554－580.

［8］Wang Z J.Spectral（finite）volume method for conservation laws on unstructured grids：Basic formulation［J］.Journal of Computational Physics，2002，178（1）：210－251.

［9］Wang Z J，Liu Y.Spectral（finite）volume method for conservation laws on unstructured gridsⅢ：One dimensional systems and partition optimization［J］.Journal of Scientific Computing，2004，20（1）：137－157.

［10］Karni S.Multicomponent flow calculations by a consistent primitive algorithm［J］.Journal of Computational Physics，1994，112（1）：31－43.

High order spectral volume method for multi－component flows

Liu Na1，2，3，Chen Yibing1，3
（1.Institute of Applied Physics and Computational Mathematics，Beijing100088，China；2.CAEP Software Center for High Performance Numerical Simulation，Beijing100088，China；3.Laboratory of Computational Physics，Institute of Applied Physics and Computational Mathematics，Beijing100088，China）

Numerical simulation of multi－material flows has been an important issues in CFD，and most CFD production codes used for multi－material flow simulation is of either first or second order accuracy，too inefficient and costly with its grid refinement for high accuracy required problems.In this paper，a high－order，efficient，compact method，called the spectral volume method，was developed for the simulation of the multi－material flow as an extension of the spectral volume method for the conservation laws.It has been pointed out that the conservative spectral volume method for the multi－material flow will cause oscillation，and the reason for this has been analyzed.So the idea of quasi－conservative scheme was borrowed to prevent the spurious oscillations in the vicinity of a material contact discontinuity.Several numerical experiments proved that there is no oscillation near the material interface and the result also demonstrates the accuracy，the efficiency and the high performance of the scheme for the multi－material flow simulation.

fluid mechanics；spectral volume method；high order；stiffened gas；multi－component flows；compact stencil

O357.4國標(biāo)學(xué)科代碼：1302511

A

10.11883／1001－1455（2017）01－0114－06

（責(zé)任編輯 丁 峰）

2015－05－11；

2015－09－20

國家自然科學(xué)基金項目（11101047，11501043，11671050，91430218，U1630247）；

中國工程物理研究院科學(xué)技術(shù)發(fā)展基金項目（2013A0202011）

劉 娜（1986— ），女，博士，助理研究員，liu＿na＠iapcm.ac.cn。