0 引　言

在水下航行体外形设计中，考虑到航行体运动过程中水动特性，存在多种航行体尾部构型方案^[1]。航行体水下运动过程中一般采用高压燃气推动或依靠尾部火箭发动机推动^[2]，不同推进方式产生的气体均会在一定程度上附着在航行体尾部形成尾空泡^{[3 − 4]}。不同尾部外形航行体尾空泡的附着状态不同，对尾空泡上游航行体物面压力分布的影响也有所差异，进而会影响航行体的受力和运动状态^{[5 − 7]}，是水下航行体工程设计的重要关注点。

在航行体尾部附着空泡影响研究方面，目前已有学者利用试验和仿真计算等手段分析了稳态或瞬态空泡演化特征^[8]。在试验方面，结合水洞提供的稳定来流条件，研究了不同来流和通气条件的尾空泡形态变化特征^[9]，也有利用减压水箱，开展变速运动条件下的尾空泡动态演化特征研究^[10]。在数值仿真方面，主要结合仿真计算获得的丰富流动信息，分析航行体水下发射过程中尾空泡收缩膨胀、颈缩拉断和出水溃灭中形态^[11]和压力变化^[12]，并研究尾空泡演化对航行体尾部流场变化以及运动状态的影响^[13]。

针对航行体尾部外形研究，在全湿流条件下，已有学者研究了不同尾部外形航行体发动机水下点火尾部压力分布和推力特性^[14]。在尾部附着空泡状态下，通过数值仿真计算了水下运动过程中不同尾部外形航行体尾空泡形态变化^[15]，分析了尾拍运动时航行体尾部流场结构^[16]，重点关注了尾部外形对航行体阻力^[17]和水动特性变化的影响^[18]。研究大多针对尾部外形对航行体受力和运动状态的影响，对于不同尾部外形航行体尾空泡振荡对上游物面压力影响规律的研究较少。

针对航行体水下运动尾部附着空泡演化问题，本文基于势流理论形成了可计算尾空泡流场的二维轴对称边界元数值计算方法，并利用已有试验数据进行校验。对不同尾部外形航行体尾空泡稳态和瞬态过程进行数值计算，重点关注了尾空泡形态和空泡上游航行体物面压力分布的变化，探究了空泡振荡影响航行体物面压力变化的主要原因，并分析了空泡振荡过程中航行体沾湿区阻力的变化情况。

1 数值计算方法 1.1 控制方程

假设空泡附着在航行体尾部低压扰流区，在无限大的流体域中收缩膨胀，且流动是理想无旋、无粘且不可压缩的，则流体速度势${\varPhi }$满足Laplace方程：

${\nabla}^2{\varPhi } = 0。$

(1)

流体速度：

$V = \nabla {\varPhi }。$

(2)

式中：${{V}}$为速度向量；${{\varPhi }} = {U_\infty }x + \varphi$为总速度势；${{{U}}_\infty }$为来流速度（设来流为$x$轴正向）；$x$为位置坐标；$\varphi$为扰动速度势。

航行体物面满足不可穿透条件，即壁面法向速度为0。

${n}{V}= \frac{\partial \varPhi }{\partial {n}}=0。$

(3)

式中：$n$为航行体物面法向方向向量。

在空泡壁面，流体速度势${{\varPhi }}$、速度$V$及压力$p$满足伯努利方程：

$\frac{{\partial \varPhi }}{{\partial t}} + \frac{1}{2}{V^2} + p/\rho = \frac{1}{2}U_\infty ^2 + {p_\infty }/\rho 。$

(4)

式中：${p_\infty }$为环境压力；$\rho$为流体密度。

1.2 数值方法

在计算方法上本文采用边界元法对速度势进行求解，以扰动速度势$\varphi$为例，航行体物面和尾空泡壁面组成边界Γ，对于边界Γ上的点$p$，可得边界积分方程：

$C(p)\varphi (p) = \int_{{\Gamma }} {\left( {\varphi *\frac{{\partial \varphi }}{{\partial n}} - \varphi \frac{{\partial \varphi *}}{{\partial n}}} \right)} {{\mathrm{d}}} {\mathit{\Gamma }} 。$

(5)

式中：$C(p)$为与边界几何形状相关的常数，本文的边界问题中$C(p) = 1/2$；$\varphi *$为拉普拉斯方程方程基本解，三维情况下$\varphi *(p,q) = - 1/4\text{π} r$，其中$r = \left| {p - q} \right|$。

将${\mathit{\Gamma }}$分成$n$个边界单元${\Gamma } _j$$(j = 1,2, \cdots ,n)$，假设航行体物面上有${n_1}$个单元，尾空泡壁面上有${n_2}$个单元，则${n_1} + {n_2} = n$。对于边界${{\mathit{\Gamma }} _i}$上的节点${p_i}(i = 1,2, \cdots ,n)$，上述积分方程可离散为：

$\sum\limits_{j = 1}^n {{\varphi _j}} {H_{ij}} = \sum\limits_{j = 1}^n {\frac{{\partial {\varphi _j}}}{{\partial n}}} {G_{ij}}。$

(6)

式中：${H_{ij}} = \int_{{\Gamma _j}} {(\partial {\varphi ^*}({p_i},q)/\partial n){{\mathrm{d}}} {{\mathit{\Gamma }}_q}} + 1/2{\delta _{ij}}$，${\delta _{ij}}$为克罗内克函数，仅当$i = j$时${\delta _{ij}} = 1$，其他情况下${\delta _{ij}} = 0$；${G_{ij}} = \int_{{{\mathit{\Gamma }} _j}} {{\varphi ^*}({p_i},q){{\mathrm{d}}} {{\mathit{\Gamma }}_q}}$；$q$为边界单元${\Gamma _j}$$(j = 1,2, \cdots ,n)$上的任意点。

式(6)中需求解$n$个边界节点的速度势和速度势法向导数值，共$2n$个未知数。本文的边界为混合边界，航行体物面满足不可穿透条件，即第二类边界条件，空泡壁面满足第一类边界条件，故边界${{\mathit{\Gamma }} _i}$上$\bar q = \partial {\varphi _i}/\partial n$$(i = 1,2, \cdots ,{n_1})$和${\bar \varphi _i}$$(i = {n_1} + 1,{n_1} + 2, \cdots ,n)$共$n$个值已知，剩余$n$个未知量，可由式(6)离散方程组求得，由此可得各边界单元节点的速度势及其法向导数值，计算可求得各节点的速度和压力，从而确定整个流场的状态。

1.3 空泡尾部闭合模型

实际流动中航行体尾部流场复杂多变，计算过程中需要对尾空泡进行简化处理，现阶段针对航行体尾部空泡演化问题，已有多种尾空泡简化闭合模型。考虑到数值计算方法的稳定性，本文采用了圆弧壁面条件对空泡尾部进行了简化处理，尾空泡闭合模型示意图如图1所示。

当航行体尾空泡尾部切线与水平方向夹角$\eta$小于某个特定角度时，采用圆弧对空泡后续形状进行简化处理，从而使尾空泡闭合。在圆弧段尾空泡仍然满足运动学边界条件，但是空泡壁面上的速度势法向导数值不再是${q_c}$，而是满足给定的分布规律：

$\frac{{\partial \varphi }}{{\partial n}} = \frac{{\partial \varphi (\theta )}}{{\partial n}}\left( {1 - \sin \left(\frac{\text{π} }{2}\left(\frac{{{\theta _0} - \theta }}{{{\theta _0}}}\right)\right)} \right) 。$

(7)

式中：${\theta _0}$为空泡尾部整个圆弧段对应的夹角；$\theta$为圆弧上的点与闭合$B$点之间对应圆弧的夹角。

1.4 数值计算有效性验证

为验证本文采用数值计算方法的有效性，航行体尾部附着空泡状态下，对水下运动航行体物面压力进行计算，并将计算结果与文献[13]中试验数据进行对比，结果如图2所示。图中，压力系数${C_p} = \left( {p - {p_\infty }} \right)/ (1/2\rho U_\infty ^2)$，$p$为航行体表面压力，${U_\infty }$为来流速度。

由图可知，在附着空泡状态下，航行体物面压力的数值计算结果与水洞试验数据拟合程度均较好，说明本文所采用的数值计算方法对航行体尾空泡的相关计算可靠。

2 不同外形附着空泡稳态分析 2.1 稳态空泡形态

采用上述数值方法，对不同尾部外形航行体进行计算，航行体外形示意图如图3所示。图中，A点为航行体尾锥顶点，$\lambda$为航行体尾部半锥角，$\lambda = {0^ \circ }$为航行体尾部无锥角情况，$L$为航行体长度，$D$为航行体直径。

航行体尾部半锥角$\lambda = {6^ \circ }$时，对不同来流情况下航行体尾部附着空泡进行迭代计算，获得空泡收敛后的形态如图4所示。图中，$X/L$和$Y/L$分别为无量纲的轴向和径向坐标。其中，空化数$\sigma = ({p_\infty } - {p_c})/(1/2\rho U_\infty ^2)$，${p_c}$为尾空泡泡内压力。

稳态情况下空泡形态与空化数相关，对于不同初始条件，只要空化数相同，收敛后的空泡形态一致。当环境压力${p_\infty }$减小、空泡内压力${p_c}$增大或来流速度${U_\infty }$增大时，空化数降低，空泡初始膨胀能力增强，导致空泡轴向长度增加，空泡几何尺寸增大。

航行体尾部外形也会影响附着空泡的形态，在空化数$\sigma = 0.338$条件下，对不同尾部外形航行体尾空泡进行计算，获得空泡形态如图5所示。

较大半锥角$\lambda$条件下，航行体尾部截面尺寸更大，尾空泡初始径向尺寸明显增加。同时随着半锥角$\lambda$增加，流体在尾部与壁面分离后轴向流动速度增大，尾空泡具有更大初始膨胀能量。以上两方面因素导致，在相同空化数条件下，随着航行体尾部半锥角$\lambda$增大，尾空泡几何尺寸更大。

2.2 物面压力变化

在尾空泡稳定状态下，不同尾空泡形态会导致航行体尾部流场速度分布产生差异。不可压缩流动中，航行体尾部流动结构改变会对空泡上游产生影响，直接表现是航行体物面压力发生显著变化。不同航行体尾部外形条件下，计算空化数变化时航行体物面压力变化情况如图6所示。图中压力变化量$\Delta {C_p} = {\left. {{C_p}} \right|_{\sigma = 0.226}} - {\left. {{C_p}} \right|_{\sigma = 0.667}}$。

空化数变化时，航行体尾空泡形态改变，航行体物面压力随之发生变化。计算结果表明，空化数减小时，尾空泡几何尺寸增大，导致航行体物面压力增大，越靠近空泡的航行体物面区域压力变化越大，航行体头部物面压力基本不变。不同尾部外形也会影响航行体物面压力变化，相同空化数变化条件下，随着尾部半锥角$\lambda$的增大，航行体物面压力变化量增大。

为分析航行体尾部外形对空泡上游物面压力的影响，仅关注尾部锥角引起的压力变化，在空化数$\sigma = 0.338$条件下，计算得到不同尾部外形航行体物面压力相对航行体尾部无锥角情况的变化量如图7所示，图中压力变化量$\Delta {C_p} = {\left. {{C_p}} \right|_\lambda } - {\left. {{C_p}} \right|_{\lambda = 0}}$，${\left. {{C_p}} \right|_{\lambda = 0}}$为航行体无锥角时物面压力。

相同空化数条件下，相比无锥角情况，航行体尾部锥角会使空泡几何尺寸变大，导致临近尾空泡的航行体物面区域压力增大，尾锥顶点A处压力变化最大。半锥角$\lambda$较大的航行体物面压力变化量较大，说明半锥角增大会导致尾空泡对上游航行体物面压力分布影响程度更大。

3 不同外形附着空泡瞬态分析 3.1 尾空泡振荡现象

尾空泡动态演化过程会影响航行体尾部流场，从而导致空泡上游航行体物面压力变化，为研究不同航行体尾部外形空泡振荡影响，在稳态空泡的计算基础上，引入扰动压力使空泡振荡，从而计算尾空泡的动态演化过程。在空化数$\sigma = 0.338$的初始条件下，减小尾空泡内压力，使空化数增大为$\sigma ' = 0.564$，空泡稳定状态被破坏开始振荡。空泡振荡过程中，不同尾部外形航行体尾空泡泡内压力变化如图8所示。图中压力系数${C_p} = ({p_c} - {p_\infty })/(1/2\rho U_\infty ^2)$，无量纲时间$t^* = {U_\infty }t/D$。

尾空泡振荡过程中泡内压力也随之周期性变化，相同初始空化数条件下，航行体尾部半锥角$\lambda$越大，泡内压力振荡峰值越大，振荡周期越长，说明半锥角增大会导致空泡振荡泡内压力变化更剧烈。

图8中$\Delta T$为尾空泡泡内压力从振荡开始到达到第一个峰值所需时间，空泡振荡过程中空泡形态也随之变化，故$\Delta T$同时也为空泡第一次收缩至最小所需时间。为研究尾空泡振荡过程中空泡形态变化，针对空泡第一次收缩过程，计算半锥角$\lambda = {6^ \circ }$时典型时刻空泡形态如图9所示，其中${t_0}$为初始时刻。

由图可知，尾空泡收缩过程中收缩速度是变化的，${t_0} \sim {t_0} + 0.25\Delta T$时间段，尾空泡收缩初期，尾空泡收缩速度较小，尾空泡形态变化不大，随着尾空泡壁面收缩速度增加，尾空泡形态差异逐渐明显。

3.2 物面压力变化

尾空泡振荡会改变航行体尾部流场结构，进而影响空泡上游航行体物面压力，导致航行体物面压力随空泡振荡而变化。以航行体尾锥顶点A点为例，在相同初始空化数条件下，计算不同尾部外形航行体尾空泡振荡过程中A点压力变化如图10所示。图中压力变化量$\Delta {C_p} = {\left. {{C_p}} \right|_{t^*}} - {\left. {{C_p}} \right|_{t^* = 0}}$。

尾空泡振荡过程中，A点压力也随之变化，航行体尾部半锥角$\lambda$越大，A点压力变化越剧烈，振荡峰值越大周期越长。对比图10与图8可知，A点压力变化与空泡振荡过程中泡内压力变化趋势大致相同，进一步说明了空泡振荡会导致航行体物面压力变化。

为进一步研究尾空泡振荡对上游航行体物面压力的影响，针对空泡第一次收缩过程，计算不同尾部外形航行体尾空泡收缩过程中压力变化$\Delta {C_p} = {\left. {{C_p}} \right|_{t^* = \Delta T}} - {\left. {{C_p}} \right|_{t^* = 0}}$，结果如图11所示。其中${\left. {{C_p}} \right|_{t^* = \Delta T}}$为空泡第一次收缩到最小时航行体物面压力。

尾空泡收缩过程中，泡内压力增大，空泡上游航行体物面压力随之增大，越靠近空泡的航行体物面区域压力变化越大。随着半锥角$\lambda$增大，航行体物面压力变化量增加，与前文结论一致。尾空泡收缩过程中形态变化会使周围流体质点速度和加速度变化，影响航行体尾部流场结构，进而导致上游航行体物面压力改变。尾空泡振荡过程中，航行体尾部流场压力变化满足式(4)伯努利方程，可得：

${C_p} = 1 - \frac{q}{{{q_\infty }}} - \frac{{\rho \dot \Phi }}{{{q_\infty }}}。$

(8)

式中：$q = 1/2\rho {v^2}$为流体的动压；${q_\infty } = 1/2\rho U_\infty ^2$为来流动压；$\dot \Phi = \partial \Phi /\partial t$为流体速度势变化率。

流场中流体质点速度改变会导致动压变化。流体质点加速度改变会导致流体速度势变化率发生变化，故可用动压${C_{p1}} = 1 - {q_\infty }/q$和速度势变化率${C_{P2}} = - \rho \dot \Phi / {q_\infty }$分别表征流场中流体质点的速度与加速度变化，流场中压力${C_P} = {C_{P1}} + {C_{P2}}$。为分析尾空泡振荡影响上游航行体物面压力变化的主要原因，分别计算空泡第一次收缩过程中航行体物面动压和速度势变化率的变化情况，结果如图12和图13所示。图中$\Delta {C_{p1}} = {\left. {{C_{p1}}} \right|_{t^* = \Delta T}} - {\left. {{C_{p1}}} \right|_{t^* = 0}}$，$\Delta {C_{p2}} = {\left. {{C_{p2}}} \right|_{t^* = \Delta T}} - {\left. {{C_{p2}}} \right|_{t^* = 0}}$。

尾空泡收缩过程中，由图11可知，空泡上游航行体物面压力增大，图12中航行体物面动压减小，与物面压力变化趋势相反，图13中航行体物面速度势变化率增大，与物面压力变化趋势一致，且不同半锥角情况下变化规律一致，说明航行体物面压力变化的主要原因是速度势变化率改变。空泡振荡过程中空泡形态变化导致周围流体质点加速度变化，流体加速度变化使流场速度势变化率改变，进而导致航行体物面压力变化。随着航行体尾部半锥角$\lambda$增大，尾空泡振荡更剧烈，对空泡上游航行体物面压力变化的影响也更大。

3.3 沾湿区阻力变化

尾空泡振荡导致上游航行体物面压力产生变化，压力扰动会影响航行体沾湿区的受力情况，进而影响航行体的运动状态。针对尾空泡振荡过程，在相同初始空化数条件下，对不同尾部外形航行体尾空泡流场进行计算，获得航行体沾湿区阻力变化情况如图14所示。图中$\Delta {C_d} = {\left. {{C_d}} \right|_{t^*}} - {\left. {{C_d}} \right|_{t^* = 0}}$，${C_d}$为航行体沾湿区阻力系数：

${C_d} = \frac{{\displaystyle\int_{{S_w}} {(p - {p_\infty }} ){n_x}{\mathrm{d}}s}}{{\dfrac{1}{2}\rho U_\infty ^2{S_{ref}}}} 。$

(9)

式中：${S_{ref}} = 1/4{\text π} {D^2}$为航行体沾湿区参考面积：${n_x}$为法向向量在X轴方向上的分量。

尾空泡振荡过程中，航行体沾湿区阻力也随之变化，相同初始条件下，航行体尾部半锥角$\lambda$越大，航行体沾湿区阻力变化越剧烈，振荡峰值越大振荡周期越长。对比图14与图8、图10可知，空泡振荡过程中，航行体沾湿区阻力变化趋势与航行体物面A点压力变化、空泡泡内压力变化基本一致。

由此可得尾空泡振荡影响航行体运动的过程，空泡振荡泡内压力变化的同时空泡形态随之改变，空泡形态变化会使周围流体质点加速度改变化，从而影响空泡上游航行体物面压力，航行体物面压力变化导致航行体沾湿区阻力变化，进而影响航行体的受力和运动状态。航行体尾部半锥角$\lambda$越大，尾空泡振荡越剧烈，对航行体运动的影响越大。

4 结　语

1）稳态空泡计算结果表明，相同空化数条件下，随着航行体尾部半锥角$\lambda$增大，尾空泡几何尺寸增加，空泡对上游航行体物面压力影响更剧烈。

2）尾空泡振荡过程中，随着航行体尾部半锥角$\lambda$的增大，空泡振荡加剧，对航行体物面压力影响程度更大。尾空泡振荡影响空泡上游航行体物面压力的主要原因是空泡振荡形态变化使周围流体质点加速度变化。

3）尾空泡振荡过程中，航行体尾部半锥角$\lambda$越大，航行体阻力变化越剧烈，对航行体运动影响越大。空泡振荡导致航行体物面压力变化进而使航行体沾湿区阻力改变，空泡振荡泡内压力变化、航行体物面压力变化和沾湿区阻力变化趋势大致相同。