应用ALE流固耦合算法，使用LS-DYNA进行仿真，其基本操作设置流程如图1所示。

图 1 LS-DYNA求解入水问题操作流程

计算单位制设定

LS-DYNA中并不限定使用何种单位制，鉴于所建立模型的尺度，选择kg- cm-ms为基本单位，注意到：
1力单位=1质量单位×1加速度单位
1加速度单位=1长度单位/1时间单位的平方
按照上述换算关系，在本文中选择的单位与国际单位的换算关系如表1所示。
表 1 LS-DYNA物理量单位与国际单位换算关系

材料定义MAT及状态方程选择EOS

材料模型是用来描述材料状态变量(如应力、应变、温度及时间)的相互关系，主要是应力与应变之间的关系。状态方程是表征流体内压力、密度和温度三个热力学参数的关系式。在本次研究中涉及的材料模型有五种，状态方程两种，对应关系如表 2所示。
表 2 结构材料模型与状态方程

空气材料模型选自MAT_NULL，这种材料允许在不计算偏应力的情况下考虑状态方程。在此模型中，空气仅需要定义密度RO=1.29kg/m2即可，其余物理量由状态方程EOS_LINEAR_POLYNOMIAL定义。此模型可以输入线性多项式状态方程的系数，方程定义如下。

其中ρ⁄ρ_0 代表密度比例，ρ是MAT_NULL定义的密度，当μ<0时，C2=C6=0。
针对于空气，C0=C1=C2=C3=C6=0,C4=C5=0.4。
水材料模型与空气一致，密度为RO=1000kg/m2，状态方程本文选择EOS_GRUNEISEN，数学表达式如下，相关系数如表 3所示。

其中C为声音在水中的传播速度，a是对GRUNEISEN系数γ_0的一阶修正，S_1、S_2、S_3是无量纲系数。
表 3 水状态方程系数取值

超空泡航行体采用材料模型为*MAT_PLASTIC_KINEMATIC，该模型适用于模拟各向同性和运动硬化塑性，航行体选择45号钢金属材料，这类钢中有害杂质及非金属夹杂物含量较少，化学成分控制得也较严格，其塑性、韧性较好,运用于制造较重要的机械零件，相关参数如表 4所示。
表 4超空泡航行体材料属性系数

缓冲头罩使用材料模型和超空泡航行体一致，材料类别选择不饱和聚酯树脂复合材料，不饱和聚酯树脂是热固性树脂中发展较快的品种之一[52]，其力学性能优良，强度和脆性较好，入水时易发生脆断，相关参数设置如表 5所示。
表 5 缓冲头罩材料属性系数

最后缓冲垫材料模型为MAT_CRUSHABLE_FOAM，用于模拟泡沫，可设置阻尼和截止拉应力。缓冲垫材料选择硬质聚氨酯泡沫，硬质聚氨酯泡沫是一种吸能显著的材料，广泛应用于航空、航海等领域。硬质聚氨酯泡沫属于弹塑性材料，针对其力学性能，参数如表 6所示。在本构关系中，应力应变曲线表现出明显的三阶段特征：弹性变形、塑性屈服和压缩密实，典型的应力应变曲线如图 2所示。在LS-DYNA中可以通过关键字DEFINE_CURVE自定义本构曲线。
表 6 缓冲垫材料属性系数

图 2 弹塑性材料典型本构曲线

沙漏控制*HOURGLASS

LS-DYNA在计算本文模型时使用的是单点高斯积分，这种方法极大的提高了运算效率，也非常有利于计算大变形问题。但是单点积分很容易引起零能模式，也称沙漏模式，沙漏是一种比结构全局响应高得多的频率震荡的零能变形模式，这是单元刚度矩阵中秩不足导致的，由于积分点的不足，沙漏模式是一种在数学上是稳定的、但在物理上无法实现的状态。它们通常没有刚度，变形呈现锯齿形网格。单点实体单元的沙漏模式如图 3所示。

图3单点沙漏模式变形状态
在分析中沙漏变形的出现会使结果无效，所以应尽量避免和减小，如果总体沙漏能超过总体内能的10%，那么分析可能是无效的，有时5%也是不被允许的，所以非常必要对其进行控制。在LS-DYNA中可以使用关键字*HOURGLASS进行控制。在该关键字中比较重要的两个属性是沙漏控制类型IHQ与沙漏系数QM，在LS-DYNA的前处理器LS-PrePost 4.7版本中，有11种沙漏控制类型，针对于本文的实体单元，提供了6种类型，在本模型中使用到了两种沙漏类型，如表 7所示。
表 7 沙漏属性及参数设置

IHQ=1是标准LS-DYNA粘性类型，此种类型适用于高速变形问题，当沙漏系数QM>0.15时，会带来不稳定，对于此种工况QM取值一般小于0.1；IHQ=6是Belytschko-Bindeman 沙漏控制，该类型适用于二维和三维实体图元，针对航行体和缓冲头罩这种剪切弹性模量远小于弹性模量的材料来说，QM取值在0.001~0.1，针对于缓冲垫材料，其QM值应控制在0.5到1.0之间，经测试选择0.6可得到较为理想的结果。

单元属性设置*SECTION

在该关键字中可以设置单元属性，包括单元算法、积分准则、节点厚度和截面特性，本模型使用到的是*SECTION_SOLID，其可以设定实体单元的类型。在此关键字中比较重要的属性是单元公式选项ELFORM，该属性提供了43种可选类型，本文模型选择如表 8所示。
表8 结构单元属性设置

ELFORM=5/7/11/12均可适用于ALE算法，11最为常用，1为恒定应力图元，是LS-DYNA缺省值。

流固耦合作用定义*CONSTRAINED

CONSTRAINED中提供了一种约束自由度的方法，在本文ALE流固耦合使用关键字CONSTRAINED_LAGRANGE_IN_SOLID定义，耦合方式采用罚函数耦合。本关键字涉及的属性较多，对算例结果的影响很大，介绍如下。
SLAVE：从结构，选择结构域航行体、缓冲头罩和缓冲垫。
MASTER：主结构，选择流体域空气和水。
NQUAD：每一个耦合Lagrangian面上有NQUAD×NQUAD个耦合点，对于泄露控制十分有效，按照算例所取网格尺寸关系，针对研究模型取NQUAD=3。
CEYPE：流固耦合方法属性，本文为体单元采取罚函数，对应编号为4。
DIREC：耦合方向，高速入水以压缩变形为主，取DIREC=2。
MCOUP：耦合物质，由于本文问题涉及气液固三相变化，所以需要全部物质进行耦合，这会消耗较多的计算资源，取MCOUP=0。
PFAC：罚函数因子，用于计算分布在从段SLAVE与主段MASTER之间的力。此属性关乎整个算例的结果，是十分重要的参数，当PFAC过大时，会使计算网格出现负体积，会造成计算意外终止；而过小会使主从段之间的耦合力太小，出现大面积的体积泄露。此参数合适的数值会根据模型的不同而发生变化，所以需要反复调整，本文选取0.0003、0.001、0.01和0.05四个工况，在不加装缓冲头罩和缓冲垫的情况下（即只有超空泡航行体入水）监测的加速度曲线如图4所示。

图 4 不同PFAC工况下航行体加速度曲线
从结果中可以看出，不同的PFAC产生了不同的结果，在PFAC=0.05时出现了负体积，在时间为2.24ms时程序终止运算，而PFAC=0.0003时出现了网格泄露，泄露状态如图5所示，可以看出有部分流体已经进入了封闭的超空泡航行体内。选择合适的PFAC不能只看加速度峰值状态，而是需要比较滑动界面能，在无摩擦或小摩擦情况下，滑动界面能是接触弹簧保持的势能，在碰撞过程中，能量的转化应该是接触弹簧的势能转化为动能，动能转化为变形能，所以计算的滑动界面能是非物理的，应当控制在0左右，尤其是不能出现过大负值，以上对应工况的滑动界面能统计如图 6所示。

图 5 PFAC=0.0003工况下造成的体积泄露

图 6不同PFAC工况下滑动界面能曲线
可以看出在PFAC=0.05时滑动界面能为负值，此时取值过大，应调小处理，在PFAC=0.0003时滑动界面能稳定在300左右，此时取值过小，应放大处理，经过反复调试，最终在PFAC=0.001时滑动界面能稳定在0左右，说明此值较为合适。
FRCMIN：缺省值为0.5，减小此值会更早的激活耦合，适用于在高速碰撞情况下避免耦合泄露，取FRCMIN=0.15。
ILEAK：耦合泄露控制标识，取ILEAK=1，代表弱耦合控制。
ALE算法参数定义ALE
此节中引用ALE_MULTI-MATERIAL_GROUP，简称AMMG，此关键字可以进行界面重构，根据物质间能否混合将各种材料定义在不同的多物质组AMMG中，每一AMMG相当于一种单独的“流体”。定义该关键字后将允许一个单元中容纳多种ALE物质材料，本文中定义了空气和水两种AMMG。
接触类型CONTACT
关键字CONTACT提供了一种处理不相交部分之间交互的方法。本文中使用到CONTACT_AUTOMATIC_SURFACE_TO_SURFACE，算例中需要涉及三个接触，即超空泡航行体-缓冲头罩、超空泡航行体-缓冲垫、缓冲头罩-缓冲垫。
时间步/ALE/能量控制CONTROL
此关键字可以对多种参数进行控制，算例中启用了CONTROL_ALE、CONTROL_TERMINATION和CONTROL_TIMESTEP，这些设置较为简单，关键字及属性值如表9所示。
表 9*CONTROL关键字属性数值

DCT：激活交替输运逻辑的标识，DCT=-1代表开启交替输运逻辑。
NADV：输运步之间的循环数，通常设为1。
METH：输运方法，METH=2表示带有HIS的Van Leer方法，二阶精度。
AFAC：ALE光滑加权因子，简单平均。
ENDTIM：计算终止时间，取4ms。
TSSFAC：时步缩放系数，高速冲击问题中取0.667。
HGEN、SLNTEN：沙漏能、滑动界面能处于开启监测状态。

初始速度INITIAL、荷载LOAD、边界条件*BOUNDARY

初始速度INITIAL_VELOCITY_GENERATION设定超空泡航行体、缓冲垫和缓冲头罩为300m/s，荷载LOAD定义所有结构重力加速度g，边界条件设定流体域四周为无反射区，下表面固定。

后处理*DATABASE

此关键字用来输出关注的物理量，本算例使用了能量监测GLSTAT、MATSUM和界面力RCFORC，其余后处理诸如运动特性、应力应变等均在LS-PrePost中可以直接得到。

能量校核

在每个算例结束后都要进行能量校核，以保证能量变化在可接受的范围内，具体监测结果如表 10所示。总能量变化率在1.2%，处于可接受的范围。另外，对于航行体的沙漏能/内能需要小于10%，否则会出现不合理的变形，如图7所示为两者的曲线走向，沙漏能与内能峰值为5.1%，符合要求。
表 10航行体能量变化汇总

图 7航行体沙漏能与内能