论文导读：：基于无网格自然单元法，提出了结构动力弹塑性响应分析的一条新途径。自然单元法是一种新兴的无网格数值计算方法，其实质是基于自然邻近插值的伽辽金法。自然单元法在本质边界条件的施加上较采用移动最小二乘法的无网格法具有明显的优势。在空间域上采用自然单元法离散，并运用加权余量法推导了动力弹塑性分析的离散控制方程。然后，采用预校正形式的Newmark法在时间域上进行求解。最后给出了数值算例，并验证了所提方法的有效性和正确性。

　　论文关键词：无网格法，自然单元法，动力响应，弹塑性

　　0引言

　　在结构的动力响应过程中，通常总是既有弹性变形，又有塑性变形，而且这两种变形以及它们之间的分界面都随时间变化[1].正是由于弹塑性动力响应问题的复杂性，往往得不到结构弹塑性动力响应问题的解析解，需要借助数值方法进行求解。 在结构弹塑性动力响应的模拟方面，不论是有限元法[2]还是边界元法[3, 4]，都是基于网格的数值方法。不合适的网格形状严重影响计算的精度，而且复杂结构的网格生成也是极具挑战性的问题。 相反，近年来发展迅速的无网格法具有近似函数不依赖于网格，且计算精度高，前后处理简单的优点[5].

　　目前发展的无网格方法主要有无单元Galerkin法[6, 7]、无网格局部Petrov-Galerkin法[ 8]、边界无单元法[9, 10]、杂交边界点法[11, 12]、无网格流形方法[13]以及自然单元法[14, 15]等。 与大多数无网格法不同，自然单元法的近似函数具有插值性且在边界结点间是线性变化的无网格法，从而可以方便地施加本质边界条件。 此外，自然单元法的形函数计算不仅具有不涉及矩阵求逆运算以及计算量较小的优点，而且也没有任何人为参数的选择问题。 因此，自然单元法是一种非常有发展前景的数值方法，并且已引起了国内外许多学者[16~22]的极大关注。

　　近年来，许多学者都致力于弹性动力学的无网格法研究。 程玉民等[10]将弹性动力学的边界积分方程方法与改进的移动最小二乘法结合，提出了弹性动力学的边界无单元法。 苗雨等[12]将双互易法同杂交边界点法相结合，提出了求解弹性动力学问题的双互易杂交边界点法。 李树忱等[13]运用无网格流形方法求解了弹性动力学问题。 刘应华等[23]发展了弹性动力学分析的基于Voronoi结构的无网格局部Petrov-Galerkin法。 虽然朱合华等[19]，张英新等[20]和江涛等[21]在无网格自然单元法应用于弹塑性力学问题方面做了一些有益的工作，但目前尚未见到动力弹塑性分析中无网格自然单元法的研究成果。

　　本文尝试将自然单元法应用于动力弹塑性分析问题的求解计算。 首先采用加权残值法详细推导了动力弹塑性分析的自然单元法理论公式，然后给出了其详细的数值实现过程。 最后，通过典型算例的计算和对比分析验证了方法的有效性和合理性。

　　1 控制方程

　　考虑二维动力弹塑性分析问题，其计算域为，边界为，则平衡方程为：

　　弹塑性（1）

　　式中，为质量密度；为结点加速度，且有，为结点位移；为应力；为应变；为体力矢量；和分别为位移和面力已知的边界。

　　应力和位移满足如下的边界条件

　　（2a）

　　弹塑性（2b）

　　和初始条件

　　弹塑性（3a）

　　（3b）

　　式中，为上点处的外法线的方向余弦；和分别为已知的位移和面力分量；和分别为初始位移和初始速度。

　　弹塑性材料在进入塑性状态后，继续加载时的应力应变关系可以表示为

　　（4）

　　式中，为总应变增量。 采用关联的流动法则以及各向同性的加工强化模型，可得到弹塑性刚度张量为[ 2]

　　（5）

　　其中

　　（6a）

　　（6b）

　　（6c）

　　这里，和分别为等效应力和等效塑性应变，为单向屈服应力，而和分别为剪切模量和泊松比。

　　2 动力弹塑性分析的无网格自然单元法

　　2.1 自然邻近插值[14~23]

　　自然邻近插值是一种多变量的插值方案，目前应用于无网格法的有Sibson插值和Laplace插值。本文采用Sibson插值。

　　对区域内任一结点，其Voronoi结构定义为：

　　（7）

　　式中，为点与结点的距离。此外无网格法，二次Voronoi结构可定义为：

　　（8）

　　图1所示为平面7个结点的Voronoi结构和待插值点的二次Voronoi结构。

　　对于Sibson插值，计算点对结点的形函数定义为与结点的二次Voronoi结构的面积与的一次Voronoi结构的面积之比，即

　　（9）

　　定义了各结点的插值函数后，点的位移函数类似于有限元法可写为

　　（10）

　　式中，是点周围自然邻结点的结点位移，为对应结点的形函数。

　　2.2 控制方程的弱形式及其离散化

　　平衡方程式（1）及力的边界条件式（2b）的等效积分形式的伽辽金提法可以表示为

　　（11）

　　对式（11）进行分部积分，则有

　　（12）

　　为便于进行数值计算，我们把上式改写成矩阵形式，有

　　（13）

　　式中

　　（14）

　　由于只对空间域进行离散，求解域内的试函数可由式（10）表示为

　　（15）

　　将空间离散后的位移表达式（15）代入式（14），并注意到结点位移变化的任意性，最终得到系统的运动方程如下：

　　（16）

　　式中，为结点的加速度向量；，和分别为质量矩阵、结点内力向量和结点外力向量，且它们各元素可具体表示为：

　　（17a）

　　（17b）

　　（17c）

　　其中

　　, （18）

　　3 时间积分方案

　　本文对时间域的离散采用应用较为广泛的Newmark方法[2].Newmark 方法是一种隐式算法，因此时刻的运动方程应得到满足。此外，由于塑性变形的非线性特性，在每个时间步都必须进行迭代计算。应用Newton-Raphson迭代，时刻的运动方程可以改写为

　　（19a）

　　（19b）

　　式中，切线刚度矩阵可以表示为

　　（20）

　　在和时间步内，Newmark方法采用的位移和速度的递归关系为

　　（21a）

　　（21b）

　　其中和是按积分精度和稳定性要求决定的参数。 由式（19b）和式（21b）可得

　　（22）

　　将式（22）代人式（19a），则得到如下的静力等效问题

　　（23）

　　式中

　　（24a）

　　（24b）

　　顺便指出，给定初始的位移和初始速度后，我们就可以从下式求出初始加速度：

　　（25）

　　下面给出时刻静力等效问题（23）的具体迭代求解过程：

　　（1）令迭代计算变量。

　　（2）在开始预估阶段，令

　　（26a）

　　（26b）

　　（26c）

　　式中

　　（27）

　　（3）利用下述方程计算残余力

　　（28）

　　（4）如有需要，应用下式形成等效刚度矩阵

　　（29）

　　否则，就应用前面已计算出的。

　　（5）求解方程

　　（30）

　　得到结点位移增量，然后利用几何关系可计算得到应变增量。在此基础上，利用第1节中所述的弹塑性本构关系以及切向预测径向返回方法[24]就可以确定应力。

　　（6）进入修正阶段无网格法，令

　　（31a）

　　（31b）

　　（31c）

　　（7）如果不满足收敛条件，令并转到第3步，否则继续进行下去。

　　（8）为下一个时间步使用，令

　　（32a）

　　（32b）

　　（32c）

　　4 数值算例

　　为了验证以上提出的数值方法的有效性，本文对一些典型算例进行了动力弹塑性分析的计算和比较。 如果不作特别说明，积分时取Delaunay三角形内的3个高斯点。 为方便起见，本文算例所有参数采用无量纲形式。

　　算例1：受突加载荷的悬臂梁

　　一个长，高度的悬壁梁，端部受到剪切荷载的作用，如图2所示。材料为理想弹塑性，并服从von Mises屈服

　　准则。 材料的弹性模量，泊松比，屈服应力，质量密度。

　　计算分析中采用如图3所示的两种结点布置方案，且时间步长取为。 为了检验方法的有效性，本文使用Owatsiriwong和Park[4]的边界元法计算结果及有限元软件ABAQUS的计算结果进行对比。图4给出了悬臂梁右端中点A的竖向位移随时间的变化曲线，从图上可以看出本文采用规则结点和不规则结点的计算结果与边界元法的计算结果以及ABAQUS的计算结果都吻合很好，从而表明了本文方法的有效性，同时也表明了本文计算结果对结点的不规则分布不敏感。

　　算例2：受瞬态载荷的带孔方板

　　如图5所示，长为36，宽为20的矩形板，中央有一直径为10的圆形孔洞，在平面应力情况下两侧作用瞬态载荷，其中的表达式为

　　材料采用双线性各向同性强化本构关系，并服从von Mises屈服准则。 材料参数为：弹性模量，泊松比，初始屈服应力无网格法，切线模量，质量密度。由于对称性，取结构的四分之一作为计算模型，其结点布置方案如图6所示，共589个结点。 采用时间步长计算得到了点A处的竖向位移随时间的变化曲线，如图7所示。 为了进行对比，图7中还给出了有限元软件ABAQUS的计算结果。 从图7可以看出，本文的计算结果与ABAQUS的计算结果吻合很好，从而进一步验证了本文方法的有效性。

　　5 结论

　　自然单元法作为一种无网格数值方法，不需要对求解域进行网格划分，前处理简单。从形函数和近似函数的构造过程来看，自然单元法更像是介于有限元法与无网格法之间的一种数值方法，兼具无网格的特性和有限元的优点。本文建立了进行结构弹塑性动力响应分析的无网格自然单元法，弹塑性动力学积分弱形式的推导采用加权残数法，空间离散采用无网格自然单元法，时间离散采用预校正形式的Newmark方法。 数值算例表明，本文提出的方法能够有效用于结构的弹塑性动力响应分析问题的求解。

　　参考文献

　　[1]杨桂通。弹塑性动力学基础[M]. 北京：科学出版社， 2008.（Yang Guitong. The basis of elastoplastic dynamics[M]. Beijing: Science Press,2008 （in Chinese）。）

　　[2]Zienkiewicz O C,Taylor R L. The Finite Element Method [M]. 5th ed. Singapore: Elsevier Pte Itd,2005.

　　[3]Hatzigeorgiou GD, Beskos D E. Dynamic elastoplastic analysis of 3-D structures by thedomain/boundary element method [J]. Computers & Structures, 2002, 80:339~347.

　　[4]Owatsiriwong A,Park K H. A BEM formulation for transient dynamic elastoplastic analysis viaparticular integrals [J]. International Journal of Solids and Structures, 2008,45: 2561~2582.

　　[5]张雄，刘岩，马上。无网格法的理论及其应用[J]. 力学进展， 2009, 39（1）： 1~36. （Zhang Xiong, Liu Yan, Ma Shang. Meshfreemethods and their applications [J]. Advances in Mechanics, 2009, 39（1）： 1~36（in Chinese）。）

　　[6]Belytschko T, LuY Y, Gu L. Element-free Galerkin method [J]. International Journal for NumericalMethods in Engineering, 1994, 37: 229~256.

　　[7]王东东，李凌，张灿辉。稳定节点积分伽辽金无网格法的应力计算方法研究[J]. 固体力学学报， 2009,30（6）： 586~591. （Wang Dongdong, Li Ling, Zhang Canhui. On stress evaluation inGalerkin meshfree methods with stabilized conforming nodal integration [J].Chinese Journal of Solid Mechanics, 2009, 30（6）： 586~591（in Chinese）。）

　　[8]Atluri S N, ZhuT L. A new meshless local Petrov-Galerkin（MLPG） approach in computationalmechanics [J]. Computational Mechanics, 1998, 22（2）： 117~127.

　　[9]程玉民，陈美娟。弹性力学的一种边界无单元法[J]. 力学学报， 2003,35（2）： 181~186. （Cheng Yumin, Chen Meijun. A boundary element-free method forlinear elasticity [J]. Acta Mechanica Sinica, 2003, 35（2）： 181~186 （in Chinese）。）

　　[10]程玉民，彭妙娟。弹性动力学的边界无单元法[J]. 中国科学， 2005,35（4）： 435~448. （Cheng Yumin, Peng Miaojuan. Boundary element-free method forelastodynamics [J]. Science in China, 2005, 35（4）：435~448 （in Chinese）。）

　　[11]Zhang J M, Yao Z H, Masatakz Tanaka. The meshless regular hybrid boundary node method for 2D linearelasticity [J]. Engineering Analysis with Boundary Elements, 2003, 27: 259~262.

　　[12]苗雨，王元汉。弹性动力学的双互易杂交边界点法[J]. 力学学报， 2009,41（1）： 122~128. （Miao Yu, Wang Yuanhan. Dual reciprocity hybrid boundary nodemethod for elastodynamics [J]. Acta Mechanica Sinica, 2009, 41（1）： 122~128 （inChinese）。）

　　[13]李树忱，李术才，朱维申。弹性动力学的无网格流形方法[J]. 岩石力学与工程学报， 2006,25（1）： 141~145 （Li Shuchen, Li Shucai, Zhu Weishen. Study on meshless manifoldmethod in elastodynamics [J]. Chinese Journal of Rock Mechanics and Engineering,2006, 25（1）： 141~145 （in Chinese）。）

　　[14]Braun J,Sambridge M. A numerical method for solving partial differential equations onhighly irregular evolving grids [J]. Nature, 1995, 376: 655~660.

　　[15]Sukumar N,Moran B, Belytschko T. The natural elements method in solid mechanics [J]. InternationalJournal for Numerical Methods in Engineering, 1998, 43: 839~887.

　　[16]Calvo B, Martinez MA, Doblare M. On solving large strain hyperelastic problems with the naturalelement method [J]. International Journal for Numerical Methods in Engineering,2005, 62: 159~185.

　　[17]Bueche D,Sukumar N, Moran B. Dispersive properties of the natural element method [J]. ComputationalMechanics, 2000, 25: 207~219.

　　[18]江涛，章青。中厚板弯曲问题的自然单元法[J]. 固体力学学报， 2009,30（4）： 424~431. （Jiang Tao, Zhang Qing. Natural element method for the bendingproblem of moderately thick plate [J]. Chinese Journal of Solid Mechanics, 2009,30（4）： 424~431 （in Chinese）。）

　　[19]朱合华，杨宝红，蔡永昌，徐斌。无网格自然单元法在弹塑性分析中的应用[J]. 岩土力学， 2004,25（4）： 671~674. （Zhu Hehua, Yang Baohong, Cai Yongchang, Xu Bin. Application ofmeshless natural element method to elastoplastic analysis [J].Rock and Soil Mechanics, 2004, 25（4）： 671~674 （inChinese）。）

　　[20]张英新，王建华，高绍武。二维弹塑性自然单元法算法实现[J]. 上海交通大学学报， 2005,39（5）： 727~730. （Zhang Yingxin, Wang Jianhua, Gao shaowu. The computationalmethods of natural element method in two-dimensional elasto-plastic analysis[J]. Journal of Shanghai Jiaotong University, 2005, 39（5）： 727~730 （inChinese）。）

　　[21]江涛，章青。线性强化材料弹塑性分析的自然单元法[J]. 力学季刊， 2010,31（2）： 288~296. （Jiang Tao, Zhang Qing. Natural element method for elasto-plasticanalysis of linear strain hardening material [J]. ChineseQuarterly of Mechanics, 2010, 31（2）： 288~296 （in Chinese）。）

　　[22]王兆清，冯伟。自然单元法研究进展[J]. 力学进展， 2004,34（4）： 437~445. （Wang Zhaoqing, Feng Wei. Advances in natural element method [J]. Advances in Mechanics, 2004, 34（4）： 437~445 （in Chinese）。）

　　[23]Liu Y H, ChenS S, Li J, Cen Z Z. A meshless local natural neighbour interpolation methodapplied to structural dynamic analysis [J]. Computer Modeling in Engineering& Sciences, 2008, 31（3）： 145~156.

　　[24]Owen D J R, HintonE. Finite Element in Plasticity: Theory and Practice [M]. Swansea, Pineridge Press, UK, 1980.