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

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

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

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

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

其中ρ⁄ρ_0 代表密度比例，ρ是MAT_NULL定义的密度，当μ0.15时，会带来不稳定，对于此种工况QM取值一般小于0.1；IHQ=6是Belytschko-Bindeman 沙漏控制，该类型适用于二维和三维实体图元，针对航行体和缓冲头罩这种剪切弹性模量远小于弹性模量的材料来说，QM取值在0.001~0.1，针对于缓冲垫材料，其QM值应控制在0.5到1.0之间，经测试选择0.6可得到较为理想的结果。

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

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

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

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

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

图 7航行体沙漏能与内能