引言

混合式CRP (contra-rotating propeller)是由两个对转的螺旋桨和一个吊舱组成，一般而言，吊舱位于螺旋桨的下游，为后桨提供电力推进.混合式CRP具有较高的推进效率及良好的系统冗余性，且船体所需推力由两个螺旋桨分担，改善了推进系统的振动、噪声性能.正因为其各方面优良的性能，被挪威船级社与德国劳氏船级社展望^[1]为2020年"绿色船舶"的一大亮点.混合式CRP最早由ABB (asea brown boveri)公司于2001年推出，随后应用于日本的两艘客滚船^[2].但混合式CRP组成构件较多的特点，也使得对它的试验及数值研究变得困难.在混合式CRP推出之初，大部分研究着重于试验研究^[3-5]，而研究的主要内容是混合式CRP的试验方法及试验步骤，比较典型的有Cheng等^[6]、Sasaki等^[7]和Quereda等^[8]的工作，正是基于这些学者的研究成果，2014年的第27届国际拖曳水池会议推进委员会提出了混合式CRP的初步试验规程^[9].在数值方法方面，混合式CRP的研究还不算深入，并且大部分是基于黏流的手段^[10-11]，这主要是由于势流方法求解组合推进器大多是基于迭代方法^[12-13]，构件之间的影响通过各自在对方位置处的诱导速度或诱导速度势进行考虑.这种迭代方法广泛应用于由两个构件组成的推进器，例如吊舱推进器^[14]和对转桨^[15-16]等，但若应用于混合式CRP这类由3个构件组成的推进器就会比较复杂.因此，一直以来基于势流方法对混合式CRP的研究比较少.尽管黏流方法能够更精细地捕捉流场^[17]，但势流方法具有较高的求解效率，在推进器理论设计方面具有重要意义.为了有效的开展混合式CRP设计分析，促进混合式CRP的工程应用，有必要对混合式CRP的势流计算方法进行研究.

本文选用基于势流理论的面元法作为数值手段，开展混合式CRP的面元法计算方法研究，研究中建立了迭代面元法和整体面元法，对比了两者的计算效果，为混合式CRP的设计分析，提供高效可靠的数值手段.

1 数值方法 1.1 迭代面元法

已有研究表明^[18-20]，基于诱导速度势的低阶面元法能够满足螺旋桨计算精度的需要，因此研究中采用基于诱导速度势的低阶迭代面元法.对于均匀来流中的单个螺旋桨问题，由于对称性，面元法求解时只需要求解一个桨叶^[21]，这大大缩短了螺旋桨的求解时间.对于吊舱推进器、对转桨以及混合式CRP等，各个构件的相互影响打破了这种对称性，即使是均匀来流, 各个桨叶上的奇点强度也是不相同的，因此对于两个构件的组合推进器一般采用迭代方法，即是先计算各个构件在对方位置处的诱导速度或者诱导速度势，取周向平均，这样就可以把各个构件拆开，分别求解，并且每一个构件可以认为是定常问题，通过相互之间的诱导速度或诱导速度势进行迭代求解，直至两者受力达到稳定，即求得结果.

由两个构件的迭代面元法的广泛应用，很容易想到也采用迭代方法求解混合式CRP，但为了简化迭代处理，需要把混合式CRP拆分为两个部分，这也是迭代面元法处理的关键所在.由于熊鹰等^[22-23]建立了吊舱推进器定常整体面元法，因此将吊舱推进器作为整体求解变得可能，于是本文的处理是将前桨当作一个整体，后桨与吊舱组成的吊舱推进器当作一个整体.前桨采用单个螺旋桨面元法进行计算，吊舱推进器采用熊鹰等$^{[22\hbox{-} 23]}$所建立的定常整体面元法进行计算.同样，螺旋桨和吊舱推进器之间的影响通过诱导速度来计算，对螺旋桨及吊舱推进器在对方位置处的诱导速度取周向平均，于是可以写出混合式CRP迭代面元法中前桨以及吊舱推进器的诱导速度势表达式^[21-22]

$ \phi _{i{\rm F}} = \sum\limits_{k_1 = 1}^{Z_1 } \sum\limits_{j = 1}^M \phi _j C_{ij}^{k_1 } - \sum\limits_{k_1 = 1}^{Z_1 } \sum\limits_{j = 1}^M [({\boldsymbol U} +{\boldsymbol V}_{\rm FP} ) \cdot {\boldsymbol n}]_j {\boldsymbol B}_{ij}^{k_1 } + \\ \qquad \sum\limits_{k_1 = 1}^{Z_1 } \sum\limits_{j = 1}^{M_W } \Delta \phi_j W_{ij}^{k_1 }, \ \ i = 1, 2, \cdots, M $

(1)

$ \phi _{i{\rm P}} = \sum\limits_{k_2 = 1}^{Z_2 } \sum\limits_{j = 1}^N \bar\phi _j C_{ij}^{k_2} + \sum\limits_{k_2 = 1}^{Z_2 } \sum\limits_{j= 1}^N [({\boldsymbol U} +{\boldsymbol V}_{\rm PF} ) \cdot {\boldsymbol n}]_j {\boldsymbol B}_{ij}^{k_2 } + \\ \;\;\;\;\;\;\sum\limits_{k_2 = 1}^{Z_2 } \sum\limits_{j = 1}^{N_W } \overline{\Delta \phi}_j W_{ij}^{k_2}+\dfrac 1{Z_2} \sum\limits^{Z_2}_{\delta=1} \sum\limits^Y_{l=1} \bar \phi_{l\delta} C_{il\delta}+\\ \qquad \dfrac 1{Z_2} \sum\limits^{Z_2}_{\delta=1} \sum\limits^Y_{l=1} [({\boldsymbol U} +{\boldsymbol V}_{\rm PF} ) \cdot {\boldsymbol n}]_{l\delta} {\boldsymbol B}_{il\delta}, \\ \qquad i = 1, 2, \cdots, N, \ N+1, \cdots, N+Y $

(2)

其中，由于物面不可穿透条件，式(1)和式(2)中的源汇影响项$\partial \phi / \partial n$用$({\boldsymbol U} +{\boldsymbol V}_{\rm FP} ) \cdot {\boldsymbol n}$或$({\boldsymbol U} +{\boldsymbol V}_{\rm PF} ) \cdot {\boldsymbol n}$代替. $\phi_{i{\rm F}} $表示前桨场点面元的速度势，$\overline \phi _{i{\rm P}} $表示吊舱推进器场点面元的速度势. $Z_1 $是前桨桨叶数，$Z_2 $是后桨桨叶数，$M$表示前桨上面元数目，$M_{\rm W} $表示前桨尾涡面元数目，$N$表示吊舱桨桨叶上面元数目，$N_{\rm W} $表示吊舱桨尾涡面上面元数目，$Y$表示吊舱上的面元数目. ${\boldsymbol U}$表示来流速度矢量，${\boldsymbol V}_{\rm FP} $表示吊舱推进器在前桨桨盘面处的周向平均诱导速度矢量，${\boldsymbol V}_{\rm PF} $表示前桨在吊舱螺旋桨盘面及吊舱位置周向平均诱导速度矢量. ${\boldsymbol n}$为面元的法向量，$C_{ij}^k $, $B_{ij}^k $, $W_{ij}^k $为影响系数，它们的定义如下

$ C_{ij}^k = \dfrac{1}{2\pi }\mathop{{\iint}\mkern-22.5mu \bigcirc}\limits_{S_{B_j } } \dfrac{\partial }{\partial n_{q_j } } \Big(\dfrac{1}{r(p_i, q_{jk} )} \Big ) {\rm{d}} S_{q_j } \\ B_{ij}^k = - \dfrac{1}{2\pi }\mathop{{\iint}\mkern-22.5mu \bigcirc}\limits_{S_{B_j } } \Big(\dfrac{1}{r(p_i, q_{jk} )} \Big) {\rm{d}} S_{q_j } \\ W_{ij}^k = \dfrac{1}{2\pi }\mathop{{\iint}\mkern-22.5mu \bigcirc}\limits_{S_{W_j } } \dfrac{\partial }{\partial n_{q_j } } \Big(\dfrac{1}{r(p_i, q_{jk} )} \Big){\rm{d}} S_{q_j } $

其中，$p_{i}$表示场点，$q_{jk}$表示$k$号桨叶上控制点，$n$表示控制点面元法向.另外，式(2)中下标$\delta $表示吊舱支柱与主桨叶错开的角度是吊舱螺旋桨相邻桨叶间隔角度的$\delta $倍.求解过程中式(1)与式(2)通过$V_{\rm FP} $和$V_{\rm PF} $迭代运行，直至前桨与吊舱推进器的受力达到稳定值即可认为求解完毕，本文在计算的过程中一般迭代5次左右即可收敛.

1.2 整体面元法

迭代面元法虽然能够解决混合式CRP的性能预报问题，但若结合有关数值方法进行混合式CRP的设计，由于这中间繁杂的迭代过程，仍旧会使得设计工作变得困难.为了探究更高效的混合式CRP面元法计算方法，研究中试图建立一种定常整体面元法，即计算中将混合式CRP视为整体，且前桨及后桨(吊舱螺旋桨)均只求解一个桨叶.在混合式CRP的分析中，仍旧按文献[22]中的处理将吊舱推进器当成一个整体，来考虑与前桨的相对位置.为了方便说明，以前桨4叶，后桨5叶为例进行混合式CRP定常整体计算公式的分析.结构位置示意图如图 1所示，图 1中实线代表后桨桨叶，虚线代表前桨桨叶，前桨桨叶用罗马数字标注桨叶编号，后桨桨叶用阿拉伯数字标注桨叶编号.吊舱以T型实线表示支柱中线.前、后桨桨叶的位置以螺旋桨参考线为准.

由于文献[22]所选取的求解位置中，吊舱的位置只与后桨各桨叶重合，于是考虑吊舱推进器整体与前桨的相对位置时，可以始终认为吊舱与后桨某一桨叶重合.如图 1所示，位置1处前桨Ⅰ号桨叶与后桨1号桨叶重合，而位置2是在位置1的基础上旋转了一个$\theta $角度使得前桨Ⅳ号桨叶与后桨5号桨叶重合.位置3则是前桨Ⅳ号桨叶与后桨1号桨叶重合.以$\theta $角度为间隔，前桨(或后桨)旋转一周的过程中，吊舱与各个后桨桨叶重合对前桨的影响效果是一致的，于是如图 1所示，在选定的各个位置下，吊舱始终与1号桨叶重合.由于文献[22]所建立的吊舱推进器整体面元法求解式(2)是一种平均状态，在此状态下各个桨叶的奇点强度保持对称性，因此在考虑以式(2)所建立的吊舱推进器系统与前桨的相互影响时，可以认为，在均匀来流下，前桨桨叶的受力是周期性变化的，变化的周期是旋转经过后桨相邻桨叶间隔角度的时间；后桨桨叶的受力变化周期是旋转经过前桨相邻桨叶间隔角度的时间.尽管前、后桨的桨叶受力变化周期不一致，但整个前、后螺旋桨的受力变化周期是一致的，均为螺旋桨(前桨或后桨)旋转角度$\theta $ (图 1中的标注)经历的时间.因此位置1、位置2、位置3处前桨、后桨受力均是相同的.在计算混合式CRP的敞水性能时，受力是平均意义上的.因此，在这里试图将受力的平均处理通过速度势来解决，从而选取前、后桨错开某些角度作为求解位置，并将各个求解位置的速度势平均处理.分别以图 1中前桨Ⅰ号桨叶与后桨1号桨叶作为研究对象(主桨叶)进行说明，从位置1开始，以$\theta $角度(图 1标示的)为间隔，将Ⅰ号桨叶或1号桨叶旋转一周所经历的所有位置作为求解位置.由这里所选取的螺旋桨桨叶数可以知道，$\theta$的大小为18$^\circ$，定义

$ S_0 = \dfrac{360^\circ}{\theta }, \ \ S_1 = \dfrac{360^\circ }{\theta \cdot Z_1 }, \ \ S_2 = \dfrac{360^\circ }{\theta \cdot Z_2 } $

(3)

于是可以写出在每一求解角度位置处，主桨叶上面元$i$上的速度势离散表达式(3)

$ \phi _{iS} = \sum\limits_{k_1 = 1}^{Z_1 } \sum\limits_{j = 1}^M \phi _{jS}^{k_1 } C_{ijS}^{k_1 } + \sum\limits_{k_1 = 1}^{Z_1 } \sum\limits_{j = 1}^M \Big(\dfrac{\partial \phi }{\partial n}\Big)_{jS}^{k_1 } B_{ijS}^{k_1 } +\\ \qquad \sum\limits_{k_1 = 1}^{Z_1 } \sum\limits_{j = 1}^{M_W } \Delta \phi _{jS}^{k_1 } W_{ijS}^{k_1 } + \sum\limits_{k_2 = 1}^{Z_2 } \sum\limits_{l = 1}^N \phi _{lS}^{k_2 } C_{ilS}^{k_2 } +\\ \qquad \sum\limits_{k_2 = 1}^{Z_2 } \sum\limits_{l = 1}^N \Big (\dfrac{\partial \phi }{\partial n}\Big)_{lS}^{k_2 } B_{ilS}^{k_2 } + \sum\limits_{k_2 = 1}^{Z_2 } \sum\limits_{l = 1}^{N_W } \Delta \phi _{lS}^{k_2 } W_{ilS}^{k_2 } +\\ \qquad \dfrac{1}{Z_2 }\sum\limits_{k_3 = 1}^{Z_2 } \sum\limits_{x = 1}^Y \phi _{xS}^{k_3 } C_{ixS}^{k_3 } + \dfrac{1}{Z_2 }\sum\limits_{k_3 = 1}^{Z_2 } \sum\limits_{x = 1}^Y \Big (\dfrac{\partial \phi }{\partial n}\Big)_{xS}^{k_3 } B_{ixS}^{k_3 }, \\ \qquad i = 1, 2, \cdots, M, M + 1, \cdots, M + N, \cdots, \\ \qquad \qquad M + N + Y ; \qquad S = 1, 2, \cdots, S_0 $

(4)

其中，$S_0 $如式(3)定义；下标$S$表示各求解的角度位置，根据角度间隔的选取，$S = 1$对应图 1中的位置1，即前、后主桨叶错开零度，$S = 2$对应图 1中的位置2，即前、后主桨叶错开一倍$\theta $，依此类推.影响系数$C_{ijS}^{k_1 } $表示在$S$角位置下，$k_1 $号桨叶上$j$号控制点对场点$i$的影响；$\phi _{jS}^{k_1 } $表示$S$角位置下，$k_1 $号桨叶上$j$号控制点的速度势；其余各影响系数项和奇点强度项的上、下标与此类似.但需要注意的是式(4)中的吊舱影响项的$k_3 $与式(2)中的$\delta $意义一致.将$S$取不同值时的式(4)叠加并取平均

$ \dfrac{1}{S_0 }\sum\limits_{S = 1}^{S_0 } \phi _{iS} = \dfrac{1}{S_0 }\left\{ {\left[{\sum\limits_{S = 1}^{S_0 } \sum\limits_{k_1 = 1}^{Z_1 } \sum\limits_{j = 1}^M \phi _{jS}^{k_1 } C_{ijS}^{k_1 } + }\right.}\right. \\ \quad \left.{\sum\limits_{S = 1}^{S_0 } \sum\limits_{k_1 = 1}^{Z_1 } \sum\limits_{j = 1}^M \Big (\dfrac{\partial \phi }{\partial n}\Big)_{jS}^{k_1 } B_{ijS}^{k_1 } + \sum\limits_{S = 1}^{S_0 } \sum\limits_{k_1 = 1}^{Z_1 } \sum\limits_{j = 1}^{M_W } \Delta \phi _{jS}^{k_1 } W_{ijS}^{k_1 } } \right]_1 +\\ \quad \left[{\sum\limits_{S = 1}^{S_0 } \sum\limits_{k_2 = 1}^{Z_2 } \sum\limits_{l = 1}^N \phi _{lS}^{k_2 } C_{ilS}^{k_2 } + \sum\limits_{S = 1}^{S_0 } \sum\limits_{k_2 = 1}^{Z_2 } \sum\limits_{l = 1}^N \Big(\dfrac{\partial \phi }{\partial n}\Big)_{lS}^{k_2 } B_{ilS}^{k_2 } + }\right. \\ \quad \left.{\sum\limits_{S = 1}^{S_0 } \sum\limits_{k_2 = 1}^{Z_2 } \sum\limits_{l = 1}^{N_W } \Delta \phi _{lS}^{k_2 } W_{ilS}^{k_2 } } \right]_2 + \left[{\sum\limits_{S = 1}^{S_0 } \dfrac{1}{Z_2 }\sum\limits_{k_3 = 1}^{Z_2 } \sum\limits_{x = 1}^Y \phi _{xS}^{k_3 } C_{ixS}^{k_3 } + }\right. \\ \quad \left.{\left.{\sum\limits_{S = 1}^{S_0 } \dfrac{1}{Z_2 }\sum\limits_{k_3 = 1}^{Z_2 } \sum\limits_{x = 1}^Y \Big (\dfrac{\partial \phi }{\partial n}\Big)_{xS}^{k_3 } B_{ixS}^{k_3 } } \right]_3 } \right\}, \\ \quad i = 1, 2, \cdots, M, M + 1, \cdots , M + N, \cdots, M + N + Y ; \\ \quad S = 1, 2, \cdots, S_0 $

(5)

为了能够运用定常面元法求解式(5)，需要对式(5)进行相关的数值处理和分析.式(5)中标号为1的项是前桨的影响项，标号为2的项为后桨影响项，标号为3的项为吊舱影响项.显然，在$S$从1到$S_0 $的变化过程中，标号1, 2和3中各自项的变化特点是一致的，因此将式(5)拆分为前桨以及吊舱推进器的表达形式，且选取标号1, 2和3中各一项进行分析.

$ \dfrac{1}{S_0 }\sum\limits_{S = 1}^{S_0 } \phi _{iS} = \dfrac{1}{S_0 }\left( {\sum\limits_{S = 1}^{S_0 } \sum\limits_{k_1 = 1}^{Z_1 } \sum\limits_{j = 1}^M \phi _{jS}^{k_1 } C_{ijS}^{k_1 } +}\right.\\ \qquad \sum\limits_{S = 1}^{S_0 } \sum\limits_{k_2 = 1}^{Z_2 } \sum\limits_{l = 1}^N \phi _{lS}^{k_2 } C_{ilS}^{k_2 } + \left. {\dfrac{1}{Z_2 }\sum\limits_{S = 1}^{S_0 } \sum\limits_{k_3 = 1}^{Z_2 } \sum\limits_{x = 1}^Y \phi _{xS}^{k_3 } C_{ixS}^{k_3 } + \cdots } \right), \\ \qquad i = 1, 2, \cdots, M; \ S = 1, 2, \cdots, S_0 $

(6)

$ \dfrac{1}{S_0 }\sum\limits_{S = 1}^{S_0 } {\phi _{iS} } = \dfrac{1}{S_0 }\left( {\sum\limits_{S = 1}^{S_0 } {\sum\limits_{k_1 = 1}^{Z_1 } {\sum\limits_{j = 1}^M {\phi _{jS}^{k_1 } } } C_{ijS}^{k_1 } } + } \right.\\ \qquad \sum\limits_{S = 1}^{S_0 } {\sum\limits_{k_2 = 1}^{Z_2 } {\sum\limits_{l = 1}^N {\phi _{lS}^{k_2 } } } C_{ilS}^{k_2 } } + \\ \qquad {\left. {\dfrac{1}{Z_2 }\sum\limits_{S = 1}^{S_0 } {\sum\limits_{k_3 = 1}^{Z_2 } {\sum\limits_{x = 1}^Y {\phi _{xS}^{k_3 } } } C_{ixS}^{k_3 } } + \cdots } \right) } \\ \qquad i = M + 1, \cdots, M + N, \cdots, M + N + Y; \\ \qquad S = 1, 2, \cdots, S_0 $

(7)

其中，省略号代表了式(5)中未列出的项.可以知道，式(6)和式(7)均选取了式(5)中偶极子影响项.首先对式(6)进行分析，式(6)中右端第一项是前桨上控制点对前桨自身场点的影响，因此$S$取不同值，$k_1 $取相同值时，影响系数$C_{ijS}^{k_1 } $是相同的，于是可以将其合并

$ \dfrac{1}{S_0 }\sum\limits_{S = 1}^{S_0 } {\sum\limits_{k_1 = 1}^{Z_1 } {\sum\limits_{j = 1}^M {\phi _{jS}^{k_1 } } } C_{ijS}^{k_1 } } = \sum\limits_{k_1 = 1}^{Z_1 } {\sum\limits_{j = 1}^M {\Big(\dfrac{1}{S_0 }} } \sum\limits_{S = 1}^{S_0 } {\phi _{jS}^{k_1 }\Big )C_{ijS}^{k_1 } } $

(8)

根据前述吊舱推进器作为整体的处理，参照图 1可以知道，式(8)中$S$取不同值，$\phi _{jS}^{k_1 } $是不同的，但$\phi _{jS}^{k_1 } $是周期性变化的，且不同桨叶上的奇点强度在不同$S$值下存在对应关系.以图 1中位置1和位置2为例，$S = 2$时，Ⅰ号桨叶上面元奇点强度与$S = 1$时，Ⅱ号桨叶面元奇点强度是相同的，同样，其他桨叶也有类似的对应关系，且$S$达到位置3时，又重复上述的等效关系.显然，对于前桨，这种等效关系在$S$从1到$S_0 $的取值过程会重复$S_1 $次.于是式(8)中的右端速度势项可以写成

$ \dfrac{1}{S_0 }\sum\limits_{S = 1}^{S_0 } \phi _{jS}^{k_1 } = \dfrac{1}{S_1 S_2 }S_1 (\phi _j^1 + \phi _j^2 + \phi _j^3 + \phi _j^4 ) = \\ \qquad \dfrac{1}{Z_1 }\sum\limits_{k_1 = 1}^{Z_1 } {\phi _j^{k_1 } } $

(9)

显然式(8)中，$k_1 $取不同值时，右端速度势项均可写成式(9)形式，故对于不同$k_1 $值，这一项是相等的.令$\overline \phi _j = \dfrac{1}{Z_1 }\sum\limits_{k_1 = 1}^{Z_1 } {\phi _j^{k_1 } } $，于是式(9)可以写为

$ \dfrac{1}{S_0 }\sum\limits_{S = 1}^{S_0 } {\sum\limits_{k_1 = 1}^{Z_1 } {\sum\limits_{j = 1}^M {\varphi _{jS}^{k_1 } } } C_{ijS}^{k_1 } } = \sum\limits_{k_1 = 1}^{Z_1 } {\sum\limits_{j = 1}^M {\overline \varphi _j C_{ij}^{k_1 } } } $

(10)

于是得到式(6)右端第1项的最终表达式(10).下面分析式(6)右端第2项，第2项是后桨面元对前桨面元的影响项，因此$S$取不同值时，即使$k_2 $取相同值，$C_{ilS}^{k_2 } $也是不断变化的，且$C_{ilS}^{k_2 } $所对应的$\phi _{lS}^{k_2 } $也是不断变化的.因此式(6)右端第2项的处理不能简单按第1项处理办法.但仔细分析式(6)右端第2项可以知道，$C_{ilS}^{k_2 } $和$\phi _{lS}^{k_2 } $在S取不同值时仍旧存在一定的对应关系.以图 1中位置1, 2, 3为例予以说明：一方面，对照位置1和位置3可以知道，后桨1桨叶面元上的奇点强度$\phi _{lS}^{k_2 } $在位置1和位置3处是相等的，但对前桨Ⅰ桨叶上面元的影响系数是不同的，因此，可以将奇点强度相同的项的影响系数进行合并.另一方面，对照位置2和位置3可以知道，位置2中后桨1号桨叶与位置3中后桨2号桨叶对前桨Ⅰ号桨叶上面元的影响系数完全一致，且两者的奇点强度也相同.显然，$S$从1到$S_0 $的变化过程中，后桨各桨叶都重复这两方面的特征.于是式(6)右端第2项可以写为

$ \dfrac{1}{S_0 }\sum\limits_{S = 1}^{S_0 } {\sum\limits_{k_2 = 1}^{Z_2 } {\sum\limits_{l = 1}^N {\phi _{lS}^{k_2 } } } C_{ilS}^{k_2 } } = \dfrac{1}{S_1 }S_1 \sum\limits_{k_2 = 1}^{Z_2 } {\sum\limits_{l = 1}^N {\phi _{l1}^{k_2 } E_{il}^{k_2 } } } =\\ \qquad \sum\limits_{k_2 = 1}^{Z_2 } {\sum\limits_{l = 1}^N {\phi _{l1}^{k_2 } E_{il}^{k_2 } } } $

(11)

其中，$k_2 $从1到$Z_2 $的变化，$E_{il}^{k_2 } $满足

$ \left. \begin{array}{l} E_{il}^1 = (C_{il1}^1 + C_{il6}^1 + C_{il11}^1 + C_{il16}^1)/{S_2}\\ E_{il}^2 = (C_{il2}^1 + C_{il7}^1 + C_{il12}^1 + C_{il17}^1)/{S_2}\\ E_{il}^3 = (C_{il3}^1 + C_{il8}^1 + C_{il13}^1 + C_{il18}^1)/{S_2}\\ E_{il}^4 = (C_{il4}^1 + C_{il9}^1 + C_{il14}^1 + C_{il19}^1)/{S_2}\\ E_{il}^5 = (C_{il5}^1 + C_{il10}^1 + C_{il15}^1 + C_{il20}^1)/{S_2} \end{array} \right\} $

(12)

为了更形象地说明问题，式(12)中影响系数的下标$S$用具体的数值代替了，显然，这里是以本文选取的后桨桨叶为5叶进行下标$S$的标注.以$C_{il1}^1 $为例，表示在角位置$S = 1$处，后桨1桨叶面元上控制点$l$对前桨主桨叶上场点$i$的影响系数.另外，注意式(12)中后桨其他桨叶对场点的影响系数都用后桨1号桨叶代替了，这是由于前述讨论了不同$S$位置处，桨叶上面元影响系数的等效关系.对照图 1，观察式(12)影响系数项可以知道，不同$k_2 $，对应的是将$E_{il}^{k_2 } $进行了旋转改变，这一旋转改变使得每一个影响系数$C_{ilS}^1 $的值是变化的.但从数值求解的角度，可以认为$E_{il}^{k_2 } $是近似相等的，这在后文的误差分析中会有详细说明，于是可将式(11)右端项改为

$ \sum\limits_{k_2 = 1}^{Z_2 } \sum\limits_{l = 1}^N \phi _{l1}^{k_2 } E_{il}^{k_2 } = \sum\limits_{l = 1}^N \Big (\sum\limits_{k_2 = 1}^{Z_2 } \phi _{l1}^{k_2 } \Big)E_{il}^1 =\\ \qquad \sum\limits_{k_2 = 1}^{Z_2 } \sum\limits_{l = 1}^N \Big ( \sum\limits_{k_2 = 1}^{Z_2 } \phi _{l1}^{k_2 } \Big / {Z_2 } \Big )E_{il}^1 = \sum\limits_{k_2 = 1}^{Z_2 } \sum\limits_{l = 1}^N \Big( \sum\limits_{k_2 = 1}^{Z_2 } \phi _{l1}^{k_2 } \Big / {Z_2 } \Big )E_{il}^{k_2 } $

(13)

显然，对于不同的$k_2 $，$\sum\limits_{k_2 = 1}^{Z_2 } \phi _{l1}^{k_2 } \Big / {Z_2 } $是相同的，令$\sum\limits_{k_2 = 1}^{Z_2 } \phi _{l1}^{k_2 } / {Z_2 } = \overline {\phi _l } $，结合式(11)和式(13)，得到式(6)右端第2项的最终表达式为

$ \dfrac{1}{S_0 }\sum\limits_{S = 1}^{S_0 } {\sum\limits_{k_2 = 1}^{Z_2 } {\sum\limits_{l = 1}^N {\varphi _{lS}^{k_2 } } } C_{ilS}^{k_2 } } = \sum\limits_{k_2 = 1}^{Z_2 } {\sum\limits_{l = 1}^N {\overline \varphi _l } } E_{il}^{k_2 } $

(14)

对于式(6)中的右端第3项是吊舱对前桨面元的影响，由于分析中将吊舱与后桨当作一个整体，因此它的影响特性与后桨影响项保持一致.于是式(6)中右端第3项可以写为

$ \dfrac{1}{Z_2 }\sum\limits_{S = 1}^{S_0 } \sum\limits_{k_3 = 1}^{Z_2 } \sum\limits_{x = 1}^Y \phi_{xS}^{k_3 } C_{ixS}^{k_3 } = \\ \qquad \sum\limits_{x = 1}^Y \Big(\dfrac{1}{S_1 }\sum\limits_{S = 1}^5 \phi _{xS} \Big ) \sum\limits_{k_3 = 1}^{Z_2 } (C_{ix1}^{k_3 } +\\ \qquad C_{ix6}^{k_3 } + C_{ix11}^{k_3 } + C_{ix16}^{k_3 } ) / Z_2 S_2 = \sum\limits_{x = 1}^Y \overline \phi_x E_{ix} $

(15)

其中，$\overline \phi _x = \dfrac{1}{S_1 }\sum\limits_{S = 1}^5 \phi_{xS} $，$E_{ix} $满足

$ E_{ix} = \dfrac{1}{S_0 }\sum\limits_{S = 1}^{S_0 } {C_{ixS} } $

(16)

由于$k_3 $取不同值时，式(15)中$C_{ixS}^{k_3 } $项与$C_{ixS} $存在对应关系，于是式(16)中的$C_{ixS} $的上标$k_3 $可以去掉.由式(10)、式(14)和式(15)得到式(6)的最终表达式

$ \overline \phi _i = \sum\limits_{k_1 = 1}^{Z_1 } \sum\limits_{j = 1}^M \overline \phi _j C_{ij}^{k_1 } + \sum\limits_{k_2 = 1}^{Z_2 } \sum\limits_{l = 1}^N \overline \phi _l E_{il}^{k_2 } + \\ \qquad \sum\limits_{x = 1}^Y \overline \phi _x E_{ix} + \cdots, \qquad i = 1, 2, \cdots, M $

(17)

下面对式(7)进行分析，对于式(7)的右端第1项是前桨桨叶面元对后桨桨叶面元的影响，因此对它的分析与式(6)中右端第2项的分析方法类似，于是可以写出式(7)中右端第1项的表达式

$ \sum\limits_{S = 1}^{S_0 } {\sum\limits_{k_1 = 1}^{Z_1 } {\sum\limits_{j = 1}^M {\varphi _{jS}^{k_1 } } } C_{ijS}^{k_1 } } = \sum\limits_{k_1 = 1}^{Z_1 } {\sum\limits_{j = 1}^M {\overline \varphi _j E_{ij}^{k_1 } } } $

(18)

其中，$E_{ij}^{k_1 } $满足

$ \left. \begin{array}{l} E_{ij}^1 = (C_{ij1}^1 + C_{ij5}^1 + C_{ij9}^1 + C_{ij13}^1 + C_{ij17}^1)/{S_1}\\ E_{ij}^2 = (C_{ij2}^1 + C_{ij6}^1 + C_{ij10}^1 + C_{ij14}^1 + C_{ij18}^1)/{S_1}\\ E_{ij}^3 = (C_{ij3}^1 + C_{ij7}^1 + C_{ij11}^1 + C_{ij15}^1 + C_{ij19}^1)/{S_1}\\ E_{ij}^4 = (C_{ij4}^1 + C_{ij8}^1 + C_{ij12}^1 + C_{ij16}^1 + C_{ij20}^1)/{S_1} \end{array} \right\} $

(19)

由于式(7)的右端第2、3项是吊舱推进器上面元对自身面元的影响项，而式(7)是在吊舱推进器整体面元法的基础上建立的，因此这两项与式(6)中右端第1项前桨对自身面元影响项的分析类似，于是式(7)中右端第2、3项可以写为

$ \sum\limits_{S = 1}^{S_0 } {\sum\limits_{k_2 = 1}^{Z_2 } {\sum\limits_{l = 1}^N {\phi _{lS}^{k_2 } } } C_{ilS}^{k_2 } } + \dfrac{1}{Z_2 }\sum\limits_{S = 1}^{S_0 } {\sum\limits_{k_3 = 1}^{Z_2 } {\sum\limits_{x = 1}^Y {\phi _{xS}^{k_3 } } } C_{ixS}^{k_3 } } =\\ \qquad \sum\limits_{k_2 = 1}^{Z_2 } {\sum\limits_{l = 1}^N {\overline \phi _l } } C_{il}^{k_2 } + \sum\limits_{x = 1}^Y {\overline \phi _x } E_{ix} $

(20)

其中，$E_{ix} $满足

$ E_{ix} = \dfrac{1}{Z_2 }\sum\limits_{k_3 = 1}^{Z_2 } {C_{ix}^{k_3 } } $

(21)

显然，当式(21)中的面元$i$表示吊舱上场点时，根据式(4)中有关$k_3 $意义的描述可以知道，不同$k_3 $下，$C_{ix}^{k_3 } $其实是相同的.

根据式(18)和式(20)得到式(7)的最终表达式

$ \overline \phi _i = \sum\limits_{k_1 = 1}^{Z_1 } {\sum\limits_{j = 1}^M {\overline \phi _j E_{ij}^{k_1 } } } + \sum\limits_{k_2 = 1}^{Z_2 } {\sum\limits_{l = 1}^N {\overline \phi _l } } C_{il}^{k_2 } + \sum\limits_{x = 1}^Y {\overline \phi _x E_{ix} } + \cdots, \\ i = M + 1, M + 2, \cdots, M + N + 1, \cdots, M + N + Y $

(22)

为了统一，将式(17)和式(22)中的影响系数$C$均用$E$来表达，可以写出式(5)的最终表达式

$ \overline {\phi _i } = \sum\limits_{k_1 = 1}^{Z_1 } {\sum\limits_{j = 1}^M {\overline \phi _j E_{ij}^{k_1 } } } + \sum\limits_{k_1 = 1}^{Z_1 } {\sum\limits_{j = 1}^M {\Big (\overline {\dfrac{\partial \phi }{\partial n}} \Big )_j F_{ij}^{k_1 } } } + \\ \qquad \sum\limits_{k_1 = 1}^{Z_1 } \sum\limits_{j = 1}^{M_W } \overline {\Delta \phi } _j G_{ij}^{k_1 } + \sum\limits_{k_2 = 1}^{Z_2 } \sum\limits_{l = 1}^N \overline \phi _l E_{il}^{k_2 } + \\ \qquad \sum\limits_{k_2 = 1}^{Z_2 } \sum\limits_{l = 1}^N \Big(\overline {\dfrac{\partial \phi }{\partial n}} \Big)_l F_{il}^{k_2 } + \sum\limits_{k_2 = 1}^{Z_2 } \sum\limits_{l = 1}^{N_W } \overline {\Delta \phi } _l G_{il}^{k_2 } +\\ \qquad \sum\limits_{x = 1}^Y \overline \phi _x E_{ix} + \Big(\overline {\dfrac{\partial \phi }{\partial n}} \Big)_x F_{ix}, \\ \qquad i = 1, 2, \cdots, M, M + 1, \cdots, M + N, \cdots, \\ \qquad \qquad M + N + Y $

(23)

其中，影响系数$E$满足式(24)，影响系数$F$和$G$的表达形式与$E$类似，只需将$E$中的影响系数$C$项换成$B$项或$W$项即可.

$ E = \left\{ \begin{array}{l} \begin{array}{*{20}{l}} {{E_{ij}} = {C_{ij}}, }&{i = {P_{{\rm{FP}}}}}\\ {{E_{il}} = {\rm{Eq}}.({\rm{12}}), }&{i = {P_{{\rm{FP}}}}}\\ {{E_{ix}} = {\rm{Eq}}.({\rm{16}}), }&{i = {P_{{\rm{FP}}}}}\\ {{E_{ij}} = {\rm{Eq}}.({\rm{19}}), }&{i = {P_{{\rm{POD}}}}}\\ {{C_{il}}, }&{i = {P_{{\rm{POD}}}}} \end{array}\\ \begin{array}{*{20}{l}} {{E_{ix}} = {\rm{Eq}}.({\rm{21}}), }&{i = {P_{{\rm{POD}}}}} \end{array} \end{array} \right. $

(24)

其中，$E_{il} = {\rm Eq.(12)}$表示$E_{il} $按式(12)定义，其他项表达与此类似，$i = P_{\rm FP} $表示$i$是前桨上面元，$i = P_{\rm POD} $表示$i$是吊舱上面元.下标$j$表示控制点在前桨桨叶上，$l$表示控制点在后桨桨叶上，$x$表示控制点在吊舱上.式(23)即是整体求解混合式CRP的面元法计算公式，观察式(23)，可以知道，不同桨叶的$\overline \phi $, $ \overline {\dfrac{\partial \phi }{\partial n}} $以及$\overline {\Delta \phi } $是相等的，因此，求解中前后桨只需各求解一个桨叶，这将大大减小整体计算混合式CRP的计算网格，同时简化了数值处理.

需要注意的是，式(23)是在吊舱推进器定常整体面元法计算式(2)的基础上得到的，因此对于吊舱桨与吊舱的求解间隔是后桨相邻桨叶的间隔(对于后桨为5叶时即是72$^\circ$).这样的角度间隔表示忽略角度间隔之间螺旋桨受力的脉动值，对于不同螺旋桨类型的适应性有待于进一步研究，为了使得求解的结果更准确，可以在混合式CRP的求解中，将后桨与吊舱推进器的求解间隔也选为图 1标示的$\theta $(即18$^\circ$)，显然这样的选取对求解的受力信息更全面.于是，将式(23)推广到后桨与吊舱求解角度间隔也是$\theta$，分析得到混合式CRP的整体面元法计算公式仍旧与式(23)保持一致，而影响系数的表达会有所改变，见式(25)

$ E = \left\{ {\begin{array}{*{20}{l}} {{E_{ij}} = {C_{ij}}, }&{i = {P_{{\rm{FP}}}}}\\ {{E_{il}} = {\rm{Eq}}.({\rm{12}}), }&{i = {P_{{\rm{FP}}}}\;({\rm{or}}\;{P_{{\rm{POD}}}})}\\ {{E_{ix}} = {\rm{Eq}}.({\rm{16}}), }&{i = {P_{{\rm{FP}}}}\;({\rm{or}}\;{P_{{\rm{AP}}}})}\\ {{E_{ij}} = {\rm{Eq}}.({\rm{19}}), }&{i = {P_{{\rm{AP}}}}\;({\rm{or}}\;{P_{{\rm{POD}}}})}\\ {{E_{il}} = {C_{il}}, }&{i = {P_{{\rm{AP}}}}}\\ {{E_{ix}} = {C_{ix}}, }&{i = {P_{{\rm{POD}}}}} \end{array}} \right. $

(25)

其中，$i = P_{\rm FP} $表示$i$是前桨上面元，$i = P_{\rm POD} $表示$i$是吊舱上面元，$i = P_{\rm AP} $表示$i$是后桨上面元.式(23)及其影响系数表达式(25)即是本文建立的混合式CRP定常整体面元法计算公式.另外需要说明的是，本文在推导混合式CRP的定常整体面元法公式时，以前桨4叶、后桨5叶为例，目的是为了使推导过程更加具体，对于其他形式的桨叶数目，本文的方法同样适用，只是在影响系数$E$, $F$和$G$的表达有所区别.

与单个螺旋桨的面元法求解相似，整体面元法计算式(23)中的$\overline {\Delta \phi } $也需要通过库塔条件确定.在这里仍然采用等压库塔条件^[24]，但由于前后桨是整体求解，因此需要将前后桨随边面元的等压条件联合成同一方程组，见式(26)所示^[25-26]

$ \left. \begin{array}{l} \Delta {p_m} = {p_{m{\rm{U}}}}-{p_{m{\rm{L}}}}\\ m = 1, 2, \cdots, {m_{\rm{F}}}, \\ \qquad {m_F} + 1, \cdots, {m_{\rm{F}}} + {m_A} \end{array} \right\} $

(26)

其中，$m_{\rm F} $为前桨径向尾涡数目，$m_{\rm A }$表示后桨径向尾涡数目.联立式(23)与式(26)即组成了整体面元法的求解方程组.但$\Delta \phi $与$\Delta p$之间是非线性关系，直接求解会有一定困难，于是建立雅可比矩阵，采用牛顿迭代法得到$\Delta \phi $与$\Delta p$的关系如式(27)所示^[25]

$ \Delta {\boldsymbol\phi}^{(k + 1)} = \Delta {\boldsymbol\phi}^{(k)}-{\boldsymbol J}^{ - 1}\Delta {\boldsymbol p}^{(k)} $

(27)

其中，上标$k + 1$表示第$k + 1$次迭代结果，$k$表示第$k$次迭代结果. ${\boldsymbol J}^{ - 1}$是雅可比矩阵的逆.雅可比矩阵的元素可由式(28)确定

$ \left. \begin{array}{l} {J_{ij}} = \frac{{\partial {{(\Delta p)}_i}}}{{\partial {{(\Delta \phi )}_j}}}\\ i({\rm{or}}\;j) = 1, 2, \cdots, {m_F}, {m_F} + 1, \cdots, {m_F} + {m_A} \end{array} \right\} $

(28)

基于式(23)、式(26)及式(27)即建立了整体面元法封闭求解方程组.

另外，尾涡面作为面元法求解的理论假设，对计算结果有较大影响.在求解之前，需要对混合式CRP的尾涡面进行说明.本文在确定混合式CRP的尾涡面时，主要依据已有对吊舱推进器和对转桨的尾流场研究结果以及混合式CRP的水动力性能研究结果.文献[27-28]的研究结果表明，对于吊舱推进器，吊舱支柱对螺旋桨的尾流既有阻塞的效果，也有排挤加速的效果.文献[29]的研究结果表明，与单独存在前桨的情况相比，对转桨前桨的尾流在流经后桨之前变化不大，流经后桨之后将会有较大幅度的加速；同时，后桨的尾流会得到前桨的加速，前桨造成的旋转尾流会与后桨旋转流进行一定程度的抵消.而已有对混合式CRP敞水性能的研究结果表明[6, 30\hbox{-}31]，混合式CRP的前桨与单独前桨的推力、扭矩相差不大.结合以上研究结果，本文主要对混合式CRP的尾涡螺距角作出式(29)的修改

$ \left. {\begin{array}{*{20}{l}} {{\beta _{{\rm{HF}}}} = \left\{ {\begin{array}{*{20}{l}} {{\omega _{{\rm{F1}}}}{\beta _{{\rm{SF}}}}, }&{{x_{{\rm{wake}}}} \le {x_{{\rm{A-plane}}}}}\\ {{\omega _{{\rm{F2}}}}{\beta _{{\rm{SF}}}}, }&{{x_{{\rm{wake}}}} > {x_{{\rm{A-plane}}}}} \end{array}} \right.}\\ {{\beta _{{\rm{HA}}}} = {\omega _{\rm{A}}}{\beta _{{\rm{SA}}}}} \end{array}{\rm{ }}} \right\} $

(29)

其中，$x_{\rm wake} $表示前桨尾涡面轴向坐标，$x_{\rm A - plane} $表示后桨桨盘面轴向坐标. $\beta _{\rm HF} $表示混合式CRP中前桨尾涡螺距角，$\beta _{\rm SF} $表示单独前桨相应轴向位置处的尾涡螺距角；$\beta _{\rm CA} $表示混合式CRP中后桨尾涡螺距角，$\beta _{\rm SA} $表示单独后桨相应轴向位置处的尾涡螺距角；$\omega _{\rm F1} $, $\omega _{\rm F2} $和$\omega _{\rm A } $为3个修正系数，它随混合式CRP的工作工况而变化，且具有一定的经验性.根据式(29)，前桨尾涡面分为两部分：后桨盘面之前和后桨盘面之后，这主要是由于后桨对前桨的抽吸影响较大.后桨的尾涡面修正主要考虑前桨旋转尾流与吊舱的阻塞作用.另外，前后桨的尾涡面在支柱内部时，认为其对场点的影响系数为零，从前、后桨叶根处泻出的尾涡面在经过吊舱段时，紧贴吊舱表面.

2 数值结果

前一小节分别给出了混合式CRP的迭代面元法和整体面元法求解公式，并对混合式CRP的面元法求解尾涡面进行了修正，随后按照求解公式编译完成迭代面元法和整体面元法的计算程序.为了对比求解精度与求解效果，运用两种面元法对文献[10]中的混合式CRP进行计算并与文献[10]中的试验值进行对比.混合式CRP的主要参数如表 1和表 2所示，吊舱的示意图如图 2所示.前、后桨间距为0.454 5倍前桨直径.

定义$K_{\rm TF} $，$K_{\rm TA} $，$K_{\rm TU} $，$K_{\rm QF} $和$K_{\rm QA} $如式(30)所示. $K_{\rm TF} $是前桨推力系数；$K_{\rm TA} $是后桨推力系数；$K_{\rm TU} $是吊舱推进器的整体推力系数; $K_{\rm QF} $是前桨扭矩系数；$K_{\rm QA} $是后桨扭矩系数.

$ \left. \begin{array}{l} J = \frac{V}{{{n_{\rm{F}}}{D_{\rm{F}}}}}\\ {K_{{\rm{TF}}}} = \frac{{{T_{\rm{F}}}}}{{\rho n_{\rm{F}}^2D_{\rm{F}}^4}}, {K_{{\rm{QF}}}} = \frac{{{Q_{\rm{F}}}}}{{\rho n_{\rm{F}}^2D_{\rm{F}}^5}}\\ {K_{{\rm{TA}}}} = \frac{{{T_{\rm{A}}}}}{{\rho n_{\rm{A}}^2D_{\rm{A}}^4}}, {K_{{\rm{QA}}}} = \frac{{{Q_{\rm{A}}}}}{{\rho n_{\rm{A}}^2D_{\rm{A}}^5}}\\ {K_{{\rm{TU}}}} = \frac{{{T_{\rm{U}}}}}{{\rho n_{\rm{A}}^2D_{\rm{A}}^4}} = \frac{{{T_{\rm{A}}} + {T_{{\rm{pod}}}}}}{{\rho n_{\rm{A}}^2D_{\rm{A}}^4}} \end{array} \right\} $

(30)

其中，$\rho $是流体密度；$T_{\rm F}$是前桨轴向推力；$Q_{\rm F}$是前桨扭矩；$T_{\rm A}$后桨轴向推力；$T_{pod} $吊舱轴向推力；$Q_{\rm A}$后桨扭矩；$V$是来流速度；$n_{\rm F}$是前桨转速；$D_{\rm F}$是前桨直径；$n_A $是后桨转速；$D_{\rm A}$是后桨直径.计算结果及与试验值的对比如图 3和图 4所示.计算中前桨及桨毂的网格数目为928，后桨及桨毂的网格数目为694，吊舱上网格数目为1 065，整个混合式CRP总的网格数目为2 687.

图 3和图 4中，iterative method表示迭代面元法计算结果，integral method表示整体面元法计算结果，experimental data表示试验值.通过图 3和图 4看出，整体面元法和迭代面元法的数值结果与试验值相比均有较好的一致性，在设计点$J = 0.781$处，相对误差均控制在5%以内.对于前桨的推力、扭矩系数，整体定常面元法的计算精度要优于迭代面元法，计算精度基本控制在2%以内.对于后桨的推力、扭矩系数，迭代面元法与整体面元法的计算精度相当.因此从工程应用的角度，迭代面元法与整体面元法均能满足计算精度的需要.研究中还对两种面元法的计算耗时进行了对比，两种方法的计算网格完全相同，计算所用的计算机也完全相同(4核Intel Core i5-4440 CPU @ 3.10GHz)，计算耗时对比如表 3所示.

根据表 3的计算耗时对比以及图 3和图 4计算结果的对比可以知道，在保证相同计算精度的条件下，整体面元法的计算时间要比迭代面元法的计算时间少.一般而言，推进器的设计会结合数值算法进行优化，这中间的计算过程需要进行大量的迭代处理，而本文建立的混合式CRP的整体计算面元法省去了面元法中的迭代，这对于混合式CRP的性能分析和设计具有较大意义.同时，根据混合式CRP整体面元法计算公式的推导可以知道，本文所建立的整体面元法可推广到其他多构件组成的组合推进器，但在应用的过程中需要对相关影响系数进行修正.

3 误差分析

基于势流理论的面元法在求解升力体问题时，将下泄的非势流区当作无厚度的势流区边界，这是面元法的理论缺陷，然而大量的螺旋桨计算结果表明^[18-21]，这一理论缺陷，可以通过建立合适的尾涡模型进行消除，从而使得数值结果与试验值有较好的一致性.本文所建立的定常整体面元法计算结果误差，同样受到尾涡模型的影响.另外，在定常整体面元法的推导过程中还作了以下近似处理，这里对这些近似处理的合理性作出说明.

第1项近似处理是：式(13)、式(18)的处理中，认为对于不同的$k_1 $(或$k_2 $), 影响系数项是相等的.以式(12)为例，认为$E_{il}^1 \approx E_{il}^2 \approx E_{il}^3 \approx E_{il}^4 \approx E_{il}^5 $.这里并不从理论公式上推导其误差，只给出数值结果的对比，以前桨对后桨影响系数$E_{ij}^{k_1 } $和后桨对前桨影响系数$E_{il}^{k_2 } $为例，通过数值结果的对比进行说明. $E_{il}^{k_2 } $值的对比如图 5和图 6所示，以$l$取为100和500，$i$取为50, 100, 150, 200, 250, 300, 350, 400, 450, 500, 550和600为例. $E_{ij}^{k_1 } $值的对比如图 7和图 8所示，以$j$取为200和600，$i$取为1 000, 1 050, 1 100, 1 150, 1 200, 1 250, 1 300, 1 350, 1 400, 1 450, 1 500和1 550为例(根据本文计算中网格的划分可以知道，$i < 929$表示前桨上的面元，$i>928$且$ < 1 622$表示后桨上的面元).

从图 5 $\sim $图 8可以看出，影响系数$E_{il}^{k_2 } $及$E_{ij}^{k_1 } $数值的量级为10$^{ - 4}$，而不同$k_2 $值，$E_{il}^{k_2 } $的差别大部分在1%以内，个别较大点的差值会达到7%，$E_{ij}^{k_1 } $的计算结果也具有类似规律.研究中，对其他面元的影响系数以及前后桨与吊舱之间的影响系数也进行了计算，得到数值结果的差别与本文中展示的数值结果差别类似.因此，可以认为式(13)及式(18)的处理所造成的数值误差是合理的.

第2项近似处理是：选取的求解角度间隔问题.选取一定的角度间隔表示忽略了角度间隔之间的螺旋桨受力的脉动值.根据文献[11-12, 23]的结果可以知道，就本文所选取的角度间隔，螺旋桨或吊舱的受力脉动值在1%以内，因此从工程计算的角度考虑，这一处理也是合理的.另外，求解角度可以选取更小值，但会增加影响系数的计算时间.

第3项近似处理是：以诱导速度势的平均来等效替代螺旋桨受力的平均.由于螺旋桨受力由伯努利方程通过速度平方项求得，即是通过诱导速度势偏导数的平方求得.因此这一平均处理可以简单用下式来描述

$ \sum\limits_{i = 1}^n \dfrac{\phi_i^2 }{n} \approx \sum\limits_{i = 1}^n \Big(\dfrac{\phi _i }{n} \Big )^2 $

(31)

其中，$n$表示个数.尽管$\phi _i $的取值不相同，但他们是同一量级的，从工程应用和数值计算的角度，可以认为这一近似处理是合理的.

综合以上的误差分析可以知道，本文提出的定常整体面元法计算公式所得到的数值结果是合理可信的.

4 总结

本文首先建立了混合式CRP的迭代面元法，随后分析混合式CRP的几何特点及桨叶和吊舱上奇点强度的变化特性，提出了将混合式CRP作为整体进行定常面元法求解的可行方法，建立了整体求解面元法的控制方程，编译完成计算程序.运用整体面元法和迭代面元法对一混合式CRP的敞水性能进行了分析，数值结果表明：与试验值相比，整体面元法与迭代面元法的敞水性能计算值在设计点的相对误差均在5%以内，计算精度能够满足工程应用的需要，同时，整体面元法比迭代面元法的计算时间要少且不需要迭代处理，这为混合式CRP的优化设计提供了高效可靠地数值手段，也为其他多对象组合式推进器的计算分析提供思路.