基于二维线弹性带孔板的虚拟元方法

A Virtual Element Method Based on 2D Linear Elastic Plate with Hole

摘要: 本文应用虚拟元方法研究二维线弹性带孔板问题，该方法的应用克服了线弹性方程数值格式的强制性、数值解的闭锁性以及应力张量的对称性等困难，即不需要显式构造基函数，仅通过单元内部及边界的自由度来构造虚拟元函数空间，进而求出数值解，并给出了半离散和全离散格式的误差估计。通过二维带有圆孔的无限大板线弹性方程的数值计算，证明了理论分析结果的正确性，且相比于传统的有限元法，该方法在解的精确性和收敛性方面具有显著优势。

Abstract: In this paper, the virtual element method is used to study the two-dimensional linear elastic plate with hole problem, the application of this method overcomes the difficulties such as the compulsion of linear elastic equation’s numerical scheme, the locking of numerical solution, and the symmetry of stress tensor, which is no need to explicitly construct the basis function, only through unit interi-or and boundary of freedom to construct the virtual function space, then calculate the numerical solution. The error estimates for semi-discrete and fully discrete schemes are given. Through the numerical calculation of linear elastic equation of two-dimensional infinite plate with circular holes, the correctness of the theoretical analysis result is proved. Compared with traditional finite element method, this method has significant advantages in solution’s accuracy and convergence.

1. 引言

近年来，随着科技的迅猛发展，数值计算应用广泛，在现代学科中发挥着愈发重要的作用，例如材料科学、工程力学及生物科学等方面，而在科学工程和数值计算中，偏微分方程是最重要的研究内容。常见的数值求解方法主要包括有限体积法 [1]、有限差分法 [2] 和有限元法 [3] 等，其中有限元方法仅限于二维的三角形和四边形，或三维的四面体和六面体的网格剖分，这种标准形状的使用简化了单元上形状函数的构造，然而，这可能需要复杂的网格算法来生成高质量的网格或捕捉拓扑变化，使得其在处理复杂问题时存在一定的困难。2013年，F. Brezzi等人提出了Virtual Element Method (简称VEM)，即虚拟元方法 [4] [5]，该方法由有限差分法演变而来，成为求解偏微分方程的一种新的离散化方法。虚拟元方法坚实的数学理论基础、实现的简单性以及计算的效率和准确性引起了许多学者的关注，逐渐开始应用于各类问题，例如对流扩散反应方程 [6]、Quasilinear ellipse问题 [7]、二维弹性问题 [8]、二维半线性sine-Gordon方程 [9]、一般二阶椭圆问题 [10] 等等。首先，虚拟元方法可用于任意多边形和多面体网格剖分，网格处理高度灵活，并能够很好地处理网格的悬挂节点和自适应问题。与传统的有限元法等方法相比，其对网格的要求更低、适应性更好。其次，该方法无需给出形状函数的显式表达，而是通过构造虚拟元函数空间，引入相应的投影算子。最后，根据单元内部及其边界的自由度来计算刚度矩阵和右端项，进而在不计算单元内部基函数的前提下求得数值解，克服了传统有限元方法在求解线弹性方程时遇到的很多障碍。由于虚拟元方法不需要形状函数的显式表达，故实现时计算量较小。

带孔板在各类工程结构物中较为常见，在工程实际中由于某种用途，不可避免地需要对板结构开孔。而由于这些孔的存在，使得带孔板在受载荷作用下应力局部增高，从而产生了应力集中现象。众所周知，应力集中区域的应力在结构强度分析中十分关键，因此获取精准、稳定、可靠的解至关重要。而在处理开孔问题时，由于开孔区域较为复杂的几何构型将导致存在形态较差的单元，使用传统的有限元方法等都无法得到良好的计算精度，因此找到一种有效的数值求解方法至关重要。

本文在第2节给出了虚拟元方法在二维线弹性带孔板问题中的理论分析。第3节通过数值算例验证理论分析结果。第4节对本文工作进行了简要的总结及对未来的展望。

2. 线弹性孔板的VEM方法

2.1. 预备知识

在本文中，我们将遵循通常的Sobolev空间和范数的常用符号。特别地，给定 $\omega \in R^2$， $k\in \mathbb{N}$，令 $P_k\left(\omega \right)$ 表示次数不超过k次的多项式空间， $P_k\left(\omega \right)$ 表示多项式向量空间。对于一个开的有界域 $\Omega$，我们将使用 ${\left|\cdot \right|}_{s,\Omega }$ 和 ${\Vert \cdot \Vert }_{s,\Omega }$ 分别表示Sobolev空间中的 $H^s\left(\Omega \right)$ 半范数和范数，而 ${\left(\cdot ,\cdot \right)}_{0,\Omega }$ 将表示 $L^2\left(\Omega \right)$ 内积， $\otimes$ 表示向量之间的笛卡尔积。令 ${\Omega }_h$ 为区域的一个剖分， $\Gamma$ 表示 ${\Omega }_h$ 的所有单元的边f的集合， $h_E$ 表示单元的直径，h表示所有单元的直径的最大值， $\left|E\right|$ 表示单元的面积。存在一个正实数 $C_s$，使得对每个h及任意的 $E\in {\Omega }_h$，都满足以下条件：

1) 剖分是由有限数量的简单多边形构成；

2) E是相对于半径 $\rho$ 不小于 $C_s{h}_E$ 的球的每一点都是星型的；

3) E的最短边与最大直径 $h_E$ 之比不小于 $C_s$。

此外，令 $a\lesssim b$，即给定两个实数 $a,b$，存在C使得 $a\le Cb$，其中正常数C与网格大小无关，且在每次出现时不一定相同。

2.2. 连续问题

我们首先考虑含孔洞的二维线弹性平面应力问题

$\sigma =\left[\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{33}\end{array}\right]={\sigma }_0\left[\begin{array}{c}1-\frac{a^2}{r^2}\left(\frac{3}{2}\cos 2\theta +\cos 4\theta \right)+\frac{3a^4}{2r^4}\cos 4\theta \\ -\frac{a^2}{r^2}\left(\frac{1}{2}\cos 2\theta -\cos 4\theta \right)+\frac{3a^4}{2r^4}\cos 4\theta \\ -\frac{a^2}{r^2}\left(\frac{1}{2}\sin 2\theta +\sin 4\theta \right)+\frac{3a^4}{2r^4}\sin 4\theta \end{array}\right],$ (1)

其中弹性体 $\Omega \subset {ℝ}^2$ 是多边形域， ${\sigma }_r$ 表示正应力， ${\sigma }_{\theta }$ 表示剪应力，a表示圆孔半径。等式(1)的变分形式为： $\exists u\in M:H_0^1\left(\Omega \right)$，使得

${\displaystyle {\int }_{\Omega }\sigma \left(x,\nabla u\left(x\right)\right):\nabla v\left(x\right)\text{d}x}={\displaystyle {\int }_{\Omega }f\left(x\right)\cdot v\left(x\right)\text{d}x},\forall v\in W,$ (2)

其中定义矩阵 $\sigma =\sigma \left(x,\nabla u\left(x\right)\right)\in {ℝ}_{symm}^2$ 表示 $\Omega$ 上应力与应变张量的本构方程，向量函数 $u:\Omega \to {ℝ}^2$ 表示物体变位场， $\left(\cdot ,\cdot \right)$ 表示定义在 $L^2\left(\Omega \right)$ 的标量积，冒号表示矩阵之间的Frobenius内积，M表示容许位移的空间，W表示其变化的空间， $f$ 表示体载荷项。由Lax-Milgram定理可以证得连续问题(2)的解存在且唯一。

2.3. 离散问题

下面对于 $\forall E\in {\Omega }_h$，我们定义虚拟元空间

$\begin{array}{l}{V}_{h,E}:=\left\{v\in {\left[H^1\left(E\right)\cap C^0\left(E\right)\right]}^2:\Delta v=0\text{in}E,{v|}_f\in P_1\left(f\right),\forall f\in \partial E\right\},\\ V_h:=\left\{v\in V:{v|}_E\in V_{h,E},\forall E\in {\Omega }_h\right\},\end{array}$

其中 $\Delta$ 表示分量方向上的Laplace算子， $V_{h,E}$ 是一个调和函数空间，且在单元边界上是分段线性和连续的。

为此可以得到该模型问题的离散双线性形式

$a\left(v,w\right)={\displaystyle {\sum }_{E\in {\Omega }_h}{a}_E}\left(v,w\right),$ (3)

其中

$a_E\left(v,w\right)={\displaystyle {\int }_E\sigma \left(x,\nabla v\left(x\right)\right):\nabla w\left(x\right)\text{d}x},$ (4)

$a\left(v,w\right)={\displaystyle {\int }_{\Omega }{\sigma }_E}\left(\nabla v\left(x\right)\right):\nabla w\left(x\right)\text{d}x,$ (5)

则连续问题(2)可以转化为： $\exists u\in M$，使得对于 $\forall v\in W$，有

$a\left(v,w\right)={\displaystyle {\int }_{\Omega }f\left(x\right)v\left(x\right)\text{d}x}.$

对于任意的 $v\in V_{h,E}$ 中，我们可以选择以下自由度

· 在单位E的顶点上：顶点处 ${\left\{v\left(v\right)\right\}}_{v\in \partial E}$ 的值；

· 在单位E的边的中点上： $v$ 在每条边f上中点的值；

· 在单位E的内部： $\frac{1}{\left|E\right|}{\displaystyle \underset{E}{\int }\nabla v\text{d}x}$ 的值。

由于 $V_{h,E}$ 包含没有显示表达式的非多项式函数，无法直接计算每个单元的局部刚度矩阵，所以我们定义多项式空间的 $L^2$ 投影算子 ${\Pi }^0,{\Pi }_E^0$，对于 $\forall G\in L^2\left(\Omega \right)$ 定义局部算子

${{\Pi }_E^{0}G|}_E=\frac{1}{\left|E\right|}{\displaystyle {\int }_{E}G\text{d}x},\forall E\in {\Omega }_h.$ (6)

定义多项式空间的 $H^1$ 投影算子 ${\Pi }^{\nabla }:V_h\to P_k\left(E\right)$，对于 $\forall v\in V_h$，有 ${\left({\Pi }^{\nabla }v\right)}_E={\Pi }_E^{\nabla }\left({v|}_E\right)\in P_1\left(E\right)$，定义局部算子

$\{\begin{array}{l}\nabla \left({\Pi }_E^{\nabla }\left({v|}_E\right)\right)={{\Pi }_E^0\nabla v|}_E,\\ \overline{{\Pi }_E^{\nabla }\left({v|}_E\right)}=\overline{{v|}_E},\end{array}\forall E\in {\Omega }_h,$ (7)

其中，对于每个单元E，都存在唯一的线性函数 ${\Pi }_E^{\nabla }\left({v|}_E\right)$，且 $w_h=\frac{1}{N^E}{\displaystyle \underset{i=1}{\overset{N^E}{\sum }}{w}_h{v}_i}$ 表示顶点处的平均值；第

二个条件是为了处理 ${\Pi }_E^{\nabla }\left({v|}_E\right)$ 中的常数部分，它可以通过 ${v|}_E$ 的自由度来计算，所以结合(6)和(7)式可以得出以下等式

$\nabla \left({\Pi }_E^{\nabla }\left({v|}_E\right)\right)={{\Pi }_E^0\nabla v|}_E=\frac{1}{\left|E\right|}{\displaystyle \underset{E}{\int }\nabla v\text{d}x}=\frac{1}{\left|E\right|}{\displaystyle \underset{f\in \partial E}{\sum }v\otimes n_f\text{d}s}.$

另外，为保证双线性形式的稳定性及收敛性，引入带有稳定项的离散双线性形式 $S_{h,E}:V_{h,E}\times V_{h,E}\to ℝ$， $S_{h,E}\left(v_h,w_h\right)=h_E^0{\displaystyle {\sum }_{v\in \partial E}{v}_h\left(v\right)w_h}\left(v\right),\forall v_h,w_h\in V_{h,E},$ 对于 $\forall s_h,v_h,w_h\in V_{h,E}$，定义 $a_{h,E}\left(\cdot ,\cdot \right)$ 的离散形式

$\begin{array}{c}{a}_{h,E}\left(s_h;v_h,w_h\right)={\tilde{a}}_{h,E}\left(v_h,w_h\right)+{\alpha }_E\left(s_h\right)S_{h,E}\left(v_h-{\Pi }_E^{\nabla }{v}_h,w_h-{\Pi }_E^{\nabla }{w}_h\right)\\ =\left|E\right|{\sigma }_E\left({\Pi }_E^0\nabla v_h\right):{\Pi }_E^0\nabla w_h+{\alpha }_E\left(s_h\right)S_{h,E}\left(v_h-{\Pi }_E^{\nabla }{v}_h,w_h-{\Pi }_E^{\nabla }{w}_h\right),\end{array}$ (8)

此时双线性形式 $a_{h,E}\left(\cdot ;\cdot ,\cdot \right)$ 仍然是 $P_1$ 一致的，即对于每一个 $q\in P_1$ 有 $q-{\Pi }_E^{\nabla }q=0$，那么对 $\forall E\in {\Omega }_h$，我们有

$a_{h,E}\left(s_h;q,v_h\right)={\displaystyle {\int }_E{\sigma }_E\left(\nabla q\right):\nabla v_h\text{d}x},\forall s_h,v_h\in V_{h,E},\forall q\in P_1\left(E\right).$ (9)

2.4. 右端项的构造

这部分我们考虑右端项 ${\langle f,v_h\rangle }_h$，参考文献 [1] 中的虚拟元方法右端项的构造，即对于每个单元E， $\forall s_h,v_h,w_h\in V_h$，全局双线性形式 $a_h\left(\cdot ;\cdot ,\cdot \right)$ 表示为

$a_h\left(s_h;v_h,w_h\right)={\displaystyle {\sum }_{E\in {\Omega }_h}{a}_{h,E}\left(s_h;v_h,w_h\right)},$ (10)

给定一个 $s_h$，对于 $\forall v_h\in V_h$， $\exists u_h\in V_h$，使得

$a_h\left(s_h;u_h,v_h\right)={\langle f,v_h\rangle }_h,$ (11)

上述右端项展开为

${\langle f,v_h\rangle }_h={\displaystyle \underset{v\in \partial E}{\sum }{\omega }_{v}f\left(v\right)v\left(v\right)},$

当处理线性问题时，通过选择权函数 ${\omega }_v$ 来提供在E上的精确积分。此外， $s_h$ 的合理选择是 $s_h=u_h$。

2.5. 误差分析

本节中我们通过引理给出了投影、插值以及右端项三部分的误差估计，进而对离散问题(10)进行误差分析。

首先对于每一个h， $\forall E\in {\Omega }_h$ 作出以下假设：

1) 函数 $\tau \mapsto {\sigma }_E\left(\tau \right)\in C^1\left({ℝ}^{2\times 2}\right)$ ；

2) 差分 $\frac{\partial {\sigma }_E\left(\tau \right)}{\partial \tau }$ 满足

i) 对于任意的 $s\in {ℝ}^{2\times 2}$ 存在 $C_{\alpha }>0$，使得 $\frac{\partial {\sigma }_E}{\partial \tau }\left(\tau \right)s:s\ge C_{\alpha }{\Vert s\Vert }^2$ ；

ii) 对于任意的 $s,t\in {ℝ}^{2\times 2}$ 存在 $C_M>0$，使得 $\frac{\partial {\sigma }_E}{\partial \tau }\left(\tau \right)s:t\le C_M\Vert s\Vert \cdot \Vert t\Vert$。

引理1 [11]：双线性形式 $a_E\left(\cdot ,\cdot \right),a\left(\cdot ,\cdot \right),a_{h,E}\left(\cdot ,\cdot \right)$ 及 $a_h\left(\cdot ,\cdot \right)$ 已经由上述(4)，(5)，(8)和(9)式给出，如果满足上述假设，则有下列不等式成立

${\left|v_h-w_h\right|}_{1,\Omega }^2\lesssim a_h\left(s_h;v_h,v_h-w_h\right)-a_h\left(s_h;w_h,v_h-w_h\right),\forall s_h,v_h,w_h\in V_h.$ (12)

$a_E\left(v,r\right)-a_E\left(w,r\right)\lesssim {\left|v-w\right|}_{1,E}{\left|r\right|}_{1,E},\forall v,w,r\in V.$ (13)

$a_{h,E}\left(s_h;v_h,r_h\right)-a_{h,E}\left(s_h;w_h,r_h\right)\lesssim {\left|v_h-w_h\right|}_{1,E}{\left|r_h\right|}_{1,E},\forall s_h,v_h,w_h,r_h\in V_h.$ (14)

引理2 [4]：对于每一个光滑的函数w，存在一个正常数C，对每个 $s\left(1\le s\le k+1\right)$ 和每个 $w\in H^s\left(E\right)$，都会存在一个 $w_{\pi }\in P_k\left(E\right)$，使得

${\Vert w-w_{\pi }\Vert }_{0,E}+h_E{\left|w-w_{\pi }\right|}_{1,E}\le C{h}_E^s{\left|w\right|}_{s,E}.$

引理3 [4]：对每个 $s\left(2\le s\le 3\right)$ 和h，以及任意的 $E\in {\Omega }_h$ 都会存在一个正常数C和一个插值函数 $w_I\in V_{h,E}$，使得

${\Vert w-w_I\Vert }_{0,E}+h_E{\left|w-w_I\right|}_{1,E}\le C{h}_E^s{\left|w\right|}_{s,E},$

且满足

${\chi }_i\left(w-w_I\right)=0,i=1,2,\cdots ,N^{dof}.$

引理4 [4]：连续问题(2)的右端项 $\left(f,v_h\right)$ 和离散问题(10)的右端项 ${\langle f,v_h\rangle }_h$ 满足

${\langle f,v_h\rangle }_h-\left(f,v_h\right)\le {\zeta }_h{\left|v\right|}_{1,\Omega },$

其中 ${\zeta }_h$ 是满足 ${\zeta }_h\le Ch{\left|f\right|}_{0,\Omega }$ 的最小常数。

在得到以上误差估计之后，我们开始对虚拟元离散的二维线弹性带孔板问题进行误差分析。

定理1：若 $u\in M$ 是连续问题(2)的解，考虑对于任意 $s_h\in V_h$，使得

$a_h\left(s_h;u_h,v_h\right)={\langle f,v_h\rangle }_h,\forall v_h\in V_h,$

有唯一解 $u_h$，则有以下误差估计成立

${\left|u-u_h\right|}_{1,\Omega }\lesssim h\left({\left|u\right|}_{2,\Omega }+{\zeta }_h\right).$

证明：对于 $u$ 的每个插值 $u_I\in V_{h,E}$ 和 $u_{\pi }\in L^2\left(\Omega \right)$，都有 $u_{\pi |E}\in P_1\left(E\right)$，令 ${\delta }_h=u_h-u_I$，结合(12)式将离散双线形式展开，再插入插值算子 $u_{\pi }$ 可以得出

$\begin{array}{c}{\left|u_h-u_I\right|}_{{}_{1,\Omega }}^2\lesssim a_h\left(s_h;u_h,{\delta }_h\right)-a_h\left(s_h;u_I,{\delta }_h\right)\\ ={\langle f,{\delta }_h\rangle }_h-{\displaystyle \underset{E\in {\Omega }_h}{\sum }{a}_{h,E}\left(s_h;u_I,{\delta }_h\right)}\\ ={\langle f,{\delta }_h\rangle }_h-{\displaystyle \underset{E\in {\Omega }_h}{\sum }\left\{a_{h,E}\left(s_h;u_I,{\delta }_h\right)-a_{h,E}\left(s_h;u_{\pi },{\delta }_h\right)+a_{h,E}\left(s_h;u_{\pi },{\delta }_h\right)\right\}},\end{array}$

结合(4)式和(9)式我们可以推出

$a_{h,E}\left(s_h;u_{\pi },{\delta }_h\right)=a_E\left(u_{\pi },{\delta }_h\right),$ (15)

将上述(15)式代入，这样就可以得到

$\begin{array}{c}{\left|u_h-u_I\right|}_{1,\Omega }^2\lesssim {\langle f,{\delta }_h\rangle }_h-{\displaystyle \underset{E\in {\Omega }_h}{\sum }\left\{\left[a_{h,E}\left(s_h;u_I,{\delta }_h\right)-a_{h,E}\left(s_h;u_{\pi },{\delta }_h\right)\right]+a_{h,E}\left(s_h;u_{\pi },{\delta }_h\right)\right\}}\\ ={\langle f,{\delta }_h\rangle }_h-{\displaystyle \underset{E\in {\Omega }_h}{\sum }\left[a_{h,E}\left(s_h;u_I,{\delta }_h\right)-a_{h,E}\left(s_h;u_{\pi },{\delta }_h\right)\right]}\\ -{\displaystyle \underset{E\in {\Omega }_h}{\sum }\left[a_E\left(u_{\pi },{\delta }_h\right)-a_E\left(u,{\delta }_h\right)\right]}-a\left(u,{\delta }_h\right)\\ =\left[{\langle f,{\delta }_h\rangle }_h-\left(f,{\delta }_h\right)\right]-{\displaystyle \underset{E\in {\Omega }_h}{\sum }\left[a_{h,E}\left(s_h;u_I,{\delta }_h\right)-a_{h,E}\left(s_h;u_{\pi },{\delta }_h\right)\right]}\\ -{\displaystyle \underset{E\in {\Omega }_h}{\sum }\left[a_E\left(u_{\pi },{\delta }_h\right)-a_E\left(u,{\delta }_h\right)\right]},\end{array}$

接下来根据(13)和(14)式并结合插值估计可以得到

${\left|u_h-u_I\right|}_{1,\Omega }^2\lesssim \left(\underset{v_h\in V_h}{\sup }\frac{{\langle f,v_h\rangle }_h-\left(f,v_h\right)}{{\left|v_h\right|}_{1,\Omega }}+{\left|u_I-u_{\pi }\right|}_{1,\Omega }+{\left|u_{\pi }-u\right|}_{1,\Omega }\right){\left|{\delta }_h\right|}_{1,\Omega },$

回顾 ${\delta }_h=u_h-u_I$，通过化简可以得到关于 $u$ 的误差估计

${\left|u_h-u_I\right|}_{1,\Omega }\lesssim \underset{v_h\in V_h}{\sup }\frac{{\langle f,v_h\rangle }_h-\left(f,v_h\right)}{{\left|v_h\right|}_{1,\Omega }}+{\left|u_I-u_{\pi }\right|}_{1,\Omega }+{\left|u_{\pi }-u\right|}_{1,\Omega },$

由三角不等式可以得到

${\left|u-u_h\right|}_{1,\Omega }\lesssim \underset{v_h\in V_h}{\sup }\frac{{\langle f,v_h\rangle }_h-\left(f,v_h\right)}{{\left|v_h\right|}_{1,\Omega }}+{\left|u-u_I\right|}_{1,\Omega }+{\left|u-u_{\pi }\right|}_{1,\Omega },$

最后结合引理2~4并将插值估计不等式带入可以得到

${\left|u-u_h\right|}_{1,\Omega }\lesssim h\left({\left|u\right|}_{2,\Omega }+{\zeta }_h\right).$ $\square$

3. 数值算例

本节我们通过数值算例来验证上节中的理论分析结果。假设研究对象的材料是各向同性的且承受均匀的远场应力，并添加适当的约束以消除刚体模态。考虑一个长为2L、宽为L的矩形平板，其中心有一个半径为a的圆孔。该板在沿x轴方向的单轴拉伸 ${\sigma }_0$ 作用下，其弹性模量和泊松比均沿径向变化。由于该模型的载荷分布具有对称性，只对其1/2进行建模求解，如图1所示

Figure 1. Geometric model of plate with hole on the left and 1/2 model on the right

图1. 左为带孔板问题的几何模型，右为1/2模型

在极坐标 $\left(r,\theta \right)$ 中 ${\sigma }_0=1$ 的精确应力分布为

$\begin{array}{l}{\sigma }_{11}\left(r,\theta \right)=1-\frac{a^2}{r^2}\left(\frac{3}{2}\cos 2\theta +\cos 4\theta \right)+\frac{3a^4}{2r^4}\cos 4\theta ,\\ {\sigma }_{22}\left(r,\theta \right)=-\frac{a^2}{r^2}\left(\frac{1}{2}\cos 2\theta -\cos 4\theta \right)+\frac{3a^4}{2r^4}\cos 4\theta ,\\ {\sigma }_{33}\left(r,\theta \right)=-\frac{a^2}{r^2}\left(\frac{1}{2}\sin 2\theta +\sin 4\theta \right)+\frac{3a^4}{2r^4}\sin 4\theta ,\end{array}$

其中a是圆孔的半径，由弹性力学可推出位移的解析解为

$\begin{array}{l}{u}_1\left(r,\theta \right)=\frac{a}{8\mu }\left[\frac{r}{a}\left(\kappa +1\right)\cos \theta +2\frac{a}{r}\left(\left(1+\kappa \right)\cos \theta +\cos 3\theta \right)-2\frac{a^3}{r^3}\cos 3\theta \right],\\ u_2\left(r,\theta \right)=\frac{a}{8\mu }\left[\frac{r}{a}\left(\kappa -3\right)\sin \theta +2\frac{a}{r}\left(\left(1-\kappa \right)\sin \theta +\sin 3\theta \right)-2\frac{a^3}{r^3}\sin 3\theta \right],\end{array}$

其中 $\mu$ 中为剪切模量，对于 $\kappa$ Kolosov)常数可定义为：在平面应变状态下取 $\kappa =3-4\nu$，在平面应力状态下取 $\kappa =\frac{3-\nu }{1+\nu }$。

在数值实验中，参数设置为： $\Omega$ 为矩形区域 $\left[1,1\right]\times \left[1,1\right]$，圆孔半径 $a=0.2$，弹性模量 $E=2\times {10}^8\text{kpa}$，泊松比 $\nu =0.3$，在右边缘施加牵引力，左、上、下边缘应用Dirichlet齐次边界条件和Neumann边界条件，接下来对该区域进行剖分，取N = 64、256、1024、4096、16384、65536六组网格剖分，如图2所示，

为了进行误差估计，下面给出H¹和L²范数下的相对误差公式

${\Vert e\Vert }_{L^2}=\frac{{\Vert u-\Pi u_h\Vert }_{L^2}}{{\Vert u\Vert }_{L^2}}=\sqrt{\frac{{\displaystyle \underset{E\in {\Omega }_h}{\sum }{\displaystyle \underset{E}{\int }{\left(u-\Pi u_h\right)}^{\text{T}}\left(u-u_h\right)\text{d}E}}}{{\displaystyle \underset{E\in {\Omega }_h}{\sum }{\displaystyle \underset{E}{\int }{\Vert u\Vert }^2\text{d}E}}}},$

${\Vert e\Vert }_{H^1}=\frac{{\Vert u-\Pi u_h\Vert }_{H^1}}{{\Vert u\Vert }_{H^1}}=\sqrt{\frac{{\displaystyle \underset{E\in {\Omega }_h}{\sum }{\displaystyle \underset{E}{\int }{\left(\epsilon \left(u\right)-\epsilon \left(\Pi u_h\right)\right)}^{\text{T}}C\left(\epsilon \left(u\right)-\epsilon \left(\Pi u_h\right)\right)\text{d}E}}}{{\displaystyle \underset{E\in {\Omega }_h}{\sum }{\displaystyle \underset{E}{\int }\epsilon {\left(u\right)}^{\text{T}}C\epsilon \left(u\right)\text{d}E}}}}.$

Figure 2. Grid division diagram

图2. 网格剖分图

通过计算我们给出了不同网格规模下误差函数的H¹和L²范数下的相对误差，如图3所示，直观展示了相应的H¹误差和L²误差的收敛趋势。从图中我们可以看出，随着网格剖分单元尺寸的不断减小，即当N越大时，数值解的收敛性越好。当采用最细密的网格(N = 65,536)时，其收敛效果最佳。这表明利用虚拟元方法求得的误差结果与我们上述收敛性分析是一致的，进而证明了该方法的准确性。另外，如图4所示，通过对比有限元方法离散线弹性带孔板问题 [12]，可以看出虚拟元方法对问题的收敛速度更快、数值精度更好。

Figure 3. Relative error in H¹ and L² norm

图3. H¹和L²范数下的相对误差

Figure 4. Error convergence trend of finite element method for plate with hole

图4. 带孔板的有限元法误差收敛趋势图

4. 总结与展望

在本文的工作中，我们应用虚拟元方法求解二维线弹性带孔板问题，通过定义相应的虚拟元空间、选取空间中函数的自由度来计算投影算子，并基于位移场的低阶近似以及对数值位移的适当处理，得到该问题的虚拟元离散格式，给出了离散格式解的先验误差和后验误差，最后通过数值算例验证该方法的稳定性和准确性，同时也证明了该方法不仅适用于各种不同的网格类型和网格剖分，而且适用于不同的区域，相较于传统的FEM具有更优的收敛性和精确度。对于弹性问题的虚拟元方法仍存在许多问题有待解决，而非线性弹性问题虚拟元方法的构造与研究将是我们今后探究的课题。

NOTES

^*通讯作者。