DOI: 10.12677/aam.2025.144145, PDF, HTML, XML,  被引量    科研立项经费支持
作者: 王亚博：太原理工大学数学学院，山西 晋中；刘红霞^*：太原理工大学人工智能学院，山西 晋中
关键词: 有限差分方法；对流扩散方程；HWENO格式；Hermite插值；Finite Difference Method； Convection-Diffusion Equations； HWENO Schemes； Hermite Interpolation

文章引用：王亚博, 刘红霞. 修正有限差分HWENO方法求解对流扩散方程[J]. 应用数学进展, 2025, 14(4): 122-134. https://doi.org/10.12677/aam.2025.144145

1. 引言

在本文中，我们根据求解双曲守恒律的数值方法[1]和文献[2]中的思想从而提出了一种修正的高阶Hermite加权本质非振荡(HWENO)方法来求解有限差分框架下对流占优的对流扩散方程。HWENO方案的构造基于从本质非振荡 (ENO)方案推导的加权本质非振荡(WENO)方案。它们是求解双曲守恒律的较有效和流行的数值方案。Harten和Osher [3] [4]开创性的基于总变分递减(TVD)方法构造了有限体积ENO方案来求解一维问题。Harten还在[5]中完成了二维问题的有限体积ENO方案的研究。在此之后，Casper和Atkins在[6] [7]中开发了一种多维、高阶精度的有限体积ENO方案。随后，Shu和Osher [2] [8]于1988年和1989年提出了一类用于解决多维问题的有限差分ENO方案。1994年，Liu等[9]在有限体积框架内提出了第一个WENO方案，这是第r阶ENO方案的改进版本，其数值精度可以达到(r + 1)阶。后来，Jiang和Shu [10]提出了一个有限差分WENO方案，使用相同的模板可以达到(2r − 1)级精度，并给出了方案的一般框架。与ENO方案相比，WENO方案的重构过程也采用了自适应模板的思想。区别在于，ENO方案仅使用所有候选模板中的一个，而WENO方案使用所有候选模板的非线性凸组合，并为每个模板分配与该模板上数值解的光滑度相关的非线性权重。WENO方案与ENO方案相比具有一些如[11]中所述的优势。

由于高阶WENO方案需要使用有关目标单元及其相邻单元中的解的信息，因此高阶WENO方案需要使用更宽的模板。为了同时实现高阶精度和紧凑性，Qiu和Shu [12] [13]提出了一类基于WENO方案的有限体积HWENO方案，并将其用作Runge-Kutta间断Galerkin方法的限制器。HWENO方案不仅具有像WENO方案一样随时间演变的数值推进格式，而且这些解及其一阶导数值都随时间演变并参与空间重构。因此，HWENO方案以增加额外的工作为代价实现紧凑性。然而，这些方案使用不同的模板来获得不同的多项式来重构解及其一阶导数，因此这些方案不够稳定和鲁棒，无法解决一些问题，例如这些方案不能很好地模拟双马赫和向前步长问题[14]。后来，Liu和Qiu [1] [15]提出了一维和二维双曲守恒律的高阶有限差分HWENO方案，并解决了这个问题。在此方案中，与一维方案的五阶精度相比，二维方案的精度降低为四阶。随后，Zhao等[16]通过使用Hermite限制器修改了解的导数，以提高方案的稳定性，并使有限差分HWENO方案在一维和二维情况下都达到五阶精度。对HWENO方案的更多参考可以在[17]-[20]及其参考文献中找到。

双曲守恒律与对流占优的对流扩散方程密切相关，因此我们希望应用HWENO方案在有限差分框架上求解对流扩散方程。在本文方法中，我们使用HWENO重构作为空间离散化，使用Runge-Kutta方法作为时间离散化，使用中心差分公式作为扩散项离散化。二维HWENO重构时我们使用逐维重构的思想，这既提高了计算效率，又使其易于推进到更高的维度。与[1]相比，对流项的解有几个共同点：首先，可以使用任意单调通量；第二种是直接使用解及其一阶导数点值插值进行计算，而不是使用函数及其一阶导数值进行计算。也有一些差异：首先，我们使用中心差分公式而不是通量分裂来计算高阶通量[2]；其次，我们使用Hermite插值而不是WENO重构来计算二维混合导数。改进的HWENO方法在一维情况下可以达到五阶精度，但由于混合导数的影响，在二维情况下只能达到四阶精度。综上所述，修改后的 HWNEO方法大大简化了计算并提高了效率，同时也保持了之前方法的优点。

2. 一维数值格式的构造

为简单起见，我们首先考虑一维标量对流扩散方程：

$\{\begin{array}{l}\frac{\partial u}{\partial t}+\frac{\partial f\left(u\right)}{\partial x}=\epsilon \frac{{\partial }^{2}u}{\partial x^2},\\ u\left(x,0\right)=u_0\left(x\right),\end{array}$ (1)

其中t是时间变量，x是一维情况下的空间变量。 $\epsilon$ 则是扩散系数。此外，u是守恒量，f(u)是空间x方向上的通量函数。所提出的数值方法均定义在均匀网格中。计算区域[a, b]被N等分，也就是说

$a=x_{\frac{1}{2}}<x_{\frac{3}{2}}<\cdots <x_{N+\frac{1}{2}}=b$ ，其中 $x_{i+\frac{1}{2}}$ 和 $x_i$ 分别表示单元 $I_i=\left[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}\right]$ 的单元边界和单元中心。且有 $x_i=x_{i-\frac{1}{2}}+\frac{\Delta x}{2}$ ，其中 $\Delta x=\frac{b-a}{N}$ 是空间步长。函数 $u\left(x,t\right)$ 关于空间变量的一阶导数记为 $v\left(x,t\right)$ 。那么，从方程(1)和它关于空间变量的导数方程有

$\{\begin{array}{l}\frac{\partial u}{\partial t}+\frac{\partial f\left(u\right)}{\partial x}=\epsilon \frac{{\partial }^{2}u}{\partial x^2},u\left(x,0\right)=u_0\left(x\right),\hfill \\ \frac{\partial v}{\partial t}+\frac{\partial h\left(u,v\right)}{\partial x}=\epsilon \frac{{\partial }^{2}v}{\partial x^2},v\left(x,0\right)=v_0\left(x\right),\hfill \end{array}$ (2)

其中 $v_0\left(x\right)={\left(u_0\right)}_x$ ， $h\left(u,v\right)=f^{\prime }\left(u\right)v$ 。我们用以下半离散有限差分方案近似(2)

$\{\begin{array}{l}\frac{\partial u_i}{\partial t}=-\frac{1}{\Delta x}\left({\hat{f}}_{i+\frac{1}{2}}-{\hat{f}}_{i-\frac{1}{2}}\right)+\epsilon {\frac{{\partial }^{2}u}{\partial x^2}|}_{x_i},\hfill \\ \frac{\partial v_i}{\partial t}=-\frac{1}{\Delta x}\left({\hat{h}}_{i+\frac{1}{2}}-{\hat{h}}_{i-\frac{1}{2}}\right)+\epsilon {\frac{{\partial }^{2}v}{\partial x^2}|}_{x_i},\hfill \end{array}$ (3)

其中 ${\hat{f}}_{i+\frac{1}{2}}$ 和 ${\hat{h}}_{i+\frac{1}{2}}$ 分别是原始方程及其空间导数方程的数值通量。它们满足数值通量的一般条件，例如单调

性、Lipschitz连续性以及与物理通量 $f\left(u\right)$ 和 $h\left(u,v\right)$ 的一致性。

如果直接使用HWENO重构来计算数值通量，则只能获得低阶精度。为了提高精度，经典的有限差分HWNEO方法使用通量分裂来重构函数及其一阶导数。本章我们没有使用这种方法，而是使用[1]中求解双曲守恒律的想法来求解对流扩散方程。在一维有限差分框架中，五阶精度的离散化为

$\{\begin{array}{l}\frac{1}{\Delta x}\left({\hat{f}}_{i+\frac{1}{2}}-{\hat{f}}_{i-\frac{1}{2}}\right)={\frac{\partial f\left(u\right)}{\partial x}|}_{x_i}+O\left(\Delta x^6\right),\hfill \\ \frac{1}{\Delta x}\left({\hat{h}}_{i+\frac{1}{2}}-{\hat{h}}_{i-\frac{1}{2}}\right)={\frac{\partial h\left(u,v\right)}{\partial x}|}_{x_i}+O\left(\Delta x^5\right),\hfill \end{array}$ (4)

在[2]中，数值通量 ${\hat{f}}_{i+\frac{1}{2}}$ 和 ${\hat{h}}_{i+\frac{1}{2}}$ 的构造如下

$\{\begin{array}{l}{\hat{f}}_{i+\frac{1}{2}}={\hat{f}}_{i+\frac{1}{2}}^L+{\hat{f}}_{i+\frac{1}{2}}^H,\hfill \\ {\hat{h}}_{i+\frac{1}{2}}={\hat{h}}_{i+\frac{1}{2}}^L+{\hat{h}}_{i+\frac{1}{2}}^H,\hfill \end{array}$ (5)

其中由Taylor展开得到的高阶修正项 ${\hat{f}}_{i+\frac{1}{2}}^H$ 和 ${\hat{h}}_{i+\frac{1}{2}}^H$ 可以来提高数值通量的精度。低阶的通量项 ${\hat{f}}_{i+\frac{1}{2}}^L={\hat{f}}^L\left(u_{i+\frac{1}{2}}^{-},u_{i+\frac{1}{2}}^{+}\right)$ 和 ${\hat{h}}_{i+\frac{1}{2}}^L={\hat{h}}^L\left(u_{i+\frac{1}{2}}^{-},u_{i+\frac{1}{2}}^{+},v_{i+\frac{1}{2}}^{-},v_{i+\frac{1}{2}}^{+}\right)$ 是任意的单调通量。 $u_{i+\frac{1}{2}}^{\pm }$ 和 $v_{i+\frac{1}{2}}^{\pm }$ 是 $u\left(x,t\right)$ 和 $v\left(x,t\right)$ 在半点 $x=x_{i+\frac{1}{2}}$ 处由偏左或偏右模板的HWENO重构得到的数值近似。在这里，我们使用以下Lax-Friedrichs通量近似 ${\hat{f}}_{i+\frac{1}{2}}^L$ 和 ${\hat{h}}_{i+\frac{1}{2}}^L$

$\{\begin{array}{l}{\hat{f}}_{i+\frac{1}{2}}^L=\frac{1}{2}\left[f\left(u_{i+\frac{1}{2}}^{-}\right)+f\left(u_{i+\frac{1}{2}}^{+}\right)-\alpha \left(u_{i+\frac{1}{2}}^{+}-u_{i+\frac{1}{2}}^{-}\right)\right],\hfill \\ {\hat{h}}_{i+\frac{1}{2}}^L=\frac{1}{2}\left[h\left(u_{i+\frac{1}{2}}^{-},v_{i+\frac{1}{2}}^{-}\right)+h\left(u_{i+\frac{1}{2}}^{+},v_{i+\frac{1}{2}}^{+}\right)-\alpha \left(v_{i+\frac{1}{2}}^{+}-v_{i+\frac{1}{2}}^{-}\right)\right],\hfill \end{array}$ (6)

其中 $\alpha$ 定义为 ${\max }_u{f}^{\prime }\left(u\right)$ 。高阶修正项 ${\hat{f}}_{i+\frac{1}{2}}^H$ 和 ${\hat{h}}_{i+\frac{1}{2}}^H$ 的系数至少是 $\Delta x^2$ ，再加上一维情况下要达到五阶精度，

因此整体截断误差应达到 $O\left(\Delta x^6\right)$ 。即有

$\{\begin{array}{l}{\hat{f}}_{i+\frac{1}{2}}^H=-\frac{1}{24}\Delta x^2{\frac{{\partial }^{2}f}{\partial x^2}|}_{x_{i+\frac{1}{2}}}+\frac{7}{5760}\Delta x^4{\frac{{\partial }^{4}f}{\partial x^4}|}_{x_{i+\frac{1}{2}}}+O\left(\Delta x^6\right),\hfill \\ {\hat{h}}_{i+\frac{1}{2}}^H=-\frac{1}{24}\Delta x^2{\frac{{\partial }^{2}h}{\partial x^2}|}_{x_{i+\frac{1}{2}}}+\frac{7}{5760}\Delta x^4{\frac{{\partial }^{4}h}{\partial x^4}|}_{x_{i+\frac{1}{2}}}+O\left(\Delta x^6\right),\hfill \end{array}$ (7)

2.1. ${\hat{f}}_{i+\frac{1}{2}}^L$ 和 ${\hat{h}}_{i+\frac{1}{2}}^L$ 的构造

要构造 ${\hat{f}}_{i+\frac{1}{2}}^L$ 和 ${\hat{h}}_{i+\frac{1}{2}}^L$ 需要先重构 $u_{i+\frac{1}{2}}^{\pm }$ 和 $v_{i+\frac{1}{2}}^{\pm }$ 。现在我们详细描述由给定点值 $\left\{u_i,v_i\right\}$ 得到点值 $\left\{u_{i+\frac{1}{2}}^{\pm },v_{i+\frac{1}{2}}^{\pm }\right\}$ 的五阶HWENO插值过程。u和v在半点 $x=x_{i+\frac{1}{2}}$ 处由偏左模板 $S^{-}=\left\{x_{i-1},x_i,x_{i+1}\right\}$ 从左侧重构得到的值即为 $u_{i+\frac{1}{2}}^{-}$ 和 $v_{i+\frac{1}{2}}^{-}$ 。同理，u和v在半点 $x=x_{i+\frac{1}{2}}$ 处由偏右模板 $S^{+}=\left\{x_i,x_{i+1},x_{i+2}\right\}$ 从右侧重构得到的值即为 $u_{i+\frac{1}{2}}^{+}$ 和 $v_{i+\frac{1}{2}}^{+}$ 。以下是 $u_{i+\frac{1}{2}}^{-}$ 和 $v_{i+\frac{1}{2}}^{-}$ 的构造算法，而 $u_{i+\frac{1}{2}}^{+}$ 和 $v_{i+\frac{1}{2}}^{+}$ 的值可以用镜像对称的方法得到。

完整的五阶HWENO重构过程可以表示为

$U_{i+\frac{1}{2}}^{-}\approx {\displaystyle \underset{j=0}{\overset{2}{\sum }}{\omega }_j}{p}_j\left(x_{i+\frac{1}{2}}\right),$ (8)

其中 $U={\left(u,v\right)}^{\text{T}}$ 。 $p_j\left(x_{i+\frac{1}{2}}\right)$ 是多项式在半点 $x=x_{i+\frac{1}{2}}$ 处的函数值，该多项式根据重构对象满足在小模板上的对应插值性质。这些插值多项式满足的性质取决于重构对象，具体所满足的性质在[1]中已经给出。此外， ${\omega }_j$ 是非线性权重，它满足 ${\omega }_j\ge 0$ 和 ${\displaystyle \underset{j=0}{\overset{2}{\sum }}{\omega }_j}=1$ 。

非线性权 ${\omega }_j$ 的具体表达式如下

${\omega }_j=\frac{{\bar{\omega }}_j}{{\displaystyle \underset{k}{\sum }{\bar{\omega }}_k}},{\bar{\omega }}_k=\frac{{\gamma }_k}{{\left({\beta }_k+\epsilon \right)}^2},j,k=0,1,2,$ (9)

其中 ${\gamma }_j$ 是线性权，它满足 $q\left(x_{i+\frac{1}{2}}\right)={\displaystyle \underset{j=0}{\overset{2}{\sum }}{\gamma }_j}{p}_j\left(x_{i+\frac{1}{2}}\right)$ 和 ${\displaystyle \underset{j=0}{\overset{2}{\sum }}{\gamma }_j}=1$ 。 $q\left(x_{i+\frac{1}{2}}\right)$ 则类似于满足某些性质的插值多项式 $p_j\left(x_{i+\frac{1}{2}}\right)$ ，不同之处在于它的模板是在大模板 $S^{-}$ 上而不是小模板上，所满足的具体性质也可见[1]。

而 $\epsilon$ 可取任意能避免分母为零的正实数，在本文均取为10⁻⁶。 ${\beta }_j$ 是光滑指示子，我们使用[10]中如下的定义

${\beta }_j={\displaystyle \underset{l=k}{\overset{r}{\sum }}{\displaystyle {\int }_{I_i}^{}\Delta }}{x}^{2l-1}{\left(\frac{{\partial }^l}{\partial x^l}{p}_j\left(x\right)\right)}^2\text{d}x.$ (10)

求和上限r是多项式 $p_j\left(x\right)$ 的次数，求和下限k取决于重构对象，即重构 $u_{i+\frac{1}{2}}^{-}$ 时取k = 1，重构 $v_{i+\frac{1}{2}}^{-}$ 时取k = 2。

2.1.1. $u_{i+\frac{1}{2}}^{-}$ 的构造算法

插值多项式的函数值如下

$\begin{array}{l}{p}_0\left(x_{i+\frac{1}{2}}\right)=\frac{1}{4}\left(-5u_{i-1}+9u_i-3\Delta x{v}_{i-1}\right)\\ p_1\left(x_{i+\frac{1}{2}}\right)=\frac{1}{4}\left(u_i+3u_{i+1}-\Delta x{v}_{i+1}\right)\\ p_2\left(x_{i+\frac{1}{2}}\right)=\frac{1}{8}\left(-u_{i-1}+6u_i+3u_{i+1}\right)\\ q\left(x_{i+\frac{1}{2}}\right)=\frac{1}{64}\left(-8u_{i-1}+36u_i+36u_{i+1}-3\Delta x{v}_{i-1}-9\Delta x{v}_{i+1}\right)\end{array}$ (11)

而线性权的值为 ${\gamma }_0=\frac{1}{16},{\gamma }_1=\frac{9}{16},{\gamma }_2=\frac{3}{8}$ 。

光滑指示子 ${\beta }_j$ 的表达式如下

$\begin{array}{l}{\beta }_0={\left(-2u_{i-1}+2u_i-\Delta x{v}_{i-1}\right)}^2+\frac{13}{3}{\left(-u_{i-1}+u_i-\Delta x{v}_{i-1}\right)}^2,\\ {\beta }_1={\left(-2u_i+2u_{i+1}-\Delta x{v}_{i+1}\right)}^2+\frac{13}{3}{\left(-u_i+u_{i+1}-\Delta x{v}_{i+1}\right)}^2,\\ {\beta }_2=\frac{1}{4}{\left(-u_{i-1}+u_{i+1}\right)}^2+\frac{13}{12}{\left(-u_{i-1}+2u_i-u_{i+1}\right)}^2.\end{array}$ (12)

2.1.2. $v_{i+\frac{1}{2}}^{-}$ 的构造算法

插值多项式的函数值如下

$\begin{array}{l}{p}_0\left(x_{i+\frac{1}{2}}\right)=\frac{1}{4\Delta x}\left(18u_{i-1}-18u_i+7\Delta x{v}_{i-1}+15\Delta x{v}_i\right),\\ p_1\left(x_{i+\frac{1}{2}}\right)=\frac{1}{4\Delta x}\left(-6u_i+6u_{i+1}-\Delta x{v}_i-\Delta x{v}_{i+1}\right),\\ p_2\left(x_{i+\frac{1}{2}}\right)=\frac{1}{8\Delta x}\left(u_{i-1}-8u_i+7u_{i+1}+2\Delta x{v}_i\right),\\ q\left(x_{i+\frac{1}{2}}\right)=\frac{1}{64\Delta x}\left(3u_{i-1}-96u_i+93u_{i+1}+\Delta x{v}_{i-1}-12\Delta x{v}_i-15\Delta x{v}_{i+1}\right).\end{array}$ (13)

则线性权的值为 ${\gamma }_0=\frac{1}{112},{\gamma }_1=\frac{15}{16},{\gamma }_2=\frac{3}{56}$ 。

光滑指示子如下

$\begin{array}{l}{\beta }_0=\frac{13}{12}{\left(12u_{i-1}-12u_i+6\Delta x{v}_{i-1}+6\Delta x{v}_i\right)}^2+{\left(6u_{i-1}-6u_i+2\Delta x{v}_{i-1}+4\Delta x{v}_i\right)}^2,\\ {\beta }_1=\frac{13}{12}{\left(12u_i-12u_{i+1}+6\Delta x{v}_i+6\Delta x{v}_{i+1}\right)}^2+{\left(6u_i-6u_{i+1}+4\Delta x{v}_i+2\Delta x{v}_{i+1}\right)}^2,\\ {\beta }_2=\frac{13}{12}{\left(-3u_{i-1}+3u_{i+1}-6\Delta x{v}_i\right)}^2+{\left(u_{i-1}-2u_i+u_{i+1}\right)}^2.\end{array}$ (14)

那么这就完成了 $u_{i+\frac{1}{2}}^{\pm }$ 和 $v_{i+\frac{1}{2}}^{\pm }$ 的HWENO重构，然后我们将其代入到(8)就完成了 ${\hat{f}}_{i+\frac{1}{2}}^L$ 和 ${\hat{h}}_{i+\frac{1}{2}}^L$ 的计算。

2.2. ${\hat{f}}_{i+\frac{1}{2}}^H$ 和 ${\hat{h}}_{i+\frac{1}{2}}^H$ 的构造

从(7)可以看出，我们只需要分别构造通量函数 $f\left(u\right)$ 和 $h\left(u,v\right)$ 的二阶和四阶导数就完成了高阶修正

项 ${\hat{f}}_{i+\frac{1}{2}}^H$ 和 ${\hat{h}}_{i+\frac{1}{2}}^H$ 的计算。为了保持数值方法的紧凑性，我们在此使用以下Hermite方法推导的中心差分法对

高阶导数进行近似处理

$\begin{array}{c}{\frac{{\partial }^{2}f}{\partial x^2}|}_{x_{i+\frac{1}{2}}}=-\frac{27}{64\Delta x^2}{f}_{i-1}+\frac{27}{64\Delta x^2}{f}_i+\frac{27}{64\Delta x^2}{f}_{i+1}-\frac{27}{64\Delta x^2}{f}_{i+2}\\ -\frac{19}{192\Delta x}{h}_{i-1}-\frac{99}{64\Delta x}{h}_i+\frac{99}{64\Delta x}{h}_{i+1}+\frac{19}{192\Delta x}{h}_{i+2}\\ {\frac{{\partial }^{4}f}{\partial x^4}|}_{x_{i+\frac{1}{2}}}=\frac{45}{4\Delta x^4}{f}_{i-1}-\frac{45}{4\Delta x^4}{f}_i-\frac{45}{4\Delta x^4}{f}_{i+1}+\frac{45}{4\Delta x^4}{f}_{i+2}\\ +\frac{11}{4\Delta x^3}{h}_{i-1}+\frac{57}{4\Delta x^3}{h}_i-\frac{57}{4\Delta x^3}{h}_{i+1}-\frac{11}{4\Delta x^3}{h}_{i+2},\end{array}$

以及公式

$\begin{array}{l}{\frac{{\partial }^{2}h}{\partial x^2}|}_{x_{i+\frac{1}{2}}}=\frac{281}{288\Delta x^3}{f}_{i-1}+\frac{513}{32\Delta x^3}{f}_i-\frac{513}{32\Delta x^3}{f}_{i+1}-\frac{281}{288\Delta x^3}{f}_{i+2}\\ +\frac{19}{96\Delta x^2}{h}_{i-1}+\frac{297}{32\Delta x^2}{h}_i+\frac{297}{32\Delta x^2}{h}_{i+1}+\frac{19}{96\Delta x^2}{h}_{i+2},\\ {\frac{{\partial }^{4}h}{\partial x^4}|}_{x_{i+\frac{1}{2}}}=-\frac{785}{18\Delta x^5}{f}_{i-1}-\frac{345}{2\Delta x^5}{f}_i+\frac{345}{2\Delta x^5}{f}_{i+1}+\frac{785}{18\Delta x^5}{f}_{i+2}\\ -\frac{55}{6\Delta x^4}{h}_{i-1}-\frac{285}{2\Delta x^4}{h}_i-\frac{285}{2\Delta x^4}{h}_{i+1}-\frac{55}{6\Delta x^4}{h}_{i+2},\end{array}$

其中 $f_i$ 和 $h_i$ 分别是 $f\left(u_i\right)$ 和 $h\left(u_i,v_i\right)$ 的简写。

2.3. 扩散项的离散

由于对流占优的对流扩散方程与双曲守恒律密切相关，因此以上就完成了半离散格式(3)中对流项的离散化。在这里，我们仍然使用Hermite插值的思想来离散扩散项。扩散项中的高阶导数的离散表达式如下

$\{\begin{array}{l}{\frac{{\partial }^{2}u}{\partial x^2}|}_{x_i}=\frac{1}{2\Delta x^2}\left(4u_{i-1}-8u_i+4u_{i+1}+\Delta x{v}_{i-1}-\Delta x{v}_{i+1}\right),\\ {\frac{{\partial }^{2}v}{\partial x^2}|}_{x_i}=\frac{1}{2\Delta x^3}\left(-15u_{i-1}+15u_{i+1}-3\Delta x{v}_{i-1}-24\Delta x{v}_i-3\Delta x{v}_{i+1}\right),\end{array}$ (15)

其中该格式可以达到五阶精度。

2.4. 时间项的离散

完成上述步骤后，我们只需要离散时间项就能完成数值求解过程。在此，我们使用经典的三阶TVD Runge-Kutta方法来离散时间项

$\{\begin{array}{l}{U}^{\left(1\right)}=U^{\left(n\right)}+\Delta tL\left(U^{\left(n\right)}\right),\hfill \\ 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),\hfill \\ 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),\hfill \end{array}$ (16)

这样，我们就完成了一个完整的对流扩散方程的一维数值格式。

3. 二维数值格式的构造

然后我们考虑二维标量对流扩散方程

$\{\begin{array}{l}\frac{\partial u}{\partial t}+\frac{\partial f\left(u\right)}{\partial x}+\frac{\partial g\left(u\right)}{\partial y}=\epsilon \left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right),\hfill \\ u\left(x,y,0\right)=u_0\left(x,y\right),\hfill \end{array}$ (17)

其中 $u=u\left(x,y,t\right)$ 是守恒量， $f\left(u\right)$ 和 $g\left(u\right)$ 分别是空间变量x和y方向上的通量函数。为简化表述，数值

方法定义在均匀网格中。也就是说计算域被均匀网格点 $\left\{\left(x_{i+\frac{1}{2}},y_{j+\frac{1}{2}}\right)\right\}$ 划分，且 $i=0,\cdots ,N_x$ ， $j=0,\cdots ,N_y$ 。 $\Delta x=x_{i+\frac{1}{2}}-x_{i-\frac{1}{2}}$ 和 $\Delta y=y_{j+\frac{1}{2}}-y_{j-\frac{1}{2}}$ 分别是x和y方向上的空间步长。那么，单元和单元中心分别定义为 $I_{i,j}=\left[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}\right]\times \left[y_{j-\frac{1}{2}},y_{j+\frac{1}{2}}\right]$ 和 $\left(x_i,y_j\right)=\left(\frac{1}{2}\left(x_{i-\frac{1}{2}}+x_{i+\frac{1}{2}}\right),\frac{1}{2}\left(y_{j-\frac{1}{2}}+y_{j+\frac{1}{2}}\right)\right)$ 。我们将守恒量u关于x和y的一阶导数记为v和w。然后，从方程(17)及其导数方程得到

$\{\begin{array}{ll}\frac{\partial u}{\partial t}+\frac{\partial f\left(u\right)}{\partial x}+\frac{\partial g\left(u\right)}{\partial y}=\epsilon \left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right),\hfill & u\left(x,y,0\right)=u_0\left(x,y\right),\hfill \\ \frac{\partial v}{\partial t}+\frac{\partial h\left(u,v\right)}{\partial x}+\frac{\partial r\left(u,v\right)}{\partial y}=\epsilon \left(\frac{{\partial }^{2}v}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y^2}\right),\hfill & v\left(x,y,0\right)=v_0\left(x,y\right),\hfill \\ \frac{\partial w}{\partial t}+\frac{\partial q\left(u,w\right)}{\partial x}+\frac{\partial s\left(u,w\right)}{\partial y}=\epsilon \left(\frac{{\partial }^{2}w}{\partial x^2}+\frac{{\partial }^{2}w}{\partial y^2}\right),\hfill & w\left(x,y,0\right)=w_0\left(x,y\right),\hfill \end{array}$ (18)

其中

$\{\begin{array}{l}h\left(u,v\right)=f^{\prime }\left(u\right)v,r\left(u,v\right)=g^{\prime }\left(u\right)v,v_0\left(x,y\right)={\left(u_0\right)}_x,\hfill \\ q\left(u,w\right)=f^{\prime }\left(u\right)w,s\left(u,w\right)=g^{\prime }\left(u\right)w,w_0\left(x,y\right)={\left(u_0\right)}_y.\hfill \end{array}$

方程(18)的半离散有限差分格式为

$\begin{array}{l}\frac{\partial u_{i,j}}{\partial t}=-\frac{1}{\Delta x}\left({\hat{f}}_{i+\frac{1}{2},j}-{\hat{f}}_{i-\frac{1}{2},j}\right)-\frac{1}{\Delta y}\left({\hat{g}}_{i,j+\frac{1}{2}}-{\hat{g}}_{i,j-\frac{1}{2}}\right)+\epsilon {\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right)|}_{\left(x_i,y_j\right)},\\ \frac{\partial v_{i,j}}{\partial t}=-\frac{1}{\Delta x}\left({\hat{h}}_{i+\frac{1}{2},j}-{\hat{h}}_{i-\frac{1}{2},j}\right)-\frac{1}{\Delta y}\left({\hat{r}}_{i,j+\frac{1}{2}}-{\hat{r}}_{i,j-\frac{1}{2}}\right)+\epsilon {\left(\frac{{\partial }^{2}v}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y^2}\right)|}_{\left(x_i,y_j\right)},\\ \frac{\partial w_{i,j}}{\partial t}=-\frac{1}{\Delta x}\left({\hat{q}}_{i+\frac{1}{2},j}-{\hat{q}}_{i-\frac{1}{2},j}\right)-\frac{1}{\Delta y}\left({\hat{s}}_{i,j+\frac{1}{2}}-{\hat{s}}_{i,j-\frac{1}{2}}\right)+\epsilon {\left(\frac{{\partial }^{2}w}{\partial x^2}+\frac{{\partial }^{2}w}{\partial y^2}\right)|}_{\left(x_i,y_j\right)},\end{array}$ (19)

其中 $u\left(x_i,y_j,t\right),v\left(x_i,y_j,t\right)$ 和 $w\left(x_i,y_j,t\right)$ 的值近似表示为 $u_{i,j},v_{i,j}$ 和 $w_{i,j}$ 。需要注意的是，在二维对流扩散方

程中的对流项和扩散项中都存在混合导数。这里，对流项中的非混合数值通量 ${\hat{f}}_{i\pm \frac{1}{2},j},{\hat{h}}_{i\pm \frac{1}{2},j},{\hat{g}}_{i,j\pm \frac{1}{2}}$ 和 ${\hat{s}}_{i,j\pm \frac{1}{2}}$ 以及扩散项中的非混合导数 $\frac{{\partial }^{2}u}{\partial x^2},\frac{{\partial }^{2}v}{\partial x^2},\frac{{\partial }^{2}u}{\partial y^2}$ 和 $\frac{{\partial }^{2}w}{\partial y^2}$ 都是以类似的一维方式逐维重构的。唯一的例外是对流项中混合数值通量 ${\hat{q}}_{i+\frac{1}{2},j},{\hat{r}}_{i,j+\frac{1}{2}}$ 和扩散项中混合导数 $\frac{{\partial }^{2}v}{\partial y^2},\frac{{\partial }^{2}w}{\partial x^2}$ 的构造，它们由简单的三点Hermite插值得出的

四阶中心差分公式近似，而不是通常的通量分裂方法。再加上虽然对流项的处理可以使用[19]的最新研究方法达到五阶数值精度，但是由于扩散项差分近似和该方法对于复杂问题的计算比较困难等原因使得本文提出的方法在二维情况下只能达到四阶精度。

混合数值通量 ${\hat{q}}_{i+\frac{1}{2},j},{\hat{r}}_{i,j+\frac{1}{2}}$ 的差分公式如下

$\{\begin{array}{l}{\hat{q}}_{i+\frac{1}{2},j}=\frac{1}{64\Delta x}\left(3f_{i-1,j}-96f_{i,j}+93f_{i+1,j}+\Delta x{q}_{i-1,j}-12\Delta x{q}_{i,j}-15\Delta x{q}_{i+1,j}\right),\hfill \\ {\hat{r}}_{i,j+\frac{1}{2}}=\frac{1}{64\Delta y}\left(3g_{i,j-1}-96g_{i,j}+93g_{i,j+1}+\Delta y{r}_{i,j-1}-12\Delta y{r}_{i,j}-15\Delta y{r}_{i,j+1}\right).\hfill \end{array}$ (20)

混合导数 $\frac{{\partial }^{2}v}{\partial y^2},\frac{{\partial }^{2}w}{\partial x^2}$ 的数值近似如下

$\{\begin{array}{l}{\frac{{\partial }^{2}v}{\partial y^2}|}_{\left(x_i,y_j\right)}=\frac{3}{2\Delta y^3}\left(-5u_{i,j-1}+5u_{i,j+1}-\Delta y{v}_{i,j-1}-8\Delta y{v}_{i,j}-\Delta y{v}_{i,j+1}\right),\hfill \\ {\frac{{\partial }^{2}w}{\partial x^2}|}_{\left(x_i,y_j\right)}=\frac{3}{2\Delta x^3}\left(-5u_{i-1,j}+5u_{i+1,j}-\Delta x{w}_{i-1,j}-8\Delta x{w}_{i,j}-\Delta x{w}_{i+1,j}\right).\hfill \end{array}$ (21)

4. 数值算例

本节提供了上一节所构造的求解对流扩散方程数值方法的经典数值算例，以证明所提方法的准确性和有效性。一维和二维均匀网格点的数量分别用N和 $N_x\times N_y$ 表示。三阶TVD Runge-Kutta方法用于一维和二维数值方法的时间项离散，CFL条件数为0.6。

4.1. 一维线性方程精度测试

$\{\begin{array}{l}\frac{\partial u}{\partial t}+\frac{\partial u}{\partial x}=\epsilon \frac{{\partial }^{2}u}{\partial x^2},\hfill \\ u\left(x,0\right)=\text{sin}\left(x\right),0\le x\le 2\pi .\hfill \end{array}$ (22)

该问题选取周期边界条件。该问题的精确解为 $u\left(x,t\right)=\exp \left(-\epsilon t\right)\text{sin}\left(x-t\right)$ 。终止时刻设置为 $t=1$ 且 $\epsilon =0.01$ 。表1给出了用HWENO数值方法得到的数值精度。我们能够发现数值方法在一维情况下达到了五阶精度。

Table 1. One-dimensional accuracy test

表1. 一维精度测试

N	$L^1$ 误差	精度	$L^{\text{2}}$ 误差	精度	$L^{\infty }$ 误差	精度
10	2.84E−03		3.18E−03		4.73E−03
20	8.45E−04	5.07	9.87E−05	5.01	1.55E−04	4.93
40	2.58E−06	5.03	2.93E−06	5.07	4.55E−06	5.09
80	7.98E−08	5.02	8.89E−08	5.04	1.35E−07	5.07
160	2.47E−09	5.01	2.74E−09	5.02	3.97E−09	5.09
320	7.58E−11	5.02	8.41E−11	5.02	1.19E−10	5.06
640	2.42E−12	4.97	2.69E−12	4.97	3.78E−12	4.98

4.2. 一维Burgers方程

$\{\begin{array}{l}\frac{\partial u}{\partial t}+\frac{\partial \left(u^2/2\right)}{\partial x}=\epsilon \frac{{\partial }^{2}u}{\partial x^2},\hfill \\ u\left(x,0\right)=-\text{sin}\left(\pi x\right),-1\le x\le 1.\hfill \end{array}$ (23)

该问题选取Dirichlet边界条件，即有 $u\left(-1,t\right)=u\left(1,t\right)=0$ 。该问题没有直接的精确解，需要通过Hopf-Cole变换[21]来获得。扩散系数取为 $\epsilon =0.01/\pi$ ，终止时刻为 $t=\text{0}\text{.69}$ 。当 $N=100$ 时我们在图1中给出数值解和精确解的拟合图像。可以看见数值解较好地拟合了精确解，这意味着所设计的数值方法是有效的。

Figure 1. One-dimensional Burgers equation

图1. 一维Burgers方程

4.3. 一维Buckley-Leverett方程

$\frac{\partial u}{\partial t}+\frac{\partial f\left(u\right)}{\partial x}=\epsilon \frac{\partial \left(v\left(u\right)u_x\right)}{\partial x},\epsilon v\left(u\right)\ge 0.$ (24)

该算例的终止时刻设置为 $t=0.2$ ，扩散系数取为 $\epsilon =0.01$ 。且

$v\left(u\right)=\{\begin{array}{ll}4u\left(1-u\right),\hfill & 0\le u\le 1,\hfill \\ 0,\hfill & 其他.\hfill \end{array}$

该算例不考虑重力影响，故通量函数可以表示为

$f\left(u\right)=\frac{u^2}{u^2+{\left(1-u\right)}^2},$

至于守恒量u的连续初始条件为

$u\left(x,0\right)=\{\begin{array}{ll}1-3x,\hfill & 0\le x\le 1/3,\hfill \\ 0,\hfill & 1/3\le x\le 1.\hfill \end{array}$

空间计算区域 $\left[-1,1\right]$ 上的边界条件为 $u\left(0,t\right)=1,u\left(1,t\right)=0$ 。该问题目前尚未有精确的解析解，故我们用HWENO数值方法在网格点 $N=500$ 时得到的数值解作为参考的精确解。图2显示了 $N=100$ 的数值解与参考精确解之间的对比。我们发现即使在不连续处也没有振荡，数值解与研究中所述方法得出的参考解非常吻合。

4.4. 二维线性方程精度测试

$\{\begin{array}{l}\frac{\partial u}{\partial t}+\frac{\partial u}{\partial x}+\frac{\partial u}{\partial y}=\epsilon \left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right),\hfill \\ u\left(x,y,0\right)=\text{sin}\left(2\pi \left(x+y\right)\right).\hfill \end{array}$ (25)

精确解为 $u\left(x,y,t\right)=\exp \left(-8{\pi }^2\epsilon t\right)\text{sin}\left(2\pi \left(x+y-t\right)\right)$ 。此时计算域 $\left(x,y\right)\in \left[0,1\right]\times \left[0,1\right]$ 上的所有边界均选取周期边界条件。终止时刻设置为 $t=0.1$ 且扩散系数取为 $\epsilon =0.001$ 。表2列出了使用修正的HWENO方法得到的数值结果。可以看出，所构造的数值方法达到了四阶精度。

Figure 2. One-dimensional Buckley-Leverett equation

图2. 一维Buckley-Leverett方程

Table 2. Two-dimensional accuracy test

表2. 二维精度测试

$N_x\times N_y$	$L^1$ 误差	精度	$L^{\text{2}}$ 误差	精度	$L^{\infty }$ 误差	精度
10 * 10	3.43E−03		1.08E−02		3.42E−02
20 * 20	8.15E−05	5.39	3.65E−04	4.89	1.63E−03	4.39
40 * 40	2.75E−06	4.89	1.74E−05	4.39	1.10E−04	3.89
80 * 80	1.22E−07	4.49	1.09E−06	3.99	9.75E−06	3.49
160 * 160	5.09E−09	4.58	6.44E−08	4.08	8.14E−07	3.58
320 * 320	2.55E−10	4.32	4.55E−09	3.82	8.15E−08	3.32

4.5. 二维Buckley-Leverett方程

$\frac{\partial u}{\partial t}+\frac{\partial f\left(u\right)}{\partial x}+\frac{\partial g\left(u\right)}{\partial y}=\epsilon \left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right).$ (26)

扩散系数取 $\epsilon =0.01$ 。两个空间方向上的通量函数可以写为

$\begin{array}{l}f\left(u\right)=\frac{u^2}{u^2+{\left(1-u\right)}^2},\\ g\left(u\right)=f\left(u\right)\left(1-5{\left(1-u\right)}^2\right).\end{array}$

初始条件为

$u\left(x,y,0\right)=\{\begin{array}{ll}1,\hfill & x^2+y^2<0.5,\hfill \\ 0,\hfill & 其他.\hfill \end{array}$

需要注意的是，该算例在y方向上有重力影响。由于二维情况下的问题比较复杂，故目前没有解析解。矩形空间域 $\left[-1.5,1.5\right]\times \left[-1.5,1.5\right]$ 被 $80\times 80$ 等分，在计算域上满足周期边界条件。图3给出了数值方法在终止时刻 $t=0.\text{5}$ 时的数值模拟结果。可以发现等高线图的结果与参考文献[22]中的结果一致。

Figure 3. Two-dimensional Buckley-Leverett equation

图3. 二维Buckley-Leverett方程

5. 总结

在本文中，在有限差分框架中设计了一种修正的高阶HWENO方法，用于求解一维和二维的对流扩散方程。在所提出的数值方法中，对流项用修正的HWENO重构，扩散项用中心差分公式近似，时间项用Runge-Kutta方法近似。修正的HWENO重构直接使用解的点值及其一阶导数，这提高了计算效率。此外，可以选择任意的单调数值通量。为了保持扩散项的紧凑性，我们使用Hermite插值推出的差分公式对其进行近似。结合方程的特点，对二维情况下的高阶数值通量和混合导数进行了适当的处理，以保证方法的有效性和紧凑性。最后的数值算例验证了所提方法的高级精度和守恒性等性质。

基金项目

该工作得到了中国山西省自然科学基金(第202103021224041号)的部分支持。

NOTES

^*通讯作者。

参考文献

[1]	Liu, H. and Qiu, J. (2015) Finite Difference Hermite WENO Schemes for Conservation Laws, II: An Alternative Approach. Journal of Scientific Computing, 66, 598-624. https://doi.org/10.1007/s10915-015-0041-4
[2]	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
[3]	Harten, A. (1983) High Resolution Schemes for Hyperbolic Conservation Laws. Journal of Computational Physics, 49, 357-393. https://doi.org/10.1016/0021-9991(83)90136-5
[4]	Harten, A., Engquist, B., Osher, S. and Chakravarthy, S.R. (1987) Uniformly High Order Accurate Essentially Non-Oscillatory Schemes, III. Journal of Computational Physics, 71, 231-303. https://doi.org/10.1016/0021-9991(87)90031-3
[5]	Harten, A. (1987) Preliminary Results on the Extension of Eno Schemes to Two-Dimensional Problems. In: Morel, J.-M., Teissier, B., et al., Eds., Lecture Notes in Mathematics, Springer, 23-40. https://doi.org/10.1007/bfb0078315
[6]	Casper, J. (1992) Finite-Volume Implementation of High-Order Essentially Nonoscillatory Schemes in Two Dimensions. AIAA Journal, 30, 2829-2835. https://doi.org/10.2514/3.11625
[7]	Casper, J. and Atkins, H.L. (1993) A Finite-Volume High-Order ENO Scheme for Two-Dimensional Hyperbolic Systems. Journal of Computational Physics, 106, 62-76. https://doi.org/10.1006/jcph.1993.1091
[8]	Shu, C.-W. and Osher, S. (1989) Efficient Implementation of Essentially Non-Oscillatory Shock-Capturing Schemes, II. Journal of Computational Physics, 83, 32-78. https://doi.org/10.1016/0021-9991(89)90222-2
[9]	Liu, X.-D., Osher, S. and Chan, T. (1994) Weighted Essentially Non-Oscillatory Schemes. Journal of Computational Physics, 115, 200-212. https://doi.org/10.1006/jcph.1994.1187
[10]	Jiang, G.-S. and Shu, C.-W. (1996) Efficient Implementation of Weighted ENO Schemes. Journal of Computational Physics, 126, 202-228. https://doi.org/10.1006/jcph.1996.0130
[11]	Shu, C.-W. (2009) High Order Weighted Essentially Nonoscillatory Schemes for Convection Dominated Problems. SIAM Review, 51, 82-126. https://doi.org/10.1137/070679065
[12]	Qiu, J.X. and Shu, C.-W. (2004) Hermite WENO Schemes and Their Application as Limiters for Runge-Kutta Discontinuous Galerkin Method: One-Dimensional Case. Journal of Computational Physics, 193, 115-135. https://doi.org/10.1016/j.jcp.2003.07.026
[13]	Qiu, J.X. and Shu, C.-W. (2005) Hermite WENO Schemes and Their Application as Limiters for Runge-Kutta Discontinuous Galerkin Method II: Two Dimensional Case. Computers & Fluids, 34, 642-663. https://doi.org/10.1016/j.compfluid.2004.05.005
[14]	Zhang, M. and Zhao, Z. (2023) A Fifth-Order Finite Difference HWENO Scheme Combined with Limiter for Hyperbolic Conservation Laws. Journal of Computational Physics, 472, Article 111676. https://doi.org/10.1016/j.jcp.2022.111676
[15]	Liu, H.X. and Qiu, J.X. (2014) Finite Difference Hermite WENO Schemes for Hyperbolic Conservation Laws. Journal of Scientific Computing, 63, 548-572. https://doi.org/10.1007/s10915-014-9905-2
[16]	Zhao, Z., Zhang, Y.-T. and Qiu, J.X. (2020) A Modified Fifth Order Finite Difference Hermite WENO Scheme for Hyperbolic Conservation Laws. Journal of Scientific Computing, 85, Article No. 29. https://doi.org/10.1007/s10915-020-01347-1
[17]	Ma, Z. and Wu, S.-P. (2018) HWENO Schemes Based on Compact Differencefor Hyperbolic Conservation Laws. Journal of Scientific Computing, 76, 1301-1325. https://doi.org/10.1007/s10915-018-0663-4
[18]	Tao, Z.J., Li, F.Y. and Qiu, J.X. (2016) High-Order Central Hermite WENO Schemes: Dimension-by-Dimension Moment-Based Reconstructions. Journal of Computational Physics, 318, 222-251. https://doi.org/10.1016/j.jcp.2016.05.005
[19]	Zhao, Z., Chen, Y.B. and Qiu, J.X. (2020) A Hybrid Hermite WENO Scheme for Hyperbolic Conservation Laws. Journal of Computational Physics, 405, Article 109175. https://doi.org/10.1016/j.jcp.2019.109175
[20]	Zhu, J. and Qiu, J.X. (2008) A Class of the Fourth Order Finite Volume Hermite Weighted Essentially Non-Oscillatory Schemes. Science in China Series A: Mathematics, 51, 1549-1560. https://doi.org/10.1007/s11425-008-0105-0
[21]	Basdevant, C., Deville, M., Haldenwang, P., Lacroix, J.M., Ouazzani, J., Peyret, R., et al. (1986) Spectral and Finite Difference Solutions of the Burgers Equation. Computers & Fluids, 14, 23-41. https://doi.org/10.1016/0045-7930(86)90036-8
[22]	Kurganov, A. and Tadmor, E. (2000) New High-Resolution Central Schemes for Nonlinear Conservation Laws and Convection-Diffusion Equations. Journal of Computational Physics, 160, 241-282. https://doi.org/10.1006/jcph.2000.6459