求解非齐次跳跃界面问题的新型浸入有限元法

A New Immersed Finite Element Method for Interface Problems with Non-Homogeneous Jump Conditions

摘要: 针对非齐次跳跃条件的椭圆型界面问题，提出了新型浸入有限元方法。该方法基于非拟合网格，将分片线性多项式延拓到整个平面，用来逼近有限元函数。通过延拓点以及延拓点之间界面上的函数积分平均来求解分片线性多项式。数值算例验证了该方法的有效性。

Abstract: A new immersed finite element method is proposed to solve the elliptic interface problem with non-homogeneous jump conditions. This method is based on unfitted mesh, and the piecewise linear polynomial is extended to the whole plane to approximate the finite element function. The piecewise linear polynomial functions are solved by the extension points and the average of the integral of function on the interface between the extension points. Numerical examples verify the effectiveness of the method.

1. 引言

在很多重要科学应用中都存在界面问题，例如多介质中的电磁场问题、复合材料中的热传导问题、多相流问题等[1] [2] [3]。通常，界面问题的解在界面上需要满足非齐次跳跃条件，例如求解静电场问题中界面上存在面电荷、多相流问题中界面上存在表面张力等。除了某些特殊情况外，解析求解带有非齐次跳跃的界面问题一般比较困难[4]。因此设计快速、有效的数值方法对解决实际问题有重要的意义。

根据离散单元和界面的关系，界面问题的有限元求解分为拟合网格法和非拟合网格法[5]。拟合网格法：网格必须与界面互相匹配，即，界面附近的网格点必须落在界面上。用拟合网格法求解界面问题的好处在于直接应用传统的有限元方法，便可以得到最优误差估计[6] [7]。在实际计算中，拟合网格法有着许多缺点。例如：对于形状复杂的界面生成拟合网格往往比较困难；对于移动界面问题，在每一时间层都需要重新生成拟合网格法，而且需要构造插值算子来插值定义在不同拟合网格上的数值解，增加了计算量.为了克服上述缺点，另外一类方法是非拟合网格法[8] [9] [10] [11]。非拟合网格法由于网格生成简单，而且对移动界面问题优点显著，受到计算数学界和工程界的青睐和广泛关注。

在众多非匹配网格方法求解界面问题的方法中，浸入界面有限元方法在计算方面有着自由度少，易于编程等明显的优点[12] [13] [14] [15]。该方法通过修改界面单元的基函数来获得最优收敛性。然而，这种修改会导致基函数在界面穿过的边上不连续，进而产生相容误差。针对该问题，Lin等[12]提出了求解齐次跳跃问题的部分罚浸入有限元方法。该方法利用间断有限元思想，在界面穿过的边上添加额外的惩罚项来惩罚基函数的不连续性。随后，Ji等[16]通过构造满足非齐次跳跃条件的修正函数，提出了求解具有非齐次跳跃条件界面问题的部分罚浸入有限元方法。该方法局限性在于构造修正函数时，用到了界面上函数的点值，但在正则性较低的情况下，函数的点值无意义。He等[17]在求解非齐次跳跃问题的浸入有限元函数时利用了积分平均的方式。该方法的局限性在于积分的区间会趋向于0，导致积分平均趋于无穷。

针对二阶非齐次跳跃条件下的椭圆界面问题。在求解齐次跳跃问题的浸入有限元方法基础上，本文对每个界面元素添加分片线性多项式，即修正函数，使其能满足非齐次跳跃条件，从而达到丰富IFE空间的目的，最后利用部分罚的浸入有限元方法求解非齐次跳跃界面问题。求解修正函数过程中，针对界面上点的函数值无意义以及积分平均的积分区间趋于0的情况。本文扩大积分的区间，保证积分区间始终大于常数C，继而利用积分平均来求解修正函数。最后通过数值算例来说明该方法的正确性。

2. 椭圆界面问题

考虑如下椭圆界面问题：

$-\nabla \cdot \left(\beta \left(x\right)\nabla u\left(x\right)\right)=f\left(x\right),x\in \Omega \backslash \Gamma ,$ (1)

$u\left(x\right)=0,x\in \partial \Omega ,$ (2)

其中解满足如下界面跳跃条件：

${\left[u\right]}_{\Gamma }=u^{+}-u^{-}=w,$ (3)

${\left[\beta \frac{\partial u}{\partial n}\right]}_{\Gamma }={\beta }^{+}\nabla u^{+}\cdot n-{\beta }^{-}\nabla u^{-}\cdot n=Q,$ (4)

其中：$\frac{\partial u}{\partial n}=\nabla u\cdot n$ 为u的法向导数，$n$ 为$\Gamma$ 单位法向量并指向${\Omega }^{+}$ ，${[\cdot ]}_{\Gamma }$ 表示界面$\Gamma$ 上的跳跃。如图1，区域$\Omega$ 是${ℝ}^2$ 上的一个矩形域，$\Gamma$ 为区域$\Omega$ 中的一条封闭光滑曲线，将区域$\Omega$ 分为${\Omega }^{+}$ 和${\Omega }^{-}$ 两个子域，$\Omega ={\Omega }^{+}\cup {\Omega }^{-}\cup \Gamma$ 。假设子域${\Omega }^{-}$ 严格包含在区域$\Omega$ 内部，其中$u^{\pm }\left(x\right)={u\left(x\right)|}_{{\Omega }^{\pm }}$ ，系数函数$\beta \left(x\right)$ 在区域$\Omega$ 上为分片函数，即

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

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

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

3. 浸入有限元方法

如图2所示，考虑$\Omega$ 上的一致三角剖分${\mathcal{T}}_h$ ，其网格尺寸$h=\sqrt{{\left(\Delta x\right)}^2+{\left(\Delta y\right)}^2}$ ，其边的集合记为${\epsilon }_h$ 。当h充分小时，可以假设界面$\Gamma$ 和剖分单元的边至多有两个交点，且这两个交点落在不同的边上。这些与界面相交的边的集合记为${\epsilon }_h^{\mathrm{int}}:=\left\{e\in {\epsilon }_h:e\cap \Gamma \ne \varnothing \right\}$ 。界面单元的集合记为${\mathcal{T}}_h^{\mathrm{int}}:=\left\{T\in {\mathcal{T}}_h:T\cap \Gamma \ne \varnothing \right\}$ ，非界面单元的集合为${\mathcal{T}}_h^{non}:={\mathcal{T}}_h\backslash {\mathcal{T}}_h^{\mathrm{int}}$ 。

Figure 2. Mesh generation

图2. 网格剖分

3.1. 非齐次跳跃条件下的浸入有限元函数

定义界面$\Gamma$ 的水平集函数$\phi \left(x\right)$ 为：

$\phi \left(x\right)\{\begin{array}{l}<0,x\in {\Omega }^{-},\\ =0,x\in \Gamma ,\\ >0,x\in {\Omega }^{+}.\end{array}$ (6)

设$I_h$ 为基于节点的线性插值算子，定义近似界面为${\Gamma }_h=\left\{x\in \Omega |I_h\phi \left(x\right)=0\right\}$ 。设$T=\Delta A_1{A}_2{A}_3$ ，考虑$T\in {\mathcal{T}}_h^{\mathrm{int}}$ ，如图3所示，${\Gamma }_h$ 与$A_1{A}_2,A_1{A}_3$ 相交于$D,E$ 两点，线段DE将T分成了$T_h^{+}$ 和$T_h^{-}$ ，$n_h$ 和$t_h$ 是线段DE的单位法向量和单位切向量。$D,E$ 的中点为$B_0=\left(D+E\right)/2$ 。$B_1$ 点为$B_0-h\cdot t_h$ ，$B_2$ 点为$B_0+h\cdot t_h$ 。定义映射$p:B_1{B}_2\to \Gamma$ ，形式为$p\left(x\right)=x+\alpha n_h$ ，其中$\alpha \in ℝ$ 。

Figure 3. Interface element

图3. 界面单元

设${\mathbb{P}}_1$ 为${ℝ}^2$ 上次数不超过一次的多项式的集合，定义界面单元T上的分片线性多项式形式为：

$v_{h,T}\left(x\right)=\{\begin{array}{l}{v}^{+}\left(x\right)\in {\mathbb{P}}_1,x\in T_h^{+},\\ v^{-}\left(x\right)\in {\mathbb{P}}_1,x\in T_h^{-},\end{array}$ (7)

满足：

$\begin{array}{l}{v}_{h,T}\left(A_1\right)=V_1,v_{h,T}\left(A_2\right)=V_2,v_{h,T}\left(A_3\right)=V_3,\\ v^{+}\left(B_1\right)-v^{-}\left(B_1\right)=w\left(p\left(B_1\right)\right),v^{+}\left(B_2\right)-v^{-}\left(B_2\right)=w\left(p\left(B_2\right)\right),\\ {\left[\beta \frac{\partial v_{h,T}}{\partial n_h}\right]}_{{\Gamma }_h\cap T}=\frac{{\displaystyle {\int }_{p\overline{\left(B_1{B}_2\right)}}Q\text{d}s}}{\left|p\left(\overline{B_1{B}_2}\right)\right|},\end{array}$ (8)

其中$V_1,V_2,V_3$ 为节点变量。

显然$v_{h,T}$ 有如下分解：$v_{h,T}=u_{h,T}^{\mathrm{hom}}+u_{h,T}^J$ ，其中

$u_{h,T}^{\mathrm{hom}}\left(x\right)=\{\begin{array}{l}{u}^{\mathrm{hom},+}\left(x\right)\in {\mathbb{P}}_1,x\in T_h^{+},\\ u^{\mathrm{hom},-}\left(x\right)\in {\mathbb{P}}_1,x\in T_h^{-},\end{array}$ (9)

满足：

$\begin{array}{l}{u}_{h,T}^{\mathrm{hom}}\left(A_1\right)=V_1\text{, }{u}_{h,T}^{\mathrm{hom}}\left(A_2\right)=V_2,u_{h,T}^{\mathrm{hom}}\left(A_3\right)=V_3,\\ u^{\mathrm{hom},+}\left(B_1\right)-u^{\mathrm{hom},-}\left(B_1\right)=0,u^{\mathrm{hom},+}\left(B_2\right)-u^{\mathrm{hom},-}\left(B_2\right)=0,\\ {\left[\beta \frac{\partial u_{h,T}^{\mathrm{hom}}}{\partial n_h}\right]}_{{\Gamma }_h\cap T}=0,\end{array}$ (10)

$u_{h,T}^J\left(x\right)=\{\begin{array}{l}{u}^{J,+}\left(x\right)\in {\mathbb{P}}_1,x\in T_h^{+},\\ u^{J-}\left(x\right)\in {\mathbb{P}}_1,x\in T_h^{-},\end{array}$ (11)

满足：

$\begin{array}{l}{u}_{h,T}^J\left(A_1\right)=0,u_{h,T}^J\left(A_2\right)=0,u_{h,T}^J\left(A_3\right)=0,\\ u^{J,+}\left(B_1\right)-u^{J,-}\left(B_1\right)=w\left(p\left(B_1\right)\right),u^{J,+}\left(B_2\right)-u^{J,-}\left(B_2\right)=w\left(p\left(B_2\right)\right),\\ {\left[\beta \frac{\partial u_{h,T}^J}{\partial n_h}\right]}_{{\Gamma }_h\cap T}=\frac{{\displaystyle {\int }_{p\left(\overline{B_1{B}_2}\right)}Q\text{d}s}}{\left|p\left(\overline{B_1{B}_2}\right)\right|}.\end{array}$ (12)

定义局部浸入有限元空间：

$S_h\left(T\right)=\left\{u_{h,T}^{\mathrm{hom}}|u_{h,T}^{\mathrm{hom}}满足\left(9\right)-\left(10\right),其中V_1,V_2,V_3为任意实数\right\}.$

网格为一致三角剖分，函数$u_{h,T}^{\mathrm{hom}}\in S_h\left(T\right)$ 由$u_{h,T}^{\mathrm{hom}}\left(A_i\right)\left(i=1,2,3\right)$ 唯一确定[18]，其中$A_i\left(i=1,2,3\right)$ 为T的顶点，函数$u_{h,T}^J$ 可以根据(11)~(12)唯一确定[17]。

注：He等[17]采用积分平均$\frac{{\displaystyle {\int }_{\Gamma \cap T}Q\text{d}s}}{\left|\overline{DE}\right|}$，其中$Q\in L^2$ ，若$\left|\overline{DE}\right|\to 0$ ，会使$\frac{{\displaystyle {\int }_{\Gamma \cap T}Q\text{d}s}}{\left|\overline{DE}\right|}\to \infty$。Ji等[16]使用修正函数中$Q\in L^2$ ，$x^{\ast }\in \Gamma$ 时，$Q\left(x^{\ast }\right)$ 可能会没有意义。

3.2. 部分罚的浸入有限元方法

设e是单元T₁和T₂的公共边，$n_e$ 是e单位法向量并指向T₂。定义e上的平均值和跳跃：

$\left\{\nabla u\right\}:=\frac{1}{2}\left({\nabla u|}_{T_1}+{\nabla u|}_{T_2}\right),〚u〛n_e:={u|}_{T_1}{n}_e-{u|}_{T_2}{n}_e.$ (13)

定义映射$p_h:{\Gamma }_h\to \Gamma$ ，形式为$p_h\left(x\right)=x+{\alpha }_h{n}_h$ ，其中${\alpha }_h\in ℝ$ 。定义在界面$\Gamma$ 上的函数Q通过$Q_h=Q\circ p_h$ 转换到${\Gamma }_h$ 。

对于非界面单元$T\in {\mathcal{T}}_h^{non}$ ，我们采用传统线性有限元空间，也记为$S_h\left(T\right)$ 。那么全局浸入有限元空间可以定义为：

$S_h\left(\Omega \right)=\left\{v_h|{v_h|}_T\in S_h\left(T\right),\forall T\in {\mathcal{T}}_h;{v_h|}_{\partial \Omega }=0;v_h在节点连续\right\}.$ (14)

考虑的有限元方法是求$u_h=u_h^{\mathrm{hom}}+u_h^J,u_h^{\mathrm{hom}}\in S_h\left(\Omega \right)$ ，使得

$a_h\left(u_h^{\mathrm{hom}},v_h\right)={\displaystyle {\int }_{\Omega }f}{v}_h\text{d}x-{\displaystyle {\int }_{{\Gamma }_h}{Q}_h}{v}_h\text{d}s-a_h\left(u_h^J,v_h\right),\forall v_h\in S_h\left(\Omega \right),$ (15)

其中

$a_h\left(w,v_h\right)={\displaystyle \sum _{T\in T_h}{\displaystyle \sum _{i=\pm }{\displaystyle {\int }_{T\cap {\Omega }^i}\beta \nabla w\cdot \nabla v_h\text{d}x}}}-{\displaystyle \sum _{e\in {\epsilon }_h^{\mathrm{int}}}{\displaystyle {\int }_e\left(\left\{\beta \nabla w\right\}\cdot n_e{\left[v_h\right]}_e+\left\{\beta \nabla v_h\right\}\cdot n_e{\left[w\right]}_e-\frac{\eta }{h}{\left[w\right]}_e{\left[v_h\right]}_e\right)\text{d}s}}.$ (16)

本文部分罚浸入有限元方法与标准的线性协调有限元方法相比，不同点在于(15)的右端项增加了

${\displaystyle {\int }_{{\Gamma }_h}{Q}_h}{v}_h\text{d}s$ 和$a_h\left(u_h^J,v_h\right)$ ，但左端项$a_h\left(u_h^{\mathrm{hom}},v_h\right)$ 和标准的线性协调有限元方法的左端项类似且都为刚度矩

阵，根据标准的线性协调有限元方法是有解，可知(15)同样有解。当本文的式(3)~(4)中$w=0,Q=0$ 以及${\beta }^{+}={\beta }^{-}$ 时，此时(15)中的$Q_h=0$ ，$u_h^J=0$ ，(16)中的第二项也变为0，$S_h\left(\Omega \right)$ 将变为${\mathbb{P}}_1$ 协调元，则本文部分罚浸入有限元方法将转化为标准的线性协调有限元方法。

4. 数值算例

以下讨论几个数值算例来说明该方法的有效性。考虑矩形区域$\Omega$ 为$\left[-1,1\right]\times \left[-1,1\right]$ ，已知解析解$u^{+}\left(x_1,x_2\right),u^{-}\left(x_1,x_2\right)$ ，系数${\beta }^{+}\left(x_1,x_2\right),{\beta }^{-}\left(x_1,x_2\right)$ 以及水平集函数$\phi \left(x_1,x_2\right)$ 。将区域$\Omega$ 划分为$2N^2$ 个直角

三角形，网格尺寸为h。在下述所有算例中，理论上需要将$\eta$ 设置充分大，本文发现将$\eta$ 设置为0时，仍然可以得到正确解。

定义如下范数：

$\begin{array}{l}{\Vert u-u_h\Vert }_{L^{\infty }}=\underset{x\in N_h}{\max }\left|u\left(x\right)-u_h\left(x\right)\right|,N_h为\Omega 内所有的顶点，\\ {\Vert u-u_h\Vert }_{L^2}={\displaystyle \sum _{T\in T_h}{\left({\displaystyle {\int }_T{\left|u\left(x\right)-u_h\left(x\right)\right|}^2\text{d}x}\right)}^{\frac{1}{2}}},{\Vert u-u_h\Vert }_{H^1}={\displaystyle \sum _{T\in T_h}{\left({\displaystyle \sum _{\left|\alpha \right|<1}{\displaystyle {\int }_T{\left|D^{\left|\alpha \right|}\left(u\left(x\right)-u_h\left(x\right)\right)\right|}^2\text{d}x}}\right)}^{\frac{1}{2}}}.\end{array}$

例1 解析解$u^{+}\left(x_1,x_2\right),u^{-}\left(x_1,x_2\right)$ ,系数${\beta }^{+}\left(x_1,x_2\right),{\beta }^{-}\left(x_1,x_2\right)$ 以及水平集函数$\phi \left(x_1,x_2\right)$ 如下

$\begin{array}{l}{u}^{+}\left(x_1,x_2\right)=\ln \left(x_1^2+x_2^2\right),u^{-}\left(x_1,x_2\right)=\sin \left(x_1+x_2\right),\\ {\beta }^{+}\left(x_1,x_2\right)=\sin \left(x_1+x_2\right)+2,{\beta }^{+}\left(x_1,x_2\right)=\cos \left(x_1+x_2\right)+2,\\ \phi \left(x_1,x_2\right)=x_1^2+x_2^2-{0.5}^2.\end{array}$

Table 1. Error and convergence order of example 1

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

N	${\Vert u-u_h\Vert }_{L^{\infty }}$	阶数	${\Vert u-u_h\Vert }_{L^2}$	阶数	${\Vert u-u_h\Vert }_{H^1}$	阶数
64	1.1142E−03	--	1.0285E−03	--	9.8450E−02	--
128	2.9693E−04	1.91	2.6053E−04	1.98	4.9469E−02	0.99
256	7.7437E−05	1.94	6.6200E−05	1.98	2.4791E−02	1.00
512	2.0057E−05	1.95	1.6688E−05	1.99	1.2412E−02	1.00
1024	5.0503E−06	1.99	4.2222E−06	1.98	6.2099E−03	1.00

Figure 4. Approximate solution

图4. 近似解

Figure 5. Error distribution

图5. 误差分布

由表1可知，该算例通过此方法得到的数值解的$L^2$ 误差为2阶，$H^1$ 误差为1阶，从而验证了方法的有效性。图4和图5为此算例$N=64$ 时该方法的近似解和误差分布。

例2 解析解$u^{+}\left(x_1,x_2\right),u^{-}\left(x_1,x_2\right)$ ，系数${\beta }^{+}\left(x_1,x_2\right),{\beta }^{-}\left(x_1,x_2\right)$ 以及水平集函数$\phi \left(x_1,x_2\right)$ 如下

$\begin{array}{l}{u}^{+}\left(x_1,x_2\right)=4-x_1^2-x_2^2,u^{-}\left(x_1,x_2\right)=x_1^2+x_2^2,\\ {\beta }^{+}\left(x_1,x_2\right)=x_1{x}_2+2,{\beta }^{+}\left(x_1,x_2\right)=x_1^2-x_2^2+3,\\ \phi \left(x_1,x_2\right)=x_1^2-x_2-1.\end{array}$

Table 2. Error and convergence order of example 2

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

N	${\Vert u-u_h\Vert }_{L^{\infty }}$	阶数	${\Vert u-u_h\Vert }_{L^2}$	阶数	${\Vert u-u_h\Vert }_{H^1}$	阶数
64	1.3103E-03	--	1.0268E-03	--	5.3567E-02	--
128	3.7122E-04	1.82	2.5350E-04	2.02	2.6108E-02	1.04
256	1.2264E-04	1.60	6.3428E-05	2.00	1.2899E-02	1.02
512	3.0989E-05	1.98	1.5833E-05	2.00	6.4130E-03	1.01
1024	8.7151E-06	1.83	3.9418E-06	2.01	3.1978E-03	1.00

Figure 6. Approximate solution

图6. 近似解

Figure 7. Error distribution

图7. 误差分布

由表2可知，该算例通过此方法得到的数值解的$L^2$ 误差为2阶，$H^1$ 误差为1阶，从而验证了方法的有效性。图6和图7为此算例$N=64$ 时该方法的近似解和误差分布。

5. 结论

针对具有非齐次跳跃条件的二阶椭圆型界面问题，本文提出了一种新型浸入有限元方法。该方法在齐次跳跃浸入有限元方法的基础上，构造满足非齐次跳跃条件的修正函数。本文通过求解分片线性多项式来得到所需的修正函数，在求解分片线性多项式过程中使用了界面与单元的交点的延拓点以及延拓点之间界面上的函数积分平均，最后使用部分罚浸入有限元方法求解具有非齐次跳跃条件的界面问题。数值算例的结果显示$L^2$ 误差为2阶收敛，$H^1$ 的误差为1阶收敛，验证了该方法的有效性。

NOTES

^*通讯作者。