0 引言

运动体从空气中穿越自由液面进入水中所引起的一系列问题统称为入水流动问题。对于刚性运动体的入水问题而言，包括入水喷溅^[1-2]、入水空泡^[3]、入水抨击^[4-6]等问题，已经有100多年的研究历史。但是对于超弹性运动体的入水问题而言，由于涉及材料大变形的流固耦合问题，具有相当的研究难度，近年来才开始有相关研究。

球体作为一种典型的运动体，由于其具有对称性且运动稳定性较好，已成为入水问题的首选研究对象。Abraham等通过数值方法分析了小球垂直入水空泡流动特性^[7]，以及阻力系数随表面张力、流动状态和入水速度的变化规律。Aristoff等研究了流体属性、球体密度以及疏水性对球体入水空泡的影响^[8-10]，获得了入水空泡形态与韦伯数、邦德数的关系，发现了空泡深闭合时间与密度无关。Johnson^[11-12]、Belden^[13]、Moxnes^[14]开展了球体入水运动参数和液体介质属性对入水跳弹的影响研究，获得了运动参数对跳弹现象的影响规律，以及跳弹现象的临界入射角度及其阻力系数。李佳川等开展了加热球体入水过程的实验和数值仿真研究^[15]，分析了球体温度对空泡形态以及空泡的形成、发展闭合溃灭的影响规律。Truscott等采用实验方法研究了旋转小球入水后的流场特性^[16-17]，发现了旋转小球入水后呈现马格努斯效应的曲线轨迹。夏维学等基于VOF多相流模型分析了低弗劳德数下高速旋转的小球入水空泡及流体动力特性^[18]。

大多数经典的入水问题研究都忽略了运动体入水过程的流固耦合作用，然而近年来由于弹性材料的发展，人们也开始逐渐开始重视入水过程中的流固耦合问题。Yan等^[19]采用实验方法研究了平板入水后的抨击压力和水动力弹性振动，发现冲击载荷作用下，水动力弹性振动使紊流气-水场的压力周期性下降。Xie等^[20]结合实验和流固耦合数值方法对复合材料船体结构入水开展了研究，通过对结构固有频率的分析，讨论了其对水动力压力和应力响应的影响。Michael等^[21]采用数值模拟和实验方法，研究了弹性球体在水面的跳跃问题，分析了弹性球体水面上跳跃的内在机理。

对于超弹性球体的研究问题，任九生等^[22]基于有限变形动力学理论，分析了受到均匀拉伸载荷作用下可压缩超弹性材料球壳的有限振动问题，得到了球壳振动的时程和相图，并且证明了可压缩超弹性球壳的振动具有拟周期性和非线性。袁学刚等研究了横向各向同性的不可压缩超弹性球表面受周期性阶梯拉伸载荷时孔洞的形成和运动问题^[23]，结果表明拉伸载荷超过某一临界值时在球体表面会形成孔洞，并且在周期性的阶跃载荷下，给出了周期振荡形成孔洞的存在条件。Peyraut等开展了受静水压力下超弹性球体的数值仿真研究^[24]，发现了牛顿-拉弗森算法在无方向保持下容易发散，且提出最优时间步载荷算法。Fanning和Hurd等采用高速摄像方法开展了可变形弹性体球入水实验^[25]，研究发现弹性球体入水后形成了独特的嵌套空泡现象，并通过材料性能和冲击条件预测了球体的变形和振动。

通过上述分析发现，目前球体入水研究主要集中于刚性球体的入水问题，同时，超弹性球体的研究主要集中于均匀或者静止载荷下超弹性球体的结构动态研究，而涉及超弹性球体入水流固耦合问题的研究相对较少。同时，虽然少数学者开展了可变形弹性体球入水实验，但是由于实验条件的限制，对于空泡轮廓的准确捕捉、流场参数的分布和应力变形的测量等问题暂时都无法很好地解决，因此结合实验开展数值模拟研究具有重要意义。本文采用实验与数值模拟相结合的方法，开展超弹性球体垂直入水空泡流动研究，分析超弹性球体入水过程的流动特性。

1 实验系统及实验方法

本文的超弹性球体入水空泡实验在室温下进行，实验系统的整体示意图如图 1所示。水箱尺寸为1.5 m×0.8 m×0.9 m，与球体相比有足够大的尺度，故可以忽略边壁效应对入水空泡形态的影响。水箱材料均为钢化玻璃，其中侧壁厚度为10 mm，底面厚度为15 mm，各面之间由酸性硅酮玻璃胶胶合。超弹性球体由释放装置在水槽中心位置的正上方无干扰地释放，以实现对超弹性球体入水速度的控制。

实验采用Photron FASTCAM SA-X型高速摄像机对超弹性球体垂直入水过程进行拍摄，拍摄帧率为1000 fps。实验时在水箱后方设置平面LED灯阵作为光源，以保证采光要求。同时在后方光源和水箱之间竖直布置柔光屏，以提供较柔和的光线，从而能够较好地捕捉流场细微结构，得到优质的图像。本文实验所捕捉的空泡形态、位移及变形量都是基于四次以上重复性实验结果获得的。

实验采用的超弹性球体是由不可压缩硅橡胶制作而成。超弹性球体的材料参数为：直径D= 60 mm，密度ρ_s=988 kg/m³，剪切模量G=12.7 kPa。

2 数值计算方法 2.1 流场模拟理论

本文采用Navier-Stocks方程描述入水空泡流动，通过VOF多相流模型捕捉运动体入水过程中的气液界面。

运动方程:

$ \frac{{\partial \left( {\rho {u_i}} \right)}}{{\partial t}} + \frac{{\partial \left( {\rho {u_i}{u_j}} \right)}}{{\partial {x_j}}} = - \frac{{\partial p}}{{\partial {x_i}}} + \mu \frac{{\partial {u_i}}}{{\partial {x_j}\partial {x_j}}} + \rho {g_i} $

(1)

连续性方程:

$ \frac{{\partial {u_i}}}{{\partial {x_i}}} = 0 $

(2)

其中，u_i为流体速度，ρ为流体密度，μ为流体动力黏度，p为流体压强，g_i为重力加速度。

$ \rho = {\rho _l}{\alpha _l} + {\rho _g}\left( {1 - {\alpha _l}} \right) $

(3)

$ \mu = {\mu _l}{\alpha _l} + {\mu _g}\left( {1 - {\alpha _l}} \right) $

(4)

其中，ρ_l为水的密度；ρ_g为空气密度；μ_l为水的动力黏度；μ_g为空气的动力黏度；α_l为液相体积分数；α_g为气相体积分数。根据VOF多相流模型：α_l=1代表流域网格单元完全为液相；α_l=0代表流域网格单元完全为气相；0 < α_l < 1代表流域网格单元将出现气液交界面。

2.2 固体结构模拟理论

基于虚功原理，采用非线性有限元法分析结构的变形和水动力荷载，必须使得物体每个单元的力和力矩得到平衡。

离散化的平衡方程可以写成:

$ {\mathit{\boldsymbol{M}}^{NM}}{{\mathit{\boldsymbol{\ddot u}}}^M} - {\mathit{\boldsymbol{P}}^N} + {\mathit{\boldsymbol{I}}^N} = 0 $

(5)

其中，P^N为外力；I^N为内力；M^NM为质量矩阵；u^M为单元加速度；M^NMü^M为由材料惯性引起的力。

内力方程可以表示为:

$ {\mathit{\boldsymbol{I}}^N} = \int\limits_V {{\mathit{\boldsymbol{\beta }}^N}} :\mathit{\boldsymbol{\sigma }}{\rm{d}}{V_0} $

(6)

其中，V₀为球体的未变形体积；σ为应力张量；β^N为弹性体本构模型定义的应变率。

外力方程可以表示为:

$ {\mathit{\boldsymbol{P}}^N} = \int\limits_S {{\mathit{\boldsymbol{N}}^N}} \cdot t{\rm{d}}S + \int\limits_S {{\mathit{\boldsymbol{N}}^N}} \cdot \mathit{\boldsymbol{F}}{\rm{d}}V $

(7)

其中N^N为单元形状函数；t为曲面上的牵引力矢量；F为球面上的力。

由于结构计算运行是隐式的，因此需要对矩阵求逆，并在每一增量中求解存在的非线性平衡方程。用牛顿法求解了非线性平衡方程的内解。

2.3 数值模型

本文基于径向基函数(RBF)的变形插值法和虚功原理的力映射法，开展超弹性球体入水流固耦合问题的数值模拟。流场计算通过在StarCCM+软件中采用重叠网格技术来实现，流体计算域尺度为45D×20D×20D，重叠域尺度约为3D×3D×3D。重叠域跟随超弹性球体在流场计算域中运动。计算域顶端为滞止入口、侧面为壁面，底端为压力入口，其压力值及液相体积分数通过自定义函数进行指定。同时，为保证计算的准确性，计算域的网格划分方式均采用六面体结构网格，在球体3D范围内进行渐进加密，加密区外围的网格稍稀疏。流场计算总体网格数量约为280万，流场网格流场计算域如图 2所示。其中，球体的初始入水速度为u₀，自由液面下方球体最低点的深度为y_b，初始直径为D，变形后直径为${\tilde D}$(其定义为通过几何中心平行于x轴与球面交点之间的距离)，定义球体变形系数为$\lambda = \tilde D/D$。

结构计算在Abaqus模拟软件中采用隐性算法来实现。在结构场中定义了弹性体材料、运动方式。为了与流场网格节点相对应，结构场计算也采用六面体结构网格，结构单元采用C3D8RH，固体结构网格约130万，固体结构网格如图 3所示。

在流固耦合计算中，SIMULIA协同仿真引擎负责Abaqus和StarCCM+之间的双向耦合代码通信，利用求解器之间的这种通信，即可以获得整个流体/固体交界面的完全求解。

3 结果及讨论 3.1 超弹性球体入水实验研究

在初始入水速度u₀=4.7 m/s条件下，针对直径D=60 mm，密度ρ_s=988 kg/m³，剪切模量G=12.7 kPa的不可压缩硅胶超弹性球体开展入水实验。入水空泡和球体变形的实验结果如图 4所示。在入水撞击阶段(t=0 ms)，超弹性球体由于经历了从空气到水的跨介质突变而引起空泡与喷溅的形成。同时，球体自身也受到强烈的冲击，由于超弹性球体刚度较小，与水面发生冲击时发生了明显的变形(t=6 ms)。随着球体入水深度增大，流动分离点上方已获得动能的流体在惯性作用下继续向外排开，导致空泡直径变大，空泡演化形成敞开趋势。此时球体变形为一个与自由液面近似平行的横向扁状球体，将水大幅度地排开，进而也形成了图中所示又宽又短的空泡(t=8 ms)。但是由于超弹性球体内部弹性力的存在，使得超弹性球体由原来近似平行于自由液面的横向椭球体反弹变形为一个纵向扁平椭球体(t=18 ms)。然后该纵向椭球体径向膨胀穿透第一个空泡壁面，同时，在水动力的作用下，纵向扁平椭球体又变形回原来的横向扁平椭球体，排开附近流体，进而出现了第二个空泡，由于第二个空泡没有大量空气进入，因此第二个空泡相对第一个空泡要小很多，第二个小空泡被称作嵌套空泡^[22](t=24~40 ms)。该嵌套空泡产生的主要原因是超弹性球体入水后发生材料震荡导致球体形态发生大变形而产生。随着超弹性球体不断下降，球体将穿透嵌套空泡，由于此时超弹性球体下降过程被空泡包裹，其所受阻力小，变形不大，小幅度地排开水，如水珠一般下降(t=66 ms)，直至最后发生空泡闭合(t=110 ms)。

3.2 超弹性球体入水数值模拟研究

在初始入水速度u₀=4.7 m/s条件下，针对直径D=60 mm，密度ρ_s=988 kg/m³, 剪切模量G=12.7 kPa的超弹性球体入水过程开展数值模拟研究。入水空泡的数值模拟结果如图 5所示。从图中可以看出，数值模拟较好地模拟出了实验中的入水空泡演化过程。在相同的入水时刻下，入水空泡形态、球体的变形行为与实验结果具有较好的一致性。

为了进一步验证数值模拟方法的有效性，选取弹性球体最低点位移作为研究对象，将位移的数值结果与实验数据进行了对比，如图 6所示。从图中可以看出球体位移的数值结果与实验数据具有较好的一致性。

为了验证数值结果中球体变形模拟的准确性，将球体变形系数λ的数值结果与实验数据进行了对比，如图 7所示。从图中可以看出二者在变形周期上稍有差异，但是在变形结构方面符合得很好，最大误差约为3.5%。

通过上述关于入水空泡形态、球体位移和变形系数的分析对比，验证了本文数值计算方法的有效性。

3.3 超弹性球体入水流动特性

为了深入研究超弹性球体入水特有的入水流动特性，在入水速度u₀=4.7 m/s的条件下，本文采用流固耦合数值模拟方法，针对直径D=60 mm，剪切模量G=12.7 kPa的超弹性球体与具有相同密度和直径、剪切模量G=337 MPa的刚性球体的入水流动特性开展对比研究。

3.3.1 超弹性球体入水空泡特性

为了定量地对比超弹性球体与刚性球体空泡的尺寸，分别定义空泡特征长度和特征宽度。其中空泡长度l取自由液面到空泡与球体接触点的距离来表征，空泡宽度分别取空泡口横向宽度D_z和纵向宽度D_y来表征，如图 8所示。

图 9给出了超弹性球体与刚性球体无量纲空泡长度的数值结果。从图中可以看出，在同一时刻，超弹性球体的空泡长度要小于刚性球体的空泡长度。其原因主要是超弹性球体在入水过程中由于变形会导致能量损耗，超弹性球体的变形能量损耗明显大于刚性球体(如图 10所示)，因此入水过程中超弹性球体速度降低幅度大于刚性球体，进而导致空泡长度缩短。

图 11给出了超弹性球体和刚性球体无量纲空泡横向宽度与纵向宽度的数值结果。从图中可以看出，超弹性球体入水空泡的横向宽度和纵向宽度均大于刚性球体的空泡宽度。其原因主要是超弹性球体入水后变形严重，从而导致流体向外排开液体的幅度越大，故入水后超弹性球体形成的空泡宽度变大。

其次，本文对比分析了某一时刻不同剪切模量条件下，多个截面处超弹性球体的入水空泡口轮廓，如图 12所示。从图中可以看出在h/D=0~3的多个不同截面处，剪切模量为G=12.7 kPa球体空泡口轮廓最大，G=25.3 kPa球体空泡口轮廓其次，G=42.0 kPa球体空泡口轮廓最小。故球体的剪切模量越小，入水空泡轮廓越大，其主要原因是剪切模量越小，入水后变形量就越大，故入水后空泡口轮廓越大。

3.3.2 超弹性球体表面压力特性

为了对比超弹性球体与刚性球体的表面压力特性，提取了xoy平面与球体表面交线的下半部分圆弧线压力，从而获得典型时刻下球体表面的压力变化规律，如图 13所示。

从图中可以看出刚性球体入水后球体表面的压力随y轴位置的变化趋势都呈现开口向下的抛物线形式，球体顶端受到的正向压力最大，随着位移向轴线两侧移动压力逐渐降低，在球体表面与空泡接触的位置处压力达到零。

从图中还可以看出超弹性球体入水后球体表面压力随y轴位置的变化在不同时刻呈现不同的变化趋势。当t=20 ms时，球体内部弹性力使得超弹性球体反弹形成纵向扁平椭球体，此时球体下表面呈现回弹收缩趋势，底部形成低压区，所以较刚性球面的压力低，且球体赤道处向内收缩变形剧烈，该变形会产生反向运动与空气尾流发生碰撞，使得弹性球体附近的流体流动受到较大的扰动，球体两侧形成局部扰流，故出现局部负压区域。当t=30 ms时，超弹性球体受到较刚性球体大2倍以上的水动压力，使得超弹性球体由纵向扁平椭球体变形为横向扁平椭球体，进而也形成了第一个嵌套空泡，同时，局部扰流导致球体赤道两侧出现两处局部负压区。当t=50 ms时, 超弹性球体受到不均匀的表面压力，此时，球体处于第二个嵌套空泡的形成阶段，逐渐由纵向扁平椭球体向横向扁平椭球体变形，同时，由于局部扰流球体赤道两侧出现两处局部负压区。同理，当t=70 ms时，超弹性球体受到局部高压，但球体两侧因湍流出现负压区。

图 14给出了超弹性球和刚性球体入水流场的压力分布云图。从图中可以看出，刚性球体入水后表面压力的变化范围约为0 ~ 7 kPa，而超弹性球体表面压力的变化范围约为-8~10 kPa。超弹性球体的表面压力变化范围要比刚性球体大，这是由于局部扰流与球体变形作用导致的。

4 结论

本文采用实验与数值模拟相结合的方法对超弹性球体低速垂直入水问题进行了研究，获得的主要结论如下：

1) 通过实验研究发现了超弹性球体入水后特有的“横向扁平球体-纵向扁平球体-横向扁平球体”的周期性变形行为，同时入水空泡的发展过程中会形成嵌入空泡。

2) 由于变形产生的能量损耗，从而造成超弹性球体入水空泡长度要小于刚性球体的入水空泡长度。同时，由于超弹性球体的变形行为，导致超弹性球体入水后形成的空泡宽度要大于刚性球体的空泡宽度。超弹性球体入水后，球体剪切模量越小，入水空泡口轮廓越大。

3) 超弹性球体不同时刻的表面压力分布与其变形情况密切相关，由于局部扰流与球体变形导致超弹性球体表面压力的变化范围较大。