0 引　言

水下航行体通过高压工质气体快速运动，实现水下垂直弹射^[1]。当航行体离开运载器筒口后，高温高压气体涌出筒口形成气团，并在航行体后续运动过程中与水环境相互作用发生膨胀、收缩、缩聚和断裂等现象^{[2 − 3]}，形成筒口压力场，对邻近设备造成显著的冲击载荷^[4]。筒盖作为直接承受筒口压力场的机构，其结构强度和布置形式直接影响了运载器的安全性、稳定性和隐蔽性^[5]，同时也通过干扰压力场在一定程度上影响了航行体的运动姿态。为降低上述风险，优化筒盖布置形式成为一种迫切需求。

为精确预测筒口气团作用下的筒盖受载，方国强等^[6]基于试验测试数据采用最小二乘法获得了计算筒盖最大载荷的经验公式，计算结果与试验结果对比平均误差在13%以内，取得了较好的一致性。王汉平等^[7]通过神经网络和遗传算法完成了复杂数据的处理和模型优化，提出了一种基于遗传算法进化神经网络的筒盖最大压力预测方法，预测结果与测试结果最大相对误差为9.54%，具备较强的实用价值。数值方法上，程载斌等^[8]采用Ls-dyna中的ALE方法建立了水、空气、航行体和筒盖等结构的耦合模型，真实反映了航行体水下垂直发射过程。祁晓斌等^[9]采用计算流体力学软件Fluent通过重叠网格中的局部重构和弹性光顺滑构建了具有柔性支撑的筒盖模型，研究了筒口气团膨胀-收缩周期性过程，获得了筒盖的运动规律和受载特性。在此基础上，姚熊亮等^[10]研究了不同运载体横向速度、筒口气团初体积、筒盖弯曲系数和筒盖开盖角度对筒盖承载的影响程度，揭示了筒盖载荷特性。当前，由于水下运载器对性能、尺度等因素有较高的要求，筒盖的许多设计参数受到限制，给相关优化工作带来了很大挑战。

本文针对筒盖设计面临的限制，提出了一种滑动筒盖的布置形式，采用CFD方法，建立了水下航行体离开筒口后筒内高温高压气体形成筒口气团作用于筒盖的计算模型，通过分析筒口气团演变过程、筒盖承载及其测点压力，对筒口气团影响下不同滑动高度的筒盖载荷特性进行了研究，为验证滑动筒盖优化设计提供参考。

1 数值方法 1.1 控制方程

水下航行体离开筒口后，高温高压气体与水环境两相流体间相互作用并发生对流传热，研究该问题需遵循3个基本控制方程。

质量守恒方程：

$ \frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \mathit{v}\right)=0 。$

(1)

式中：$ \rho $为流体密度；$ \mathit{v} $为流体速度；$ t $为时间。

动量守恒方程：

$ \frac{\partial \left(\rho \mathit{v}\right)}{\partial t}+\nabla \cdot \left(\rho \mathit{v}\otimes \mathit{v}\right)=\nabla \cdot \mathit{\sigma }+{\mathit{f}}_{{b}}。$

(2)

式中：$ \mathit{\sigma } $为应力张量，对于流体而言，通常写为法向应力和剪切应力的总和；$ {\mathit{f}}_{{b}} $为作用于流体单位体积的体积力合力。

能量守恒方程：

$ \frac{\partial \left(\rho E\right)}{\partial t}+\nabla \cdot \left(\rho E\mathit{v}\right)= {\mathit{f}}_{{b}}\cdot \mathit{v}+\nabla \cdot \left(\mathit{v}\cdot \mathit{\sigma }\right)-\nabla \cdot \mathit{q}+{S}_{E} 。$

(3)

式中：$ E $为流体单位质量的总能量；$ \mathit{q} $为热通量；$ {S}_{E} $为单位体积的能量源。

1.2 多相流传热模型

为较好地捕捉筒口气团在水中的膨胀、收缩、溃灭等现象，采用VOF多相流模型描述两相流体交界面处相i的分布和位置，其相体积分数$ {\alpha }_{i} $为：

$ {\alpha }_{i}=\frac{{V}_{i}}{V}。$

(4)

式中：$ {V}_{i} $为相i在网格单元中的体积；$ V $为网格单元的总体积。则空气相$ {\alpha }_{{g}} $和水相$ {\alpha }_{{w}} $在网格单元中的体积分数之和遵循：

$ {\alpha }_{{g}}+{\alpha }_{{w}}=1 。$

(5)

同一网格单元中流体的密度$ \rho $和动力粘度$ \mu $表示为：

$ \rho ={\rho }_{{g}}{\alpha }_{{g}}+{\rho }_{{w}}{\alpha }_{{w}}，$

(6)

$ \mu ={\mu }_{{g}}{\alpha }_{{g}}+{\mu }_{{w}}{\alpha }_{{w}} 。$

(7)

式中：$ {\rho }_{{g}} $和$ {\mu }_{{g}} $分别为空气相的密度和动力粘度；$ {\rho }_{{w}} $和$ {\mu }_{{w}} $分别为水相的密度和动力粘度。

高温气体与水环境间的自然对流以及航行体运动影响下的强制对流影响着筒口气团的演化，对流热传递遵循符合牛顿冷却定律：

$ {\dot{\mathit{q}}}^{{'}{'}}=h({T}_{{g}}-{T}_{{w}}) 。$

(8)

式中：$ {\dot{\mathit{q}}}^{{'}{'}} $为局部热通量；$ h $为局部对流热传递系数；$ ({T}_{{g}}-{T}_{{w}}) $为两相间的温度差。

1.3 湍流模型

湍流模型的应用既符合航行体出筒的流场状态，同时会影响两相流体间的热传递系数。筒盖载荷主要来源于筒口整体流场演化形成的筒口压力场，近壁面处的流场影响不大。因此，采用适用性较好的Realizable k−ε两层湍流模型开展数值计算。该模型关于湍动能$ {k}_{t} $和湍动耗散率$ {\varepsilon }_{t} $的传输方程如下：

$ \frac{\partial }{\partial t}\left(\rho {k}_{t}\right) + \nabla \cdot \left(\rho {k}_{t}\stackrel-{\mathit{v}}\right) = \nabla \cdot \left[\left(\mu + \frac{{\mu }_{t}}{{\sigma }_{k}}\right)\nabla {k}_{t}\right] + {P}_{k} - \rho \left({\varepsilon }_{t} - {\varepsilon }_{0}\right) + {S}_{k} ，$

(9)

$ \begin{split} &\frac{\partial }{\partial t}\left(\rho {\varepsilon }_{t}\right)+\nabla \cdot \left(\rho {\varepsilon }_{t}\bar{v}\right)=\nabla \cdot \left[\left(\mu +\frac{{\mu }_{t}}{{\sigma }_{\varepsilon }}\right)\nabla {\varepsilon }_{t}\right]+\\ &\quad\frac{1}{{T}_{e}}{{C}_{\varepsilon 1}P}_{\varepsilon }-{C}_{\varepsilon 2}{f}_{2}\rho \left(\frac{{\varepsilon }_{t}}{{T}_{e}}-\frac{{\varepsilon }_{0}}{{T}_{0}}\right)+{S}_{\varepsilon }。\end{split}$

(10)

式中：$ {\mu }_{t} $为粘涡系数，是方程的求解项；$ \bar{v} $为平均速度；$ {\sigma }_{k} $和$ {\sigma }_{\varepsilon } $分别为与湍动能$ {k}_{t} $和湍动耗散率$ {\varepsilon }_{t} $对应的普朗数（Pr）；$ {T}_{e}={k}_{t}/{\varepsilon }_{t} $为大涡时间尺度；$ {P}_{k} $和$ {P}_{\varepsilon } $为湍动能产生项；$ {C}_{\varepsilon 1} $和$ {C}_{\varepsilon 2} $均为经验参数；$ {f}_{2} $为阻尼函数；$ {S}_{k} $和$ {S}_{\varepsilon } $均为源项；$ {\varepsilon }_{0} $为源项中抵消湍流衰减的环境湍流值；$ {T}_{0} $为施加源项的同时产生的单位时间尺度。

2 计算模型及验证 2.1 模型假定

筒盖载荷主要受航行体离开筒口后形成的筒口气团的影响，因此，本文关注水下航行体离开发射筒口，高温高压燃气与水环境开始接触后的过程，并针对筒口气团作用下不同滑动距离的筒盖响应情况进行数值建模并进行一定简化：

1） 简化了航行体尾部、筒口等复杂结构，优化了整体型线，以利于数值计算收敛和提高计算效率；

2） 航行体出筒后横向运动对筒口压力场影响轻微，简化为方向向上的单自由度运动，每个工况运动速度一致；

3） 以航行体出筒瞬间状态为基准，筒内燃气压力和温度视为均匀分布，燃气为理想气体；

4） 空化汽化模型对筒口压力场影响较小^[11]，本文不考虑空化汽化现象；

5） 为研究不同滑动距离的筒盖布置形式在筒口气团作用下的响应情况，忽略实际设计中筒盖因不同滑动距离而在角度方向等方面的变化。

2.2 计算模型

图1为初始状态水下航行体、筒体及筒盖的布置情况和计算流体域的边界条件。航行体直径为D，长度为L，与筒体及筒盖均为壁面。计算水域为正方体，边长为2L，空气域分布在筒内。计算域底面为壁面边界条件，以模拟运载体壁面，与筒口存在一定距离；其余面为压力出口边界条件；航行体在流体域中的运动通过重叠网格实现。

图2为网格划分情况，分别对航行体运动路径和筒内空间进行加密，并针对筒口气团演化区域及筒盖部分进行局部加密，采用自适应时间步提高计算的收敛性和效率。

对于计算模型的初始状态，航行体出筒速度为V，筒内高温高压气体压力为P₀，温度为T₀。为研究筒口气团对不同滑动距离筒盖的影响规律，在筒盖内表面设置了7个压力测点，如图3所示。

2.3 模型验证

为衡量不同产品的载荷特性，本文对数据进行了无量纲化处理：

$ m=\frac{{m}_{t}}{{m}_{\mathrm{c}\mathrm{o}\mathrm{v}\mathrm{e}\mathrm{r}}} ，$

(11)

$ t=\frac{{t}_{0}}{V/g}，$

(12)

$ h=\frac{{h}_{\mathrm{c}\mathrm{o}\mathrm{v}\mathrm{e}\mathrm{r}}}{D}，$

(13)

$ P=\frac{{P}_{t}}{{m}_{\mathrm{c}\mathrm{o}\mathrm{v}\mathrm{e}\mathrm{r}}g/{D}^{2}}，$

(14)

$ T=\frac{{T}_{t}}{{T}_{\mathrm{R}\mathrm{T}}} 。$

(15)

式中：$ m $、$ t $、$ h $、$ P $、$ T $分别为无量纲化的筒盖载荷、时间、筒盖高度、压力和温度；$ {m}_{t} $、$ {P}_{t} $、$ {T}_{t} $分别为航行体出筒过程中筒盖的实际载荷、环境压力和温度；$ {t}_{0} $为实际时间历程；$ g $为重量加速度；$ {m}_{\mathrm{c}\mathrm{o}\mathrm{v}\mathrm{e}\mathrm{r}} $和$ {h}_{\mathrm{c}\mathrm{o}\mathrm{v}\mathrm{e}\mathrm{r}} $分别为筒盖的质量和距筒口的实际高度；$ {T}_{\mathrm{R}\mathrm{T}} $为室温。

为验证计算模型的准确性，参考了单次试验测试结果，并将试验参数代入到式(16)中获得筒盖最大载荷的经验值进行对比验证^[6]。

$ F=74.82+0.069{V}^{{2}}。$

(16)

式中：$ F $为筒盖承受的最大载荷。当前算例筒盖最大载荷的对比结果见表1，可知，数值结果与试验结果和经验公式有较好一致性，当前计算模型能够较好地模拟筒口气团作用下的筒盖受载情况。

3 结果及分析 3.1 网格无关性分析

针对本文滑动筒盖优化设计研究模型，首先考虑了3种不同的网格尺度来验证模型的网格无关性，网格总数分别为185万、325万和600万，如图4所示。可知，在载荷波峰区间内，大网格尺度、低网格数时的载荷曲线峰值明显偏大，整体曲线有所偏移。为保证计算精度和效率，选取325万网格数下的网格尺度开展后续计算和分析。

3.2 流场演化分析

图5为h = 0工况下，航行体离开筒口过程中筒口附近的水体积分数、压力场和温度场演化云图，可以清楚地观察到高温高压气体与水环境作用过程。结合图4，t=0时，航行体出筒，筒内气体迅速涌出筒口，高温高压气团快速膨胀，并作用于筒盖下端，施加正向载荷。期间，筒口气团温度不断降低，压力开始快速降低且低于水环境，逐渐产生相对于水环境的负压力场，筒盖开始承受反向载荷，如t = 0.085时所示。随后，在t = 0.153～0.170过程中，航行体尾部气团与筒口气团分离，重力场作用逐渐显著，筒口气团开始发生收缩，根据理想气体状态方程，气团压力和温度开始回升，与此同时，筒盖外侧压力也快速恢复，使筒盖承受的反向载荷继续增大，气团开始出现断裂现象。当t = 0.183状态时，气团体积继续收缩，温度和压力也不断升高，并在断裂处发生小范围溃灭，反向载荷达到峰值开始减小，气团发生完全断裂。在t = 0.194时，筒口气团温度和压力达到阈值，发生大范围溃灭，筒口上方压力骤升，筒盖正向载荷快速升高并达到峰值。能量释放完毕后，多余的气团在重力场作用下上浮，筒口压力和温度恢复至环境值，如t = 0.256所示。

3.3 参数化分析

为研究不同滑动距离的筒盖布置形式承受筒口气团作用的载荷特性，以筒盖原始高度作为基准工况h = 0，同时考虑了筒盖高度等差下降后的5个工况，即h = −0.215、−0.430、−0.645、−0.860、−1.075，开展后续计算和分析。

图6为筒盖在不同高度下的承载曲线。针对t=0～0.05的过程，筒盖受到筒口气团膨胀造成的正向压力随筒盖高度h的递减先增大后减小，在h = −0.430和h = −0.645状态时达到最大。

随后，筒盖开始承受负压力场后，筒盖高度h越小，反向载荷增长越快，当h≤−0.645后，增长速度趋于一致，但随着h的减小，气团开始溃灭的时间提前，峰值逐步减小，如图7所示。

筒口气团收缩溃灭阶段，筒盖承受正向载荷，载荷峰值随h的递减先增大后逐步减小，当h = −1.075时，又突然增大，如图8所示。相对于h = 0的工况，h = −0.860时的筒盖位置降载效果最好，最大载荷降低了28.0%，最大反向载荷降低了29.0%。因此，可以认为当前范围内筒盖与筒口的相对位置对筒盖最大载荷有较大影响。

3.4 机理分析

对2个关键时刻不同筒盖高度下筒口附近的水体积分数进行分析（见图8）。h = −0.430和h = −0.645时膨胀后的气团几乎完全正面作用于筒盖，导致这2个位置的筒盖在气团膨胀阶段承受了较大的正向载荷。

下一阶段，随着筒盖位置的递减，为筒口水体让出了更多空间。气团结束膨胀后在重力场作用下加快了对筒口气团的压缩，进而使气团溃灭时间提前，降低了反向载荷峰值，这在图8中可以明显捕捉到该现象。

在气团溃灭阶段，当h = 0和h = −0.215时，筒盖基本位于筒口以上，从图8和图9中对比2个工况可以看到，尽管h = −0.215时筒盖下端略低于筒口减弱了承载，但其余部分与筒口更为接近，使影响更加明显，造成了更大的正向载荷。如图9所示，随着筒盖高度的进一步递减，筒口气团溃灭形成压力场的影响逐渐减弱，筒盖峰值测点逐步上移，整体载荷降低。

在t = 0.194时不同筒盖高度下筒口附近的水体积分数云图中可以注意到，部分工况在筒口气团收缩溃灭阶段观察到了局部溃灭现象。当h ≥ −0.215时，筒口气团只发生了单次整体溃灭；而当h为−0.430、−0.645、−0.860时，由于筒盖的存在延缓了筒盖侧气团的收缩，同时，没有筒盖阻碍的水体在重力场作用下对该侧气团进行了分割，导致气团整体溃灭的同时发生局部溃灭。

对于h = −1.075时的工况，表现出了值得关注的结果。局部溃灭先于整体溃灭单独发生，能量未经耗散且靠近筒盖，造成了较大的筒盖承载，而后的整体溃灭由于部分能量释放和筒盖位置较低，未对筒盖承载造成明显影响。

针对部分工况气团收缩溃灭阶段出现第二波峰的现象，在图9中可以观察到h = −1.075时筒盖上端测点压力的急剧上升，本文认为可能是气团溃灭后的残余气团在筒盖上端部附近与水环境、筒盖三者相互作用导致的。

4 结　语

本文为验证一种滑动筒盖的优化设计，基于CFD方法建立并验证了筒口气团作用下的筒盖承载计算模型，对不同筒盖高度的承载特性进行了研究，得出以下结论：

1） 向下滑动的筒盖布置形式很大程度上避免了筒口气团的正面作用，能有效减小筒盖载荷，在现有设计限制下，当筒盖高度降低至h = −0.860时效果最好，最大载荷降低了28.0%，最大反向载荷降低了29.0%；

2） 筒盖高度的降低会加快筒口气团的溃灭，在筒盖尚未完全远离筒口的条件下，例如h = −1.075时，筒口气团甚至提前发生靠近筒盖的局部溃灭，增加筒盖承载。

滑动筒盖的布置形式不仅可以降低筒口气团作用下的筒盖承载，并且与相比原先的位置相比，包络空间更小，重心更低，有利于运载体噪声、稳定性等总体性能指标的提升。