基于LS_DYNA的弹体撞水过程流固耦合动力分析

第 22 卷第 6 期
系 统 仿 真 学 报?
Vol. 22 No. 6
2010 年 6 月
Journal of System Simulation Jun., 2010
? 1498 ?
基于LS-DYNA的弹体撞水过程流固耦合动力分析
孙 琦, 周 军，林 鹏
（西北工业大学精确制导与控制研究所，西安 710072）
摘 要：基于有限元软件 ANSYS/LS-DYNA，建立了三维弹性弹丸体在多物质流体介质中运动的
有限元动力分析模型，采用 ALE 方法对弹体撞水过程进行了流固耦合数值模拟，探讨了其中关键
控制参数对弹体入水过程的影响。仿真结果表明，该仿真模型可以解决弹体撞水流固耦合计算中流
体形变引起的数值计算困难，求解出的水—气交界面变化情况比较合理，求解出的不同撞水初态下
弹体姿态和位置变化情况，可以为弹体撞水初态的选择提供依据。
关键词：流固耦合；撞水响应；LS-DYNA；数值模拟
中图分类号：TP391.9 文献标识码：A 文章编号：1004-731X (2010) 06-1498-04
Dynamic Analysis of Fluid-Structure Interaction
for Water Impact of Projectile Using LS-DYNA
SUN Qi, ZHOU Jun, LIN Peng
(Institute of Precision Guidance and Control, Northwestern Polytechnical University, Xi’an 710072, China)
Abstract: Multi-material ALE method was used to fluid-structure interaction dynamic analysis on water impact of
three-dimensional elastic projectile, by the explicit formulation dynamic analysis finite element software LS-DYNA. The
influences of key control parameters were discussed. Numerical simulation results indicate that, this model can solve the
numerical computation difficulty originated from fluid great shape change in fluid-structure interaction dynamic analysis,
and the water-gas interface is reasonable. Pressure applied from fluid on projectile is analyzed, important references for the
design of projectile structure. Projectile attitude and position variety under different initial state is compared, which is
helpful to design appropriate initial attitude in water impact of projectile.
Key words: fluid-structure interaction; water impact of projectile; LS-DYNA; numerical simulation
引 言
1
结构物撞水的研究，对于解决鱼雷、反潜导弹、深水炸
弹、破障炮弹入水的运动，宇宙飞船返回舱的水上回收，水
上飞机的起降等工程实际问题具有重要意义
[1-2]
。当结构物
与流体发生撞击作用时，不但流体会对结构物产生强大的撞
击作用力，结构物的运动和形变也会对流场产生显著的影
响，因此两者之间存在着强烈的流固耦合作用。而高速弹体
的冲击入水，还需要考虑水—气交界面的变化情况，对于这
样具有复杂边界条件的实际工程问题，很难进行解析求解，
借助于计算机的数值分析方法进行近似求解是比较有效的
方法。随着计算机硬件和软件的高速发展，有限元方法在流
固耦

合动力学分析中得到了越来越多的应用
[3-4]
，尤其是基
于显式积分有限元软件 LS-DYNA 的 ALE 方法，被国内外
学者认为是一种分析流固耦合的有效方法
[5-8]
。
1 基于 ALE 网格的流固耦合方法
固体结构物的有限元计算通常采用 Lagrange 列式，其
收稿日期：2009-10-14 修回日期：2010-03-01
作者简介：孙琦(1982-), 男, 湖南益阳人, 博士生, 研究方向为飞行器结
构设计、导弹控制系统设计与仿真；周军(1966-), 男, 江苏常州人, 教授,
博导，研究方向为飞行器导航、制导与控制；林鹏(1978-), 男, 辽宁本
溪人, 讲师, 博士后，研究方向为飞行器结构设计、导弹先进控制理论
与应用。
单元网格与结构是重合的，网格随着结构的变形而变形。但
是，对于流固耦合问题，材料流动将造成有限元网格的严重
畸变，引起数值计算困难，导致仿真计算无法完成。为此，
我们采用 LS-DYNA 程序的 ALE (Arbitrary Lagrange-Euler )
列式来解决这个问题。ALE 列式先执行一个或几个 Lagrange
时步计算，此时单元网格随材料流动而产生变形，然后执行
ALE 时步计算：首先，保持变形后的物体边界条件，对内
部单元进行重分网格，网格的拓扑关系保持不变，称为
Smooth Step；然后，将变形网格中的单元变量（密度、能量、
应力张量等）和节点速度矢量输运到重分后的新网格中，称
为 Advection Step。
对于 ALE 网格，可以理解为有两层网格重叠在一起，
一层是空间点网格，另一层是附着在材料上并随着材料在空
间网格中运动的网格，设流体质点 K ( a , b , c )和网格点
I ( i , j , k )的运动可用下式分别表示：
( )
,
L
X = X K t， ( )
,
A
X =X I t
取 f 为流场内某一个物理量，它是坐标与时间变量的连续函
数，则 f 对网格 I 与质点 K 的随体导数可分别表示为
d
d
A
A
X
f f f X
t t X t
I
? ? ?
= + ?
? ? ?
(1)
d
d
L
L
K
X
f f f X
t t X t
? ? ?
= + ?
? ? ?
(2)
式中，? f ?t 和 ?f ?X 分别表示当地导数和位变导数，
A
?X ?t第 22 卷第 6 期 Vol. 22 No. 6
2010 年 6 月 孙琦, 等：基于 LS-DYNA 的弹体撞水过程流固耦合动力分析 Jun., 2010
? 1499 ?
和
L
?X ?t 分别表示网格点和流体质点的速度。
基于 ALE 网格结构，在 LS-DYNA 程序中可以很方便
的实现欧拉—拉格朗日耦合算法，从而实现流体—结构的耦
合计算。采用流固耦合算法模拟问题时，往往要对拉格朗日
结构进行约束，将结构的相关参量传递给流体单元。我们采
用罚函数约束的方法来实现耦合。罚函数耦合系数追踪拉格
朗日节点（弹体结构，即从物质）和欧拉流体（水和空气，
即主物质）位置间的相对位移 d ，

检查每一个从节点对主物
质表面的贯穿，如果从节点不出现贯穿，就不进行任何操作；
如果发生从节点对主物质表面的贯穿，界面力 F 就会分布
到欧拉流体的节点上。界面力的大小与发生的贯穿数量成正
比，即
i
F = k ? d(3)
式中
i
k 表示基于主从节点质量模型特性的刚度系数。
2 弹体撞水计算模型
在进行弹体撞水动力分析时，我们所主要关心的是弹体
结构受到的流体反作用的压力峰值，以便为提高弹体结构的
抗撞强度提供一个参考依据；以及在入水过程中，弹体位置
和姿态的变化情况，通过选择合适的入水初态，使整个撞水
过程中，弹体姿态满足一定的要求。
2.1 弹体材料模型
在材料本构模型中，主要有弹性、塑性、粘性以及三者
的组合。考虑到弹体高速入水时承受的巨大撞击力，有时需
要通过一定的结构来进行缓冲，存在结构的变形及失效问
题，因此，我们采用弹塑性材料来对弹体进行建模。
弹塑性材料在加载段应力与应变保持线性，当应力大于
屈服应力
s
σ 时，材料进入塑性，此后如果继续加载，应力
—应变关系仍然为线性，但是斜率发生变化。卸载曲线与加
载段曲线斜率相同，这样当完全卸载（应力等于零）后，材
料中将保留永久的塑性变形
p
ε 。一维情况下弹塑性材料本
构的描述，首先判断结构的应力状态是否达到屈服应力，如
果没有达到，则按照线弹性材料本构进行处理。如果材料内
的应力已经超过屈服应力，则按照塑性变形本构计算结构中
的应力—应变。在三维条件下，判断材料是否进入塑性可以
使用 V.Mises 屈服准则，即
2
1 1
0
2 3
ij ij s
s s ? σ= (4)
式中，
ij ij m ij
s = σ ?σ δ
为斜应力张量；
11 22 33
1
( )
3
m
σ = σ + σ + σ为平均应力。
2.2 流体材料模型
仿真模型中的流体材料包括水和空气两种介质，建模时
将二者间处理为共面，划分网格时两种材料间即为共节点。
水和空气均采用 Null 材料模型，水的状态方程采用
Gruneisen 状态方程，即
2 02
0
0
2
2 3
1 2 3
2
1 1
2 2
( )
1 ( 1)
1 ( 1)
a
C
a E
S S S
p
γ
ρ μ μ μ
γ μ
μ μ
μ
μ μ
? ? ??
+ ? ?
? ? ??
? ? ??
= + +
? ?
? ? ? ?
? ?
+ +
? ?
(5)
式中， p 为压力，E 为单位体积内能，C 、
1
S 、
2
S 、
3
S 为
材料常数。空气模型采用线性多项式状态方程，其压力计算
公式为
2 3 2
0 1 2 3 4 5 6
p = c + c μ + c μ + c μ + ( c + c μ + c μ)E(6)
2.3 显式分析求解
在建立弹体撞水的有限元动力分析模型时，外流场需包
含整个弹体运动区域，我们将流场边界到弹体的最近距离取
为弹体相应尺寸的 3~5 倍

，流场外界采用无反射边界条件，
考虑到弹体的高速运动，这样得到的流场区域往往很大，为
了将求解时间控制在可以接受的范围内，需要对时间步长进
行严格控制。
对于 LS-DYNA 程序的显式时间积分算法，若给定了模
型的材料特性，则仿真的最小时间步长由最小单元尺寸控
制：对于给定的网格划分，
min
Δt 取决于材料的纵波速，它
是材料性质的函数（密度、弹性模量和泊松比）。由于仿真
模型的外流场区域远大于弹体体积，在划分网格时，弹体网
格单元往往会比流场网格单元小很多，为了避免流场中的大
单元也采用很小的时间步长，我们通过子循环（混合时间积
分）的方法来加快分析速度。子循环方法根据计算时间步长
的大小把单元分成许多组，每一组的时间步长值可以是单元
最小时间步长的整数倍。这样，较大单元组的时间步长可能
为最小时间步长的 2 倍、3 倍或 4 倍等。通过仿真可以发现，
对于单元大小差别较大的模型，使用子循环技术可以显著加
快分析速度。
仿真模型由弹体、水和空气三部分组成，均采用
LS-DYNA 提供的 3D SOLID164 单元进行有限元分析，其离
散化结构运动方程为
Mx = P - F +H -Cx (7)
其时间积分采用显式中心差分方法，基本格式如下：
1
1/2
( ) [ ( ) ( ) ( ) ( )]
n n n n n
x t M P t F t H t Cx t
?
?
= ? + ? (8)
1/2 1/2 1
( ) ( ) ( )( ) / 2
n n n n n
t t t t t
+ ? ?
x = x + x Δ + Δ(9)
1 1/2
( ) ( ) ( )
n n n n
t t t t
? +
x = x + x Δ(10)
式中，M 为总体质量矩阵，P 为总体节点载荷向量，F 为
应力散度矢量，H 为沙漏粘性阻尼力向量， ( )
n
x t为
n
t 时刻
的节点位置坐标向量。
考虑到弹体撞水过程中涉及到水和空气两种流体介质，
采用多物质 ALE 算法进行求解，单元算法采用单点
Euler/ALE 多物质单元算法。多物质单元是指这种单元划分
的网格中，允许多种物质的流动，在同一个网格中，可以包
含多种材料的物质，因此，该算法需要求解物质在网格中的
输运过程。
此外，由于高速碰撞会在材料内部形成冲击波，使压力、第 22 卷第 6 期 Vol. 22 No. 6
2010 年 6 月 系 统 仿 真 学 报 Jun., 2010
? 1500 ?
密度、质点加速度和能量等各种物理量出现跳跃间断，给动
力学微分方程组的数值求解带来困难。所以，需要在计算中
引入人工体积粘性项来修正静水压力项。人工体积粘性 q 的
加入，可以将冲击波的强间断平滑为在相当窄的区域内急剧
变化但却是连续变化的情况。对于人工体积粘性项，
LS-DYNA 程序的标准算法是：
2
0 1
( ) 0
0 0
kk kk kk
kk
q lC l C a
q
ρ ε εε
ε
?
?
?= <
?
?= ≥
?
(11)
其中，

特征长度
3
l = V， a 为局部声速， ρ 为当前质量密
度， ( )
kk11 22 33
ε = ε + ε + ε是应变率张量的迹，
0
C 和
1
C 为无
量纲常数，通常取为 1.5 和 0.06。
引入人工体积粘性 q 后，应力计算公式被修正为：
( )
ij ij ij
σ = S + p + qδ(12)
式中， p 为压力，
ij
S 为偏应力张量。
2.4 基于 ANSYS/LS-DYNA 的建模方法
由于 LS-DYNA 本身只是一个求解器，仅能对由前处理
器输出的 K 文件进行求解。所以在进行仿真计算时，需要
通过前处理软件来建立对象几何模型、划分网格以及设置大
部分的仿真条件，通常采用 ANSYS 软件作为前处理器。但
是，由于 ANSYS 不支持 LS-DYNA 的某些关键字，在采用
ALE 算法进行分析计算时，需要对 ANSYS 软件输出的 K
文件进行文本编辑，增加相关关键字，这样才能应用
LS-DYNA 软件的多物质 ALE 算法进行流固耦合动力分析
求解。采用 ANSYS/LS-DYNA 软件对三维弹性弹丸体撞水
过程进行建模和仿真分析的主要步骤包括：
(1) 建立对象实体模型
建立仿真对象实体模型的过程在 ANSYS/LS-DYNA 的
图形用户界面中完成，包括定义单元属性、定义材料模型、
构建几何模型、划分网格、定义约束和设置初始条件等。为
了能够应用 LS-DYNA 求解器进行求解，需要将对象单元属
性指定为显示单元类型，这里我们采用 LS-DYNA Explicit>
3D SOLID164 单元来建立对象的 3D 模型。在构建好对象的
几何模型后，通过分别划分弹体结构和流场的网格来建立有
限元模型，在划分网格时需要指定各部件（PART）的单元
属性、材料属性、网格类型和网格密度。这里，为了采用多
物质 ALE 算法进行求解，需要将空气流场和水流场定义为
不同的 PART，采用不同的材料模型参数。
(2) 设置求解分析基本参数
在 ANSYS/LS-DYNA 图形用户界面的求解器 Main
Menu>Solution>中，可以设置大部分的求解分析参数，包括
能量选项、人工体积粘性选项、时间步长因子、求解时间、
结果文件输出等。根据各次仿真的具体要求，设置完所有选
项后，通过菜单项 Main Menu>Solution>Write Jobname.k 输
出模型关键字 K 文件，该文件是一个 ASCII 格式的文本文
件，其格式为 LS-DYNA 的关键字命令格式。
(3) 修改关键字文件
用文本编辑器打开前面得到的模型关键字文件，在其中
修改和添加 ANSYS 未写入而 LS-DYNA 求解器所需要的信
息字段。这里，K 文件中需要修改的主要内容包括：将
*SECTION_SOLID 关键字改为*SECTION_SOLID_ALE，其
单元算法选项 ELFORM 取 11，即 ALE 多物质单元；添加
ALE 算法控制关键字*CONTROL_ALE；添加定义多物质单
元的关键字段*ALE_MULTI-MATERIAL_GROUP，允许流
体网格中同时存在空气和水两

种物质；添加流固耦合关键字
段*CONSTRAINED_LAGRANGE_IN_SOLID，将 Lagrange
描述的固体结构部件的 ID 编号作为从属（SLAVE），将 ALE
描述的流体部件组合作为主部件（MASTER），通过罚函数
约束算法实现流固耦合。
(4) 递交 LS-DYNA 程序求解
将修改完成的关键字K文件递交LS-DYNA程序进行求
解分析，计算完成后，根据所设置的结果文件输出类型可以
得到各种后处理软件能识别的结果文件。
3 数值模拟结果与分析
基于 LS-DYNA 软件的多物质 ALE 算法，对弹塑性材
料弹丸体撞水过程（弹体入水初速 0.6 Ma）进行了三维有限
元模拟，采用 LS-PREPOST 软件对 LS-DYNA 求解出的数
据进行后处理。仿真过程中为了减少计算量，采用面对称假
设，只对一半的流场和弹丸体进行建模，并对弹体外形进行
了简化，以下是部分模拟仿真结果。
3.1 流场状态分析
图 1 所示为入水过程中某时刻的流体密度等值线图，图
图 1 流体密度等值线图
图 2 压力等值线图第 22 卷第 6 期 Vol. 22 No. 6
2010 年 6 月 孙琦, 等：基于 LS-DYNA 的弹体撞水过程流固耦合动力分析 Jun., 2010
? 1501 ?
中上半部分是空气，下半部分是水，没有考虑海浪的影响，
从图中可以看出因弹丸撞击而引起的水面起伏，以及流体密
度的变化。图 2 所示为入水过程中某时刻的压力等值线图，
为了更好的显示弹丸和水的压力等值线在三维空间的分布
情况，图中在显示时已将空气流场隐藏掉了。
3.2 冲击响应分析
图 3 所示为某初始撞水条件下弹丸体垂向加速度随时
间变化曲线，仿真模型中采用的时间单位是 ms，图中纵坐
标为弹丸体垂向加速度，单位是
3 2
10 m/ms
?
。从仿真结果来
看，弹体受到的过载（加速度）是一个急剧上升然后逐渐减
小的过程。弹体不同初始条件下撞水过程的仿真结果表明，
撞水初速、入水姿态和弹体材料是影响过载峰值和过载响应
曲线的主要因素。通常情况下，入水时弹丸体垂向速度越大、
弹体纵轴与水平面的夹角越小、弹体材料的弹性模量越大，
则过载峰值越大，过载响应越剧烈。
图 3 弹丸体垂向加速度随时间变化曲线
图 4 所示为撞水过程中某时刻弹丸体的应力分布等值
线图，可以看出弹体的应力主要集中在头部，同时图中还列
出了弹体的网格结构，并标出了所选定用来进行对比分析的
三个网格节点 No:116 、 No:984 和 No:1106。
图 4 弹丸体应力分布等值线图
3.3 弹体运动情况
图 5 所示为图 4 中所标出的节点 No:116 、 No:984 和
No:1106 的垂向速度随时间变化曲线，速度的单位是 m/ms 。
从图中可以看出三个节点的垂向速度并非完全一致，这说明弹
体在入水过程中

姿态发生了变化，根据弹体上不同位置的节点
的相对速度变化曲线可以得到弹体姿态变化曲线。为了使整个
入水过程中弹体的姿态都满足要求，需要设计合适的入水初
态，因此，需要对不同初态下的弹体入水过程进行仿真分析。
图 5 节点垂向速度随时间变化曲线
4 结论
基于 LS-DYNA 程序的 ALE 有限元方法处理弹体撞水
问题具有较大的优越性，它几乎不受几何外形、边界条件以
及载荷情况的限制，可以对弹体结构、空气流场和水流场进
行非定常耦合计算，可以克服许多解析方法无法解决的困
难，如撞水引起的液面变化，弹体结构的变形耦合，撞水的
空气层影响等等。
仿真结果表明，该仿真模型可以很好的解决弹体撞水流
固耦合计算中流体形变引起的数值计算困难，对于亚音速弹
体撞水的三维流固耦合数值模拟计算，通过合理的划分网格
和选择算法参数，可以使撞水过程的沙漏能量小于总能量的
10% ，从而满足显式时间积分算法对沙漏控制的要求。采
用该仿真模型求解出的水—气交界面变化情况比较合理，求
解出的弹体受到的流体压力峰值，可以为弹体结构设计提供
理论依据；求解出的不同撞水初态下弹体姿态和位置变化情
况，可以为弹体撞水初态的选择提供依据。
参考文献:
[1] 尹莉, 钱勤, 罗宁, 等. 结构撞水的流固耦合动力分析[J]. 华中科
技大学学报, 2006, 23(2): 5-7.
[2] 马春生, 黄世霖, 张金换, 等. LS-DYNA 的 ALE 方法在飞船返回
舱着落仿真中的应用[J]. 清华大学学报, 2006, 46(8): 1455-1457.
[3] 张馨, 王善, 陈振勇, 等. 水下接触爆炸作用下加筋板的动态响应
分析[J]. 系统仿真学报, 2007, 19(2): 257-260.
[4] Jinwook Lee, Reuben R Rohrschneider, Stephen M Ruffin, Robert D
Braun. Fluid-Structure Analysis of a Clamped Ballute in Titan's
atmosphere [C]// AIAA Computational Fluid Dynamics Conference,
Miami, USA, 2007. USA: AIAA, 2007-4308.
[5] 陈向东, 金先龙, 李政. 波浪作用下软体浮囊动态响应仿真方法
研究[J]. 系统仿真学报, 2007, 19(20): 4616-4619.
[6] Benjamin A Tutt, Anthony P Taylor. The use of LS-DYNA to
Simulate the Inflation of a Parachute Canopy [C]// AIAA
Aerodynamic Decelerator Systems Technology Conference and
Seminar, Munich, Germany. USA: AIAA, 2005-1608.
[7] Edwin L Fasanella, Karen E Jackson. Dynamic Impact Tolerance of
Shuttle RCC Leading Edge Panels Using LS-DYNA [C]// AIAA Joint
Propulsion Conference and Exhibit, Tucson, USA, 2005. USA: AIAA,
2005-3631.
[8] Vijay K Goyal, Carlos A Huertas, Jose R Borrero, Tomas R Leutwiler.
Robust Bird-Strike Modeling Based On ALE Formulation Using
LS-DYNA [C]// AIAA Structures, Structural Dynamics, and
Materials Con. Newport, USA, 2006. USA: AIAA, 2006-1759.
Y
-
elocityv
lesuRtant
igidR

odyB
rceleAation
E-03)(
0 50 100
150
Time
LS-DYNA user input
0
0.005
0.01
0.015
0.02
0.03
0.025
0.035
Material No:
A1
0 50 100 150
Time
LS-DYNA user input
-0.016
-0.015
-0.014
-0.012
-0.013
-0.011
Node no.
1106
A
C984
B
116