光滑粒子法与有限元的耦合算法及其在冲击动力学中

第２７卷第６期

２００７年１１月爆炸与冲击ＥＸＰＬＯＳｌ０ＮＡＮＤＳＨＯＣＫＶ０１．２７．Ｎｏ．６Ｎｏｖ．．２００７ＷＡＶＥＳ

文章编号：ｌＯＯｌ一１４５５（２００７）０６一０５２２一０７

光滑粒子法与有限元的耦合算法

及其在冲击动力学中的应用。

王吉，王肖钧，卞梁

（中国科学技术大学力学与机械工程系，安徽合肥２３００２６）

摘要：扼要讨论了光滑粒子法的离散思想，充分利用光滑粒子法和有限元方法各自的优点，提出了一种

初始时刻用有限元建模，计算过程中大变形单元自动转换为光滑粒子的耦合算法。高速碰撞的系列算例说

明，耦合算法不但适宜于计算大变形冲击动力学问题，而且由于集两种方法的优点于一身，可以更高效地模拟

一些高速碰撞问题，提高计算效率。

关键词：固体力学；耦合算法；光滑粒子法；高速碰撞；有限元

中图分类号：０３４７；０４１３．１国标学科代码：１３０?１５７０文献标志码：Ａ

１引言

高速碰撞是典型的冲击力学问题，高速碰撞数值模拟计算一直是冲击工程数值计算的难点。当采用Ｌａｇｒａｎｇｅ有限元或有限差分方法计算高速碰撞问题时，必然遇到网格大畸变和滑移面处理等一系列数值计算中的关键问题，他们常常成为Ｌａｇｒａｎｇｅ有限元或有限差分法在工程应用中的制约因素。虽然Ｅｕｌｅｒ方法不存在计算大畸变问题的困难，但是难以准确描述各类界面。因此寻找适宜于高速碰撞数值模拟的新的计算方法，一直是计算力学工作者的努力目标。光滑粒子法（ｓｍｏｏｔｈｅｄ

ｎａｍｉｃｓｐａｒｔｉｃｌｅｈｙｄｒｏｄｙ—ｍｅｔｈｏｄ）便是其中的一种典型方法。光滑粒子法是一种纯Ｌａｇｒａｎｇｅ无网格算法，通过带质量的粒子离散计算域，粒子本身便代表材料，不同材料的粒子自然地构成界面，不同材料粒子的相对运动便形成所谓界面的滑移。因此从理论上说，光滑粒子法可以比较“自然地”模拟高速碰撞、侵彻贯穿等物理现象，因此是冲击力学中一种比较理想和有前途的数值方法，但其计算效率要低于有限元法。为了提高计算效率，Ｇ．Ｒ．Ｊｏｈｎｓｏｎ等［１屯］将光滑粒子法与有限元方法耦合起来，发挥他们各自的长处，在变形大的地方采用光滑粒子方法，其他区域使用有限元算法，开展侵彻贯穿方面的数值计算，取得了一些有意义的结果。

本文中在Ｇ．Ｒ．Ｊｏｈｎｓｏｎ等［１吨］的工作基础上，提出一种初始时刻用有限元方法建模、随着变形的增大、大变形区域的有限元单元自动转换成光滑粒子的耦合算法。并详细讨论耦合算法中的几个关键问题，如粒子生成的条件与方式、滑移计算与界面计算等。通过经典Ｔａｙｌｏｒ碰撞、高速和超高速碰撞等问题的计算实例，分析耦合算法的长处和特点。

２光滑粒子法

Ｊ．Ｋ．Ｃｈｅｎ等口］从Ｔａｙｌｏｒ展开式出发，引入核函数的概念，建立了光滑粒子法的基本算式。设函数，（工’），对粒子ｉ所在位置工；进行Ｔａｙｌｏｒ展开，并引入核函数ｗ；＝Ⅳ（一一藏Ｉ，＾），其中＾为光滑长度。在解域ｎ中积分，便可得到函数，（ｘ’）一阶偏导数核估计的求和形式

?收稿日期：２００６一０５—２２；修回日期：２００６－０９一０４

基金项目：国家自然科学基金项目（１０５７２１３４）作者简介：王吉（１９７９一），男，博士研究生。

万方数据　

第６期王吉等：光滑粒子法与有限元的耦合算法及其在冲击动力学中的应用５２３

怠§竺！厶二五！鲨幽一型§竺』！茎二兰！！里ｇ！￡ｌＤ』ａｚ？台ｆ１、”７

１０』

式中：Ｎ表示粒子数，由于核函数ｗ是强尖峰函数，因此实际计算时，只需考虑ｗ某一邻域里的粒子，ｍＪ是粒子ｚＪ的质量，ＷⅢ表示离散核函数Ｗｄ对空间坐标ｚ４的偏导数。（１）式是以空间坐标ｚ｝为方向的代数方程组，求解该方程组，便可获得函数厂（ｘ’）沿三个方向的一阶偏导数。可以证明由此获得的一阶偏导数的核估计，对于内点精度为二阶，而对于边界点则为一阶。采用（１）式对守恒方程中的空间导数进行离散，便可得到光滑粒子法的基本算式。

３光滑粒子法和有限元法的耦合

在高速和超高速碰撞的数值模拟计算中，Ｌａｇｒａｎｇｅ有限元常常遭遇局部区域一些网格畸变过大，以及侵彻面推进过程中两侧单元滑移计算的困难。虽然光滑粒子法可以自然地避开这些问题，但是其计算效率远低于有限元方法，因此其工程应用受到了很大限制。考虑到网格畸变过大和滑移计算只发生在局部区域的一些个别单元上，倘若能将出现计算困难的单元自动转换为光滑粒子，按光滑粒子法计算，其余区域依然用有限元方法计算，则可充分发挥两者的优势，切实解决高速碰撞中的计算难题。然而耦合算法也提出了新的课题：（１）如何实现有限元单元向光滑粒子的自动转换；（２）光滑粒子与有限元单元间的界面计算问题。

３．１有限元单元向光滑粒子的转换

Ｌａｇｒａｎｇｅ有限元在高速侵彻数值计算中，滑移面两侧单元常常是大畸变单元，因此需要依据网格具体变形情况设计光滑粒子生成的条件。在耦合算法中，我们借用有限元计算中处理大畸变单元的思想，以等效塑性应变作为单元转换的判据，即当位于滑移面两侧Ｌａｇｒａｎｇｅ单元的等效塑性应变达到某一特定值时，将此单元转换成光滑粒子。此时，光滑粒子的质量、体积和应力都等同于被置换Ｌａｇｒａｎｇｅ单元相应的物理量，而粒子的几何量，如坐标和速度则取为该Ｌａｇｒａｎｇｅ单元节点的平均值，再利用光滑粒子的速度分布，求出该时刻粒子的应变率张量和旋转率张量。光滑粒子的光滑长度＾可以近似通过质量守恒条件获得

矗＝忍ｏ（风／ｌＤ）１７”

间的维数，轴对称情况下取ｕ一３。

图１为位于滑移面上的Ｌａｇｒａｎｇｅ单元向光滑

粒子转换的示意图，虚线Ｎ为初始滑移面，实线Ｍ

为单元转换以后的滑移面。以该时刻为例，滑移面

上有４个单元满足转换条件而转换成光滑粒子（黑

色的圆）。

应该指出，有限元算法里，等效塑性应变是作为

单元因畸变过大而破坏所设立的判据，但在耦合算

法里，等效塑性应变是单元转换为粒子的阈值条件，

两者的物理意义并不相同。耦合算法里的等效塑性

应变直接关系到光滑粒子的个数。由于粒子数的多

寡对计算结果的影响并不是本文研究的重点，因此

在下面的算例里将参考Ｇ．Ｒ．Ｊｏｈｎｓｏｎ的工作‘¨，将Ｆｉｇ．１

阈值取为０．５，以说明本文方法的有效性。

３．２光滑粒子与有限元交界面的滑移计算

界面处理也是高速碰撞计算中的关键问题，在耦合算法中，同样也存在着光滑粒子与有限元界面

（滑移面）之间的接触问题。如图２所示，图中每一个圆都代表一个光滑粒子，当粒子侵入滑移面Ｍ时图１（２）式中：＾。表示初始设定的光滑长度，它的选取带有经验性，并和初始的Ｌａｇｒａｎｇｅ单元尺寸相关，ｕ为空Ｌ。ｇ，。ｎｇ。单元与光滑粒子的转换示意图ＤｅｓｃｒｉｐｔｉｏｎｏｆｔｈｅｓＰＨｎｏｄｅｇｅｎｅｒａｔｉｏｎａｌｇｏｒｉｔｈｍ万方数据　

５２４爆炸与冲击第２７卷

ｒｍｌ口ｌ＋ｍ２耽＋ｍｆ口ｆ＝ｍｌｕｌ＋ｍ２钝＋ｍ∥ｆ

ｌ优ｌ可１ｒ１＋ｍ２忱ｒ２＋ｍｆ可ｉｒ‘＝ｍｌ西ｒ１＋ｍ２谠ｒ２＋ｍ∥：ｒｉ

ｌｕ：一西＝（ｕ：一西）（ｒ｛一，．。）／（ｒ。一ｒ。）（３）

Ａｒｅａｏｆｓｍｏｏｔｈｅｄｐａｒｔｉｃｌｅｓ

ＳｌｉｄｉｎｇｉｎｔｅｒｆＡｃｅ

ＡｒｅａｏｆＬａｇｒａｎｇｅｆｊｎｉｔｅｅｌｅｍｅｍｓ

图２Ｌａｇｒａｎｇｅ单元与光滑粒子的界面模型

Ｆｉｇ．２图３交界面附近粒子的计算示意图Ｆｉｇ．３ＣｏｍｐｕｔａｔｉｏｎｓｏｆｔｈｅＳＰＨｎｏｄｅｓ

ｎｅａｒＳＰＨｎｏｄｅｓｌｉｄｉｎｇｏｎａｓｔａｎｄａｒｄｇｒｉｄｔｈｅｓｌｉｄｉｎｇｉｎｔｅｒｆａｃｅ

４计算实例和讨论

４．１Ｔａｙｌｏｒ碰撞

Ｔａｙｌｏｒ碰撞是典型的高速碰撞问题。１９４８年，Ｇ．Ｉ．Ｔａｙｌｏｒ研究了一平头柱形杆弹法向撞击刚性靶板问题。从一维弹塑性波理论出发，采用刚塑性本构模型，解析获得了杆弹撞击后柱体变形与材料动态屈服强度间的近似关系。由于Ｔａｙｌｏｒ碰撞体现了材料大变形高应变率等典型的冲击力学特性，因此通过Ｔａｙｌｏｒ碰撞研究材料力学性能一直受到关注。ｗ．Ｈ．Ｇｕｓｔ［４］和Ｇ．Ｒ．Ｊｏｈｎｓｏｎ等［２３都先后通过数值方法研究了Ｔａｙｌｏｒ碰撞问题，通过与实验结果的比较，验证本构模型和计算方法的合理性和准确性。为此我们采用光滑粒子法及耦合算法对Ｔａｙｌｏｒ碰撞开展数值模拟，以便与已有的理论研究和实验资料相比较，确定本文所述方法和所编程序的可靠性。

设直径Ｄｏ＝７．６ｍｍ、长度Ｌｏ＝２５．４ｍｍ的钢柱，以速度矾一２２１ｍ／ｓ法向撞击刚性平板，材料参数分别为［５］：密度胁一７．８９０Ｔ／ｍ３，冲击压缩实验中测得的Ｄ“线性关系中的常数ｃ。一３．６ｋｍ／ｓ，ｓ＝１。８０，剪切模量户＝８０ＧＰａ，简单拉伸条件下的屈服应力ｙ０＝ｏ．５ＧＰａ，Ｇｒｕｎｅｉｓｅｎ系数ｒ＝１．８１。设材万方数据　

第６期王吉等：光滑粒子法与有限元的耦合算法及其在冲击动力学中的应用５２５料为理想弹塑性体。图４给出了碰撞后的变形示意图，定义变形后柱体的长度为Ｌ，撞击端直径为Ｄ，离固壁ｏ．２Ｌ０处的直径为Ｗ。张刚明等‘６３给出了同样算例的实验结果，Ｌ＝１９．８ｍｍ，Ｄ一１３．７ｍｍ，Ｗ

、

＝８．６ｍｍ。

分别采用纯光滑粒子法和耦合算法对Ｔａｙｌｏｒ碰撞开展数值模

拟计算。初始时刻，耦合算法以三角形单元建模，共６５８个单元，单

元数与纯光滑粒子法模型的６３０个粒子相近，计算过程中，随着有限

元单元变形的发展，部分单元自动按上述方法转换成光滑粒子。

表１计算和实验结果的比较

ＴａｂＩｅ１Ｃ伽叩ｎ幅ｉｏｎｂｅｔｗｅｅｎｃＯｍｐｕｔａｔｉＯ璐ａｎｄｅｘｐｅｒｉｍｅｎｔ

图４ＴａｙＩｏｒ碰撞示意图

Ｆｉｇ．４Ｔａｙｌｏｒｉｍｐａｃｔ

图５～６分别给出了纯光滑粒子法和耦合算法在ｌｏｐｓ和４０ｐｓ时的变形图，耦合算法中杆弹上部仍为有限元，下部已转换为光滑粒子。表１给出了计算值与实验结果的比较。虽然几个特征尺寸耦合算法要比光滑粒子法更接近实验值，但是从撞击端的总体构形上看，两者仍有差异，这说明耦合算法中粒子的生成和界面算法对结果有一定影响。通过与实验结果的比较，进一步改进相关算法，将是本文后续工作的重点。然而Ｔａｙｌｏｒ碰撞的初步计算结果表明，采用本文给出的思想建立耦合算法是成功的。

图５１０肛ｓ时变形图

Ｆｉｇ．５Ｄｅｆｏｒｍａｔｉｏｎａｔｚ—ｌＯｐｓ

图６４０ｐｓ时变形图

Ｆｉｇ．６Ｄｅｆｏｒｍａｔｉｏｎａｔ‘一４０

ｐｓ

万　方数据

爆炸与冲击第２７卷

４．２有限元、光滑粒子法和耦合算法的对｝匕计算

钨合金平头杆弹，弹长２０ｍｍ、直径１０ｍｍ，以初速‰＝１．毒ｋｍ／ｓ正撞击厚为１０ｍｍ的钢靶，分别采用有限元法、光滑粒子法和耦合算法进行对比计算，靶板和杆弹视为理想弹塑性体。表２给出了计算中用到的两种材料的参数。

表２两种材料的参数

ＴａｂＩｅ２Ｐａ甩ｍｅｔｅ伟０ｆｔｗ０ｍａｔｅｒｉａｌｓ

采用有限元程序计算时，单元取为三角形单元，共２０４８个单元。采用光滑粒子法计算时，共１０２４个粒子。而在耦合计算中，初始模型只取了８００个三角形单元。

图７～８分别给出了三种方法在不同时刻（３弘ｓ和１０肛ｓ）的计算结果。对于畸变过大的单元，有限元方法用消蚀法处理，即当单元等效塑性应变超过某一阈值时，该单元被认为破坏，只能承受静水压力。计算中钨合金和钢靶的极限等效应变分别取为１．ｏ和ｏ．９。光滑粒子法由于不存在单元，等效塑性应变增大，相当于粒子距离减小，因此不会出现有限元方法中因单元变形太大导致计算无法进行下去的情形。而在耦合算法中，由于只有部分单元转换成粒子（如３ｐｓ时为９８个粒子，１０弘ｓ时为３１０个粒子），因此既能有效处理畸变过大的计算问题，又保留了有限元计算效率高的优势。从图中可以看到，三种方法的变形图像相差不是很大。由于在有限元方法中，破坏的单元对子弹不再存在阻力作用，因此相同时刻杆弹的速度要比光滑粒子法计算的结果大。有限元方法中子弹前部速度比较小的单元由于达到了破坏条件，成为破坏单元，在图上没有画出。光滑粒子法中，由于子弹前部的变形较大，形成的蘑菇头要比有限元方法的大，从而开坑半径要大一些。而耦合算法计算结果更趋近于光滑粒子法。

图７３肚ｓ时变形图

Ｆｉｇ．７Ｄｅｆｏｒｍａｔｉｏｎａｔ￡＝３ｐｓ

图８１０ｐｓ时变形图

Ｆｉｇ．８Ｄｅｆｏｒｍａｔｉｏｎａｔｔ＝ｌＯ肛ｓ

４．３超高速碰撞计算

考虑高速平头弹丸侵彻半无限厚靶板的情况，子弹长２０ｍｍ、弹径１０ｍｍ，初始弹速４ｋｍ／ｓ，子弹和靶板的材料与算例２相同。光滑粒子法共用了１０４９６个粒子。耦合算法中共用了８２００个单元（粒

子）。图９～１１分别给出了两种算法在几个不同时刻的变形图。图上颜色淡一点的是子弹的粒子，而颜万方数据