圆柱体垂直入水三维数值模拟及影响因素研究

邹田春， 高 飞， 魏家威， 解 江

(中国民航大学 适航学院，天津 300300)

结构物入水冲击涉及气-固-液三者之间的相互作用，是典型的流固耦合过程，在结构入水过程中，结构与自由液面冲击产生的较大冲击力极易导致外部结构破损，内部构件的失效，返回舱入水、鱼雷入水、水上飞机的起降等都属于结构入水的范畴。因此结构物入水的研究在军事、船舶、海洋工程、航空航天等领域都有重要的工程应用背景。

近十年来研究人员针对结构入水冲击问题的分析研究大多以楔形体[1-4]、圆柱体[5-8]、球体[9-10]等作为研究对象，采取试验和仿真方法开展研究。Shah进行了刚性楔形体入水冲击试验，利用高速摄像成功捕捉到了整个二维流体流动动态过程，并运用光滑粒子流体动力学(smoothed particle hydrodynamics，SPH)方法建立相应模型，对SPH方法模拟结构入水的准确性进行了验证。夏维学等以圆柱体为研究对象开展入水试验研究，对入水空泡壁面运动特性与空泡演化过程之间的关系进行分析，并基于空泡壁面运动特性划分空泡演化各阶段，分析了空泡演化在各阶段的流动机制。Luo等[11]对三维钢质楔形体进行了一系列自由落体入水冲击试验，研究楔形体入水冲击加速度、压力及应力响应，并建立了楔形体试验件的有限元模型，通过在模型上施加试验测得的入水冲击压力来预测结构的瞬态响应，并对水弹性效应进行讨论。施书文等[12]针对入水抨击问题提出了一种弱可压光滑粒子水动力与光滑点插值流固耦合法，运用该方法计算了弹性斜板低速入水和楔形体高速入水问题，计算结果与试验值和半解析解吻合良好，证明了该方法对模拟弹性体入水问题的有效性。结构入水仿真方法主要有有限体积法(finite volume method，FVM)[13-17]、SPH方法[18-20]以及任意拉格朗日欧拉(arbitrary Lagrangian-Eulerian，ALE)方法[21-26]等。ALE方法兼具拉格朗日法和欧拉法的优点，在结构入水问题研究领域采用较多且较为成熟。在传统ALE理论和算法不变的情况下，S-ALE(structured-arbitrary Lagrange-Euler)方法则利用网格的规律性，将ALE模型规则网格通过输入三维点位坐标方式由程序本身自行创建网格，极大降低数值模拟的计算成本。

针对结构入水问题的研究，在国防、军事上大多以楔形体为研究对象，当前我国正在大力发展民用客机且目前民用大型客机都是典型筒段机身结构，考虑到民机机身水上迫降的工程背景，本文选择圆柱体为研究对象。基于圆柱体垂直入水冲击试验研究结果，本文建立了圆柱体三维数值模型，采用S-ALE算法对圆柱体垂直入水冲击过程进行了数值模拟，分析其入水冲击响应特性，并通过对比圆柱体入水位移及自由液面动态变化来验证模型有效性；基于验证的数值模型，进一步研究圆柱体入水速度、质量和长度对其垂直入水冲击的影响规律。

1 基本理论

1.1 S-ALE基础方程

S-ALE在实现上是独立于ALE而全新开发的，程序简洁高效，但在理论方面S-ALE求解器与ALE完全相同。在ALE方法的描述中，引入了除Lagrange和Euler坐标系之外的第三个任意参照坐标系，与参照坐标相关的材料微商可以采用式(1)描述

(1)

式中：Xi为Lagrange坐标；xi为Euler坐标；wi为相对速度。

假设物质速度为v，网格速度为u，引入相对速度w=v-u。则ALE方法的控制方程可以由如下守恒方程给出

(1)质量守恒方程

(2)

(2)动量守恒方程

控制固定域内的牛顿流体流动问题的增强形式由控制方程和对应的初始及边界条件组成，控制流体问题的方程是Navier-Stokes方程的ALE描述

(3)

应力张量为

σij=-pδij+μ(vi,j+vj,i)

(4)

式中：p为压强；bi为单位质量的体积力；δij为Kronecker δ函数；μ为动力黏性系数。该方程可与下列边界条件和初始条件联立求解

(5)

σijnj=0,Γ2边界

(6)

Γ1∪Γ2=Γ,Γ1∩Γ2=0

(7)

式中：Г为计算域的完整边界，分为Г1和Г2两个部分，其中，Г1为固体边界，Г2为自由液面边界[27]；上标0为该参数的初始指定值；nj为边界的外法线单位向量。假设整个计算域在t=0时刻的速度场是静止的，即

vi(xi,0)=0

(8)

(3)能量守恒方程

(9)

采用算子分离法求解ALE方程，该方法将每个时间步分为两个阶段。首先执行拉格朗日过程，此时网格随物质一起运动。该过程计算速度、内外力引起的内能改变量，平衡方程为

(10)

(11)

单元边界没有物质流过，所以满足质量守恒。第二个计算阶段为对流阶段，在对流阶段之前采用VOF(volume of fluid)方法捕捉不同材料之间的界面位置，然后对穿过单元边界的质量运输、内能以及动量进行计算，即将拉格朗日阶段的移位网格重映射回初始位置或者任意位置。

对于每个节点，速度和位移按下列等式进行更新

(12)

xn+1=xn+1+Δtun+1/2

(13)

外力矢量又分为体积力、边界力、非反射边界条件及接触力。本文采用标准罚函数算法计算流体与结构之间的接触力。罚函数算法中允许从节点穿透到主段，所以在每一个时间步，系统首先检查从节点是否穿透主面，未穿透前，不做任何处理，当发生穿透，则引入一个与穿透距离d成比例的界面接触力F并将其分配到流体单元节点上，计算公式为

F=ki·d

(14)

式中,ki为基于主从节点质量模型特性的刚度系数。

1.2 水和空气状态方程

圆柱体垂直入水模型中，水和空气采用空材料模型，通过一个截断压力值来定义压力的更小边界值，式(15)为任意方向黏性应力的表达式。通过状态方程来定义压力，其中水采用Gruneisen状态方程，空气采用Polynomial状态方程。

(15)

Gruneisen状态方程对于压缩材料由式(16)计算，膨胀材料由式(17)计算。

(γ0+au)E

(16)

P=ρ0C2u+(γ0+au)E

(17)

式中：C为冲击波速度；E为流体的比内能；ρ0为流体密度；γ0为Gruneisen伽马；a为对γ0的一阶体积修正；S1，S2，S3为曲线斜率的系数；由于材料不受初始载荷，初始能量项E和修正因子a均设为1。

线性多项式状态方程

P=C0+C1μ+C2μ2+C3μ3+(C4+C5μ+C6μ2)E

(18)

式中：C0，C1，C2，C3，C4，C5，C6为自定义的常数；E为单位初始体积的内能。

水的Gruneisen状态方程参数，如表1所示。空气的Polynomial状态方程参数，如表2所示。

表1 水的Gruneisen状态方程参数Tab.1 Parameters of Gruneisen state equation for water

表2 空气的Polynomial状态方程参数Tab.2 Parameters of Polynomial state equation for air

2 有限元建模

圆柱体垂直入水三维模型，如图1所示。圆柱体直径D=160 mm，长L=295 mm；水和空气域长L=810 mm，宽W=660 mm，高H=700 mm其中水深h=480 mm。模型采用实体单元建模，圆柱体采用Lagrange单元，空气和水采用Euler单元。通过定义ALE_STRUCTURED_MESH_CONTROL_POINTS和ALE_STRUCTURED_MESH两个关键字在计算过程中由程序自动生成水和空气的网格。

为防止圆柱体在试验过程中发生变形，将PVC圆柱体内部用泡沫材料进行实心填充，所以圆柱体采用刚体进行建模。水和空气采用NULL材料模型和流体状态方程来模拟，并对模型施加向下的重力加速度g，可以产生水的压力梯度。

约束流体域底面单元节点的所有自由度来模拟水池底面边界条件；在流体域四周设置无反射边界条件，模拟无限水域，可减小流体域大小，降低计算成本，还可防止边界产生的反射波与结构发生流固耦合作用。运用S-ALE罚函数耦合算法计算圆柱体与流体之间的耦合作用。

基于Russo等开展的圆柱体入水试验结果，选取质量为m0(2.214 kg)试验件提升1 m高度的入水试验工况，如图2所示。作为本文圆柱体入水冲击仿真的初始工况。初始时刻圆柱体距液面高度为零，圆柱体初始入水速度为v0，通过运用GetData Graph Digitizer确定初始入水速度v0为3.6 m/s。将圆柱体入水冲击仿真结果与试验结果进行对比验证。

图1 数值模拟模型Fig.1 Numerical simulation model

图2 不同质量试验件提升高度与入水时刻速度关系图Fig.2 Relation diagram of lifting height and water entry velocities of test pieces with different masses

3 有限元仿真及验证

3.1 S-ALE与ALE计算对比

采用10 mm网格分别建立ALE模型和S-ALE模型，圆柱体与空气、水网格比例为1 ∶1，设置计算总时长为0.25 s，使用LS-DYNA 971求解器在主频为2.8 GHz的工作站上采用8核进行计算。两种模型在相同设置和工况下计算时间与最大位移值，如表3所示。S-ALE模型计算时间相对于ALE模型减少了20.2%。两种模型计算得到的圆柱体位移-时间曲线对比，如图3所示。从图3可知，变化趋势完全一致。其中，ALE模型得到最大位移值为0.231 6 m，S-ALE模型得到最大位移值为0.231 8 m，两者相差可以忽略不计。因此，S-ALE模型在获得相同仿真精度条件下可以大大减少计算时间。

表3 计算时间及最大位移值Tab.3 Calculation time and maximum displacement

图3 S-ALE模型与ALE模型仿真结果对比Fig.3 Comparison of simulation results between S-ALE model and ALE model

3.2 模型验证

圆柱体与空气、水的网格尺寸相同，分别选取7.5 mm、10 mm、15 mm和20 mm进行计算。采用S-ALE模型获得的圆柱体入水位移-时间曲线与试验结果对比，如图4所示。

随着网格的细化，仿真获得的圆柱体最大位移值增加，且更接近试验最大位移值(0.241 m)；网格细化到10 mm时，仿真获得的圆柱体最大位移值趋于稳定。当网格尺寸为10 mm时，圆柱体最大位移值为0.231 m，与试验值相差4.15%；网格尺寸为7.5 mm时，圆柱体最大位移值为0.235 m，与试验值相差2.49%。因此，选用7.5 mm网格尺寸，模型总网格数为898 372，圆柱体不同入水时刻获得的仿真位移值与试验值对比情况，如表4所示。位移误差均在3%以内，因此S-ALE模型能够精确模拟圆柱体入水位移时间历程。

图4 不同网格大小计算结果与试验结果对比Fig.4 Comparison of calculation results of different grid sizes and test results

表4 不同时刻圆柱体仿真位移值与试验结果对比Tab.4 The simulated displacement value of cylinder at different times compared with the experimental results

圆柱体入水仿真结果与试验结果对比,如图5所示。当圆柱体接触到自由液面时，圆柱体动能开始转移为水的动能，圆柱体所挤占空间处的水在接触处被喷射出来形成射流；圆柱体在下沉过程中高速的自由液面射流不断增长与扩大(T<τ/5)(τ为圆柱体达到最大向下位移的时刻)；在(τ/5)

图5 仿真结果和试验结果对比Fig.5 Comparison of simulation and experimental results

S-ALE模型计算得到的液面运动情况与试验之间有微小的差异，对于射流的飞溅效果，模拟结果不如试验明显，但仿真模拟基本上再现了圆柱体撞击水面、自由液面射流形成、空气腔打开、空气腔闭合的典型现象。整体来看自由液面动态仿真模拟结果与试验基本一致。

3.3 圆柱体入水速度分析

7.5 mm尺寸网格模型仿真得到的圆柱体入水速度变化曲线,如图6所示。从图6可知，圆柱体入水整个运动过程可分为撞击水面、减速下移和上浮三个阶段。圆柱体在撞击水面时受到较大冲击载荷，速度迅速下降，由于圆柱体密度小于水，在冲击波衰减过后，其受水的作用力仍大于重力，因此圆柱体在撞击水面后做加速度相对较小的减速运动，并在0.198 s时速度减小到0，此时圆柱体在水中处于悬浮状态，随后在浮力的作用下圆柱体开始上浮。

图6 圆柱体入水速度-时间曲线Fig.6 Velocity-time history curve of cylinder water entry

4 不同因素影响研究

4.1 圆柱体入水速度影响

圆柱体(D=160 mm，L=295 mm)不同入水速度下的加速度-时间历程曲线，如图7所示。从图7可知，三种速度下加速度曲线都有一较大峰值，此时段为圆柱体冲击水面时段，在较大冲击载荷作用下会对结构产生较大破坏，也是结构物入水研究关注的重点，对比加速度曲线可知，入水速度对圆柱体加速度峰值影响非常明显，随着速度的增大加速度峰值增大，而对加速度脉冲持续时间几乎没有影响。

图7 圆柱体不同入水速度加速度-时间曲线Fig.7 Acceleration-time history curves of cylinder with different water entry velocities

圆柱体(D=160 mm，L=295 mm)在不同速度下的入水位移-时间历程曲线，如图8所示。从图8可知，三种速度下圆柱体位移趋势一致，随着入水速度增大，圆柱体向下运动的越快且最大位移增大，但达到最大位移的时刻基本保持不变。

图8 圆柱体不同入水速度位移-时间曲线Fig.8 Displacement-time history curves of cylinder with different water entry velocities

4.2 圆柱体质量影响

通过在圆柱体上方增加集中质量的方式改变圆柱体质量从而获得2.214 kg，3.214 kg，4.214 kg，5.214 kg圆柱体入水冲击加速度。图9为圆柱体(D=160 mm，L=295 mm)入水速度3.6 m/s，不同质量工况下加速度曲线，如图9所示。从图9可知，随着质量增大，加速度峰值有减小趋势，但质量增加未改变加速度脉冲持续时间。圆柱体受力示意图如图10所示。由牛顿第二定律式(19)可以得到圆柱体入水冲击在垂直方向上所受水的合外力F(包括浮力、动量传输引起的附加力、形状阻力)，不同质量圆柱体入水所受阻力曲线，如图11所示。由图11可知，在圆柱体外形不变的条件下不同质量圆柱体入水冲击所受水的合外力基本不变，则由式(20)可知，外形相同的圆柱体入水冲击时，质量m越大，加速度a越小。

F-mg=ma

(19)

a=F/m-g

(20)

式中：m为圆柱体质量；g为重力加速度；a为圆柱体加速度。

图9 不同质量圆柱体入水加速度-时间曲线Fig.9 The water entry acceleration-time history curves of cylinders with different masses

不同质量圆柱体(D=160 mm，L=295 mm)以2.2 m/s速度入水时位移-时间历程曲线，如图12所示。从图12可知，由于加速度大小与质量呈负相关，则质量越大圆柱体入水减速越缓，向下运动的越快；随着质量增加，圆柱体入水最大位移及其到达时刻都增大；当质量增加到圆柱体重力大于浮力时，圆柱体入水过程最终将会达到稳定状态，保持匀速一直下落。

图10 圆柱体入水冲击时的受力图Fig.10 The force diagram of cylinder impact water

图11 不同质量圆柱体入水所受阻力曲线Fig.11 The water entry drag force curves of the cylinders with different masses

图12 不同质量圆柱体入水位移-时间历程曲线Fig.12 The water entry displacement-time history curves of cylinders with different masses

4.3 圆柱体长度影响

为研究圆柱体入水三维效应，针对给定直径D的圆柱体，保持其密度ρ不变，只改变长度L，研究长度对圆柱体入水特性的影响。

选取D=160 mm圆柱体，不同长度圆柱体入水加速度-时间历程曲线，如图13所示。不同长度(D=160 mm)圆柱体加速度峰值大小，如图14所示。从图14可知，加速度峰值随着长度的增大而增大，且当L/D>1时，增加趋势逐渐变缓，当L/D>3时，L/D的增加对加速度峰值几乎没有影响。L/D对加速度脉冲持续时间几乎没有影响。相对于速度，长度改变所造成的加速度变化幅度明显要小。

图13 不同长度(D=160 mm)圆柱体入水加速度-时间曲线Fig.13 The water entry acceleration-time history curves of cylinders(D=160 mm)of different lengths

图14 不同长度(D=160 mm)圆柱体入水加速度峰值Fig.14 The water entry peak accelerations of cylinders (D=160 mm) with different lengths

不同长度(D=80 mm)圆柱体入水位移-时间历程曲线，如图15所示。从图15可知，随着L/D增大，圆柱体最大入水位移与其达到最大位移的时刻都增大；当L/D>2时，L/D的增大对最大位移的影响逐渐降低；当L/D>4时，L/D的增加不再影响最大位移。

图15 不同长度(D=80 mm)圆柱体入水位移-时间曲线Fig.15 The water entry displacement-time history curves of cylinders(D=80 mm)with different lengths

不同长度(D=120 mm)圆柱体入水位移-时间历程曲线，如图16所示。当L/D<3时，随着L/D的增加，圆柱体最大入水位移与其达到最大位移的时刻都增大；当L/D>3时，随着L/D的增加，最大入水位移开始减小。

图16 不同长度(D=120 mm)圆柱体入水位移-时间曲线Fig.16 The water entry displacement-time curves of cylinders(D=120 mm)with different lengths

图17给出了不同长度(D=160 mm)圆柱体入水位移-时间历程曲线，当L/D<2时，圆柱体最大入水位移与其达到最大位移的时刻都增大；当L/D>2时，随着L/D的增加，最大入水位移逐渐减小。

图17 不同长度(D=160 mm)圆柱体入水位移-时间曲线Fig.17 The water entry displacement-time curves of cylinders(D=160 mm)with different lengths

L/D对圆柱体入水的影响在开始时间阶段较小，且L/D越小向下运动的越快，随着向下运动的进行，L/D对圆柱体位移的影响变的更加明显，且因直径的不同而呈现不同的规律。

圆柱体最大入水位移与不同L/D的关系曲线，如图18所示。从图18可知，随着直径的增大，最大位移增大；三种直径圆柱体，在L/D<2时，L/D对位移的影响较大；当L/D>2时，影响逐渐变缓，且不同直径呈现不同的影响规律。

图18 不同长度、不同直径圆柱体入水位移曲线Fig.18 The water entry displacement curves of cylinders with different lengths and diameters

5 结 论

(1)相对ALE模型，S-ALE模型在保证计算精度的前提下，可大大提升计算效率。S-ALE模型可以完整的再现圆柱体冲击水面、自由液面射流形成、空气腔打开、空气腔闭合的整个过程。通过圆柱体入水过程典型时刻自由液面动态以及位移的数值计算结果与试验结果对比，验证了本文所建立模型的有效性和准确性。

(2)随着圆柱体入水冲击速度的增大，圆柱体加速度峰值及最大位移值都增大，但对加速度脉冲持续时间和达到最大位移时刻没有影响。

(3)加速度大小与质量呈负相关，但圆柱体质量改变，加速度脉冲持续时间基本不变，质量增大最大位移及其到达时刻增大。

(4)长度增加，圆柱体入水冲击加速度增加，但对加速度脉冲持续时间没有影响，且相对于速度，长度改变引起的加速度变化幅度明显要小；入水位移在开始时间段受长度的影响较小，且L/D越小圆柱体向下运动的越快，随着圆柱体向下运动的进行，L/D对圆柱体位移的影响变的更加明显，且因直径的不同而呈现不同的规律。