黏弹性问题的插值型无单元Galerkin方法*

张鹏轩 彭妙娟

(上海大学土木工程系, 上海 200444)

1 引 言

无网格方法作为一种新的数值方法, 因其只需要节点信息, 从而可以避免网格的划分以及重构[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δ函数性质, 使得基于其形成的无网格方法不能直接施加本质边界条件.Lancaster和Salkauskas[40]提出了移动最小二乘插值法, 该方法得到的逼近函数满足Kroneckerδ函数性质, 从而可以直接施加本质边界条件.

Ren和Cheng[41,42]提出了改进的移动最小二乘插值法, 该方法证明了Lancaster的形函数公式中的一些内积可以为零, 可得到更为简单的形函数计算公式, 从而减少了方程个数以及未知量个数,提高了计算效率.Ren和Cheng[41,42]基于改进的移动最小二乘插值法, 建立了势问题和弹性力学的插值型无单元Galerkin方法(interpolating elementfree 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方法以及有限元的解或解析解比较, 说明了该方法具有提高计算效率的优点.

2 改进的移动最小二乘插值法

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

其中参数α是正整数,ρI是影响域xI的影响半径.

定义任意函数f(x) 和g(x) 的内积为

首先对基函数pi(x) 做标准化处理.在空间span(p1,p2,...,pm)中, 在x点将基函数p1(x)≡1单位化,

再将p2(x),p3(x),...,pm(x)与q1(x) 正交化,可得

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

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

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

3 黏弹性问题的基本方程

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

平衡方程

其中L是微分算子矩阵

σ是域内任意一点的应力

b是域内任意一点的体力

几何方程

其中ε和u分别为域内任意点的应变和位移,

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

对Kelvin 模型,

对三参数模型,

其中G,G1和G2分别是各模型中弹簧的剪切弹性模量,η是各模型中黏壶的黏性系数,K为体积弹性模量.

边界条件

n1和n2是边界点的外法线方向余弦.

由于可以直接施加本质边界条件, 因此黏弹性问题的Galerkin积分弱形式为

4 黏弹性问题的插值型无单元Galerkin方法

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

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

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

是正应变矩阵;

是偏应变矩阵.

对平面应力情况, 即σ33=0 , 可得

其中β1和β2是常数,

则(53)式可写为

对平面应变情况, 即ε33=0 , 可得

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

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

由δUT的任意性, 得到最后的离散系统方程为

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

对于平面应变问题,

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

5 数值算例

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

定义相对误差公式为

其中L2范数定义为

定义方差公式为

算例1受均布荷载的悬臂梁

图1所示的受均布荷载作用的悬臂梁, 梁长L=8m , 梁高D=2m , 取单位厚度, 按平面应力计算.材料体积变化是弹性的,E=1.0×108Pa ,v=0.25, 剪切变形流变性质满足Kelvin黏弹性模型, 其参数为G=2.0×108Pa ,η=6.0×108Pa·s ,均布荷载为p=3.0×104Pa , 不计自重.

图1 受均布荷载的悬臂梁Fig.1.A cantilever beam subjected to a distributed loading.

为保证有限元法(FEM)数值解的精确性, 进行不同节点布置以得到较为精确的数值解.采用四边形4节点单元, 具体节点布置为: 1 7×9 ,21×11, 2 9×11 , 3 3×11 , 3 3×13 , 3 7×17 , 图2给出了不同节点分布下有限元法解的方差.从图2可以看出, 节点分布越密方差越小, 说明有限元法得到的解是收敛的.在此选用 3 3×13 的节点布置.

图2 不同节点分布下有限元法解的方差Fig.2.The variances of the solutions of FEM under different node distributions.

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

图3 不同节点分布下的相对误差Fig.3.The relative error under different node distributions.

图4 节点布置Fig.4.Node distribution.

不同影响域比例参数对本文方法数值解的精度具有一定影响, 如图5所示, 当dmax=1.7-1.9时, 该方法数值解的精度较高.本算例选用dmax=1.7.

图5 不同影响域比例参数下的相对误差Fig.5.The relative error for different scale parameters of influence domains.

采用插值型无单元Galerkin方法求解时, 节点布置为 1 7×7 , 背景积分网格选取 1 6×6 ,dmax=1.7, 此时中轴线各点挠度和右端中点挠度随时间变化的总相对误差为 1.29% , 计算时间为31.19 s.

采用无单元Galerkin方法求解时, 节点布置为 1 7×7 , 背景积分网格选取 1 6×6 ,dmax=3.0 ,α=2.3×1010, 从而可以得到与插值型无单元Galerkin方法相近的计算精度, 计算时间为34.67 s.

上述两种方法所得的数值解与有限元法的数值解的对比如图6和图7所示, 可以看出, 两种方法得到的数值解与有限元法的数值解较为符合, 但插值型无单元Galerkin方法能够提高计算效率.

图6 t=20s 时悬臂梁中轴线上各点的挠度Fig.6.Vertical displacements of nodes on the neutral axis of the beam when t=20s.

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

算例2受纯弯曲的梁

图8为受纯弯曲的梁, 几何参数为L=5m ,H=2m, 单位厚度.材料体积变化是弹性的,E=1.0×106Pa ,ν=0.3 , 剪切变形的流变性质满足三参数模型, 其参数为G1=5.0×105Pa ,G2=1.0×106Pa 和η=2.0×106Pa·s.其所受三角形分布荷载大小如图8所示, 不计体力, 按平面应力问题计算.

图8 纯弯曲的梁Fig.8.A beam subjected to simple bending.

由于对称性和反对称性, 只取1/4梁进行计算.在x1轴和x2轴上, 位移u1均为 0, 为了消除结构在x2方向上刚体位移, 可令原点处u2为0.有限元法采用四边形4节点单元, 节点布置为 2 1×13.无单元Galerkin方法和插值型无单元Galerkin方法均采用图9所示的 1 5×7 的节点布置, 背景积分网格选取 1 4×6.

图9 节点分布Fig.9.Node distribution.

采用插值型无单元Galerkin方法求解时,dmax=2.0, 此时中轴线各点挠度和右端中点挠度随时间变化的总相对误差为 0.30% , 计算时间为21.88 s.

时采用无单元Galerkin方法求解时,dmax=3.0 ,α=1.5×108, 从而可以得到与插值型无单元Galerkin方法相近的计算精度, 计算时间为28.76 s.

上述两种方法所得数值解与有限元法数值解的对比如图10和图11所示.可以看出, 插值型无单元Galerkin方法和无单元Galerkin方法都可以达到较高的精度, 但插值型无单元Galerkin方法能够提高计算效率.

算例3受均布内压的圆环

受均布内压的圆环如图12所示, 其内表面承受均匀分布的压力p=30kPa.几何参数为a=1m ,b=5m.材料体积变化是弹性的,E=1.0×107Pa ,ν=0.25, 剪切变形的流变性质采用Kelvin模型.其参数为G=5.0×105Pa ,η=2.0×106Pa·s.不计体力, 按平面应变问题计算.

图10 t=30s 时梁中轴线上的节点挠度Fig.10.Vertical displacements of nodes on the neutral axis of the beam when t=30s.

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

图12 受均布内压的厚壁圆筒Fig.12.Circular ring under a distributed inner pressure.

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

其中

(70)式表明Kelvin模型加载瞬时位移为零, 当时间足够长, 位移将趋于稳定.

如图13, 由于对称性, 只取1/4区域为研究对象.无单元Galerkin方法和插值型无单元Galerkin方法均选用如图14所示的的节点布置, 背景积分网格选取.

图13 受均布内压1/4圆筒Fig.13.A quarter of the circular ring under a distributed inner pressure.

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

采用插值型无单元Galerkin方法求解时,dmax=1.01, 此 时t=30s 时 沿x2=0 线上节点的位移和点 (2,0) 的径向位移随时间变化的总相对误差为 1.10% , 计算时间为23.38 s.

采用无单元Galerkin方法求解时,dmax=1.5 ,α=8.0×107, 此时可得到与上述插值型无单元Galerkin方法相近的相对误差, 计算时间为28.66 s.

两种方法得到的数值解与解析解的对比如图15和图16所示.从以上分析和图15和图16的对比可看出, 插值型无单元Galerkin方法和无单元Galerkin方法的数值解均与解析解符合得较好,但插值型无单元Galerkin方法能够提高计算效率.

图15 t=30s 时沿 x2=0 线上节点的位移Fig.15.Radial displacements at x2=0 when t=30s.

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

算例4工程算例: 受静水压力的混凝土水坝

为了证明本方法在较复杂几何求解域的可行性, 本算例计算一个受到静水压力的混凝土水坝,水坝的几何参数如图17所示.坝高100 m, 上游水位90 m, 下游无水.坝底部完全固定, 仅考虑坝体受到的上游水荷载.坝体材料为混凝土, 所以需要考虑建造完工后期混凝土的流变效应.假定混凝土材料体积变化是弹性的,E=2.0×1010Pa ,ν=0.2 ,混凝土的早期流变性质满足三参数模型, 其参数为G1=8.33×109Pa ,G2=4.0×1010Pa 和η=8.0×1016Pa·s.不计体力, 按平面应变问题计算.

图17 受静水压力的混凝土水坝Fig.17.A concrete dam under hydrostatic pressure.

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

图18 混凝土水坝的节点分布Fig.18.Node distribution of a concrete dam.

采用插值型无单元Galerkin方法求解时,dmax=1.2,t=500 d时沿x1=15 方向上节点的水平位移和点 (15,50) 的水平位移与时间关系的总相对误差为 0.27% , 计算时间为50.42 s.

采用无单元Galerkin方法求解时,dmax=3.0 ,α=2.0×1012, 其总相对误差为 1.36% , 计算时间为94.44 s.

两种方法的数值解与有限元法的数值解的对比如图19和图20所示.从以上分析和图19和图20的对比可以看出, 与无单元Galerkin方法相比, 插值型无单元Galerkin方法具有更高的精度和效率, 说明该方法在解决复杂工程问题时能够提高计算精度和计算效率.

图19 t=500 d时沿 x1=15 方向上节点的水平位移Fig.19.Horizontal displacements at x1=15 when t=500d.

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

6 结 论

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