无网格Galerkin法求解二维Helmholtz方程

An Element-Free Galerkin Method for Solving Two-Dimensional Helmholtz Equation

摘要: 无网格Galerkin法本来是对固体力学问题进行数值模拟的一个重要方法，在本文我们用此方法来对Helmholtz方程进行求解，并对求解所得的数值模拟结果进行了验证。在求解Helmholtz方程的过程中，首先用移动最小二乘近似来构造形函数，然后边值条件由拉格朗日乘子法引入，最后对平面规则域上的一般Helmholtz方程进行了数值模拟，所得到的结果趋向于精确解，且随着节点的增加，其精确度越来越高，验证了该方法具有良好的收敛性。

Abstract: The element-free Galerkin method is originally an important method for numerical simulation of solid mechanics problems. In this paper, we use this method to solve the Helmholtz equation and verify the numerical simulation results. In the process of solving Helmholtz equation, the shape function is constructed by moving least squares approximation, then the boundary condition is introduced by Lagrange multiplier method. Finally, the general Helmholtz equation on the plane rule domain is numerically simulated. The obtained results tend to be exact solutions. Furthermore, high precision and good convergence could be contained as the nodes increased.

1. 引言

无网格法是一种通过节点信息来求解偏微分方程问题的数值方法，相比于有限元法和有限差分法等传统的数值方法，无网格法可减少或完全消除对网格的依赖，还保证了计算的精度，降低了计算的难度。Belytschko等人在1994年提出了无网格Galerkin法 [1] ，并将其运用到对固体力学问题进行数值模拟的方面上。虽然全局弱式的无网格Galerkin法具有精度高、稳定性好等优点，但仍需要背景网格或背景单元进行积分。

定义1：无网格法是在建立的整个问题域的系统代数方程时，不需要利用预定义的网格信息进行域离散的方法 [2] 。

在Helmholtz问题的无网格法求解中，李美香等使用了点插值法 [3] 、姜欣荣等使用基本解法(MFS)和边界节点法(BKM) [4] 以及杨子乐等使用无网格介点法(MIP) [5] [6] 来求解此类问题。对于无网格Galerkin法，Mehdi，Dehghan等人用这种方法来较好的模拟了肿瘤细胞的扩散现象 [7] 。

2. 移动最小二乘近似(MLS)

假设 $\Omega \subset R^d$ 是一个非空有界集，设 $X=\left\{x_1,x_2,\cdots ,x_N\right\}$ 为在求解区域 $\Omega$ 的N个点的集合， $u:\Omega \to R$ 是一个连续函数。

那么对于 $\forall x\in \Omega$ 的移动最小二乘近似可以被写为：

$u\left(x\right)\approx u^h\left(x,x_i\right)={\displaystyle \underset{i\in I\left(x\right)}{\sum }{\varphi }_i\left(x\right)u\left(x_i\right)},i=1,2,\cdots ,I\left(x\right)$ (1)

(1)式中 ${\varphi }_i\left(x\right)$ 是形函数， $u\left(x_i\right)$ 是在节点 $x_i$ 上的节点参数， $I\left(x\right)$ 表示在支持域 $x$ 中的节点指标集，定义如下所示：

$I\left(x\right)=\left\{i\in \left\{1,2,\cdots ,N\right\}:{\Vert x-x_i\Vert }_2\le \delta \right\}$

$\delta$ 为支持域的半径。

(1)式可由近似函数 $u^h\left(x\right)$ 在这些节点 $x=x_i$ 处误差的加权平方和推导可得，如下：

$u^h\left(x,x_i\right)={\displaystyle \underset{i\in I\left(x\right)}{\sum }{p}_i\left(x_i\right)a_i\left(x\right)}=p^{\text{T}}\left(x_i\right)a\left(x\right),p\left(x\right)\in P_m\left(R^d\right)$ (2)

$J\left(a\right)={\displaystyle \underset{i\in I\left(x\right)}{\sum }{\left(u^h\left(x,x_i\right)-u\left(x_i\right)\right)}^{\text{2}}w\left(x,x_i\right)}$ (3)

$p\left(x\right)$ 为基函数， $P_m\left(R^d\right)$ 中m指的是基函数的最高次幂，d指的是空间的维数。则函函数的个数为 $Q=\left(\begin{array}{c}m+d\\ d\end{array}\right)$ 。

其中(2)式为u的局部近似，(3)式为精确解与近似解的误差加权平方和。

通过(3)式求关于 $a\left(x\right)$ 的驻点值，可得：

$\frac{\partial J}{\partial a}=0$

并对该式进行运算和整理可得形函数为：

${\Phi }^{\text{T}}\left(x\right):=\left[{\varphi }_1\left(x\right),{\varphi }_2\left(x\right),\cdots ,{\varphi }_{\left|I\left(x\right)\right|}\left(x\right)\right]=p^{\text{T}}\left(x\right)A^{-1}\left(x\right)B\left(x\right)$ (4)

其中 $A\left(x\right)$ 和 $B\left(x\right)$ 为：

$A\left(x\right)={\displaystyle \underset{i\in I\left(x\right)}{\sum }w\left(x,x_i\right)p\left(x_i\right)p^{\text{T}}\left(x_i\right)}$ (5)

$B\left(x\right)=\left[w\left(x,x_1\right)p\left(x_1\right)w\left(x,x_2\right)p\left(x_2\right)\cdots w\left(x,x_{I\left(x\right)}\right)p\left(x_{I\left(x\right)}\right)\right]$ (6)

在(3) (4) (5)式中的 $w\left(x,x_i\right)$ 为权函数，在本文中我们选择用Gauss权函数，如下所示：

$w\left(x,x_i\right)=\{\begin{array}{l}\frac{\exp \left(-{\epsilon }^2{r}^2\right)-\exp \left(-{\epsilon }^2\right)}{1-\exp \left(-{\epsilon }^2\right)},0\le r\le 1\\ 0,r>1\end{array}$ (7)

$r:=\Vert x-x_i\Vert /\delta$ ，其中 $\delta$ 为节点 $x_i$ 的支撑域半径，且 $\epsilon$ 是一个正参数。

3. Helmholtz方程的离散格式

考虑二维的Helmholtz方程 [5] ：

$\{\begin{array}{l}{\nabla }^{2}u+k^{2}u=f,x在\Omega 内\\ u=\bar{u},x在{\Gamma }_u上\\ \frac{\partial u}{\partial n}=q,x在{\Gamma }_t上\end{array}$ (8)

其中k为波数，u为解函数，q为解函数的外法向导数， $n_x$ 和 $n_y$ 分别为边界 ${\Gamma }_u$ 上的外法线方向余弦， ${\Gamma }_t$ 和 ${\Gamma }_u$ 分别为Neumann边界条件和Dirichlet边界条件，两种边界条件构成了求解域围成的边界 $\partial \Omega$ 。那么接下来运用能量法之一的最小势能原理引出平衡方程的另一种表达方式，然后将求解问题转化为变分问题，并引入拉格朗日乘子法来对原有方程施加边界条件，最后可通过所得的离散方程求得近似数值解。

首先通过最小势能原理 [8] 可得：

$E=\frac{\text{1}}{\text{2}}{\displaystyle {\int }_{\Omega }{u}_x^2+u_y^2+k^2{u}^2\text{d}\Omega }-{\displaystyle {\int }_{\Omega }{u}^{\text{T}}f\text{d}\Omega }-{\displaystyle {\int }_{{\Gamma }_t}{u}^{\text{T}}q\text{d}\Gamma }$ (9)

对(9)式进行变分可表示为：

$\delta E=\delta \left(\frac{\text{1}}{\text{2}}{\displaystyle {\int }_{\Omega }{u}_x^2+u_y^2+k^2{u}^2\text{d}\Omega }-{\displaystyle {\int }_{\Omega }{u}^{\text{T}}f\text{d}\Omega }-{\displaystyle {\int }_{{\Gamma }_t}{u}^{\text{T}}q\text{d}\Gamma }\right)$ (10)

根据变分的定义以及链式法则整理可得：

$\delta E={\displaystyle {\int }_{\Omega }{\left(\delta u\right)}_x{u}_x+{\left(\delta u\right)}_y{u}_y+\left(\delta u\right)k^{2}u\text{d}\Omega }-{\displaystyle {\int }_{\Omega }\delta u^{\text{T}}f\text{d}\Omega }-{\displaystyle {\int }_{{\Gamma }_t}\delta u^{\text{T}}q\text{d}\Gamma }$ (11)

因为总势能的最小要求为 $\delta E=0$ ，故可得到：

${\displaystyle {\int }_{\Omega }{\left(\delta u\right)}_x{u}_x+{\left(\delta u\right)}_y{u}_y+\left(\delta u\right)k^{2}u\text{d}\Omega }-{\displaystyle {\int }_{\Omega }\delta u^{\text{T}}f\text{d}\Omega }-{\displaystyle {\int }_{{\Gamma }_t}\delta u^{\text{T}}q\text{d}\Gamma }=0$ (12)

令 $\delta u=\Phi \left(x\right)$ ，将(1)式代入(12)式可得：

${\displaystyle \underset{i\in I\left(x\right)}{\sum }{\displaystyle \underset{j=1}{\overset{N}{\sum }}\left(-{\displaystyle {\int }_{\Omega }\nabla {\varphi }_i\cdot \nabla {\varphi }_j\text{d}\Omega }+{\displaystyle {\int }_{\Omega }{k}^2{\varphi }_i\cdot {\varphi }_j}\text{d}\Omega \right)u_i}}={\displaystyle \underset{i\in I\left(x\right)}{\sum }{\displaystyle {\int }_{\Omega }{\varphi }_i^{\text{T}}f\text{d}\Omega }}+{\displaystyle \underset{i\in I\left(x\right)}{\sum }{\displaystyle {\int }_{{\Gamma }_t}{\varphi }_i^{\text{T}}q\text{d}\Gamma }}$ (13)

上式可简写为：

$K{U}_N=F$ (14)

其中 $K={\left[-{\displaystyle {\int }_{\Omega }\nabla {\varphi }_i\cdot \nabla {\varphi }_j\text{d}\Omega }+{\displaystyle {\int }_{\Omega }{k}^2{\varphi }_i\cdot {\varphi }_j\text{d}\Omega }\right]}_{1\le i,j\le N}$ ， $F={\left[{\displaystyle {\int }_{\Omega }{\varphi }_i^{\text{T}}f\text{d}\Omega }+{\displaystyle {\int }_{{\Gamma }_t}{\varphi }_i^{\text{T}}q\text{d}\Gamma }\right]}_{1\le i\le N}$

这样通过上式就可求解出n个在该积分节点的近似参数。但在该部分未引入边值条件，所以求得的解因为在边界条件的不同会对数值解产生不同的影响。我们在下一部分引入了拉格朗日乘子法来对边界条件进行处理，求得较为精确的数值解。

接着使用Lagrange乘子法约束本质边界条件，则可利用拉格朗日乘子 $\lambda$ 写出边界条件的积分形式为：

${\displaystyle {\int }_{{\Gamma }_u}{\lambda }^{\text{T}}\left(u-\bar{u}\right)\text{d}\Gamma }$ 。

则式(12)可重写为：

$\begin{array}{l}{\displaystyle {\int }_{\Omega }{\left(\delta u\right)}_x{u}_x+{\left(\delta u\right)}_y{u}_y+\left(\delta u\right)k^{2}u\text{d}\Omega }-{\displaystyle {\int }_{\Omega }\delta u^{\text{T}}f\text{d}\Omega }-{\displaystyle {\int }_{{\Gamma }_t}\delta u^{\text{T}}q\text{d}\Gamma }\\ -{\displaystyle {\int }_{{\Gamma }_u}\delta {\lambda }^{\text{T}}\left(u-\bar{u}\right)\text{d}\Gamma }-{\displaystyle {\int }_{{\Gamma }_u}\delta u^{\text{T}}\lambda \text{d}\Gamma }=0\end{array}$ (15)

在上式中的最后两项的目的是为了使得在我们把 $u^h$ 近似为u时会产生微小误差，为减小误差对数值模拟结果的影响，则可通过拉格朗日乘子使得边界条件 $u-\bar{u}=0$ 在任意的边界处能够成立。

其中拉格朗日乘子的表达式为： $\lambda =N\left(s\right)v$ 。

式中的 $N\left(s\right)$ 为本质边界条件上的形函数，可设矩阵 $N\left(s\right)$ 的插值节点个数为m，s为沿本质边界条件上的弧长，v为本质边界条件上的节点向量。

令 $\delta \lambda =N\left(s\right)$ ，以及将(1)代入式(15)中，可得：

$\left[\begin{array}{cc}K& G\\ G^{\text{T}}& 0\end{array}\right]\left[\begin{array}{c}u\\ v\end{array}\right]=\left[\begin{array}{c}F\\ Q\end{array}\right]$ (16)

其中的K，F，u与式(14)中的相同，而另外的G和Q分别为：

$G={\left[-{\displaystyle {\int }_{{\Gamma }_u}{N}_i^{\text{T}}{\varphi }_j\text{d}\Gamma }\right]}_{1\le j\le N,1\le i\le m}$ ， $Q={\left[-{\displaystyle {\int }_{{\Gamma }_u}{N}_i^{\text{T}}\bar{u}\text{d}\Gamma }\right]}_{1\le j\le N,1\le i\le m}$

对于上式的积分可采用高斯积分数值方法进行计算，接着通过对于所得的(16)式进行求解，最终可得所求数值解u。

4. 数值算例

为验证本文的数值模拟方法，我们考虑二维Helmholtz方程在确定的规则域中 $\Omega ={\left[0,1\right]}^2:=\left\{x:={\left(x,y\right)}^{\text{T}}\in {\text{R}}^{\text{2}}:\text{0}\le x,y\le 1\right\}$ ，具体的方程如下：

${\nabla }^{\text{2}}u+k^{2}u=f$ (17)

式中的Dirichlet边界条件可由解析解方程可得到，其解析解方程为

$u=\sin \sqrt{k^2/2}x\cos \sqrt{k^2/2}y$ (18)

在进行数值模拟时采用均匀离散节点的方案， $h=1/40$ 的离散节点图(见图1)，以及我们取=25 $h=1/40$ 时的数值k模拟图像(见图2)。

Figure 1. Discrete nodes in regular domain

图1. 规则域的离散节点

Figure 2. Meshless method for numerical solution of cloud maps

图2. 无网格法数值解云图

为更好的分析数值模拟结果的精确度和稳定性，我们将函数的结果误差erru定为：

$erru=\frac{\sqrt{{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\left(u^h\left(x_i\right)-u\left(x_i\right)\right)}^2}}}{\sqrt{{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\left(u\left(x_i\right)\right)}^2}}}$

从表1中的数据可以看出，随着节点的增多，此时的相应误差在不断的减少，我们可以看出数值模拟解在不断的趋向于解析解，说明该数值模拟方法具有良好的收敛性。

Table 1. Error analysis of different nodes

表1. 不同节点的误差结果分析

5. 结语

在本文中，运用了无网格Galerkin法来对Helmholtz方程进行离散和求解的过程中，首先通过移动最小二乘近似来构造形函数，边值条件由拉格朗日乘子法来进行引入，最后对该方程进行了数值模拟，给出的结果随着所选取节点的增加，其精确度越来越高，说明了该方法具有良好的收敛性。虽然此方法在处理Helmholtz方程时表现出较好的结果，但仍存在着许多问题，需要做进一步的研究和探讨。