一种光滑粒子流体动力学-有限元法转换算法及其在冲击动力学中的应用

第４５卷第１期２０１１年１月＝。－——

西安交通大

学学

报

Ｖ０１．４５Ｎｏ．１

ＪＯＵＲＮＡＬＯＦＸＩ’ＡＮＪＩＡＯＴ（）ＮＧＵＮＩＶＥＲＳＩＴＹ

－－——－．，——＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

Ｊａｎ．２０１１

一种光滑粒子流体动力学一有限元法转换算法

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

张志春，强洪夫，高巍然

（第二炮兵工程学院２０１室，７１００２５，西安）

摘要：为了解决冲击动力学中的大变形问题，提出了一种光滑粒子流体动力学一有限元法（ＳＰＨ—ＦＥＭ）转换算法，以等效Ｍｉｓｅｓ应力作为转换判据，将冲击过程中局部大变形区域的有限元网格转换为ＳＰＨ粒子．该算法在大变形区域使用具有优势的ＳＰＨ，在小变形区域使用精度和效率更高的ＦＥＭ，为冲击动力学问题的数值计算提供了一条有效途径．使用ＳＰＨ—ＦＥＭ转换算法对圆柱形钢弹正冲击钢板发生冲塞破坏的过程进行了三维数值计算，计算结果与实验吻合较好，显示了该算法在计算精度方面的优势．在实际工程中，需要根据具体材料的失效模式，选择更加合适的转换判据．关键词：光滑粒子流体动力学；有限元法；转换算法；冲击

中图分类号：０３４７．３文献标志码：Ａ文章编号：０２５３—９８７Ｘ（２０１１）０１—０１０５—０６

Ｃｏｎｖｅｒｓｉｏｎｏｆ３ＤＤｉｓｔｏｒｔｅｄＦｉｎｉｔｅＥｌｅｍｅｎｔｓｉｎｔｏＳＰＨＰａｒｔｉｃｌｅｓ

ＤｕｒｉｎｇＩｍｐａｃｔＤｙｎａｍｉｃＤｅｆｏｒｍａｔｉｏｎ

ＺＨＡＮＧＺｈｉｃｈｕｎ，ＱＩＡＮＧＨｏｎｇｆｕ，ＧＡＯＷｅｉｒａｎ

（ＳｔａｆｆＲｏｏｍ２０１，Ｘｉ’ａｎＨｉ—ＴｅｃｈＩｎｓｔｉｔｕｔｅ，Ｘｉ’ａＤ７１００２５，Ｃｈｉｎａ）

Ａｂｓｔｒａｃｔ：Ａｎａｌｇｏｒｉｔｈｍｔｏａｕｔｏｍａｔｉｃａｌｌｙｃｏｎｖｅｒｔ３ＤｄｉｓｔｏｒｔｅｄｆｉｎｉｔｅｅｌｅｍｅｎｔｓｉｎｔｏＳＰＨ（ｓｍｏｏｔｈｅｄ

ｐａｒｔｉｃｌｅｈｙｄｒｏｄｙｎａｍｉｃｓ）ｐａｒｔｉｃｌｅｓｉｓｐｒｏｐｏｓｅｄ，ｗｈｉｃｈｈａｎｄｌｅｓｔｈｅｐｒｏｂｌｅｍｓｗｉｔｈｅｘｔｒｅｍｅｄｉｓｔｏｒ—

ｔｉｏｎｓｉｎｉｍｐａｃｔｄｙｎａｍｉｃｓ．Ｍｉｓｅｓ

ｓｔｒｅｓｓ

ｃｒｉｔｅｒｉｏｎｉｓｔａｋｅｎ，ａｎｄｔｈｅ

ａｃｃｕｒａｔｅ

ＳＰＨ—ＦＥＭ（ｆｉｎｉｔｅｅｌｅｍｅｎｔｍｅｔｈ—

ｏｄ）ｃｏｎｖｅｒｓｉｏｎａｌｇｏｒｉｔｈｍａｌｌｏｗｓ

ｔｏｕｓｅ

ａｎｄｅｆｆｉｃｉｅｎｔＦＥＭｉｎｔｈｅｌｏｗｅｒｄｉｓｔｏｒｔｉｏｎｒｅｇｉｏｎｓ，

ａ

ａｎｄＳＰＨｉｎｔｈｅｈｉｇｈｅｒｄｉｓｔｏｒｔｉｏｎｒｅｇｉｏｎｓ．Ｔｈｅｐｅｒｆｏｒａｔｉｏｎｏｆ

ｔｉｎｇ

ａ

ｃｙｌｉｎｄｒｉｃａｌｓｔｅｅｌｐｒｏｊｅｃｔｉｌｅｉｍｐａｃ—

ｐｌａｔｅ

ｔａｒｇｅｔ

ｉｓｓｉｍｕｌａｔｅｄｉｎ３Ｄｕｓｉｎｇｔｈｅ

ＳＰＨ—ＦＥＭ

ｃｏｎｖｅｒｓｉｏｎａｌｇｏｒｉｔｈｍ．Ｔｈｅｇｏｏｄａｇｒｅｅ—

ａｃｃｕｒａｃｙ

ｍｅｎｔｂｅｔｗｅｅｎｔｈｅｃｏｍｐｕｔｅｄｒｅｓｕｌｔａｎｄｔｈｅｅｘｐｅｒｉｍｅｎｔａｌｏｂｓｅｒｖａｔｉｏｎｓｈｏｗｓｔｈｅ

ｓｕｐｅｒｉｏｒ—

ｉｔｙｏｆｔｈｉｓａｌｇｏｒｉｔｈｍ．ＣｏｒｒｅｓｐｏｎｄｉｎｇｃｏｎｖｅｒｓｉｏｎｃｒｉｔｅｒｉｏｎｓｈｏｕｌｄｂｅｃｈｏｓｅｎａｃｃｏｒｄｉｎｇｔＯｔｈｅｍａｔｅ—ｒｉａｌｆａｉｌｕｒｅｍｏｄｅｉｎｅｎｇｉｎｅｅｒｉｎｇａｐｐｌｉｃａｔｉｏｎ．

Ｋｅｙｗｏｒｄｓ：ｓｍｏｏｔｈｅｄｐａｒｔｉｃｌｅｈｙｄｒｏｄｙｎａｍｉｃｓ；ｆｉｎｉｔｅｅｌｅｍｅｎｔｍｅｔｈｏｄ；ｃｏｎｖｅｒｓｉｏｎ

ｉｍｐａｃｔ

ａｌｇｏｒｉｔｈｍ；

冲击侵彻过程中遇到的大变形问题是计算力学的难点，有限元中的拉格朗日方法在计算此类问题时不可避免地会遇到网格畸变和缠绕问题，造成计算无法进行下去．通常的解决方案是采用单元侵蚀

算法将严重畸变的有限元从模型中删除．尽管该方案将单元质量分配到单元节点，保证了质量守恒，但单元的体积却被删除，降低了计算精度‘１１．欧拉方法可以避免网格缠绕，但要精确确定模型的自由表面、

收稿日期：２０１０—０７—１５．作者简介：张志春（１９８２一），男，博士生；强洪夫（联系人），男，教授，博士生导师．基金项目：国家

“９７３计划”资助项目（６１３１０２０３０３０２）；教育部跨世纪优秀人才培养计划资助项目（ＮＣＥＴ一４１３８Ｃ２ＸＢ）；第二炮兵工程学院创新性探索研究基金资助项目（ＫＸ２００８１７２）．网络出版时间：２０１０—１０—１９

网络出版地址：ｈｔｔｐ：／／ｗｗｗ．ｃｎｋｉ．ｎｅｔ／ｋｃｍｓ／ｄｅｔａｉｌ／６１—１０６９一Ｔ．２０１０１０１９．２０５８．００５．ｈｔｍｌ

西安交通大学学报第４５卷

变形边界和运动交界面等则存在困难．光滑粒子流体动力学（ＳＰＨ）［２’３］是一种无网格方法，通过使用一系列任意分布的粒子来求解相应的偏微分方程组，得到精确稳定的数值解．由于这是一种纯拉格朗日方法，并且计算中不需要借助网格，所以在处理涉及大变形、自由表面等问题时显示出明显的优势，目前在冲击侵彻问题的数值模拟研究中得到了广泛的应用．但是，对于冲击动力学问题中的小变形部分，如靶板的边缘区域，使用ＳＰＨ方法就存在明显劣势：计算精度和效率均较有限元法（ＦＥＭ）低．同时，ＳＰＨ还具有拉伸不稳定和施加边界困难等问题，在工程应用中受到了很大的限制．考虑到冲击侵彻过程中大变形仅发生在局部区域，将发生畸变的有限元转换为ＳＰＨ粒子，在大变形区域使用具有优势的ＳＰＨ，在小变形区域使用精度和效率更高的ＦＥＭ，成为计算冲击动力学问题的一种有效途径．Ｊｏｈｎｓｏｎ

ｈ＝（走）∑ｅ‰（尝）”５

㈣

式中：Ｎ。表示与粒子相关联的单元数；ｒｏ是初始单元尺寸；ｐｏ和ｐ分别表示初始密度和当前密度．ＳＰＨ—ＦＥＭ转换算法的计算流程如图２所示：对更新后的有限元进行转换依据判断，满足转换要求的有限元转换为ＳＰＨ粒子，未转换的有限元与生成的ＳＰＨ粒子分别按照各自的积分方式进行求解．继续参与计算的ＳＰＨ粒子与有限元之间存在２类问题：①ＳＰＨ粒子与有限元的固结；②ＳＰＨ粒子与有限元的接触．

等人［４｛］在解决高速冲击问题时，将严重畸变的有限元转换为一种无网格粒子算法——广义粒子算法

（ｇｅｎｅｒａｌｉｚｅｄｐａｒｔｉｃｌｅａｌｇｏｒｉｔｈｍ，ＧＰＡ），采用主面一

从点的方式，通过动量守恒、角动量守恒等条件，调整相关节点的速度，实现ＧＰＡ粒子与有限元之间的固结和接触，取得了一些有意义的成果．

本文为解决冲击侵彻过程中的大变形问题，提

计算时间步长

图１畸变有限元转换为ＳＰＨ粒子

出了一种ＳＰＨ—ＦＥＭ转换算法，以等效Ｍｉｓｅｓ应力作为转换判据，将畸变有限元转换为ＳＰＨ粒子，并且计算了ＳＰＨ粒子和有限元之间的固结和接触问题．采用ＳＰＨ－ＦＥＭ转换算法，对圆柱形钢弹正冲击钢板导致冲塞破坏的过程进行了三维数值模拟，并且与ＳＰＨ方法、ＬＳ－ＤＹＮＡ的计算结果及实验结果进行了对比，显示了该转换算法在计算精度方面的优势．１

动力学方程

ｉ点．粒亍｝触判帮

更新有限元节点位移

节固一／ｅ；抟

篓

二［

Ｉ

接触力

算法流程

ＳＰＨ—ＦＥＭ转换算法的思路是将畸变有限元转

１．１有限元向ＳＰＨ粒子的转换

换为ＳＰＨ粒子，如图１所示，图中虚线网格表示畸变的有限元，实线圆表示Ｉ型ＳＰＨ粒子（周围部分有限元发生畸变），虚线圆表示Ⅱ型ＳＰＨ粒子（周围所有有限元均发生畸变）．采用Ｍｉｓｅｓ应力作为单元转换判据，即当单元的Ｍｉｓｅｓ应力达到用户定义的某一参考值时，单元自动转换为ＳＰＨ粒子．生成ＳＰＨ粒子的质量、位置、速度、应力均与畸变单元相应节点的物理量保持一致，粒子的光滑长度ｈ可以按照单元质量守恒关系获得

ｈｔｔｐ

垂

图２

位移

徊ｉ｛；圈ｉ

ｌ回ｉ

更新ＳＰＨ粒子

质量、速度、位置、能量

一一－ｊ＿－…ｉ

…………－＿…………一Ｊ

ＳＰＨ—ＦＥＭ转换算法流程

１．２

ＳＰＨ粒子与有限元固结

冲击过程中靶板中心部分单元转换为ＳＰＨ粒

子后，为确保ＳＰＨ粒子与未畸变单元之间物理量的连续性，需要考虑两者之间的固结问题．本文采用的方法是将ＳＰＨ粒子固结在有限元节点，并且将有限元节点纳入到ＳＰＨ粒子的邻近搜索范围，即任何位于邻近搜索范围内的ＳＰＨ粒子和有限元节点都被加入到ＳＰＨ粒子的临近列表．将ＳＰＨ粒子的积分

ｔ？｝嘲、删．、ｄｘｂ．ｃｎ

第１期张志春。等：一种光滑粒子流体动力学一有限元法转换算法及其在冲击动力学中的应用

１０７

方式作如下修正

Ｎ

Ｎｎ

分法是条件稳定的，即时间步长必须小于由该问题求解方程的性质所决定的某个临界值．ＳＰＨ部分采用跳蛙格式求解Ｎａｖｉｅｒ－Ｓｔｏｋｅｓ方程，时间步长与最小光滑长度成比例．在拉格朗日框架下，从有限元到ＳＰＨ粒子的转换要求两者的积分必须是同步的，数据的传输必须是在同一时间点，这就要求在每一步计算中，ＳＰＨ和ＦＥＭ采用相同的时间步长，在此选择ＳＰＨ与ＦＥＭ时间步长中的较小者

（ｆ（ｘｉ）卜蚤（－竺，，－）ｆ（ｘｊ）Ｗ。＋蚤（，。）ｆ（ｘｊ）Ｗ。

（２）

式中：厂（ｚ）是被积函数；Ｎ和Ｎ。分别是ＳＰＨ粒子ｉ支持域内的ＳＰＨ粒子数和有限元节点数；Ⅳ是ＳＰＨ光滑核函数；ｍｊ和』Ｄ，分别表示ＳＰＨ粒子歹的质量和密度；ｍ。，和ＩＤ。，分别表示有限元节点ｊ‘的质量和密度；＜＞表示核近似算子．

使用该方案处理ＳＰＨ粒子与有限元的耦合界面时，对于ＳＰＨ部分，由于将界面附近的有限元节点纳入到ＳＰＨ粒子的临近列表，积分时ＳＰＨ的支持域不会被问题域截断，因此消除了ＳＰＨ的边界效应［６１；对于有限元部分，耦合界面的物理量由Ｉ型ＳＰＨ粒子确定，相当于对有限元施加了本质边界条件，由于有限元的形函数满足Ｋｒｏｎｅｃｈｅｒｄｅｌｔａ张量性质，故该方案是可行的．

１．３

＆阱删一ｍｉｎ｛△￡ｓＰＨ，＆洲｝

ｔ

（４）

式中：ＳＰＨ时间步长＆ｓＰＨ—ａ＿ｎ，ＦＥＭ时间步长

‘

ｒ

＆吲≤导，其中ｈ是光滑长度，Ｌ是最小单元尺寸，

Ｌ

ｃ是材料的声速，ａ是时间步长比例系数．

２圆柱形钢弹正冲击钢板的

三维数值计算

冲击侵彻是涉及到几何、材料和接触非线性的复杂动力学问题，不仅包含大应变以及不同应变率和温度条件下的应力一应变关系，而且涉及损伤和失效，很难获得精确的解析解．本文采用ＳＰＨ－ＦＥＭ转换算法对圆柱形平头钢弹正冲击钢板发生冲塞破坏的过程进行三维数值计算，其中子弹采用Ａｒｎｅｔｏｏｌ钢，靶板采用Ｗｅｌｄｏｘ

４６０

ＳＰＨ粒子与有限元接触

本文采用施加接触力的方式计算冲击过程中

ＳＰＨ粒子和有限元之间的接触问题．接触力的计算借鉴无网格接触算法的思想［７３，将有限元节点视为粒子，任何位于有限元节点支持域内的ｓＰＨ粒子都会对节点产生接触力，反之任何位于ＳＰＨ支持域内的有限元节点也都会对粒子产生接触力．施加在ＳＰＨ粒子和有限元节点上的接触力的计算形式为

Ｅ钢．计算模型的尺寸如

图３所示，子弹和靶板均采用六面体有限元离散，单

㈦一‘葶（署）（著）Ｋ竹（％落）Ｖ川州

（３）

式中：Ｎ。。表示粒子ｉ支持域内属于不同体的临近粒子数（此处将有限元节点也称为粒子）；ｒ０是粒子间距；Ａｐ。。是粒子间光滑长度的平均值；Ｋ和卵是用户自定义参数．接触力的方向由ＳＰＨ光滑核函数的梯度确定．接触力分别以外力的形式加入到ＳＰＨ动量方程和有限元动力学方程中求解．

使用该方法计算ＳＰＨ粒子和有限元之间的接触，可以充分利用ＳＰＨ粒子的临近搜索列表，并且与ＳＰＨ的求解格式保持一致．只要粒子和节点间距达到２倍光滑长度时，就在相关粒子和节点上施加接触力．与有限元处理接触问题的方法不同，该方法在每个时间步不需要定义接触界面和界面法线方向，可以方便地扩展到三维．１．４时间步长控制

如图２所示，在ＳＰＨ—ＦＥＭ转换算法中，有限元部分采用中心差分法求解显式动力学方程。中心差

ｈｔｔｐ

ｔ

图３冲击模型尺寸

ｆ｝ＶＣ＇ＩＮＶｖ＇ｏｊｄｘｂ．ｃａ

西安交通大学学报第４５卷

元总数为４２９４４．为减少计算量，仅在子弹前端和靶板中心部位进行单元加密，加密区域的单元尺寸为

１０ｑＩｎ．子弹与靶板初始间距为２×１０一ｍ，靶板四

式中：Ｐａ是损伤阈值；Ｐｒ是断裂等效塑性应变，与材料的应力三轴度、应变率和温度相关．Ｊｏｈｎｓｏｎ等人［１２］提出的剪切损伤演化模型中将Ｐｆ描述为

Ａ一［Ｄ１＋Ｄｚｅｘｐ（Ｄ３ｄ。）］

（１＋Ｄ４ｌｎｐ。）（１＋Ｄ５Ｔ’）

（７）

周采用固支边界条件，计算时间为１６０×１０一Ｓ．在ＳＰＨ的计算过程中，采用Ｍｏｎａｇｈａｎ［８］提出的人工应力法消除ＳＰＨ拉伸不稳定造成的数值断裂，采用强洪夫［９－１０］提出的完全变光滑长度法解决变光滑长度带来的影响．

靶板中心区域在冲击过程中出现大变形，采用Ｂｏｒｖｉｋ等人Ｄ１］提出的修正Ｊｏｈｎｓｏｎ－Ｃｏｏｋ本构模型和Ｇｒｕｎｅｉｓｅｎ状态方程；子弹变形和损伤较小，采用线弹性模型．修正的Ｊｏｈｎｓｏｎ－Ｃｏｏｋ本构模型包含线性热弹性、Ｍｉｓｅｓ屈服准则、各向同性应变硬化、应变率硬化、各向同性损伤引起的软化以及失效准则，其等效Ｍｉｓｅｓ屈服应力表示如下

Ｏ＇ｅｑ一（１一Ｄ）（Ａ＋Ｂｒ”）（１＋Ｃ１ｎｒ＋）（１一Ｔ”）

（５）

式中：Ｄ，～Ｄ。为材料常数；口’＝‰／％为应力三轴度，其中Ｏ＇ｍ一（以＋民＋以）／３为平均正应力．靶板的材料参数见表１．子弹的材料参数为：Ｅ一２００

ｖ＝０．３３；ｐ一７８３８ｋｇ／ｍ３．

ＧＰａ

当子弹初速Ｖｏ一２８５．４ｍ／ｓ时，数值计算结果如图４所示，相关的实验结果如图５、图６所示ｎ１｜．从中可以看出，子弹、靶塞和靶板冲塞型破坏的计算结果均与实验结果吻合较好．子弹变形很小，表明数值计算中子弹选择线弹性模型是可行的．在￡一７×

１０－６

Ｓ时，子弹与靶板接触并开始减速，子弹正前方

靶板开始加速．靶板与子弹边缘接触部位有限元的Ｍｉｓｅｓ应力达到了ＳＰＨ—ＦＥＭ转换阈值，转换为ＳＰＨ粒子．靶板被挤压变形，背部出现凸起，子弹周围的靶板出现剪切区．在子弹正前方的压缩区和子弹周围的剪切区内，单元应力继续增大，转换为ＳＰＨ粒子继续参与计算，并且部分粒子出现损伤．靶板变形继续增大，粒子损伤迅速演化，部分粒子的损伤变量达到临界值并开始失效，裂纹逐渐向靶板背部扩展．在冲击的后期，失效模式同时包含靶板内部的剪切断裂和靶板背面的拉伸损伤，最终形成ＳＰＨ粒子形式的靶塞并完全从靶板脱离．

表２给出了数值计算结果与实验结果的对比，其中‰，表示子弹余速，勘。表示靶塞速度．数值计算采用了３种算法：①ＳＰＨ—ＦＥＭ转换算法；②ＳＰＨ算法（子弹和靶板均由ＳＰＨ粒子离散，采用ＳＰＨ粒子接触算法）；③ＬＳ－ＤＹＮＡ算法［１１］（子弹和靶板均由有限元离散）．对比结果显示，ＳＰＨ—ＦＥＭ转换算

式中：Ａ、Ｂ、Ｃ、咒、ｍ是材料常数；Ｄ为损伤变量，Ｄ一０表示材料没有损伤，Ｄ一１表示材料完全失效；Ｔ。

一≠｝＝畿是一致性温度，其中Ｔｏ是室温，Ｔｍ是

材料熔点；ｒ是累积损伤塑性应变，ｒ＝（１一Ｄ）Ｐ，Ｐ

是累积塑性应变；｝。一（１一Ｄ）多＋，多。＝ｐｌｂ。，多。是

用户定义的参考应变率．实际上当材料出现宏观裂纹时，损伤变量的临界值小于１，失效准则可以简化为Ｄ＝Ｄ。≤ｌ，其中Ｄｃ是损伤变量临界值．式（５）的第１项描述由于损伤造成的材料性能退化，第２项描述应变硬化，第３项描述应变率效应，第４项描述温度软化效应．

损伤演化通常与残余应变相关，可以表示为累积塑性应变Ｐ的函数，并且只有当Ｐ大于某一阈值时，损伤才开始演化，因此Ｂｏｒｖｉｋ等人［１１］提出如下损伤演化规律

Ｌ由ｆ一加…。肛１篇，ｐ≥加

．

ｆ０，Ｐ＜Ｐｄ

∞’

法的计算结果与实验值吻合较好，并且比ＳＰＨ算法和ＬＳ－ＤＹＮＡ算法具有更高的精度．

表１靶板材料参数

ｈｔｔｐ：？｛瑚删．、ｄｘｂ．ｖｉｉ

第１期

张志春，等：一种光滑粒子流体动力学一有限元法转换算法及其在冲击动力学中的应用

１０９

表２数值计算结果与实验结果的对比

ｍ ｓ叫

１２

３

２８５４２４４２２２４

７

１７０７ｌ２９６ｌ１９

２

２２ｒＯ

』ｔ

Ⅲ

ｇ｝

４４４

２３３６１３１１０６

９８

１０９９４

１４２１２３

１８１．１１３２．６１１３。７

２２４．７１８７．７１６９．０

１８１１ｌ７

７●ｏｏ

眠

图４子弹侵彻钢板的模拟结果（ＳＰＨ—ＦＥＭ转换算法）

图５冲击后的子弹和靶塞

图６冲击后的靶板横截面

图７为ＳＰＨ—ＦＥＭ转换算法模拟的不同初速情况下子弹的速度一时间历程曲线．由曲线的斜率可知，在冲击侵彻过程中子弹和靶塞进行了多次碰撞．初次接触时，子弹前方的靶板急剧加速，子弹则急剧减速，此时靶塞速度高于子弹，之后由于靶板剪切区的限制，靶塞速度降低，与子弹再次发生碰撞，最后当靶塞从靶板完全脱离时，子弹与靶塞完成了一次２个自由体之间的碰撞，造成靶塞的速度高于子弹的残余速度．由不同的子弹初速对比可知：子弹初速越高，则侵彻靶板需要的时间越短，子弹损失的动能越多．

图７子弹的速度一时间历程曲线

；

｛

３结论

本文考虑到冲击侵彻过程中大变形仅发生在局

部区域，提出了一种ＳＰＨ—ＦＥＭ转换算法，在计算过程中将严重畸变的有限元自动转换为ＳＰＨ粒子，并且同时计算ＳＰＨ粒子与有限元之间的固结和接触．

ｈｔｔｐ：／／ｗｗⅣ．ｊｄｘｂ．ｃｎ