摘要:基于改进的移动最小二乘插值法, 提出了黏弹性问题的插值型无单元Galerkin方法. 采用改进的移动最小二乘插值法建立形函数, 根据黏弹性问题的Galerkin弱形式建立离散方程, 推导了相应的计算公式. 与无单元Galerkin方法相比, 本文提出的黏弹性问题的插值型无单元Galerkin方法具有直接施加本质边界条件的优点. 通过数值算例讨论了影响域、节点数对计算精确性的影响, 说明了该方法具有较好的收敛性; 将计算结果与无单元Galerkin方法和有限元方法或解析解比较, 说明了该方法具有提高计算效率的优点.
关键词: 无网格方法/
改进的移动最小二乘插值法/
插值型无单元Galerkin方法/
黏弹性问题

English Abstract

全文HTML

--> --> -->

无网格方法作为一种新的数值方法, 因其只需要节点信息, 从而可以避免网格的划分以及重构^[1], 因此具有前处理简单、计算精度高等特点, 已经成为科学和工程计算发展的重要趋势之一. 当前已经发展的无网格方法有无单元Galerkin方法(element-free Galerkin method, EFG)^[1-3]、重构核粒子方法(reproducing kernel particle method, RKPM)^[4-6]、有限点方法(finite point method, FPM)^[7]、无网格局部Petrov-Galerkin方法(meshless local Petrov-Galerkin method, MLPG)^[8]、单位分解法^[9,10]、复变量无网格方法^[11-15]、维数分裂无网格方法^[16-19]和无网格的边界积分方程方法^[20-26]等.
无单元Galerkin方法是目前应用最广泛的无网格方法之一, 该方法基于移动最小二乘法构造逼近函数, 所以具有较高的计算精度. 但该方法的不足之处在于形成方程组时, 有时是病态的或者是奇异的. 因此, 有可能难以求解或者得到不正确的解. 为了改进该方法的不足, 程玉民和陈美娟^[20]以及Cheng和Peng^[21]提出了改进的移动最小二乘法. 该方法选取正交函数作为基函数, 从而致使法方程既不病态也不会奇异, 同时也不需要求解矩阵的逆, 可以直接得到方程组的解, 降低了计算量. Zhang等^[27-29]采用改进的移动最小二乘法建立形函数, 提出了瞬态热传导、波动方程和弹性动力学等问题改进的无单元Galerkin方法. Cheng和Liew^[30]以及Cheng和Wei^[31]利用改进的无单元Galerkin方法求解波动方程和广义Camassa-Holm方程. 文献[32-35]建立了三维黏弹性、三维弹塑性和弹塑性大变形等问题改进的无单元Galerkin方法, 并将改进的无单元Galerkin方法用于机场复合道面的断裂力学分析. Wu等^[36]建立了弹性力学拓扑优化问题改进的无单元Galerkin方法. Meng等^[37-39]将维数分裂法与改进的无单元Galerkin方法相结合, 提出了三维势问题、瞬态热传导问题和波动方程的杂交无单元Galerkin方法.
由于改进的移动最小二乘法构造的逼近函数不满足Kronecker $\delta $函数性质, 使得基于其形成的无网格方法不能直接施加本质边界条件. Lancaster和Salkauskas^[40]提出了移动最小二乘插值法, 该方法得到的逼近函数满足Kronecker $\delta $函数性质, 从而可以直接施加本质边界条件.
Ren和Cheng^[41,42]提出了改进的移动最小二乘插值法, 该方法证明了Lancaster的形函数公式中的一些内积可以为零, 可得到更为简单的形函数计算公式, 从而减少了方程个数以及未知量个数, 提高了计算效率. Ren和Cheng ^[41,42]基于改进的移动最小二乘插值法, 建立了势问题和弹性力学的插值型无单元Galerkin方法(interpolating element-free Galerkin method, IEFG). Cheng等^[43,44]建立了弹塑性和非线性大变形等问题的插值型无单元Galerkin方法. Deng等^[45]建立温度场问题的插值型复变量无单元Galerkin方法.
为了克服移动最小二乘插值法因权函数奇异导致的计算困难并避免产生截断误差, Wang等^[46]和Sun等^[47,48]提出了基于非奇异权函数改进的移动最小二乘插值法, 建立了势问题、弹性问题和弹塑性问题改进的插值型无单元Galerkin方法. 该方法形函数的待定系数比传统的移动最小二乘法少一个, 求逆矩阵的阶数比移动最小二乘法少一阶, 提高了计算效率和计算精度. 基于非奇异权函数改进的移动最小二乘插值法, Wang等^[49,50]研究了两点边值问题改进的插值型无单元Galerkin方法的误差估计和插值型移动最小二乘法的误差估计. Sun等^[51]分析了在n维空间的插值型移动最小二乘法的误差估计. Liu等^[52-54]建立了弹性大变形、弹塑性大变形和凝胶非均匀溶胀大变形等问题改进的插值型无单元Galerkin方法.
黏弹性问题在力学和工程中具有十分重要的应用. 由于其具有非线性的特点, 目前国内外主要是运用有限元法或者边界元法等方法来求解黏弹性问题. 但是, 近年来无网格方法发展快速, 越来越多的****运用此方法来解决黏弹性问题. Yang和Liu^[55]把无单元Galerkin方法和时域精细算法结合, 用于求解黏弹性问题. Canelas和Sensale^[56]运用无网格边界节点法分析了弹性和黏弹性材料中的简谐波问题. 彭妙娟和程玉民等建立了改进的无单元Galerkin方法^[32]、黏弹性问题复变量无单元Galerkin方法^[57]和改进的复变量无单元Galerkin方法^[58].
本文基于改进的移动最小二乘插值法构建形函数, 根据Galerkin积分弱形式建立方程, 提出了黏弹性问题的插值型无单元Galerkin方法. 通过数值算例讨论了影响域、节点数对计算精确性的影响, 说明了该方法具有较好的收敛性; 与无单元Galerkin方法以及有限元的解或解析解比较, 说明了该方法具有提高计算效率的优点.

对于改进的移动最小二乘插值法^[41,42], 定义奇异权函数

$w({{x}},{{{x}}_I}) = w({{x}} - {{{x}}_I}) = \left\{ \begin{aligned}& {\left| {{{x}} - {{{x}}_I}} \right|^{ - \alpha }}, \left\| {{{x}} - {{{x}}_I}} \right\| \leqslant {\rho _I},\\& 0,\;\;\;\;\;\;\;\;{\text{其他}},\end{aligned} \right.$

其中参数$\alpha $是正整数, ${\rho _I}$是影响域${{{x}}_I}$的影响半径.
定义任意函数$f({{x}})$和$g({{x}})$的内积为

${(f,g)_x} = \sum\limits_{I = 1}^n {w({{x}},{{{x}}_I})} f({{{x}}_I})g({{{x}}_I}), \forall f,g \in {C^0}(\bar \varOmega ),$

其中参数$\bar \varOmega = \varOmega \cup \partial \varOmega $, ${C^0}(\bar \varOmega )$是$\bar \varOmega $域中所有连续函数的集合, $n$是影响域覆盖场点$x$的节点数. 这样可得

${\left\| f \right\|_x} = {\left[ {{{(f,f)}_x}} \right]^{\frac{1}{2}}}.$

首先对基函数${p_i}({{x}})$做标准化处理. 在空间${\rm{span}} ({p_1},{p_2}, \ldots,{p_m})$中, 在${{x}}$点将基函数${p_1}({{x}}) \equiv 1$单位化,

${q_1}({{x}}) = \frac{{{p_1}}}{{{{\left\| {{p_1}} \right\|}_x}}} = \frac{1}{{{{\left[ {\sum\limits_{I = 1}^n {w({{x}} - {{{x}}_I})} } \right]}^{\frac{1}{2}}}}}. $

再将${p_2}({{x}}),{p_3}({{x}}), \ldots,{p_m}({{x}})$与${q_1}({{x}})$正交化, 可得

$\begin{split}{q_i}({{x}})\, & = {p_i}({{x}}) - {({p_i}, {q_1})_{{x}}}{q_1}({{x}})\\ &= {p_i}({{x}}) - \frac{{\sum\limits_{I = 1}^n {{p_i}({{{x}}_I})w({{x}},{{{x}}_I})} }}{{\sum\limits_{I = 1}^n {w({{x}},{{{x}}_I})} }},{\rm{ }}(i = 2,3, \cdots,m).\end{split}$

将新的基函数${q_i}({{x}})$应用于移动最小二乘法可得

${u^h}({{x}}) = {{{v}}^{\rm{T}}}({{x}}){{u}} + {{{b}}^{\rm{T}}}({{x}}){{a}}({{x}}),$

其中

${{{u}}^{\rm{T}}} = ({u_1},{u_2}, \cdots,{u_n}),$

${{{v}}^{\rm{T}}}({{x}}) = ({v_1}({{x}}),{v_2}({{x}}), \cdots,{v_n}({{x}})),$

${{{b}}^{\rm{T}}}({{x}}) = ({q_2}({{x}}),{q_3}({{x}}), \cdots,{q_n}({{x}})),$

$\begin{split}{{\alpha}} ({{x}})\,& = ({\alpha _1}({{x}}),{\alpha _2}({{x}}), \cdots,{\alpha _{m - 1}}({{x}})) \\ &= {{A}}_x^{ - 1}({{x}}){{{B}}_x}({{x}}){{u}},\end{split}$

${{{v}}_I}({{x}}) = \frac{{w({{x}},{{{x}}_I})}}{{\sum\limits_{I = 1}^n {w({{x}},{{{x}}_I})} }}, (I = 1,2, \cdots,n),$

${{{A}}_x}({{x}}) = {{{C}}^{\rm{T}}}({{x}}){{W}}({{x}}){{C}}({{x}}),$

${{{B}}_x}({{x}}) = {{{C}}^{\rm{T}}}({{x}}){{W}}({{x}}),$

${{C}}({{x}}) = \left[ {\begin{array}{*{20}{c}} {{q_2}({{{x}}_1})}&{{q_3}({{{x}}_1})}& \cdots &{{q_m}({{{x}}_1})} \\ {{q_2}({{{x}}_2})}&{{q_3}({{{x}}_2})}& \cdots &{{q_m}({{{x}}_2})} \\ \vdots & \vdots & \ddots & \vdots \\ {{q_2}({{{x}}_n})}&{{q_3}({{{x}}_n})}& \cdots &{{q_m}({{{x}}_n})} \end{array}} \right],$

${{W}}({{x}}) = \left[\!\!\!{\begin{array}{*{20}{c}} {w({{x}} - {{{x}}_1})}&0& \cdots &0 \\ 0&{w({{x}} - {{{x}}_1})}& \cdots &0 \\ \vdots & \vdots & \ddots & \vdots \\ 0&0& \cdots &{w({{x}} - {{{x}}_n})} \end{array}}\!\!\!\right].$

将(10)式代入(6)式可得

${u^h}({{x}}) = {{{v}}^{\rm{T}}}({{x}}){{u}} + {{{b}}^{\rm{T}}}({{x}}){{A}}_x^{ - 1}({{x}}){{{B}}_x}({{x}}){{u}}.$

(16)式可以写成

${u^h}({{x}}) = {{\varPhi}} ({{x}}){{u}} = \sum\limits_{I = 1}^n {{\varPhi _I}({{x}}){u_I}},$

其中

$\begin{split}{{\varPhi}} ({{x}})\, &= ({\varPhi _1}({{x}}), {\varPhi _2}({{x}}), \cdots,{\varPhi _n}({{x}})) \\ &= {{{v}}^{\rm{T}}}({{x}}) + {{{b}}^{\rm{T}}}(x){{A}}_x^{ - 1}({{x}}){{{B}}_{{x}}}({{x}})\end{split}$

为形函数.
由于改进的移动最小二乘插值法建立的形函数满足Kronecker $\delta $函数性质, 因此在建立离散方程时, 可以直接施加本质边界条件, 从而减少了计算时间, 提高了计算效率.

黏弹性问题的应力-应变关系与时间相关, 其材料参数是时间的函数. 假定${{\sigma }}(t)$在整个加载过程中恒定, 然而${{\varepsilon}} ({{t}})$和位移${{u}} ({{t}})$随着时间变化. 如果给定求解域$\varOmega $内的体力${{b}}$, 应力边界${\varGamma _t}$上的面力${\bar {{t}}}$及位移边界${\varGamma _u}$上的位移${\bar {{u}}}$(边界条件与时间无关), 下面给出二维黏弹性问题的基本方程和边界条件.
平衡方程

${{{L}}^{\rm{T}}}{{\sigma}} + {{b}} = {\rm{0}}\;\;\;\;({\text{在}}\varOmega{\text{域}}),$

其中${{L}}$是微分算子矩阵

${{L}} = \left[ {\begin{array}{*{20}{c}} {{\partial _{,1}}}&0 \\ 0&{{\partial _{,2}}} \\ {{\partial _{,2}}}&{{\partial _{,1}}} \end{array}} \right],$

${{\sigma}} $是域内任意一点的应力

${{{\sigma}} ^{\rm{T}}} = ({\sigma _{11}}, {\sigma _{22}}, {\sigma _{12}}),$

${{b}}$是域内任意一点的体力

${{{b}}^{\rm{T}}} = ({b_1}, {b_2}).$

几何方程

${{\varepsilon}} = {{Lu}}\;\;({\text{在}}\varOmega{\text{域}}),$

其中${{\varepsilon}} $和${{u}}$分别为域内任意点的应变和位移,

${{{\varepsilon}} ^{\rm{T}}} = ({\varepsilon _{11}}, {\varepsilon _{22}}, {\varepsilon _{12}}),$

${{{u}}^{\rm{T}}} = ({u_1}, {u_2}).$

物理方程

${{\sigma}} (x) = {{{M}}_1}\left\{ {\begin{aligned} e \\ e \\ 0 \end{aligned}} \right\} + {{{M}}_2}\left\{ {\begin{aligned} {{e_{11}}} \\ {{e_{22}}} \\ {{e_{12}}} \end{aligned}} \right\}, $

其中

${{{M}}_1} = \left[ {\begin{array}{*{20}{c}} {3K}&0&0 \\ 0&{3K}&0 \\ 0&0&{3K} \end{array}} \right], $

$ {{{M}}_2} = \left[ {\begin{array}{*{20}{c}} {\dfrac{1}{{{J_1}(t)}}}&0&0 \\ 0&{\dfrac{1}{{{J_1}(t)}}}&0 \\ 0&0&{\dfrac{1}{{{J_1}(t)}}} \end{array}} \right]. $

$J(t)$为不同流变模型的蠕变柔量, 对Maxwell模型,

$J(t) = \frac{1}{G} + \frac{t}{\eta }; $

对Kelvin 模型,

$J(t) = \frac{1}{G}(1 - {{\rm{e}}^{ - \frac{{Gt}}{\eta }}}); $

对三参数模型,

$J(t) = \frac{1}{{{G_2}}} + \frac{1}{{{G_1}}}(1 - {{\rm{e}}^{ - \frac{{{G_1}t}}{\eta }}}). $

其中$G$, ${G_{\rm{1}}}$和${G_{\rm{2}}}$分别是各模型中弹簧的剪切弹性模量, $\eta $是各模型中黏壶的黏性系数, $K$为体积弹性模量.
边界条件

${{u}} = \bar {{u}}\;\;\;\;({\text{在边界}}{\varGamma _u}{\text{上}}), $

${{n}} \cdot {{\sigma}} = {\bar {{t}}}\;\;\;\;({\text{在边界}}{\varGamma _t}{\text{上}}), $

其中

${{n}} = \left[ {\begin{array}{*{20}{c}} {{n_1}}&0&{{n_2}} \\ 0&{{n_2}}&{{n_1}} \end{array}} \right], $

${n_{\rm{1}}}$和${n_{\rm{2}}}$是边界点的外法线方向余弦.
由于可以直接施加本质边界条件, 因此黏弹性问题的Galerkin积分弱形式为

$\int_{ \varOmega } { \delta {{{\varepsilon}} ^{\rm{T}}}{{\sigma}} \rm{d} \varOmega - } \int_{ \varOmega } { \delta {{{u}}^{\rm{T}}}{{b}} {\rm{d}} \varOmega } - \int_{ {\varGamma _t}} {\delta {{{u}}^{\rm{T}}}\bar {{t}} {\rm{d}} {\varGamma _t} = 0}. $

设问题求解域$\varOmega $内配有$M$个节点, 节点的影响域${\varOmega _I}$($I = 1, 2, \cdots, M$)的并集覆盖了整个域$\varOmega $. 由改进的移动最小二乘插值法, 位移向量

${{u}} = {{u}}({{x}}) = {({u_1}({{x}}),{u_2}({{x}}))^{\rm{T}}} $

的插值函数为

${u_1}({{x}}) = \sum\limits_{I = 1}^n {{\varPhi _I}({{x}})} {u_1}({{{x}}_I}), $

${u_2}({{x}}) = \sum\limits_{I = 1}^n {{\varPhi _I}({{x}})} {u_2}({{{x}}_I}). $

(36)式可以写为

$\begin{split}{{u}}({{x}})\, & = \left[ {\begin{array}{*{20}{c}} {{u_1}({{x}})} \\ {{u_2}({{x}})} \end{array}} \right] \\ & = \sum\limits_{I = 1}^n {\left[ {\begin{array}{*{20}{c}} {{\varPhi _I}({{x}})}& 0 \\ 0&{{\rm{ }}{\varPhi _I}({{x}})} \end{array}} \right]} \left[ {\begin{array}{*{20}{c}} {{u_1}({{{x}}_I})} \\ {{u_2}({{{x}}_I})} \end{array}} \right]\\& = \sum\limits_{I = 1}^n {{{{N}}_I}({{x}})} {{{u}}_I} = {{N}}({{x}}){{U}},\end{split}$

其中

${{N}}({{x}}) = ({{{N}}_1}({{x}}), {{{N}}_2}({{x}}), \cdots , {{{N}}_n}({{x}})), $

${{{N}}_I}({{x}}) = \left[ {\begin{array}{*{20}{c}} {{\varPhi _I}({{x}})}&{{\rm{ }}0} \\ 0&{{\rm{ }}{\varPhi _I}({{x}})} \end{array}} \right], $

${{{U}}^{\rm{T}}} = ({{{u}}^{\rm{T}}}({{{x}}_1}), {{{u}}^{\rm{T}}}({{{x}}_2}) , \cdots , {{{u}}^{\rm{T}}}({{{x}}_n})), $

${{{u}}_I} = {{u}}({{{x}}_I}) = \left[ {\begin{array}{*{20}{c}} {{u_1}({{{x}}_I})} \\ {{u_2}({{{x}}_I})} \end{array}} \right]. $

黏弹性平面问题的应变张量可表示为

$\begin{split}{{\varepsilon}} ({{x}})\, & = \left\{ {\begin{array}{*{20}{c}} {{\varepsilon _{11}}} \\ {{\varepsilon _{22}}} \\ {{\varepsilon _{12}}} \end{array}} \right\} = {{Lu}} = \sum\limits_{I = 1}^n {{{{B}}_I}({{x}}){{u}}({{{x}}_I})} = {{B}}({{x}}){{U}}\\ & = \left[\!\!{\begin{array}{*{20}{c}} {{{\left( \sum\limits_{I = 1}^n {{\varPhi _I}({{x}}) u({{{x}}_I})}\right)}_{, 1}}} \\ {{{\left( \sum\limits_{I = 1}^n {{\varPhi _I}({{x}}) u({{{x}}_I})}\right)}_{, 2}}} \\ {{{\left(\sum\limits_{I = 1}^n {{\varPhi _I}({{x}}) u({{{x}}_I})}\right)}_{, 2}} + {{\left( \sum\limits_{I = 1}^n {{\varPhi _I}({{x}}) u({{{x}}_I})}\right)}_{, 1}}} \end{array}}\!\!\right],\end{split} $

其中

${{B}}({{x}}) = ({{{B}}_1}({{x}}),{{{B}}_2}({{x}}), \cdots ,{{{B}}_n}({{x}})), $

${{{B}}_{I}}({{x}})= \left[ {\begin{array}{*{20}{c}} { {\varPhi _{I,1}}({{x}})}&0 \\ 0&{{\varPhi _{I,2}}({{x}})} \\ {{\varPhi _{I,2}}({{x}})}&{ {\varPhi _{I,1}}({{x}})} \end{array}} \right], $

$ \begin{split}{{U}} =\, & ({u_1}({{ x}_1}), {u_2}({{ x}_1}), {u_1}({ x}_2),{u_2}({ x}_2), \\ &\cdots , {u_1}({{{x}}_n}), {u_2}({{{x}}_n}))^{\rm{T}}.\end{split}$

黏弹性平面问题的应力张量可表示为

$\begin{split}\sigma ({{x}}) \, & = \left\{ {\begin{aligned} {{\sigma _{11}}} \\ {{\sigma _{22}}} \\ {{\sigma _{12}}} \end{aligned}} \right\} = \left\{ {\begin{aligned} {{\sigma _v}} \\ {{\sigma _v}} \\ 0\;\end{aligned}} \right\} + \left\{ {\begin{array}{*{20}{c}} {{\sigma _{11}} - {\sigma _v}} \\ {{\sigma _{22}} - {\sigma _v}} \\ {{\sigma _{12}}} \end{array}} \right\} \\ & = {M_1}\left\{ {\begin{aligned} e \\ e \\ 0 \end{aligned}} \right\} + {M_2}\left\{ {\begin{aligned} {{e_{11}}} \\ {{e_{22}}} \\ {{e_{12}}} \end{aligned}} \right\},\end{split}$

其中

$\left\{ {\begin{aligned} e \\ e \\ 0 \end{aligned}} \right\} = \frac{1}{3}\left\{ {\begin{array}{*{20}{c}} {{\varepsilon _{11}} + {\varepsilon _{22}} + {\varepsilon _{33}}} \\ {{\varepsilon _{11}} + {\varepsilon _{22}} + {\varepsilon _{33}}} \\ 0 \end{array}} \right\} = {{{B}}^{(1)}}({{x}}){{U}}, $

$\left\{ {\begin{aligned} {{e_{11}}} \\ {{e_{22}}} \\ {{e_{12}}} \end{aligned}} \right\} = \left\{\!\!{\begin{array}{*{20}{c}} {{\varepsilon _{11}} - \dfrac{1}{3}({\varepsilon _{11}} + {\varepsilon _{22}} + {\varepsilon _{33}})} \\ {{\varepsilon _{22}} - \dfrac{1}{3}({\varepsilon _{11}} + {\varepsilon _{22}} + {\varepsilon _{33}})} \\ {\dfrac{1}{2}{\varepsilon _{12}}} \end{array}}\!\!\right\} = {{{B}}^{(2)}}({{x}}){{U}}, $

${{{B}}^{(1)}}({{x}}) = ({{B}}_1^{(1)}({{x}}),{{B}}_2^{(1)}({{x}}), \cdots ,{{B}}_n^{(1)}({{x}})) $

是正应变矩阵;

${{{B}}^{(2)}}({{x}}) = ({{B}}_1^{(2)}({{x}}),{{B}}_2^{(2)}({{x}}), \cdots,{{B}}_n^{(2)}({{x}})) $

是偏应变矩阵.
对平面应力情况, 即${\sigma _{33}} = 0$, 可得

${{B}}_I^{(1)}({{x}}) = \frac{1}{{3K{J_1}(t) + 2}}\left[ {\begin{array}{*{20}{c}} {{\varPhi _{I,1}}({{x}})}&{{\varPhi _{I,2}}({{x}})} \\ {{\varPhi _{I,1}}({{x}})}&{{\varPhi _{I,2}}({{x}})} \\ 0&0 \end{array}} \right], $

${{B}}_I^{(2)}({{x}}) = \left[ {\begin{array}{*{20}{c}} {{\beta _1}{\varPhi _{I,1}}({{x}})}&{ - {\beta _2}{\varPhi _{I,2}}({{x}})} \\ { - {\beta _2}{\varPhi _{I,1}}({{x}})}&{{\beta _1}{\varPhi _{I,2}}({{x}})} \\ {\dfrac{1}{2}{\varPhi _{I,2}}({{x}})}&{\frac{1}{2}{\varPhi _{I,1}}({{x}})} \end{array}} \right], $

其中${\beta _1}$和${\beta _2}$是常数,

${\beta _1} = \frac{{3K{J_1}(t) + 1}}{{3K{J_1}(t) + 2}}, $

${\beta _2} = \frac{1}{{3K{J_1}(t) + 2}}. $

令

${{B}}_I^{(1*)}({{x}}) = \left[ {\begin{array}{*{20}{c}} {{\varPhi _{I,1}}({{x}})}&{{\varPhi _{I,2}}({{x}})} \\ {{\varPhi _{I,1}}({{x}})}&{{\varPhi _{I,2}}({{x}})} \\ 0&0 \end{array}} \right] ,$

则(53)式可写为

${{B}}_I^{(1)}({{x}}) = \frac{1}{{3K{J_1}(t) + 2}}{{B}}_I^{(1*)}({{x}}). $

对平面应变情况, 即${\varepsilon _{{\rm{33}}}} = {\rm{0}}$, 可得

${{B}}_I^{(1)}({{x}}) = \frac{1}{3}\left[ {\begin{array}{*{20}{c}} {{\varPhi _{I,1}}({{x}})}&{{\varPhi _{I,2}}({{x}})} \\ {{\varPhi _{I,1}}({{x}})}&{{\varPhi _{I,2}}({{x}})} \\ 0&0 \end{array}} \right], $

${{B}}_I^{(2)}({{x}}) = \left[ {\begin{array}{*{20}{c}} {\dfrac{2}{3}{\varPhi _{I,1}}({{x}})}&{ - \dfrac{1}{3}{\varPhi _{I,2}}({{x}})} \\ { - \dfrac{2}{3}{\varPhi _{I,1}}({{x}})}&{\dfrac{1}{3}{\varPhi _{I,2}}({{x}})} \\ {\dfrac{1}{2}{\varPhi _{I,2}}({{x}})}&{\dfrac{1}{2}{\varPhi _{I,1}}({{x}})} \end{array}} \right]. $

将(49)式和(50)式代入(48)式得到

${{\sigma}} ({{x}}) = \left\{ {\begin{aligned} {{\sigma _{11}}} \\ {{\sigma _{22}}} \\ {{\sigma _{12}}} \end{aligned}} \right\} = {{{M}}_1}{{{B}}^{(1)}}{{U}} + {{{M}}_2}{{{B}}^{(2)}}{{U}}. $

将(44)式和(61)式代入(35)式得到

$\begin{split} &\int_{ \varOmega } {\delta {{{U}}^{\rm{T}}}({{{B}}^{\rm{T}}}} {{{M}}_1}{{{B}}^{(1)}}){{U}}{\rm{d}}\varOmega \\+\, & \int_{ \varOmega } {\delta {{{U}}^{\rm{T}}}({{{B}}^{\rm{T}}}} {{{M}}_2}{{{B}}^{(2)}}){{U}}{\rm{d}}\varOmega \\ -\,& \int_{ \varOmega } { \delta {{{U}}^{\rm{T}}}{{{N}}^{\rm{T}}}{{b}}{\rm{d}}\varOmega } - \int_{ {\varGamma _t}} {\delta {{{U}}^{\rm{T}}}{{{N}}^{\rm{T}}}\bar {{t}}{\rm{d}} {\varGamma _t}} = 0.\end{split}$

由$\delta {{{U}}^{\rm{T}}}$的任意性, 得到最后的离散系统方程为

${{KU}} = {{F}}, $

(63)式中, 对于平面应力问题,

$\begin{split}{{K}} =\, & \frac{1}{{3K{J_1}(t) + 2}}\int_{ \varOmega } { {{{B}}^{\rm{T}}}} {{{M}}_1}{{{B}}^{(1*)}} {\rm{d}} \varOmega \\ &+ \int_{ \varOmega } { {{{B}}^{\rm{T}}}} {{{M}}_2}{{{B}}^{(2)}} {\rm{d}} \varOmega ;\end{split}$

对于平面应变问题,

${{K}} = \int_{ \varOmega } { {{{B}}^{\rm{T}}}} {{{M}}_1}{{{B}}^{(1)}} {\rm{d}} \varOmega + \int_{ \varOmega } { {{{B}}^{\rm{T}}}} {{{M}}_2}{{{B}}^{(2)}} {\rm{d}} \varOmega,$

${{F}} = \int_{ \varOmega } {{{{N}}^{\rm{T}}}(x){{b}} {\rm{d}} \varOmega } + \int_{ \varGamma _t} {{{{N}}^{\rm{T}}}({{x}})\bar {{t}}{\rm{d}} {\varGamma _t}} .$

以上即为黏弹性问题的插值型无单元Galerkin方法.

为了验证上述理论的准确性, 以下通过该方法对4个算例进行了计算, 并将计算结果与无单元Galerkin方法以及有限元方法的计算结果或解析解进行对比. 在算例中, 构造移动最小二乘插值法的形函数采用线性基函数和奇异权函数. 在每个积分单元, 采用$4 \times 4$的Gauss积分.
定义相对误差公式为

$E = \left\| {u - {u^h}} \right\|_{{L^2}(\varOmega )}^{{\rm{rel}}} = \frac{{{{\left\| {u - {u^h}} \right\|}_{{L^2}(\varOmega )}}}}{{{{\left\| u \right\|}_{{L^2}(\varOmega )}}}}, $

其中${L^2}$范数定义为

${\left\| {u - {u^h}} \right\|_{{L^2}(\varOmega )}} = {\left( {\int_\varOmega {{{(u - {u^h})}^2}{\rm{d}}\varOmega } } \right)^{{1 / 2}}}. $

定义方差公式为

${E_{IJ}} = \frac{1}{M}\sum\limits_{k = 1}^M {[{{(u_k^I - \bar u_k^{IJ})}^2} + {{(u_k^J - \bar u_k^{IJ})}^2}]}, $

其中$u_k^I$为第$I$种节点分布时变量$u$在节点${x_k}$处的数值解, $\bar u_k^{IJ}$为第$I$种和$J$种节点分布时变量$u$在节点${x_k}$处数值解的平均值.
算例1　受均布荷载的悬臂梁
图1所示的受均布荷载作用的悬臂梁, 梁长$L = {\rm{8\;}}{\rm{m}}$, 梁高$D = {\rm{2\;}}{\rm{m}}$, 取单位厚度, 按平面应力计算. 材料体积变化是弹性的, $E = 1.0 \times 1{0^8}\; {\rm{Pa}}$, $v = 0.25$, 剪切变形流变性质满足Kelvin黏弹性模型, 其参数为$G = 2.0 \times 1{0^8} \;{\rm{Pa}}$, $\eta = 6.0 \times 1{0^8} \;{\rm{Pa}} \cdot {\rm{s}}$, 均布荷载为$p = 3.0 \times 1{0^4} \;{\rm{Pa}}$, 不计自重.
图 1 受均布荷载的悬臂梁
Figure1. A cantilever beam subjected to a distributed loading

为保证有限元法(FEM)数值解的精确性, 进行不同节点布置以得到较为精确的数值解. 采用四边形4节点单元, 具体节点布置为: $17 \times 9$, $21 \times $11, $29 \times 11$, $33 \times 11$, $33 \times 13$, $37 \times 17$, 图2给出了不同节点分布下有限元法解的方差. 从图2可以看出, 节点分布越密方差越小, 说明有限元法得到的解是收敛的. 在此选用$33 \times 13$的节点布置.
图 2 不同节点分布下有限元法解的方差
Figure2. The variances of the solutions of FEM under different node distributions.

在不同节点分布$9 \times 3$, $11 \times 4$, $13 \times 5$, $15 \times 6$, $17 \times 7$和$19 \times 8$的情况下, 利用本文提出的插值型无单元Galerkin方法计算, 图3给出了不同节点分布下数值解的相对误差. 从图3可以看出, 节点分布越密相对误差越小, 即插值型无单元Galerkin方法具有较好的收敛性. 在考虑计算效率的情况下, 采用图4所示的$17 \times 7$节点布置, 背景积分网格选取$16 \times 6$.
图 3 不同节点分布下的相对误差
Figure3. The relative error under different node distributions.

图 4 节点布置
Figure4. Node distribution.

不同影响域比例参数对本文方法数值解的精度具有一定影响, 如图5所示, 当${d_{{\rm{max}}}} = $1.7—1.9时, 该方法数值解的精度较高. 本算例选用${d_{{\rm{max}}}} = 1.7$.
图 5 不同影响域比例参数下的相对误差
Figure5. The relative error for different scale parameters of influence domains.

采用插值型无单元Galerkin方法求解时, 节点布置为$17 \times 7$, 背景积分网格选取$16 \times 6$, ${d_{{\rm{max}}}} = 1.7$, 此时中轴线各点挠度和右端中点挠度随时间变化的总相对误差为$1.29\% $, 计算时间为31.19 s.
采用无单元Galerkin方法求解时, 节点布置为$17 \times 7$, 背景积分网格选取$16 \times 6$, ${d_{{\rm{max}}}} = 3.0$, $\alpha = 2.3 \times 1{0^{10}}$, 从而可以得到与插值型无单元Galerkin方法相近的计算精度, 计算时间为34.67 s.
上述两种方法所得的数值解与有限元法的数值解的对比如图6和图7所示, 可以看出, 两种方法得到的数值解与有限元法的数值解较为符合, 但插值型无单元Galerkin方法能够提高计算效率.
图 6 $t = 20 \;{\rm{s}}$时悬臂梁中轴线上各点的挠度
Figure6. Vertical displacements of nodes on the neutral axis of the beam when $t = 20 \;{\rm{s}}$.

图 7 梁右端中点的挠度随时间的变化
Figure7. Time history of vertical displacement of midpoint in the right end of the beam.

算例2　受纯弯曲的梁
图8为受纯弯曲的梁, 几何参数为$L = 5\;{\rm{m}}$, $H = 2\;{\rm{m}}$, 单位厚度. 材料体积变化是弹性的, $E = 1.0 \times 1{0^6}\; {\rm{Pa}}$, $\nu = 0.3$, 剪切变形的流变性质满足三参数模型, 其参数为${G_1} = 5.0 \times 1{0^5}\; {\rm{Pa}}$, G₂ = $ 1.0 \times 1{0^6} \;{\rm{Pa}}$和$\eta = 2.0 \times 1{0^6} \;{\rm{Pa}} \cdot {\rm{s}}$. 其所受三角形分布荷载大小如图8所示, 不计体力, 按平面应力问题计算.
图 8 纯弯曲的梁
Figure8. A beam subjected to simple bending

由于对称性和反对称性, 只取1/4梁进行计算. 在${x_{\rm{1}}}$轴和${x_{\rm{2}}}$轴上, 位移${u_{\rm{1}}}$均为0, 为了消除结构在${x_{\rm{2}}}$方向上刚体位移, 可令原点处${u_{\rm{2}}}$为0. 有限元法采用四边形4节点单元, 节点布置为$21 \times 13$. 无单元Galerkin方法和插值型无单元Galerkin方法均采用图9所示的$15 \times 7$的节点布置, 背景积分网格选取$14 \times 6$.
图 9 节点分布
Figure9. Node distribution.

采用插值型无单元Galerkin方法求解时, ${d_{{\rm{max}}}} = 2.0$, 此时中轴线各点挠度和右端中点挠度随时间变化的总相对误差为$0.30\% $, 计算时间为21.88 s.
时采用无单元Galerkin方法求解时, ${d_{{\rm{max}}}} = 3.0$, $\alpha = 1.5 \times 1{0^8}$, 从而可以得到与插值型无单元Galerkin方法相近的计算精度, 计算时间为28.76 s.
上述两种方法所得数值解与有限元法数值解的对比如图10和图11所示. 可以看出, 插值型无单元Galerkin方法和无单元Galerkin方法都可以达到较高的精度, 但插值型无单元Galerkin方法能够提高计算效率.
图 10 $t = 30 \;{\rm{s}}$时梁中轴线上的节点挠度
Figure10. Vertical displacements of nodes on the neutral axis of the beam when $t = 30 \;{\rm{s}}$.

图 11 梁右端中点的挠度随时间$t$的变化
Figure11. Time history of vertical displacement of midpoint in the right end of the beam.

算例3　受均布内压的圆环
受均布内压的圆环如图12所示, 其内表面承受均匀分布的压力$p = 30 \;{\rm{kPa}}$. 几何参数为$a = 1\;{\rm{m}}$, $b = 5\;{\rm{m}}$. 材料体积变化是弹性的, $E = 1.0 \times 1{0^7}\; {\rm{Pa}}$, $\nu = 0.25$, 剪切变形的流变性质采用Kelvin模型. 其参数为$G = 5.0 \times 1{0^5}\; {\rm{Pa}}$, $\eta = 2.0 \times 1{0^6} \;{\rm{Pa}} \cdot {\rm{s}}$. 不计体力, 按平面应变问题计算.
图 12 受均布内压的厚壁圆筒
Figure12. Circular ring under a distributed inner pressure

在极坐标系下, 采用Kelvin模型, 圆筒径向位移随时间变化的解析解为

$\begin{split}{u_r}(r,t) =\, &\frac{{p{a^2}}}{{{b^2} - {a^2}}}\left[\left( {\frac{{3 r}}{{6K + {G_1}}}} \right) \cdot \left( {1 - {{\rm{e}}^{ - {\lambda _1}t}}} \right)\right.\\ &+\left. \left( {\frac{{{b^2}}}{{{G_1}r}}} \right)\left( {1 - {{\rm{e}}^{ - \frac{{{G_1}t}}{\eta }}}} \right)\right],\end{split}$

其中

$\lambda { _{\rm{1}}} = \frac{{{\rm{6}}K + {G_{\rm{1}}}}}{\eta }. $

(70)式表明Kelvin模型加载瞬时位移为零, 当时间足够长, 位移将趋于稳定.
如图13, 由于对称性, 只取1/4区域为研究对象. 无单元Galerkin方法和插值型无单元Galerkin方法均选用如图14所示的${\rm{1}}7 \times 12$的节点布置, 背景积分网格选取$16 \times 11$.
图 13 受均布内压1/4圆筒
Figure13. A quarter of the circular ring under a distributed inner pressure

图 14 1/4圆筒的节点分布
Figure14. Node distribution of a quarter of the circular ring.

采用插值型无单元Galerkin方法求解时, ${d_{{\rm{max}}}} = 1.01$, 此时$t = 30 \;{\rm{s}}$时沿${x_2} = {\rm{0}}$线上节点的位移和点$(2,0)$的径向位移随时间变化的总相对误差为$1.10\% $, 计算时间为23.38 s.
采用无单元Galerkin方法求解时, ${d_{{\rm{max}}}} = 1.5$, $\alpha = 8.0 \times 1{0^7}$, 此时可得到与上述插值型无单元Galerkin方法相近的相对误差, 计算时间为28.66 s.
两种方法得到的数值解与解析解的对比如图15和图16所示. 从以上分析和图15和图16的对比可看出, 插值型无单元Galerkin方法和无单元Galerkin方法的数值解均与解析解符合得较好, 但插值型无单元Galerkin方法能够提高计算效率.
图 15 $t = 30 \;{\rm{s}}$时沿${x_2} = {\rm{0}}$线上节点的位移
Figure15. Radial displacements at${x_2} = {\rm{0}}$ when $t = 30 \;{\rm{s}}$.

图 16 点$(2,0)$的径向位移随时间$t$的变化
Figure16. Time history of radial displacement at point $(2,0)$.

算例4　工程算例: 受静水压力的混凝土水坝
为了证明本方法在较复杂几何求解域的可行性, 本算例计算一个受到静水压力的混凝土水坝, 水坝的几何参数如图17所示. 坝高100 m, 上游水位90 m, 下游无水. 坝底部完全固定, 仅考虑坝体受到的上游水荷载. 坝体材料为混凝土, 所以需要考虑建造完工后期混凝土的流变效应. 假定混凝土材料体积变化是弹性的, $E = 2.0 \times 1{0^{10}}\; {\rm{Pa}}$, $\nu = 0.2$, 混凝土的早期流变性质满足三参数模型, 其参数为${G_1} = 8.33 \times 1{0^9}\; {\rm{Pa}}$, ${G_2} = 4.0 \times 1{0^{10}}\; {\rm{Pa}}$和$\eta = 8.0 \times $$ 10^{16} \;{\rm{Pa}} \cdot {\rm{s}}$. 不计体力, 按平面应变问题计算.
图 17 受静水压力的混凝土水坝
Figure17. A concrete dam under hydrostatic pressure.

有限元法采用四边形4节点单元, 共布置234个节点. 无单元Galerkin方法和插值型无单元Galerkin方法均在求解域内布置如图18所示的255个节点, 背景积分网格选取231个.
图 18 混凝土水坝的节点分布
Figure18. Node distribution of a concrete dam.

采用插值型无单元Galerkin方法求解时, ${d_{{\rm{max}}}} = 1.2$, $t = 500$ d时沿${x_1} = 15$方向上节点的水平位移和点$(15, 50)$的水平位移与时间关系的总相对误差为$0.27\% $, 计算时间为50.42 s.
采用无单元Galerkin方法求解时, ${d_{{\rm{max}}}} = 3.0$, $\alpha = 2.0 \times 1{0^{12}}$, 其总相对误差为$1.36\% $, 计算时间为94.44 s.
两种方法的数值解与有限元法的数值解的对比如图19和图20所示. 从以上分析和图19和图20的对比可以看出, 与无单元Galerkin方法相比, 插值型无单元Galerkin方法具有更高的精度和效率, 说明该方法在解决复杂工程问题时能够提高计算精度和计算效率.
图 19 $t = 500$ d时沿${x_1} = 15$方向上节点的水平位移
Figure19. Horizontal displacements at ${x_1} = 15$ when $t = 500\; {\rm d}$.

图 20 混凝土坝上点$(15, 50)$的水平位移与时间的关系
Figure20. Time history of horizontal displacement of the point $(15, 50)$.

本文基于改进的移动最小二乘插值法构造形函数, 建立了黏弹性问题的插值型无单元Galerkin方法. 由于改进的移动最小二乘插值法的形函数满足Kronecker$\delta $函数性质, 避免了传统无单元Galerkin方法要利用Lagrange乘子法或者罚函数法来施加本质边界条件, 大大减少了计算量, 从而有效地提高了在求解黏弹性问题时的计算效率. 通过数值算例讨论了影响域大小、节点数对计算精确性的影响, 说明了该方法具有较好的收敛性; 将计算结果和有限元解或解析解进行对比, 说明了插值型无单元Galerkin方法较无单元Galerkin方法在求解黏弹性问题上具有提高计算效率的优点.