对结构高速入水的研究，在导弹、鱼雷入水、航行器水上回收、船舶砰击和水上飞机迫降等实际问题上具有重要的工程意义。鱼雷等武器多采用薄壳结构形式，在高速入水过程将面临着结构破坏变形、内部仪器失灵等问题，鱼雷是否能进入正常的作战状态决定于鱼雷的毁伤性能。回转体入水过程涉及气-液-固三相耦合作用^[1]，在高速入水过程中流固耦合现象更复杂。回转体在高速入水瞬间受到极大的冲击载荷，结构可能发生较大变形甚至破坏。因此，在研究过程中需考虑流固耦合作用，分析结构的动响应特征，为结构强度设计提供依据。在对高速入水过程的流固耦合研究中，理论和实验的方法往往难以实施。随着计算机和数值技术的发展，数值方法在解决入水问题上不受几何外形和边界条件等因素的限制，具有一定的优越性。

对入水问题，已有了大量的理论和实验研究。1929年，von Karman^[2]提出并采用了附加质量和动量守恒的方法研究入水冲击问题，在理论上求解了二维楔形体入水的冲击载荷。在此基础上，1932年，Wagner^[3]采用伯努利方程，得到了小攻角楔形体浸水面上的压力分布。他们为后来入水问题的理论研究奠定了基础。Gavrilenko等^[4]考虑流体的可压缩性，对边界条件进行Laplace变化，在理论上分析了不同形状的轴对称刚体结构的入水。Worthington等^[5]最早采用高速摄影技术，完整地观察了小球入水的完整过程，对过程中的自由液面变化和空泡演变过程进行了分析。May等^[6-7]研究了不同结构参数的钢球入水空泡特性，并分析了不同速度入水下钢球的载荷参数。Abelson^[8]开展了锥头圆柱体以速度76.2 m/s入水的实验，分析了入水角度与压力分布的关系。随着数值计算能力的发展，数值模拟方法逐步成为解决入水问题的重要研究方法。Panahi^[9]基于有限体积方法和动网格技术，分析了结构物入水过程中的运动规律和动力特征。Erfanian等^[10]考虑流固耦合现象，对球鼻艏入水过程进行了数值模拟，提高了数值精度。Neaves等^[11]引入空化模型和Tait状态并考虑流体的可压缩性，开展了高速射弹垂直入水过程的数值模拟，较好地分析了流域的压力场分布情况。

目前，在我国，对于入水问题也有了广泛的研究。孙辉等^[12]通过实验研究了楔形体入水流固耦合动响应参数特征和不同参数的结构响应规律。张伟等^[13]、郭子涛等^[14]进行了速度35～160 m/s、不同头型参数的入水实验，研究了头型参数对弹道稳定性的影响并得到了圆柱形结构入水速度衰减的理论公式。陈诚等^[15]开展了超空泡航行器的小角度倾斜入水实验，研究了入水速度和空化器面积对冲击载荷的影响。Yan等^[16]开展了AUV航行器高速入水实验，并基于有限元和SPH粒子法耦合的方法对航行器的冲击载荷进行了分析。钱铖铖等^[17]采用动网格技术对射弹高速入水进行了数值模拟，得到了阻力、压力和入水速度的关系。张佳悦等^[18]采用VOF方法和动网格技术对航行体高速垂直入水过程进行了数值模拟，得到了航行器阻力系数与入水速度的关系及空泡发展规律。何春涛等^[19]用有限体积方法结合动网格技术，对结构垂直入水进行了数值模拟，通过与实验结果对比验证了数值算法的有效性，进一步分析了空泡闭合特性及空泡演变规律。侯昭等^[20]基于雷诺平均方法并采用重叠网格技术，建立了六自由度圆柱倾斜入水的模拟方法，深入讨论了空泡形态、弹道特性和受力情况。黄志刚等^[21]采用ALE方法，深入分析了不同壁厚回转体100 m/s入水时冲击特性和结构强度，指出结构形式是影响入水回转体结构强度的重要因素。

由以上研究可以看出，入水冲击问题受到广泛关注并有了大量有价值的研究成果。目前，研究重点多集中在低速入水问题，作为一个冲击时间极短的强非线性问题，高速入水过程中涉及流固耦合作用的实验和数值研究开展较少。因此，采用流固耦合数值方法，开展高速入水过程中弹塑性材料结构的动响应规律、流体动力特性和空泡演变规律的研究，可以加深对高速入水流场特征的理解，为高速入水结构物设计提供依据，同时为下一步开展高速入水数值计算提供参考。

1.   数值计算模型

1.1   流体基本控制方程

在流体求解器中，采用均质平衡流模型中的VOF模型作为多相流模型。VOF方法将气、汽、液组成的混合介质看作一种可变密度的单一流体，在整个流域内实现速度场，压力场共享，各相介质成分由各相体积分数α来描述。该方法可以有效地捕捉自由液面处的变形，在模拟多相流界面之间运动方面应用广泛。

连续性方程为：

式中：ρ_m为混合流体密度，各相介质密度由体积分数α决定，u为混合物平均速度，u_q为第q相速度。

动量方程为：

式中：p为压力，${\tau _{ij}} = \mu \left[ {\left( {\dfrac{{\partial {u_i}}}{{\partial {x_j}}}{\rm{ + }}\dfrac{{\partial {u_j}}}{{\partial {x_i}}}} \right) - \dfrac{2}{3}{\delta _{ij}}\dfrac{{\partial {u_k}}}{{\partial {x_k}}}} \right]$为黏性剪切应力，$\mu = \displaystyle\sum\limits_{q = 1}^n {{\alpha _q}{u_q}}$为动力黏性系数，S为源项。

对于第q相，体积分数方程为：

式中：ρ_q为第q相的密度，${\dot m}_{pq}$为第p相到第q相的质量传递率，$S_{\alpha q}$为质量源项。

体积分数守恒方程为：

湍流模型采用了应用范围较广的SST k-w模型，该模型综合了k-w和k-ε模型的优势，在近壁面处采用k-w模型而在边界层外采用k-ε模型。SST k-w湍流模型包含了修正的湍流黏性公式，考虑了湍流的剪切应力效应，能更精确地模拟反压力梯度引起的分离点和分离区大小。其输运方程为：

式中：k为湍动能，w为湍流耗散率，G_k为层流速度梯度产生的湍动能，G_w为湍动能速度梯度，Y_k为波动扩张引起的湍动能耗散，Г_k、Г_w分别为k、w的有效扩散项，D_w为正交扩散项，S_k、S_w分别为用户自定义源项。Г_k、Г_w分别为：

式中：σ_k、σ_w分别为k和w的普朗特数，F₁和F₂为混合函数，μ_t为湍流黏性系数，a₁=0.31为模型参数，S为应变率，对于高雷诺数模型$a^*=1$。

1.2   固体求解器单元和材料模型

鱼雷壳体等水下武器一般采用薄壳结构，为了使回转体在高速入水过程中的结构动响应特性更明显，本文中在倾斜入水数值模拟算例的固体求解器中采用壳单元模型。固体求解器中采用数值积分方法，分别计算厚度方向每个积分点的应力和应变，并允许非线性材料行为。采取通用壳单元S4R形式，该单元模型在多数情况下能提供比较精确的模拟结果。

材料本构关系等材料性能数据对数值模拟的精确性和有效性的影响很大，尤其对高速入水冲击等非线性问题。Johnson-Cook本构关系是一个经验性的黏塑性模型，用于高应变率和高温情况^[22]。Johnson-Cook材料模型形式简单、使用方便，能够很好地反映金属材料的性能，在冲击动力学上有广泛的应用。其应力应变关系为：

式中：σ为等效应力，ε_eq为等效塑性应变，A、B、C、n、m为模型参数，可通过参考温度和参考应变率下的拉伸实验获得，$\varepsilon_{\rm{eq}}^*$为无量纲等效应变率，${T^*} = \dfrac{{T - {T_{\rm{t}}}}}{{{T_{\rm{m}}} - {T_{\rm{t}}}}}$为相对温度，T_t为转化温度一般取室温，T_m为熔化温度。

45钢的主要材料参数^[23]分别为：密度7 850 kg/m³，杨氏模量210 GPa，泊松比0.29；Johnson-Cook本构关系参数分别为：A=595 MPa，B=580 MPa，C=0.023，n=0.133，m=2.03，T_m=1 492 ℃，T_t=25 ℃。

1.3   耦合求解模式

STAR-CCM+和ABAQUS的双向耦合是一种强流固耦合形式，双向耦合考虑流体的压力和剪切力对固体作用引起变形，也同时考虑变形导致流体流动的变化，当固体求解器和流体求解器都采用隐式格式时，该耦合模式允许STAR-CCM+与ABAQUS在一个时间步内多次交换数据。因此，需增加一个内部迭代过程。在内迭代耦合计算时，固体求解器先根据流场中的初始条件计算结构的变形量和位移量估值，耦合交界面处STAR-CCM+流体域界面上的载荷和ABAQUS求解器中耦合界面节点处的载荷需要满足载荷边界条件：

式中：${{\sigma}} _{{\varGamma }}^{\rm{s}}$、${{\sigma}} _{{\varGamma }}^{\rm{f}}$为固体和流体表面载荷。固体求解器把计算结果传递给流体求解器后，流体求解器对边界进行一致性修正，根据耦合界面的位移条件进行计算：

式中：d_s、d_f为耦合交界面处固体和流体位移。

流体求解器采用网格变形技术实现网格的变形和移动，完成一步耦合计算并继续更新压力和剪切力，重复至最大时间或满足位移收敛条件后进行下一步耦合。在双向耦合过程中，信息在流体求解器和固体求解器中进行更新和传递，并在一定间隔内进行多次交换，称为耦合时间步。耦合时间步的选择取决于问题的耦合强度，对于复杂的流固耦合问题，耦合时间步一般与固体及流体计算的时间步相同，以实现收敛和稳定性。过程如图1所示。

图  1  双向耦合求解过程

Figure  1.  Computing process of fluid-structure interaction method

下载: 全尺寸图片 幻灯片

在流固耦合过程中，固体结构被施加形变和刚体运动，因此采用网格变形技术且需适应刚体运动影响。STAR-CCM+采用morphing网格变形技术，在固体网格形变结束后，流体网格通过映射到固体网格节点上的方式产生变形移动。在有限元分析中，此方法称为任意拉格朗日-欧拉(ALE)方法。在计算流体力学中，当流体网格与固体网格共形时，流体求解器通过额外求解输运方程获得固体的运动信息，通过保角映射来获取固体的形变信息。

1.4   计算模型设置

为了进一步研究倾斜入水过程中结构动响应情况和流场变化，采用双向流固耦合方法，开展回转体倾斜60°入水过程的数值模拟分析。如图2所示，计算区域分为空气域和水域，空气域尺寸为1.5 m×1.0 m×0.5 m，水域尺寸为1.5 m×1.0 m×1.2 m。回转体入水速度为60 m/s，与水平面夹角为60°，弹体长0.197 m、直径0.05 m、壁厚3 mm，采用45钢Johnson-Cook材料模型。

图  2  计算域

Figure  2.  Computational domain

下载: 全尺寸图片 幻灯片

设计算域四周和底部为速度入口，顶部为压力出口，避免壁面反射的影响模拟无限水域；设圆柱体表面为壁面，设重叠网格与流体域的界面为重叠网格边界。如图3所示，计算时在水面运动区域和圆柱体周围进行网格加密，其中流体网格数量为3 285 858，固体网格数量为16 432。

图  3  流体网格模型

Figure  3.  Mesh of fluid domain

下载: 全尺寸图片 幻灯片

2.   计算结果和分析

2.1   方法有效性验证

为了验证本次数值模拟方法的有效性，选取文献[24]中的实验模型进行计算，并与实验结果进行对比。实验系统主要包括发射系统、水箱、高速摄影机、照明装置等，水箱尺寸为3.3 m×3.3 m×6.0 m。实验模型为45钢实心圆柱体，直径6 mm，实测质量4.88 g，入水速度106.6 m/s。数值模拟中，模型尺寸、质量和入水速度与实验一致，45钢的弹塑性阶段使用Johnson-Cook材料模型。通过对比数值模拟和入水实验结果中的速度衰减、位移和不同时刻的空泡演变规律，验证流固耦合算法的有效性。

图4为速度曲线的流固耦合模拟和实验结果。在回转体高速撞水瞬间，头部会受到一个极大的冲击力，体现为在触水瞬间回转体存在一个剧烈的速度衰减，双向流固耦合模拟可以清楚地捕捉到该现象。同时也可以看出，速度曲线的模拟结果和实验结果吻合较好。由图5可以看出，侵蚀位移的数值模拟结果和实验结果趋近一致，证明了流固耦合数值模拟方法的有效性。

图  4  速度衰减曲线

Figure  4.  Velocity attenuation curves

下载: 全尺寸图片 幻灯片

图  5  侵蚀位移曲线

Figure  5.  Erosion displacement curves

下载: 全尺寸图片 幻灯片

图6给出了回转体以速度106.8 m/s入水时，不同时刻实验和数值模拟的空泡形态对比。由图6可以清楚地观察到撞水开空泡、空泡发展、空泡闭合、空泡溃灭等过程。在触水阶段，回转体撞击水面对水产生冲击波，空泡开始形成。在0.8 ms时，回转体完全入水，空泡完全包裹回转体，水表面收到冲击产生飞溅现象。在2.0 ms时，随着回转体继续运动，空泡逐渐拉长、体积不断扩展、直径将要达到最大，并开始产生回射流。可以看出，数值模拟的空泡演变结果与实验结果具有较高的一致性，进一步验证了流固耦合算法的可行有效。

图  6  不同时刻的空泡形态

Figure  6.  Cavity features at different times

下载: 全尺寸图片 幻灯片

2.2   速度60 m/s倾斜入水

先进行入水速度60 m/s下回转体的结构动响应分析和流场分析。图7为回转体轴向、底部径向单元分布示意图，图8为底部、侧面和顶部中心单元的Mises应力曲线。由图8可以看出：回转体底部应力峰值远大于顶部应力和侧面中心处峰值且出现时间早于二者，因此在回转体倾斜入水过程中，应力主要集中在回转体底部；在入水过程中，回转体的有效应力曲线在17～19 ms时会出现一个脉动峰值，该现象出现的原因将在下面分析。

图  7  回转体轴向、底部径向的单元分布

Figure  7.  Element distribution of revolution body

下载: 全尺寸图片 幻灯片

图  8  回转体的应力沿轴向分布曲线

Figure  8.  Equivalent stress curves of elements along axis

下载: 全尺寸图片 幻灯片

图9为不同时刻回转体底部的应力分布，可以清楚地观察到应力主要集中在回转体底部中心和圆周边缘。在触水瞬间，圆周边缘部分与水接触，受到极大的冲击载荷，应力集中在触水圆周边缘处；随后，应力向底部中心处转移，扩展至单元中心处和整个圆周边缘；最终，应力始终集中在底部中心处和圆周处，保持稳定状态。图10～11为回转体底部中心单元、半径中心单元和圆周处单元的有效应力、应变曲线，可以看出：底部中心应力的峰值为171 MPa，大于半径中心处和圆周处的；底部中心应力衰减频率大于半径中心和圆周处的，在18.8 ms左右应力发生波动产生脉冲峰值，三者应变情况与应力情况趋势保持一致。

图  9  不同时刻回转体底部应力分布

Figure  9.  Stress distributions of revolution body bottom at different times

下载: 全尺寸图片 幻灯片

图  10  回转体底部的应力沿径向分布曲线

Figure  10.  Equivalent stress curves of bottom elements along radius

下载: 全尺寸图片 幻灯片

图  11  回转体底部的应变沿径向分布曲线

Figure  11.  Strain curves of bottom elements along radius

下载: 全尺寸图片 幻灯片

图12～13为不同时刻空泡轮廓图和流场压力云图，可进一步分析回转体入水过程中结构响应的脉动情况。在约17.0 ms时，空泡发生表面闭合，在约2.0 ms时，尾部空泡收缩使空泡内部产生一个短暂的高压区；由于回转体阻力主要为压差阻力，空泡内部的高压区使回转体首尾压差降低，回转体加速度呈下降趋势；空泡内部压力升高，导致回转体整体受压变大，应力、应变相应增加。在19.5 ms时，空泡尾部收缩完成，空泡内部高压区不明显，回转体整体加速度变大且应力应变减小。在回转体的有效应力曲线上，这个现象体现为17.0～19.0 ms出现的脉动峰值。

图  12  不同时刻的空泡演变

Figure  12.  Cavity features at different times

下载: 全尺寸图片 幻灯片

图  13  不同时刻的应力分布

Figure  13.  Stress distributions at different times

下载: 全尺寸图片 幻灯片

采用流体固体求解器协同模拟方法，可以清楚地观察到回转体倾斜入水过程的完整空泡演变过程，如图14所示。回转体经历撞水开空泡、空泡形成、空泡发展、空泡颈缩和空泡闭合的过程。初始时刻，回转体撞水后，水面发生流动分离向回转体两侧运动，撞水部分右侧的水面隆起，回转体轴向空泡呈现不对称形态，轴向右侧空泡大于左侧，随着入水时间增加，空泡不对称性变弱；在4 ms时，空泡完全包裹回转体；在12 ms前，空泡沿径向和回转体轴向扩张，空泡充分发展；在17 ms左右，空表面闭合，产生回射流；然后，空泡尾部开始收缩，空泡即将开始溃灭。

图  14  速度60 m/s倾斜入水的空泡演变

Figure  14.  Cavity features for oblique water-entry at v₀=60 m/s

下载: 全尺寸图片 幻灯片

2.3   入水速度的影响

为了研究不同入水速度对于回转体的载荷特性和流场演变的影响，基于已建立的流固耦合方法，同时开展了入水速度60、80、100 m/s下的数值模拟。在流固耦合数值模拟前，先比较刚体和变形体底部中心压力和速度的变化情况。图15为刚体和变形体底部中心点压力曲线。表1为不同入水速度下圆柱底部中心点的压力峰值，相比刚体结果，流固耦合方法的回转体底面中心点压力曲线的峰值较小，且随着时间会有明显的脉动现象。这是因为，弹性体圆柱入水时其底部触水部分会发生变形来起到缓冲作用。因为缓冲变形呈反复状态，弹性回转体的压力曲线产生脉动，一段时间后，随着变形量的减少，脉动幅值也逐渐减小。由图16可见：初始时刻，弹性回转体速度和刚体回转体速度变化基本相同；随后，弹性回转体在入水过程中产生变形来缓冲受到的冲击力，反映在速度曲线上，弹性回转体速度衰减将慢于刚体回转体。

表  1  不同入水速度的压力峰值

Table  1.  Peak pressures at different initial water entry velocities

入水速度/(m·s−1)	压力峰值/MPa
	CFD	FSI
60	9.93	9.42
80	16.03	15.59
100	27.92	25.12

下载: 导出CSV 

| 显示表格

图  15  入水速度80 m/s时回转体底部中心点压力曲线

Figure  15.  Pressure curve of bottom central point at v₀=80 m/s

下载: 全尺寸图片 幻灯片

图  16  速度80 m/s入水时速度曲线

Figure  16.  Velocity curves at v₀=80 m/s

下载: 全尺寸图片 幻灯片

图17为采用流固耦合方法的回转体以不同速度入水的归一化速度衰减曲线。在初始时刻，由于受到撞水瞬间产生的巨大冲击力，回转体速度衰减剧烈且速度越大衰减程度越剧烈；随后，由于冲击力迅速降低，在短时间内回转体速度衰减程度放缓且衰减程度仍与入水速度有关，入水速度越大衰减程度越大；最后，在阻力、惯性力和重力的联合作用下，速度呈类似线性下降，速度衰减状态趋于稳定且不同入水速度的衰减程度近乎一致。

图  17  归一化的速度衰减曲线

Figure  17.  Normalized velocity attenuation curves

下载: 全尺寸图片 幻灯片

图18为不同入水速度下回转体底部中心单元的等效应力曲线。回转体以不同速度入水，其底部中心单元的等效应力变化趋势保持一致，入水瞬间回转体应力急剧升高，形成瞬时峰值，随后应力峰迅速下降，最后呈震荡衰减状态。应力变化趋势与入水速度有关，等效应力峰值与入水速度平方相关，入水速度越大其峰值出现时间越早，且应力衰减震荡现象越明显，衰减幅度越大。回转体100 m/s入水时，底部中心单元等效应力峰值为596.4 MPa，略微超出材料的屈服应力595 MPa，回转体底部产生了小幅度的塑性应变。由此可见，在入水速度超过100 m/s时，该回转体模型会导致结构产生塑性形变甚至失效，所以在设计回转体模型时，还应该考虑入水瞬间的冲击力作用下的结构强度问题。

图  18  回转体中心单元的有效应力曲线

Figure  18.  Equivalent stress curves of central elements

下载: 全尺寸图片 幻灯片

图19为不同入水速度下流场不同时刻的空泡形态曲线。在4 ms时，空泡处于敞开阶段，空泡后部分呈外凸状态，空泡完全包裹回转体；在8 ms时，空泡充分发展，且入水速度越高空泡长度越长，半径也越大；在12 ms时，空泡发展完毕，空泡径向尺寸开始收缩，随着入水位移加深，在水压力和空泡表面张力作用下，入水速度80、100 m/s的回转体开始出现颈缩现象，60 m/s的回转体颈缩现象明显；在16 ms时，空泡下半部分直径较小，上半部分直径较大，因为回转体速度已经大幅度降低，回转体底部附近流体获得的动能较低，可见60 m/s的回转体即将发生表面闭合，而80、100 m/s的回转体表面闭合的时间更晚。因为入水速度较高时，在径向水表面会获得较高的扩张速度，直径也会扩展更大，表面闭合在压差作用下向轴向收紧所需要的时间也更长。

图  19  空泡形态对比

Figure  19.  Comparisons of cavity features

下载: 全尺寸图片 幻灯片

3.   结　论

采用双向流固耦合方法，对回转体倾斜入水过程中结构载荷特性和流场演变进行了数值模拟，得到了如下结论。

（1）回转体倾斜高速入水过程中，载荷先集中在触水部分边缘处，后集中于底部中心处，因空泡表面闭合后收缩，会产生载荷脉动峰值。

（2）采用流固耦合方法，转体入水载荷结果小于刚体计算的，因弹性体圆柱入水时其底部触水部分会发生变形起到缓冲作用，而缓冲变形的反复，使弹性回转体的载荷曲线产生波动。

（3）回转体撞水后，轴向空泡先呈现不对称形态，随着入水时间延长，空泡不对称性变弱；入水速度60 m/s下，空泡先发生表面闭合，回转体入水越快，空泡表面闭合越晚。

（4）回转体入水载荷变化与入水速度有关，入水速度越大其峰值出现越早，且震荡越明显，速度超过100 m/s时回转体会产生塑性形变，在结构设计上需考虑回转体的结构强度。

本文中，入水速度远低于声速，未来将对亚声速和超声速下的入水问题展开流固耦合数值模拟，分析流体可压缩性对入水过程的影响；同时，可进一步考虑结构的模态、材料参数和结构参数等对高速入水过程中结构动响应的影响及结构的损伤破坏情况。期望为结构物的强度设计提供参考。