求解具有间断系数波动方程的浸入有限元方法

Immersed Finite Element Method for Solving Wave Equations with Discontinuous Coefficients

摘要: 数值求解具有间断系数的波动方程在实际工程中有着非常广泛的应用。本文针对具有间断系数的波动方程，提出一种基于界面非拟合网格的浸入有限元方法。该方法具有无需人为调整罚参数的优点。本文给了一系列数值试验，实验结果验证了该方法的有效性。

Abstract: The numerical solution of wave equation with discontinuous coefficient is widely used in practical engineering. In this paper, an immersed finite element method based on interface unfitted meshes is proposed for wave equations with discontinuous coefficients. This method has the advantage of not needing to adjust penalty parameters artificially. In this paper, a series of numerical experiments are given, and the experimental results verify the effectiveness of the method.

1. 引言

波动方程用来描述各种波动现象，在科学和工程学中有重要应用。在实际问题中，所考虑的计算区域通常由许多不同介质组成，波在不同介质中的传播系数不同，在介质的交界面上需要满足一定的界面条件，这便产生了界面问题。具有间断系数的波动方程模拟了波在充满不同均匀介质的界面上的传播，该问题可以用来研究声波在地球上的传播，在不同地质材料的界面处，材料性质会发生突然变化，具有间断系数的波动方程可以用作地球地震建模的简化模型[1]。还有许多应用包括声学透镜或波导的研究、海洋声学的研究，以及超声波成像问题等。求解具有间断系数的波动方程，不仅需要考虑通常的初始条件和边界条件，还需要考虑界面上满足的跳跃条件，这使得方程的解析求解十分困难，因此设计一个快速有效的数值求解方法对于解决实际界面问题非常重要。

求解界面问题的一种传统的做法是采用界面拟合网格离散[2]，但该方法在应用中通常会受到限制，即网格中的每一个单元必须包含在一种介质中，任何一个单元都不能有界面从中穿过。这使得该方法在处理移动界面时需要多次重新划分网格，这大大的增加了存储空间和计算量。为了克服这样的缺点，学者们提出了一些与界面无关的非拟合网格数值计算方法，浸入有限元方法就是其中之一[3]。

浸入有限元方法通过在界面单元上根据界面跳跃条件设计分片多项式，在非界面单元上使用标准的多项式有限元函数来解决非拟合网格上的界面问题。该方法已应用于求解椭圆型问题[4] [5]、抛物型问题[4]以及双曲型问题[6] [7]等。由于浸入有限元空间在与界面相交的边上不连续，因此数值格式会产生相容误差，从而得不到最优收敛性。为了克服这一缺点，Lin等人在浸入有限元方法中引入了部分惩罚思想，在界面边上添加了惩罚项，提出了部分惩罚浸入有限元方法[4]。然而该方法在计算时需要人为选择足够大的罚参数来保证方法的稳定性，针对这一缺陷，Ji等人提出了改进的部分罚浸入有限元方法[5]。该方法通过构造提升算子来保证方法的稳定性，从而避免了人为选取罚参数的麻烦。

本文采用改进的无需罚参数的部分罚浸入有限元方法研究了具有间断系数的波动方程，提出了时间二阶精度的全离散数值格式，接着，通过一系列算例，验证了该方法的有效性。

2. 模型问题

设$\Omega$ 是$R^2$ 中的一个有界域，光滑曲线$\Gamma \subset \Omega$ 将$\Omega$ 分为${\Omega }^{-}$ 和${\Omega }^{+}$ 两个子区域，如图1所示，考虑以下具有间断系数的波动方程：

$\frac{{\partial }^{2}u}{\partial t^2}-\nabla \cdot \left(\beta \nabla u\right)=f\left(x,y,t\right),\left(x,y\right)\in \Omega ,t\in (0,T],$ (1)

$u\left(x,y,t\right)=0,\left(x,y\right)\in \partial \Omega ,t\in (0,T],$ (2)

$u\left(x,y,0\right)=u_0\left(x,y\right),\left(x,y\right)\in \Omega ,$ (3)

$\frac{\partial u}{\partial t}\left(x,y,0\right)=u_1\left(x,y\right),\left(x,y\right)\in \Omega ,$ (4)

其解满足跳跃条件

${\left[u\right]}_{\Gamma }=0,$ (5)

${\left[\beta \nabla u\cdot n_{\Gamma }\right]}_{\Gamma }=0.$ (6)

Figure 1. Region $\Omega ={\Omega }^{-}\cup {\Omega }^{+}\cup \Gamma$

图1. 区域$\Omega ={\Omega }^{-}\cup {\Omega }^{+}\cup \Gamma$

上式中$n_{\Gamma }$ 为界面$\Gamma$ 上从${\Omega }^{+}$ 指向${\Omega }^{-}$ 的单位法向量，${\left[\cdot \right]}_{\Gamma }$ 表示在界面$\Gamma$ 处的跳跃，即${\left[v\right]}_{\Gamma }={v^{+}|}_{\Gamma }-{v^{-}|}_{\Gamma }$ ，其中$v^{+}={v|}_{{\Omega }^{+}}$ ，$v^{-}={v|}_{{\Omega }^{-}}$ ，系数$\beta$ 在$\Omega$ 上为分片常数函数，即

$\beta \left(x,y\right)=\{\begin{array}{l}{\beta }^{+}>0,\left(x,y\right)\in {\Omega }^{+},\\ {\beta }^{-}>0,\left(x,y\right)\in {\Omega }^{-}.\end{array}$

3. 全离散浸入有限元方法

设${\mathcal{T}}_h$ 是$\Omega$ 上的一个一致三角剖分。对于任意的$K\in {\mathcal{T}}_h$ ，定义$h_K$ 是K的直径，网格参数$h={\max }_{K\in {\mathcal{T}}_h}{h}_K$ 。记${\mathcal{T}}_h^i=\left\{K|K\in {\mathcal{T}}_h,K\cap \Gamma \ne \varnothing \right\}$ 是界面单元的集合，${\mathcal{T}}_h^n=\left\{K|K\in {\mathcal{T}}_h,K\cap \Gamma =\varnothing \right\}$ 是非界面单元的集合。

Figure 2. Interface element

图2. 界面单元

下面以图2的界面单元$K=\Delta P_1{P}_2{P}_3$ 为例定义浸入有限元空间。设界面$\Gamma$ 与K的边交于DE两点，点$P_1\in {\Omega }^{-}$ ，点$P_2$ 和点$P_3\in {\Omega }^{+}$ ，直线DE将K分为两个子单元$K_h^{-}=\Delta P_{1}DE$ 和$K_h^{+}=K\backslash K_h^{-}$ 。定义局部线性基函数${\varphi }_i,i=1,2,3$ 为：

${\varphi }_i\left(x,y\right)=\{\begin{array}{c}{\varphi }_i^{+}\left(x,y\right)=a_i^{+}x+b_i^{+}y+c_i^{+},\left(x,y\right)\in K_h^{+},\\ {\varphi }_i^{-}\left(x,y\right)=a_i^{-}x+b_i^{-}y+c_i^{-},\left(x,y\right)\in K_h^{-},\end{array}$

满足：

$•$ ${\varphi }_i\left(P_j\right)={\delta }_{ij},i,j=1,2,3$ 。

$•$ ${\varphi }_i\left(x,y\right)$ 在DE上连续，即${\varphi }_i^{+}\left(D\right)={\varphi }_i^{-}\left(D\right)$ ，${\varphi }_i^{+}\left(E\right)={\varphi }_i^{-}\left(E\right)$ 。

$•$ ${\varphi }_i\left(x,y\right)$ 的通量在DE上连续，即${\displaystyle {\int }_{\overline{DE}}\beta \nabla {\varphi }_i^{+}\cdot n_h\text{d}s}={\displaystyle {\int }_{\overline{DE}}\beta \nabla {\varphi }_i^{-}\cdot n_h\text{d}s}$。

其中${\delta }_{ij}=\{\begin{array}{l}1i=j\\ 0i\ne j\end{array}$ ，$n_h$ 是DE的单位法向量，方向从$K^{+}$ 指向$K^{-}$ 。局部浸入有限元空间定义为$S_h\left(K\right)=Span\left\{{\varphi }_i,i=1,2,3\right\}$ 。

区域$\Omega$ 上的全局有限元空间定义为：

$S_h\left(\Omega \right)=\left\{v\in L^2\left(\Omega \right):{v|}_K\in S_h\left(K\right),\forall K\in {\mathcal{T}}_h,v在每个顶点上连续,{v|}_{\partial \Omega }=0\right\}$ .

接下来我们给出问题(1)~(4)的空间半离散格式。记${\epsilon }_h^i$ 和${\epsilon }_h^n$ 分别为内部界面边和内部非界面边的集合。对于每条边$e\in {\epsilon }_h^i$ ，设K₁和K₂是相邻的两个单元，$n_e$ 为e的单位法向量，方向从K₁指向K₂。对于函数v，定义e上的平均和跳跃为

${\left\{v\right\}}_e=\frac{1}{2}\left({\left({v|}_{K_1}\right)|}_e+{\left({v|}_{K_2}\right)|}_e\right),{\left[v\right]}_e={\left({v|}_{K_1}\right)|}_e-{\left({v|}_{K_2}\right)|}_e.$ (7)

令$W_e=\left\{w_h\in {\left(L^2\left(\Omega \right)\right)}^2:{w_h|}_{\Omega \backslash \left(K_1\cup K_2\right)}=0\right\}$ ，定义局部提升算子

$r_e:L^2\left(e\right)\to W_e,$

满足

${\displaystyle {\int }_{\Omega }{\beta }_h{r}_e\left(\phi \right)\cdot w_h\text{d}x}={\displaystyle {\int }_e{\left\{\beta w_h\cdot n_e\right\}}_e\phi \text{d}s},\forall w_h\in W_e.$ (8)

方程(1)~(4)的空间半离散格式为：求$U\left(\cdot ,t\right)\in S_h\left(\Omega \right)$ 满足

$\left(\frac{{\partial }^{2}U}{\partial t^2},V\right)+A_h\left(U,V\right)=\left(f,V\right),\forall V\in S_h\left(\Omega \right);$ (9)

$U\left(\cdot ,0\right)={\tilde{u}}_0,$ (10)

$\frac{\partial U}{\partial t}\left(\cdot ,0\right)={\tilde{u}}_1.$ (11)

其中

$A_h\left(U,V\right)=a_h\left(U,V\right)+s_h\left(U,V\right),$ (12)

$a_h\left(U,V\right)={\displaystyle \sum _{K\in {\mathcal{T}}_h}{\displaystyle {\int }_K{\beta }_h\nabla U\cdot \nabla V\text{d}x}}-{\displaystyle \sum _{e\in {\epsilon }_h^i}{\displaystyle {\int }_e({\left\{\beta \nabla U\cdot n_e\right\}}_e{\left[V\right]}_e}}+{\left\{\beta \nabla V\cdot n_e\right\}}_e{\left[U\right]}_e)\text{d}s,$ (13)

$s_h\left(U,V\right)=4{\displaystyle \sum _{e\in {\epsilon }_h^i}{\displaystyle {\int }_{\Omega }{\beta }_h{r}_e\left({\left[U\right]}_e\right)\cdot r_e\left({\left[V\right]}_e\right)\text{d}x}},$ (14)

$\left(f,V\right)={\displaystyle {\int }_{\Omega }fV\text{d}x}.$ (15)

${\tilde{u}}_0$ 和${\tilde{u}}_1$ 分别表示$u_0\left(x,y\right)$ 和$u_1\left(x,y\right)$ 的插值，${\beta }_h\left(x,y\right)=\{\begin{array}{l}{\beta }^{+}\left(x,y\right)\in {\Omega }_h^{+}\\ {\beta }^{-}\left(x,y\right)\in {\Omega }_h^{-}\end{array}$ 。

最后，我们将上述半离散格式进行时间离散。给定一个正整数N，令$\Delta t=T/N$ ，记$U_n=U\left(\cdot ,t_n\right)$ ，其中$t_n=n\Delta t,n=0,\cdots ,N$ 。定义

$U_{n,\frac{1}{4}}=\frac{1}{4}{U}_{n+1}+\frac{1}{2}{U}_n+\frac{1}{4}{U}_{n-1},{\partial }_t^2{U}_n=\frac{U_{n+1}-2U_n+U_{n-1}}{{\left(\Delta t\right)}^2},$

则可以得到如下全离散格式：求$U_n\in S_h\left(\Omega \right)$ ，对于$n=2,\cdots ,N$ 满足

$\left({\partial }_t^2{U}_n,V\right)+A_h\left(U_{n,\frac{1}{4}},V\right)=\left(f_{n,\frac{1}{4}},V\right),\forall V\in S_h\left(\Omega \right)$ (16)

$U_0={\tilde{u}}_0,$ (17)

$U_1={\tilde{u}}^{*},$ (18)

其中${\tilde{u}}^{*}$ 是$u^{*}$ 的插值，$u^{*}=u_0+\Delta t{u}_1+\frac{\Delta t^2}{2}{u}_{tt}\left(x,y,0\right)$ ，且$u_{tt}\left(x,y,0\right)$ 可由(1)计算得出。

4. 数值算例

例1 考虑计算区域为$\Omega =\left(-1,1\right)\times \left(-1,1\right)$ ，界面为$\Gamma =\left\{\left(x,y\right)|x^2+y^2-r_0^2=0\right\}$ ，真解为

$u\left(x,y,t\right)=\{\begin{array}{l}\frac{1}{{\beta }^{-}}{r}^{\alpha }\cos \left(2t\right),\left(x,y,t\right)\in {\Omega }^{-}\times \left[0,1\right],\\ \left(\frac{1}{{\beta }^{+}}{r}^{\alpha }+\left(\frac{1}{{\beta }^{-}}-\frac{1}{{\beta }^{+}}\right)r_0^{\alpha }\right)\cos \left(2t\right),\left(x,y,t\right)\in {\Omega }^{+}\times \left[0,1\right].\end{array}$

其中$r_0=\pi /6.28$ ，$\alpha =3$ ，$r=\sqrt{x^2+y^2}$ ，${\beta }^{-}=1$ ，${\beta }^{+}=10$ 。

我们将计算区域划分为$N\times N$ 个矩形，再沿着同一方向把每个矩形由对角线分为两个三角形。考虑时间$T=1$ ，在时间方向选取步长$\Delta t=T/N$ 。表1给出了$L^2$ 误差以及误差的收敛阶，可以看出数值解的$L^2$ 误差为2阶收敛。

Table 1. Error and convergence order of example 1

表1. 例1的误差及收敛阶

N	L^2误差	收敛阶
4	0.78991E−00
8	0.22660E−00	1.80157
16	0.61090E−01	1.89111
32	0.15483E−01	1.98028
64	0.39185E−02	1.98229
128	0.98193E−03	1.99661
256	0.24564E−03	1.99911

例2 考虑某种声波在传播过程中遇到两种不同均匀介质的传播情况，假设传播区域为

$\Omega =\left(-2.1,2.1\right)\times \left(-2.1,2.1\right)$ ，椭圆形界面$\Gamma =\left\{\left(x,y\right)|\frac{{\left(x-x_0\right)}^2}{r_x^2}+\frac{{\left(y-y_0\right)}^2}{r_y^2}-1=0\right\}$ 将区域$\Omega$ 划分为两个子域：

${\Omega }^{-}=\left\{\left(x,y\right)|\frac{{\left(x-x_0\right)}^2}{r_x^2}+\frac{{\left(y-y_0\right)}^2}{r_y^2}<1\right\},{\Omega }^{+}=\Omega \backslash \overline{{\Omega }^{-}}.$

其中$\left(x_0,y_0\right)=\left(1.15,0\right)$ ，$r_x=\pi /4.84$ ，$r_y=\pi /1.97$ ，设${\beta }^{-}=1$ ，${\beta }^{+}=4$ ，$f=0$ ，$u_1=0$ ，$T=0.6$ ，并且初始值$u_0\left(x,y\right)$ 设为：

$u_0\left(x,y\right)=\{\begin{array}{l}0,\left(x,y\right)\notin S,\\ 600\cdot \exp \left(\frac{-1}{1-\sqrt{{\left(x-0.3\right)}^2+y^2}-\sqrt{{\left(x+0.3\right)}^2+y^2}}\right),\left(x,y\right)\in S.\end{array}$

其中$S=\left\{\left(x,y\right)|\frac{x^2}{0.25}+\frac{y^2}{0.16}-1<0\right\}$ 。

使用$N=256$ 的网格，设置$\Delta t=T/N$ ，选取$T=0,T=0.3,T=0.6$ 三个时刻，计算出了图3所示波在传播过程中冲击界面前后的数值结果。

(a) $t=0$ (b) $t=0.3$

(c) $t=0.6$

Figure 3. The solution $u_h$ at $t=0,t=0.3,t=0.6$

图3. 在$t=0,t=0.3,t=0.6$ 时的解$u_h$

例3 考虑某种在两种不同均匀介质中平行向右传播的波，假设传播区域为$\Omega =\left(0,5\right)\times \left(0,3\right)$ ，圆形界面$\Gamma =\left\{\left(x,y\right)|{\left(x-3\right)}^2+{\left(y-1.5\right)}^2=1\right\}$ 将区域$\Omega$ 划分为两个子域：

${\Omega }^{-}=\left\{\left(x,y\right)|{\left(x-3\right)}^2+{\left(y-1.5\right)}^2<1\right\},{\Omega }^{+}=\Omega \backslash \overline{{\Omega }^{-}}.$

设$f=0,T=3$ 并且初始值$u_0\left(x,y\right)$ ，$u_1\left(x,y\right)$ 分别设为：

$\begin{array}{l}{u}_0\left(x,y\right)=\exp \left(-\frac{3}{2}{x}^2\right),\left(x,y\right)\in \Omega \\ u_1\left(x,y\right)=3x\cdot \exp \left(-\frac{3}{2}{x}^2\right),\left(x,y\right)\in \Omega \end{array}$

另外我们将(2)中的零边界条件换成非齐次边界条件$u\left(x,y,t\right)=g\left(x,y,t\right)$ ，其中

$g\left(x,y,t\right)=\exp \left(-\frac{3}{2}{\left(x-t\right)}^2\right),\left(x,y,t\right)\in \Omega \times \left[0,3\right]$

再使用传统的处理方法，令$N=256$ ，设置$\Delta t=T/N$ ，选取$T=0,T=2,T=3$ 三个时刻，分别计算出${\beta }^{-}=0.2,{\beta }^{+}=1$ 和${\beta }^{-}=5,{\beta }^{+}=1$ 两种系数条件下波在传播过程中冲击界面前后的数值结果，模拟出了波在这两种均匀介质中的传播，如图4、图5所示。

(a) $t=0$

(b) $t=2$

(c) $t=3$

Figure 4. The solution $u_h$ at $t=0,t=2,t=3$ when ${\beta }^{-}=0.2,{\beta }^{+}=1$

图4. ${\beta }^{-}=0.2,{\beta }^{+}=1$ 时在$t=0,t=2,t=3$ 时的解$u_h$

(a) $t=0$

(b) $t=2$

(c) $t=3$

Figure 5. The solution $u_h$ at $t=0,t=2,t=3$ when ${\beta }^{-}=5,{\beta }^{+}=1$

图5. ${\beta }^{-}=5,{\beta }^{+}=1$ 时在$t=0,t=2,t=3$ 时的解$u_h$

5. 结论

本文提出了一种基于界面非拟合网格的浸入有限元方法。该方法能够数值求解具有间断系数波动方程，与经典的部分罚浸入有限元方法相比，该方法具有无需人为调整罚参数的优点。本文给出了全离散格式，并且通过一系列数值试验验证了该方法的有效性。然而目前只是针对二维问题提出了本文的方法，接下来我们希望推广该方法来处理三维的问题，并且理论分析也是接下来的研究工作之一。