弹性体力学是固体力学的一个分支,是研究弹性体在外力作用或者温度变化等的影响下,弹性体内部所引起的应力、应变和位移方面的变化,广泛应用于材料力学、机械工程、航空航天等领域中,对其理论及数值方法的研究具有重要的理论和实际价值.求解弹性问题最常用的方法是有限元方法,参见文献[1-5]等.但在一些实际工程的计算中,如橡胶、塑料等相当多的材料会呈现出几乎不可压缩的情况,即泊松比$\nu\rightarrow{0.5}$,采用有限元近似求解无法达到最优收敛阶,甚至不再收敛,即发生所谓的“闭锁”现象. Babu$\check{\rm s}$ka与Suri在文献[6-7]中详述了弹性问题中的这一现象. "闭锁"现象的产生与否与$\lambda$值的选取密切相关,文献[3]等引入一个新的变量将弹性问题转化为一个鞍点问题,不仅可以得到较好的位移近似,而且还可以得到压力或应力的较好近似,但与此同时,需要选取的有限元空间满足LBB条件.我们将Uzawa迭代和多重网格方法耦合,建立多重网格和套迭代多重网格Uzawa型混合有限元方法,不仅可以克服"闭锁"现象,而且可以达到优化初值选择、减少迭代步数、有限元空间选择自由、提高运算效率等效果.

线弹性模型是一般的非线性弹性理论的简化模型,是研究固体在受到附加外力的时候发生形变的变化情况.当弹性体发生微小形变时, Hooke定律给出了应力张量$\sigma(\vec{u})=\sigma_{ij}(\vec{u})$与应变张量$\varepsilon(\vec{u})=\varepsilon_{ij}(\vec{u})$在笛卡尔坐标系下的本构关系为

其中$c_{ijkl}$称为胡克张量,具有对称性,且满足$c_{ijkl}=c_{klij}=c_{jikl}=c_{ijlk}$.

若$c_{ijkl}(x)$为常量(即弹性体是各向同性且同质的),则本构方程为

其中$\delta_{ij}$是Kronecker符号函数, $\lambda, \mu$是描述固体力学性质的Lamé常数, $\lambda=\frac{E\nu}{(1+\nu)(1-2\nu)},$$\mu=\frac{E}{2(1+\nu)}$, $E$是杨氏模量, $\nu\in[0, 0.5)$为泊松比, $\varepsilon_{ij}(\vec{u})=\frac{1}{2}(\frac{\partial{u_i}}{\partial{x_j}} +\frac{\partial{u_j}}{\partial{x_i}}).$

设弹性体所在的多边形区域$\Omega\subset{\Bbb R}^d(d=2, 3)$,由力的平衡法则知

其中$\sigma(\vec{u})=2\mu\varepsilon(\vec{u}) +\lambda\nabla\cdot{\vec{u}} {\cal I}$表示应力张量, $\vec{u}:\Omega\rightarrow{{\Bbb R}^d}$表示位移, $\vec{f}:\Omega\rightarrow{{\Bbb R}^d}$表示给定的载荷密度, ${\cal I}$表示$d$阶单位矩阵, $\varepsilon(\cdot):{\Bbb R}^d \rightarrow{\Bbb R}^{d\times{d}}$表示应变张量

假设$0<\mu_0\leq{\mu}\leq{\mu_1}<+\infty$, $0<\lambda_0\leq{\lambda}<+\infty$,一般情况下,当泊松比$0\leq\nu\leq0.4$时,称对应的弹性问题(2.3)-(2.4)为可压缩弹性问题;当$0.4<\nu<0.5$时,称相应的问题为几乎不可压缩弹性问题.

令${\bf V}=(H_0^1(\Omega))^d$, $B(\vec{u}, \vec{v}) =\int_\Omega(2\mu{\varepsilon(\vec{u})}:\varepsilon(\vec{v}) +\lambda{\nabla\cdot\vec{u}}\nabla\cdot\vec{v})$,其中"$:$"定义为$A:B=\sum\limits_{i, j=1}^n a_{ij}b_{ij}$, $A=(a_{ij})n\times n, B=(b_{ij})n\times n$,则问题(2.3)-(2.4)的变分形式为:求$\vec{u}\in{{\bf V}}$,使得

设${\bf V}_h$为${\bf V}$的有限维子空间,则问题(2.6)的离散格式可以表示为:求$\vec{u}_h\in{\bf V}_h$,使得

由Céa引理,可知

定理2.1  设$\vec{u}, \vec{u}_h$分别为问题(2.6)和(2.7)的解,则

其中$|\cdot|=B(\cdot, \cdot)^{\frac{1}{2}}$.

由(2.8)式可知:当$\lambda\rightarrow+\infty(\nu\rightarrow0.5)$时,即便$h$很小,也得不到理想的收敛结果.为了克服这一困难,引入"压力"变量: $p=\lambda{\rm div(\vec{u})}$,将问题(2.3)-(2.4)转化为如下的边值问题

因此,当弹性体为完全不可压缩时, (2.10)式就变成${\rm div}\vec{u}=0$,这样,问题(2.9)-(2.11)就变成广义Stokes问题.

问题(2.9)-(2.11)的变分形式为:求$(\vec{u}, p)\in{{\bf V}\times(P=L_0^2(\Omega))}$,使得

其中$a(\cdot, \cdot)$, $b(\cdot, \cdot)$和$c(\cdot, \cdot)$分别为

定义算子${\cal A}: {\bf V}\rightarrow{\bf V}^\ast, {\cal B}: {\bf V}\rightarrow{P^\ast}$, ${\cal B}^\ast: P\rightarrow{\bf V}^\ast$以及${\cal C}:P\rightarrow{P^\ast}$满足

定义等距逆算子${\cal A}^{-1}: {\bf V}^\ast\rightarrow{\bf V}$和${\cal C}^{-1}: P^\ast\rightarrow{P}$分别满足

因此,变分形式(2.12)-(2.13)可以改写成算子形式

设${\bf V}_h$和$P_h$分别为${\bf V}$和$P$的子空间,对任意的$\vec{u}_h, \vec{v}_h\in{\bf V}_h$和$p_h, q_h\in{P_h}$,定义算子${\cal A}_h=-2\mu{\rm div(\varepsilon(\cdot))}: {\bf V}_h\rightarrow{\bf V}_h^\ast, {\cal B}_h=\nabla\cdot: {\bf V}_h\rightarrow{P_h^\ast}$, ${\cal B}_h^\ast=-\nabla: P_h\rightarrow{\bf V}_h^\ast$以及${\cal C}_h={\cal I}: P_h\rightarrow{P_h^\ast}$分别为

设$\Omega\subset{\Bbb R}^d(d=2, 3)$为有界多边形区域, ${\cal T}_h$为其上的拟一致三角剖分,问题(2.12)-(2.13)的混合离散格式为:求$(\vec{u}_h, p_h)\in{\bf V}_h\times{P_h}$,使得

令双线性形式$B(\vec{u}_h, p_h;\vec{v}_h, q_h)=a(\vec{u}_h, \vec{v}_h) +b(\vec{v}_h, p_h)-b(\vec{u}_h, q_h)+\lambda^{-1}c(p_h, q_h)$,则离散格式(3.1)-(3.2)等价于

易知$B(\vec{u}_h, p_h;\vec{v}_h, q_h)$满足

由Lax-Milgram定理可知, (3.3)式存在唯一解等价于$a(\cdot, \cdot), c(\cdot, \cdot)$在${\bf V}_h\times{\bf V}_h$与$P_h\times{P_h}$上的有界性和强制性,以及$b(\cdot, \cdot)$在${\bf V}_h\times{P_h}$上的有界性,对于可压缩弹性问题这一理论是成立的,但对于几乎不可压缩弹性问题,当$\lambda$值大到一定程度时, (3.2)式接近自由散度约束条件${\rm div}\vec{u}_h=0$,原问题就转化为一个广义Stokes问题,这一理论就不再成立了,从后面的数值算例可以看到这一现象.

关于Uzawa型混合有限元方法有如下收敛性结果,详细证明参见文献[8].

定理3.1  若$0<\alpha<\frac{2}{M^2}$,设$(\vec{u}_j, p_j)$为Uzawa型混合有限元方法的数值解,则

其中$M=\max \{ M_h\}$, $M_h^2$为${\cal C}_h^{-1}{\cal B}_h {\cal A}_h^{-1}{\cal B}_h^\ast$的谱半径的上界, $\gamma_h=\|{\cal I}-\alpha{\cal C}_h^{-1} {\cal B}_h{\cal A}_h^{-1}{\cal B}_h^\ast -\alpha\lambda^{-1}\|$, $C$与$h$无关,而与$\lambda$有关, $j$为迭代次数.

令${\bf V}_{k-1}\subset{\bf V}_k, P_{k-1}\subset{P_k}$,其中, ${\bf V}_k, P_k$和${\bf V}_{k-1}, P_{k-1}$分别是基于三角剖分${\cal T}_k$和${\cal T}_{k-1}$建立的有限元子空间,且$h_{k}=\frac{h_{k-1}}{2},$$k=1, 2, \cdots$.

定义$L^2$投影算子$Q_{k-1}:{\bf V}_{k}\rightarrow{\bf V}_{k-1}$满足$(Q_{k-1}\vec{v}_k, \vec{w}_{k-1})=(\vec{v}_k, \vec{w}_{k-1}),$$\forall\vec{v}_k\in{\bf V}_k,$$\vec{w}_{k-1}\in{\bf V}_{k-1}$, $R_{k-1}:P_k\rightarrow{P_{k-1}}$满足$(R_{k-1}p_k, q_{k-1})=(p_k, q_{k-1}), \forall{p_k}\in{P}_k, q_{k-1}\in{P}_{k-1}$.

算法3.1  给定初值$\vec{u}_0, p_0$及松弛系数$\alpha$,令$j=1, 2, \cdots, m_1+m_2+1$,问题(3.1)-(3.2)在有限元子空间对${\bf V}_k\times{P_k}$上的多重网格迭代解$MG(k, (\vec{u}_0, p_0), (\vec{f}, 0))$需以下几步求得

(1)前光滑过程:对$1\leq{j}\leq{m_1}$,求解

(2)校正过程:令$\vec{r}_{j}=\vec{f}-{\cal A}\vec{u}_{j} -{\cal B}^{\ast}p_{j}, s_{j}=\lambda^{-1}{\cal C}p_{j} -{\cal B}\vec{u}_{j}$,将前光滑过程得到的结果分别投影到${\bf V}_{k-1}$和$P_{k-1}$上,即

并在${\bf V}_{k-1}\times{P_{k-1}}$上,求解

其解记为$(\vec{u'}_i, p'_i)=MG(k-1, (\vec{u'}_{i-1}, p'_{i-1}), (\vec{\widehat{f}}, q)),$其中$\vec{u}'_0={\bf 0}, p'_0=0, 1\leq{i}\leq{l}.$再令

(3)后光滑过程:迭代步数为$m_2$步,对$m_1+1\leq{j}\leq{m_1+m_2}$,有

该算法称为多重网格Uzawa型混合有限元方法(MUMFEM).

取$m_1=m, m_2=0$,多重网格Uzawa型混合有限元方法的收敛性结果:

定理3.2  设问题(2.12)-(2.13)的真解为$(\vec{u}, p)\in(H^{\kappa+1}(\Omega)\cap{\bf V}, H^{\kappa+1}(\Omega)\cap{P})$,若$\alpha\in(0, \frac{2}{M^2})$,光滑次数$m$充分大,则

证  记问题(3.1)-(3.2)在$({\bf V}_3, P_3)$上的有限元近似解为$(\vec{u}_h, p_h),$由于前光滑过程使用Uzawa型混合有限元方法,所以$(\vec{u}_m, p_m)$满足

第二层网格$({\bf V}_{2}, P_{2})$上也使用Uzawa型混合有限元方法,设问题(3.7)-(3.8)在$({\bf V}_{2}, P_{2})$上的有限元近似解为$(\vec{u}'_h, p'_h)$,也满足

由(3.7)-(3.8)式可知

因此,有

上式联立(3.12)和(3.13)式并利用定理3.1,得

若光滑次数$m$充分大时,多重网格Uzawa型混合有限元方法收敛且满足估计式(3.11).证毕.

比较(3.4)式与(3.11)式这两个估计式,可以看出将Uzawa型混合有限元方法的结果在粗网格进行校正,虽然无法克服迭代初值$p_0$对结果的影响($k\geq3$),但因多了$\gamma_h^{m-1}h^{k+1}|p|_{k+1}$这一交叉项可以大大改善$p$的误差结果.

为了进一步消除初值选择对于误差结果的影响,基于算法3.1建立套迭代多重网格Uzawa型混合有限元方法:

算法3.2  给定迭代初始值$\vec{u}_0, p_0$及松弛系数$\alpha$,对问题(3.1)-(3.2),给定一个初始三角剖分${\cal T}_1$,并由此建立有限元子空间${\bf V}_1, P_1$,每次网格加细策略是将每个三角形单元剖分成四个同样的三角形得到剖分${\cal T}_k, k=2, 3, \cdots, r$,且记由${\cal T}_k$建立的有限元子空间为${\bf V}_k, P_k$,则循环过程如下:

(1)对第一层网格${\cal T}_1$,求解鞍点问题

得到第一层网格上的解$(\vec{\widehat{u}}_{1}, \widehat{p}_{1})$;

(2)对$k\geq2$,以上层迭代解作为迭代初值,并按照算法3.1进行$r$次迭代,也即按照如下循环过程定义$(\vec{\widehat{u}}_{k}, \widehat{p}_{k})$:

(ⅰ) $(\vec{u}_k^0, p_k^0)=(\vec{\widehat{u}}_{k-1}, \widehat{p}_{k-1})$;

(ⅱ) $(\vec{u}_k^i, p_k^i)=MG(l, (\vec{u}_k^{i-1}, p_k^{i-1}), (\vec{f}, 0)), 1\leq{i}\leq{r}$;

(ⅲ) $(\vec{\widehat{u}}_{k}, \widehat{p}_{k})=(\vec{u}_k^r, p_k^r)$;

(3)令$k=k+1$,按照步骤(2)进行循环.

对算法3.2,有如下的收敛结果:

定理3.3  记问题(2.12)-(2.13)的真解为$(\vec{u}, p)\in(H^{\kappa+1}(\Omega)\cap{\bf V}, H^{\kappa+1}(\Omega)\cap{P})$,网格层数$k\geq2$,则

证  由于在${\cal T}_1$上精确求解,所以对第一层网格, $(\vec{\widehat{u}}_1, \widehat{p}_1)$满足

而对第二层网格,若$i=1,$由于$(\vec{u}_k^1, p_k^1)=MG(l, (\vec{u}_k^{0}, p_k^{0}), (\vec{f}, 0))$,所以满足误差估计式(3.11),即

对$i=2,$有$(\vec{u}_k^2, p_k^2)=MG(l, (\vec{u}_k^2, p_k^2), (\vec{f}, 0))$,

这样继续下去,由于$(\vec{\widehat{u}}_2, \widehat{p}_2)$是多重网格Uzawa型有限元算法进行$r$次得到的迭代解,所以$(\vec{\widehat{u}}_2, \widehat{p}_2)$满足估计式(3.14).

对第$k$层网格${\cal T}_k(k>2),$将初值设为$(\vec{\widehat{u}}_{k-1}, \widehat{p}_{k-1}),$采用同样的递推方式,可以得到多重网格迭代解$(\vec{\widehat{u}}_k, \widehat{p}_k)$满足误差估计式(3.14).证毕.

取$\Omega$为单位方形$(0, 1)\times{(0, 1)}$,位移精确解为$\vec{u}=(u_1, u_2)^T$,其中

加载外力项$\vec{f}=(f_1, f_2)^T$,其中$f_1$和$f_2$分别为

取线性协调有限元空间${\bf V}_h=\{{\bf v}_h\in{[C(\overline{\Omega})]^2} |{\bf v}_h|_K\in{[{\cal P}_1]^2}, \forall{K\in{\cal T}_h}\}$,取网格大小为$h=h_K=\frac{\sqrt{2}}{2^4},$另设参数$\mu=1, \ \lambda= 0.1, 1, 10^3$, $10^5$, $10^{10}$,可以得到如表 4.1中的结果. 表 4.1验证了定理2.1中关于误差的相关结论,当$\lambda$较大时,会产生"闭锁"现象.

下面给出选取${\cal P}_1-{\cal P}_1$元, $\lambda=10^{17}$时几乎不可压缩弹性问题的误差结果.

取$\lambda=10^{20}$,其他参数不变,迭代初值取为$p_0=0,$迭代次数$j=10,$使用Uzawa型混合有限元方法,计算结果见表 4.3.

值得一提的是:从定理3.1中估计式(3.4)可以看出,对Uzawa型混合有限元方法(参见文献[8]),若不能选取一个合适的初值$p_0$,在同样的情况下$\vec{u}_h$仍然可以达到最优收敛阶,但$p_h$就无法达到最优收敛阶,甚至不再收敛.如在不改变其它参数的情况下,取初值为$p_0=0.01$, $\lambda=10^{10}$分别迭代10, 30, 100, 500次时的平均收敛阶情况如表 4.4所示.

针对Uzawa型混合有限元方法在求解几乎不可压缩弹性问题时对迭代初值$p_0$限制,分别使用多重网格Uzawa型混合有限元方法(MUMFEM) (迭代20次,取$m_1=10, m_2=0, l=1$)和套迭代多重网格Uzawa型混合有限元方法(IMUMFEM),其他参数不变,可以去掉对$p_0$的取值限制,计算结果如表 4.5所示.

仍使用套迭代多重网格Uzawa型混合有限元方法对表 4.5中出现的问题进一步改进:设迭代层数为5,取$m_1=10, m_2=0, l=1$,并令$r=1$,每层有校正过程和无校正过程的结果如下.

通过对比表 4.5与表 4.6可以看出,有校正过程的套迭代算法对位移$\vec{u}$和压力$p$的误差可以达到最优收敛阶,与定理3.2中的误差估计一致.

本文对几乎不可压缩线性弹性问题的数值方法进行了深入研究,提出了多重网格和套迭代多重网格Uzawa型混合有限元方法,这两类方法对所有有限元空间对都稳定,可以克服"闭锁"现象.并用数值算例验证了理论结果.