abcd Boussinesq系统的高阶中心间断Galerkin-有限元方法

High Order Central Discontinuous Galerkin-Finite Element Methods for the abcd Boussinesq System

摘要: Boussinesq类型的方程是一类重要的非线性色散浅水波方程，广泛应用于研究小振幅浅水波的传播现象。我们考虑一个含有四个参数的Boussinesq系统：abcd Boussinesq系统，并针对该系统设计了一个高阶中心间断Galerkin-有限元方法。该方法首先将abcd Boussinesq系统改写为守恒律方程与椭圆型方程的一个耦合系统，然后采用中心局部间断Galerkin方法求解守恒律方程，使用有限元方法求解椭圆型方程，最后我们通过一系列的数值算例验证所提出方法的高阶精度和有效性。

Abstract: Boussinesq type equations are an important class of nonlinear dispersive shallow water wave equa-tions, it has been widely studied to model the propagation phenomena of small amplitude shallow water waves. In this work, we consider a Boussinesq system with four parameters: abcd Boussinesq system, and design a high order central discontinuous Galerkin-finite element methods for solving this system. In our numerical approach, we first reformulate the abcd Boussinesq system into a coupled system of conservation law and elliptic equations. Then we propose a family of high order numerical methods which discretize the conservation law with central local discontinuous Galerkin methods and the elliptic equations with continuous finite element methods. Different types of nu-merical tests are provided to illustrate the accuracy and effectiveness of the proposed schemes.

1. 引言

Boussinesq方程是一类重要的非线性偏微分方程，它是模拟近海岸浅水波传播现象的重要工具之一。关于水波的Boussinesq近似，最初是由Boussinesq [1] 在1871年提出，从欧拉方程出发推导出了第一个Boussinesq方程。此后，不同类型的Boussinesq方程相继被提出，例如“good” Boussinesq方程 [2] 与improved Boussinesq方程 [3] 。直到2002年，Bona等人 [4] 提出了一个更一般化的Boussinesq方程，即abcd Boussinesq系统：

$\{\begin{array}{l}{\eta }_t+u_x+{\left(\eta u\right)}_x+a{u}_{xxx}-b{\eta }_{xxt}=0,\\ u_t+{\eta }_x+{\left(\frac{u^2}{2}\right)}_x+c{\eta }_{xxx}-d{u}_{xxt}=0,\end{array}$ (1)

这里， $\eta$ 表示水的自由表面与其静止位置的垂直偏差，u表示水平方向上的速度，a，b，c，d是四个常数，满足以下的关系式：

$a+b=\frac{1}{2}\left({\theta }^2-\frac{1}{3}\right)$ ， $c+d=\frac{1}{2}\left(1-{\theta }^2\right)\ge 0$ ， $\theta \in \left[0,1\right]$ 。 (2)

四个参数a，b，c，d的不同取值对应着不同的浅水波模型。例如，当 $a=b=c=0,d=1/3$ 时方程(1)对应着经典的Boussinesq方程， $a=c=0,b=d=1/6$ 时相应地变成了耦合的Benjamin-Bona-Mahony (BBM)方程， $a=c=1/6,b=d=0$ 时也就是耦合的Korteweg-de Vries (KdV)方程。

目前关于Boussinesq类型方程的数值方法已经得到了很好的发展。例如Peregrine等人 [5] 在1966年针对经典的Boussinesq方程提出了一个有限差分近似方法。Bona等人 [6] 在2007年采用带有周期光滑样条基函数的Galerkin方法数值求解了耦合的KdV-KdV系统。2019年，Burtea和Courtès [7] 关于abcd Boussinesq系统提出了一个完全离散的有限体积方法，并证明了该方法在时间方向上的一阶收敛速度和空间方向上的二阶收敛速度。2022年，Jiawei Sun，Shusen Xie与Yulong Xing [8] 针对abcd Boussinesq系统发展了一个局部间断Galerkin方法，通过选取不同的数值通量，对于四个参数a，b，c，d不同的取值情况，给出了一个严格的先验误差估计，并通过一系列的数值实验验证了理论结果。采用数值方法数值求解abcd Boussinesq系统一般会遇到以下两个问题：一、abcd Boussinesq系统是一个高阶的非线性偏微分方程组，包含二阶、三阶偏导数，所以不能采用一般的低阶数值方法去数值求解该偏微分方程组。二、方程组里面包含有对时间变量的导数以及对时间变量与空间变量的混合偏导数，一样成为了数值求解该偏微分方程组的一个难点。

中心间断Galerkin方法 [9] [10] [11] 是一种基于中心方法 [12] 与间断Galerkin方法 [13] [14] [15] 发展出来的一种新型数值计算方法，该方法定义在两套重叠网格上，交错的计算两组近似解，从而避免了间断Galerkin方法里面需要定义数值通量的问题，因此它很适用于方程特征系统未知的问题。本文的基本结构如下：首先将原始的abcd Boussinesq系统改写为守恒律方程与椭圆型方程的一个耦合系统，使得该耦合系统中不再含有对时间变量以及空间变量的混合导数项。然后采用中心局部间断Galerkin方法求解守恒律方程，使用连续有限元方法求解椭圆型方程，最后我们通过一系列的数值实验(精度测试、行波碰撞)验证所提出数值方法的高阶精度和有效性。

2. 数值方法

2.1. 改写abcd Boussinesq系统

为了消除abcd Boussinesq系统(1)中的混合时间与空间导数项，引入两个人工变量

$\{\begin{array}{l}p=\eta -b{\eta }_{xx},\\ q=u-d{u}_{xx},\end{array}$ (3)

那么方程组(1)就可以被改写为

$\{\begin{array}{l}{p}_t+{\left(u+\eta u+a{u}_{xx}\right)}_x=0,\\ q_t+{\left(\eta +\frac{u^2}{2}+c{\eta }_{xx}\right)}_x=0.\end{array}$ (4)

此时，方程组(4)中通量里面依然含有两个二阶导数项，所以为了将(4)进一步改写为一个一阶系统，再引入四个人工变量

$\{\begin{array}{l}\xi ={\eta }_x,\\ \gamma ={\xi }_x={\eta }_{xx},\\ v=u_x,\\ w=v_x=u_{xx},\end{array}$ (5)

至此，就得到了最终的一阶系统：

$\{\begin{array}{l}{p}_t+f_1{\left(\eta ,u,w\right)}_x=0,\\ q_t+f_2{\left(\eta ,u,\gamma \right)}_x=0,\end{array}$ (6)

其中 $f_1\left(\eta ,u,w\right)=u+\eta u+aw$ ， $f_2\left(\eta ,u,\lambda \right)=\eta +\frac{u^2}{2}+c\gamma$ 是通量函数。因此数值求解abcd Boussinesq系统(1)就等价于数值求解方程组(3)、(5)、(6)。

2.2. 中心间断Galerkin-有限元方法

假设 ${\left\{x_j\right\}}_j$ 是计算区域 $\Omega =\left[x_{\min },x_{\max }\right]$ 的一个网格剖分，定义 $x_{j+\frac{1}{2}}=\frac{1}{2}\left(x_j+x_{j+1}\right)$ ， $I_j=\left(x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}\right)$ ， $I_{j+\frac{1}{2}}=\left(x_j,x_{j+1}\right)$ ，则 ${\left\{I_{j+\frac{1}{2}}\right\}}_j$ 是与 ${\left\{I_j\right\}}_j$ 相互重叠的一套网格。基于以上的两套网格，定义以下的四个有限维离散函数空间：

，

，

，

。

为了简单起见，只给出向前Euler时间离散方法的数值格式，后面将讨论高阶时间离散的问题。所提及的数值方法涉及到两套解，假设在 $t=t_n$ 时刻 ${\eta }^{n,\ast }$ 与 $u^{n,\ast }$ 是已知的，我们需要去计算在 $t=t_{n+1}=t_n+\Delta t_n$ 时刻的解 ${\eta }^{n+1,\ast }$ 和 $u^{n+1,\ast }$ ，这里符号 $\star$ 表示C或者D。为了简便，只给出计算 ${\eta }^{n+1,C}$ 和 $u^{n+1,C}$ 的计算公式，关于 ${\eta }^{n+1,D}$ 和 $u^{n+1,D}$ 的计算公式类似可得。

已知 ${\eta }^{n,\ast }$ 与 $u^{n,\ast }$ ，为了得到 $p^{n+1,C}$ 和 $q^{n+1,C}$ ，应用中心间断Galerkin方法作为空间离散方法，向前欧拉法作为时间离散方法去离散方程组(5)、(6)，即寻找 ${\xi }^{n,C}$ ， ${\gamma }^{n,C}$ ， $v^{n,C}$ ， $w^{n,C}$ ， $p^{n+1,C}$ 和 $q^{n+1,C}\in {\mathcal{V}}^{C,k}$ ，使得对于任意的 $\varsigma$ ， $\chi$ ， $\bar{\varsigma }$ ， $\bar{\chi }$ ， $\phi$ ， $\varphi \in {{\mathcal{V}}^{C,k}|}_{I_j}$ 有

${\displaystyle {\int }_{I_j}{\xi }^{n,C}}\varsigma \text{d}x=-{\displaystyle {\int }_{I_j}{\eta }^{n,D}{\varsigma }_x}\text{d}x+{\eta }^{n,D}\left(x_{j+\frac{1}{2}}\right)\varsigma \left(x_{j+\frac{1}{2}}^{-}\right)-{\eta }^{n,D}\left(x_{j-\frac{1}{2}}\right)\varsigma \left(x_{j-\frac{1}{2}}^{+}\right)$ ， (7)

${\displaystyle {\int }_{I_j}{\gamma }^{n,C}}\chi \text{d}x=-{\displaystyle {\int }_{I_j}{\xi }^{n,D}{\chi }_x}\text{d}x+{\xi }^{n,D}\left(x_{j+\frac{1}{2}}\right)\chi \left(x_{j+\frac{1}{2}}^{-}\right)-{\xi }^{n,D}\left(x_{j-\frac{1}{2}}\right)\chi \left(x_{j-\frac{1}{2}}^{+}\right)$ ， (8)

${\displaystyle {\int }_{I_j}{v}^{n,C}}\bar{\varsigma }\text{d}x=-{\displaystyle {\int }_{I_j}{u}^{n,D}{\bar{\varsigma }}_x}\text{d}x+u^{n,D}\left(x_{j+\frac{1}{2}}\right)\bar{\varsigma }\left(x_{j+\frac{1}{2}}^{-}\right)-u^{n,D}\left(x_{j-\frac{1}{2}}\right)\bar{\varsigma }\left(x_{j-\frac{1}{2}}^{+}\right)$ ， (9)

${\displaystyle {\int }_{I_j}{w}^{n,C}}\bar{\chi }\text{d}x=-{\displaystyle {\int }_{I_j}{v}^{n,D}{\bar{\chi }}_x}\text{d}x+v^{n,D}\left(x_{j+\frac{1}{2}}\right)\bar{\chi }\left(x_{j+\frac{1}{2}}^{-}\right)-v^{n,D}\left(x_{j-\frac{1}{2}}\right)\bar{\chi }\left(x_{j-\frac{1}{2}}^{+}\right)$ ， (10)

$\begin{array}{c}{\displaystyle {\int }_{I_j}{p}^{n+1,C}}\phi \text{d}x={\displaystyle {\int }_{I_j}\left({\theta }_n{p}^{n,D}+\left(1-{\theta }_n\right)p^{n,C}\right)\phi \text{d}x}+\Delta t_n{\displaystyle {\int }_{I_j}{f}_1}\left({\eta }^{n,D},u^{n,D},w^{n,D}\right){\phi }_x\text{d}x\\ -\Delta t_n[f_1\left({\eta }^{n,D}\left(x_{j+\frac{1}{2}}\right),u^{n,D}\left(x_{j+\frac{1}{2}}\right),w^{n,D}\left(x_{j+\frac{1}{2}}\right)\right)\phi \left(x_{j+\frac{1}{2}}^{-}\right)\\ -f_1\left({\eta }^{n,D}\left(x_{j-\frac{1}{2}}\right),u^{n,D}\left(x_{j-\frac{1}{2}}\right),w^{n,D}\left(x_{j-\frac{1}{2}}\right)\right)\phi \left(x_{j-\frac{1}{2}}^{+}\right)],\end{array}$ (11)

$\begin{array}{c}{\displaystyle {\int }_{I_j}{q}^{n+1,C}}\varphi \text{d}x={\displaystyle {\int }_{I_j}\left({\theta }_n{q}^{n,D}+\left(1-{\theta }_n\right)q^{n,C}\right)\varphi \text{d}x}+\Delta t_n{\displaystyle {\int }_{I_j}{f}_2}\left({\eta }^{n,D},u^{n,D},{\gamma }^{n,D}\right){\varphi }_x\text{d}x\\ -\Delta t_n[f_2\left({\eta }^{n,D}\left(x_{j+\frac{1}{2}}\right),u^{n,D}\left(x_{j+\frac{1}{2}}\right),{\gamma }^{n,D}\left(x_{j+\frac{1}{2}}\right)\right)\varphi \left(x_{j+\frac{1}{2}}^{-}\right)\\ -f_2\left({\eta }^{n,D}\left(x_{j-\frac{1}{2}}\right),u^{n,D}\left(x_{j-\frac{1}{2}}\right),{\gamma }^{n,D}\left(x_{j-\frac{1}{2}}\right)\right)\varphi \left(x_{j-\frac{1}{2}}^{+}\right)],\end{array}$ (12)

这里 ${\theta }_n=\Delta t_n/\tau \in \left[0,1\right]$ ， $\tau$ 是由CFL条件所限制的时间步长的一个上限，对于任意的函数 $\omega$ ， $\omega \left(x^{\pm }\right)=\underset{ϵ\to 0+}{\lim }\omega \left(x\pm ϵ\right)$ 。

一旦有了 $p^{n+1,C}$ 和 $q^{n+1,C}$ ，就可以通过连续有限元方法求解方程组(3)得到 ${\eta }^{n+1,C}$ 和 $u^{n+1,C}$ ，也就是寻找 ${\eta }^{n+1,C}$ ， $u^{n+1,C}\in {\tilde{W}}^{C,k}$ 使得对于任意的 $\psi$ ， $\bar{\psi }\in {\hat{W}}^{C,k}$ 有

${\displaystyle {\int }_{\Omega }{\eta }^{n+1,C}\psi \text{d}x}+b{\displaystyle {\int }_{\Omega }{\left({\eta }^{n+1,C}\right)}_x{\psi }_x\text{d}x}={\displaystyle {\int }_{\Omega }{p}^{n+1,C}}\psi \text{d}x$ ， (13)

${\displaystyle {\int }_{\Omega }{u}^{n+1,C}\bar{\psi }\text{d}x}+d{\displaystyle {\int }_{\Omega }{\left(u^{n+1,C}\right)}_x{\bar{\psi }}_x\text{d}x}={\displaystyle {\int }_{\Omega }{q}^{n+1,C}}\bar{\psi }\text{d}x$ ， (14)

这里 ${\tilde{W}}^{C,k}$ ， ${\hat{W}}^{C,k}$ 是 ${\mathcal{W}}^{C,k}$ 根据问题边界条件变化而得的连续函数空间，例如，当 $\eta$ ，u的边界条件是Dirichlet类型时，有

${\tilde{W}}^{C,k}=\left\{w\in {\mathcal{W}}^{C,k}:w与\eta ,u有相同的边界条件\right\}$ ，

${\hat{W}}^{C,k}=\left\{w\in {\mathcal{W}}^{C,k}:w有零边界条件\right\}$ ；

当 $\eta$ ，u的边界条件是周期边界类型时，有

${\tilde{W}}^{C,k}={\hat{W}}^{C,k}=\left\{w\in {\mathcal{W}}^{C,k}:w\left(x_{\min }\right)=w\left(x_{\max }\right)\right\}$ ；

当 $\eta$ ，u的边界条件是Neumann类型( $\frac{\partial \eta }{\partial n}=\frac{\partial u}{\partial n}=0$ )时，有

${\tilde{W}}^{C,k}={\hat{W}}^{C,k}={\mathcal{W}}^{C,k}$ 。

2.3. 高阶时间离散方法

上面部分为了表述的方便，仅仅采用了向前欧拉时间离散方法，然而为了获得时间上的更高阶精度以及保持数值算法的稳定性，在实际计算过程中通常采用保持高度稳定(Strong stability preserving (SSP))的高阶Runge-Kutta时间离散方法 [16] ，该方法可以写成向前欧拉时间离散方法的一个凸组合。例如，对于求解方程

$u_t=F\left(u\right)$ ， (15)

这里， $F\left(u\right)$ 是一个空间离散算子，可以采用下面的三阶TVD Runge-Kutta时间离散方法：

$\begin{array}{l}{u}^{\left(1\right)}=u^n+\Delta tF\left(u^n\right),\\ u^{\left(2\right)}=\frac{3}{4}{u}^n+\frac{1}{4}\left(u^{\left(1\right)}+\Delta tF\left(u^{\left(1\right)}\right)\right),\\ u^{n+1}=\frac{1}{3}{u}^n+\frac{2}{3}\left(u^{\left(2\right)}+\Delta tF\left(u^{\left(2\right)}\right)\right).\end{array}$

3. 数值算例

这一部分，我们提出一系列的数值算例去验证所提出数值方法的精确性以及可靠性。首先通过四个参数a，b，c，d取不同的值，研究中心间断Galerkin-有限元方法的计算精度，然后通过模拟有限时间内的爆破波、孤立波的迎面碰撞研究中心间断Galerkin-有限元方法的可靠性。计算过程中，我们使用三阶TVD Runge-Kutta时间离散方法去离散时间变量，时间步长被动态的决定

$\Delta t_n=\frac{C_{cfl}\min \left(\Delta x_i\right)}{\max \left(\left|u\right|\right)}$ ，

当速度为零时，我们使用固定的时间步长： $\Delta t=C_{cfl}\min \left(\Delta x_i\right)$ 。空间上，我们使用( $k=1$ )线性多项式以及( $k=2$ )二次多项式近似，一般来说，对于 $k=1$ 的情况，取 $C_{cfl}=0.03$ ， $k=2$ 时取 $C_{cfl}=0.015$ 。

3.1. 精度测试一

分别取四个参数 $a=-\frac{7}{30}$ ， $b=\frac{7}{15}$ ， $c=-\frac{2}{5}$ ， $d=\frac{1}{2}$ ，此时方程组的准确解为 [17]

$\{\begin{array}{l}\eta \left(x,t\right)=\frac{3}{8}{\mathrm{sech}}^2\left(\frac{1}{2}\sqrt{\frac{5}{7}}\left(x-20-\frac{5\sqrt{2}}{6}t\right)\right),\\ u\left(x,t\right)=\frac{1}{2\sqrt{2}}{\mathrm{sech}}^2\left(\frac{1}{2}\sqrt{\frac{5}{7}}\left(x-20-\frac{5\sqrt{2}}{6}t\right)\right),\end{array}$ (16)

取计算区域为 $\Omega =\left[0,40\right]$ ，计算时间为 $T=0.8$ 。 $\eta$ 和u在 $L^1$ 、 $L^2$ 、 $L^{\infty }$ 范数意义下的误差和收敛阶分别如表1和表2所示，从表中可以看出对于空间k次多项式的近似，所提出的中心间断Galerkin-有限元方法具有最优的 $\left(k+1\right)$ 阶收敛速度。

Table 1. L 1 , L 2 , L ∞ errors and orders of accuracy of η

表1. $\eta$ 的 $L^1$ 、 $L^2$ 、 $L^{\infty }$ 误差和相应的收敛阶

Table 2. L 1 , L 2 , L ∞ errors and orders of accuracy of u

表2. u的 $L^1$ 、 $L^2$ 、 $L^{\infty }$ 误差和相应的收敛阶数

3.2. 精度测试二

取四个参数 $a=\frac{1}{6}$ ， $b=0$ ， $c=\frac{1}{6}$ ， $d=0$ ，此时的abcd Boussinesq系统退化为耦合的KdV-KdV方程组，它的准确解为

$\{\begin{array}{l}\eta \left(x,t\right)=-1+\frac{3}{2}{\mathrm{sech}}^2\left(\sqrt{\frac{3}{2}}\left(x-20-\sqrt{2}t\right)\right),\\ u\left(x,t\right)=\frac{3}{\sqrt{2}}{\mathrm{sech}}^2\left(\sqrt{\frac{3}{2}}\left(x-20-\sqrt{2}t\right)\right),\end{array}$ (17)

取计算区域为 $\Omega =\left[0,40\right]$ ，计算时间为 $T=0.01$ 。 $\eta$ 和u在 $L^1$ 、 $L^2$ 、 $L^{\infty }$ 范数意义下的误差和收敛阶分别如表3和表4所示，同样从表中可以看出对于空间k次多项式的近似，我们提出的中心间断Galerkin-有限元方法具有最优的 $\left(k+1\right)$ 阶收敛速度。

Table 3. L 1 , L 2 , L ∞ errors and orders of accuracy of η

表3. $\eta$ 的 $L^1$ 、 $L^2$ 、 $L^{\infty }$ 误差和相应的收敛阶数

Table 4. L 1 , L 2 , L ∞ errors and orders of accuracy of u

表4. u的 $L^1$ 、 $L^2$ 、 $L^{\infty }$ 误差和相应的收敛阶数

3.3. 精度测试三

取四个参数 $a=0$ ， $b=\frac{1}{3}$ ， $c=-\frac{1}{3}$ ， $d=\frac{1}{3}$ ，此时abcd Boussinesq系统的准确解为

$\{\begin{array}{l}\eta \left(x,t\right)=-1,\\ u\left(x,t\right)=1+6{\mathrm{sech}}^2\left(\frac{1}{\sqrt{2}}\left(x-20-3t\right)\right),\end{array}$ (18)

依然取计算区域为 $\Omega =\left[0,40\right]$ ，计算时间为 $T=0.01$ 。 $\eta$ 和u在 $L^1$ 、 $L^2$ 、 $L^{\infty }$ 范数意义下的误差和收敛阶分别如表5和表6所示，从表中可以看出对于空间k次多项式的近似，所给出的中心间断Galerkin-有限元方法在 $\eta$ 上面保持机器精度，在u上面具有最优的 $\left(k+1\right)$ 阶收敛速度。

Table 5. L 1 ， L 2 ， L ∞ errors and orders of accuracy of η

表5. $\eta$ 的 $L^1$ 、 $L^2$ 、 $L^{\infty }$ 误差和相应的收敛阶数

Table 6. L 1 ， L 2 ， L ∞ errors and orders of accuracy of u

表6. u的 $L^1$ 、 $L^2$ 、 $L^{\infty }$ 误差和相应的收敛阶数

3.4. 有限时间爆破

这一个算例里面，取四个参数 $a=0$ ， $b=\frac{1}{6}$ ， $c=0$ ， $d=\frac{1}{6}$ ，此时abcd Boussinesq系统退化为耦合的BBM方程组，并具有以下形式的行波解 [7]

$\{\begin{array}{l}{\eta }_{\pm }\left(x,t\right)=\frac{15}{2}{\mathrm{sech}}^2\left(\frac{3}{\sqrt{10}}\left(x\pm \frac{5}{2}t\right)\right)-\frac{45}{4}{\mathrm{sech}}^4\left(\frac{3}{\sqrt{10}}\left(x\pm \frac{5}{2}t\right)\right),\\ u_{\pm }\left(x,t\right)=\mp \frac{15}{2}{\mathrm{sech}}^2\left(\frac{3}{\sqrt{10}}\left(x\pm \frac{5}{2}t\right)\right).\end{array}$ (19)

取初始条件为

$\begin{array}{l}\eta \left(x,0\right)={\eta }_{+}\left(x-x_{+},0\right)+{\eta }_{-}\left(x-x_{-},0\right),\\ u\left(x,0\right)=u_{+}\left(x-x_{+},0\right)+u_{-}\left(x-x_{-},0\right),\end{array}$ (20)

式中 $x_{\pm }=\pm 7$ ，计算区域取为 $\Omega =\left[-14,14\right]$ ，计算时间为 $T=4.5$ 。 $\eta$ 和u在爆破前不同时刻的数值解如图1所示。

(a) (b) (c) (d)

Figure 1. The numerical solutions of $\eta$ and u at different times (before the blow-up)

图1. $\eta$ 和u在爆破前不同时刻的数值解

3.5. 孤立波的迎面碰撞

最后一个算例，考虑孤立波的迎面碰撞，取四个参数 $a=-\frac{7}{30}$ ， $b=\frac{7}{15}$ ， $c=-\frac{2}{5}$ ， $d=\frac{1}{2}$ ，初始条件为

$\begin{array}{l}\eta \left(x,0\right)={\eta }_{+}\left(x\right)+{\eta }_{-}\left(x\right),\\ u\left(x,0\right)=u_{+}\left(x\right)+u_{-}\left(x\right),\end{array}$ (21)

其中

$\{\begin{array}{l}{\eta }_{\pm }\left(x\right)=\frac{1}{\sqrt{8}}{\mathrm{sech}}^2\left(\sqrt{\frac{5}{28}}\left(x-x_{\pm }\right)\right),\\ u_{\pm }\left(x\right)=\pm \frac{3}{8}{\mathrm{sech}}^2\left(\sqrt{\frac{5}{28}}\left(x+x_{\pm }\right)\right),\end{array}$ (22)

$x_{\pm }=\pm 7$ 。计算区域取为 $\Omega =\left[-14,14\right]$ ， $\eta$ 和u在迎面碰撞前、后不同时刻的数值解如图2所示。

(a) (b) (c) (d) (e)

Figure 2. The numerical solutions of $\eta$ and u at different times (head-on collision)

图2. $\eta$ 和u在迎面碰撞前后不同时刻的数值解

4. 结论

本文设计了一个数值求解abcd Boussinesq系统的高效数值方法，该方法首先通过引入人工变量来改写最初的模型，然后再将中心间断Galerkin-有限元方法应用到改写后的方程组，最后通过一系列的数值算例验证了该方法的高阶精度和可靠性。那么基于本文的工作，可以考虑将此类方法应用于更加复杂的非线性、色散浅水波模型的数值模拟。

基金项目

重庆市教育委员会科学技术研究项目(KJQN202200725)。

NOTES

^*通讯作者。