一维Euler方程的无振荡高阶间断Petrov-Galerkin方法

Oscillation-Free High-Order Discontinuous Petrov-Galerkin Method for One-Dimensional Euler Equations

摘要

针对传统高阶数值方法在处理间断问题时易产生非物理振荡的难题，提出一种新型无振荡间断Petrov-Galerkin方法来求解一维Euler方程。采用SSP Runge-Kutta方法进行时间离散，并在每个时间步进后引入一个阻尼算子，通过自适应调节机制抑制数值计算中可能产生的非物理振荡，同时保持间断Petrov-Galerkin方法原有的高精度特性。数值实验结果表明，该方法在光滑区域保持最优收敛阶，在间断附近通过阻尼机制自动调节数值耗散，实现了稳定性和高分辨率的有效平衡。

Abstract

To address the challenge of non-physical oscillations generated by traditional high-order numerical methods when dealing with discontinuous problems， a new oscillation-free discontinuous Petrov-Galerkin method was proposed to solve one-dimensional Euler equations.The method employed the SSP Runge-Kutta scheme for temporal discretization and incorporated a damping operator at every time step， so that the possible non-physical oscillations in numerical calculations could be suppressed through the adaptive regulation mechanism， while maintaining the original high-precision characteristics of the discontinuous Petrov-Galerkin method.The results of numerical experiments show that this method maintains the optimal convergence order in smooth regions and automatically adjusts the numerical dissipation near discontinuities through the damping mechanism， achieving an effective balance between stability and high resolution.

Graphical abstract

关键词

一维Euler方程 / 间断Petrov-Galerkin方法 / SSP Runge-Kutta方法 / 阻尼算子

Key words

one-dimensional Euler equation / discontinuous Petrov-Galerkin method / SSP Runge-Kutta method / damping operator

Euler方程在流体力学中用来刻画无黏流动过程，是典型的非线性守恒律方程组。其非线性对流特征常常可能导致解在有限时间内演变为激波间断，因此，设计稳定且高效的数值计算方法具有很大的挑战性。

数值求解双曲守恒律方程的众多方法中，在有限差分或有限体积方法的基本框架下，基于Harten引理^［1］的TVD（Total variation diminishing）格式在抑制间断附近的数值振荡方面效果显著，但是其内在的极值Clipping效应^［2］导致TVD格式的空间逼近精度很难突破二阶，这限制了TVD格式成为真正的高分辨率格式（高精度且无振荡）。在有限元方法的基本框架下，Cockburn等^［3-5］将有限体积格式处理间断的基本思想与传统DGFEM（Discontinuous Galerkin finite element method）相结合，提出的Runge-Kutta DGFEM具备良好的抑制数值振荡能力和自适应提高逼近精度的灵活性，从而成为一类重要的方法。从有限元离散过程看，DGFEM在局部离散中采用的是标准有限元方法，即试探函数空间和检验函数空间是相同的多项式空间。与此相对应的不同思路是在局部离散中采用非标准有限元过程，即试探函数空间和检验函数空间是不相同的多项式空间，这类有限元常常称作Petrov-Galerkin方法。Bottasso等^［6-7］最早在DG框架下引入Petrov-Galerkin方法，提出采用间断Petrov-Galerkin（DPG）方法求解椭圆问题，并对稳定性和收敛性进行理论分析，将其拓展应用于平流扩散问题的数值模拟。随后，Demkowicz等^［8］提出一种可自动构造检验函数实现优化的DPG方法。Stipcich等^［9］结合DGFEM构建了间断控制体积有限元法（DCVFEM）处理平流扩散问题，从理论上证明该方法具有良好的收敛性和低耗散优势，且能灵活地处理强梯度和间断解问题。陈大伟等^［10］基于Runge-Kutta DG方法和广义差分法提出Runge-Kutta控制体积间断有限元法（RKCVDFEM），首次将DPG方法扩展到一维和二维双曲守恒律的数值求解^［11-14］，该方法能够轻松处理间断解和复杂的边界条件，具有保持物理量局部守恒的特性，并且比Runge-Kutta DG方法更简单。由于DPG方法具有诸多优良特性，近年来已成为偏微分方程数值求解领域中备受关注的方法之一。

与其他高阶方法类似，线性DPG方法也无法避免在解的大梯度或者间断附近产生数值振荡。现有的DPG方法常常利用非线性限制器^［15-16］来控制振荡，它的原理是在获得数值解后，使用各种指标来检测问题单元，然后在问题单元内重构解多项式。但是很多限制器依赖人工调参，这些参数缺乏足够的适应性，如果参数设置不合适，会出现限制效果不足或过度降阶的问题。Lu等^［17］提出一种可以自适应调节振荡的DG方法，直接在半离散DG格式中引入一个阻尼项实现振荡控制。然而在强间断附近，阻尼项的离散会得到一个刚性时间常微分方程组，从而导致非常大的时间计算成本^［18］。Peng等^［19］提出OEDG（Oscillation-eliminating discontinuous Galerkin）方法，它通过求解一个准线性阻尼常微分方程实现自适应调节振荡，且其对应的时间常微分方程不再具有明显刚性，大大降低了时间计算成本。

本文提出一种简单且高效的高分辨率数值方法，该方法在空间离散中采用DPG格式，时间离散采用三阶SSP Runge-Kutta（Strong stability preserving Runge-Kutta）方法，类似于OEDG方法，引入OE型阻尼算子来自动调节非物理振荡。该方法无需复杂的振荡检测，直接在DPG方法计算后统一实施阻尼处理，所构建的阻尼算子基于解的局部梯度跳跃特征，通过自适应耗散机制，调节局部多项式的高阶项系数，实现数值振荡的自动抑制，从而在保持一致高精度的同时，克服了传统限制器常常依赖人工参数的不足。

1 一维Euler方程

一维Euler方程的初值问题为

\left\{\begin{array}{ll}{\mathit{u}}_t+\mathit{f}{(\mathit{u})}_x=\mathrm{0},& (x,t)\in (a,b)\times (\mathrm{0},T]\\ \mathit{u}(x,\mathrm{0})={\mathit{u}}_{\mathrm{0}}(x) ,& x\in (a,b)\end{array}\right.（1）

式中，\mathit{u}={(\rho ,\rho u,E)}^{\mathrm{T}}，\mathit{f}(\mathit{u})=(\rho u,\rho u^{\mathrm{2}}+p,u{(E+p))}^{\mathrm{T}}，\rho为密度，u为速度，\rho u表示动量，E为总能量，p表示压力，总能量与压力的关系为

E=\frac{p}{\gamma -\mathrm{1}}+\frac{\mathrm{1}}{\mathrm{2}}\rho u^{\mathrm{2}}，

其中，无量纲参数\gamma =\mathrm{1.4}，为理想气体的比热比。

引入向量函数\mathit{q}=(q_{\mathrm{1}},q_{\mathrm{2}},q_{\mathrm{3}}{)}^{\mathrm{T}}={(\rho ,\rho u,E)}^{\mathrm{T}}，则流通量向量为

\mathit{f}(\mathit{q})=\left(\begin{array}{c}\rho u\\ \rho u^{\mathrm{2}}+p\\ u(E+p)\end{array}\right)=\left(\begin{array}{c}{q}_{\mathrm{2}}\\ \frac{\mathrm{1}}{\mathrm{2}}(\mathrm{3}-\gamma )\frac{q_{\mathrm{2}}^{\mathrm{2}}}{q_{\mathrm{1}}}+(\gamma -\mathrm{1})q_{\mathrm{3}}\\ \gamma \frac{q_{\mathrm{2}}{q}_{\mathrm{3}}}{q_{\mathrm{1}}}-\frac{\mathrm{1}}{\mathrm{2}}(\gamma -\mathrm{1})\frac{q_{\mathrm{2}}^{\mathrm{3}}}{q_{\mathrm{1}}^{\mathrm{2}}}\end{array}\right)，

它的Jacobian矩阵为

{\mathit{f}}^{\text{'}}\left(\mathit{q}\right)=\left(\begin{array}{ccc}\mathrm{0}& \mathrm{1}& \mathrm{0}\\ \frac{\mathrm{1}}{\mathrm{2}}(\gamma -\mathrm{3}){\left(\frac{q_{\mathrm{2}}}{q_{\mathrm{1}}}\right)}^{\mathrm{2}}& (\mathrm{3}-\gamma )\frac{q_{\mathrm{2}}}{q_{\mathrm{1}}}& \gamma -\mathrm{1}\\ -\gamma \frac{q_{\mathrm{2}}{q}_{\mathrm{3}}}{q_{\mathrm{1}}^{\mathrm{2}}}+(\gamma -\mathrm{1}){\left(\frac{q_{\mathrm{2}}}{q_{\mathrm{1}}}\right)}^{\mathrm{3}}& \gamma \frac{q_{\mathrm{3}}}{q_{\mathrm{1}}}-\frac{\mathrm{3}}{\mathrm{2}}(\gamma -\mathrm{1}){\left(\frac{q_{\mathrm{2}}}{q_{\mathrm{1}}}\right)}^{\mathrm{2}}& \gamma \frac{q_{\mathrm{2}}}{q_{\mathrm{1}}}\end{array}\right)，

Jacobian矩阵的特征值为

{\lambda }_{\mathrm{1}}=u-c,{\lambda }_{\mathrm{2}}=u,{\lambda }_{\mathrm{3}}=u+c,

其中c为声速，满足

c^{\mathrm{2}}=\frac{\gamma p}{\rho }=\gamma (\gamma -\mathrm{1})\left[\frac{q_{\mathrm{3}}}{q_{\mathrm{1}}}-\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{q_{\mathrm{2}}}{q_{\mathrm{1}}}\right)}^{\mathrm{2}}\right]。

2 无振荡间断Petrov-Galerkin方法

2.1 DPG方法空间离散

将求解区间\Omega =\left(a,b\right)剖分成N个小单元，a=x_{\mathrm{1}/\mathrm{2}}\le x_{\mathrm{3}/\mathrm{2}}\le \cdots \le x_{N+\mathrm{1}/\mathrm{2}}=b，记x_{j+\mathrm{1}/\mathrm{2}}(j=\mathrm{0,1},\cdots ,N)为剖分节点，第j个小单元为I_j=(x_{j-\mathrm{1}/\mathrm{2}},x_{j+\mathrm{1}/\mathrm{2}})，单元中点为x_j=\frac{\mathrm{1}}{\mathrm{2}}(x_{j-\mathrm{1}/\mathrm{2}}+x_{j+\mathrm{1}/\mathrm{2}})，单元长度为h_j=x_{j+\mathrm{1}/\mathrm{2}}-x_{j-\mathrm{1}/\mathrm{2}}。

定义有限元空间为

{\mathit{V}}_h^k:=\left\{\mathit{q}={\left(q_{\mathrm{1}},q_{\mathrm{2}},q_{\mathrm{3}}\right)}^{\mathrm{T}}:q_i\in L^{\mathrm{2}}\left(\left[a,b\right]\right)\mathrm{且}{q}_i\left|{}_{I_j}\right.\in P^k\left(I_j\right),i=\mathrm{1,2},\mathrm{3},j=\mathrm{1},\cdots ,N\right\}，

P^k(I_j)是单元I_j上次数不超过k的多项式空间，并且取P^k(I_j)=\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}{\left\{{\phi }_j^l(x)\right\}}_{l=\mathrm{0}}^k，{\phi }_j^l(x)为I_j上的Legendre正交多项式，其中

\begin{array}{l}{\phi }_j^{\mathrm{0}}\left(x\right)=\mathrm{1},{\phi }_j^{\mathrm{1}}\left(x\right)=\frac{x-x_j}{h_j/\mathrm{2}},{\phi }_j^{\mathrm{2}}\left(x\right)={\left(\frac{x-x_j}{h_j/\mathrm{2}}\right)}^{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{3}},\\ {\phi }_j^{\mathrm{3}}\left(x\right)={\left(\frac{x-x_j}{h_j/\mathrm{2}}\right)}^{\mathrm{3}}-\frac{\mathrm{3}}{\mathrm{5}}\left(\frac{x-x_j}{h_j/\mathrm{2}}\right),\end{array}

则数值解{\mathit{u}}_h\in {\mathit{V}}_h^k为

{\mathit{u}}_h\left(x,t\right)={\displaystyle \sum _{l=\mathrm{0}}^k}{\mathit{u}}_j^l\left(t\right){\phi }_j^l\left(x\right),x\in I_j（2）

其中，{\mathit{u}}_j^l(t)为线性组合的系数。

在每个单元I_j上选取x_{j-\mathrm{1}/\mathrm{2}},x_j^{\mathrm{0}},x_j^{\mathrm{1}},\cdots ,x_j^{k-\mathrm{1}},x_{j+\mathrm{1}/\mathrm{2}}作为插值节点，插值节点的取法不是唯一的，要通过对比不同节点选取方案下的数值结果，基于高斯点的数值方法展现出更高分辨率，故本文选取高斯点作为插值节点。选取I_j上的控制体积如图1所示，并选取检验函数空间为

{\mathit{W}}_h:=\left\{\mathit{w}={\left(w_{\mathrm{1}},w_{\mathrm{2}},w_{\mathrm{3}}\right)}^{\mathrm{T}}:w_i\in W_h,i=\mathrm{1,2},\mathrm{3}\right\}，

其中，W_h为分片常数空间，相应基函数为

{\psi }_{\mathit{j}}^{\mathit{l}}=\left\{\begin{array}{cc}\mathrm{1},& \forall x\in I_j^l\\ \mathrm{0},& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}\right.（3）

在方程（1）两边乘以检验函数\mathit{w}\in {\mathit{W}}_h，并在单元I_j上积分，得

{\int }_{I_j}{\mathit{u}}_t\mathit{w}\mathrm{d}x+{\int }_{I_j}\mathit{f}{\left(\mathit{u}\right)}_x\mathit{w}\mathrm{d}x=\mathrm{0}，

用近似解{\mathit{u}}_h代替精确解\mathit{u}，并取\mathit{w}为检验函数空间的基函数{\mathit{\psi }}_j^l，得

{\int }_{I_j}{\left({\mathit{u}}_h\right)}_t{\mathit{\psi }}_j^l\mathrm{d}x+{\int }_{I_j}\mathit{f}{\left({\mathit{u}}_h\right)}_x{\mathit{\psi }}_j^l\mathrm{d}x=\mathrm{0},l=\mathrm{0,1},\cdots ,k，

令{\mathit{u}}_h=(q_{\mathrm{1}h},q_{\mathrm{2}h},q_{\mathrm{3}h}{)}^{\mathrm{T}}，\mathit{f}({\mathit{u}}_h)=(f_{\mathrm{1}}({\mathit{u}}_h),f_{\mathrm{2}}({\mathit{u}}_h),f_{\mathrm{3}}({\mathit{u}}_h{))}^{\mathrm{T}}，则上式展开得

\left\{\begin{array}{ll}{\int }_{I_j}{\left(q_{\mathrm{1}h}\right)}_t{\psi }_j^l\mathrm{d}x+{\int }_{I_j}{f}_{\mathrm{1}}{\left({\mathit{u}}_h\right)}_x{\psi }_j^l\mathrm{d}x=\mathrm{0},& l=\mathrm{0,1},\cdots ,k\\ {\int }_{I_j}{\left(q_{\mathrm{2}h}\right)}_t{\psi }_j^l\mathrm{d}x+{\int }_{I_j}{f}_{\mathrm{2}}{\left({\mathit{u}}_h\right)}_x{\psi }_j^l\mathrm{d}x=\mathrm{0}, & l=\mathrm{0,1},\cdots ,k\\ {\int }_{I_j}{\left(q_{\mathrm{3}h}\right)}_t{\psi }_j^l\mathrm{d}x+{\int }_{I_j}{f}_{\mathrm{3}}{\left({\mathit{u}}_h\right)}_x{\psi }_j^l\mathrm{d}x=\mathrm{0}, & l=\mathrm{0,1},\cdots ,k\end{array}\right.，

将式（3）代入上式，可得

{\int }_{I_j^l}{\left({\mathit{u}}_h\right)}_t\mathrm{d}x+{\int }_{I_j^l}\mathit{f}{\left({\mathit{u}}_h\right)}_x\mathrm{d}x=\mathrm{0}, l=\mathrm{0,1},\cdots ,k，

分部积分得

\frac{\mathrm{d}}{\mathrm{d}t}\left({\int }_{I_j^l}{\mathit{u}}_h\mathrm{d}x\right)=-\mathit{f}\left({\mathit{u}}_h\right)\left|{}_{I_j^l}\right., l=\mathrm{0,1},\cdots ,k，

将式（2）代入上式得

\left\{\begin{array}{c}\frac{\mathrm{d}}{\mathrm{d}t}\left[\left({\int }_{I_j^{\mathrm{0}}}{\phi }_j^{\mathrm{0}}\mathrm{d}x\right){\mathit{u}}_j^{\mathrm{0}}+\left({\int }_{I_j^{\mathrm{0}}}{\phi }_j^{\mathrm{1}}\mathrm{d}x\right){\mathit{u}}_j^{\mathrm{1}}+\cdots +\left({\int }_{I_j^{\mathrm{0}}}{\phi }_j^k\mathrm{d}x\right){\mathit{u}}_j^k\right]=-\mathit{f}\left({\mathit{u}}_h\right)\left|{}_{I_j^{\mathrm{0}}}\right.\\ \frac{\mathrm{d}}{\mathrm{d}t}\left[\left({\int }_{I_j^{\mathrm{1}}}{\phi }_j^{\mathrm{0}}\mathrm{d}x\right){\mathit{u}}_j^{\mathrm{0}}+\left({\int }_{I_j^{\mathrm{1}}}{\phi }_j^{\mathrm{1}}\mathrm{d}x\right){\mathit{u}}_j^{\mathrm{1}}+\cdots +\left({\int }_{I_j^{\mathrm{1}}}{\phi }_j^k\mathrm{d}x\right){\mathit{u}}_j^k\right]=-\mathit{f}\left({\mathit{u}}_h\right)\left|{}_{I_j^{\mathrm{1}}}\right.\\ ⋮\\ \frac{\mathrm{d}}{\mathrm{d}t}\left[\left({\int }_{I_j^k}{\phi }_j^{\mathrm{0}}\mathrm{d}x\right){\mathit{u}}_j^{\mathrm{0}}+\left({\int }_{I_j^k}{\phi }_j^{\mathrm{1}}\mathrm{d}x\right){\mathit{u}}_j^{\mathrm{1}}+\cdots +\left({\int }_{I_j^k}{\phi }_j^k\mathrm{d}x\right){\mathit{u}}_j^k\right]=-\mathit{f}\left({\mathit{u}}_h\right)\left|{}_{I_j^k}\right.\end{array}\right.，

整理得

\frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{array}{c}{\mathit{u}}_j^{\mathrm{0}}\\ {\mathit{u}}_j^{\mathrm{1}}\\ ⋮\\ {\mathit{u}}_j^{\mathrm{k}}\end{array}\right)={\mathit{A}}^{-\mathrm{1}}\left(\begin{array}{c}\mathit{f}\left({\mathit{u}}_h\left(x_{j-\mathrm{1}/\mathrm{2}},t\right)\right)-\mathit{f}\left({\mathit{u}}_h\left(x_j^{\mathrm{0}},t\right)\right)\\ \mathit{f}\left({\mathit{u}}_h\left(x_j^{\mathrm{0}},t\right)\right)-\mathit{f}\left({\mathit{u}}_h\left(x_j^{\mathrm{1}},t\right)\right)\\ ⋮\\ \mathit{f}\left({\mathit{u}}_h\left(x_j^{k-\mathrm{1}},t\right)\right)-\mathit{f}\left({\mathit{u}}_h\left(x_{j+\mathrm{1}/\mathrm{2}},t\right)\right)\end{array}\right)（4）

其中

\mathit{A}=\left(\begin{array}{cccc}{\int }_{I_j^{\mathrm{0}}}{\phi }_j^{\mathrm{0}}\mathrm{d}x& {\int }_{I_j^{\mathrm{0}}}{\phi }_j^{\mathrm{1}}\mathrm{d}x& \cdots & {\int }_{I_j^{\mathrm{0}}}{\phi }_j^k\mathrm{d}x\\ {\int }_{I_j^{\mathrm{1}}}{\phi }_j^{\mathrm{0}}\mathrm{d}x& {\int }_{I_j^{\mathrm{1}}}{\phi }_j^{\mathrm{1}}\mathrm{d}x& \cdots & {\int }_{I_j^{\mathrm{1}}}{\phi }_j^k\mathrm{d}x\\ ⋮& ⋮& \ddots & ⋮\\ {\int }_{I_j^k}{\phi }_j^{\mathrm{0}}\mathrm{d}x& {\int }_{I_j^k}{\phi }_j^{\mathrm{1}}\mathrm{d}x& \cdots & {\int }_{I_j^k}{\phi }_j^k\mathrm{d}x\end{array}\right)。

数值解在单元交界处允许间断，由于\mathit{f}({\mathit{u}}_h(x_{j\pm \mathrm{1}/\mathrm{2}},t))无定义，故用相容单调的Lipschitz数值流通量{\widehat{\mathit{f}}}_{j\pm \mathrm{1}/\mathrm{2}}=\widehat{\mathit{f}}({\mathit{u}}_{j\pm \mathrm{1}/\mathrm{2}}^{-},{\mathit{u}}_{j\pm \mathrm{1}/\mathrm{2}}^{+})来代替\mathit{f}({\mathit{u}}_h(x_{j\pm \mathrm{1}/\mathrm{2}},t))，这里采用Lax-Friedrichs数值流通量^［20］

{\widehat{\mathit{f}}}_{j\pm \mathrm{1}/\mathrm{2}}=\frac{\mathrm{1}}{\mathrm{2}}\left[\mathit{f}\left({\mathit{u}}_{j\pm \mathrm{1}/\mathrm{2}}^{+}\right)+\mathit{f}\left({\mathit{u}}_{j\pm \mathrm{1}/\mathrm{2}}^{-}\right)-{\alpha }_{j\pm \mathrm{1}/\mathrm{2}}\left({\mathit{u}}_{j\pm \mathrm{1}/\mathrm{2}}^{+}-{\mathit{u}}_{j\pm \mathrm{1}/\mathrm{2}}^{-}\right)\right]，

其中{\mathit{u}}_{j+\mathrm{1}/\mathrm{2}}^{\pm }为界面值的左右极限，{\alpha }_{j\pm \mathrm{1}/\mathrm{2}}表示\mathit{f}的Jacobian矩阵的局部最大波速，具体计算公式为{\alpha }_{j\pm \mathrm{1}/\mathrm{2}}=\mathrm{m}\mathrm{a}\mathrm{x}\left\{|u|-|c|,|u|,|u|+|c|\right\}。

至此，间断Petrov-Galerkin格式空间离散完毕，方程组（4）可简记为

\frac{\mathrm{d}{\mathit{u}}_h}{\mathrm{d}t}=L_h\left({\mathit{u}}_h,t\right)。

2.2 OEDPG方法完全离散

DPG方法在间断附近可能会产生数值振荡，如果不控制这些振荡，会导致严重的数值不稳定和模拟崩溃。因此，在每个Runge-Kutta阶段后设计OE型阻尼算子抑制数值振荡，从而提高解的稳定性和准确性。当使用三阶SSP Runge-Kutta方法进行时间步进时，得到OEDPG方法从t^n到t^{n+\mathrm{1}}的实现过程为

\left\{\begin{array}{ll}{\mathit{u}}_h^{n,\mathrm{1}}={\mathit{u}}_{\sigma }^n+\mathrm{\Delta }t{L}_h\left({\mathit{u}}_{\sigma }^n\right),& {\mathit{u}}_{\sigma }^{n,\mathrm{1}}={\mathrm{ℱ}}_{\mathrm{\Delta }t}\left({\mathit{u}}_h^{n,\mathrm{1}}\right)\\ {\mathit{u}}_h^{n,\mathrm{2}}=\frac{\mathrm{3}}{\mathrm{4}}{\mathit{u}}_{\sigma }^n+\frac{\mathrm{1}}{\mathrm{4}}{\mathit{u}}_{\sigma }^{n,\mathrm{1}}+\frac{\mathrm{1}}{\mathrm{4}}\mathrm{\Delta }t{L}_h\left({\mathit{u}}_{\sigma }^{n,\mathrm{1}}\right),& {\mathit{u}}_{\sigma }^{n,\mathrm{2}}={\mathrm{ℱ}}_{\mathrm{\Delta }t}\left({\mathit{u}}_h^{n,\mathrm{2}}\right)\\ {\mathit{u}}_h^{n+\mathrm{1}}=\frac{\mathrm{1}}{\mathrm{3}}{\mathit{u}}_{\sigma }^n+\frac{\mathrm{2}}{\mathrm{3}}{\mathit{u}}_{\sigma }^{n,\mathrm{2}}+\frac{\mathrm{2}}{\mathrm{3}}\mathrm{\Delta }t{L}_h\left({\mathit{u}}_{\sigma }^{n,\mathrm{2}}\right),& {\mathit{u}}_{\sigma }^{n+\mathrm{1}}={\mathrm{ℱ}}_{\mathrm{\Delta }t}\left({\mathit{u}}_h^{n+\mathrm{1}}\right)\end{array}\right.，

其中，{\mathit{u}}_{\sigma }^n={\mathit{u}}_h^n={\mathit{u}}_h(x,t^n)是t^n时间层的数值解，\mathrm{\Delta }t=t^{n+\mathrm{1}}-t^n是时间步长，{\mathrm{ℱ}}_{\mathrm{\Delta }t}是从{\mathit{V}}_h^k到{\mathit{V}}_h^k的振荡消除算子，定义{\mathrm{ℱ}}_{\mathrm{\Delta }t}({\mathit{u}}_h^{n,l})={\mathit{u}}_{\sigma }(x,\mathrm{\Delta }t)，其中{\mathit{u}}_{\sigma }(x,\widehat{t})\in {\mathit{V}}_h^k是下列初值问题的解：

\left\{\begin{array}{l}\frac{\mathrm{d}}{\mathrm{d}\widehat{t}}{\int }_{I_j}{\mathit{u}}_{\sigma }\mathit{v}\mathrm{d}x+{\displaystyle \sum _{m=\mathrm{0}}^k}{\beta }_j\frac{{\sigma }_j^m\left({\mathit{u}}_h^{n,l}\right)}{h_j}{\int }_{I_j}\left({\mathit{u}}_{\sigma }-P_h^{m-\mathrm{1}}{\mathit{u}}_{\sigma }\right)\mathit{v}\mathrm{d}x=\mathrm{0}, \mathit{v}\in {\mathit{V}}_h^k\\ {\mathit{u}}_{\sigma }\left(x,\mathrm{0}\right)={\mathit{u}}_h^{n,l}\left(x\right)\end{array}\right.（5）

其中，\widehat{t}是不同于t的伪时间，{\beta }_j是Jacobian矩阵\frac{\partial f}{\partial u}({\overline{\mathit{u}}}_j^{n,l})的谱半径，{\overline{\mathit{u}}}_j^{n,l}是{\mathit{u}}_j^{n,l}在I_j上的单元积分平均值，算子P_h^m是{\mathit{V}}_h^k到{\mathit{V}}_h^m的标准L^{\mathrm{2}}投影，对任意函数\mathit{v}\in {\mathit{V}}_h^k，有P_h^m\mathit{v}\in {\mathit{V}}_h^m且满足

{\int }_{I_j}\left(P_h^m\mathit{v}-\mathit{v}\right)\mathit{w}\mathrm{d}x=\mathrm{0},\forall \mathit{w}\in {\mathit{V}}_h^m,

并且定义P_h^{-\mathrm{1}}=P_h^{\mathrm{0}}，特别地，在I_j上P_h^{\mathrm{0}}\mathit{v}={\overline{\mathit{v}}}_j。

定义阻尼系数为

{\sigma }_j^m\left({\mathit{u}}_h\right)=\underset{\mathrm{1}\le i\le Q}{\mathrm{m}\mathrm{a}\mathrm{x}}{\sigma }_j^m\left(u_h^{\left(i\right)}\right)，

其中，u_h^{\left(i\right)}是{\mathit{u}}_h的第i个分量，Q=\mathrm{3}，且{\sigma }_j^m\left(u_h^{(i)}\right)定义为

{\sigma }_j^m\left(u_h^{(i)}\right)=\left\{\begin{array}{ll}\mathrm{0},& u_h^{(i)}\equiv \mathrm{a}\mathrm{v}\mathrm{g}\left(u_h^{(i)}\right)\\ \frac{(\mathrm{2}m+\mathrm{1})h_j^m}{(\mathrm{2}k-\mathrm{1})m!}\frac{\left|{⟦{\partial }_x^m{u}_h^{(i)}⟧}_{j-\mathrm{1}/\mathrm{2}}\right|+\left|{⟦{\partial }_x^m{u}_h^{(i)}⟧}_{j+\mathrm{1}/\mathrm{2}}\right|}{\mathrm{2}{‖u_h^{(i)}-\mathrm{a}\mathrm{v}\mathrm{g}\left(u_h^{(i)}\right)‖}_{L^{\mathrm{\infty }}(\Omega )}},& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}\right.，

其中，{⟦{\partial }_x^m{u}_h^{\left(i\right)}⟧}_{j\pm \mathrm{1}/\mathrm{2}}={\partial }_x^m{u}_h^{\left(i\right)}(x_{j\pm \mathrm{1}/\mathrm{2}}^{+})-{\partial }_x^m{u}_h^{\left(i\right)}(x_{j\pm \mathrm{1}/\mathrm{2}}^{-})表示函数{\partial }_x^m{u}_h^{\left(i\right)}在x_{j\pm \mathrm{1}/\mathrm{2}}处的跳跃，\mathrm{a}\mathrm{v}\mathrm{g}(u_h^{(i)})=\frac{\mathrm{1}}{\left|\Omega \right|}{\int }_{\Omega }{u}_h^{(i)}(x)\mathrm{d}x表示u_h^{\left(i\right)}在\Omega上的平均值。

初值问题（5）是线性的，可以精确求解，设问题的解为

{\mathit{u}}_{\sigma }\left(x,\widehat{t}\right)={\displaystyle \sum _{l=\mathrm{0}}^k}{\mathit{u}}_j^l\left(\widehat{t}\right){\phi }_j^l\left(x\right), x\in I_j（6）

根据算子P_h^m的定义，有

\left({\mathit{u}}_{\sigma }-P_h^{m-\mathrm{1}}{\mathit{u}}_{\sigma }\right)\left(x,\widehat{t}\right)={\displaystyle \sum _{l=\mathrm{m}\mathrm{a}\mathrm{x}\{\mathrm{1},m\}}^k}{\mathit{u}}_j^l\left(\widehat{t}\right){\phi }_j^l\left(x\right)（7）

将式（6）和式（7）代入式（5），并取\mathit{v}={\mathit{\phi }}_j^l，得

\frac{\mathrm{d}}{\mathrm{d}\widehat{t}}{\mathit{u}}_j^l\left(\widehat{t}\right){\int }_{I_j}{\left({\phi }_j^l\left(x\right)\right)}^{\mathrm{2}}\mathrm{d}x+{\displaystyle \sum _{m=\mathrm{0}}^l}{\beta }_j\frac{{\sigma }_j^m}{h_j}{\mathit{u}}_j^l\left(\widehat{t}\right){\int }_{I_j}{\left({\phi }_j^l\left(x\right)\right)}^{\mathrm{2}}\mathrm{d}x=\mathrm{0}, l=\mathrm{1,2},\cdots ,k（8）

化简得

\frac{\mathrm{d}}{\mathrm{d}\widehat{t}}{\mathit{u}}_j^l\left(\widehat{t}\right)+\left(\frac{{\beta }_j}{h_j}{\displaystyle \sum _{m=\mathrm{0}}^l}{\sigma }_j^m\right){\mathit{u}}_j^l\left(\widehat{t}\right)=\mathrm{0}, l=\mathrm{1,2},\cdots ,k\mathrm{。}

从0到\widehat{t}进行积分，并取\widehat{t}=\mathrm{\Delta }t得

{\mathit{u}}_j^l\left(\mathrm{\Delta }t\right)={\mathit{u}}_j^l\left(\mathrm{0}\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{{\beta }_j\mathrm{\Delta }t}{h_j}{\displaystyle \sum _{m=\mathrm{0}}^l}{\sigma }_j^m\right), l=\mathrm{1,2},\cdots ,k\mathrm{。}

当l=\mathrm{0}时，{\int }_{I_j}\left({\mathit{u}}_{\sigma }-P_h^{m-\mathrm{1}}{\mathit{u}}_{\sigma }\right){\phi }_j^{\mathrm{0}}\left(x\right)\mathrm{d}x\equiv \mathrm{0}，式（8）化为

\frac{\mathrm{d}}{\mathrm{d}\widehat{t}}{\mathit{u}}_j^{\mathrm{0}}\left(\widehat{t}\right){\int }_{I_j}{\left({\phi }_j^{\mathrm{0}}\left(x\right)\right)}^{\mathrm{2}}\mathrm{d}x=\mathrm{0}，

解得

{\mathit{u}}_j^{\mathrm{0}}\left(\mathrm{\Delta }t\right)={\mathit{u}}_j^{\mathrm{0}}\left(\mathrm{0}\right)。

综上，

{\mathit{u}}_{\sigma }^{n,l}={\mathrm{ℱ}}_{\mathrm{\Delta }t}\left({\mathit{u}}_h^{n,l}\right)={\mathit{u}}_{\sigma }\left(x,\widehat{t}\right)\left|{}_{\widehat{t}=\mathrm{\Delta }t}\right.={\mathit{u}}_j^{\mathrm{0}}\left(\mathrm{0}\right){\phi }_j^{\mathrm{0}}\left(x\right)+{\displaystyle \sum _{l=\mathrm{1}}^k}{\mathit{u}}_j^l\left(\mathrm{0}\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{{\beta }_j\mathrm{\Delta }t}{h_j}{\displaystyle \sum _{m=\mathrm{0}}^l}{\sigma }_j^m\right){\phi }_j^l\left(x\right)，

{\mathit{u}}_j^l\left(\mathrm{0}\right)=\frac{{\int }_{I_j}{\mathit{u}}_h^{n,l}{\phi }_j^l\left(x\right)\mathrm{d}x}{{‖{\phi }_j^l\left(x\right)‖}_{L^{\mathrm{2}}\left(I_j\right)}^{\mathrm{2}}}。

阻尼系数基于解在界面处的导数跳跃设计，在光滑区域，导数跳跃小，阻尼项自动弱化，保持高阶多项式的高精度逼近；在间断区域，导数跳跃大，阻尼项增强，自适应耗散高阶多项式系数，从而构建起振荡消除算子{\mathrm{ℱ}}_{\mathrm{\Delta }t}的自适应振荡抑制机制。

3 数值算例

本节通过5个数值算例证明所提算法的有效性和稳定性。空间网格均匀剖分，对于光滑问题，采用k+\mathrm{1}阶SSP Runge-Kutta时间离散来验证基于P^k元的OEDPG方法的最优收敛阶；对于间断问题，使用三阶SSP Runge-Kutta时间离散。CFL条件数设为C_{\mathrm{C}\mathrm{F}\mathrm{L}}=\frac{\mathrm{0.95}}{\mathrm{2}k+\mathrm{1}}^［19］，k是有限元空间{\mathit{V}}_h^k上多项式的最高次数，时间步长\mathrm{\Delta }t=C_{\mathrm{C}\mathrm{F}\mathrm{L}}\frac{h}{\beta }，其中\beta为局部最大波速。

3.1 精度测试

算例1 对于一维Euler方程（1），首先对一个具有周期边界条件的光滑算例进行精度测试，给定初始条件为

\left(\rho ,u,p\right)\left(x,\mathrm{0}\right)=\left(\mathrm{1}+\mathrm{0.2}\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{2}\mathrm{\pi }x),\mathrm{1,1}\right)，

精确解为

\left(\rho ,u,p\right)\left(x,t\right)=\left(\mathrm{1}+\mathrm{0.2}\mathrm{s}\mathrm{i}\mathrm{n}[\mathrm{2}\mathrm{\pi }(x-t)],\mathrm{1,1}\right),

求解区域为[\mathrm{0,1}]，最终计算时刻t=\mathrm{1}，表1给出了OEDPG方法分别在P^{\mathrm{1}}元、P^{\mathrm{2}}元和P^{\mathrm{3}}元下密度的L^{\mathrm{1}}、L^{\mathrm{2}}误差及收敛阶^［21］。注意，在粗网格上，高阶阻尼效应在数值误差中占主导地位，导致OEDPG方法计算得到的收敛阶数高于k+\mathrm{1}，表1的结果证明了基于P^k元的OEDPG方法在两种不同范数下达到k+\mathrm{1}阶的最优收敛阶数。表2是在相同条件下没有OE过程的DPG方法的误差及收敛阶。

3.2 无振荡测试

算例2 Lax激波管问题的初始条件为

\left(\rho ,u,p\right)\left(x,\mathrm{0}\right)=\left\{\begin{array}{ll}\left(\mathrm{0.445,0.698,3.528}\right),& x\le \mathrm{0}\\ \left(\mathrm{0.5,0},\mathrm{0.571}\right),& x>\mathrm{0}\end{array}\right.，

计算区间为x\in [-\mathrm{1,1}]，最终计算时刻t=\mathrm{0.28}，采用无反射边界条件以避免波在边界反射。图2是基于P^k元的OEDPG方法在空间剖分为200份时密度\rho的计算结果，从图2中可以看出在激波间断处没有明显振荡，提出的新方法能够很好地模拟准确解。

算例3 Sod激波管问题描述了一个在封闭管道中理想气体的运动，它的初始条件为

\left(\rho ,u,p\right)\left(x,\mathrm{0}\right)=\left\{\begin{array}{ll}\left(\mathrm{1,0},\mathrm{1}\right),& x\le \mathrm{0}\\ \left(\mathrm{0.125,0},\mathrm{0.1}\right),& x>\mathrm{0}\end{array}\right.，

计算区间为x\in [-\mathrm{1,1}]，最终计算时刻t=\mathrm{0.28}，采用无反射边界条件，图3是空间剖分为200份时密度ρ的计算结果。从图3可以看出OEDPG方法能够有效消除数值振荡。

算例4 Shu-Osher问题描述了正弦波与激波的相互作用，这是一个具有复杂结构的基准测试，可用于检验数值方法处理激波、接触间断和复杂波系的能力。初始条件为

，

计算区域为[-\mathrm{5,5}]，最终计算时刻t=\mathrm{1.8}，采用无反射边界条件，图4显示了基于P^k元的OEDPG方法在空间剖分为400份时密度ρ的计算结果，可以发现P^{\mathrm{2}}元和P^{\mathrm{3}}元的数值解明显优于P^{\mathrm{1}}元的数值解，表明高阶数值方案能够更好地解析波。

算例5 Woodward-Colella blast wave问题描述了两个冲击波的相互作用，初始条件是

，

计算区域为x\in [\mathrm{0,1}]，在x=\mathrm{0}和x=\mathrm{1}处采用反射边界条件，图5给出了在400个剖分网格下t=\mathrm{0.038}时密度ρ的数值结果。图5再次表明新方法在间断附近没有虚假振荡，在固定自由度的情况下，P^{\mathrm{3}}元的数值解更接近准确解，表明高阶方法具有更高的分辨率。

4 结论

针对一维Euler方程，通过在DPG格式中引入OE型阻尼算子，构建了一种具备自适应振荡调节机制的高分辨率数值方法。光滑初值算例的计算结果表明，该方法可以保持DPG格式的高阶精度；经典间断解算例的计算结果表明，该方法能够有效抑制非物理振荡。在相同网格条件下，随着格式精度阶的提升，新方法对激波间断的捕捉更精准，分辨率和计算精度显著提高，数值耗散降低。另外，高阶格式在网格加密时仍能保持良好的计算效率。

参考文献

[1] SHU C W.Total-variation-diminishing time discretizations［J］.SIAM Journal on Scientific and Statistical Computing，1988，9（6）：1073-1084.

[2] ZHANG X X， SHU C W.A genuinely high order total variation diminishing scheme for one-dimensional scalar conservation laws［J］.SIAM Journal on Numerical Analysis，2010，48（2）：772-795.

[3] COCKBURN B， HOU S， SHU C W.The Runge-Kutta local projection discontinuous galerkin finite element method for conservation laws Ⅳ：The multidimensional case［J］.Mathematics of Computation，1990，54（190）：545-581.

[4] COCKBURN B， LIN S Y， SHU C W.TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws Ⅲ：One-dimensional systems［J］.Journal of Computational Physics，1989，84（1）：90-113.

[5] COCKBURN B， SHU C W.TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws Ⅱ：General framework［J］.Mathematics of Computation，1989，52（186）：411-435.

[6] BOTTASSO C L， MICHELETTI S， SACCO R.The discontinuous Petrov-Galerkin method for elliptic problems［J］.Computer Methods in Applied Mechanics and Engineering，2002，191（31）：3391-3409.

[7] BOTTASSO C L， MICHELETTI S， SACCO R.A multiscale formulation of the Discontinuous Petrov-Galerkin method for advective-diffusive problems［J］.Computer Methods in Applied Mechanics and Engineering，2005，194（25/26）：2819-2838.

[8] DEMKOWICZ L， GOPALAKRISHNAN J.A class of discontinuous Petrov-Galerkin methods.Ⅱ.Optimal test functions［J］.Numerical Methods for Partial Differential Equations，2011，27（1）：70-105.

[9] STIPCICH G， PILLER M， PIVETTA M，et al.Discontinuous control-volume/finite-element method for advection-diffusion problems［J］.Computers & Fluids，2011，52：33-49.

[10] 陈大伟，蔚喜军.一维双曲守恒律的龙格-库塔控制体积间断有限元方法［J］.计算物理，2009，26（4）：501-509.

[11] CHEN D W， YU X J， CHEN Z X.The Runge-Kutta control volume discontinuous finite element method for systems of hyperbolic conservation laws［J］.International Journal for Numerical Methods in Fluids，2011，67（6）：771-786.

[12] ZHAO G Z， YU X J， GUO P Y.The discontinuous Petrov-Galerkin method for one-dimensional compressible Euler equations in the Lagrangian coordinate［J］.Chinese Physics B，2013，22（5）：96-103.

[13] ZHAO G Z， YU X J， LI Z Z.Runge-Kutta control volume discontinuous finite element method for multi-medium fluid simulations［J］.计算物理，2014，31（3）：271-284.

[14] ZHAO G Z， YU X J， GUO H M.A discontinuous Petrov-Galerkin method for two-dimensional compressible gas dynamic equations in Lagrangian coordinates［J］.计算物理，2017，34（3）：294-308.

[15] QIU J X， SHU C W.Runge-Kutta discontinuous Galerkin method using WENO limiters［J］.SIAM Journal on Scientific Computing，2005，26（3）：907-929.

[16] ZHONG X H， SHU C W.A simple weighted essentially nonoscillatory limiter for Runge-Kutta discontinuous Galerkin methods［J］.Journal of Computational Physics，2013，232（1）：397-415.

[17] LU J F， LIU Y， SHU C W.An oscillation-free discontinuous Galerkin method for scalar hyperbolic conservation laws［J］.SIAM Journal on Numerical Analysis，2021，59（3）：1299-1324.

[18] LIU Y， LU J F， SHU C W.An essentially oscillation-free discontinuous Galerkin method for hyperbolic systems［J］.SIAM Journal on Scientific Computing，2022，44（1）：A230-A259.

[19] PENG M T， SUN Z， WU K L.OEDG：Oscillation-eliminating discontinuous Galerkin method for hyperbolic conservation laws［J］.Mathematics of Computation，2025，94（353）：1147-1198.

[20] QIU J X， KHOO B C， SHU C W.A numerical study for the performance of the Runge-Kutta discontinuous Galerkin method based on different numerical fluxes［J］.Journal of Computational Physics，2006，212（2）：540-565.

[21] 高巍，李德茂.非线性奇异椭圆问题的有限元误差分析［J］.内蒙古大学学报（自然科学版），2004，35（1）：16-19.

基金资助

内蒙古自治区人才开发基金项目(12000-1300020240)