1.2 无网格法的研究进展
1977年,Lucy[1]和Gingold 等[2]分别提出光滑质点流体动力学方法(Smoothed particle hydrodynamics,SPH),并成功应用于求解天体和流体问题,成为最早的无网格方法。SPH 的缺点在于对于非均匀分布的节点不满足一致性条件,对于均匀分布的节点在边界上不满足一致性条件,从而导致计算精度的下降。再加之当时FEM 正处于巨大成功的阶段,该类方法并未受到学术界的关注,发展相对缓慢。直到1992年,Nayroles 等[3]基于移动最小二乘法(Moving least square method,MLS),提出了散射元法(Diffuse element method,DEM),才真正拉开了无网格法研究的序幕。DEM 存在计算精度不高,边界条件不易施加等缺点。1994年,美国西北大学的著名学者Belytschko 等[4]对DEM 进行了细致的改进,提出了著名的无网格伽辽金法(Element free Galerkin method,EFGM),在计算力学界掀起了研究无网格法的热潮。较SPH,EFGM 的计算费用较高,但计算精度和稳定性更好。因此,该方法被广泛应用于裂纹扩展分析[5-11]、板壳分析[12]、三维撞击和流体晃动分析[13-14]、相变分析[15]、节理岩体变形分析[16]、扩散问题[17]、拱坝开裂[18]、边坡开挖[19]、固结问题[20]以及弹性地基板[21]等领域。
1995年,Liu 等[22,23]利用再生核函数思想和函数积分变换法,基于Galerkin 弱形式提出了再生核质点法(Reproducing kernel particle method,RKPM)。该方法中的近似函数由窗函数(Windows function)和傅里叶变换建立。由于窗函数具有平移、缩放等优良特性,非常适合应用于弹性、弹塑性和动力学等问题。之后,Liu 等[24]又借鉴小波函数理论中的多尺度分析技术,提出了多尺度再生核质点法(Multiple scale reproducing kernel particle method,MRKPM)。RKPM 及其改进方法在许多问题的分析中得到了很好的应用,如延性断裂、声学分析、动力学分析、中厚度板、金属加工成形等。1995年,Duarte 等[25]利用MLS建立单位分解近似函数,并由此构造形函数和权函数,系统方程由Galerkin 弱形式建立,提出了hp 云团法(hp-clouds meshless method,hp)。目前,该方法已被应用于求解悬臂梁问题、厚板弯曲以及二维裂纹扩展的自适应分析等问题。1998年Oden等[26,27]选取有限元插值函数作为单位分解函数,提出了基于云团法的新型hp 有限元(Clouds based on hp FEM,hp-FEM)。Liszka 等[28]采用配点法取代Galerkin 法,提出hp 无网格云团法(hp meshless clouds method)。
1996年,Onate 等[29]采用MLS 构造近似函数,并用配点格式离散系统方程,提出有限点法(Finite point method,FPM)。该方法无需使用背景网格,效率较高,目前主要应用于流体力学问题分析[30]和弹性力学问题[31]。同年,Babuska 等[32]提出了单位分解法(Partition of unity method,PUM),并对其进行了严格的数学证明;同时还指出可以使用特征解析解来增强基函数。Menlenk等[33]和Babuska 等[34]又先后将PUM 与有限元方法结合,提出了单位分解有限元(Partition unity of unity finite element method,PUFEM) 和广义有限元(Generalized finite element method,GFEM)。这些方法由于网格划分独立于结构几何形态,因此在处理强、弱不连续问题时显示了良好的性能,可以处理任意形状裂纹。
1998年,Zhu 等[35]和Atluri 等[36]分别采用MLS 构造试函数和加权残值法中的权函数,用局部Petrove-Galerkin 法建立系统离散方程,提出了局部边界积分方程方法(Local boundary integral equation method,LBIE) 和无网格局部Petrove 伽辽金法(Meshless local petrove Galerkin methods,MLPG)。实际上,LBIE是MLPG 的一种特殊形式,需要奇异积分。后者较前者具有更高的计算效率和精度。两种方法均不需要插值网格和背景积分网格,积分在局部子域或局部边界完成,是真正意义上的无网格法。刘桂荣等[37]将MLPG 与有限元和边界元方法耦合,以发挥它们各自的优势。目前,该方法已被用于求解弹性力学问题、二维位势问题、梁弯曲和裂纹扩展等问题。
2000—2001年,清华大学的张雄教授先后提出了配点型无网格法[38]、最小二乘配点无网格法[39]等新方法。该类方法解决了配点型无网格法精度低、稳定性差的问题,较Galerkin 法计算量要少很多。2001年,张建明、姚振汉[40]结合MLS 和修正变分原理,提出了杂交边界点法(Hybrid boundary node method,HBNM)。HBNM 是一种仅需在边界布点的无网格法,计算精度高,收敛性好,已被成功应用于求解弹性力学问题。同年,陈文等[41,42]提出边界节点法(Boundary knot method,BKM)和边界粒子法(Boundary particle method,BPM)。BKM 和BPM 是径向基函数和边界元法结合的产物,只需在问题域的边界上离散节点,而不需要离散整个计算区域。也是在同年,刘桂荣等[43]提出了无网格点插值法(Point interpolation method,PIM)。PIM 采用点插值法建立无网格形函数,具有Kronecker-delta 函数特性,可以像有限元那样方便地施加本质边界条件。但由于插值函数受节点分布的影响,插值函数的稳定性比较差。为了改善PIM 的性能,2002年刘桂荣等[44-48]又提出和发展了径向基点插值法(Radial basis point interpolation method,RPIM)和局部径向基点插值无网格法(Local radial point interpolation method,LRPIM)。RPIM和LPIM 都是在PIM 的多项式基函数中加入径向基函数,大大改善了插值函数不稳定的缺点,而且还保留了PIM 的所有优点。虽然其形函数破坏了全局的相容性条件,但较PIM 具有更好的收敛性和精度。目前RPIM 和LRPIM 已被应用于二维弹性力学、弹塑性问题、智能材料、桥墩和地基等问题中。
Venture 等[49]受XFEM 的启发,提出了单位分解扩展无网格法。当扩展发生时,通过向量水平集方法确定裂纹路径。该方法为裂纹的产生、分岔和多裂纹连接等问题提供了非常有效的处理手段,不足之处在于水平集更新算法十分复杂。粒子断裂法[50]也是基于阶跃函数的,在裂纹线上离散节点,这些点的集合即可描述裂纹。在该方法的基础上,Rabczuk 等[51]提出扩展无网格伽辽金法。
程玉民、李九红等[52,53]将复变量MLS 和边界无网格法结合,提出求解弹性力学问题的复变量边界无网格法。该方法的形函数具有Kronecker-delta 性质,很容易施加本质边界条件,提高了精度和计算效率。目前,该方法主要应用于断裂问题[54]和大变形问题[56]。
张雄等[56]以冲击爆炸问题为背景,对物质点法进行了系统研究,研发了三维物质点法仿真软件MPM3D 软件,可用于模拟超高速碰撞、冲击侵切和爆炸等强冲击荷载作用下材料与结构的力学行为。物质点法采用拉格朗日和欧拉双重描述,将物体离散为一组在空间网格中运动的质点。质点携带了所有的物质信息,如质量、速度、应力和应变等,很容易跟踪材料界面和引入与变形历史有关的材料模型。质点在空间网格中运动,运动方程在空间网格上求解,避免了网格畸变,非常适合分析特大变形和流动问题。
考虑配点型和局部弱式无网格法的优、缺点,刘桂荣和顾元通[57]提出了无网格强弱式(Meshless weak-strong form,MWS)。MWS 的基本思想是如果某个节点位于问题域内或本质边界上,采用配点法建立系统方程;如果某个节点位于或非常接近导数边界上则采用局部弱式建立系统方程。然后利用无网格形函数,从而得到离散的系统方程。目前,该方法已被应用于二维弹性力学问题和不可压缩流体问题,取得了令人满意的计算结果。MWS 不但继承了Galerkin 法计算精度高、稳定性好的特点,还大大提高了计算效率。
为了更好地完善EFGM,许多研究工作者致力于对其性质和改进方面的研究。Belytschko 等[58,59]给出了EFGM 的误差估计,对其数值积分方案进行了深入研究,奠定了该方法的数学基础。刘桂荣等[60]将EFGM 与边界元耦合,并应用于固体力学问题。Belytschko 等[61]和Hegem[62]将EFGM 与有限元耦合,可以很好地发挥它们各自的优势。Beissel 等[63]为了避免在EFGM 中使用背景积分网格,提出节点积分方案,但计算稳定性较差。Chen 等[64,65]利用应变光滑技术,根据线性分片实验推导出Galerkin 无网格法对应的积分约束,提出了稳定一致节点积分(Stabilized conforming nodal integration,SCNI)。SCNI 很好地继承了节点积分的优势,同时完全避免了形成整体刚度矩阵时复杂形函数导数计算。另外,SCNI 能保证线性精度,而与无网格离散模式无关。为了弱化场函数近似的一致性要求,刘桂荣及其团队[66,67]推广了应变光滑技术,提出了G 空间理论和广义光滑伽辽金(GS-Galerkin)弱形式,建立了两类光滑数值方法(weakened weak,W2)的理论基础。这两类方法分别称为光滑有限元法(smoothed finite element method,S-FEM) 和光滑点插值法(smoothed point interpolation method,S-PIM)。S-FEM[68-71]的基本思想是:在有限元网格上创建的光滑区域内光滑化场函数导数,进而获得比传统有限元“更软”的光滑刚度矩阵,提高解的精度和收敛性。由于不需要计算形函数导数,避免了等参映射过程。因此,S-FEM 对畸形网格具有极强的适应性。S-PIM 与SFEM 的基本原理是相同的,不同之处在于S-PIM 使用多项式或径向基函数建立节点形函数。刘桂荣[72]指出,S-FEM 仅仅是S-PIM 的特殊形式。Cui 等[73]将MLS 与GS-Galerkin 弱形式结合,提出了光滑伽辽金法。但与Galerkin 无网格相比,其计算精度较低。段庆林等[74-76]将SCNI 中的线性积分约束推广到二阶情况,提出了二阶一致节点积分(Quadratially consistent nodal integration,QCI)。QCI 满足二阶精度,但为了求解形函数导数,必须在每个光滑区域内额外求解2 个3×3 的代数方程,计算代价并不小。王东东等[77]将三角形背景网格划分为两水平嵌套光滑区域,利用梯度光滑技术和Richardson 外推法,提出了嵌套子域梯度光滑积分(NSGSI)算法。NSGSI 算法将SCNI 的精度提高到4 阶。与高斯积分相比,NSGSI 仅需要6 个积分点,极大地提高了计算效率。为了得到二阶精度的数值算法,NSGSI 必须使边界积分、外力项积分与刚度积分保持一致。显然,NSGSI 的一致性要求大大降低了其灵活性。马文涛[78]将Richardson 外推法与广义光滑梯度技术相结合,提出了嵌套光滑无网格法(NSMM)。NSMM 不仅保持了无网格法高精度的特点,还极大地提高了计算效率。