NUMERICAL INVESTIGATION ON THE WATER-ENTRY CAVITY FEATURE AND FLOW STRUCTURE OF A SPINNING SPHERE BASED ON LARGE-EDDY SIMULATION
-
摘要: 圆球旋转入水过程对于基于先导物投放的新型入水降载方式具有重要研究价值. 采用大涡模拟方法结合均质多相流模型和VOF界面捕捉算法, 对低弗劳德数条件下疏水圆球高速旋转入水的自由运动过程进行了数值模拟, 研究了转速对入水空泡演化、流场结构和水动力特性的影响. 采用动网格与滑移网格技术实现圆球的自由运动, 并基于试验结果对比验证了数值模拟的可靠性与正确性. 旋转运动的升力效应导致圆球入水弹道发生偏转并从水面携带横向楔形射流侵入空泡内部. 采用入水砰击速度与转速进行归一化分析, 结果表明入水转速的增加显著改变了圆球的动力特性: 水平方向的速度和位移以及升力峰值都随入水转速的增加而变大, 但升力峰值受到入水速度的限制; 而垂直方向的速度和加速度以及空泡断裂深度几乎不受转速增加的影响, 并且空泡深闭合发生前圆球转速变化不大. 入水转速的增加也使液面飞溅环和空泡断裂的非对称性增强, 在较低转速时发生空泡表面闭合, 而在较高转速时则发生空泡深闭合. 对于空泡深闭合模式, 入水转速的增加带来更强的横向楔形射流, 并且抑制了空泡断裂产生的高压以及相应涡结构的生成, 致使圆球在入水砰击阶段承受更低的侧向压力.Abstract: The water-entry process of spinning sphere has great significance to the research of the up-to-date load reduction method of water-entry based on pre-launched object. In the present work, large-eddy simulation method is used together with the homogeneous multiphase flow model and VOF algorithm of interface capturing, to simulate the water-entry free motion of a fast-spinning sphere with hydrophobic coating at low Froude number, thus to investigate the water-entry cavity evolution, the flow structure and the hydrodynamic features. The free motion of the sphere is achieved through the dynamic mesh and sliding mesh techniques. The reliability and accuracy of the numerical simulation results are validated by comparison with previously published experimental results with good agreement on the transient cavity shape and the motion of the sphere. The spinning motion induces a lift force on the sphere and the trajectory of the sphere has significant curvature along its descent. A persistent wedge of fluid is emerged across the center of the cavity due to the fluid along the surface dragged by the sphere. The velocity and spin rate were normalized with the impact velocity and spin rate to analyze the numerical results. It shows that the spin rate has significant influence on the cavity evolution and hydrodynamic characteristics. Both of those cavity shapes have asymmetrical splash curtain and collapse asymmetrically. As spin rate increases, the horizontal velocity and the maximum lift force increase, while the maximum lift force is also limited by the impact velocity. The spin rate increase also leads to a stronger wedge of fluid forming. As a result, the pinch-off pressure maximum decreases and less vortex structures are observed. And also, the spin rate increase leads to lower side pressure during the initial impact phase. However, the vertical dynamic characteristics of spheres, like vertical velocity, acceleration and immersion depth of pinch-off, are less affected by the spin rate. Moreover, the sphere spin rate is less affected by the impact spin rate increase before cavity pinch-off.
-
Keywords:
- cavity shape /
- cavity evolution /
- water-entry trajectory /
- large eddy simulation /
- spinning sphere
-
引 言
物体穿越水面的入水问题与许多军事应用直接相关, 如: 舰载发射航行体入水[1]、船舶砰击[2]、飞机迫降[3]等. 物体入水问题的研究已经持续了近一个世纪, 这个经典问题的影响因素有很多: 物体几何形状[4-6]、冲击速度和角度[5]、物体密度和表面材质[6]等. 当引入变化的物体密度[7]和物体入水带有旋转时[8], 问题变得更加复杂. 研究圆球高速旋转入水过程的空泡特性和精细流场结构, 可为将来基于多物体入水[9]的新型航行体入水降载方式提供关键性的力学机理认识.
入水问题的关注点主要包括入水冲击和入水空泡两个方面. Von Karman[10]提出的附加质量方法为入水冲击研究提供了先驱性理论指导, 其后国际上的理论研究大多围绕该理论展开. 对于圆球入水问题, Watanabe[11]开展的相关试验为后人理论模型的发展提供了关键的验证数据. 入水空泡的研究最早始于19世纪末Worthington和Cole[12]的工作, 其采用电子闪光摄像捕捉圆球的垂直入水过程, 首次系统地阐述了水面飞溅、空泡闭合和射流等现象. Truscott等[8,13-14]对圆球高速旋转入水作了系统性的试验探究, 包括采用不同入水速度和入水转速, 球面采用不同接触角的涂层以及采用不同密度的圆球进行试验. 他们发现由于圆球表面的无滑移边界, 使得圆球旋转会在空泡内部形成横向楔形射流, 并导致圆球弹道会发生明显的横向偏转.
近年来, 高性能计算机的发展使得数值模拟成为求解入水问题越来越重要的手段. 数值模拟在湍流场细节和涡结构时空演变的精细捕捉方面, 相比试验技术有天然的优势. 而对于入水问题, 数值模拟方法有雷诺平均(RANS)类方法, 大涡模拟(LES)类方法, 也有部分直接数值模拟方法(DNS)[15], 还有一些基于拉格朗日体系的无网格方法[16-18]等.
RANS类方法是最常用的湍流模拟方法, Yu等[19]采用RANS方法模拟了圆球的入水过程, 发现无量纲化之后的冲击力与入水速度和圆球半径都无关, 只和球的质量有关. 周波等[20]使用SST模型对不同旋转方向强制运动圆球的倾斜入水过程进行数值模拟, 发现球体绕不同轴旋转倾斜入水时的运动受力特性的差异显著.
然而传统RANS方法采用的涡黏湍流模型对于空泡流这种包含大范围流动分离和涡旋流动的情形难以获得较好的计算结果. 相对而言大涡模拟类的高精度数值模拟方法(LES, DES, LES/RANS混合等) 在湍流空泡流领域的应用在近年来取得了长足的进步. Li等[21]建立了六自由度的圆盘撞击水面空泡演化的数值计算模型, 他们发现冲击生成的空泡在圆盘前进方向存在非对称与旋转的现象. Zhao等[22]对低弗劳德数下圆柱斜入水问题进行了数值模拟和试验探究, 他们利用LES与VOF (volume of fluid)相结合的方法求解湍流场和气液界面. 并且成功捕捉了多种尺度下圆柱周围的涡结构, 同时对多种入水参数进行了数值模拟.
综上所述, 目前入水问题的研究方法众多, 而旋转圆球入水还有许多复杂的现象值得深入探究, 在这一方面采用LES类方法进行模拟的相关研究较少. 本文采用LES结合动网格与滑移网格技术, 基于壁面边界层八叉树各向同性网格, 对低弗劳德数下旋转圆球入水的自由运动过程开展数值模拟, 并与试验结果[7-8]进行对比, 着重分析空泡形态演化和水中弹道随入水旋转速度的变化规律, 揭示圆球高速旋转入水过程的空泡流场精细结构, 为多物体入水降载问题提供基础力学机理的认识.
1. 数学模型
1.1 基本控制方程
本文采用均质多相流的VOF模型, 气液混合介质由各相体积分数比例混合而成, 并将这种多相混合物看作一种单一的变密度流体. 各相介质共享相同的速度、压力、湍流量等各种流动物理量, 混合物的物性参数则通过各相加权平均获得.
混合介质的流场控制方程如下
$$ \frac{{\partial \rho }}{{\partial t}}{\text{ + }}\frac{{\partial (\rho {u_j})}}{{\partial {x_j}}} = 0 $$ (1) $$ \begin{split} & \frac{{\partial (\rho {u_i})}}{{\partial t}} + \frac{{\partial (\rho {u_i}{u_j})}}{{\partial {x_j}}} =\\ &\qquad - \frac{{\partial p}}{{\partial {x_i}}} + \mu \frac{\partial }{{\partial {x_j}}}\left( {\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_i}}}} \right) + {g_i} + F_i^{st} \end{split} $$ (2) 其中ui为速度矢量, ρ为混合介质密度, 即
$ \rho = $ $ {\rho _l}{\alpha _l} + {\rho _g}{\alpha _g} $ , αl为液相体积分数, αg=1−αl为气相体积分数, μ为混合介质的动力学黏性系数, 即$ \mu = {\mu _l}{\alpha _l} + $ $ {\mu _g}{\alpha _g} $ , gi为重力加速度,$F_i^{st}$ 为单位体积的表面张力.1.2 大涡模拟滤波的控制方程
在大涡模拟(LES)中, 物理量的求解是通过滤波核函数
$G({\boldsymbol{x}},{{\boldsymbol{x}}^\prime })$ 在空间上对上述方程进行滤波来得到的. 滤波过程能很好的拆分尺度大于滤波长度的涡与尺度小于滤波长度的涡. 滤波之后的变量如下定义$$ \bar \phi ({\boldsymbol{x}}) = \int_D \phi ({{\boldsymbol{x}}^\prime })G({\boldsymbol{x}},{{\boldsymbol{x}}^\prime }){\rm{d}}{{\boldsymbol{x}}^\prime } $$ (3) 滤波后小尺度变量为
$\phi ' = \phi - \bar \phi $ ,而在有限体积法的数值求解过程中, 网格尺度V是隐式的滤波尺度, 即$G({\boldsymbol{x}},{{\boldsymbol{x}}^\prime }) = 1/V$ 通过对式(1)和式(2)进行滤波运算之后, 能得到以下滤波之后的控制方程$$ \frac{{\partial \rho }}{{\partial t}}{\text{ + }}\frac{{\partial (\rho {{\bar u}_j})}}{{\partial {x_j}}} = 0 $$ (4) $$ \begin{split} & \frac{{\partial (\rho {{\bar u}_i})}}{{\partial t}} + \frac{{\partial (\rho {{\bar u}_i}{{\bar u}_j})}}{{\partial {x_j}}} = - \frac{{\partial \bar p}}{{\partial {x_i}}} +\\ &\qquad \mu \frac{\partial }{{\partial {x_j}}}\left( {\frac{{\partial {{\bar u}_i}}}{{\partial {x_j}}} + \frac{{\partial {{\bar u}_j}}}{{\partial {x_i}}}} \right) + {g_i} + F_i^{st} - \frac{{\partial {\tau _{ij}}}}{{\partial {x_j}}} \end{split} $$ (5) 其中τij为亚格子应力, 定义如下
$$ {\tau _{ij}} \equiv \overline {{u_i}{u_j}} - {\bar u_i}{\bar u_j} $$ (6) 1.3 亚格子尺度应力模型(SGS)
对N-S方程滤波之后的亚格子应力项是未知的, 所以需要对其进行建模来求解. 采用Boussinesq 假定[23], 亚格子应力计算方法如下
$$ {\tau _{ij}} - \frac{1}{3}{\tau _{kk}}{\delta _{ij}} = - 2{\mu _t}{\bar S_{ij}} $$ (7) 其中μt 为亚格子涡黏系数. 采用常规滤波方法, 所以各向同性部分τkk可以与滤波之后的压力项合并 [24],
${\bar S_{ij}}$ 为滤波之后的应变率张量定义如下$$ {\bar S_{ij}} = \frac{1}{2}\left( {\frac{{\partial {{\bar u}_i}}}{{\partial {x_j}}} + \frac{{\partial {{\bar u}_j}}}{{\partial {x_i}}}} \right) $$ (8) 而对于亚格子涡黏系数μt, 采用Nicoud与Ducros[25]提出的WALE (wall-adapting local eddy-viscosity)模型. 该模型优势在于剪切层流的μt=0, 相比传统的Smagorinsky-Lilly 模型[26]能较好地返回流体区域中层流的特性, 并且在壁面拥有更好的对数律y3 性质. 该模型具体形式如下
$$ {\mu _t} = \rho L_s^2\frac{{{{\left( {S_{ij}^dS_{ij}^d} \right)}^{3/2}}}}{{{{\left( {{{\bar S}_{ij}}{{\bar S}_{ij}}} \right)}^{5/2}} + {{\left( {S_{ij}^dS_{ij}^d} \right)}^{5/4}}}} $$ (9) 其中在WALE模型中的亚格子尺度混合长度Ls与
$ S_{ij}^d $ 定义如下$$ \;\;\;\;{L_s} = \min \left( {\kappa d,{C_w}{V^{1/3}}} \right) $$ (10) $$ \left.\begin{array}{l} S_{ij}^d = \dfrac{1}{2}\left( {\bar g_{ij}^2 + \bar g_{ji}^2} \right) - \dfrac{1}{3}{\delta _{ij}}\bar g_{kk}^2\\ {\bar g_{ij}} = \dfrac{{\partial {{\bar u}_i}}}{{\partial {x_j}}} \end{array} \right\}$$ (11) 其中冯卡门常数κ =0.4187, d为单元中心到壁面距离, ΔV为单元体积. WALE常数采用Cw=0.325, 这个值对于大多数流动拥有更好的适应性.
1.4 表面张力模型
当气液固三相接触时, 在三相交界处产生接触线与接触角. 当温度和压强给定, 在平衡状态下这个接触角被称为静态接触角. 然而通常在试验中能观察到移动的接触线, 这个接触线对应的接触角为动态接触角(如图1所示). 而在入水过程中动态接触角对空泡的产生与形态有很重要的影响[27].
针对这个问题, 本文采用 Brackbill等[28]提出的CSS (continuum surface stress)表面张力模型描述球体表面的润湿性, 将表面张力作为VOF动量方程的源项, 具体模型如下
$$ {\boldsymbol{T}} = \sigma ({\boldsymbol{I}} - \hat {\boldsymbol{n}} \otimes \hat {\boldsymbol{n}})|{\boldsymbol{n}}| $$ (12) $$ {\boldsymbol{ n}} = \nabla \alpha $$ (13) $$ \hat {\boldsymbol{n }}= \frac{{{\boldsymbol{ n}}}}{{\left| {{\boldsymbol{ n}}} \right|}} $$ (14) 其中T为毛细应力矩阵, I为单位张量, σ 为表面张力系数,
$ \otimes $ 代表张量积, α 代表体积分数,${\boldsymbol{ n}}$ 体积分数梯度, 得到的表面张力为$$ {F^{st}} = \nabla \cdot {\boldsymbol{T}} $$ (15) CSS模型避免显式计算界面曲率, 将表面张力影响作为各向异性的毛细应力来计算, 从而能对界面上的尖角得到较好的模拟结果. 而对于气液界面与球体壁面的接触线, 也引入了壁面依附约束模型 (wall adhesion) 来处理. 将流体与壁面的接触角用于计算距离壁面的1个流体单元界面法向的条件, 如下所示
$$ \hat {\boldsymbol{n}} = {\hat {\boldsymbol{n}}_w}\cos \left( {{\text{π}} - {\theta _w}} \right) + {\hat {\boldsymbol{t}}_w}\sin \left( {{\text{π}} - {\theta _w}} \right) $$ (16) 其中θw为壁面依附角,
${\hat {\boldsymbol{n}}_w}$ 为壁面法向单位矢量,${\hat{\boldsymbol{ t}}_w}$ 为壁面切向单位矢量.2. 数值算法与验证
本文采用有限体积法对流动控制方程进行时空离散, 同时采用基于压力的分离式求解器(SIMPLE)对控制方程中耦合的非线性项进行解耦迭代求解. LES需要解析流场中大部分高波数的湍流脉动, 从而需要高精度的数值格式以及高精度的网格, 而偶阶次的数值格式天生不具有耗散性, 很好的规避了对SGS模型在耗散上的影响. 所以采用中心差分格式能很好的适应LES方法, 但是中心差分格式会对流场引入数值震荡, 所以最终我们选择有界中心差分格式[29]作为动量的离散格式. 而对于方程中变量的梯度项, 采用最小二乘格式来进行计算. 由于入水过程为低速所以忽略空化效应的影响. 本文采用了VOF多相流模型对气液界面进行捕捉, 而迎风格式对耗散会过估计, 中心差分则有数值震荡都无法很好的对界面进行捕捉. 为了克服这些问题, 采用一种高分辨率界面捕捉算法(HRIC)进行空泡面捕捉[30], 相比QUICK与其他二阶格式拥有更高的计算精度, 相比传统的几何重构格式计算速度能更快. 圆球的三自由度运动采用动网格方法中的动框架来实现计算域平移运动控制, 圆球带空泡在水中平移和旋转运动采用相对于背景计算域的滑移网格来控制.
2.1 计算域及边界条件
为了校验当前数值模拟方案的可靠性, 与Aristoff等[3]的垂直入水试验进行了对比. 球直径为D=0.0254 m, 半径为R, 分别对比了钢和聚丙烯两种不同密度比m*=ρl /ρs=7.86, 0.86的情况. 如图2 所示, 整个计算域采用球直径R进行标准化, 垂直长度为100R, 水平长度为40R, 宽度为20R, 球心距离顶部压力边界设置为80R, 距离侧面速度边界为20R, 在宽度方向选取了对称面, 同时在球面上设置为无滑移边界. 水面位置距离球为3R. 坐标轴原点同水面保持一致如图3 所示, 其中Vl代表圆球入水时左侧壁面速度大小, Vr代表圆球入水时右侧壁面速度大小, 圆球各项初始化参数如表1所示.
Parameters Symbol Value impact velocity Vo 2.17 m/s impact spin rate ωo 0 rad/s diameter D 0.0254 m water density ρl 998.2 kg/m3 density ratio m* 7.86, 0.86 surface tension coefficient σ 0.0717 N/m contact angle θw 160° gravity acceleration g 9.81 m/s2 Froude number Fr ${V_o}/\sqrt {gD} $ non-dimensional spin parameter S ωoR/Vo 具体计算网格采用结构化与非结构化相结合的方式, 尽可能采用结构化网格以减小数值误差(如图4所示). 为了实现精度较高的流场大涡模拟, 近壁面边界层附近采用尽量各向同性和均匀的网格, 以减小大涡模拟低通滤波算法中滤波函数的交换误差. 计算网格采用尽量合理的网格疏密分布. 在物体表面布置致密网格, 网格的精细度达到能够刻画湍流边界层结构和典型湍流拟序结构的程度. 除了满足数值模拟对网格质量的一般性要求, 即保证无量纲壁面垂向距离y + <1之外, 还要尽量使壁面流向和展向的无量纲网格单元大小满足Δx + ≤100和Δz + ≤30. 对于当前模拟工况Δx + , Δz + 列在表2中, 选取的时间为圆球刚刚沉入水面以下的时刻(t=10 ms 左右).
表 2 流向与展向平均无量纲网格尺度表(t=10 ms)Table 2. Averaged nondimensional grid sizes of the finest mesh in the streamwise and spanwise (t=10 ms)Computational cases Δx + Δz + Vo=5.45 m/s, ωo=199 rad/s 12 11 Vo=5.45 m/s, ωo=150 rad/s 25 25 Vo=5.45 m/s, ωo=100 rad/s 25 25 Vo=5.45 m/s, ωo=50 rad/s 12 11 而另一种普遍采用的衡量LES解析度的方式是Benard等[31]建议整体网格量能计算出流动区域内至少80%的湍动能. 将未解析湍动能比率定义为
$$ {Q_c} = \frac{{{E_{sgs}}}}{{{E_{sgs}} + {E_R}}} \leqslant 0.2 $$ (17) 其中
${E_R} = \overline {{u_i}^\prime {u_i}^\prime }/2$ 为解析的湍动能, Esgs为亚格子湍动能, 可以通过亚格子运动黏性系数νsgs来估计$$ {E_{sgs}} = C{\left( {\frac{{{\nu _{sgs}}}}{{{L_s}}}} \right)^2} $$ (18) 其中参数C选取为100.
虽然湍流脉动速度
${u_i}^\prime $ 和湍动能ER在LES中不是显式求解获得的物理量, 但是对于那些宏观时均流动几乎定常的流场, 可以通过数据采样统计的方法获得${u_i}^\prime $ 和ER. 虽然圆球入水空泡流动并不是宏观定常流场, 但是可以采用同样的网格对圆球在水中重力场中无空泡垂直降落过程这样类似的时均定常流场进行LES模拟, 进行数据采样统计从而获得脉动速度并计算湍动能, 以检验所采用的网格分辨率的可靠性. 所以模拟了圆球初始速度Vo=5.45 m/s不带空泡在水中下沉, 其中t∈[0.07, 0.18] (1100时间步数据采样), 开启数据采样, 得到均方根(RMS)的$ \sqrt {\overline {{{\left( {u'} \right)}^2}} } $ ,$ \sqrt {\overline {{{\left( {v'} \right)}^2}} } $ 与$ \sqrt {\overline {{{\left( {w'} \right)}^2}} } $ 来计算ER与Esgs.图5展示了对称面上Qc与ER大小, 可以发现在ER较高的圆球尾流区Qc远远小于0.2. 所以可以认为在流场中绝大多数区域Qc<0.2的要求都能被满足.
根据经典湍流理论, 湍流场能量在惯性子区满足−5/3次方律(Kolmogorov的K41理论)[32]
$$ E\left(k\right)={C}_{K}{\varepsilon }^{2/3}{k}^{-5/3} $$ (19) 而湍流压力脉动的能谱满足文献[33]的−7/3次方律
$$ {E}_{pp}\left(k\right)={C}_{P}{\varepsilon}^{4/3}{k}^{-7/3} $$ (20) 这也是LES研究论文中经常采用的验证方法. 为了进一步验证所采用的网格分辨率的可靠性, 实时监测获取稳定之后流场中的压力测点(在随球运动的动坐标系内, 测点坐标x = 0, y = 3R, z = 0)的压力信息, 采用快速傅里叶变换得到归一化的功率谱密度估计PSD (见图6). 其中采样频率Fs = 104, NFFT = 1024. 可以直接测量如下的双对数坐标PSD (power spectral density)图中的峰值包络线斜率, 发现与上述经典湍流理论中的−7/3次方律非常接近.
2.2 垂直入水与试验对比验证结果
图7给出了圆球垂直入水数值模拟结果与文献[7]的试验结果的对比. 从图中可以看出, 当前模拟方案在总体上对圆球深度有些低估, 但总体上模拟与试验结果是非常接近的. 表3选择了几个特定时间, 对比了圆球模拟深度与试验深度, 可以发现相对误差不大于10%, 总体来说误差水平是可以接受的.
图8展示了数值模拟和试验泡型的对比. 此处初始时刻定义同文献[7]维持一致, 即将圆球球心通过水面位置作为初始时刻. 虽然在图8(a) 68.9 ms与图8 (b) 62.2 ms接近空泡断裂时上半部分空泡略有偏小. 但总体对于定性的泡形, 当前的模拟方案也能和试验吻合较好. 特别是在空泡断裂之后(图8 (a) 75.9 ms与图8 (b) 68.7 ms), 还能捕捉到空泡面的收缩与向上的“Worthington”射流.
表 3 垂直入水圆球深度试验[7]与模拟对比表(Vo=2.17 m/s)Table 3. Comparison between the experimental[7] and simulation vertical depths for the case at Vo=2.17 m/sTime/ms Material Experiment depth/cm Simulation depth/cm Absolute error/cm Relative error/% 40.2 polypropylene −5.73 −5.31 −0.42 7.33 50.1 −6.62 −6.14 −0.48 7.25 60.2 −7.28 −6.76 −0.52 7.14 40.3 steel −9.61 −8.90 −0.71 7.39 50.3 −11.98 −11.14 −0.84 7.01 60.3 −14.32 −13.38 −0.94 6.56 3. 不同转速对圆球旋转入水影响
本文主要考察相同入水速度Vo=5.45 m/s下不同入水转速ωo=50, 100, 150, 199 rad/s对于入水空泡演化、圆球弹道、流场各个物性参数的影响. 所有初始时刻t=0 ms选取在圆球刚刚接触自由液面的时刻.
3.1 入水角速度对圆球运动学性质影响
本节将详细讨论圆球高速旋转入水的水动力学特性, 主要关注在圆球的弹道(图9)、速度(图10)、加速度(图11)以及圆球受力情况(图12)随时间t的变化.
图9展示了当圆球入水速度为Vo=5.45 m/s 时, 不同入水旋转角速度ωo=50, 100, 150, 199 rad/s 对于圆球弹道的影响. 其中为了和文献[8]的试验结果进行对比, 使用密度比m*=1.74, 半径R=28.6 mm的圆球进行模拟. 可以发现ω0=199 rad/s的模拟结果与文献[8]的试验结果吻合较好. 圆球入水后在马格努斯力的作用下弹道发生偏转, 在入水的初始阶段不同转速ωo对弹道的影响并不是很大, 但是随着圆球进入水面以下(t=10 ms之后, 见图13), 入水转速越大的球, 相应的弹道偏转也越大. 当空泡发生深闭合时, 空泡断裂时圆球所在深度都非常近似.
为了进一步探究弹道偏转与入水转速的关系, 在图10中展示了圆球质心速度与角速度随时间变化关系图, 其中采用入水速度Vo与角速度ωo对速度和角速度进行归一化. 通过比较归一化之后的结果可以发现, 入水转速的改变对水平方向速度vx/Vo的影响很大(见图10 (a)). 相应的转速越大, 水平方向速度大小也越大, 从而使得弹道偏转也越大. 对于ωo = 50 rad/s, vx/Vo近似于线性下降, 而对于ωo = 100, 150, 199 rad/s 都随时间t, 先快速下降(t = 10~20 ms), 然后趋于常数(t > 60 ms). 通过对比x方向加速度ax随时间变化(见图11(a))可以发现, 入水砰击阶段(t < 5 ms) ωo = 50 rad/s的加速度峰值不到 ωo = 100, 150, 199 rad/s 的30%, 并且当砰击阶段结束, ax也趋于平稳, 这就导致了vx/Vo是近似的线性下降. 而对于垂直方向速度vy/Vo, 并不受到ωo改变的影响(见图10 (b)), 所有工况下都呈现在入水砰击阶段(t < 5 ms), 下降较快, 之后下降趋势逐渐放缓的情况. 而圆球转速ω/ωo随时间t变化总体不大(见图10 (c), min (ω/ωo) = 0.94).
在图11中展示了圆球加速度与角加速度随时间t的变化, 在入水砰击阶段 (见图11(a)和图11(b) t = 0~4 ms), y方向加速度ay到达峰值之后快速下降, 然后趋于不变并且不同ωo对于 ay影响不大. 但是对于ωo = 199, 150, 100 rad/s时, ax继续变大直到圆球大部分进入自由液面以下(t = 6 ms)而到达峰值. 峰值的大小随ωo变大而变大, 而ωo = 199, 150 rad/s 两个情况峰值大小非常接近, 说明入水砰击阶段ax峰值大小受到Vo大小的限制. 在t = 20 ms之后ax趋于稳定直到空泡断裂(t > 100 ms). 圆球角加速度 β (见图11 (c)), 在砰击阶段震荡过后都基本维持不变, β的大小也是随着ωo的变大而变大. 结合圆球表面未沾湿面积之比Aair/Asphere随时间t变化(见图14), 可以发现ωo越大, 未沾湿比例越小, 水依附在球面的面积越大, 从而导致球面上受到黏性力也而越大, 所以 β 随着 ωo的变大而变大. 过了入水砰击阶段(t > 5 ms) Aair/Asphere基本维持不变, 使得β维持不变.
基于上述运动学的讨论, 圆球质心的受力可以做如下分析
$$ {m_s}{a_x} = {F_{hx}} $$ (21) $$ {m_s}{a_y} = - {m_s}g + {F_{hy}} $$ (22) 其中ms为小球质量, Fhx, Fhy为球面上所受黏性力与压力的合力, 球面上还应有表面张力Fst对球质心的作用, 根据式(15)表面张力大小的最大值不会超过Fst=πDσ, 其大小不到重力1%: 6σ/(D2ρsg) ≈ 10−3, ρs为小球密度. 所以为了简化分析, 忽略了表面张力对球受力的影响[15,34].
圆球所受升力Fl, 阻力Fd的方向定义如图15 所示, 故根据式(21)和式(22)可以将升阻力化简为如下形式
$$ {F_l} = {F_{hx}}\cos \gamma + \left( {{F_{hy}} - {m_s}g} \right)\sin \gamma $$ (23) $$ {F_d} = - {F_{hx}}\sin \gamma + \left( {{F_{hy}} - {m_s}g} \right)\cos \gamma $$ (24) 其中弹道切线角度γ=arctan(−vx/vy), 从而得到相应的升阻力系数
$$ \left. \begin{array}{l}{C_d} = \dfrac{{{2F_d}}}{{{\rho _l}V_o^2A}}\\ {C_l} = \dfrac{{{2F_l}}}{{{\rho _l}V_o^2A}} \end{array} \right\}$$ (25) 其中A=0.5πR2为圆球截面面积.
根据上述定义在图12中展示了圆球所受升力系数Cl, 阻力系数Cd随时间t变化. 可以发现Cl与Cd的变化趋势分别与ax , ay比较接近. Cd在度过入水砰击阶段之后(t>5 ms)从Cd=0.2开始随时间线性下降直到空泡断裂, 并且ωo对其影响并不是很明显. 但是对于升力系数Cl, ωo越大入水初始阶段(t=0~10 ms)的峰值越大.
3.2 入水角速度对空泡形态的影响
本节展示了ωo = 199, 150, 100, 50 rad/s 的空泡泡形(图13), 对称面上压力、速度云图, 以及Q准则等值面, 其中Q准则为速度梯度张量的第二不变量定义如下[35].
$$ {{Q}} = \frac{1}{2}\left( {{{\left\| {\boldsymbol{\varOmega }} \right\|}^2} - {{\left\| {{{\boldsymbol{S}}}} \right\|}^2}} \right) = \frac{1}{2}\left( {{{\bar \varOmega }_{ij}}{{\bar \varOmega }_{ji}} - {{\bar S}_{ij}}{{\bar S}_{ji}}} \right) $$ (26) $$ {\bar \varOmega _{ij}} = - \frac{1}{2}\frac{{\partial {{\bar u}_i}}}{{\partial {x_j}}}\frac{{\partial {{\bar u}_j}}}{{\partial {x_i}}} $$ (27) 图13展示了空泡随时间的演化序列, 其中球面上为水气混合密度的云图, 即不同颜色的接触线代表了圆球的水气接触线, 在图13(a) 中定性比较了文献[8]的试验与当前的模拟结果. 可以发现模拟泡形相比试验结果略小一些, 这在空泡断裂之后较为明显(图13 (a)实验111 ms与模拟110.8 ms). 这可能是由于随着圆球的下沉网格逐渐稀疏导致的. 同时在试验中飞溅冠在液面以上先行闭合, 而在当前模拟中无法很好的模拟水面以上的飞溅冠使得这一现象无法观察到. 在图13和图16可以看到有横向楔形劈尖在入水的初始阶段(t =10 ms附近)开始生成侵入空泡内部, 结合图11中ax随t变化可以发现, 横向劈尖在球面上的攀升是ax变大的主要原因. 当t>20 ms时, 横向劈尖在球面上不再继续延展 (见图13(a)~图13(c)) ax也趋于稳定. 随着圆球的下沉, 横向劈尖逐渐发展直到击穿泡面(ωo =199, 150 rad/s). 而相同时刻, ωo越大横向劈尖也越大(见图16). 而对于ωo = 100 rad/s, 上半部分楔形劈尖还未到达另一边空泡面空泡已经发生断裂.
对于ωo = 50 rad/s 的情况(图13, 图16 (d) ), 能观察到空泡在自由液面上发生表面闭合, 也没有明显的横向劈尖产生(图16 (d)), 这是由于ωo =50 rad/s对应的无量纲旋转参数S = 0.262<1, 对应Vr<0, 也就是接触线大部分都在圆球赤道以下, 所以没有横向劈尖的产生, 与文献[8]的结论保持一致.
入水初始阶段(t=10 ms附近)与空泡断裂后(t = 110 ms附近)不同ωo的压力云图展示在图17中, 其中实线代表混合密度ρm = 400 kg/m3自由液面, 可以看出由于马格努斯效应圆球一侧形成低压区并且随着ωo变大, 低压区也在变大, 这也使得弹道的横向位移也越大. 同时ωo的变化对于空泡断裂后的压力场有显著的影响. 对于ωo = 100 rad/s (见图17(c)), 在空泡的断裂位置有明显的高压区, 但是对于ωo = 199 rad/s (见图17(a)), 并不能观察到明显的高压区. 这可以如下解释: 横向劈尖侵入泡内从而影响空泡深闭合的断裂, 而ωo越大横向楔形越大, 从而使得水进入泡内更快. 最终结果就是空泡因静水压夹断时, 产生的高压区更小乃至几乎没有.
图18展示了ωo = 199, 150, 100 rad/s 在100 ms附近时的对称面速度云图与矢量场, 对比图13可以发现虽然此时空泡还未断裂, 但在对称面上由于横向楔形劈尖的存在, 已经接近断裂. 观察速度云图可以发现ωo=100 rad/s 在自由液面上没有出现明显的横向气流射出与ωo=199, 150 rad/s有所不同, 这可以被解释为ωo较小, 横向劈尖速度小, 使得自由液面变形速度较慢从而抑制了横向气流的生成.
为了探究涡结构对空泡泡形的影响, 图19展示了Q = 100000 s−2的等值面, 可以发现涡结构生成的时间点有两个: 入水初始阶段(t = 11 ms附近)与空泡断裂阶段(t = 90~111 ms). 这一点与垂直入水的情况保持一致, 但横向劈尖的存在使得空泡上半部分始终能观察到相应的涡结构. ωo的增加使入水初始阶段的涡结构更加趋向于非对称, 但是到了空泡断裂阶段, ωo越大涡结构越少. 这可以被解释为小球横向位移的增加, 使得泡内空气与外界交换受到阻碍, 抑制了涡结构的生成.
4. 结论
本文采用大涡模拟对圆球高速旋转入水建立了数值模型. 研究了圆球在相同入水速度
$V_o $ =5.45 m/s下不同入水转速ωo=50, 100, 150, 199 rad/s的空泡演化与水动力学特性, 主要结论如下.(1)数值模拟结果通过与垂直入水, 高速旋转入水试验对比, 在圆球弹道, 空泡泡形上与试验结果有较好的一致性, 验证了数值模型的正确性与可靠性.
(2)入水转速的改变对圆球的水动力学特性有显著的影响. 随着入水转速的增加, 圆球弹道水平位移更大, 圆球水平方向获得的速度更大, 在入水砰击阶段圆球所获得的升力也更大, 但升力的峰值受到入水速度的限制. 而入水转速对于圆球垂直方向特性影响不大, 包括垂直方向速度, 加速度, 空泡断裂时圆球所在深度.
(3)入水转速对于空泡形态与流场特性得到详细的探究. 在模拟中观察到ωo=50 rad/s 的表面闭合ωo =100, 150, 199 rad/s的深闭合. 对于深闭合, 入水转速的增加带来了更强的横向楔形射流并且抑制了空泡断裂带来的高压以及相应涡结构的生成.
-
Parameters Symbol Value impact velocity Vo 2.17 m/s impact spin rate ωo 0 rad/s diameter D 0.0254 m water density ρl 998.2 kg/m3 density ratio m* 7.86, 0.86 surface tension coefficient σ 0.0717 N/m contact angle θw 160° gravity acceleration g 9.81 m/s2 Froude number Fr ${V_o}/\sqrt {gD} $ non-dimensional spin parameter S ωoR/Vo 表 2 流向与展向平均无量纲网格尺度表(t=10 ms)
Table 2 Averaged nondimensional grid sizes of the finest mesh in the streamwise and spanwise (t=10 ms)
Computational cases Δx + Δz + Vo=5.45 m/s, ωo=199 rad/s 12 11 Vo=5.45 m/s, ωo=150 rad/s 25 25 Vo=5.45 m/s, ωo=100 rad/s 25 25 Vo=5.45 m/s, ωo=50 rad/s 12 11 表 3 垂直入水圆球深度试验[7]与模拟对比表(Vo=2.17 m/s)
Table 3 Comparison between the experimental[7] and simulation vertical depths for the case at Vo=2.17 m/s
Time/ms Material Experiment depth/cm Simulation depth/cm Absolute error/cm Relative error/% 40.2 polypropylene −5.73 −5.31 −0.42 7.33 50.1 −6.62 −6.14 −0.48 7.25 60.2 −7.28 −6.76 −0.52 7.14 40.3 steel −9.61 −8.90 −0.71 7.39 50.3 −11.98 −11.14 −0.84 7.01 60.3 −14.32 −13.38 −0.94 6.56 -
[1] Shi Y, Pan G, Yan GX, et al. Numerical study on the cavity characteristics and impact loads of auv water entry. Applied Ocean Research, 2019, 89: 44-58 doi: 10.1016/j.apor.2019.05.012
[2] Xie H, Liu F, Tang H, et al. Numerical study on the dynamic response of a truncated ship-hull structure under asymmetrical slamming. Marine Structures, 2020, 72: 102767
[3] Aristiff I , Grizzi S. Cavitation and ventilation modalities during ditching. Physics of Fluids, 2019, 31(5): 052101
[4] Louf JF, Chang B, Eshraghi J, et al. Cavity ripple dynamics after pinch-off. Journal of Fluid Mechanics, 2018, 850: 611-623 doi: 10.1017/jfm.2018.459
[5] Yang C, Zhang HX. Improved moving particle semi-implicit method with large eddy simulation for determining water-entry impact and damage to flat-bottomed structures. Journal of Ship Mechanics, 2019, 23(9): 1070-1085
[6] 路中磊, 孙铁志, 魏英杰等. 开放空腔壳体倾斜入水运动特性试验研究. 力学学报, 2018, 50(2): 263-273 (Lu Zhonglei, Sun Tiezhi, Wei Yingjie, et al. Experimental investigation on the motion feature of inclined water-entry of a semi-closed cylinder. Chinese Journal of Theoretical and Applied Mechanics, 2018, 50(2): 263-273 (in Chinese) doi: 10.6052/0459-1879-17-191 [7] 黄超, 翁翕, 刘谋斌. 超疏水小球低速入水空泡研究. 力学学报, 2019, 51(1): 36-45 (Huang Chao, Wen Xi, Liu Moubin. Study on low-speed water entry of super-hydrophobic small spheres. Chinese Journal of Theoretical and Applied Mechanics, 2019, 51(1): 36-45 (in Chinese) doi: 10.6052/0459-1879-18-310 [8] Truscott TT, Techet AH. Water entry of spinning spheres. Journal of Fluid Mechanics, 2009, 625: 135-165 doi: 10.1017/S0022112008005533
[9] 卢佳兴, 魏英杰, 王聪等. 圆柱体并联入水过程空泡演化特性实验研究. 力学学报, 2019, 51(2): 450-461 (Lu Jiaxing, Wei Yingjie, Wang Cong, et al. Experimental study on cavity evolution characteristics in the water-entry process of parallel cylinders. Chinese Journal of Theoretical and Applied Mechanics, 2019, 51(2): 450-461 (in Chinese) doi: 10.6052/0459-1879-18-288 [10] Von Karman T. The Impact on Seaplane Floats During Landing, Washington DC: National Advisory Committee for Aeronautics, 1929: 1-8
[11] Watanabe S. Resistance of impact on water surface, part V-sphere. Institute of Physical and Chemical Research, 1930, 23(484): 202-208
[12] Worthington AM, Cole RSV. Impact with a liquid surface, studied by the aid of instantaneous photography. Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character, 1897, 189: 137-148
[13] Techet AH, Truscott TT. Water entry of spinning hydrophobic and hydrophilic spheres. Journal of Fluids and Structures, 2011, 27(5-6): 716-726 doi: 10.1016/j.jfluidstructs.2011.03.014
[14] Truscott TT, Techet AH. A spin on cavity formation during water entry of hydrophobic and hydrophilic spheres. Physics of Fluids, 2009, 21(12): 121703
[15] Li D, Zhao X, Kong D, et al. Numerical investigation of the water entry of a hydrophobic sphere with spin. International Journal of Multiphase Flow, 2020, 126: 103234 doi: 10.1016/j.ijmultiphaseflow.2020.103234
[16] Molteni D, Colagrossi A. A simple procedure to improve the pressure evaluation in hydrodynamic context using the SPH. Computer Physics Communications, 2009, 180(6): 861-872 doi: 10.1016/j.cpc.2008.12.004
[17] Adami S, Hu XY, Adams NA. A generalized wall boundary condition for smoothed particle hydrodynamics. Journal of Computational Physics, 2012, 231(21): 7057-7075 doi: 10.1016/j.jcp.2012.05.005
[18] Fourtakas G, Jose MD, Vacondio R, et al. Local uniform stencil (lust) boundary condition for arbitrary 3-D boundaries in parallel smoothed particle hydrodynamics (SPH) models. Computers & Fluids, 2019, 190: 346-361
[19] Yu P, Shen C, Zhen C, et al. Parametric study on the free-fall water entry of a sphere by using the rans method. Journal of Marine Science and Engineering, 2019, 7(5): 122
[20] 周波, 刘辉, 王杰等. 旋转球倾斜入水空泡演变和运动特性数值模拟. 华中科技大学学报(自然科学版), 2020, 48(10): 92-98 (Zhou Bo, Liu hui, Wang Jie, et al. Numerical simulation on cavity evolution and motion characteristics for oblique water entry of a rotation sphere. Journal of Huazhong Universiy of Science and Technology (Natural Sciences Edition) , 2020, 48(10): 92-98 (in Chinese) [21] Li C, Wang C, Wei Y, et al. Three-dimensional numerical simulation of cavity dynamics of a stone with different spinning velocities. International Journal of Multiphase Flow, 2020, 129: 103339
[22] Zhou H, Sun T, Quan X, et al. Large eddy simulation and experimental investigation on the cavity dynamics and vortex evolution for oblique water entry of a cylinder. Applied Ocean Research, 2018, 81: 76-92 doi: 10.1016/j.apor.2018.10.008
[23] Hinze JO. Turbulence. 2nd ed. New York: McGraw-Hill, 1975: 92-96
[24] Erlebacher G, Hussaini MY, Speziale CG, et al. Toward the large-eddy simulation of compressible turbulent flows. Journal of Fluid Mechanics, 1992, 238: 155-185 doi: 10.1017/S0022112092001678
[25] Nicoud F, Ducros F. Subgrid-scale stress modelling based on the square of the velocity gradient tensor. Flow Turbulence and Combustion, 1999, 62(3): 183-200 doi: 10.1023/A:1009995426001
[26] Smagorinsky J. General circulation experiments with the primitive equations I. The basic experiment. Mon. Weather Rev., 1963, 91: 99-164 doi: 10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2
[27] Duez C, Ybert C, Clanet C, et al. Making a splash with water repellency. Nature Physics, 2007, 3(3): 180-183 doi: 10.1038/nphys545
[28] Brackbill JU, Kothe DB, Zemach C. A continuum method for modeling surface tension. Journal of Computational Physics, 1992, 100(2): 335-354 doi: 10.1016/0021-9991(92)90240-Y
[29] Leonard BP. The ultimate conservative difference scheme applied to unsteady one-dimensional advection. Computer Methods in Applied Mechanics and Engineering, 1991, 88(1): 17-74 doi: 10.1016/0045-7825(91)90232-U
[30] Muzaferija S, Peric M, Sames P, et al. A two-fluid navier-stokes solver to simulate water entry//Proc. 22nd Symposium on Naval Hydrodynamics, Washington, 1998: 277–289
[31] Benard P, Balarac G, Moureau V, et al. Mesh adaptation for large-eddy simulations in complex geometries. International Journal for Numerical Methods in Fluids, 2016, 81(12): 719-740 doi: 10.1002/fld.4204
[32] Lesieur M. Turbulence in Fluids: Fourth Revised and Enlarged, Edition 4. Dordrecht: Springer Netherlands, 2008: 197-199
[33] Batchelor G K. Pressure fluctuations in isotropic turbulence. Math. Proc. Camb. Phil. Soc. 1951, 47(2): 359-374.
[34] Mansoor MM, Vakarelski IU, Marston JO, et al. Stable-streamlined and helical cavities following the impact of leidenfrost spheres. Journal of Fluid Mechanics, 2017, 823: 716-754 doi: 10.1017/jfm.2017.337
[35] Hunt JCR, Wray AA, Moin P. Eddies, stream, and convergence zones in turbulent flows. in: Studying Turbulence Using Numerical Simulation Databases-II, Center for Turbulence Research 2, 1988: 193-208
-
期刊类型引用(4)
1. 张兴权,刘航轩,段士伟,裴善报,左立生,张晖. 激光轰击金属液面溅射过程的数值模拟. 力学学报. 2024(01): 285-297 . 本站查看
2. 王连安,徐艳,王尊策,李森,张井龙,刘海水. 风琴管喷嘴空化水射流流场的大涡模拟. 机械科学与技术. 2024(09): 1514-1521 . 百度学术
3. 石晓,蒋勤,钟振宇. 波浪对水平板冲击作用水气二相流数值模拟研究. 力学学报. 2024(10): 2878-2887 . 本站查看
4. 祁晓斌,施瑶,刘喜燕,潘光. 阶梯式圆柱射弹小角度入水弹道特性研究. 力学学报. 2023(11): 2468-2479 . 本站查看
其他类型引用(0)