泵喷推进器由导管、定子和转子构成，定转子的组合可以实现机械能的高效转化，在导管的保护和整流的作用下，转子的效率会得到进一步的提高，同时导管可以有效地屏蔽和降低辐射噪声。泵喷具有体积小、效率高、噪声低等特点，受到了国内外学者的关注和重视，研究人员在工作机理和流场特性等方面展开了大量的研究，但主要为试验研究^[1-2]和CFD数值仿真研究^[3-7]。

基于势流的水动力预报技术在船舶及水下航行体的推进器理论研究和设计工作中占有重要地位，但由于泵喷独特的结构形式和复杂多变的内部流场环境，使得基于势流方法的泵喷水动力性能预报精度并不理想，在相关的研究文献中大多是以结构相似、外形简化的带定子导管桨^[8-13]为主，并没有针对泵喷推进器的水动力预报方法。对于装配在鱼雷或潜艇的泵喷推进器，由于受到尾锥体倾斜壁面的限制，转子和定子的端面会随之倾斜，因而以往研究中基于圆柱面的网格划分方法将无法用于这种几何外形的叶片。另外，泵喷的结构更为紧凑，转子和导管的缝隙区域存在强烈的间隙流动和泄露涡^{[14, 15]}，而以往对带定子导管桨的研究中并没有考虑间隙泄露涡带来的影响。

本文以定常面元法为基础，建立了泵喷推进器各部分组件间的干扰数值模型，并提出了适用于平顶式叶片的梢部泄露涡模型，通过计算结果的对比验证方法的准确性。并基于此方法，从推进性能指标、组件的压力分布和流场的速度分布等方面，系统地探讨导管翼型剖面的拱度对泵喷水动力性能的影响。

1 泵喷计算模型

本文采用以速度势为基础的低阶面元法，物体表面的离散单元为双曲面元，面元上布置的奇点为源汇和偶极的混合分布。采用求解有升力体问题的基本方程建立泵喷推进器的数值模型，并在原模型的基础上加入梢部泄露涡模型。

1.1 干扰数值模型

将泵喷推进器的部件分为在同一个坐标系下的两个独立存在的有升力体，分别为转子-桨毂系统和定子-导管系统，如图 1所示。两个系统分别组成基本方程，系统间的相互影响通过对彼此系统周围流场的扰动进行考虑，而干扰作用则是以诱导速度的形式来表达。依次求解各系统速度势方程和诱导速度方程，通过诱导速度循环迭代的方式，使两个系统的数值解都达到收敛标准，得到泵喷推进器整体的计算结果。以下是速度势和诱导速度的表达方法。

假定泵喷工作于速度为V₀的无旋、非粘性、不可压缩来流中，根据格林公式^[16]，定子-导管和转子-桨毂在边界面上的速度势的求解方程分别为

$ \begin{array}{*{20}{c}} {2{\rm{\pi }}{\varphi ^k}\left( P \right) = \iint\limits_{{S_{\text{s}}} + {S_{\text{d}}}} {{\varphi ^k}\left( Q \right)\frac{\partial }{{\partial {n_Q}}}\left( {\frac{1}{{{R_{PQ}}}}} \right){\text{d}}S} + } \\ {\iint\limits_{{S_{{\text{ws}}}} + {S_{{\text{wd}}}}} {\Delta {\varphi ^k}\left( {{Q_1}} \right)\frac{\partial }{{\partial {n_{{Q_1}}}}}\left( {\frac{1}{{{R_{P{Q_1}}}}}} \right){\text{d}}S} + } \\ {\iint\limits_{{S_{\text{s}}} + {S_{\text{d}}}} {\left( {\left( {{\mathit{\boldsymbol{V}}_0} + \mathit{\boldsymbol{V}}_{{\text{sd}},{\text{rh}}}^k} \right) \cdot {\mathit{\boldsymbol{n}}_Q}} \right)\left( {\frac{1}{{{R_{PQ}}}}} \right){\text{d}}S}} \end{array} $

(1)

$ \begin{array}{*{20}{c}} {2\pi {\varphi ^k}\left( P \right) = \iint\limits_{{S_{\text{r}}} + {S_{\text{h}}}} {{\varphi ^k}\left( Q \right)\frac{\partial }{{\partial {\mathit{\boldsymbol{n}}_Q}}}\left( {\frac{1}{{{R_{PQ}}}}} \right){\text{d}}S} + } \\ {\iint\limits_{{S_{{\text{wr}}}}} {\Delta {\varphi ^k}\left( {{Q_1}} \right)\frac{\partial }{{\partial {n_{{Q_1}}}}}\left( {\frac{1}{{{R_{P{Q_1}}}}}} \right){\text{d}}S} + } \\ {\iint\limits_{{S_{\text{r}}} + {S_{\text{h}}}} {\left( {\left( {{\mathit{\boldsymbol{V}}_0} + r{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_\theta } + \mathit{\boldsymbol{V}}_{{\text{rh,sd}}}^k} \right) \cdot {\mathit{\boldsymbol{n}}_Q}} \right)\left( {\frac{1}{{{R_{PQ}}}}} \right){\text{d}}S}} \end{array} $

(2)

式中：S_s、S_d、S_r和S_h分别为定子表面、导管表面、转子表面及桨毂表面，S_ws、S_wd和S_wr分别为对应部件的尾涡面；P和Q为边界面S上的任意控制点；R_PQ为P、Q两点间的直线距离；n为相应物面的单位法向量，指向流场内；Δφ为通过尾涡面的速度势跳跃；Ω_θ为螺旋桨的旋转角速度；上角标k表示第k次迭代运算；V_{sd, rh}表示转子-桨毂系统在定子和导管的表面产生的诱导速度，作为定子-导管系统相对进流速度的一部分，V_{rh, sd}则相反，其表达方程为

$ \begin{array}{*{20}{c}} {4{\rm{\pi }}\mathit{\boldsymbol{V}}_{{\text{sd}},{\text{rh}}}^k\left( P \right) = \iint\limits_{{S_r} + {S_h}} {{\varphi ^k}\left( Q \right){\nabla _P}\frac{\partial }{{\partial {n_Q}}}\left( {\frac{1}{{{R_{PQ}}}}} \right){\text{d}}S} + } \\ {\iint\limits_{{S_{wr}}} {\Delta {\varphi ^k}\left( {{Q_1}} \right){\nabla _P}\frac{\partial }{{\partial {n_{{Q_1}}}}}\left( {\frac{1}{{{R_{P{Q_1}}}}}} \right){\text{d}}S} + } \\ {\iint\limits_{{S_r} + {S_h}} {\left( {\left( {{\mathit{\boldsymbol{V}}_0} + r{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_\theta } + \mathit{\boldsymbol{V}}_{{\text{rh,sd}}}^{k - 1}} \right) \cdot {\mathit{\boldsymbol{n}}_Q}} \right){\nabla _P}\left( {\frac{1}{{{R_{PQ}}}}} \right){\text{d}}S}} \end{array} $

(3)

$ \begin{array}{*{20}{c}} {4{\rm{\pi }}\mathit{\boldsymbol{V}}_{{\text{rh}},{\text{sd}}}^k\left( P \right) = \iint\limits_{{S_s} + {S_d}} {{\varphi ^k}\left( Q \right){\nabla _P}\frac{\partial }{{\partial {n_Q}}}\left( {\frac{1}{{{R_{PQ}}}}} \right){\text{d}}S} + } \\ {\iint\limits_{{S_{{\text{ws}}}} + {S_{{\text{wd}}}}} {\Delta {\varphi ^k}\left( {{Q_1}} \right){\nabla _P}\frac{\partial }{{\partial {n_{{Q_1}}}}}\left( {\frac{1}{{{R_{P{Q_1}}}}}} \right){\text{d}}S} + } \\ {\iint\limits_{{S_{\text{s}}} + {S_{\text{d}}}} {\left( {\left( {{\mathit{\boldsymbol{V}}_0} + \mathit{\boldsymbol{V}}_{{\text{sd,rh}}}^{k - 1}} \right) \cdot {\mathit{\boldsymbol{n}}_Q}} \right){\nabla _P}\left( {\frac{1}{{{R_{PQ}}}}} \right){\text{d}}S}} \end{array} $

(4)

初次求解时，如果首先计算转子和桨毂，V_{rh, sd}⁰为0；否则V_{sd, rh}⁰为0。式(1)~(4)为各系统速度势和诱导速度的积分方程，根据图 1中划分的网格单元，离散积分方程，形成由偶极和源汇构成的离散方程组。再将各系统的离散方程组满足相应部件在随边处的等压库塔条件来进行求解，等压库塔条件的表达式如下

$ {\left( {\Delta p} \right)_{{\rm{TE}}}} = P_{{\rm{TE}}}^ + - P_{{\rm{TE}}}^ - = 0 $

(5)

1.2 梢部泄露涡模型

文献[14-15]表明，泄露涡的涡线通常是从叶面顶端发出，经过间隙区域的壁面约束，在叶背一侧发展完全。涡线的螺距角随着轴向距离的增加而不断减小，直至在叶片随边的下游某处，与尾涡涡线的螺距角相近。基于上述结论，针对所采用的泵喷推进器，本文初步建立了转子梢部泄露涡的几何模型，如图 2所示。

涡线从转子顶端面与叶面的连接点发出，各涡线的初始螺距角为涡线起始点的水动力螺距角β₀，随着轴向距离的增加螺距角β不断减小，直到下游某处等于转子梢部叶剖面的几何螺距角β_g，然后以固定螺距β_g继续沿下游延伸。为简化梢涡面元的划分方法，本文假定各涡线的螺距角在同一轴向位置T点趋于相同，设T点的轴向坐标为x_t，如图 2所示。涡线的径向坐标变化趋势与导管内壁面一致，可通过插值获得，则涡线上第i个单元的几何坐标可表达为

$ \left\{ \begin{array}{l} {x_i} = {x_{i - 1}} + \Delta x\\ {r_i} = f\left( {{x_i}} \right)\\ {\theta _i} = {\theta _{i - 1}} + \left( {{x_i} - {x_{i - 1}}} \right)/{r_i}/\tan \left( {{\beta _i}} \right) \end{array} \right. $

(6)

$ {\beta _i} = \left\{ {\begin{array}{*{20}{c}} \begin{array}{l} {\beta _0} + \left( {{\beta _g} - {\beta _0}} \right)\left( {i - 1} \right)/{\rm{nu}}{{\rm{m}}_T},\\ {\beta _g}, \end{array}&\begin{array}{l} 0 \le {x_i} \le {x_T}\\ {x_i} \ge {x_T} \end{array} \end{array}} \right. $

(7)

式中num_T为从起始点至T点的单元个数。

将梢部泄露涡模型作为转子尾涡面的一部分，即转子叶面和顶端面的相接面元、叶面和叶背在随边处的相接面元同时满足等压库塔条件。

1.3 计算迭代过程

初始计算时选取的系统对计算过程的收敛性和结果精度没有影响，这里以转子-桨毂系统为例，为简化描述，两个系统分别为转子系统和定子系统。首先在不考虑定子系统影响的情况下，计算转子系统的敞水水动力性能，然后将转子系统对定子系统产生的诱导速度作为定子和导管来流的一部分，并计算定子系统的水动力性能和对转子和桨毂的诱导速度，再重新计算转子系统的水动力性能(计及定子系统的干扰)，如此迭代计算直至泵喷整体的推力系数收敛为止，图 3为计算迭代流程。

泵喷推进器各部件和整体的性能指标表达式如下所示:

$ {K_{{\rm{Tr}}}} = \frac{{{T_{\rm{r}}}}}{{\rho {n^2}{D_{\rm{r}}}^4}};{K_{{\rm{Qr}}}} = \frac{{{Q_{\rm{r}}}}}{{\rho {n^2}{D_{\rm{r}}}^5}};{K_{{\rm{Ts}}}} = \frac{{{T_{\rm{s}}}}}{{\rho {n^2}{D_{\rm{r}}}^4}}; $

$ {K_{{\rm{Qs}}}} = \frac{{{Q_{\rm{s}}}}}{{\rho {n^2}{D_{\rm{r}}}^5}};{K_{{\rm{Td}}}} = \frac{{{T_{\rm{d}}}}}{{\rho {n^2}{D_{\rm{r}}}^4}};\eta = \frac{{{K_{T - {\rm{all}}}}}}{{{K_{{\rm{Qr}}}}}}\frac{J}{{2{\rm{ \mathsf{ π} }}}} $

式中:K_T和K_Q分别为推力系数和转矩系数，下角标r、s、d分别对应转子、定子和导管，总推力系数分别为K_T-all=K_Tr+K_Ts+K_Td，平衡力矩系数K_Q-all=K_Qr+K_Qs；η为泵喷的效率；进速系数J=V₀/nD_r；D_r为转子最大直径；n为转子转速。

2 方法的准确性验证

由于泵喷推进器的用途较为特殊，缺乏公开的试验数据，因此利用较为成熟的CFD仿真技术代替试验手段，对本文提出的泵喷预报方法进行验证和计算结果分析。

2.1 CFD方法可行性分析

为了验证CFD方法用于计算带定子的导管类推进器数值仿真的适用性，本文另外选取具有试验数据的前置定子导管桨进行CFD方法的精度验证。几何模型如图 4所示，包含9个定子叶片和4个转子叶片，转子直径D为0.254 m，详细信息参见文献[8]。

2.1.1 计算域及网格划分

计算流场分为模拟无限远边界和模拟转子旋转的内域两个柱形计算域。外域直径取15D以消除壁面效应对旋转区域的影响，同时圆柱面设置为对称面，速度入口及压力出口位置距导管前缘均为10D。内域仅包含转子，外边界为导管内壁面，上游边界取导管入口处，下游边界位于定转子之间，内外域通过交互面传递数据，如图 5所示。

计算域和近壁面分别采用笛卡尔网格和棱柱网格进行离散。考虑到边界层的影响，定转子表面第一层厚度设置为0.05 mm，对应的Y⁺值在60左右。为了控制Y⁺值和网格数量，布置至少3层棱柱网格，空间网格以1.2倍倍率增长。各部件的边线、转子顶面以及梢部间隙区域进行加密处理，如图 6所示。选取RNGk-ε湍流模型封闭计算方程，定常性能预报采用多重旋转坐标系(MRF)模型。扩散相采用中心差分格式，对流相采用一阶迎风格式，压力速度耦合采用SIMPLEC算法。

2.1.2 网格无关性分析

依据Roy^[17]对笛卡尔网格在粘流计算中的不确定度研究方法，建立多套不同数量的网格。网格尺寸细化率r与网格总数满足如下关系：

$ r = {\left( {\frac{{{N_{{\rm{fine}}}}}}{{{N_{{\rm{coarse}}}}}}} \right)^{1/d}} $

式中：N表示网格总数，d表示问题的空间维数。各套网格在整个计算域中沿三个方向的网格细化率保持一致。本文将r取值为1.2，得到了网格总数为200万、350万、611万的网格。由表 1列出的结果显示，随着网格数的增加，推力和转矩的计算值与试验值的误差减小，当网格数大于350万时，具有良好的精度，K_t-all和K_q的误差均低于4%。

图 7为利用611万网格预报的全工况水动力性能结果，CFD的计算结果在设计工况范围内与试验值基本吻合，仅在进速系数较高时存在偏差。以上结果表明CFD方法在带定子导管类推进器的数值模拟中具有较高的精度和可信度，可用于泵喷推进器的敞水性能计算和后续的势流方法验证。

2.3 势流方法准确性验证

本文选取的势流方法研究对象为后置定子泵喷推进器，计算模型包含11个转子叶片和9个定子叶片，转子最大直径D_r为0.238 m，转速n为2 000 r/m，转子顶端与导管的最小间隙为3 mm。计算模型的网格划分见图 1，转子和定子单个叶片的网格数均为32×10，单通道的桨毂和导管的面元数分别为56×3、60×9。本文共计算了4种工况，进速系数J分别为1.29、1.94、2.59、3.24。将泵喷推进器定子、转子以及整体的性能计算值与CFD计算结果进行对比，如图 8和图 9所示。

由图 8和图 9计算结果分析可得：梢涡模型的引入提高了定转子以及整体性能的预报精度，推力、转矩以及效率随进速的变化趋势均与粘流结果保持一致，吻合良好。在进速系数较低时，梢涡模型对结果的改善程度更加明显，间接说明当推进器处于重载时，转子梢部存在较为剧烈的泄露涡，梢涡模型可以准确地模拟泄露涡对周围流场的扰动情况。在研究中发现，泄露涡模型仅对接触面元的影响较大，转子和定子均与梢涡相接或相交，因此水动力性能变化明显，而导管面元并不与梢涡模型接触，因此导管推力的变化较小。从泵喷整体的效率曲线图可以发现，有无梢涡模型的结果相近。这是因为整体效率的变化趋势主要取决于转子的效率曲线。但从转子的结果来看，梢涡模型对转子效率的影响较小。另外，在进速较高时，转子推力系数和效率的计算结果偏高，这是由于本文采用的梢涡和尾涡均为线性结构，简化了泄露涡的真实流动。

3 导管拱度的变尺度分析

本节沿用导管原翼型剖面进行分析，为了降低导管几何变化对间隙流动、转子外形带来的影响，在改变拱度的过程中尽可能使转子间隙区域保持不变，在保留当前翼型剖面的攻角、弦线沿轴向的长度、沿弦长方向的厚度分布以及最大厚度处的空间坐标的前提下，按照各弦向位置的厚度t的百分比来改变在相应位置的拱度分布，如图 10所示。当拱度f≤0时，导管内壁面会加速内部流体的通过，因此表现为加速型导管；反之当拱度大于零时，表现为减速型导管。

由图 11泵喷效率和转子效率可看出，导管拱度的增加拓宽了泵喷的有效做功范围，提高泵喷整体和转子的效率峰值，并且使两者的最高效率点向高速工况偏移。当泵喷采用减速型导管时，拱度变化对η_r和η曲线的影响要小于加速导管，特别在1.29≤J≤2.59时，整体效率的变化幅度很小，只在J=3.24时随着拱度的增加而略有提升。由于泄露涡结构的线性假设，转子效率在高进速时的计算结果偏高，但从不同拱度的变化趋势可以发现，相同转子在不同拱度的导管内运转，不仅最佳效率点相差较大，而且效率峰值存在明显差异。这是因为拱度变化导致转子来流和过流的径向速度分布发生变化，进而影响了转子的径向载荷分布。

以来流速度为例，图 12为在J=3.24工况下导管-定子系统对导管入口盘面各径向位置的轴向诱导速度V_{x-in, sd}，转子来流速度的轴向分量为相应径向位置的V₀和V_{x-in, sd}的总和。从图中可知，V_{x-in, sd}在整个导管入口的分布方式大致相同，但由于拱度越小，导管的入口直径越大，转子的来流区域占整个入口盘面的比例越小，转子来流的径向速度分布发生改变，例如转子叶梢位置的来流速度在f=1.0 t的结果中处于急剧下降阶段，而当f=-0.5 t时，处于缓慢减小阶段。增加拱度降低了入口和转子绕流的速度，由伯努利方程可知，转子的压力系数随之增加，有助于延迟叶片的空化发生。

结合J=3.24的不同拱度的转子径向环量分布如图 13所示，转子的来流分布对叶片径向的载荷分布影响较大，随着拱度的增加，转子在各径向处的环量值增大，并且径向位置越高提升幅度越大，梢部由此逐渐变为做功的主要部位。

拱度变化同样对导管性能影响较大，如图 14所示，推力系数为负，说明导管表现为阻力，当进速较低时，低拱度的导管阻力偏小，随着进速系数的增加，所有拱度的K_Td逐渐减小，但加速导管的曲线斜率要大于减速导管，当J>3.24时，加速导管的阻力将大于减速导管。并且加速导管的拱度越小，阻力增长速度越快；而对于减速导管，拱度变化对K_Td曲线的下降趋势影响较小。加速导管和减速导管的阻力曲线存在差异是因为影响轴向力变化的主要贡献区域不同。

图 15和图 16分别为f=-0.5 t和f=1.0 t导管在θ=20°处纵向剖面的压力系数C_p及其轴向力-Cp·n_x，其中n_x为物面单位法向量的轴向分量。由图 15可知，加速导管的内壁面压力系数明显低于外壁面，但只有转子范围的内壁面对轴向力的贡献度较大，而定子范围的内壁面的n_x接近于零，如图 10的轮廓图所示，因此轴向力很小。轴向力随进速的主要变化区域为导管前缘区域和内壁面的转子区域，且整体的-C_p·n_x均与进速呈正比关系。由图 16的结果可得，减速导管的压力系数在前缘和定子范围的内壁面对进速的敏感度要高于加速导管。轴向力的主要变化区域与加速导管不同，为外壁面的前缘以及内壁面的后半段，且两个区域随进速的变化趋势相反，因此在当前攻角下，减速导管的阻力增长速度要远小于加速导管。

由以上的变拱度的计算结果可以发现，增加拱度有利于提高泵喷的推进效率和有效做功范围，f≥0.5 t的导管在高速工况下的水动力性能表现要优于f < 0.5 t的导管，但是当f≥0.5 t时，导管阻力持续增大，导致泵喷的效率增幅并不明显。因此对于装备在高速航行体倾斜尾锥面的泵喷推进器，应选取拱度为f≥0.5 t的减速导管，以保证较高的推进效率和良好的高速适应性，同时在推进性能满足设计要求的前提下，应尽量减小导管内壁面的弧度以兼顾导管加工、定转子叶型匹配以及流道宽度变化引起的轴系激振力等问题。

4 结论

1) 泵喷内部流场存在较为强烈的间隙泄露涡，影响转子梢部以及下游定子的绕流场，梢涡模型可以模拟间隙泄露涡对周围流场的扰动作用，从而提高转子和定子的推力和转矩等性能指标的预报精度，计算结果的准确性验证表明了本文提出的泵喷水动力性能预报方法的有效性。

2) 导管拱度的增加拓宽了泵喷的有效做功范围，提高泵喷整体和转子的效率峰值，并且使两者的最高效率点向高速工况偏移。当导管表现为减速型时，拱度变化对效率的影响减弱。

3) 拱度变化导致入口及转子绕流的径向速度分布发生变化，影响转子的径向载荷分布。随着拱度的增加，降低了内流场的轴向速度，转子叶片的压力值随之增大，空化性能提升；同时转子梢部的做功能力增强，从而扩大了叶片的有效做功面积。

4) 拱度变化影响导管的水动力特性。拱度越大，导管的阻力越大，但由于加速导管的阻力性能对进速系数的变化比减速导管更为敏感，因此随着进速的增加，加速导管的阻力将逐渐超过减速导管。

5) 对于装备在高速航行体倾斜尾锥面的泵喷推进器，应选取减速型导管，并在满足设计要求的前提下，尽量减小导管内壁面的弧度以兼顾导管加工、定转子叶型匹配以及流道宽度变化引起的轴系激振力等问题。

由于篇幅和时间所限，本文只研究了导管的拱度对泵喷性能和流场特征的影响，导管的其他几何参数，如间隙尺寸、攻角、长细比等，将作为后续的研究工作，为导管几何参数的选择提供更为详细的指导性建议和理论支撑。