权晓波， 尤天庆， 张晨星， 王凡瑜， 孔德才

(1.中国运载火箭技术研究院， 北京 100076； 2.北京宇航系统工程研究所， 北京 100076)

0 引言

航行体水下垂直弹射出筒后，发射筒内燃气会附着在航行体尾部形成尾空泡[1]。在水下垂直向上运动过程中，尾空泡呈现显著的瞬态振荡演化特性，对周围流动产生明显影响，是航行体水下垂直发射过程的主要多相流动现象[2-3]。

尾空泡发展演化会对水下发射流体动力产生重要影响，对其研究目前主要在压力振荡、演化形态以及气体与水(简称气水)掺混流态3方面进行。在压力振荡研究方面，鉴于尾空泡对航行体轴向弹道的严重影响，研究者基于气泡动力方程对尾空泡压力振荡进行模拟并耦合航行体运动求解，获得与试验吻合较好的压力振荡规律[4-5]。在空泡演化形态方面，研究者基于弹射试验获得的尾空泡演化物理特征认识，基于空泡独立膨胀原理或势流理论，建立尾空泡形态计算模型，获得了出筒速度等因素对空泡演化的影响规律[6-8]。此外基于Navier-Stokes方程的多相流模拟技术发展，使尾空泡形态演化模拟考虑的影响因素更加丰富[2,9]。在尾空泡气水多相流态方面，研究多数聚焦于多相流态对水下喷气射流影响，以及多相流的膨胀断裂等复杂流动现象[10-11]。

水下发射水介质不可压缩流动特性决定了尾空泡演化对流体动力的影响并不局限于航行体尾部。除上述尾空泡3方面影响外，航行体上游物面附近流动也应考虑尾空泡扰动。针对此问题研究，势流算法可抓住不可压流动主要特征，在假设无黏、无旋条件下，以较小计算资源消耗进行影响因素研究。同时文献[12-15]表明，势流算法在进行空泡流计算时，体现了较好的精度。

针对水下发射尾空泡演化对上游流动的影响问题，本文主要基于二维轴对称势流边界元方法，开展了针对水下垂直发射航行体尾部附着空泡的数值计算研究。重点关注航行体物面压力计算，通过与水洞试验数据对比，验证了数值计算方法的有效性。在此基础上，针对尾空泡瞬态振荡影响，分别从上游物面压力、沾湿区阻力以及空泡发展3方面进行了数值计算研究。

1 数值计算方法

1.1 控制方程

假设空泡附着在物面扰流低压区，在无限大流体域内振荡。并假设流体无黏、无旋且不可压，则Navier-Stokes方程可简化为拉普拉斯方程。

(1)

(2)

式中:Φ为总速度势，Φ=ux+φ，u为航行体运动速度(沿x轴负向)，x为位置坐标，φ为扰动速度势；v为速度向量。

在物面，流体需满足不可穿透条件。即壁面法向速度为0 m/s，

(3)

式中：n为物面法向方向向量。

在空泡壁面，速度势Φ需满足伯努利方程：

(4)

式中：p0为水面压力；ρ为水的密度；g为重力加速度；pc为空泡内压力。空泡壁面坐标x需满足：

(5)

航行体表面压力变化是本文研究的重点关注参数，流场中任意一点压力p，可由伯努利方程计算得到：

(6)

式中：∂Φ/∂t直接采用速度势时间差分进行求解是比较困难的。本文求解中利用∂Φ/∂t也满足拉普拉斯方程的特点，采用与φ求解相同的数值方法。

1.2 数值算法

本文在数值方法上采用间接边界元法求解，以求解扰动速度势φ为例，介绍求解方法。设空间场点为p，物面和空泡壁面上分布的点源位置为q，其强度为λ(q)，则其在p点产生的速度势和速度分别为

(7)

(8)

式中：G(p,q)为格林函数；S为面积；np为p点速度方向矢量；dSq为q点附近微元面积。对于拉普拉斯方程，三维情况下格林函数G(p,q)=-1/4πr，其中r=|p-q|.

数值计算中，航行体与空泡表面被划分为N个单元，对航行体表面的第i个单元，考虑在边界j(j=1，2，…，N)上源汇自身对速度的影响，扰动势所满足的边界元积分控制方程如(9)式所示：

(9)

式中：Sj为第j个单元面积；u为航行体运动速度向量；nn为单元法向向量；λi为第i个单元的强度；λj为第j个单元的强度；dSqj为qj点附近微元面积。

空泡分离点为位于第d个单元，其扰动势所满足的边界元积分控制方程如(10)式所示：

(10)

(9)式和(10)式即为边界元计算的基本方程，通过此方程即可求得物面及空泡壁面的λ分布，进而确定流场中速度势、速度及压力。

为节省计算资源提高计算效率，本文采用二维轴对称计算模型，在航行体及空泡表面的边界单元为圆环状。与速度势及速度相关的源汇均匀沿圆环分布。对于子午线长度为l圆环边界单元，其对空间点p(xp,yp,0)的影响速度势和速度由圆环上源汇q(xq，yqcosθ，yqsinθ)积分获得，θ为圆柱坐标系下的方位角，影响速度势和速度表达式分别为

(11)

(12)

(13)

由于第一类完整椭圆积分函数的奇异性，当R趋近于1时，(11)式～(13)式所示的边界单元速度势和速度影响函数值趋近于无穷大，但其积分是有限的。本文利用偶数积分点的Gauss积分进行边界单元积分计算。

为验证数值计算代码的正确性，在无黏、无旋理想流体条件下，计算了A和B两种轴对称旋成体零攻角物面压力分布，并与水洞测压试验数据对比，如图1所示。图1中：Cp为压力系数，Cp=(p-p0)/0.5ρu2；X为航行体压力监测点的位置坐标，D为航行体直径；A为流线性外形，旋成体表面为全沾湿流动，对比表明数值计算与试验数据吻合较好；B为锥柱组合外形，旋成体锥肩部附着有低压空泡，数值计算结果较准确地描述了空泡低压区范围以及空泡末端回射压力变化趋势。

图1 数值计算与试验数据对比Fig.1 Comparison between calculated and experimental results

需要说明的是，空泡末端实际流动中存在较强的湍流现象，本文计算对此进行了简化处理，导致数值计算的回射压力量值与试验结果存在一定差异，但这并不妨碍后续肩空泡回射压力变化趋势分析。

2 尾空泡振荡影响分析

水下垂直发射航行体向水面运动过程中环境压力逐渐减小，尾空泡发展演化以及航行体物面流动参数不断变化。此过程流场的计算即涉及航行体相对固定的物面流场计算，又涉及空泡动态发展过程的计算。基于上述混合边界条件的边界元计算方法，可求得航行体物面和空泡壁面的流场速度、加速度以及压力等参数。在一个时间步内求得各物理量后，需通过边界条件和空泡形态，对求解变量进行更新，本文采用4阶龙格- 库塔法进行数值求解。在时间步推进过程中，由于求解误差积累，空泡壁面会产生锯齿状形态最终导致计算发散并终止，因此本文在空泡壁面动网格计算中进行了数值光顺处理。

本文将从尾空泡振荡过程模拟、对上游物面压力影响、对沾湿区阻力影响以及对肩空泡影响这四方面进行论述，参见图2.图2中：坐标X/D≤0部分为尾空泡上游航行体物面沾湿区；坐标X/D>0部分为尾空泡轮廓，其内部压力为pct；当航行体肩部物面压力较低时，会附着有肩空泡，其内部压力为pcn.

任何脱离了所在环境文化的能力都是不成熟的，新生代企业家的优秀的胜任力，一定是要能在现有的中国情境下，仍然能够出色发挥，而不是被动地等待环境来适应新生代企业家管理理念的胜任力。要融入中国情境，就必须先从多个角度来了解中国文化，符合大众的价值取向。同时，笔者也试着从中西方科学文化的共同、共通之处为起点，对相异之处进行对比论证，为存在的问题做出合理的解释。

图2 尾空泡振荡影响示意图Fig.2 Sketch map of oscillation influence of tail cavity

2.1 尾空泡振荡的模拟

为模拟航行体由一定深度向水面运动过程中尾空泡振荡过程，首先进行一定深度条件下空泡流场稳态计算，获得σ0=(p0+ρgh-pct)/0.5ρu2条件下空泡流场，其中σ0为此稳态计算获得流场的初始空泡数，h为深度。以此稳态计算结果作为后续瞬态计算初始条件，改变空泡内初始压力，使空泡数为σ′0，引入尾空泡振荡扰动，尾空泡压力振荡如图3所示，其中t1、t2、t3、t4、t5为空泡收缩过程中等时间间隔点。由于在计算过程中航行体垂直向上运动，周围环境压力逐渐降低，因此尾空泡压力在振荡中逐渐降低。此外，本文数值算法对空泡壁面运动速度和空泡形态进行了平滑处理，在计算过程中引入了数值黏性，因此压力振荡体现出了衰减特征。

图3 尾空泡内压力振荡数值计算结果Fig.3 Numerically calculated results of pressure oscillation of tail cavity

尾空泡振荡过程典型时刻空泡形态如图4所示。对应空泡瞬态收缩泡内压力升高的过程，t1～t5等时间间隔，具体时刻如图3中标注。由空泡振荡形态变化可见，空泡壁面收缩或膨胀为变速过程。例如针对收缩过程，t1～t2时间段空泡壁面具有较大加速度，但收缩速度较小，整体表现t1和t2时刻空泡形态变化不大。随着空泡壁面加速度逐渐累加，空泡壁面加速度收缩，空泡形态差异逐渐明显。

图4 空泡收缩过程形态变化Fig.4 Profile change of cavity during shrinking

2.2 尾空泡振荡对上游物面压力影响

尾空泡瞬态振荡使周围水介质运动加速度和速度产生显著振荡，进而对空泡上游航行体物面压力产生明显影响，如图5所示。图5中：T为压力波动周期，ΔCp，max为各位置压力系数变化量最大值。随着航行体向水面运动，环境压力逐渐降低，物面压力逐渐下降，尾空泡振荡使上游航行体物面压力下降过程叠加了一定幅值的压力波动。尾空泡起始位置为X/D=0，相应位置压力系数变化曲线表明，测点位置越接近尾空泡压力波动幅度越大。

图5 尾空泡振荡引起的沾湿区物面压力波动Fig.5 Pressure oscillation on wet surface caused by tail cavity

为分析物面压力波动规律，滤去压力变化时程曲线中上升运动引起的静压变化，仅关注流体运动速度和加速度变化所引起的压力系数变化量ΔCp=Cp-ρgh/0.5ρu2，采用ΔCp，max和T对压力系数曲线进行无量纲化。取水下垂直发射弹射试验航行体物面典型测点压力，与数值计算结果对比，如图6所示。二者在压力变化规律上吻合较好，一方面体现在压力随时时间变化趋势上，另一方面体现在压力波动在空间上的衰减趋势上。

图6 尾空泡压力振荡计算与试验数据对比Fig.6 Comparison between calculated and experimental results

为进一步分析尾空泡振荡对航行体上游物面压力影响在空间上的分布规律，提取了典型时刻物面压力波动量值沿空间分布情况，如图7所示。数值计算结果表明，尾空泡振荡引起的物面压力扰动以近似指数函数在空间分布，因此基于上述分析认为，当航行体长径比较小时，尾空泡振荡对航行体整体受力及上游流动会有较为明显的影响。

图7 不同初始压力尾空泡收缩引起的物面压力变化Fig.7 Pressure change of tail cavity during shrinking

2.3 尾空泡振荡对航行体沾湿区阻力影响

计算长径比L/D为6情况下，空泡振荡条件下航行体垂直向水面运动过程受力情况，如图8所示。其中，航行体阻力系数Cd由物面压力积分求得

图8 尾空泡振荡对沾湿区阻力影响Fig.8 Effect of tail cavity oscillation on resistance in wet area

(14)

式中：Sw为航行体沾湿区域；Sr为参考面积；nx为微元法向向量x轴分量。

尾空泡振荡引起的航行体上游物面压力波动明显地表现在阻力系数上，如图8所示。其中无尾空泡振荡影响情况，即在向水面运动过程中空泡形态保持不变的计算结果。由于向水面运动过程中，随着环境压力逐渐降低，尾空泡振荡过程中形态逐渐增长，导致周围流场加速度和速度变化区域增加。因此，相比尾空泡无变化工况，尾空泡振荡对阻力系数影响，体现为振荡中的缓慢变化。

2.4 尾空泡振荡对肩空泡流场影响

尾空泡振荡引起的上游物面压力波动，不仅对航行体受力产生明显影响，在上游存在空泡情况下，也会对上游空泡流场演化产生影响。针对此问题，本文数值计算获得了一定泡内压力条件下肩尾空泡稳态流场，并以此稳态计算结果为初始条件，通过降低尾空泡泡内压力，引入流场扰动开始瞬态计算。计算一定时间段后，获得肩尾空泡流场，如图9所示。为方便对比分析，在计算过程中忽略尾空泡壁面速度，使尾空泡形态保持不变，进而获得忽略尾空泡振荡影响的计算结果。

由于在瞬态计算初始条件给定中降低了尾空泡内压力，进而导致航行体周围压力显著下降，肩空泡逐渐开始膨胀。由于尾空泡收缩使流场上游压力提高，因此肩空泡发展相对缓慢。这种影响体现在受尾空泡影响下肩空泡尺寸相对较小，如图9所示，同时也体现在肩空泡泡内压力变化上，如图10所示。在肩空泡内同样气体含量条件下，尾空泡收缩过程明显使肩空泡泡内压力提升。

图9 尾空泡收缩对肩空泡发展的影响Fig.9 Effect of tail cavity shrinking on evolution of shoulder cavity

图10 尾空泡收缩对肩空泡内部压力影响Fig.10 Effect of tail cavity shrinking on pressure of shoulder cavity

在肩空泡末端，空泡壁面水流凹向壁面，使其损失径向速度，形成回射压力峰。尾空泡振荡在改变肩空泡发展速度情况下，使航行体沾湿区压力显著提升，进而使肩空泡末端回射压力峰值增大，如图11所示。

图11 尾空泡收缩对肩空泡流场压力系数影响Fig.11 Effect of tail cavity shrinking on pressure distribution of shoulder cavity

3 结论

本文通过求解二维轴对称航行体水下垂直运动过程空泡流场瞬态演化，分析了尾空泡振荡对航行体上游物面压力、受力以及空泡演化的影响。得到主要结论如下：

1) 通过与水洞测压数据比对，验证了本文数值方法对航行体上游物面压力计算的准确性，同时航行体附着尾空泡的混合边界元计算可模拟尾空泡泡内压力以及空泡形态的瞬态振荡过程。

2) 在尾空泡振荡对上游物面压力扰动计算方面，计算结果取得与弹射试验数据相一致的变化规律；同时计算分析表明，随着关注位置远离尾空泡，压力扰动呈指数函数衰减。

3) 针对航行体上游物面整体受力，随着航行体垂直向水面运动环境压力逐渐降低，尾空泡振荡过程中长度逐渐增加，其振荡对阻力系数影响，体现为振荡中的缓慢变化。

4) 针对航行体上游肩空泡流场演化，肩尾空泡耦合计算表明，尾空泡收缩使流场上游肩部空泡发展速度减缓，同时使上游肩空泡泡内压力和末端回射压力峰值偏大。