一种可求解混合型双曲守恒律方程的高阶数值格式

A High Order Numerical Scheme for Solving Conservation Laws of Mixed Type

摘要: 在计算流体力学等领域中，双曲守恒方程及其对流占优问题可以通过加权基本无振荡(WENO)方法进行高精度的数值求解。本文旨在通过建立一个修正高阶有限差分WENO格式来求解双曲–椭圆混合型方程。通过引入特殊的通量分裂方法将通量分解为两部分，在这两个分量中分别应用双曲WENO算子，并在守恒律方程数值通量中加入高阶修正项，获得了一种可求解混合型守恒律方程的高精度有限差分WENO格式。该离散格式主要用于求解viscosity-capillarity容许性条件下的双曲–椭圆型范德华方程。数值测试验证了该算法的高精度和有效性，结果表明，该格式不仅能在强间断区域保持无振荡，在解的光滑部分保持高阶数值精度，还可以有效描述复杂波结构。

Abstract: Hyperbolic conservation law equations and convection-dominant situations in computational fluid dynamics and other areas can be solved numerically with great precision using the weighted essentially non-oscillatory (WENO) techniques. In this paper, we attempt to address the hyperbolic-elliptic mixed equations by developing a corrected high-order finite difference WENO scheme. By introducing a special flux splitting method, the flux is decomposed into two parts. We then apply the hyperbolic WENO operator to these two components, and add the higher order correction term to the numerical flux of the conservation laws, and finally a high-precision finite-difference WENO scheme is obtained. The discretization scheme is mainly used to solve the hyperbolic-elliptic van der Waals equations under the viscosity-capillarity admissibility criterion. It can be shown by numerical examples that the scheme not only can preserve the necessary no oscillation in the discontinuous region, but also retain high order numerical accuracy in the smooth part of the solution, and can effectively describe the complex wave structure.

1. 引言

若守恒律方程$u_t+f{\left(u\right)}_x=0$ 的雅可比矩阵$\frac{\partial f\left(u\right)}{\partial u}$ 可以通过实特征值相似对角化，则方程是双曲型的。

在众多激波捕捉方法中，ENO/WENO格式可以在解的光滑区域保持高精度，在强间断附件保持本质无震荡，已经被广泛用于求解双曲型守恒律方程[1] [2]。若雅可比矩阵包含复特征值，则方程是椭圆型。但在实际应用中，更多方程是混合双曲–椭圆类型，如弹性波动方程[3]，与洛伦兹方程相关的热对流方程[4]以及固体中的传播相边界[5]。在本文中，我们主要求解一类具有非单调应力–应变关系的混合型守恒律方程，即低温下的范德华流体[6]：

$\{\begin{array}{l}{v}_t+p{\left(w\right)}_x=0,w_t-v_x=0\hfill \\ v\left(x,0\right)=v^0\left(x\right),w\left(x,0\right)=w^0\left(x\right)\hfill \end{array}$ (1)

其中w是比容，v是速度，p是自变量为w的非单调压力函数。该方程的雅可比矩阵为$\pm \sqrt{p^{\prime }\left(w\right)}$ ，在理想气体状态下，压力p是递减函数，构成双曲方程，但在气液共存状态下，p可能在w的某个区域变为正值，使该方程呈现椭圆特点。为了方便讨论，我们主要分黎曼初始条件问题：

$u\left(x,0\right)=\{\begin{array}{ll}{u}_L,\hfill & x<0\hfill \\ u_R,\hfill & x>0.\hfill \end{array}$ (2)

Shearer研究发现该初始条件下范德华方程的解可由激波、稀疏波、相边界(波速为0且连接不同相解的激波)等基本波构成。在双曲区域可能形成激波[7]，在跨越椭圆区域时会产生相变，即跨越椭圆区域的间断，这揭示了解在椭圆区域不稳定的物理现象。由Truskinovsky [8]和Slemrod [9]提出的viscosity-capillarity容许性条件可用于选择具有物理意义的容许性弱解。该条件同样适用于选择黎曼初始条件下范德华方程的容许性弱解。方程(2)的viscosity-capillarity容许性条件为

$\{\begin{array}{l}{v}_t+p{\left(w\right)}_x=\mu v_{xx}-\omega {\mu }^2{w}_{xxx},v\left(x,0\right)=v_0\left(x\right),\hfill \\ w_t-v_x=0,w\left(x,0\right)=w_0\left(x\right),\hfill \end{array}$ (3)

其中$\mu >0$ 为粘度，$\omega {\mu }^2$ 为毛细界面系数，通常情况下$\omega \in \left[0,\frac{1}{4}\right]$ 。我们一般认为在$\omega =\frac{1}{4}$ 且$\mu \to 0$ 时选择的弱解$u={\left(v,w\right)}^{\text{T}}$ 为符合物理意义的容许性弱解。目前研究者们已经设计出很多一阶或二阶的格式来求解混合型范德华方程，其中比较经典的是二阶Lax-Friedrichs离散[10]，但该格式求出的数值解在相边界附近存在伪尖峰和过耗散。Shi Jin确定了数值尖峰的实际来源，并设计出一种二阶方法来消除尖峰[11]。适用于双曲方程的激波捕获方法可能不再适用于椭圆区域，因为椭圆区域无法进行局部逐场分解。基于这种情况，Shu设计了一种特殊的通量分裂方法[12]，在椭圆区域直接重构原始变量。如果黎曼初始条件包含椭圆区域，即使在理论上椭圆区域中没有精确解，这些不连续也很可能会出现。本文使用的通量分裂方法，将带有高阶修正项的有限差分WENO算子分别作用在两个分裂通量中，求解范德华方程并观察其容许性弱解跨越椭圆区域的行为。

2. 修正有限差分WENO格式

2.1. 有限差分WENO格式

我们首先对WENO格式进行简单描述。为了方便计算，等分计算区间$\left[x_a,x_b\right]$ ：

$x_a=x_{\frac{1}{2}}<x_{\frac{3}{2}}<\cdots <x_{N-\frac{1}{2}}<x_{N+\frac{1}{2}}=x_b.$ (4)

其中$x_{i+\frac{1}{2}}$ 为单元边界值，$\Delta x=\frac{b-a}{N}$ 为单元长度，$x_i$ 为单元中心。给出函数值$S=\left\{f_{i-2},f_{i-1},f_i,f_{i+1},f_{i+2},f_{i+3}\right\}$ ，f在单元边界的左侧重构值$f_{i+\frac{1}{2}}^{-}$ 可由S的偏左模板$S_5^{-}=\left\{f_{i-2},f_{i-1},f_i,f_{i+1},f_{i+2}\right\}$ 得到，右侧重构值$f_{i+\frac{1}{2}}^{+}$ 可由偏右模板$S_5^{+}=\left\{f_{i-1},f_i,f_{i+1},f_{i+2},f_{i+3}\right\}$ 得到。以下仅给出$f_{i+\frac{1}{2}}^{-}$ 的重构过程，通过镜像对称，可得到$f_{i+\frac{1}{2}}^{+}$ 的重构值。

$S_5^{-}=\left\{f_{i-2},f_{i-1},f_i,f_{i+1},f_{i+2}\right\}$ 可分为三个小模板$s^k=\left\{f_{i-2+k},f_{i-1+k},f_{i+k}\right\}\left(k=0,1,2\right)$ ，在每个小模板上分别构造二次拉格朗日插值多项式，并计算在$x=x_{i+\frac{1}{2}}$ 的值，得到：

$f^{\left(0\right)}\left(x_{i+\frac{1}{2}}\right)=\frac{3}{8}{f}_{i-2}-\frac{5}{4}{f}_{i-1}+\frac{15}{8}{f}_i$

$f^{\left(1\right)}\left(x_{i+\frac{1}{2}}\right)=-\frac{1}{8}{f}_{i-1}+\frac{3}{4}{f}_i+\frac{3}{8}{f}_{i+1}$

$f^{\left(2\right)}\left(x_{i+\frac{1}{2}}\right)=\frac{3}{8}{f}_i+\frac{3}{4}{f}_{i+1}-\frac{1}{8}{f}_{i+2}$ (5)

$f_{i+\frac{1}{2}}^{-}$ 可表示为$f^k\left(x_{i+\frac{1}{2}}\right)$ 的线性组合

$f_{i+\frac{1}{2}}^{-}={\displaystyle \sum _{k=0}^2{\omega }_{\left(k\right)}}{f}^{\left(k\right)}\left(x_{i+\frac{1}{2}}\right)$ (6)

${\omega }_{\left(k\right)}$ 为非线性权。WENO格式通过调整模板上非线性权的比重，使间断处的模板比重尽量小，实现了间断附近高精度的数值稳定。非线性权为

${\omega }_k=\frac{{\alpha }_k}{{\displaystyle \sum _{s=0}^2{\alpha }_s}},k=0,1,2$ (7)

其中${\alpha }_k=\frac{d_k}{ϵ+{\beta }_k}$ ，$\left\{d_0=\frac{1}{16},d_1=\frac{5}{8},d_2=\frac{5}{16}\right\}$ 为线性权，参数$\epsilon$ 为正实数以防分母为0，在我们的算例中，$\epsilon ={10}^{-8}$ ，${\beta }_k$ 为光滑指示子，满足

${\beta }_r={\displaystyle \sum _{l=1}^{k-1}{\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}\Delta }}{x}^{2l-1}{\left(\frac{{\partial }^l{f}^{\left(r\right)}\left(x\right)}{{\partial }^{l}x}\right)}^2\text{d}x$ (8)

具体有：

${\beta }_0=\frac{13}{12}{\left(f_{i-2}-2f_{i-1}+f_i\right)}^2+\frac{1}{4}{\left(f_{i-2}-4f_{i-1}+3f_i\right)}^2,$

${\beta }_1=\frac{13}{12}{\left(f_{i-1}-2f_i+f_{i+1}\right)}^2+\frac{1}{4}{\left(f_{i-}-f_{i+1}\right)}^2,$

${\beta }_2=\frac{13}{12}{\left(f_i-2f_{i+1}+f_{i+2}\right)}^2+\frac{1}{4}{\left(3f_i-4f_{i+1}+f_{i+2}\right)}^2.$ (9)

2.2. 修正有限差分WENO格式

考虑到迎风性质和混合型守恒律方程的特殊性，即其雅可比矩阵$\frac{\partial \left(f\left(u\right)\right)}{\partial u}$ 存在复特征值，我们将数值通量分为两部分：

$f\left(u\right)=f^{+}\left(u\right)+f^{-}\left(u\right),$ (10)

引入[12]中的通量分裂方法：

$f^{\pm }\left(u\right)=\frac{1}{2}\left(f\left(u\right)\pm \overline{\alpha }u\right),\overline{\alpha }=\left(\begin{array}{cccc}{\alpha }_1& & & \\ & {\alpha }_2& & \\ & & \ddots & \\ & & & {\alpha }_m\end{array}\right)$.(11)

对范德华方程，有如下计算结果：

${\alpha }_2=\underset{w}{\max }\left(\max \left(0,\frac{\sqrt{M^2-4p^{\prime }\left(w\right)}-M}{2}\right)\right),$

${\alpha }_1={\alpha }_2+M,$

$M=a\underset{w}{\max }\sqrt{\max \left(0,p^{\prime }\left(w\right)\right)},a>2.$ (12)

该分裂思想大致是构造逆矩阵$g\left(u\right)=M{\left(v,0\right)}^{\text{T}}$ ，然后计算找到最小的M使得雅可比矩阵$\frac{\partial \left(f\left(u\right)\pm g\left(u\right)\right)}{\partial u}$ 有不同的实特征值，$a=2$ 为临界值，当$a>2$ ，雅可比矩阵可通过实特征值相似对角化。一旦做到这一点，就很容易使用Lax-Friedrich通量分裂思想，即用合适的$\alpha$ 对$\alpha u$ 进行加减，以完成分裂[12]。

守恒律方程的空间项$-f{\left(u\right)}_x=-f^{+}{\left(u\right)}_x-f^{-}{\left(u\right)}_x$ ，可由以下守恒半离散形式逼近：

$L{\left(u\right)}_i=L^{+}{\left(u\right)}_i+L^{-}{\left(u\right)}_i=-\frac{1}{\Delta x}\left({\widehat{f}}_{i+\frac{1}{2}}^{+}-{\widehat{f}}_{i-\frac{1}{2}}^{+}\right)-\frac{1}{\Delta x}\left({\widehat{f}}_{i+\frac{1}{2}}^{-}-{\widehat{f}}_{i-\frac{1}{2}}^{-}\right),$ (13)

${\widehat{f}}_{i+\frac{1}{2}}$ 为数值通量。在有限差分的框架内，五阶离散意味着

$\frac{1}{\Delta x}\left({\widehat{f}}_{i+\frac{1}{2}}^{+}-{\widehat{f}}_{i-\frac{1}{2}}^{+}\right)={\frac{\partial f^{+}\left(u\right)}{\partial x}|}_{x_i}+O\left(\Delta x^r\right),\frac{1}{\Delta x}\left({\widehat{f}}_{i+\frac{1}{2}}^{-}-{\widehat{f}}_{i-\frac{1}{2}}^{-}\right)={\frac{\partial f^{-}\left(u\right)}{\partial x}|}_{x_i}+O\left(\Delta x^r\right).$ (14)

数值通量${\widehat{f}}_{i\pm \frac{1}{2}}^{\pm }$ 可由一个低阶通量项和泰勒展开式构成：

${\widehat{f}}_{i+\frac{1}{2}}^{\pm }=f_{i+\frac{1}{2}}^{\pm }+{\displaystyle \sum _{k=1}^{m-1}{a}_{2k}}\Delta x^{2k}{\left(\frac{{\partial }^{2k}}{\partial x^{2k}}{f}^{\pm }\right)}_{i+\frac{1}{2}}+O\left(\Delta x^{2m+1}\right),$ (15)

则该格式有2m阶精度。在本文中，取$m=2$ ，则可获得一个五阶数值通量：

${\widehat{f}}_{i+\frac{1}{2}}^{\pm }=f_{i+\frac{1}{2}}^{\pm }-{\frac{1}{24}\Delta x^2{f}_{xx}^{\pm }|}_{i+\frac{1}{2}}+{\frac{7}{5760}\Delta x^4{f}_{xxxx}^{\pm }|}_{i+\frac{1}{2}}+O\left(\Delta x^6\right)=\mathcal{D}\left[f_{i+\frac{1}{2}}^{\pm }\right]+{\mathcal{T}}_{i+\frac{1}{2}}\left[f^{\pm }\right],$ (16)

其中，$\mathcal{D}[\cdot ]$ 为低阶通量项，可由上节有限差分WENO格式逼近。${\mathcal{T}}_{i+\frac{1}{2}}[\cdot ]$ 为保证整个格式为五阶精度的高阶中心差分算子，可由函数值的线性组合计算：

${\mathcal{T}}_{i+\frac{1}{2}}\left[f^{\pm }\right]={\frac{1}{24}\Delta x^2{f}_{xx}^{\pm }|}_{i+\frac{1}{2}}+{\frac{7}{5760}\Delta x^4{f}_{xxxx}^{\pm }|}_{i+\frac{1}{2}},$

${\Delta x^2{f}_{xx}^{\pm }|}_{i+\frac{1}{2}}=\frac{1}{48}\left(-5f_{i-2}^{\pm }+39f_{i-1}^{\pm }-34f_i^{\pm }-34f_{i+1}^{\pm }+39f_{i+2}^{\pm }-5f_{i+3}^{\pm }\right)+O\left(\Delta x^6\right),$

${\Delta x^4{f}_{xxxx}^{\pm }|}_{i+\frac{1}{2}}=\frac{1}{2}\left(f_{i-2}^{\pm }-3f_{i-1}^{\pm }+2f_i^{\pm }+2f_{i+1}^{\pm }-3f_{i+2}^{\pm }+f_{i+3}^{\pm }\right)+O\left(\Delta x^6\right).$ (17)

在我们的方法中，空间项的离散可由以下步骤组成。

(1) 由$S_5^{-}=\left\{f_{i-2},f_{i-1},f_i,f_{i+1},f_{i+2}\right\}$ 和上述WENO重构方法，计算得到所有单元边界值$f_{i+\frac{1}{2}}^{-}$ ；

(2) 由$S_5^{+}=\left\{f_{i-1},f_i,f_{i+1},f_{i+2},f_{i+3}\right\}$ 和WENO重构方法，计算得到所有单元边界值$f_{i+\frac{1}{2}}^{+}$ ；

(3) 由$S=\left\{f_{i-2}^{+},f_{i-1}^{+},f_i^{+},f_{i+1}^{+},f_{i+2}^{+},f_{i+3}^{+}\right\}$ 计算所有高阶中心差分算子${\mathcal{T}}_{i+\frac{1}{2}}\left[f^{+}\right]$ ；

(4) 由$S=\left\{f_{i-2}^{-},f_{i-1}^{-},f_i^{-},f_{i+1}^{-},f_{i+2}^{-},f_{i+3}^{-}\right\}$ 计算所有高阶差分中心算子${\mathcal{T}}_{i+\frac{1}{2}}\left[f^{-}\right]$ ；

(5) 数值通量为：

${\widehat{f}}_{i+\frac{1}{2}}^{+}=\mathcal{D}\left[f_{i+\frac{1}{2}}^{+}\right]+{\mathcal{T}}_{i+\frac{1}{2}}\left[f^{+}\right],{\widehat{f}}_{i+\frac{1}{2}}^{-}=\mathcal{D}\left[f_{i+\frac{1}{2}}^{-}\right]+{\mathcal{T}}_{i+\frac{1}{2}}\left[f^{-}\right];$ (18)

(6) 空间算子离散格式为：

$\frac{\partial u_i}{\partial t}=-\frac{1}{\Delta x}\left({\widehat{f}}_{i+\frac{1}{2}}^{+}-{\widehat{f}}_{i-\frac{1}{2}}^{+}\right)-\frac{1}{\Delta x}\left({\widehat{f}}_{i+\frac{1}{2}}^{-}-{\widehat{f}}_{i-\frac{1}{2}}^{-}\right)$ . (19)

将守恒律方程写成ODE形式$u_t=L\left(u\right)$ ，时间项由三阶TVD Runge Kutta方法离散[13]：

$u^{\left(1\right)}=u^{\left(n\right)}+\Delta tL\left(u^{\left(n\right)}\right),$

$u^{\left(2\right)}=\frac{3}{4}{u}^{\left(n\right)}+\frac{1}{4}{u}^{\left(1\right)}+\frac{1}{4}\Delta tL\left(u^{\left(1\right)}\right),$

$u^{\left(n+1\right)}=\frac{1}{3}{u}^{\left(n\right)}+\frac{2}{3}{u}^{\left(2\right)}+\frac{2}{3}\Delta tL\left(u^{\left(2\right)}\right).$ (20)

2.3. 稳定性分析

通过对双曲方程使用式(10)~(12)的通量分裂方法，我们恢复了经典的Lax-Friedrichs格式，该格式可以改写为一个中心格式加上一个等于$\frac{1}{2}\alpha \Delta x{u}_{xx}$ 的耗散项。如果使用式(11)的分裂，我们仍可以把它看成一个中心格式加上形式为$\frac{1}{2}\overline{\alpha }\Delta x{u}_{xx}$ 的耗散项。因此，参考[14]，可以合理的期望该格式收敛于可容许性弱解。我们可以进行大量的数值实验来评估高阶格式的收敛性和可容许性，第3节的数值算例仅是该方向初步的计算结果。

如果对分裂式(10)使用分步法，我们最终得到如下形式的格式(式(13))：

$u^{n+1}=\left(I+\Delta t{L}^{+}\right)\left(I+\Delta t{L}^{-}\right)u^n$ (21)

对于具有光滑解的线性或非线性问题，由于式(10)中$f^{\pm }\left(u\right)$ 的双曲性，选择稳定算子$\left(I+\Delta t{L}^{\pm }\right)$ 不太困难。但这并不代表格式(21)是稳定的，因为$\left(I+\Delta t{L}^{\pm }\right)$ 不能相互交换，也不能同时对角化。对任意一致的范数，如果算子满足更严格的条件如式(22)，则分步格式(21)是稳定的。

$\Vert I+\Delta t{L}^{+}\Vert \le 1+\text{Ο}\left(\Delta t\right),\Vert I+\Delta t{L}^{-}\Vert \le 1+\text{Ο}\left(\Delta t\right)$ (22)

3. 数值算例

本节先使用标量方程进行精度测试，随后计算范德华方程以验证格式的有效性。

3.1. 精度测试

例1. 一维无粘标量方程，使用周期边界条件，

$u_t+u_x=0;u\left(x,0\right)=0.5+\sin \left(\pi x\right),0\le x\le 2.$

精确解为$u\left(x,t\right)=0.5+\text{sin}\left(\pi \left(x-\frac{1}{2\pi }\right)\right)$ 。使用最后时刻的精确解和数值解计算误差范数：

$L^2=\frac{1}{N+1}\sqrt{{\displaystyle \sum _{i=0}^N{\left|u_i^n-u\left(x_i,t\right)\right|}^2}},L^{\infty }=\underset{0\le i\le N}{\max }\left|u_i^n-u\left(x_i,t\right)\right|,$ (23)

计算到$T=\frac{1}{2\pi }$ ，CFL = 0.2。表1、表2分别给出修正前后WENO格式的计算结果，可看出修正后的格式计算误差更小，且均可达到五阶精度。

例2. 有粘标量方程，使用周期边界条件：

$u_t+u_x=r{u}_{xx};u\left(x,0\right)=\sin \left(x\right),0\le x\le 2\pi$

精确解为$u\left(x,t\right)={\text{e}}^{-rt}\sin \left(x-t\right)$ 。扩散项$r{u}_{xx}$ 由简单的中心差分近似：

$r{u}_{xx}=\frac{-u_{i-2}+16u_{i-1}-30u_i+16u_{I+1}-u_{i+2}}{12{\left(\Delta x\right)}^2}$ (24)

扩散系数$r=0.01$ 。表3、表4给出$T=1$ 修正前后WENO格式的计算结果，精度有较明显的提高。

Table 1. Accuracy of corrected fifth WENO on $u_t+u_x=0,u\left(x,0\right)=0.5+\sin \left(\pi x\right)$

表1. $u_t+u_x=0,u\left(x,0\right)=0.5+\sin \left(\pi x\right)$ ，五阶修正WENO计算精度

N	$L^2$ error	$L^2$ order	$L^{\infty }$ error	$L^{\infty }$ order
10	1.11E−03	-	3.53E−03	-
20	2.59E−06	5.42	1.03E−04	5.09
40	5.39E−07	5.59	3.02E−06	5.10
80	1.13E−08	5.57	8.93E−08	5.08
160	2.43E−10	5.54	2.72E−09	5.04
320	5.30E−12	5.52	8.45E−11	5.01

Table 2. Accuracy of fifth WENO on $u_t+u_x=0,u\left(x,0\right)=0.5+\sin \left(\pi x\right)$

表2. $u_t+u_x=0,u\left(x,0\right)=0.5+\sin \left(\pi x\right)$ ，五阶WENO计算精度

N	$L^2$ error	$L^2$ order	$L^{\infty }$ error	$L^{\infty }$ order
10	1.40E−03	-	4.49E−03	-
20	3.33E−05	5.40	1.33E−04	5.07
40	6.92E−07	5.59	3.88E−06	5.10
80	1.46E−08	5.57	1.45E−07	5.08
160	3.13E−10	5.54	3.49E−09	5.03
320	6.81E−12	5.52	1.09E−10	5.01

Table 3. Accuracy of corrected fifth WENO on $u_t+u_x=r{u}_{xx},u\left(x,0\right)=\sin \left(x\right)$

表3. $u_t+u_x=r{u}_{xx},u\left(x,0\right)=\sin \left(x\right)$ ，五阶修正WENO计算精度

N	$L^2$ error	$L^2$ order	$L^{\infty }$ error	$L^{\infty }$ order
10	1.27E−03	-	5.43E−03	-
20	3.39E−05	5.23	1.86E−04	4.87
40	7.03E−07	5.59	5.53E−06	5.10
80	1.48E−08	5.57	1.65E−07	5.04
160	3.18E−10	5.54	5.11E−09	5.01
320	6.81E−12	5.54	1.55E−10	5.04

Table 4. Accuracy of fifth WENO on $u_t+u_x=r{u}_{xx},u\left(x,0\right)=\sin \left(x\right)$

表4. $u_t+u_x=r{u}_{xx},u\left(x,0\right)=\sin \left(x\right)$ ，五阶WENO计算精度

N	$L^2$ error	$L^2$ order	$L^{\infty }$ error	$L^{\infty }$ order
10	2.27E−03	-	7.51E−03	-
20	6.19E−05	5.20	2.53E−04	4.88
40	1.29E−06	5.58	7.20E−06	5.13
80	2.72E−08	5.57	2.18E−07	5.04
160	5.57E−10	5.53	6.72E−09	5.02
320	1.27E−11	5.53	2.05E−10	5.04

3.2. 范德华方程

使用第二节所获得的WENO5格式以及通量分裂方法，$a=2.2$ ，时间步长$\text{Δ}t=\frac{\text{CFLΔ}x}{\rho \left(\frac{\partial f^{+}\left(u\right)}{\partial u}+\frac{\partial f^{-}\left(u\right)}{\partial u}\right)}$ ，谱半径$\rho (\ast )$ 为

$\rho \left(\frac{\partial f^{+}\left(u\right)}{\partial u}+\frac{\partial f^{-}\left(u\right)}{\partial u}\right)=\underset{w}{\max }\left(\frac{1}{2}\left(2{\alpha }_2+M+\sqrt{M^2-4p^{\prime }\left(w\right)}\right)\right).$ (23)

例3 $p\left(w\right)=\frac{1}{w-0.25}-\frac{0.9}{w^2}$ ，图1是该压力状态函数的图像。$w<c$ 或$w>d$ 时为理想气体状态。

$p^{\prime }\left(w\right)<0$ ($c=0.574912$ , $d=1.036251$ )，构成双曲方程。$c<w<d$ 时，$p^{\prime }\left(w\right)$ 为正，呈现椭圆特点。图1中BE为Maxwell线(与曲线相交于点$m=0.494273$ 和$M=1.405065$ ，使得阴影部分面积相等)，画水平线AD，CF得到点$a=0.483100$ 和$f=1.918618$ ，为了观察解的状态，在以下图中均用两条水平虚线表示椭圆区域。CFL = 0.6，使用特征边界条件，我们先计算几个初始黎曼条件问题。

Figure 1. The image of $p\left(w\right)$

图1. $p\left(w\right)$ 的图像

(1) $\left(v_L,w_L\right)=\left(1,m\right)$ ，$\left(v_R,w_R\right)=\left(1,M\right)$ ，$T=25$ 。m和M为上述的Maxwell值。该初始条件满足热力学理论和viscosity-capillarity容许性条件[9]。图2给出了该条件下的数值弱解，我们发现数值解在双曲区域保持稳定，在跨越椭圆区域时伴随有大幅度的跳跃。此外，在发生相变的椭圆区域内，数值解中最多包含两个过渡点，但它们并不影响数值解的稳定性。

(2) $\left(v_L,w_L\right)=\left(1,0.54\right)$ ，$\left(v_R,w_R\right)=\left(1,1.8517\right)$ ，$T=4$ 。该初始条件位于两端的双曲区域，满足平稳跳

跃的Rankine-Hugoniot条件，但物理原理和viscosity-capillarity容许性条件[9]表明这种跳跃是不容许的。数值结果由图3(a)表示，我们观察到解逐渐演变为更复杂的跳跃结构，这显然是由2.3节所述的该方案固有的数值粘性所导致。除此之外，在相边界附近出现了振荡行为，这是因为在椭圆区域中，我们的格式是逐分量估计的，由于无法识别相应特征场，波分裂为两个或多个波(有些是双曲的)，产生了数值振荡。但是与[12]中ENO格式的计算结果相比，修正的WENO格式所产生的数值振荡明显有所控制且没有出现较为明显的数值尖峰，这表明我们的格式更能抑制伪振荡。当我们继续细化网格时，振荡减小且更受限制，解接近稳定状态。图3(b)给出了2000网格时的收敛解，实验表明，即使存在这些振荡，该方案仍具有收敛性。

图3(a)使用修正的WENO格式，由2000网格点计算。它与4000点的计算结果一致，因此可以认为是一个收敛解。为检验弱解是否满足viscosity-capillarity容许性条件，对流项使用修正WENO方法和通量分裂，扩散项用五阶差分格式近似，使用龙格库塔方法离散时间项。我们通过反复细化网格，对每个固定的$\mu$ 验证其分辨率，直至解呈稳定状态(使用最大的网格点为6000)，图4给出了$\omega =\frac{1}{4}$ 和$\mu \to 0$ ($\mu =0.1,0.01,0.001$ 时的计算结果，我们可以清楚看到$\mu \to 0$ 时的极限解收敛于我们的数值解。

Figure 2. $\left(v_L,w_L\right)=\left(1,m\right)$ , $\left(v_R,w_R\right)=\left(1,M\right)$ , $T=25$ , numerical solution (+, 200 points), convergent solution (solid line, 2000 points)

图2. $\left(v_L,w_L\right)=\left(1,m\right)$ ，$\left(v_R,w_R\right)=\left(1,M\right)$ ，$T=25$ ，数值解(+，200点)，收敛解(实线，2000点)

(a)

(b)

Figure 3. $\left(v_L,w_L\right)=\left(1,0.54\right)$ , $\left(v_R,w_R\right)=\left(1,1.8517\right)$ , $T=4$ , (a) Numerical solution (+, 200 points), convergent solution (solid line, 2000 points); (b) convergent solution (2000 points)

图3. $\left(v_L,w_L\right)=\left(1,0.54\right)$ ，$\left(v_R,w_R\right)=\left(1,1.8517\right)$ ，$T=4$ ，(a) 数值解(+，200点)，收敛解(实线，2000点)；(b) 收敛解(2000点)

Figure 4. Numerical solutions of viscosity-capillarity equations with $\mu =0.1,0.01,0.001$ , 2000 points

图4. Viscosity-capillarity方程数值解，$\mu =0.1,0.01,0.001$ ，2000点

(3) $\left(v_L,w_L\right)=\left(1,0.45\right)$ ，$\left(v_R,w_R\right)=\left(2,1.5\right)$ ，$T=1.5$ 。这种情况比前两个更容易计算，因为该初始条件是一个稳定的可容许弱解。图5(a)给出了数值解(+，200网格)和收敛解(实线，2000网格)，图5(b)表明，$\mu \to 0$ ($\mu =0.1,0.01,0.001$ )和$\omega =\frac{1}{4}$ 时viscosity-capillarity方程的解收敛于我们计算所得的数值解。

(a)

(b)

Figure 5. $\left(v_L,w_L\right)=\left(1,0.45\right)$ , $\left(v_R,w_R\right)=\left(2,1.5\right)$ , $T=1.5$ , (a) Numerical solution (+, 200 points), convergent solution (solid line, 2000 points); (b) Numerical solutions of viscosity-capillarity equations with $\mu =0.1,0.01,0.001$ , 2000 points

图5. $\left(v_L,w_L\right)=\left(1,0.45\right)$ ，$\left(v_R,w_R\right)=\left(2,1.5\right)$ ，$T=1.5$ ，(a) 数值解(+，200点)，收敛解(实线，2000点)，(b) Viscosity-capillarity方程数值解，$\mu =0.1,0.01,0.001$ ，2000点

例2. 压力状态$p\left(w\right)=w-w^3$ ，此时椭圆区域为$w\in \left(-\frac{\sqrt{3}}{3},\frac{\sqrt{3}}{3}\right)$ ，初始条件为$\left(v_L,w_L\right)=\left(-0.327,0.7\right)$ ，$\left(v_R,w_R\right)=\left(-0.327,-0.7\right)$ ，CFL = 0.6，计算到$T=2$ 。图6(a)给出了数值解和收敛解。如图所示，解包含向左向右移动的激波，激波之间存在稳定的相边界。图6(b)给出了viscosity-capillarity方程的收敛解。

(a)

(b)

Figure 6. $\left(v_L,w_L\right)=\left(-0.327,0.7\right)$ , $\left(v_R,w_R\right)=\left(-0.327,-0.7\right)$ , $T=2$ , (a) Numerical solution (+, 200 points), convergent solution (solid line, 2000 points); (b) Numerical solutions of viscosity-capillarity equations with $\mu =0.1,0.01,0.001$ , 2000 points

图6. $\left(v_L,w_L\right)=\left(-0.327,0.7\right)$ ，$\left(v_R,w_R\right)=\left(-0.327,-0.7\right)$ ，$T=2$ ，(a) 数值解(+，200点)，收敛解(实线，2000点)；(b) viscosity-capillarity方程数值解，$\mu =0.1,0.01,0.001$ ，2000点

例3. 在例1的压力状态下，使用光滑初始条件以及周期边界条件。

(1) $\left(v^0\left(x\right),w^0\left(x\right)\right)=\left(1-0.5\cos \left(x\right),1+0.5\sin \left(x\right)\right)$ ，该条件下初始条件包含椭圆区域。如图7(a)所示，$t<1$ 时，部分数值解包含在椭圆区域中，随着时间推移，解逐渐变成由相变所连接的分段光滑解，并完全包含在两个双曲区域内。

(2) $\left(v^0\left(x\right),w^0\left(x\right)\right)=\left(1-0.5\cos \left(x\right),0.8+0.2\sin \left(x\right)\right)$ ，初始条件完全包含在椭圆区域内。如图7(b)所示，此条件下的解结构与(1)类似。与ENO格式[12]相比，我们的方法在椭圆区域没有振荡行为且跨越椭圆区域时解的演变更为光滑。当$t<1$ 时，数值解在计算区域的两端振荡更剧烈，$t>1$ 后振荡行为消失。这大概率是因为随着格式精度的提高，$t<1$ 时计算区域两端会出现龙格库塔现象，$t>1$ 时数值解的演化变得相对稳定。

(a)

(b)

Figure 7. t = 0, 0.2, 0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, 2 (400 points), (a) $\left(v^0\left(x\right),w^0\left(x\right)\right)=\left(1-0.5\cos \left(x\right),1+0.5\sin \left(x\right)\right)$ , (b) $\left(v^0\left(x\right),w^0\left(x\right)\right)=\left(1-0.5\cos \left(x\right),0.8+0.2\sin \left(x\right)\right)$

图7. t = 0，0.2，0.4，0.6，0.8，1，1.2，1.4，1.6，1.8，2 (400点)，(a) $\left(v^0\left(x\right),w^0\left(x\right)\right)=\left(1-0.5\cos \left(x\right),1+0.5\sin \left(x\right)\right)$ ；(b) $\left(v^0\left(x\right),w^0\left(x\right)\right)=\left(1-0.5\cos \left(x\right),0.8+0.2\sin \left(x\right)\right)$

4. 结论

本文通过引入特殊的通量分裂方法，构造了一种修正的五阶有限差分WENO格式来求解混合双曲–椭圆型范德华方程。我们所关心的是数值解在椭圆区域内或跨越椭圆区域的行为，测试结果表明，我们的格式可以有效模拟包含相边界的复杂波结构，并在整个椭圆区域内振荡较小。与其他方法相比，该格式可以有效减少过渡点，或者过渡点在椭圆区域内有明显的后移，这表明数值解可以更快地穿过椭圆区域，即双曲–椭圆的过渡区域出现了明显收缩。此外，我们可以将该方法应用于更多的混合型方程，例如与Lorenz系统相关的PDE的存在性，这为我们未来的研究提供了建议。

基金项目

山西省自然科学基金项目资助(202103021224041)。

NOTES

^*通讯作者。