在海上进行反水雷作战时，往往需要短时间内投放多枚射弹，使敌方防御系统短期内应接不暇.这是典型的运动体并联入水过程.该过程中，两个或多个运动体同时或在极短时间间隔内平行地穿越自由液面进入水中，由于各运动体流场和弹道存在耦合作用，力学环境更加复杂，运动体极易发生运动失稳.

国内外对单个运动体入水开展了大量的研究，包括空泡形态特性、运动体受载特性和弹道特性. Logvinovich^[1]提出了空泡独立膨胀原理，为后续的空泡形态研究提供了理论基础；Lee等^[2]提出了高速入水空泡形态发展的理论分析方法，该方法可以模拟高速入水过程中空泡的生成及溃灭；Aristoff等^[3]基于空泡压力平衡理论和势流理论提出了空泡形态的理论求解方法；王京华等^[4]基于空泡独立膨胀原理，研究了空泡的记忆效应及其对空泡形态的影响，详细计算了超空泡航行体各部分所受的流体动力，并提出了避免滑行力出现的预测和控制方法；张伟等^[5]进行了不同速度、形状弹体入水试验，得到了3种弹体的入水初期空泡形态与运动特性，建立了平头运动体入水的空泡形态和弹道预测模型；王柏秋等^[6]基于动网格技术对超空泡射弹的阻力系数开展了数值研究，结果表明同一模型长细比下，模型阻力系数随时间呈准线性变化；马庆鹏等^[7]开展了不同速度下的锥头柱体高速入水过程数值研究，发现入水速度越高，航行体头部压力峰值越大，相同时刻入水空泡最大直径也越大；何乾坤等^[8-14]对超空泡射弹的尾拍特性进行了大量数值研究.

目前国内外对并联入水的研究较少.卢佳兴等^[15]基于高速摄像技术开展了圆柱体并联入水空泡演化特性的实验研究，发现各深度空泡截面扩张和空泡中心偏移的峰值和持续时间均随弗劳德数增大而先增后减；宋武超等^[16]基于势流理论提出了回转体低速并联入水过程空泡形态发展的预测方法.综上所述，对入水问题的研究多针对单体入水，涉及并联入水的研究非常少，且未考虑空化现象对运动体并联入水流体和运动的影响.

本文基于对急动度的限制性假设，推导出描述圆柱体并联入水侧向位移和偏航角的预测公式，研究初始空化数和初始净距等初始条件对侧向和偏航运动特性的影响及其机理，提出避免回转体运动失稳的策略.

1 数值计算方法 1.1 控制方程及其求解

本文基于有限体积法对雷诺时均的N-S方程进行离散，使用VOF多相流模型描述各相分布. 混合介质的连续性方程为

$ \frac{\partial \rho_{m}}{\partial t}+\frac{\partial}{\partial x_{i}}\left(\rho_{m} u_{i}\right)=0, $

(1)

其中i=1, 2, 3，下同. 动量方程为

$ \begin{array}{c} \frac{\partial}{\partial t}\left(\rho_{m} u_{i}\right)+\frac{\partial}{\partial x_{i}}\left(\rho_{m} u_{i} u_{j}\right)= \\ -\frac{\partial p}{\partial x_{i}}+\frac{\partial}{\partial x_{i}}\left[\left(\mu_{m}+\mu_{t}\right)\left(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right)\right]. \end{array} $

(2)

式中：ρ_m=α_lρ_l+α_gρ_g+α_vρ_v，μ_m=α_lμ_l+α_gμ_g+α_vμ_v.α_l、α_g、α_v分别为水、空气和水蒸气的体积分数；ρ_l、ρ_g、ρ_v为三相的密度；μ_l、μ_g、μ_v为三相的动力黏度；μ_t=ρ_mC_μk²/ε为湍流黏性系数；C_μ为经验常数；k为湍动能；ε为湍动耗散率；u_i、u_j为速度分量；x_i、x_j为位移分量.C_μ由下式确定：

$ C_{\mu}=\frac{1}{A_{0}+\frac{A_{s} k U^{*}}{\varepsilon}}. $

(3)

式中: A₀=4.04；${A_s} = \sqrt 6 {\rm{cos}}\left[ {\frac{1}{3}{\rm{arccos}}\left( {\sqrt 6 \frac{{{S_{i, j}}{S_{j, k}}{S_{k, j}}}}{{{{\left( {{S_{i, j}}{S_{i, j}}} \right)}^{3/2}}}}} \right)} \right]$, 其中${S_{i, j}} = \frac{1}{2}\left( {\frac{{\partial {u_j}}}{{\partial {x_i}}} + \frac{{\partial {u_i}}}{{\partial {x_j}}}} \right);{U^*} = \sqrt {{S_{i, j}}{S_{i, j}}} $.

本文采用适用于大雷诺数流动的realizable k-ε模型^[17]，其湍动能和湍动耗散率的输运方程分别为:

$ \begin{array}{l} \frac{\partial(\rho k)}{\partial t}+\frac{\partial\left(\rho u_{j} k\right)}{\partial x_{j}}=\frac{\partial}{\partial x_{j}}\left[\left(\mu+\frac{\mu_{t}}{\sigma_{k}}\right) \frac{\partial k}{\partial x_{j}}\right]+ \\ G_{k}+G_{b}-\rho \varepsilon-Y_{M}+S_{k}, \end{array} $

(4)

$ \begin{array}{l} \frac{\partial(\rho \varepsilon)}{\partial t}+\frac{\partial\left(\rho \varepsilon u_{j}\right)}{\partial x_{j}}=\frac{\partial}{\partial x_{j}}\left[\left(\mu+\frac{\mu_{t}}{\sigma_{\varepsilon}}\right) \frac{\partial \varepsilon}{\partial x_{j}}\right]+ \\ \rho C_{1} S_{\varepsilon}-\rho C_{2} \frac{\varepsilon^{2}}{k+\sqrt{v \varepsilon}}+C_{1 \varepsilon} \frac{\varepsilon}{k} C_{3 \varepsilon} G_{b}+S_{\varepsilon} . \end{array} $

(5)

式中: G_k为时均速度梯度引起的湍动能；G_b为浮力引起的湍动能；Y_M为脉动扩张在湍流耗散率所占当量；${C_1} = \max \left[ {0.43, \frac{\eta }{{\eta + 5}}} \right]$；C₂=1.9；C_1ε、C_2ε、C_3ε为预测常数，一般取C_1ε=1.44, C_2ε=1.92；σ_k、σ_ε分别为k和ε的湍流普朗特常数，一般取σ_k=1, σ_ε=1.2；S_k、S_ε为源项.

本文采用Schnerr and Sauer空化模型^[18]描述空化现象.水蒸气相输运方程为

$ \begin{array}{l} \frac{\partial \alpha_{v}}{\partial t}+\frac{\partial}{\partial x_{i}}\left(\alpha_{v} u_{i}\right)=F_{\mathrm{vap}} \frac{2 \alpha_{\text {nuc }}\left(1-\alpha_{v}\right) \rho_{v}}{R_{B}}\cdot \\ \sqrt{\frac{2}{3} \frac{p_{v}-p}{\rho_{l}}}-F_{\text {cond }} \frac{3 \alpha_{v} \rho_{v}}{R_{B}} \sqrt{\frac{2}{3} \frac{p-p_{v}}{\rho_{l}}}. \end{array} $

(6)

式中: R_B=1×10^－6m为气核半径；α_nuc=5×10^－4为不可凝结气体体积分数；p为远场压力；p_v为饱和蒸气压；F_vap=50, F_cond=0.001.

1.2 方法有效性验证

文献[4]中给出了同实验结果符合良好的空泡形态模型，其表达式为

$ R(z)=\left[R_{0}^{2}+2 R_{0} \sqrt{\frac{C_{d}}{2 N}}\left(z-z_{0}\right)-\frac{\sigma}{2 N}\left(z-z_{0}\right)^{2}\right]^{1 / 2}. $

(7)

式中: R为空泡半径；R₀为头部半径；z为位移；z₀为初始位移；σ₀、σ分别为初始空化数和空化数；C_d= C₀(1+σ), 0.82≤C₀≤0.83；N为经验系数，取为2.

为验证计算方法的有效性，对半径5 mm的圆柱体以98.7 m/s入水过程进行数值计算，并和式(7)进行对比.由图 1中可知两者较吻合，说明本文数值计算方法是可信的.

2 侧向位移和偏航角预测公式

并联入水具有对称性，两回转体的侧向位移和偏航角大小相等，方向相反.本文约定质心相互远离的方向为侧向位移的正方向，头部远离而尾部靠近的旋转方向为偏航角的正方向，如图 2所示. 图 3给出了直径D=10 mm，长l=60 mm铝制圆柱体在初始空化数σ=0.02，初始净距d=0.6D条件下并联入水的侧向和偏航运动曲线. 由图 2可以看出，由于头部靠近内侧的压力较大，回转体受到负方向偏航力矩，导致负方向偏航运动；偏航运动使头部压力对回转体产生正向附加侧力，引起正方向侧向位移.

进一步考察侧向加速度和偏航角加速度可知，入水初期两者波动性均较强，随着时间增加逐渐趋于线性变化.因此若忽略波幅较小的波动，可近似认为t较小时侧向加速度和偏航角加速度曲线均有4个极值点，而t充分大时两者为线性变化，如图 4所示.因此由罗尔中值定理可知，侧向急动度和偏航角急动度在t较小时都至少有3个极值点，而当t充分大时两者都几乎为零.于是对侧向急动度和偏航角急动度作如下两个限制性假设：1)侧向急动度和偏航角急动度曲线均有3个极值点；2)侧向急动度和偏航角急动度在t充分大时均收敛于零.

据此设侧向急动度和偏航角急动度分别为:

$ j_{X}(t)=\sum\limits_{n=0}^{3} c_{X n} t^{n} f_{X n}(t), j_{\theta}(t)=\sum\limits_{n=0}^{3} c_{\theta n} t^{n} f_{\theta n}(t). $

(8)

式中：c_Xn、c_θn为待定参数；f_Xn(t)、f_θn(t)为待定的单调减函数，为使假设成立应满足:

$ t^{n} f_{X n}(\infty)=0, t^{n} f_{\theta n}(\infty)=0. $

(9)

对式(8)积分，并注意到积分第二中值定理可得：

$ a_{X}=\sum\limits_{n=0}^{3} c_{X n} t^{n} F_{X n}(t), a_{\theta}=\sum\limits_{n=0}^{3} c_{\theta n} t^{n} F_{\theta n}(t), $

(10)

其中：F_Xn(t)、F_θn(t)由下式确定：

$ \begin{array}{l} F_{X n}(t)=\int_{\xi_{X_{n}}(t)}^{t} f_{X n}(t) \mathrm{d} t, F_{\theta n}(t)=\int_{\xi_{\theta n}(t)}^{t} f_{\theta n}(t) \mathrm{d} t \\ \xi_{X n}(t), \xi_{\theta n}(t) \in[0, t]. \end{array} $

(11)

同样地，应用积分第二中值定理对式(10)连续积分两次得到侧向位移和偏航角分别为:

$ X=\sum\limits_{n=0}^{3} c_{X n} t^{n} \varLambda_{X_{n}}(t), \theta=\sum\limits_{n=0}^{3} c_{\theta n} t^{n} \varLambda_{\theta n}(t). $

(12)

令c_Xi=c_θi=0(i=0, 1, 2), 以及c_X=6c_X3, c_θ=6c_θ3；并取Λ_X3(t)=exp(－c_rXd$\sqrt t $/σD), Λ_θ3(t)= exp(－c_rθd t /σD)，得到侧向位移和偏航角预测公式:

$ \begin{array}{l} X=c_{X} t^{3} \exp \left(-c_{r X} d \sqrt{t} / \sigma D\right), \\ \theta=c_{\theta} t^{3} \exp \left(-c_{r \theta} d \sqrt{t} / \sigma D\right). \end{array} $

(13)

式中: c_X、c_θ分别为侧向和偏航角名义急动度；c_rX、c_rθ分别为侧向和偏航角修正因子.这些参数对于确定的并联入水过程是常数，而对于给定外形和材质的回转体则是初始空化数和初始净距的函数.

3 运动特性分析

取表 1中各初始条件的数值仿真数据对式(13)进行拟合，如图 5所示，拟合结果同式(13)吻合.

从图 5还可以看出，同一时刻侧向位移和偏航角随初始空化数减小而增加，随初始净距增加而先增后减.结合图 2对此解释如下：头部压力随初始空化数减小而增大，因此初始空化数越小，回转体所受偏航力矩和附加侧力越大，从而侧向位移和偏航角越大；速度驻点位置随初始净距变小而越靠近内侧，但初始净距较小时，速度驻点位置几乎不再向内侧偏移，故初始净距较大时，偏航力矩随初始净距减小而增大，使偏航角随初始净距减小而增大，进而使附加侧力随初始净距减小而增大，而初始净距较小时，内侧靠近尾部的低压区对回转体有正方向力矩，使偏航力矩和附加侧力随初始净距减小而减小，从而使侧向和偏航运动受到限制.

为定量考察初始空化数和初始净距对侧向和偏航运动的影响机理，图 6给出了名义急动度c_X、c_θ和修正因子c_rX、c_rθ随初始条件的变化.可以看出，c_X和c_θ随着σ减小而增大，而随d/D增大而先增后减；c_rX和c_rθ随d/D减小而增大.这进一步说明初始空化数越小，头部压力对侧向和偏航运动的促进作用越强，而初始净距越小，内侧低压对侧向和偏航运动的限制作用越强.

4 避免失稳的策略

侧向和偏航运动使回转体逐渐偏离平衡点，可能使回转体碰撞或过度远离，导致运动失稳.因此获得碰撞和过度远离的条件很有必要.

定义运动状态因数β为净距变化量与初始净距之比的相反数，即

$ \beta=-\frac{2 X+D(1-\sec \theta)+(l+D \tan \theta) \sin \theta}{d}, $

(14)

回转体发生碰撞时，头部内侧与对称线接触，如图 7所示.此时满足如下几何关系:

$ \begin{array}{l} X+\frac{1}{2} d+\frac{1}{2} D-\frac{1}{2} D \sec (-\theta)= \\ \left(\frac{1}{2} l-\frac{1}{2} D \tan (-\theta)\right) \sin (-\theta). \end{array} $

(15)

亦即

$ \beta=-\frac{2 X+D(1-\sec \theta)+(l+D \tan \theta) \sin \theta}{d}=1, $

(16)

上式即是碰撞状态的判据.本文约定|β|≤1，当且仅当β=1时回转体碰撞，当且仅当β=-1时为过度远离状态.式(14)可改写为

$ \beta=-\frac{2 X / D+(1-\sec \theta)+(0.25 \lambda+\tan \theta) \sin \theta}{d / D}, $

(17)

其中，λ为长细比，对于本文的回转体λ=24.

将式(13)代入式(17)即得β表达式.该式表明，运动状态不仅和初始条件(初始空化数和初始净距)有关，还同回转体的外形(长细比)相关.

图 8比较了不同初始条件对β的影响.可以看出，对同一初始净距，存在使回转体最接近初始状态的初始空化数；对同一初始空化数，存在使回转体净距最接近初始状态的初始净距.而为避免运动失稳，应当力求|β|尽量小.因此对于已确定形状和材料的回转体，应控制入水初速度和初始净距在合理的范围内.对于本文的圆柱体，若初始净距为0.4D，则入水初速度应控制在79.64 m/s左右(对应于图 8中σ=0.031, d/D=0.4)；若入水初速度为98.7 m/s，则初始净距应控制在0.4D~0.6D之间(对应于图 8中σ=0.020, d/D=0.4及σ=0.020, d/D=0.6).

5 结论

1) 预测公式与数值计算结果符合良好，具有一定的可信度.

2) 入水后同一时刻，运动体的侧向位移和偏航角均随初始空化数减小而增大，随初始净距减小而先增后减.

3) 初始空化数较低时头部压力对侧向和偏航运动具有促进作用，且初始空化数越小，侧向名义急动度和偏航角名义急动度越大，表明侧向和偏航运动所受促进作用越强；初始净距较小时，内侧压力对侧向和偏航运动具有抑制作用，且初始净距越小，侧向修正因子和偏航角修正因子越大，表明侧向和偏航运动所受抑制作用越强

4) 初始条件和回转体外形共同决定了并联入水过程中回转体的运动状态.对于本文运动体，当初始净距为0.4D，则入水初速度应控制在79.64 m/s左右；若入水初速度为98.7 m/s，则初始净距应控制在0.4D~0.6D之间.