一維低密度低壓力多介質(zhì)流動(dòng)問(wèn)題的高精度數(shù)值方法

一維低密度低壓力多介質(zhì)流動(dòng)問(wèn)題的高精度數(shù)值方法

吳招生

(南京航空航天大學(xué) 理學(xué)院數(shù)學(xué)系，南京 210016)

摘要：在爆炸波等多介質(zhì)流體力學(xué)問(wèn)題的數(shù)值模擬中，會(huì)出現(xiàn)低密度低壓力區(qū)域，由于數(shù)值誤差的影響，導(dǎo)致在計(jì)算過(guò)程中該區(qū)域的密度和壓力為負(fù)值，使得計(jì)算難以繼續(xù).文獻(xiàn)[1]給出了求解理想氣體歐拉方程的正保護(hù)RKDG(Runge-Kutta Discontinuous Galerkin)方法.在此基礎(chǔ)上，給出了多介質(zhì)流體力學(xué)問(wèn)題數(shù)值模擬的正保護(hù)RKDG方法，結(jié)合RGFM(Real Ghost Fluid Method)界面處理方法，對(duì)一維多介質(zhì)可壓縮流體進(jìn)行了數(shù)值試驗(yàn)，結(jié)果表明該方法能有效地避免壓力為負(fù)的情況，并能準(zhǔn)確地捕捉各類(lèi)間斷.

關(guān)鍵詞：RKDG方法；低密度；低壓力；正保護(hù)限制器

中圖分類(lèi)號(hào)：O241.82文獻(xiàn)標(biāo)志碼：A

文章編號(hào)：1008-5564(2015)03-0026-05

收稿日期：2015-04-20

作者簡(jiǎn)介：劉小剛(1983—)，男，陜西延安人，西北大學(xué)現(xiàn)代學(xué)院助教，碩士，主要從事微分方程與動(dòng)力系統(tǒng)研究;

High-precision Numerical Method for Solving Problems ofOne-dimensional Multi-medium Flows with Low-density and Low-pressure

WU Zhao-sheng

(Department of Mathematics, School of Science, Nanjing University of

Aeronautics and Astronautics, Nanjing 210016, China)

Abstract：In the numerical simulation of problems of the multi-medium fluid mechanics of blast waves, the domain with low density and low pressure would arise, and density and pressure in this domain would be negative during the calculation because of numerical error, so the calculation is difficult to continue. The positivity-preserving RKDG (Discontinuous Galerkin Runge-Kutta) method was presented for solving the ideal gas Euler equation in reference [1]. Based on all of above, the positivity-preserving RKDG method of numerical simulation was presented for solving problems of the multi-medium fluid mechanics. Combining with the RGFM (Ghost Fluid Method Real) interface processing method, numerical experiment was conducted for one-dimensional compressible fluid. The results show that this method can effectively avoid the negative pressure and can accurately capture all kinds of discontinuities.

Key words：RKDG method; Low density; Low pressure; the positivity-preserving limiter

間斷有限元方法是1973年由Reed和Hill提出的，用來(lái)求解中子運(yùn)輸問(wèn)題[2].20世紀(jì)80年代后期，Cockburn和C.W.Shu等構(gòu)造了求解一維標(biāo)量守恒律方程[3]的高階RKDG有限元方法，并且將該方法成功地推廣到多維標(biāo)量守恒律[4]和多維守恒律方程組[5].RKDG方法的優(yōu)點(diǎn)是具有高精度且易于處理復(fù)雜邊界.此后，RKDG方法開(kāi)始廣泛應(yīng)用于流體力學(xué)問(wèn)題的數(shù)值模擬.

高精度守恒型差分格式[6]由于其能夠準(zhǔn)確捕捉激波及其它物理現(xiàn)象，在流體力學(xué)數(shù)值計(jì)算中得到了廣泛的應(yīng)用.物理上密度和壓力應(yīng)該為正，但在爆炸波和高速射流等問(wèn)題中通常會(huì)出現(xiàn)低密度和低壓力區(qū)域，在用高精度格式求解時(shí)，由于數(shù)值誤差的影響，導(dǎo)致密度或壓力出現(xiàn)負(fù)值，使得計(jì)算難以繼續(xù)進(jìn)行.文獻(xiàn)[1]和[7]分別針對(duì)高精度RKDG格式和WENO格式構(gòu)造了正保護(hù)格式，并推廣到含有源項(xiàng)的歐拉方程組[8].本文在此工作的基礎(chǔ)上，開(kāi)展多介質(zhì)流動(dòng)問(wèn)題高精度數(shù)值計(jì)算格式正保護(hù)性質(zhì)的研究.

對(duì)于多介質(zhì)流動(dòng)問(wèn)題，界面兩邊是不同性質(zhì)的流體，它們的狀態(tài)方程可能會(huì)有很大差異，流體的密度會(huì)出現(xiàn)大梯度的變化，如果還采用單介質(zhì)流動(dòng)問(wèn)題的數(shù)值模擬方法，則會(huì)出現(xiàn)數(shù)值不穩(wěn)定的情況，在界面處會(huì)產(chǎn)生非物理震蕩，甚至使得計(jì)算難以進(jìn)行.因此，在多介質(zhì)流動(dòng)問(wèn)題的數(shù)值計(jì)算中，界面處理方法的研究是關(guān)鍵.Fedkiw等[9]提出的Ghost方法，把多介質(zhì)流體問(wèn)題轉(zhuǎn)化為多個(gè)單介質(zhì)問(wèn)題求解.在此基礎(chǔ)上提出的基于求解Riemann問(wèn)題的界面處理方法(簡(jiǎn)稱(chēng)RGFM方法[10])，給出了更加精確的界面邊界條件，能準(zhǔn)確地捕捉界面位置.本文利用RGFM方法定義界面邊界條件，將多介質(zhì)流動(dòng)問(wèn)題轉(zhuǎn)化為單介質(zhì)流動(dòng)問(wèn)題，利用所給出的具有正保護(hù)性質(zhì)的高精度RKDG格式進(jìn)行計(jì)算.

本文在不改變文獻(xiàn)[1]中正保護(hù)限制器性質(zhì)的情況下，利用壓力函數(shù)是凹函數(shù)的性質(zhì)，進(jìn)行線性插值求解，避免了求解一元二次方程，節(jié)省了計(jì)算量，并將此正保護(hù)限制器推廣應(yīng)用于水的狀態(tài)方程(Tait方程).本文采用三階具有正保護(hù)性質(zhì)的RKDG格式結(jié)合RGFM界面處理方法數(shù)值模擬一維多介質(zhì)可壓縮Euler方程組，對(duì)具有大密度比的一維氣-氣，氣-液兩種介質(zhì)可壓縮流體問(wèn)題進(jìn)行了數(shù)值實(shí)驗(yàn).

1方程

1.1控制方程

考慮一維無(wú)粘可壓縮流體力學(xué)方程組：

(1)

1.2狀態(tài)方程

理想氣體的狀態(tài)方程：

p=(γ-1)ρe

(2)

水的狀態(tài)方程：

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

(3)

其中γ和PC均為正常數(shù)，將在具體算例中給出.

2數(shù)值方法

2.1RKDG方法

考慮Euler方程組(1)的初值問(wèn)題：

(4)

其中[a,b]為計(jì)算區(qū)域.

(5)

(6)

(7)

(8)

用Uh(x,t)代替U(x,t)，得到：

(9)

(9)式可以寫(xiě)成：

(10)

方程組(10)的時(shí)間離散采用三階TVB Runge-Kutta方法[6].

2.2RGFM界面處理方法

假設(shè)流體Ⅰ與流體Ⅱ的界面位于網(wǎng)格格點(diǎn)i和i+1之間，對(duì)流體Ⅰ而言，DG格式邊界計(jì)算時(shí)需要網(wǎng)格單元i+1的狀態(tài)值，這個(gè)單元的狀態(tài)值由在界面處構(gòu)造的Riemann問(wèn)題的解重新賦值.Riemann問(wèn)題的初始狀態(tài)取值如下：

(11)

圖1　流體I界面邊界條件的定義

分別定義了流體Ⅰ與流體Ⅱ的界面邊界條件后，再利用所得到的高精度RKDG格式分別求解流體Ⅰ與流體Ⅱ，得到下一個(gè)時(shí)間步的數(shù)值解.

3正保護(hù)限制器

正保護(hù)限制器的具體實(shí)施分兩大步驟：第一步，對(duì)密度修正，保證每個(gè)網(wǎng)格單元上的3個(gè)Gauss-Lobatto點(diǎn)的密度近似函數(shù)的值大于零；第二步，當(dāng)流體是理想氣體時(shí)，對(duì)壓力修正.當(dāng)流體是液體時(shí)，對(duì)內(nèi)能修正.保證每個(gè)網(wǎng)格單元上的3個(gè)Gauss-Lobatto點(diǎn)的壓力近似函數(shù)或者內(nèi)能的近似函數(shù)的值大于零.限制器放在時(shí)間t推進(jìn)到t+Δt的每個(gè)RK子步后進(jìn)行限制，不影響原計(jì)算格式的精度和守恒性.

下面給出不同流體情形下的正保護(hù)限制器的具體實(shí)施和算法.

3.1流體為理想氣體

對(duì)于理想氣體而言，密度和壓力在物理意義下始終大于零.定義集合G：

(12)

可以證明集合G是凸集合，壓力函數(shù)p是凹函數(shù).

3.2流體為液體

對(duì)于液體而言，密度和內(nèi)能在物理意義下始終為正.定義集合G：

(13)

(1)密度修正：

(2)內(nèi)能修正：

4數(shù)值實(shí)驗(yàn)

4.1氣-氦激波管問(wèn)題

這是一個(gè)一維氣-氣多介質(zhì)問(wèn)題，初始條件如下：

由于數(shù)值誤差的影響，不使用正保護(hù)限制器，在低壓區(qū)壓力出現(xiàn)負(fù)值，程序崩潰.取網(wǎng)格數(shù)800，TVB參數(shù)取為(M1,M2,M3)=(1010,10,1010).隨著時(shí)間推進(jìn)，初始間斷分解為向右的強(qiáng)激波和向左的稀疏波，圖2～圖4分別給出了t=6時(shí)刻的密度、壓力和速度分布圖.可見(jiàn)，本文所給出的算法能夠保證低壓區(qū)的數(shù)值計(jì)算，保證壓力為正.

圖2　氣-氦激波管問(wèn)題 密度

圖3　氣-氦激波管問(wèn)題 壓力

圖4　氣-氦激波管問(wèn)題 速度

4.2氣-液激波管問(wèn)題

這是一個(gè)一維氣-液多介質(zhì)問(wèn)題，初始條件如下：

由于數(shù)值誤差的影響，不使用正保護(hù)限制器，水中內(nèi)能出現(xiàn)負(fù)值，程序崩潰.取網(wǎng)格數(shù)200，TVB參數(shù)取為(M1,M2,M3)=(50,50,50).隨著時(shí)間推進(jìn)，初始間斷分解為向左的強(qiáng)激波和向右的稀疏波，圖5～圖7分別給出了t=0.000 24時(shí)刻的密度、速度和內(nèi)能分布圖.可見(jiàn)，本文所給出的算法能夠保證內(nèi)能為正.

圖5　氣-液激波管問(wèn)題 密度

圖6　氣-液激波管問(wèn)題 速度

圖7　氣-液激波管問(wèn)題 內(nèi)能

5總結(jié)與展望

本文對(duì)原有正保護(hù)限制器進(jìn)行修正，避免了求解一元二次方程，節(jié)省了計(jì)算量，并將此正保護(hù)限制器推廣到水的狀態(tài)方程(Tait方程)，結(jié)合RGFM界面處理方法，應(yīng)用于一維多介質(zhì)流體力學(xué)問(wèn)題的數(shù)值模擬.數(shù)值試驗(yàn)表明該方法能準(zhǔn)確捕捉界面和其它間斷，體現(xiàn)了算法的保正性.

[參考文獻(xiàn)]

[1]ZHANGX,SHUCW.OnpositivitypreservinghighorderdiscontinuousGalerkinschemesforcompressibleEulerequationsonrectangularmeshes[J].J.Comput.Phys.2010(229)：8918-8934.

[2]REEDWH,HILLTR.Triangularmeshmethodsforneutrontransportequation[R].Tech.ReportLA-UR,1973:473-479.

[3]COCKBURNB,SHUCW.TVBRunge-KuttalocalprojectiondiscontinuousGalerkinfiniteelementmethodforconservationlawsII:generalframework[J].MathematicsofComputation,1989,52(186):411- 435.

[4]COCKBURNB,LINSY,SHUCW.TVBRunge-KuttalocalprojectiondiscontinuousGalerkinfiniteelementmethodforconservationlawsIII:onedimensionalsystems[J].JournalofComputationPhysics,1989,84(1):90-113.

[5]COCKBURNB,HOUS,SHUCW.TheRunge-KuttalocalprojectiondiscontinuousGalerkinfiniteelementmethodforconservationlawsIV:themultidimensionalcase[J].MathematicsofComputation,1990,54(190):545-581.

[6]GOTTLIEBS,SHUCW,TADMORE.StrongStability-Preservinghighordertimediscretizationmethods[J].SocietyforIndustrialandAppliedMathematics,1990，43(1)：89-112.

[7]ZHANGX,SHUCW.Positivity-preservinghighorderfinitedifferenceWENOschemesforcompressibleEulerequations[J].J.Comput.Phys.2012,231:2245-2258.

[8]ZHANGX,SHUCW.Positivity-preservinghighorderdiscontinuousGalerkinschemesforcompressibleEulerequationswithsourceterms[J].J.Comput.Phys.2011(230)：1238-1248.

[9]FEDKIWR,ASLAMT,MERRIMANB.Anon-oscillatoryEulerianapproachtointerfacesinmultimaterialflows(theGhostFluidMethod)[J].JournalofComputationalPhysics,1999,152:457-492.

[10]WANGCW,LIUTG,KHOOBC.Areal-ghostfluidmethodforthesimulationofmulti-mediumcompressibleflow[J].SIAMJ.Sci.Comput，2006,28:278-302.

[責(zé)任編輯王新奇]

Vol.18No.3Jul.2015

劉慧敏(1981—)，女，河南鄭州人，鄭州成功財(cái)經(jīng)學(xué)院共同學(xué)科部講師，碩士，主要從事微分方程研究.