一维低密度低压力多介质流动问题的高精度数值方法

一维低密度低压力多介质流动问题的高精度数值方法

吴招生

(南京航空航天大学 理学院数学系，南京 210016)

摘要：在爆炸波等多介质流体力学问题的数值模拟中，会出现低密度低压力区域，由于数值误差的影响，导致在计算过程中该区域的密度和压力为负值，使得计算难以继续.文献[1]给出了求解理想气体欧拉方程的正保护RKDG(Runge-Kutta Discontinuous Galerkin)方法.在此基础上，给出了多介质流体力学问题数值模拟的正保护RKDG方法，结合RGFM(Real Ghost Fluid Method)界面处理方法，对一维多介质可压缩流体进行了数值试验，结果表明该方法能有效地避免压力为负的情况，并能准确地捕捉各类间断.

关键词：RKDG方法；低密度；低压力；正保护限制器

中图分类号：O241.82文献标志码：A

文章编号：1008-5564(2015)03-0026-05

收稿日期：2015-04-20

作者简介：刘小刚(1983—)，男，陕西延安人，西北大学现代学院助教，硕士，主要从事微分方程与动力系统研究;

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提出的，用来求解中子运输问题[2].20世纪80年代后期，Cockburn和C.W.Shu等构造了求解一维标量守恒律方程[3]的高阶RKDG有限元方法，并且将该方法成功地推广到多维标量守恒律[4]和多维守恒律方程组[5].RKDG方法的优点是具有高精度且易于处理复杂边界.此后，RKDG方法开始广泛应用于流体力学问题的数值模拟.

高精度守恒型差分格式[6]由于其能够准确捕捉激波及其它物理现象，在流体力学数值计算中得到了广泛的应用.物理上密度和压力应该为正，但在爆炸波和高速射流等问题中通常会出现低密度和低压力区域，在用高精度格式求解时，由于数值误差的影响，导致密度或压力出现负值，使得计算难以继续进行.文献[1]和[7]分别针对高精度RKDG格式和WENO格式构造了正保护格式，并推广到含有源项的欧拉方程组[8].本文在此工作的基础上，开展多介质流动问题高精度数值计算格式正保护性质的研究.

对于多介质流动问题，界面两边是不同性质的流体，它们的状态方程可能会有很大差异，流体的密度会出现大梯度的变化，如果还采用单介质流动问题的数值模拟方法，则会出现数值不稳定的情况，在界面处会产生非物理震荡，甚至使得计算难以进行.因此，在多介质流动问题的数值计算中，界面处理方法的研究是关键.Fedkiw等[9]提出的Ghost方法，把多介质流体问题转化为多个单介质问题求解.在此基础上提出的基于求解Riemann问题的界面处理方法(简称RGFM方法[10])，给出了更加精确的界面边界条件，能准确地捕捉界面位置.本文利用RGFM方法定义界面边界条件，将多介质流动问题转化为单介质流动问题，利用所给出的具有正保护性质的高精度RKDG格式进行计算.

本文在不改变文献[1]中正保护限制器性质的情况下，利用压力函数是凹函数的性质，进行线性插值求解，避免了求解一元二次方程，节省了计算量，并将此正保护限制器推广应用于水的状态方程(Tait方程).本文采用三阶具有正保护性质的RKDG格式结合RGFM界面处理方法数值模拟一维多介质可压缩Euler方程组，对具有大密度比的一维气-气，气-液两种介质可压缩流体问题进行了数值实验.

1方程

1.1控制方程

考虑一维无粘可压缩流体力学方程组：

(1)

1.2状态方程

理想气体的状态方程：

p=(γ-1)ρe

(2)

水的状态方程：

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

(3)

其中γ和PC均为正常数，将在具体算例中给出.

2数值方法

2.1RKDG方法

考虑Euler方程组(1)的初值问题：

(4)

其中[a,b]为计算区域.

(5)

(6)

(7)

(8)

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

(9)

(9)式可以写成：

(10)

方程组(10)的时间离散采用三阶TVB Runge-Kutta方法[6].

2.2RGFM界面处理方法

假设流体Ⅰ与流体Ⅱ的界面位于网格格点i和i+1之间，对流体Ⅰ而言，DG格式边界计算时需要网格单元i+1的状态值，这个单元的状态值由在界面处构造的Riemann问题的解重新赋值.Riemann问题的初始状态取值如下：

(11)

图1　流体I界面边界条件的定义

分别定义了流体Ⅰ与流体Ⅱ的界面边界条件后，再利用所得到的高精度RKDG格式分别求解流体Ⅰ与流体Ⅱ，得到下一个时间步的数值解.

3正保护限制器

正保护限制器的具体实施分两大步骤：第一步，对密度修正，保证每个网格单元上的3个Gauss-Lobatto点的密度近似函数的值大于零；第二步，当流体是理想气体时，对压力修正.当流体是液体时，对内能修正.保证每个网格单元上的3个Gauss-Lobatto点的压力近似函数或者内能的近似函数的值大于零.限制器放在时间t推进到t+Δt的每个RK子步后进行限制，不影响原计算格式的精度和守恒性.

下面给出不同流体情形下的正保护限制器的具体实施和算法.

3.1流体为理想气体

对于理想气体而言，密度和压力在物理意义下始终大于零.定义集合G：

(12)

可以证明集合G是凸集合，压力函数p是凹函数.

3.2流体为液体

对于液体而言，密度和内能在物理意义下始终为正.定义集合G：

(13)

(1)密度修正：

(2)内能修正：

4数值实验

4.1气-氦激波管问题

这是一个一维气-气多介质问题，初始条件如下：

由于数值误差的影响，不使用正保护限制器，在低压区压力出现负值，程序崩溃.取网格数800，TVB参数取为(M1,M2,M3)=(1010,10,1010).随着时间推进，初始间断分解为向右的强激波和向左的稀疏波，图2～图4分别给出了t=6时刻的密度、压力和速度分布图.可见，本文所给出的算法能够保证低压区的数值计算，保证压力为正.

图2　气-氦激波管问题 密度

图3　气-氦激波管问题 压力

图4　气-氦激波管问题 速度

4.2气-液激波管问题

这是一个一维气-液多介质问题，初始条件如下：

由于数值误差的影响，不使用正保护限制器，水中内能出现负值，程序崩溃.取网格数200，TVB参数取为(M1,M2,M3)=(50,50,50).随着时间推进，初始间断分解为向左的强激波和向右的稀疏波，图5～图7分别给出了t=0.000 24时刻的密度、速度和内能分布图.可见，本文所给出的算法能够保证内能为正.

图5　气-液激波管问题 密度

图6　气-液激波管问题 速度

图7　气-液激波管问题 内能

5总结与展望

本文对原有正保护限制器进行修正，避免了求解一元二次方程，节省了计算量，并将此正保护限制器推广到水的状态方程(Tait方程)，结合RGFM界面处理方法，应用于一维多介质流体力学问题的数值模拟.数值试验表明该方法能准确捕捉界面和其它间断，体现了算法的保正性.

[参考文献]

[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.

[责任编辑王新奇]

Vol.18No.3Jul.2015

刘慧敏(1981—)，女，河南郑州人，郑州成功财经学院共同学科部讲师，硕士，主要从事微分方程研究.