1. 引言
1959年,Godunov首次提出以Riemann问题解为基础来构造网格物理量均值的Godunov格式,它有效地解决了所有内部状态变量的局部黎曼问题 [1] [2]。随后,Van Leer首次讨论了二阶Godunov格式的构造得到了MUSCL (Monotone Upwind Scheme for Conservation Laws)格式,这个格式包括了保单调性的插值和利用Riemann解得到的迎风通量,提高了数值精度和分辨率 [2]。MUSCL格式已经被许多人研究 [3] [4] [5],并被广泛用于科学和工程问题的模拟 [5]。
1959年,Wiin-Nielsen提出了一种修正轨迹方法 [6],随后Sawyer将其发展为Semi-Lagrangian方法 [7]。该方法具有固定的(欧拉)数值网格,根据Riemann不变量沿特征线不变的性质,从而推进每个时间步,使得时间步长不受CFL条件的限制,减少计算工作量和数值扩散 [8]。Semi-Lagrangian方法成功地推广到了有限差分,有限元和谱方法,并广泛应用于天气预报,海洋学模型的数值模拟和线性对流的建模 [9] [10]。1982年,Robert将半隐式(SI)和半拉格朗日(SL)这两种方法相结合,在NWP模型中实现了显著的效率增益,为许多环境应用产生了许多有效的SISL模型 [11]。1999年,Strain将水平集方法和半拉格朗日方法相结合,能有效、自适应和模块化求解几何运动界面问题 [10] [12]。2006年,Bonaventura和Restelli等人将DG方法与半拉格朗日方法相结合为SLDG,能有效地融合了两种方法优点,又避免了缺点 [13],对于shallow water方程模拟,SLDG具有较高的精度和有效性,并可自动选择局部逼近程度 [11]。2010年,Qiu等人将Strang时间分裂与空间高阶WENO重构相结合,应用于Semi-Lagrangian方法来求解Vlasov方程,比使用显式时间步长的Vlasov-Poisson系统的欧拉公式计算更高效,更灵活 [14]。
然而,对于一维或高维Euler方程组,求解特征线方程仍然是Semi-Lagrangian方法的难点。2011年,Lentine等人基于文 [15] 提出的分裂方法,用保守Semi-Lagrangian对流代替Meno对流格式求解方程组的对流项,但局部产生虚假振荡 [16]。目前尚未在现有文献中发现更好的处理方案。2015年,吴浪根据文 [17] 提出的方法,针对守恒律方程给出了几种高阶Semi-Lagrangian格式,通过数值模拟验证了该方法的高阶精度和无振荡性质 [18]。处理过程用到了一维Euler方程组的Riemann不变量,并反解出原物理量。这个过程增加了计算量,针对不同维数Euler方程组需要给出不同的Riemann不变量,不方便推广到多维问题,即便是通过交替方向法也不行。
本文给出了一种高阶数值通量的构造方法,将其应用于有限体积(差分) ENO,WENO和DG,即相应的Semi-Lagrangian方法。针对于一维Euler方程组,在特征空间上给出了一种新的处理特征线的方法。数值实验表明:高阶数值通量能应用于有限体积(差分) ENO,WENO和DG等格式,针对Euler方程组的特征线处理方法能保持空间离散格式本身的本质无震荡性,能获得更高的分辨率且具备更高的计算效率。
2. 二阶数值通量的推导
为了推导双曲型方程数值计算格式中的二阶数值通量,我们考察如下方程
(1)
(2)
已知
在第j个单元(
)端点的单侧极限分别为
和
,容易给出它在
内的插值多项式
(3)
其中
。
过点
的特征线与直线
相交于一点X
(4)
其中
。利用一阶Runge-Kutta方法,可得
(5)
将X代入
,得
(6)
其中
是近似代替
。
最后,利用梯形公式可计算出
内
处的二阶数值通量
(7)
3. 二阶数值通量在DG格式的应用
为了将二阶数值通量应用于DG格式,我们将方程(1)两边同时乘
,并在时空单元
上求积分,
(8)
其中
为区间
上的光滑函数。由格林公式,可得
(9)
我们再将方程(2)两边同乘
并在单元
上积分得
(10)
这样就得到了方程(1-2)的弱形式(9-10)。
参照经典二阶DG格式,取Legendre多项式
为局部基函数,并设
,其中
。将方程(9)中的
替换成
,则有
当l = 0时,(9)式简写为
(11)
当l = 1时,(9)式简写为
(12)
将(12)式中的二重积分项简记为K,利用如下梯形公式计算
(13)
其中
继续利用梯形公式(或其他积分公式),
由方程(7,11)得到的单元均值
来近似代替,得
(14)
再利用上一节所介绍的二阶数值通量的推导方法,可求解(12)式中的
这样,我们得到相应的单步二阶DG格式(10)~(12)。
4. 一维Euler方程组上的二阶数值通量的推导
一维Euler方程组可以表示为
(15)
其中压力
,比热比
。Jacobian矩阵
,A的特征值为
方程组(15)化成特征空间形式
其中R为右特征矩阵,
为左特征矩阵,
。
已知
在第j个单元(
)端点的单侧极限分别为
和
可以给出它在
内的插值多项式
特征线方程
与直线
相交,有
(16)
利用二阶Runge-Kutta法求解常微分方程(16),再将
代入
得到
。令
,则
。
最后,利用梯形公式可计算出
内
处的二阶数值通量
(17)
5. 数值算例
求解方程(1)的守恒型有限体积(差分)格式通常可以写成
(18)
我们将由公式(7)或(17))计算所得的通量替代(18)式中的
,这样就得到单步二阶格式.
算例1:一维Burger方程
考虑方程(1-2),其中
,采用周期边界条件。为了考
Table 1. Numerical error table for single-step second-order ENO, DG format and original format at t = 0.5
表1. 单步二阶ENO,DG格式与原格式在t = 0.5时刻的数值误差表
Table 2. Numerical error meter of third order WENO scheme in = 0.5 FVWENO3-HNF
表2. t = 0.5时三阶WENO格式的误差表
察格式的截断误差阶,我们给出
(尚未产生间断)的误差表,见表1~表2。其中后缀“HNF”表示使用公式(7)的单步二阶格式,另一个表示使用相应的二阶SSP Runge-Kutta时间离散格式。
由表1~表2可以看出:1) 本文给出的二阶数值通量能应用于ENO,WENO和DG等格式,得到相应的单步二阶格式。2) 单步二阶格式的计算误差比相应原格式更小(三阶WENO格式当网格充分细密时出现超收敛,误差更小)。3) 单步二阶格式所耗CPU时间更少,当N = 6400时,FV-ENO2-HNF格式比原格式节约了48%,FD-ENO2-HNF格式比原格式节约了56%,DG2-HNF格式比原格式节省了88%。
Figure 1. The local comparison of real and numerical solutions in t = 1.5
图1. t = 1.5真解与数值解的局部对比
当t = 1.5时,该问题已产生激波,见图1。由图可知:单步二阶格式在处理激波间断时具有更高的分辨率,模拟效果更好。
算例2:Sod激波管问题
考虑方程(15),其中
。采用紧支边界,初始条件:t=0时刻,
真解N=1600,数值解N=100时,使用二阶数值通量单步ENO格式与标准二阶FD-ENO2格式的对比(见图2~图3)。
Figure 2. The overall comparison between the true solution and the numerical solution in t = 0.15
图2. t = 0.15时真解与数值解的整体对比
由图2~图3可以看出:使用二阶数值通量的单步格式比原格式具有更高的分辨率。当N = 1000时,
(a) 密度 (b) 速度 (c) 压力
Figure 3. The local comparison between the true solution and the numerical solution in t = 0.15
图3. t = 0.15时真解与数值解的局部对比
FD-ENO2-HNF格式所消耗的CPU时间为282毫秒(ms),标准FD-ENO2格式所消耗的CPU时间为400毫秒(ms),FD-ENO-HNF格式比标准FD-ENO2格式节约了29.5%左右。
6. 总结
本文以二阶数值通量为例,给出了一种高阶数值通量的构造方法,能应用于有限体积(差分) ENO,WENO和DG格式得到单步高阶Semi-Lagrangian格式。格式天然具备守恒性、紧致性、并行性,且并行效率比时间离散应用SSP Runge-Kutta方法的相应格式更高。本文还给出了一种在特征空间处理特征线的方法,能应用于Euler方程组,比原格式具备更好的模拟效果和更高的效率。同时,这种方法能利用交替方向法推广到多维问题,并无需推导其Riemann不变量,也为Semi-Lagrangian方法在多维Euler方程组中的研究与应用提供了一种有效的途径。
参考文献
NOTES
*通讯作者