1. 引言

双曲方程作为偏微分方程中的一类重要方程，其解中往往出现间断，因此对于数值方法提出很大挑战。针对双曲守恒律方程，一系列的格式被提出和发展，例如Godunov格式 [1]，Lax-Friedrichs格式 [2]，总变差递减(Total Variation Diminishing, TVD) [3] [4]，总变差有界(Total Variation Bounded, TVB) [5]，加权本质无震荡格式(Weighted Essentially Non-Oscillatory, WENO) [6] 等。

间断有限体积元方法最早在2004年 [7] 中被YeXiu提出，作者将有限体积元方法与DG方法结合，针对椭圆问题提出了间断有限体积元方法。之后，毕春加将此方法推广到了抛物问题及双曲问题 [8] [9]。该方法不仅构造简单，还具备保持物理量间局部守恒、可以处理复杂边界、可并行计算等优点，发展前景广阔。在2009年，陈大伟等借鉴了控制体积有限元方法与龙格–库塔间断有限元方法，提出龙格-库塔控制体积间断有限元方法(Runge-Kutta Control Volume Discontinuous Finite Element Method, RKCVDFEM)，保证了数值解的局部守恒性且可以很好的处理激波间断问题。但RKCVDFEM使用了minmod限制器，由于该限制器中参数M的选取依赖于不同的问题，从而限制了方法的适用范围。

区别于传统的间断有限体积元方法 [10] [11] [12] [13] [14]，本文基于SSP Runge-Kutta时间离散方法，采用MPS (Maximum-principle-satisfying)的高阶精度限制器，构造了保极值间断有限体积元方法来求解双曲守恒律方程，并通过经典数值算例来验证方法的可行性与高效性。

2. 格式的构造

2.1. 空间离散

一维守恒律方程

$\{\begin{array}{l}{u}_t+f{\left(u\right)}_x=0\\ u\left(x,0\right)=u_0\end{array}$

这里的f是流通量函数。剖分定义域，令 $x_{i-1/2}$ 为剖分节点，第i个单元记为 $I_i=\left(x_{i-1/2},x_{i+1/2}\right)$，单元长度为 $h=x_{i+1/2}-x_{i-1/2}$。在每个控制单元 $I_i$ 上依次选取3个插值节点，分别记为 $x_{i-1/2},x_i,x_{i+1/2}$，围绕插值节点选取控制体积,三个控制单元分别记为 $I_{i-1/2}^{*-}=\left(x_{i-1/2},x_{i-1/4}\right)$， $I_i^{*}=\left(x_{i-1/4},x_{i+1/4}\right)$， $I_{i+1/2}^{*+}=\left(x_{i+1/4},x_{i+1/2}\right)$，试探函数空间 $U_h$ 选为线性元空间，相应的基函数选取为Lagrange基函数。选取的基函数记为 ${\phi }_j^{\left(l\right)}$，则单元形状函数 $u_h$ 为

$u_h\left(x,t\right)={\displaystyle \underset{l=1}{\overset{3}{\sum }}{u}_j^{\left(l\right)}{\phi }_j^{\left(l\right)}(\; x\; )}$

选取检验函数为控制单元 $I^{\ast }$ 上的分片常数函数，相应的基函数

${\psi }_j^l=\{\begin{array}{l}1,\forall x\in I^{\ast }\hfill \\ 0,其它\hfill \end{array},l=1,2,3$

从原方程出发进行空间离散。先在方程两边同时乘以检验函数v，并在一个单元 $I_i$ 上做积分，u用 $u_h$ 来代替

$\{\begin{array}{l}\frac{\text{d}}{\text{d}t}\left({\displaystyle {\int }_{I_{i-1/2}^{*+}}{\phi }_j^{\left(1\right)}\text{d}x*u_j^{\left(1\right)}}+{\displaystyle {\int }_{I_{i-1/2}^{*+}}{\phi }_j^{\left(2\right)}\text{d}x*u_j^{\left(2\right)}}+0\right)=f_{i-1/2}-f_{i-1/4}\hfill \\ \frac{\text{d}}{\text{d}t}\left({\displaystyle {\int }_{I_i^{*}}{\phi }_{{}_j}^{\left(1\right)}\text{d}x*u_j^{\left(1\right)}}+{\displaystyle {\int }_{I_i^{*}}{\phi }_j^{\left(2\right)}\text{d}x*u_j^{\left(2\right)}}+{\displaystyle {\int }_{I_i^{*}}{\phi }_j^{\left(3\right)}\text{d}x*u_j^{\left(3\right)}}\right)=f_{i-1/4}-f_{i+1/4}\hfill \\ \frac{\text{d}}{\text{d}t}\left(0+{\displaystyle {\int }_{I_{i+1/2}^{*-}}{\phi }_j^{\left(2\right)}\text{d}x*u_{{}^j}^{\left(2\right)}}+{\displaystyle {\int }_{I_{i+1/2}^{*-}}{\phi }_j^{\left(3\right)}\text{d}x*u_j^{\left(3\right)}}\right)=f_{i+1/4}-f_{i+1/2}\hfill \end{array}$

其中 $f_{i\pm 1/4}=f\left(u_h\left(x_{i\pm 1/4},t\right)\right)$，即 $f_{i\pm 1/4}$ 就是把 $x_{i\pm 1/4}$ 带入单元形状函数 $u_h$ 中求得，而又由于 $f_{i\pm 1/2}$ 无定义,故本文采用Lax-Friedrichs数值流通量来定义 $f_{i\pm 1/2}$

$f\langle u^{-},u^{+}\rangle =\frac{1}{2}\left[f\left(u^{-}\right)+f\left(u^{+}\right)-C\left(u^{+}-u^{-}\right)\right]$

其中 $C=\underset{u}{\max }\left|f^{\prime }\left(u\right)\right|$。整理方程组得

$A\ast \frac{\text{d}{u}_h}{\text{d}t}=g$， $u_h=\left(\begin{array}{c}{u}_j^{\left(1\right)}\\ u_j^{\left(2\right)}\\ u_j^{\left(3\right)}\end{array}\right)$， $g=\left(\begin{array}{c}\frac{f_{i-1/2}-f_{i-1/4}}{h}\\ \frac{f_{i-1/4}-f_{i+1/4}}{h}\\ \frac{f_{i+1/4}-f_{i+1/2}}{h}\end{array}\right)$

当基函数取线性元与二次元分别对应的矩阵A如下

$A_{\left(p^1\right)}=\left(\begin{array}{ccc}\frac{3}{16}& \frac{1}{16}& 0\\ \frac{1}{16}& \frac{6}{16}& \frac{1}{16}\\ 0& \frac{1}{16}& \frac{3}{16}\end{array}\right)$， $A_{\left(p^2\right)}=\left(\begin{array}{ccc}\frac{1}{6}& \frac{5}{48}& \frac{-1}{48}\\ \frac{1}{48}& \frac{11}{24}& \frac{1}{48}\\ \frac{-1}{48}& \frac{5}{48}& \frac{1}{6}\end{array}\right)$。

2.2. 斜率限制器

保极值的限制器

为了抑制数值振荡，引入MPS (Maximum-principle-satisfying)高阶精度限制器。

${\tilde{u}}_j\left(x\right)=\theta \left(u_j\left(x\right)-{\bar{u}}_j\right)+{\bar{u}}_j$， $\theta =\min \left\{\left|\frac{u_{\max }-{\bar{u}}_j}{M_j-{\bar{u}}_j}\right|,\left|\frac{u_{\min }-{\bar{u}}_j}{m_j-{\bar{u}}_j}\right|,1\right\}$

$M_j=\underset{x\in I_j}{\max }{u}_j\left(x\right)$， $m_j=\underset{x\in I_j}{\min }{u}_j(\; x\; )$

其中 ${\bar{u}}_j=\frac{1}{\Delta x}{\displaystyle {\int }_{I_j}{u}_j\left(x\right)\text{d}x}={\displaystyle \underset{\partial =1}{\overset{N}{\sum }}{\omega }_{\partial }{u}_j\left(x\right)}$，设 ${\omega }_{\partial }$ 是区间 $\left[-\frac{1}{2},\frac{1}{2}\right]$ 上的正交权重，可以令 ${\displaystyle \underset{\partial =1}{\overset{N}{\sum }}{\omega }_{\partial }}=1$。

在该限制器中 $u_j$ 表示在 $t_n$ 时刻单元上的数值解，其单元平均值 ${\bar{u}}_j$ 属于区间 $\left[u_{\min },u_{\max }\right]$，参考 [15] [16]。

Shu在 [17] 中证明了该限制器不会破坏高阶的精度。

2.3. 时间离散

半离散后得到常微分方程组，使用SSP Runge-Kutta格式 [18]

$\begin{array}{l}{u}^{\left(1\right)}=u^n+\Delta t\cdot L\left(u^n\right)\\ u^{\left(2\right)}=\frac{3}{4}{u}^n+\frac{1}{4}{u}^{\left(1\right)}+\frac{1}{4}\Delta t\cdot L\left(u^{\left(1\right)}\right)\\ u^{n+1}=\frac{1}{3}{u}^n+\frac{2}{3}{u}^{\left(2\right)}+\frac{2}{3}\Delta t\cdot L\left(u^{(\; 2\; )}\right)\end{array}$

3. 算例

3.1. 算例1线性对流方程

$\{\begin{array}{l}{u}_t+u_x=0\\ u\left(x,0\right)=\frac{1}{2}+\frac{1}{4}\sin (\; \pi \; x\; )\end{array}$

计算讨论在周期边界条件下的该算例在相同时间下，基函数分别为线性元、二次元时，它们的误差与收敛阶。计算结果见表1和表2。

通过例子可以看出，对于 $p^1$ 元，实际计算可以达到二阶精度，而对 $p^2$ 元，计算可以达到三阶精度。

3.2. 算例2 Burgers方程

$\{\begin{array}{l}{u}_t+{\left(\frac{u^2}{2}\right)}_x=0\\ u\left(x,0\right)=\frac{1}{2}+\frac{1}{4}\sin (\; \pi \; x\; )\end{array}$

边界条件选取周期边界条件，在 $t=\frac{1}{\pi }\approx 0.3183$ 时出现激波。由图1~2可知当 $t=0.3$ 时结果没有激波出现，由图3~6可知当 $t=\frac{2}{\pi }$ 和 $t=1.1$ 时，结果出现激波。可以看到数值解在间断附近有良好的分辨率。

3.3. 算例3 Buckely-Leverett方程

$\{\begin{array}{l}{u}_t+{\left(\frac{4u^2}{4u^2+{\left(1-u\right)}^2}\right)}_x=0\\ u_0=\{\begin{array}{l}1-0.5<x<0\\ 0其他\end{array}\end{array}$

由图7~8可知，通过比较数值解与五阶WENO格式下的参考解可以发现，保极值间断有限体积元方法下的数值解在大梯度处无震荡，而且与参考解有完全类似地良好计算效果。

3.4. 算例4 交通流问题

考虑一段20 km长的公路，入口的密度是50 veh/km。(入口密度用单位长度路段内的车辆总数来表示)对于所给出的初值条件，我们考虑入口条件变化。在两个小时段分析，最初10分钟入口关闭，10分钟到30分钟，提高到入口密度75 veh/km，30分钟后入口密度恢复到50 veh/km。初值条件见图9。

$\{\begin{array}{l}{u}_t+f{\left(u\right)}_x=0\\ u_0=\{\begin{array}{l}500\le x<10\hfill \\ 35010\le x<15\hfill \\ 1400-70x15\le x<20\hfill \end{array}\end{array}$

$f\left(u\right)=\{\begin{array}{l}-0.4u^2+100u0\le u<50\hfill \\ -0.1u^2+15u+350050\le u<100\hfill \\ -0.024u^2-5.2u+4760100\le u<350\hfill \end{array}$

计算结果见图10~13，在本算例中，分析了90分钟内因入口条件变化导致的交通流问题，通过将各个时间节点的车流密度变化图与该算例的准确解作对比，可以看出吻合情况较好。

4. 结论

本文研究了保极值的间断有限体积元方法，通过使用MPS的限制器使其保持高阶精度，守恒的同时又能满足极值原理。通过计算结果进行对比分析，证实了该保极值的间断有限体积元方法具有其理论设计的高精度和对间断高分辨率的特性。

基金

内蒙古大学科研发展基金(21100-5187133)；内蒙古自治区人才开发基金项目(12000-1300020240)。

参考文献