杂交无单元Galerkin方法施加Dirichlet边界条件研究

Research on Application of the Dirichlet Boundary Condition Using Hybrid Element⁃Free Galerkin Method

摘要

Lagrange乘子法和罚函数法是无网格方法施加边界条件常用的两种方法， 为了比较两种方法的优缺点， 本文研究了三维Helmholtz方程的杂交无单元Galerkin（Hybrid Element-Free Galerkin， HEFG）方法。引入维数分裂法将控制方程分裂为若干个二维问题， 对于每个二维问题， 分别采用Lagrange乘子法和罚函数法施加边界条件， 建立等价的泛函， 并推导相应的积分弱形式。引入改进的移动最小二乘法建立形函数， 进而推导二维问题的离散方程。在维数分裂方向采用有限差分法将这些二维离散方程进行耦合， 得到原三维Helmholtz方程的离散求解方程。数值算例中对数值解的精度和时间进行对比， 分析了两种方法施加Dirichlet边界条件的优缺点， 得出采用罚函数法施加边界条件较好的结论。

Abstract

The Lagrange multiplier method and the penalty method are the common methods when applying essential boundary conditions in meshless method. In order to compare the advantage and the disadvantage of two methods， the hybrid element-free Galerkin （HEFG） method was presented for analyzing 3D Helmholtz equation. By introducing the dimensional split method， the governing equation could be split into a few 2D forms， for every 2D problem， the Lagrange multiplier method and the penalty method were used to apply the boundary conditions， and the equivalent functional could be established， thus the corresponding integral weak forms could be derived. By introducing the improved moving least squares （IMLS） approximation to establish shape functions， the discrete equation of 2D forms could be obtained. In dimensional split direction， the finite difference method was selected to couple these 2D equations， thus the final discrete equation of 3D Helmholtz equation was obtained. In numerical examples， by comparing the computational accuracy and computational time of numerical results， the advantages and disadvantages of two methods for applying boundary conditions were analyzed， respectively. It is shown that the penalty method is better than the Lagrange multiplier method when applying essential boundary conditions.

Graphical abstract

关键词

Lagrange乘子法 / 罚函数法 / Helmholtz方程 / 杂交无单元Galerkin方法

Key words

Lagrange multiplier method / penalty method / Helmholtz equation / hybrid element-free Galerkin method

0 引　言

无网格方法^［1］基于点的近似构造逼近函数， 是继有限元之后一种重要的数值方法， 在解决力学领域复杂的非线性大变形问题时不会出现有限元所伴随的网格畸变的麻烦。如今， 许多研究者对无网格方法产生了浓厚的兴趣。无单元Galerkin（简称EFG）方法^［2］是目前无网格方法研究和应用中最常见的一种方法， 建立逼近函数时采用移动最小二乘法^［3］（简称MLS）， 此后学者们对该逼近函数进行了一系列改进， 研究了改进的移动最小二乘法^［4］（简称IMLS）、 插值型以及复变量移动最小二乘法^［5-6］等。另外， 采用这些方法构造逼近函数， 建立了改进的无单元Galerkin（简称IEFG）方法^［7-8］、 插值型以及复变量EFG方法^［6，9］。

采用传统的EFG方法和IEFG方法对三维问题求解时计算速度较慢， 其主要原因是不同点的形函数及其导数不同， 对于三维问题来说， 每个点的影响域内的高斯点要比二维问题多得多， 因此每个高斯点都需要计算形函数和导数， 在形函数的计算过程中又涉及到矩阵求逆及多个矩阵的相乘， 远比有限元法复杂。因此， 如何提高三维问题无网格方法的计算效率是目前无网格方法需要解决的问题之一。为了解决该问题， 程珩等^［10-11］将有限差分法和改进的复变量EFG方法相结合， 提出了维数分裂复变量EFG方法， 在计算精度近似的前提下， 新方法可以大幅度提高EFG和IEFG方法求解三维问题的计算速度。孟智娟等^［12-13］提出了维数分裂EFG方法和插值型维数分裂EFG方法； 彭飘飘^［14］建立了维数分裂重构核粒子法； 王诗涵^［15］提出了维数分裂插值型EFG方法。这些方法的研究成果说明了维数分裂法是提高无网格方法求解三维问题计算效率的有效途径。

传统的EFG方法、 IEFG方法、 复变量EFG方法和重构核粒子法， 其形函数都不具有插值特性， 目前边界条件的施加最常用的是Lagrange乘子法和罚函数法。程珩等^［16］将有限差分法和IEFG方法结合， 研究了三维Helmholtz方程的HEFG方法， 该研究选择了Lagrange乘子法而并未采用罚函数法施加边界条件。

通过对HEFG方法的研究发现， 边界条件的施加对公式推导和程序编写影响比较大。因此很有必要将两种方法进行对比研究， 本文以三维Helmholtz方程的杂交无单元Galerkin方法为例， 采用维数分裂法将三维Helmholtz方程分裂为若干个二维形式， 分别采用Lagrange乘子法和罚函数法对二维问题施加边界条件， 引入改进的移动最小二乘法建立形函数， 推导二维问题的离散方程， 第三个方向采用有限差分法对二维离散方程进行耦合， 得到三维Helmholtz方程的离散求解方程。通过数值算例分析了两种方法施加边界条件时数值解的精度和速度， 并分析各自的优势和不足， 从而为HEFG方法解决科学和工程领域中的三维问题提供参考。

1 改进的移动最小二乘法

对于任一点\mathit{x}， 其逼近函数可以表示为

u^h(\mathit{x})={\displaystyle \sum _{I=\mathrm{1}}^n}{\stackrel{̑}{\Phi }}_I{u}_I=\stackrel{̑}{\mathit{\Phi }}\mathit{u}， （\mathit{x}\in \mathrm{\Omega }），（1）

其中

{\mathit{u}}^{\mathrm{T}}=(u_{\mathrm{1}},u_{\mathrm{2}},\cdots ,u_n)\mathrm{。}（2）

形函数

\stackrel{̑}{\mathit{\Phi }}={\mathit{p}}^{\mathrm{T}}(\mathit{x})\mathit{A}\mathit{B}=({\stackrel{̑}{\Phi }}_{\mathrm{1}},{\stackrel{̑}{\Phi }}_{\mathrm{2}},\cdots ,{\stackrel{̑}{\Phi }}_n),（3）

式中： {\mathit{p}}^{\mathrm{T}}(\mathit{x})为基函数向量。

\mathit{A}=\left[\begin{array}{cccc}\frac{\mathrm{1}}{(p_{\mathrm{1}},p_{\mathrm{1}})}& \mathrm{0}& \cdots & \mathrm{0}\\ \mathrm{0}& \frac{\mathrm{1}}{(p_{\mathrm{2}},p_{\mathrm{2}})}& \mathrm{0}& \mathrm{0}\\ ⋮& ⋮& \ddots & ⋮\\ \mathrm{0}& \mathrm{0}& \cdots & \frac{\mathrm{1}}{(p_n,p_n)}\end{array}\right],（4）

\mathit{B}(\mathit{x})={\mathit{P}}^{\mathrm{T}}\mathit{W}(\mathit{x}),（5）

\mathit{P}=\left[\begin{array}{cccc}{p}_{\mathrm{1}}(x_{\mathrm{1}})& p_{\mathrm{2}}(x_{\mathrm{1}})& \cdots & p_m(x_{\mathrm{1}})\\ p_{\mathrm{1}}(x_{\mathrm{2}})& p_{\mathrm{2}}(x_{\mathrm{2}})& \cdots & p_m(x_{\mathrm{2}})\\ ⋮& ⋮& \ddots & ⋮\\ p_{\mathrm{1}}(x_n)& p_{\mathrm{2}}(x_n)& \cdots & p_m(x_n)\end{array}\right],（6）

\mathit{W}=

\left[\begin{array}{cccc}w(x-x_{\mathrm{1}})& \mathrm{0}& \cdots & \mathrm{0}\\ \mathrm{0}& w(x-x_{\mathrm{2}})& \cdots & \mathrm{0}\\ ⋮& ⋮& \ddots & ⋮\\ \mathrm{0}& \mathrm{0}& \cdots & w(x-x_n)\end{array}\right],（7）

式中： w(x-x_I)为权函数； x_I为影响域内覆盖\mathit{x}的节点。

以上为改进的移动最小二乘法^［4］。

2 三维Helmholtz方程的HEFG方法

控制方程为

\mathrm{\Delta }u+{\tilde{k}}^{\mathrm{2}}\mathit{u}=f(\mathit{x})， \mathit{x}=(x_{\mathrm{1}},x_{\mathrm{2}},x_{\mathrm{3}})\in \Omega \mathrm{。}（8）

边界条件为

\mathit{u}=\overline{u}(\mathit{x})， \mathit{x}\in {\Gamma }_u,（9）

q(\mathit{x})=u_{,\mathrm{1}}{n}_{\mathrm{1}}+u_{,\mathrm{2}}{n}_{\mathrm{2}}+u_{,\mathrm{3}}{n}_{\mathrm{3}}=\overline{q}(\mathit{x}),

\mathit{x}\in {\Gamma }_q,（10）

式中： \overline{u}和\overline{q}是已知的； f(\mathit{x})为给定的函数； {\tilde{k}}^{\mathrm{2}}为波数； n_i为x_i方向边界\Gamma上的外法线， \Gamma ={\Gamma }_u\bigcup {\Gamma }_q且{\Gamma }_u\bigcap {\Gamma }_q=\mathrm{\varnothing }。

为了采用HEFG方法对该问题进行求解， 需要将式（8）转化为二维形式

\frac{{\partial }^{\mathrm{2}}{u}^{(k)}}{\partial x_{\mathrm{1}}^{\mathrm{2}}}+\frac{{\partial }^{\mathrm{2}}{u}^{(k)}}{\partial x_{\mathrm{2}}^{\mathrm{2}}}=f^{(k)}-{\tilde{k}}^{\mathrm{2}}{u}^{(k)}-\frac{{\partial }^{\mathrm{2}}{u}^{(k)}}{\partial x_{\mathrm{3}}^{\mathrm{2}}},(x_{\mathrm{1}},x_{\mathrm{2}})\in {\Omega }^{(k)}， x_{\mathrm{3}}=x_{\mathrm{3}}^{(k)}\mathrm{。}（11）

将原三维问题求解域\Omega分裂为若干个二维区域， {\Omega }^{(k)}表示\Omega的k层区域。

\Omega ={\displaystyle \underset{k=\mathrm{1}}{\overset{L}{\cup }}}\left\{{\Omega }^{(k-\mathrm{1})}\times [x_{\mathrm{3}}^{(k-\mathrm{1})},x_{\mathrm{3}}^{(k)})\right\}\bigcup {\Omega }^{(L)},（12）

u^{(k)}=u(x_{\mathrm{1}},x_{\mathrm{2}},x_{\mathrm{3}}^{(k)}),（13）

f^{(k)}=f(x_{\mathrm{1}},x_{\mathrm{2}},x_{\mathrm{3}}^{(k)})\mathrm{。}（14）

每个二维区域边界条件为

u^{(k)}={\overline{u}}^{(k)}=\overline{u}(x_{\mathrm{1}},x_{\mathrm{2}},x_{\mathrm{3}}^{(k)}),

(x_{\mathrm{1}},x_{\mathrm{2}})\in {\Gamma }_u^{(k)},（15）

q^{(k)}={\overline{q}}^{(k)}=\overline{q}(x_{\mathrm{1}},x_{\mathrm{2}},x_{\mathrm{3}}^{(k)}),

(x_{\mathrm{1}},x_{\mathrm{2}})\in {\Gamma }_q^{(k)},（16）

式中： {\Gamma }_q^{(k)}和{\Gamma }_u^{(k)}分别为自然边界和本质边界； {\Gamma }^{(k)}={\Gamma }_u^{(k)}\bigcup {\Gamma }_q^{(k)}且{\Gamma }_u^{(k)}\bigcap {\Gamma }_q^{(k)}=\mathrm{\varnothing }。

采用Lagrange乘子法施加边界条件时所形成的HEFG方法的公式推导参见文献［16］。

采用罚函数法施加边界条件， 得到二维形式的等价泛函为

{\Pi }^{\mathrm{*}}={\int }_{{\Omega }^{(k)}}\left[u\left(\frac{{\partial }^{\mathrm{2}}u}{\partial x_{\mathrm{3}}^{\mathrm{2}}}+\frac{\mathrm{1}}{\mathrm{2}}{\tilde{k}}^{\mathrm{2}}u-f\right)\right]\mathrm{d}{\Omega }^{(k)}-

{\int }_{{\Omega }^{(k)}}\frac{\mathrm{1}}{\mathrm{2}}\left[{\left(\frac{\partial u}{\partial x_{\mathrm{1}}}\right)}^{\mathrm{2}}+{\left(\frac{\partial u}{\partial x_{\mathrm{2}}}\right)}^{\mathrm{2}}\right]\mathrm{d}{\Omega }^{(k)}-

{\int }_{{\Gamma }_q^{(k)}}u\overline{q}\mathrm{d}{\Gamma }^{(k)}+\frac{\alpha }{\mathrm{2}}{\int }_{{\Gamma }_u^{(k)}}(u-\overline{u})(u-\overline{u})\mathrm{d}{\Gamma }^{(k)},（17）

式中： \alpha为罚因子。

令

\delta {\Pi }^{\mathrm{*}}=\mathrm{0},（18）

可以得到积分弱形式为

{\int }_{{\Omega }^{(k)}}\delta u\cdot {\tilde{k}}^{\mathrm{2}}u\mathrm{d}{\Omega }^{(k)}+{\int }_{{\Omega }^{(k)}}\delta u\cdot \frac{{\partial }^{\mathrm{2}}u}{\partial x_{\mathrm{3}}^{\mathrm{2}}}\mathrm{d}{\Omega }^{(k)}-

{\int }_{{\Omega }^{(k)}}\delta {(\mathit{L}u)}^{\mathrm{T}}\cdot (\mathit{L}u)\mathrm{d}{\Omega }^{(k)}-{\int }_{{\Gamma }_q^{(k)}}\delta u\cdot \overline{q}\mathrm{d}{\Gamma }^{(k)}-{\int }_{{\mathrm{\Omega }}^{(k)}}\delta u\cdot f\mathrm{d}{\Omega }^{(k)}+\alpha {\int }_{{\Gamma }_u^{(k)}}\delta u\cdot u\mathrm{d}{\Gamma }^{(k)}-\alpha {\int }_{{\Gamma }_u^{(k)}}\delta u\cdot \overline{u}\mathrm{d}{\Gamma }^{(k)}=\mathrm{0},（19）

其中

\mathit{L}(\cdot )=\left[\begin{array}{c}\frac{\partial }{\partial x_{\mathrm{1}}}\\ \frac{\partial }{\partial x_{\mathrm{2}}}\end{array}\right](\cdot )\mathrm{。}（20）

在{\Omega }^{(k)}内选取M个点{\mathit{x}}_I^{(k)}， 则{\mathit{x}}_I^{(k)}的函数值为

u_I^{}=u^{(k)}({\mathit{x}}_I^{(k)})=u({\mathit{x}}_I^{(k)},x_{\mathrm{3}}^{(k)})\mathrm{。}（21）

由改进的移动最小二乘法可以得知

u({\mathit{x}}^{(k)},x_{\mathrm{3}}^{(k)})=u^{(k)}=

{\displaystyle \sum _{I=\mathrm{1}}^n}{\stackrel{̑}{\Phi }}_I({\mathit{x}}^{(k)})u_I=\stackrel{̑}{\mathit{\Phi }}({\mathit{x}}^{(k)})\mathit{u}，（22）

其中， \mathit{u}与式（2）一致， 从而可以得到

\frac{{\partial }^{\mathrm{2}}u({\mathit{x}}^{(k)},x_{\mathrm{3}}^{(k)})}{\partial x_{\mathrm{3}}^{\mathrm{2}}}=\frac{{\partial }^{\mathrm{2}}}{\partial x_{\mathrm{3}}^{\mathrm{2}}}{\displaystyle \sum _{I=\mathrm{1}}^n}{\stackrel{̑}{\Phi }}_I{u}^{(k)}=

{\displaystyle \sum _{I=\mathrm{1}}^n}{\stackrel{̑}{\Phi }}_I\frac{{\partial }^{\mathrm{2}}{u}_I}{\partial x_{\mathrm{3}}^{\mathrm{2}}}=\stackrel{̑}{\mathit{\Phi }}\mathit{u}",（23）

\mathit{L}{u}^{(k)}={\displaystyle \sum _{I=\mathrm{1}}^n}\left[\begin{array}{c}\frac{\partial }{\partial x_{\mathrm{1}}}\\ \frac{\partial }{\partial x_{\mathrm{2}}}\end{array}\right]{\stackrel{̑}{\Phi }}_I{u}_I={\displaystyle \sum _{I=\mathrm{1}}^n}{B}_I{u}_I=B\mathit{u},（24）

其中

\mathit{u}″={\left(\frac{{\partial }^{\mathrm{2}}{u}_{\mathrm{1}}}{\partial x_{\mathrm{3}}^{\mathrm{2}}},\frac{{\partial }^{\mathrm{2}}{u}_{\mathrm{2}}}{\partial x_{\mathrm{3}}^{\mathrm{2}}},\cdots ,\frac{{\partial }^{\mathrm{2}}{u}_n}{\partial x_{\mathrm{3}}^{\mathrm{2}}}\right)}^{\mathrm{T}},（25）

\mathit{B}=({\mathit{B}}_{\mathrm{1}},{\mathit{B}}_{\mathrm{2}},\cdots ,{\mathit{B}}_n),（26）

{\mathit{B}}_I=\left[\begin{array}{c}{\stackrel{̑}{\Phi }}_{I,\mathrm{1}}({\mathit{x}}^{(k)})\\ {\stackrel{̑}{\Phi }}_{I,\mathrm{1}}({\mathit{x}}^{(k)})\end{array}\right]\mathrm{。}（27）

将式（22）， 式（23）和式（24）代入式（19）可得

\mathit{C}{\mathit{u}}^{″}+\tilde{\mathit{K}}\mathit{u}=\tilde{\mathit{F}},（28）

其中

\tilde{\mathit{K}}={\mathit{K}}_{\alpha }+{\tilde{k}}^{\mathrm{2}}\mathit{C}-\mathit{K},（29）

\tilde{\mathit{F}}={\mathit{F}}_{\mathrm{1}}+{\mathit{F}}_{\mathrm{2}}+{\mathit{F}}_{\alpha },（30）

{\mathit{K}}_{\alpha }=\alpha {\int }_{{\Gamma }_u^{(k)}}{\stackrel{̑}{\mathit{\Phi }}}^{\mathrm{T}}\stackrel{̑}{\mathit{\Phi }}\mathrm{d}{\Gamma }^{(k)},（31）

\mathit{K}={\int }_{{\Omega }^{(k)}}{\mathit{B}}^{\mathrm{T}}\mathit{B}\mathrm{d}{\Omega }^{(k)},（32）

\mathit{C}={\int }_{{\Omega }^{(k)}}{\stackrel{̑}{\mathit{\Phi }}}^{\mathrm{T}}\stackrel{̑}{\mathit{\Phi }}\mathrm{d}{\Omega }^{(k)},（33）

{\mathit{F}}_{\mathrm{1}}={\int }_{{\Omega }^{(k)}}{\stackrel{̑}{\mathit{\Phi }}}^{\mathrm{T}}f\mathrm{d}{\Omega }^{(k)},（34）

{\mathit{F}}_{\mathrm{2}}={\int }_{{\Gamma }_q^{(k)}}{\stackrel{̑}{\mathit{\Phi }}}^{\mathrm{T}}\overline{q}\mathrm{d}{\Gamma }^{(k)},（35）

{\mathit{F}}_{\alpha }=\alpha {\int }_{{\Gamma }_u^{(k)}}{\stackrel{̑}{\mathit{\Phi }}}^{\mathrm{T}}\overline{u}\mathrm{d}{\Gamma }^{(k)}\mathrm{。}（36）

对于式（28）， 采用有限差分法对{\mathit{u}}^{″}进行离散， 可以得到三维Helmholtz方程的求解方程为

\mathit{E}\widehat{\mathit{u}}=\mathit{W},（37）

其中

\mathit{E}=\frac{\mathrm{1}}{(\Delta x_{\mathrm{3}}{)}^{\mathrm{2}}}\left[\begin{array}{cccccc}\mathit{H}& \mathit{C}& & & & \\ \mathit{C}& \mathit{H}& \mathit{C}& & & \\ & \mathit{C}& \mathit{H}& \mathit{C}& & \\ & & \ddots & \ddots & \ddots & \\ & & & \mathit{C}& \mathit{H}& \mathit{C}\\ & & & & \mathit{C}& \mathit{H}\end{array}\right],（38）

\mathit{H}=(\Delta x_{\mathrm{3}}{)}^{\mathrm{2}}\tilde{\mathit{K}}-\mathrm{2}\mathit{C},（39）

\mathit{W}=\left({\left({\tilde{\mathit{F}}}^{\left(\mathrm{1}\right)}-\frac{\mathit{C}{\mathit{u}}^{(\mathrm{0})}}{(\Delta x_{\mathrm{3}}{)}^{\mathrm{2}}}\right)}^T,{\tilde{\mathit{F}}}^{{\left(\mathrm{2}\right)}^T},\cdots ,{\tilde{\mathit{F}}}^{{(L-\mathrm{2})}^T},\right.

{\left.{\left({\tilde{\mathit{F}}}^{(L-\mathrm{1})}-\frac{\mathit{C}{\mathit{u}}^{(L)}}{(\Delta x_{\mathrm{3}}{)}^{\mathrm{2}}}\right)}^{\mathrm{T}}\right)}^{\mathrm{T}},（40）

\widehat{\mathit{u}}={\left({\mathit{u}}^{{(\mathrm{1})}^{\mathrm{T}}},{\mathit{u}}^{{(\mathrm{2})}^{\mathrm{T}}},\cdots ,{\mathit{u}}^{{(L-\mathrm{1})}^{\mathrm{T}}}\right)}^{\mathrm{T}}\mathrm{。}（41）

以上即为基于罚函数法施加边界条件下的三维Helmholtz方程的HEFG方法。

3 数值算例

本节采用HEFG方法对3个数值算例进行求解。基函数选择线性基函数， 每个积分网格内选择4×4个高斯积分点， 问题求解域内布置均匀分布的节点。相对误差公式为

{‖u-u^h‖}_{L^{\mathrm{2}}(\Omega )}^{rel}=\frac{{\left({\int }_{\mathrm{\Omega }}(u-u^h{)}^{\mathrm{2}}\mathrm{d}\Omega \right)}^{\frac{\mathrm{1}}{\mathrm{2}}}}{{‖u‖}_{L^{\mathrm{2}}(\Omega )}}\mathrm{。}（42）

第1个算例的控制方程为

\mathrm{\Delta }u+u=\mathrm{s}\mathrm{i}\mathrm{n}{x}_{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}{x}_{\mathrm{3}}\cdot (\mathrm{12}{x}_{\mathrm{1}}^{\mathrm{2}}-x_{\mathrm{1}}^{\mathrm{4}})，（43）

问题求解域\Omega ={[\mathrm{0},\mathrm{\pi }]}^{\mathrm{3}}， 该算例的解析解为

u=x_{\mathrm{1}}^{\mathrm{4}}\mathrm{s}\mathrm{i}\mathrm{n}{x}_{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}{x}_{\mathrm{3}}\mathrm{。}（44）

该算例的本质边界条件通过解析解可以得到。采用HEFG方法求解， 权函数选择三次样条函数， 在x₁方向进行分裂， 每个二维区域内节点数为15×15， 分裂层数为15。当采用Lagrange乘子法施加边界条件时， d_max取1.2， 可以得到较高的精度， 误差为0.225 3%， 计算时间为4.66 s。当采用罚函数法施加边界条件时， d_max取1.21， α取1.0×10⁵， 可以得到较高的精度， 误差为0.224 2%， 计算时间为4.85 s。数值解和解析解的对比如图 1~图 3 所示。

由图 1~图 3 可以看出， 两种方法的数值解和解析解吻合得都很好， 但Lagrange乘子法的计算时间略短。

第2个算例的控制方程为

\mathrm{\Delta }u+\mathrm{100}u=

(k^{\mathrm{2}}-\mathrm{3}{\mathrm{\pi }}^{\mathrm{2}})\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{\pi }{x}_{\mathrm{2}})\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{\pi }{x}_{\mathrm{3}})\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{\pi }{x}_{\mathrm{1}})\mathrm{。}（45）

问题求解域\Omega =[\mathrm{0,1}]\times [\mathrm{0,1}]\times [\mathrm{0,1}]， 边界条件为

\frac{\partial u(\mathrm{0},x_{\mathrm{2}},x_{\mathrm{3}})}{\partial x_{\mathrm{1}}}=\frac{\partial u(\mathrm{1},x_{\mathrm{2}},x_{\mathrm{3}})}{\partial x_{\mathrm{1}}}=

u(x_{\mathrm{1}},\mathrm{0},x_{\mathrm{3}})=u(x_{\mathrm{1}},\mathrm{1},x_{\mathrm{3}})=

u(x_{\mathrm{1}},x_{\mathrm{2}},\mathrm{0})=u(x_{\mathrm{1}},x_{\mathrm{2}},\mathrm{1})=\mathrm{0}\mathrm{。}（46）

解析解为

u=\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{\pi }{x}_{\mathrm{1}})\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{\pi }{x}_{\mathrm{2}})\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{\pi }{x}_{\mathrm{3}})\mathrm{。}（47）

采用HEFG方法求解， 权函数选择三次样条函数， 在x₂或x₃方向进行分裂， 每个二维区域内节点数为19×19， 分裂层数为19。当采用Lagrange乘子法施加边界条件时， d_max取1.3， 可以得到较高的精度， 误差为0.495 2%， 计算时间为8.75 s。当采用罚函数法施加边界条件时， d_max取1.24， α取1.5×10⁵， 可以得到较高的精度， 误差为0.519 2%， 计算时间为10.52 s。数值解和解析解的对比如图 4~图 6 所示。由图 4~图 6 可以看出， 两种数值解与解析解吻合得都比较好， 但采用Lagrange乘子法施加边界条件时的计算时间比罚函数法略短。

第3个算例的控制方程为

\mathrm{\Delta }u-k^{\mathrm{2}}u=\mathrm{0}\mathrm{。}（48）

问题求解域和算例2相同， 解析解为

u=\mathrm{e}{}_{}{}^{(c_{\mathrm{1}}{x}_{\mathrm{1}}+c_{\mathrm{2}}{x}_{\mathrm{2}}+c_{\mathrm{3}}{x}_{\mathrm{3}})}.（49）

该算例的本质边界条件通过解析解可以得到。选择k=5， c₁=3， c₂=2.7。采用HEFG方法求解， 权函数选择三次样条函数， 在x₁方向进行分裂， 每个二维区域内节点数为15×15， 分裂层数为15。当采用Lagrange乘子法施加边界条件时， d_max取1.15， 可以得到较高的精度， 误差为0.279 3%， 计算时间为0.95 s。当采用罚函数法施加边界条件时， d_max取1.21， α取8.5×10²， 可以得到较高的精度， 误差为0.051 4%， 计算时间为1.13 s。

数值解和解析解的对比如图 7~图 9 所示。由图 7~图 9 可以看出， 两种数值解与解析解吻合得都很好。从本算例分析可以得知， 采用罚函数法施加边界条件时精度略高， 时间稍慢一些。

4 结　论

从本文公式推导和编程的角度来看， 采用Lagrange乘子法施加边界条件， 公式推导较为复杂， 同时也增加了MATLAB程序编写的复杂性。

从调参数的角度考虑， 节点和积分网格选定后， Lagrange乘子法仅需要通过调d_max的大小来获得较小的相对误差。若采用罚函数法， 需要通过经验不断尝试调节d_max和α两个参数的大小， 花费时间较多。

从最终数值解的精度和计算效率的角度来看， 两种方法都可以得到精度较高的数值解， 而且精度相差不大， 都可以满足精度要求。Lagrange乘子法的计算效率略占优势， 但两种方法的计算时间相差并不大。

综上所述， 在杂交无单元Galerkin方法的研究中， 采用罚函数法施加边界条件比较好。主要原因是公式推导简单， 程序编写容易， 可以满足精度和效率的使用要求。相比其优势带来的便利， 通过经验调节d_max和α两个参数所花费的时间是次要矛盾。这也是目前采用罚函数法施加边界条件被广泛应用的重要原因。

由于改进的无单元Galerkin方法的形函数不具有插值性质， 必须借助于其它方法施加边界条件， 而文献［9］研究的插值型无单元Galerkin方法可以直接施加边界条件， 在今后的研究中， 将考虑采用该方法替代改进的无单元Galerkin方法对三维Helmholtz方程进行求解。

参考文献

[1] 程玉民 ．无网格方法［M］．北京： 科学出版社， 2015.

[2] BELYTSCHKO T， LU Y Y， GU L. Element-free Galerkin methods［J］. International Journal for Numerical Methods in Engineering， 1994， 37： 229-256.

[3] LANCASTER P， SALKAUSKAS K. Surfaces generated by moving least squares methods［J］. Mathematics of Computation， 1981， 37（155）： 141-158.

[4] 陈美娟， 程玉民 ．改进的移动最小二乘法［J］．力学季刊， 2003， 24（2）： 266-272.

[5] CHEN Meijuan， CHENG Yumin. The improved moving least-squares approximation［J］. Chinese Quarterly of Mechanics， 2003， 24（2）： 266-272. （in Chinese）

[6] 任红萍， 程玉民， 张武 ．改进的移动最小二乘插值法研究［J］．工程数学学报， 2010， 27（6）： 1021-1029.

[7] REN Hongping， CHENG Yumin， ZHANG Wu. Researches on the improved interpolating moving least-squares method［J］. Chinese Journal of Engineering Mathematics， 2010， 27（6）： 1021-1029. （in Chinese）

[8] CHENG Yumin， WANG Jianfei， BAI Funong. A new complex variable element-free Galerkin method for two-dimensional potential problems［J］. Chinese Physics B， 2012， 21（9）： 090203.

[9] 蔡小杰， 彭妙娟， 程玉民 ．弹塑性大变形问题的改进的无单元Galerkin方法［J］．中国科学： 物理学 力学 天文学， 2018， 48（2）： 024701.

[10] CAI Xiaojie， PENG Miaojuan， CHENG Yumin. The improved element-free Galerkin method for elastoplasticity large deformation problems［J］. Scientia Sinica Physics， Mechanics & Astronomy， 2018， 48（2）： 024701. （in Chinese）

[11] 程珩， 彭妙娟， 程玉民 ．三维Schrödinger方程的改进的无单元Galerkin方法［J］．力学季刊， 2021， 42（1）： 14-26.

[12] CHENG Heng， PENG Miaojuan， CHENG Yumin. The improved element-free Galerkin method for 3D Schrödinger equations［J］. Chinese Quarterly of Mechanics， 2021， 42（1）： 14-26. （in Chinese）

[13] REN Hongping， CHENG Yumin. The interpolating element-free Galerkin （IEFG） method for two-dimensional elasticity problems［J］. International Journal of Applied Mechanics， 2011， 3（4）： 735-758.

[14] 程珩 ．杂交复变量无单元Galerkin方法研究［D］．上海： 上海大学， 2019.

[15] CHENG Heng， LIU Yan， LIANG Dongqiong. Analyzing 3D Helmholtz equations by using the hybrid complex variable element-free Galerkin method［J］. International Journal of Computational Materials Science and Engineering， 2023， 12（3）： 2350005.

[16] 孟智娟 ．三维问题的维数分裂无单元Galerkin方法研究［D］．上海： 上海大学， 2019.

[17] MENG Zhijuan， CHI Xiaofei. An improved interpolation dimension split element-free Galerkin method for 3D wave equations［J］. Engineering Analysis with Boundary Elements， 2022， 134： 96-106.

[18] 彭飘飘 ．三维问题的杂交重构核粒子法［D］．上海： 上海大学， 2021.

[19] 王诗涵 ．三维波动方程和三维弹性力学问题的插值型维数分裂无单元Galerkin方法［D］．上海： 上海大学， 2022.

[20] CHENG Heng， ZHANG Jiao， XING Zebin. The hybrid element-free Galerkin method for 3D Helmholtz equations［J］. International Journal of Applied Mechanics， 2022， 14（9）： 2250084.

基金资助

山西省青年基金资助项目(20210302124388)

山西省创新训练项目(20230712)