基于ALE方法的楔形体入水分析

李上明 屈明

摘要：

基于ALE方法分析楔形体入水问题，楔形体采用拉格朗日网格离散，空气和水采用ALE网格离散。将楔形体视为刚体，空气和水的力学行为分别采用Gamma定律和GRUNEISEN状态方程模拟。讨论ALE流固耦合关键字中罚函数罚因子的取值方法，提出相应的建议原则；分析水域截断边界对楔形体响应的影响，给出模拟无限水域的截断边界位置的建议值；分析楔形体表面压力振荡的原因，提出楔形体表面压力获取方法。通过试验结果对比分析，验证方法的合理性。

关键词：

入水； ALE方法； 截断边界； 流固耦合； 无限水域； 罚函数

中图分类号： O347.5

文献标志码： B

Wedge water entry analysis based on ALE method

LI Shangming， QU Ming

（Institute of Systems Engineering， China Academy of Engineering Physics， Mianyang 621999， Sichuan， China）

Abstract：

The wedge water entry problem is analyzed based on ALE method. The wedge is dispersed using Lagrange meshes. The air and water are dispersed using ALE meshes. The wedge is considered as a rigid body. The mechanics behaviors of air and water are simulated by Gamma law and GRUNEISEN state equation respectively. The method to get the penalty factors of penalty function of ALE fluidsolid coupling keywords is discussed， then accordingly principle is proposed. The influence of the truncated boundary of water area on the wedge response is analyzed. The position of truncated boundary is proposed for modelling infinite water area. The reason of pressure vibration on wedge surface is explained， and a method is proposed to obtain the pressure on the wedge surface. The rationality of the method is verified by the comparison and analysis of test results.

Key words：

water entry； ALE method； truncated boundary； fluidsolid coupling； infinite water area； penalty function

0 引 言

随着人们对飞机水上迫降、返回舱水上回收、空投鱼雷入水等问题的进一步关注，越来越多的研究开始重视结构冲击入水的动态响应特性。响应评估的手段有理论计算、数值模拟和试验等。与理论计算和试验相比，数值模拟具有适用范围广、成本低等特点，经常被研究人员和工程师采用。

目前，常用于结构入水的数值模拟方法有基于VOF的CFD方法[15]、SPH方法[69]、ALE方法[1015]等。CFD方法大多集中于刚体结构入水分析方面，其获取的结构入水压力波动较小；SPH方法适合模拟大变形问题，常用于结构入水分析中水域的模拟，但其边界处理还存在较大的提升空间，且压力波动较大。ALE方法虽然在流固耦合面上的压力波動也较大，但因为其已经被集成在LSDYNA中，所以获得广泛的应用。调整相关参数可以有效控制ALE压力的波动，提高分析的合理性。

对结构入水流固耦合模拟问题，LSDYNA中的ALE方法基于统一的连续介质动量守恒方程，采用类似固体之间的接触算法模拟界面力，利用显式有限元法统一计算，固体采用本构关系模拟，流体采用状态方程模拟，适合模拟结构大变形和以层流特征为主的流体大位移等情况。根据ALE显式动力学有限元法的计算特点，该方法在模拟快速运动（即惯性效应较强的情况）比较有效，但因其采用罚函数法进行流固耦合载荷的传递，故流固耦合界面力会出现振荡。本文通过参数分析，讨论罚函数法相关罚因子的选择方法和压力波动的原因，进而确定结构入水分析的方法。

1 流固耦合建模方法

基于某二维楔形体结构入水问题，讨论LSDYNA流固耦合建模方法。该楔形体为等腰楔形体，底边长为1.20 m，底角为25°。楔形体顶点

左、右两边水和空气的宽为1.85 m，水深为1.00 m，空气高为0.40 m。楔形体视为刚体，采用平面应变单元离散；水域和空气域采用ALE单元离散。空气和水的自由液体界面网格共节点处理。楔形体网格与空气、水的离散网格重叠。楔形体附近水域最小网格为2 mm，相对远处区域的网格逐步放大，以减小网格数量。楔形体网格大小约2 mm，其网格离散和压力测点布置见图1。

1.1 物理参数

楔形体密度为467 kg/m3，弹性模量为500 MPa，泊松比为0，采用*mat_rigid模拟。水和空气采用*mat_null和状态方程模拟，其中：水采用GRUNEISEN状态方程模拟，见式（1）；空气通过线性多项式状态方程模拟，气体的Gamma定律状态方程见式（2）。

1.2 流固耦合关键字

为有效进行结构入水流固耦合分析，通常要求在LSDYNA的输入文件k文件中包括如下关键字：*ALE_MULTIMATERIAL_GROUP关键字，用以指定并跟踪水和空气材料及其边界；*CONTROL_ALE关键字，用以控制ALE单元的相关算法和部分边界处理方法；*CONSTRAINED_LAGRANGE_IN_SOLID關键字，用以描述流体与固体间的耦合作用，在流固耦合面上实现载荷传递。本文采用罚函数法进行流固耦合分析。

罚函数的物理意义相当于在结构节点（从节点）与流体单元表面（被穿透表面、主面）之间放置一个法向弹簧，以限制结构节点与流体单元面之间的穿透[13]，其耦合力

式中：K和V分别为被穿透流体单元的体积模量和体积；A为连接从节点的结构单元平均面积；p为罚函数的罚因子。p用于调整流固耦合接触刚度，其取值的大小影响流固耦合面的相互作用力。式（3）和（4）中，除p为输入量外，其余参数均由

LSDYNA软件自动计算。

2 算例分析

针对某二维楔形体入水，主要讨论左、右两侧截断边界和罚函数的罚因子对分析结果的影响。边界条件设为：水域底部固定垂直（y向）位移，左、右两侧水域固定水平（x向）位移。初始条件为：楔形体以5.05 m/s初始速度撞击水面，其余均为静止。外力条件为：考虑静水压力和重力加速度。因该问题具有对称性，故在分析中楔形体只考虑y向位移，x向位移约束为0。

截断边界讨论的主要目的是探索截断边界远近对计算结果的影响，以评估模拟无限水域时截断边界距离的选取原则。罚函数罚因子讨论的主要目的是评估其大小的选取原则，以便确保模拟结果的合理性。

2.1 截断边界讨论

考虑楔形体顶点两侧水域3种不同截断长度下的楔形体入水响应，即：工况一为水域侧面单边长为1.20 m，即楔形体半边长度的2.0倍；工况二为水域侧面单边长1.85 m，即楔形体半边长度的3.1倍；工况三为水域侧面单边长3.20 m，即楔形体半边长度的5.3倍。

在不同工况下楔形体的速度差见图2。在0.10 s内，水域侧面长从1.20 m增大到3.20 m时，最大速度差不到12 μm/s。该数据表明：当水域截断长度远大于结构尺寸时，截断边界的远近对结构入水冲击的整体速度影响较小。因此，在下面的分析中，水域截断长度均取1.85 m。

2.2 罚因子讨论

在不同罚因子下楔形体入水的计算结果见图3。由此可知：不同的罚因子p对速度和流固耦合面的合力影响很小。当p=0.01时，楔形体在0.10 s时刻的速度最大，为0.795 m/s；当p=1.00时，楔形体在0.10 s时刻的速度最小，为0.794 m/s。罚因子从0.01变化到1.00，0.10 s时刻的剩余速度只相差0.001 m/s，相对误差不到0.13%。因此，罚因子的变化对剩余速度的影响可忽略不计。p=0.10时，楔形体滑移能最小。在纯弹性碰撞中，滑移能应完全转化为动能，因此分析时应尽可能控制滑移能在很小的范围内，越接近0越好，所以建议p取0.10。

2.3 楔形体压力讨论

1 khz以上高频滤波后楔形体表面测点压力的响应曲线见图4。由于测点选择在楔形体的壁面上，其压力存在较大的数值振荡。

经反复数值试验发现，当某时间步的流固耦合界面刚好超过流体网格的边界节点且超过量较小时，该流体节点会将流固耦合面附近的流体网格“穿透”，形成向内“穿刺”现象，“穿刺”距离未超过1个单元格，见图5。这可能是引起边界压力波动的根本原因。

为进一步探索压力波动的原因，采用数值试验方法分析测点1和测点6附近不同位置处的压力。距离楔形体测点6法向距离分别为0、1和2个网格处的压力曲线见图6a）。3个位置到楔形体壁面法向的距离分别为0、2和4 mm。计算中心区域网格大小为2 mm（对角线长约为2.8 mm）。由图6a）可以看出：距测点6越远的点，压力振荡越小；与测点6距离为0和1个网格处的点，因未避免流体向内

个网格处的点压力振荡几乎消失。

距离楔形体测点1法向距离为3和4 mm处的压力响应曲线见图6b）。由此可以看出，这2个点的压力振荡很小，压力响应基本一致。这说明在法向距离大于1个网格对角线距离的位置压力较平稳（压力曲线自身小的振荡是由小体积应变引起的），可避免ALE网格“穿刺”引起的振荡。为获得相对稳定的压力值，耦合界面压力记录位置应该选择在距离流固耦合界面法向1个网格对角线处。在此位置上，尽管会存在一定的压力损失，但所获得的压力是稳定的，能相对真实地反映实际压力。

3 试验结果对比分析

为说明上述方法的合理性，将模拟结果与相应的试验结果进行对比。

文献[16]楔形体速度的试验数据与本文数值模拟结果对比见图7。由此可知：数值模拟的结果与ZHAO模型的结果一致性最好，并介于von Karman模型和Wagner模型结果之间；与试验结果相比，数值模拟结果中间一段存在一定差异，其原因可能是本文计算简化成平面应变问题，只考虑水池长度方向上的边界影响，未考虑水池宽度方向的边界影响。

试验与数值模拟的测点压力对比见图8。在压力波形上，数值模拟与试验较一致。靠近楔形体顶面的测点数值模拟压力曲线存在较小的振荡。该振荡可能是由水池长度方向的边界反射和数值计算中流固耦合界面接近自由表面引起的。

随着测点到楔形体顶点距离的增加，试验中测点压力上升的起始时间越来越早于数值模拟，在最远的测点12上升时间差达到约5 ms。这可能是由试验的三维效应和数值模拟中所记录的位置与其法向存在3 mm距离共同导致的。随着试验时间的增加，因楔形体厚度为1.2 m而水池在该方向上只有2.0 m厚，其压力反射效果逐渐明显，进而导致峰值前移。

部分测点的试验与数值模拟的压力峰值及其相对偏差见表1。由此可以看出：测点1的试验与数值模拟结果误差最大，达到21.9%；误差最小点为测点5，在给定有效数字范围内，其峰值正好一致。对冲击问题而言，该偏差在工程分析中完全可以接受。

4 结 论

基于ALE方法分析研究楔形体结构入水问题，并通过与试验数据的对比验证该方法的合理性。该方法可应用于以惯性效应为主的结构入水流固耦合问题。基于该方法，通过系列数值分析，获得以下建议。

（1） 采用罚函数模拟流固耦合时，因滑移能应完全转化为动能，建议以最小滑移能来确定罚因子的取值。本文算例根据数值模拟结果建议罚函数罚因子取0.10。

（2） 模拟无限水域时，其截断边界应远离结构，建议距离结构3倍以上结构尺寸为宜。

（3） 由于结构表面水压力振荡较大，建议以距离结构表面法向距离为1.4倍单元长度网格对角线位置处的压力为流固耦合面的表面压力，评估结构表面的压力情况。

参考文献：

[1] 邹星， 李海涛， 宗智. 基于VOF 模型的结构物出水过程数值模拟[J]. 武汉理工大学学报（信息与管理工程版）， 2012， 34（5）： 558561. DOI： 10.3963/j.issn.20953852.2012.5.007.

[2] 马庆鹏， 魏英杰， 王聪， 等. 锥头圆柱体高速入水空泡数值模拟[J]. 北京航空航天大学学报， 2014， 40（2）： 204209.

[3] JIANG C X， SHUAI Z J， ZHANG X Y， et al. Numerical study on transient behavior of waterentry supercavitating flow around a cylindrical projectile influenced by turbulent dragreducing additives[J]. Applied Thermal Engineering， 2016， 104： 450460. DOI： 10.1016/j.applthermaleng.2016.05.102.

[4] FACCI A L ， PORFIRI M ， UBERTINI S. Threedimensional water entry of a solid body： A computational study[J]. Journal of Fluids and Structures， 2016， 66： 3653. DOI： 10.1016/j.jfluidstructs.2016.07.015.

[5] ABRAHAM J P， GORMAN J M， FRANCO R， et al. Modeling and numerical simulation of forces acting on a sphere during earlywater entry[J]. Ocean Engineering， 2014， 76： 19. DOI： 10.1016/j.oceaneng.2013.11.015.

[6] OGER G， DORING M， ALESSANDRINI B， et al. Twodimensional SPH simulations of wedge water entries[J]. Journal of Computational Physics， 2006， 213（2）： 803822. DOI： 10.1016/j.jcp.2005.09.004.

[7] JI Z， XU F， TAKAHASHI A， et al. Large scale water entry simulation with smoothed particle hydrodynamics on single and multiGPU systems[J]. Computer Physics Communications， 2016， 209： 112. DOI： 10.1016/j.cpc.2016.05.016.

[8] GONG K， SHAO S D， LIU H， et al. Twophase SPH simulation of fluidstructure interactions[J]. Journal of Fluids and Structures， 2016， 65： 155179. DOI： 10.1016/j.jfluidstructs.2016.05.012.

[9] SHAO S D. Incompressible SPH simulation of water entry of a freefalling object[J]. International Journal for Numerical Methods in Fluids， 2009， 59（1）： 91115. DOI： 10.1002/fld.1813.

[10] 孙琦， 周军， 林鹏. 基于LSDYNA的弹体撞水过程流固耦合动力分析[J]. 系统仿真学报， 2010， 22（6）： 14981501.

[11] 郑金伟， 宗智. 三维刚体椭圆头结构高速倾斜入水冲击模拟[J]. 船海工程， 2012， 41（3）： 79. DOI： 10.3963/j.issn.16717953.2012.03.003.

[12] 张岳青， 徐绯， 金思雅， 等. 飞船返回舱水上回收的冲击响应和入水姿态分析[J]. 振动与冲击， 2014， 33（18）： 204208. DOI： 110.13465/j.cnki.jvs.2014.18.033.

[13] WANG S， SOARES C G. Numerical study on the water impact of 3D bodies by an explicit finite element method[J]. Ocean Engineering， 2014， 78： 7388. DOI： 10.1016/j.oceaneng.2013.12.008.

[14] 田媛， 刘均， 汪浩， 等. 砰击载荷下轻质波纹夹芯夹层板动力响应特性分析[J]. 船舶力学， 2016， 20（10）： 12991307. DOI： 10.3969/j.issn.10077294.2016.10.010.

[15] 李艷臣， 熊伟， 方佩文. 抛落式玻璃钢救生艇入水分析[J]. 船舶工程， 2016， 38（2）： 1417. DOI： 10.13788/j.cnki.cbgc.2016.02.014.

[16] YETTOU E， DESROCHERS A， CHAMPOUX Y. Experimental study on water impact of a symmetrical wedge[J]. Fluid Dynamics Research， 2006， 38（1）： 4766. DOI： 10.1016/j.fluiddyn.2005.09.003.

（编辑 武晓英）