一维双曲型守恒律方程的Lax-Wendroff型中心间断伽辽金方法
Lax-Wendroff Type Central Discontinuous Galerkin Method for One-Dimensional Hyperbolic Conservation Law Equations
摘要: 物理系统中波动、传播等现象通常用双曲型守恒律方程的数学模型来描述，特别是在流体力学领域尤为重要。针对此类方程，我们考虑了Lax-Wendroff型中心间断伽辽金方法。该方法首先采用Lax-Wendroff型时间离散方法，也就是通过泰勒级数展开处理时间导数，然后在空间上运用中心间断伽辽金方法，从而避免了传统的多步时间积分方法。最后我们对多个双曲型守恒律方程开展数值实验，验证所提出方法在计算效率和精度上的有效性。
Abstract: In physical systems, phenomena like wave fluctuation and propagation are often described using hyperbolic conservation law equations, which play a crucial role in fluid mechanics. To solve these equations, we employ the Lax-Wendroff central discontinuous Galerkin method. This approach begins with the Lax-Wendroff time discretization, where time derivatives are managed through a Taylor series expansion. It then incorporates the central discontinuous Galerkin method for spatial discretization and effectively eliminates the need for traditional multi-step time integration schemes. Finally, numerical experiments on various hyperbolic conservation law equations are constructed to validate the effectiveness of our method in terms of both computational efficiency and accuracy.
文章引用：罗艺. 一维双曲型守恒律方程的Lax-Wendroff型中心间断伽辽金方法[J]. 应用数学进展, 2025, 14(2): 376-387. https://doi.org/10.12677/aam.2025.142078

1. 引言

双曲型守恒律方程是一类描述物理系统中某些守恒量随时间和空间演化的双曲型方程[1]，广泛应用于流体力学、气体动力学和电磁学等领域。由于这些方程的解经常表现出激波等间断性，且数值模拟可能受到震荡或数值扩散的影响，目前，开发有效求解双曲型守恒律方程的数值方法在计算数学领域的研究中依然面临着很大的挑战。

传统的数值方法包括有限差分法、有限体积法和有限元方法等等，这些方法经过多年的发展和改进，已在许多实际问题中得到广泛应用。在这些方法中，间断伽辽金方法最早由Reed和Hill于1973年引入，用于中子输运问题[2]，并在1988年得到了进一步发展[3]，成为求解双曲型守恒律方程的重要工具。间断伽辽金方法因其高阶精度和处理复杂几何的灵活性而受到广泛关注，它通过在每个单元内采用分段多项式函数，能够实现高阶空间精度，并自然处理解中的间断性[4]。然而，在处理强间断性时，间断伽辽金方法仍可能出现数值震荡，且计算成本较高，尤其是在多维问题中，通常需要较小的时间步长以确保稳定性。

为了解决这些问题，近年来开发了几种改进的间断伽辽金方法。其中，中心间断伽辽金方法作为间断伽辽金方法的一种重要变型[5] [6]，通过在交错网格上处理离散项，避免了依赖于单元边界上的精确解或数值黎曼解。它通过不在单元边界处应用数值通量计算两组数值解，从而降低了计算复杂度。目前，中心间断伽辽金方法已成功应用于多个研究领域，如理想磁流体动力学方程[7] [8]、Hamilton-Jacobi方程[9]、Green-Naghdi方程[10]、Camassa-Holm方程[11]、广义Korteweg-de Vries方程[12]和波动方程[13]等等。

在时间离散化方面，三阶全变分递减(TVD) Runge-Kutta方法是一种常用的离散化方法[4] [14]，它通过三步迭代更新提供高阶精度。然而，TVD型Runge-Kutta方法的精度通常仅限于四阶[15]，即便在放宽TVD要求的情况下，高阶Runge-Kutta方法仍需要比精度阶数多更多的步数，这增加了计算复杂性，使得其在高阶时间离散化中效率较低。此外，对于具有非线性或复杂间断结构的问题，TVD型Runge-Kutta方法可能会产生高频震荡。

为克服这些局限性，我们考虑了Lax-Wendroff型时间离散方法[16] [17]。Lax-Wendroff型时间离散方法也被称为泰勒型或Cauchy-Kowalewski型，指的是在时间导数项中应用泰勒展开，或者在偏微分方程中的类似Cauchy-Kowalewski过程[1]。该方法基于经典的Lax-Wendroff方法，通过运用控制方程及其微分形式，将所有时间导数转换为空间导数，从而实现高阶时间精度并保持显式方案。与Runge-Kutta方法[4] [14]相比，Lax-Wendroff方法降低了计算复杂度，并且在波传播问题中表现出色，特别是在捕捉激波和间断方面。因此，将Lax-Wendroff型时间离散方法引入中心间断伽辽金框架，可以将高阶时间精度与增强的间断性处理能力相结合，从而得到更为稳健的数值方法。该研究旨在开发一种Lax-Wendroff型中心间断伽辽金方法，为求解双曲型守恒律方程提供一种既能解决光滑解又能处理间断解的方法，同时保持高精度和计算效率。

基于上述背景，本文探讨了Lax-Wendroff型中心间断伽辽金方法的数学公式，并展示了一系列数值实验的结果。实验结果表明，Lax-Wendroff型中心间断伽辽金方法不仅在高阶精度下保持了优异的性能，还展示了处理复杂间断解的卓越能力。本文的组织结构如下：第2部分详细描述了Lax-Wendroff型中心间断伽辽金方法在一维标量守恒律方程和Euler方程中的构造与实现；第3部分提供了多个数值实例用以验证方法的有效性和准确性；最后，在第4部分进行了总结和展望。

2. 数值方法

2.1. 一维守恒律方程的Lax-Wendroff型中心间断伽辽金方法

本节主要基于Lax-Wendroff型中心间断伽辽金方法求解一维守恒律方程(1)。首先，一维标量守恒律方程的数学表达式为：

{ u t + f ( u ) x = 0 ,   ( x , t ) ∈ Ω × ( 0 , T ] u ( x , 0 ) = u 0 ( x ) . (1)

其中 u ( x , t ) 表示关于位置x和时间t的物理量， f ( u ) 是关于u的流量函数。针对方程(1)，Lax-Wendroff型时间离散化如下，设 Δ t 为时间步长，其中 Δ t = t n + 1 − t n ，如公式(2)所示，我们对 t n + 1 时刻上的解在时间 t n 上进行泰勒展开以实现三阶时间精度，那么我们得到：

u ( x , t n + Δ t ) = u ( x , t n ) + Δ t   u t + Δ t 2 2 u t t + Δ t 3 6 u t t t + O ( Δ t 3 ) . (2)

根据一维守恒律方程(1)，时间导数项 u t 、 u t t 和 u t t t 可以转化为空间导数项，其表达形式如下：

u t = − f ( u ) x = − f ′ ( u ) u x , u t t = ( f ′ ( u ) f ( u ) x ) x , u t t t = − ( f ″ ( u ) ( f ( u ) x ) 2 + f ′ ( u ) ( f ′ ( u ) f ( u ) x ) x ) x ,

然后，我们结合上述公式，公式(2)可以被改写为：

u ( x , t n + Δ t ) = u ( x , t n ) + Δ t F ( u ) x , (3)

其中 F ( u ) = − f ( u ) + Δ t 2 f ′ ( u ) f ( u ) x − Δ t 2 6 ( f ″ ( u ) ( f ( u ) x ) 2 + f ′ ( u ) ( f ′ ( u ) f ( u ) x ) x ) .

经过Lax-Wendroff型时间离散化后，我们在空间上采用中心间断伽辽金方法。设 { x j } j 为计算域 Ω = [ a , b ] 的一个划分，记单元为 I j = [ x j − 1 / 2 , x j + 1 / 2 ] 和 I j + 1 / 2 = [ x j , x j + 1 ] ，这里 x j + 1 / 2 = 1 2 ( x j + x j + 1 ) 。此外，我们定义两个离散函数空间，分别与重叠网格 { I j } j 和 { I j + 1 / 2 } j 相关，

V h C = V h C , k = { v : v | I j ∈ P k ( I j ) , ∀ j } , V h D = V h D , k = { v : v | I j + 1 / 2 ∈ P k ( I j + 1 / 2 ) , ∀ j } ,

其中I表示在单元I上次数不超过k的多项式空间，v在单元内是连续函数。

我们发现，公式(3)中的数值通量 F ( u ) 存在高阶导数，对此采用了中心局部间断伽辽金(CLDG)方法，为了重构数值通量 F ( u ) 的高阶导数，我们引入了两个辅助变量：

p = f ( u ) x ,   q = ( f ′ ( u ) p ) x ,

并且将公式(3)改写为：

u ( x , t n + Δ t ) = u ( x , t n ) + Δ t   F ( u , p , q ) x = u ( x , t n ) + Δ t   [ − f ( u ) + Δ t 2 f ′ ( u ) p − Δ t 2 6 ( f ″ ( u ) p 2 + f ′ ( u ) q ) ] x . (4)

基于上述公式，接下来我们将阐述CLDG方法的两组形式。

CLDG方法的第一组形式：

为了简便，只给出计算 u C ( x , t n + Δ t ) 、 p C 和 q C 的计算公式，关于 u D ( x , t n + Δ t ) 、 p D 和 q D 的计算公式类似可得。在第一组网格中，我们使用空间 V h C ，其半离散形式如下：寻找 u C , p C , q C ∈ V h C , k ，使得对于任意 v C , v ˜ C , v ^ C ∈ V h C , k ，有：

∫ I j u C ( x , t n + Δ t ) v C   d x = ∫ I j u C ( x , t n ) v C   d x − Δ t ∫ I j F ( u D , p D , q D ) ∂ x v C   d x   + Δ t [ [ F ( u D , p D , q D ) ] j + 1 2 v C ( x j + 1 2 − ) − [ F ( u D , p D , q D ) ] j − 1 2 v C ( x j − 1 2 + ) ] , (5)

∫ I j p C v ˜ C   d x = − ∫ I j f ( u D ) ∂ x v ˜ C   d x + [ f ( u D ) ] j + 1 2 v ˜ C ( x j + 1 2 − ) − [ f ( u D ) ] j − 1 2 v ˜ C ( x j − 1 2 + ) , (6)

∫ I j q C v ^ C   d x = − ∫ I j f ′ ( u ) p D ∂ x v ^ C   d x + f ′ ( u ) ( p D ) j + 1 2 v ^ C ( x j + 1 2 − ) − f ′ ( u ) ( p D ) j − 1 2 v ^ C ( x j − 1 2 + ) . (7)

CLDG方法的第二组形式：

第二组形式相比第一组形式增加了一项数值耗散项，用于处理不同单元上两个解之间的差异，并且该项的引入在分片线性情况下恢复了最优收敛率，且在非光滑初始条件下提供了更小的误差。其半离散形式如下：寻找 u C , p C , q C ∈ V h C , k ，使得对于任意 v C , v ˜ C , v ^ C ∈ V h C , k ，有：

∫ I j u C ( x , t n + Δ t ) v C   d x = 1 τ m a x ∫ I j ( u C − u D )   v C   d x + ∫ I j u C ( x , t n ) v C   d x − Δ t ∫ I j F ( u D , p D , q D ) ∂ x v C   d x   + Δ t [ F ( u D , p D , q D ) ] j + 1 2 v C ( x j + 1 2 − ) − Δ t [ F ( u D , p D , q D ) ] j − 1 2 v C φ ( x j − 1 2 + ) , (8)

∫ I j p C v ˜ C   d x = − ∫ I j f ( u D ) ∂ x v ˜ C   d x + [ f ( u D ) ] j + 1 2 v ˜ C ( x j + 1 2 − ) − [ f ( u D ) ] j − 1 2 v ˜ C ( x j − 1 2 + ) , (9)

∫ I j q C v ^ C   d x = − ∫ I j f ′ ( u ) p D ∂ x v ^ C   d x + f ′ ( u ) ( p D ) j + 1 2 v ^ C ( x j + 1 2 − ) − f ′ ( u ) ( p D ) j − 1 2 v ^ C ( x j − 1 2 + ) ， (10)

其中， τ max 是由Courant-Friedrichs-Lewy (CFL)条件所决定的时间步长的上界，即 τ max = c Δ x ，其中c是由稳定性决定的CFL参数。

2.2. 一维Euler方程的Lax-Wendroff型中心间断伽辽金方法

一维Euler方程的Lax-Wendroff型中心间断伽辽金方法与一维守恒律方程的类似，仅在自变量的数量上有所不同。我们注意到，Lax-Wendroff型中心间断伽辽金方法可以扩展用于求解双曲守恒律系统，为了说明这一方法，我们考虑一维Euler方程：

u t + f ( u ) x = 0 , (11)

u = ( ρ m E ) ,   f ( u ) = ( ρ u ρ u 2 + p ( E + p ) u ) ,

其中 m = ρ u ， E = 1 2 ρ u 2 + ρ e ， p = γ ^ ρ e ， ρ 表示密度，m表示动量，u表示速度，E表示总能量，e表示

内能，p表示压力， γ ^ = γ − 1 且 γ > 1 是比热容(空气中 γ = 1.4 )。接下来，我们通过运用方程中每个变量的微分形式，展示Euler方程的三阶Lax-Wendroff型时间离散形式，得到以下结果：

ρ ( t n + Δ t ) ≈ ρ ( t n ) + Δ t ρ t + Δ t 2 2 ρ t t + Δ t 3 6 ρ t t t = ρ ( t n ) + Δ t [ ρ t + Δ t 2 ρ t t + Δ t 2 6 ρ t t t ] = ρ ( t n ) + Δ t [ − m − Δ t 2 m t − Δ t 2 6 m t t ] x , (12)

m ( t n + Δ t ) ≈ m ( t n ) + Δ t m t + Δ t 2 2 m t t + Δ t 3 6 m t t t = m ( t n ) + Δ t [ m t + Δ t 2 m t t + Δ t 2 6 m t t t ] = m ( t n ) + Δ t [ − ( γ ^ E + 3 − γ 2 m u ) − Δ t 2 ( γ ^ E t + 3 − γ 2 m t u + 3 − γ 2 m u t )   − Δ t 2 6 ( γ ^ E t t + 3 − γ 2 m t t u + 3 − γ 2 m u t t + ( 3 − γ ) m t u t ) ] x , (13)

E ( t n + Δ t ) ≈ E ( t n ) + Δ t E t + Δ t 2 2 E t t + Δ t 3 6 E t t t = E ( t n ) + Δ t [ E t + Δ t 2 E t t + Δ t 2 6 E t t t ] = E ( t n ) + Δ t [ − ( γ E u + γ ^ 2 m u 2 ) − Δ t 2 ( γ E t u + γ E u t + γ ^ 2 m t u 2 + γ ^ m u u t )   − Δ t 2 6 ( γ ^ E t t u + γ E u t t + 2 γ E t u t − γ ^ 2 m t t u 2 + γ ^ m ( u t 2 + u u t t ) + 2 γ ^ m t u u t ) ] x . (14)

上述方程中的变量 u x , u x x , u t , u t x , u t t , m t , m t x , m t t , E t , E t x , E t t 已在[18]中有表述。基于[18]中变量的微分形式，我们引入辅助变量：

p 1 = ρ x ,   p 2 = m x ,   p 3 = E x ,   p 4 = ρ x x ,   p 5 = m x x ,   p 6 = E x x ,

从而将其重写为：

u x = m x ρ − u ρ x ρ = p 2 ρ − u p 1 ρ ,

m t = − ( γ ^ E x + 3 − γ 2 m x u + 3 − γ 2 m u x ) = − ( γ ^ p 3 + 3 − γ 2 p 2 u + 3 − γ 2 m u x ) ,

E t = − ( γ E x u + γ E u x − γ ^ 2 m x u 2 − γ ^ m u u x ) = − ( γ p 3 u + γ E u x − γ ^ 2 p 2 u 2 − γ ^ m u u x ) ,

u x x = − 2 m x ρ x ρ 2 + m x x ρ + m ( 2 ρ x 2 ρ 3 − ρ x x ρ 2 ) = − 2 p 2 p 1 ρ 2 + p 5 ρ + m ( 2 p 1 2 ρ 3 − p 4 ρ 2 ) ,

u t = m t ρ + u m x ρ = m t ρ + u p 2 ρ ,

m t x = − ( γ ^ E x x + 3 − γ 2 ( m x x u + 2 m x u x + m u x x ) ) = − ( γ ^ p 6 + 3 − γ 2 ( p 5 u + 2 p 2 u x + m u x x ) ) ,

u t x = m t x ρ + m x 2 + m m x x − m t ρ x − 2 u ρ x m x ρ 2 = m t x ρ + p 2 2 + m p 5 − m t p 1 − 2 u p 1 p 2 ρ 2 ,

u t t = 2 m t m x ρ 2 + m t t ρ + u ( 2 m x 2 ρ 2 + m t x ρ ) = 2 m t p 2 ρ 2 + m t t ρ + u ( 2 p 2 2 ρ 2 + m t x ρ ) ,

E t x = − ( γ E x x u + 2 γ E x u x + γ E u x x − γ ^ 2 m x x u 2 − γ ^ ( 2 m x u u x + m u x 2 + m u u x x ) ) = − ( γ p 6 u + 2 γ p 3 u x + γ E u x x − γ ^ 2 p 5 u 2 − γ ^ ( 2 p 2 u u x + m u x 2 + m u u x x ) ) ,

m t t = − ( γ ^ E t x + 3 − γ 2 ( m t x u + m x u t + m t u x + m u t x ) ) = − ( γ ^ E t x + 3 − γ 2 ( m t x u + p 2 u t + m t u x + m u t x ) ) ,

E t t = − ( γ ( E t x u + E x u t + E t u x + E u t x ) − γ ^ 2 m t x u 2 − γ ^ ( m x u u t + m t u u x + m u t u x + m u u t x ) ) = − ( γ ( E t x u + p 3 u t + E t u x + E u t x ) − γ ^ 2 m t x u 2 − γ ^ ( p 2 u u t + m t u u x + m u t u x + m u u t x ) ) .

接下来，我们将引入的辅助变量和上述变量的微分形式代入式(12)~(14)中，改写成与公式(4)类似的守恒律形式，并构造成与一维守恒律方程类似的三阶CLDG方法。同样，通过在截断的泰勒展开式中加入额外的导数项，可以实现更高阶的精度。

3. 数值算例

这一部分，我们展示了一些数值算例用以验证所提出的数值方法的有效性以及准确性。实验涵盖了方程连续解的精度测试和间断解的数值计算，其结果采用图表的形式进行展示分析。值得注意的是，在时间离散过程中，依据一维情况下的时间步长公式 Δ t = c   Δ x ，为避免较大CFL数引发的时间精度退化问题，特别是在处理非线性方程时，我们选择了较小的CFL数，以确保数值解的稳定性和精度。通过合理的CFL数调控，使得实验能成功验证该方法在不同问题下的有效性与准确性。此外，对于具有间断解的Burgers方程和Euler方程，我们应用了总变差有界(TVB)限制器，并给出了相应的TVB常数。

3.1. 线性对流方程的精度测试

我们对一维线性对流方程进行精度测试，该方程由以下表达式给出：

u t + u x = 0 ,   x ∈ [ 0 , 2 π ] ,

其初始条件为 u ( x , 0 ) = sin ( x ) ，精确解为 u ( x , t ) = sin ( x − t ) ，并在实验中采用了周期边界条件。接着，我们计算出 t = 1.0 时刻该方程的数值解在 L 1 、 L 2 和 L ∞ 范数意义下的误差及收敛阶，结果如表1所示。根据表格中的数值结果，我们发现在 P k 情况下的收敛阶达到了 k + 1 ，这与误差估计理论一致，此结果进一步验证了Lax-Wendroff型中心间断伽辽金方法在处理线性对流方程中的有效性。

Table 1. The L 1 , L 2 and L ∞ errors and orders of accuracy of the numerical solution for the 1D linear advection equation

表1. 一维线性对流方程数值解的 L 1 、 L 2 和 L ∞ 误差及收敛阶

N x

L 1 误差

收敛阶

L 2 误差

收敛阶

L ∞ 误差

收敛阶

k = 1

10

6.02E−02

2.93E−02

2.29E−02

20

1.33E−02

2.18

6.90E−03

2.09

6.60E−03

1.79

40

3.09E−03

2.10

1.69E−03

2.03

1.77E−03

1.90

80

7.58E−04

2.03

4.22E−04

2.01

4.55E−04

1.96

160

1.89E−04

2.00

1.05E−04

2.00

1.16E−04

1.98

k = 2

10

2.69E−03

1.49E−03

1.48E−03

20

3.23E−04

3.06

1.85E−04

3.01

1.89E−04

2.98

40

4.00E−05

3.01

2.32E−05

3.00

2.38E−05

2.99

80

5.03E−06

2.99

2.90E−06

3.00

3.12E−06

2.93

160

6.77E−07

2.89

3.90E−07

2.89

5.57E−07

2.49

3.2. 无粘性Burgers方程的精度测试

我们对一维无粘性Burgers方程进行精度测试，该方程如下：

u t + ( u 2 2 ) x = 0 ,   x ∈ [ 0 , 2 π ] ,

方程的初始条件是 u ( x , 0 ) = sin ( x ) + 2.0 ，并在实验中采用了周期边界条件。值得注意的是，对于无粘性Burgers方程，我们基于初始条件运用牛顿迭代法获得了参考解。表2中展示了方程数值解在 L 1 、 L 2 和 L ∞ 范数意义下的误差及收敛阶，其中计算时间为 t = 0.3 。结果表明，尽管求解非线性方程的复杂性较高，可能会出现震荡和数值不稳定性，但Lax-Wendroff型中心间断伽辽金方法仍然在 P k 情况下的收敛阶数达到了 k + 1 ，这表明数值结果与理论预期一致。

Table 2. The L 1 , L 2 and L ∞ errors and orders of accuracy of the numerical solution for the 1D inviscid Burgers equation

表2. 一维无粘性Burgers方程数值解的 L 1 、 L 2 和 L ∞ 误差及收敛阶

N x

L 1 误差

收敛阶

L 2 误差

收敛阶

L ∞ 误差

收敛阶

k = 1

10

5.99E−02

3.23E−02

4.20E−02

20

1.44E−02

2.06

7.95E−03

2.02

1.26E−02

1.74

40

3.56E−03

2.02

1.97E−03

2.01

3.32E−03

1.92

80

8.82E−04

2.01

4.93E−04

2.00

8.42E−04

1.98

160

2.20E−04

2.00

1.23E−04

2.00

2.07E−04

2.02

k = 2

10

4.43E−03

2.93E−03

8.51E−03

20

4.77E−04

3.21

3.47E−04

3.08

1.24E−03

2.78

40

5.91E−05

3.01

4.33E−05

3.01

1.64E−04

2.91

80

7.36E−06

3.00

5.41E−06

3.00

2.04E−05

3.01

160

9.20E−07

3.00

6.76E−07

3.00

2.48E−06

3.04

3.3. Euler方程的精度测试

我们对一维Euler方程(11)进行精度测试，对此选择了区域 Ω = [ 0 , 2 π ] ，并施加初始条件：

ρ ( x , 0 ) = 1.0 + 0.99 sin ( x ) ,   u ( x , 0 ) = 1.0 ,   p ( x , 0 ) = 1.0 ,

精确解为 ρ ( x , t ) = 1.0 + 0.99 sin ( x − t ) , u ( x , t ) = 1.0 , p ( x , t ) = 1.0 ，实验中同样采用了周期边界条件。我们使用Lax-Wendroff型中心间断伽辽金方法计算了在 t = 1.0 时刻Euler方程中密度 ρ 的数值解，以及其在 L 1 、 L 2 和 L ∞ 范数意义下的误差及收敛阶，结果展示在表3中，从表中可知密度 ρ 在 P k 情况下的收敛阶数仍达到了 k + 1 。

Table 3. The L 1 , L 2 and L ∞ errors and orders of accuracy of the numerical solution for the 1D Euler equation

表3. 一维Euler方程数值解的 L 1 、 L 2 和 L ∞ 误差及收敛阶

N x

L 1 误差

收敛阶

L 2 误差

收敛阶

L ∞ 误差

收敛阶

k = 1

10

1.61E−01

7.48E−02

5.47E−02

20

3.07E−02

2.39

1.49E−02

2.33

1.16E−02

2.24

40

5.33E−03

2.52

2.68E−03

2.49

2.20E−03

2.40

80

1.04E−03

2.36

4.90E−04

2.43

4.12E−04

2.42

160

2.18E−04

2.25

1.07E−04

2.20

8.36E−05

2.30

k = 2

10

2.77E−03

1.49E−03

1.45E−03

20

3.25E−04

3.09

1.84E−04

3.02

1.85E−04

2.96

40

3.97E−05

3.03

2.29E−05

3.00

2.33E−05

2.99

80

4.94E−06

3.01

2.87E−06

3.00

3.06E−06

2.93

160

6.27E−07

2.98

3.65E−07

2.97

5.44E−07

2.49

3.4. 无粘性Burgers方程的间断解测试

在测试了具有连续解的Burgers方程的精度之后，我们还考虑了一维无粘性Burgers方程的初始条件为间断的：

u 0 ( x ) = { 1 , 1 ≤ x ≤ 2 ， 0 , 其 他 ，

计算区域为 Ω = [ 0 , 2 π ] ，左边界条件为inflow边界条件，右边界条件为outflow边界条件。在图1中，我们展示了 t = 1.0 时刻下的数值解与精确解的数据比较，其结果表明，Lax-Wendroff型中心间断伽辽金方法得到的数值解与精确解的图像曲线非常吻合。此外，我们采用了TVB限制器来抑制因高阶方法引入的非物理振荡，以确保数值解在处理间断时的稳定性，同时保留了高分辨率的解结构，从而提升了整体计算的可靠性。

Figure 1. A comparison between the numerical solution and the exact solution at t = 0.1 for N = 160 ，the constant M associated with the TVB limiter is M = 10

图1. N = 160 时在 t = 0.1 下的数值解与精确解的比较图，其中与T限制器有关的常数M的值为 M = 10

3.5. Sod问题和Lax问题

为了评估所提出的Lax-Wendroff型中心间断伽辽金方法在Euler方程间断解上的表现，我们使用一维空间中的两个极端Riemann问题进行测试。对于Sod问题，我们选择区域 Ω = [ 0 , 1 ] ，并施加初始条件：

( ρ ( x , 0 ) , u ( x , 0 ) , p ( x , 0 ) ) = { ( 1 , 0 , 1 ) , 0 ≤ x < 0.5 , ( 0.125 , 0 , 0.1 ) , 0.5 ≤ x ≤ 1.

对于Lax问题，我们同样选择区域 Ω = [ 0 , 1 ] ，并施加初始条件：

( ρ ( x , 0 ) , u ( x , 0 ) , p ( x , 0 ) ) = { ( 0.445 , 0.698 , 3.528 ) , 0 ≤ x < 0.5 , ( 0.5 , 0 , 0.571 ) , 0.5 ≤ x ≤ 1.

气体动力学中由Euler方程控制的Riemann问题的精确解见于[19]，在两个实验中都应用了outflow边界条件。与无粘性Burgers方程的间断解实验类似，我们采用TVB限制器来抑制伪振荡[20]，并在网格 N = 160 上绘制了在 t = 0.2 时刻的Sod问题和在 t = 0.14 时刻的Lax问题的密度数值解，将这些结果与精确解进行比较，如图2所示。其结果表明密度、速度和压力的正性得到了保持，并且Lax-Wendroff型中心间断伽辽金方法下的数值解与精确解非常吻合。

4. 结论

本文设计了求解一维双曲型守恒律方程的Lax-Wendroff型中心间断伽辽金方法，该方法首先运用了Lax-Wendroff型时间离散方法，也就是利用泰勒级数展开处理时间导数，避免了传统的多步时间积分，

(a) (b)

(c) (d)

(e) (f)

Figure 2. The density, velocity and pressure’s exact solution and numerical solution for Sod problem(left) at t = 0.2 and Lax problem(right) at t = 0.14 with N = 160 , the constant M associated with the TVB limiter is M = 2000 (left) for Sod problem and M = 20000 (right) for Lax problem

图2. N = 160 时在 t = 0.2 下的Sod问题(左)和在 t = 0.14 下的Lax问题(右)的密度、速度和压力的精确解和数值解比较图，其中与TVB限制器有关的常数M的值为 M = 2000 (左)和 M = 20000 (右)

然后空间上运用了中心间断伽辽金方法，增加了适用性和准确性。最后，通过对一系列一维双曲型守恒律问题开展数值实验，验证了该方法能够捕捉复杂的波动结构并实现高精度。未来的研究将集中于Lax-Wendroff型中心间断伽辽金方法的稳定性分析和误差估计。

基金项目

重庆市研究生校级科研创新项目(2024S0134)。

参考文献

[1] Bleecker, D. (2018) Basic Partial Differential Equations. Chapman and Hall/CRC.
[2] Shu, C.W. (2009) Discontinuous Galerkin Methods: General Approach and Stability. Numerical Solutions of Partial Differential Equations, 201, 149-201.
[3] Shu, C. and Osher, S. (1988) Efficient Implementation of Essentially Non-Oscillatory Shock-Capturing Schemes. Journal of Computational Physics, 77, 439-471.
https://doi.org/10.1016/0021-9991(88)90177-5
[4] Shu, C. (1988) Total-Variation-Diminishing Time Discretizations. SIAM Journal on Scientific and Statistical Computing, 9, 1073-1084.
https://doi.org/10.1137/0909073
[5] Liu, Y., Shu, C., Tadmor, E. and Zhang, M. (2007) Central Discontinuous Galerkin Methods on Overlapping Cells with a Nonoscillatory Hierarchical Reconstruction. SIAM Journal on Numerical Analysis, 45, 2442-2467.
https://doi.org/10.1137/060666974
[6] Liu, Y., Shu, C., Tadmor, E. and Zhang, M. (2011) Central Local Discontinuous Galerkin Methods on Overlapping Cells for Diffusion Equations. ESAIM: Mathematical Modelling and Numerical Analysis, 45, 1009-1032.
https://doi.org/10.1051/m2an/2011007
[7] Li, F., Xu, L. and Yakovlev, S. (2011) Central Discontinuous Galerkin Methods for Ideal MHD Equations with the Exactly Divergence-Free Magnetic Field. Journal of Computational Physics, 230, 4828-4847.
https://doi.org/10.1016/j.jcp.2011.03.006
[8] Wu, K., Jiang, H. and Shu, C. (2023) Provably Positive Central Discontinuous Galerkin Schemes via Geometric Quasilinearization for Ideal MHD Equations. SIAM Journal on Numerical Analysis, 61, 250-285.
https://doi.org/10.1137/22m1486996
[9] Li, F. and Yakovlev, S. (2009) A Central Discontinuous Galerkin Method for Hamilton-Jacobi Equations. Journal of Scientific Computing, 45, 404-428.
https://doi.org/10.1007/s10915-009-9340-y
[10] Li, M., Guyenne, P., Li, F. and Xu, L. (2014) High Order Well-Balanced CDG-FE Methods for Shallow Water Waves by a Green-Naghdi Model. Journal of Computational Physics, 257, 169-192.
https://doi.org/10.1016/j.jcp.2013.09.050
[11] Li, M. and Chen, A. (2014) High Order Central Discontinuous Galerkin-Finite Element Methods for the Camassa-Holm Equation. Applied Mathematics and Computation, 227, 237-245.
https://doi.org/10.1016/j.amc.2013.11.016
[12] Jiao, M.J., Cheng, Y.D., Liu, Y. and Zhang, M.P. (2020) Central Discontinuous Galerkin Methods for the Generalized Korteweg-De Vries Equation. Communications in Computational Physics, 28, 927-966.
https://doi.org/10.4208/cicp.oa-2019-0099
[13] Liu, Y., Lu, J., Shu, C. and Zhang, M. (2021) Central Discontinuous Galerkin Methods on Overlapping Meshes for Wave Equations. ESAIM: Mathematical Modelling and Numerical Analysis, 55, 329-356.
https://doi.org/10.1051/m2an/2020069
[14] Cockburn, B. and Shu, C. (1998) The Runge-Kutta Discontinuous Galerkin Method for Conservation Laws V: Multidimensional Systems. Journal of Computational Physics, 141, 199-224.
https://doi.org/10.1006/jcph.1998.5892
[15] Ruuth, S.J. and Spiteri, R.J. (2002) Two Barriers on Strong-Stability-Preserving Time Discretization Methods. Journal of Scientific Computing, 17, 211-220.
https://doi.org/10.1023/a:1015156832269
[16] Guo, W., Qiu, J. and Qiu, J. (2014) A New Lax-Wendroff Discontinuous Galerkin Method with Superconvergence. Journal of Scientific Computing, 65, 299-326.
https://doi.org/10.1007/s10915-014-9968-0
[17] Qiu, J., Dumbser, M. and Shu, C. (2005) The Discontinuous Galerkin Method with Lax-Wendroff Type Time Discretizations. Computer Methods in Applied Mechanics and Engineering, 194, 4528-4543.
https://doi.org/10.1016/j.cma.2004.11.007
[18] Xu, Z. and Shu, C. (2022) Third Order Maximum-Principle-Satisfying and Positivity-Preserving Lax-Wendroff Discontinuous Galerkin Methods for Hyperbolic Conservation Laws. Journal of Computational Physics, 470, Article ID: 111591.
https://doi.org/10.1016/j.jcp.2022.111591
[19] Toro, E.F. (2013) Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction. Springer Science & Business Media.
[20] Cockburn, B. and Shu, C. (1989) TVB Runge-Kutta Local Projection Discontinuous Galerkin Finite Element Method for Conservation Laws. II. General Framework. Mathematics of Computation, 52, 411-435.
https://doi.org/10.1090/s0025-5718-1989-0983311-4

为你推荐
Baidu
map