螺旋桨水动力性能CFD预报中预处理的程式化实现

杨琼方，郭 薇，王永生，黄 斌

（1海军工程大学 船舶与动力学院，武汉 430033；2海军驻武汉四三八厂军事代表室，武汉 430064）

1 引 言

基于螺旋桨型值表的桨叶几何建模和计算域空间网格离散是其水动力性能CFD预报的两个前提步骤。几何建模的效率和准确度以及网格划分的效率和网格质量都会直接影响到CFD计算的时间和精度。为了准确刻画桨叶导边、随边、叶梢点和叶尖几何细节，几何建模通常在CAD软件中完成，像UG和CATIA。为了合理利用有限的计算资源和尽可能高精度地捕捉上述几何细节的局部流场信息，网格划分和计算求解通常在单通道计算域内完成。其中，网格设置通常包括在四面体网格中嵌入三棱柱网格作为桨叶附面层的混合网格和直接建立在网格拓扑结构上的六面体结构化网格两种，六面体网格中桨叶近壁面流动直接由节点疏密来实现网格加密。无论是从计算资源的占用、网格数值耗散大小、网格对桨叶几何细节的映射效果，还是数值计算的精度来考虑，六面体结构化网格都要优于四面体非结构化网格[1-3]，但前提是单通道域的划分和网格拓扑结构的布置要合理，这是采用商用CFD软件进行求解时的核心技术所在。

另外，现有多数学者在进行螺旋桨水动力性能CFD计算时，为便于桨叶表面面元的布置（面元法）或者是提高叶梢附近区域的网格质量（RANS模拟），通常将叶尖部分削掉，以一个小端面来替代原叶梢点[4-6]，使几何模型在整个分析过程中都存在偏差。这样处理对积分量（推力和力矩）影响甚小，但是对于梢涡模拟以及螺旋桨空化初生判定来说影响明显。

本文以无侧斜桨和100%侧斜桨为分析对象，首先采用VC编程，结合UG软件，程序化实现了桨叶几何及其单通道域的建模，并实现了通用化。然后，利用ICEM专业网格划分软件完成六面体结构化网格划分并生成相应的命令流文件，调用该文件即可在保证网格质量的条件下实现不同螺距时调距桨或者是几何类似桨的计算域的快速便捷网格划分，从而实现螺旋桨水动力性能CFD预报中预处理的程式化操作，大幅缩短时间。最后，以4119、4381和4383桨为例，预报并校验了其无空化和空化水动力性能，分析了网格因素对CFD预报精度的影响，将对应精度最高的网格的操作命令作为网格拓扑数据库中的子样本。

2 螺旋桨桨叶三维建模的编程实现

螺旋桨桨叶通常由不同半径处的叶截面对应的伸张叶切面的二维型值来描述。型值的表述方式一般分为二种：一种是在给出伸张切面厚度分布和拱度分布规律的条件下，如NACA66（mod）+NACA α=0.8，再给出不同半径r处的弦长c、螺距P、侧斜角θs、纵斜Xm、最大厚度t0和最大拱度f0。另一种是在给出不同半径切面处螺距、弦长和随边至基线距离cT后直接给出该切面处沿弦长分布的二维坐标值。此时侧斜已由各切面参考点相对基线位置不同而体现。两种方式中对应伸张切面的二维坐标沿弦长方向排序均可以是从随边开始，或者是从导边开始。由二维伸张叶切面生成三维桨叶的过程即是得到伸张切面的逆过程，在CAD软件中对应为叶切面的横移、旋转、纵移和缠绕过程。

桨叶几何建模的编程实现即是将上述变化过程采用参数来进行描述。对于半径r处伸张切面上任一点p（x,y），横移、旋转和纵移后对应到叶截面上相应点为p4（X,Y,Z），其坐标为：

若型值以第二种方式给出，将（1）式中x-c/2-rθs/cosφ项替换成x-cT即可。若切面二维坐标从导边开始排序，将（1）式中 x-c/2-rθs/cosφ 项替换成 c/2-x-rθs/cosφ 或即可。若桨叶是左旋，将点p3相对基线镜像即可，即将 （x-c/2-rθs/cos φ）·cosφ-y·sinφ 项替换成-（x-c/2-rθs/cos φ）·cosφ+y·sinφ。 若纵倾是向船艏方向，将Xm取为负值即可。至此完成了桨叶几何建模的通用参数控制。程序图形界面如图1所示。

以DTMB4383（72°侧斜，首纵倾）桨为例，将程序得到的叶截面三维坐标导入UG软件后得到桨叶几何如图2所示。借助程序输出的切面上沿弦长等分点，将叶截面上占弦长一定百分点的区域作为与桨叶压力面和吸力面相独立的导边区和随边区，并且将一定半径截面（如0.95R）以上的桨叶部分独立作为叶尖区，见图2中所示，这样在网格划分时可以利用单独的拓扑块对这些面进行关联和控制，以提高局部流场细节控制的精度。

图1 桨叶及单通道域几何程序控制界面Fig.1 Graphical user interface of blade geometry and single passage modeling

图2 DTMB4383螺旋桨桨叶几何Fig.2 Blade geometry of DTMB4383 propeller

3 螺旋桨单通道域几何的编程实现

与泵类似，螺旋桨亦可以从属于旋转机械，可以借助GAMBIT网格划分软件中的Turbo网格生成模块或者是ANSYS TurboGrid旋转机械专业网格处理模块中的网格生成思想来完成高质量的六面体结构化网格划分[7]。但是，与泵不同的是，泵叶片与外壳间存在着固定的法向叶顶间隙距离，叶梢仍为一端面，而桨叶0.75R以上叶截面弦长迅速减小，直至叶梢变为一点，若直接套用Turbo模块来进行网格处理，要么在叶梢点关联拓扑时出现错误，要么是叶梢部分网格质量非常差，将影响计算的收敛和精度。

解决上述问题的方法有两种。一是将单通道域几何在径向方向从某个截面处分为上下两部分，如0.8R或者0.9R。该截面以下桨叶直接采用Turbo模块生成结构化网格，在径向分离截面处虚构一个圆柱面等价于泵外壳，而叶片与外壳间采用无叶顶间隙的方式处理。该截面以上桨叶作为叶尖区，采用带三棱柱网格的四面体非结构化网格来进行网格划分。最终单通道域内采用混合网格来进行空间离散，如图3所示。但该处理方式存在的问题是计算域中网格交界面增多，交界面两侧网格节点的不匹配会引入变量插值误差，且会影响计算求解的收敛效果。另一种方法是借鉴Turbo模块的拓扑结构布置特点，在ICEM CFD专业网格划分软件中采用手动拆分单通道域拓扑结构，立足于真实桨叶叶梢几何对桨叶进行六面体全结构化网格划分，并且保证网格质量满足计算的要求，处理效果与文献[1]中类同，见图3中所示。这两种解决方法实现的前提都是桨叶单通道域几何的合理、便捷地划分，其中周期界面的形状是核心因素，会直接影响到网格拓扑结构和网格质量。切削掉叶尖区的桨叶单通道域在Turbo模块中描述如图4所示，其周期界面的中间部分区域曲面形状与桨叶叶面与叶背形状相同，使单通道域形状能够反映桨叶形状，最大程度地便于叶面两侧拓扑层的对称布置，增强网格的正交性。以下重点阐述第二种方法所描述的带真实叶梢点的桨叶几何单通道域划分的程序实现。

图3 螺旋桨桨叶表面网格处理Fig.3 Surface mesh on propeller blade

图4 Turbo模块中桨叶单通道域几何描述Fig.4 Blade single passage domain description in Turbo block

仍以DTMB4383（100%侧斜）大侧斜桨为分析对象，研究结果可直接应用于无侧斜桨4381。利用图1所示界面生成满足要求的单通道域周期界面轮廓曲线如图5所示。它由螺旋线段、延伸段和样条曲线段的组合来灵活控制，程序中将该组合曲线称为包络线。包络线上连接叶截面导边和随边点的中间段为螺旋线段，用来捕捉桨叶表面形状，对应为交互界面上的“螺旋线所在圆柱长度”和“叶片数”空白框；包络线两端为3次样条曲线段，其外侧端点在周向方向上位置可调，对应交互界面上的“包络线角度调整”空白框，用于避免单通道域因桨叶大侧斜而引起的在周向上的扭曲过大，使网格的正交性降低；同时，包络线端点在轴向上的位置也可调，对应交互界面上的“包络线两端超过桨叶的长度”空白框，用来控制螺旋桨单通道域在轴向上的跨距，便于总体把握CFD计算时内旋转域的网格节点数量[8]；螺旋线段与样条曲线段之间为一段延伸距离，对应交互界面上的“延伸长度”空白框，用于实现螺旋线段与样条曲线段之间的曲率连续过渡。由不同叶截面处的包络线即可放样生成周期界面曲面，再结合进出口的圆形扇面以及桨毂和虚拟的径向圆柱外壳面，即可生成螺旋桨单通道域几何，见图6中所示。图5中，h0为对应半径r处叶截面的圆柱面在轴向上的长度，在交互界面中输入。h1为该叶截面处包络线上的螺旋线段在轴向上的投影长度。θ0和θ1分别为螺旋线段和圆柱面上的整条螺旋线在圆周方向的投影点所形成的夹角，且θ∈（0°～360）°。

图5 单通道域周期界面轮廓曲线参数描述Fig.5 Parametric control of baseline of periodic interface in single passage domain

图6 螺旋桨单通道域几何和桨叶面网格Fig.6 Propeller single passage domain geometry and blade surface mesh

依据螺距角的定义可得：

再由螺旋线的周向和轴向速度合成原理可得：

上式中，v和w分别为螺旋线对应的周向和轴向速度分量。

联立（2）式和（3）式，并令w=1 rad/s，则得到包络线上螺旋线段的参数表达式为：

上式中β为螺旋线轴向起点与x轴的夹角，见图5中所示。螺旋线起点由对应螺距为h0且过叶截面导边和随边点的螺旋线确定。

在螺旋线段向两侧延伸一定距离后，结合包络线角度α0（见图5中所示）的设定，在包络线端点轴向位置确定后，即可确定周期界面轮廓，再通过调整包络线的角度，从整体上控制单通道域在周向上的扭曲程度，即可确定单通道域几何，见图6中所示。

4 六面体结构化网格划分命令流操作

在单通道域几何确定后，在计算域内需要构建拓扑结构来关联实际几何，以映射生成域内的空间离散网格。因网格只与拓扑结构相关，与实际几何不再直接关联，所以拓扑结构对实际几何的正确、合理的映射与否，既是对实际几何准确描述的关键，也是网格质量的直接决定因素。

接下来利用ICEM CFD专业网格划分工具来完成单通道域的拓扑结构的布置和六面体网格划分。首先是整个单通道域外部轮廓由一个方块来进行关联，然后在桨叶径向方向不同叶截面处插入拓扑层，用于叶截面曲线的关联，并将方块在径向上分成多个子域以独立控制。桨叶表面近壁面采用O形结构来关联。每个拓扑层上对叶截面导边区和随边区采用单独的线来进行关联和加密控制。对叶尖区采用单独的面映射来进行关联和控制，如图7所示，最后通过微调控制线以生成满足所采用求解器的求解质量要求的网格，包括最小正则度、最小角和最大角等。

ICEM软件是基于拓扑映射几何的思想生成网格的。在得到面网格后，可通过单一命令操作方便地得到单通道域的体网格，然后依据计算需要决定是否复制生成全通道计算域网格即可。与其他应用软件一样，ICEM软件也具备脚本语言的命令流操作控制，能够对每一步的操作都用命令行的形式给予记录，以用于对操作过程的重新编绎，包括完整的重复操作、部分命令行的重复操作、对命令的重新排序操作等，如图8所示。

图7 螺旋桨单通道域网格拓扑Fig.7 Grid topology of propeller single passage domain

图8 ICEM软件操作过程的命令行控制框Fig.8 Control window of commands scripts in ICEM

图9 螺旋桨六面体结构化网格命令流应用Fig.9 Application of the command script for hexed mesh generation

命令流的二次操作主要用于以下两种场合中重新生成螺旋桨网格时的快速控制：一是计算工况的改变，如远离设计点工况时，因雷诺数发生明显变化，使得桨叶表面边界层厚度发生改变，随之要求桨叶近壁面区网格厚度分布和拓扑结构作相应的局部调整，此时可导入命令流文件中的全部命令行后再作细小的调整；二是当螺旋桨几何发生改变时，如调距桨在不同螺距时的水动力性能CFD分析，或者是桨叶形状与原来分析的桨相近的其它螺旋桨，此时只需要对该螺旋桨的真实几何作简单的清理操作后，便可载入原命令行文件实现拓扑结构和网格的快速便捷地生成，相当于实现拓扑结构与真实几何的分离操作后再统一进行关联映射即可，从而达到螺旋桨单通道域六面体结构化网格划分的程式化操作，大大缩短CFD计算时预处理的时间周期。以DTMB4383螺旋桨的命令流为基准，快速得到NSRDC4383（72°侧斜,无纵倾）和某大侧斜桨A的面网格如图9所示，也均满足计算要求。

5 螺旋桨水动力性能预报应用

首先以无侧斜桨DTMB4119为例，应用上述程式化操作流程，并且为了分析网格类型、网格拓扑结构、网格密度和网格空间节点分布规律对CFD预报精度的影响，生成4套桨叶面网格如图10所示。其中，G1为四面体加三棱柱混合网格，G2和G3分别为拓扑结构1、2映射六面体网格，G4由G3局部加密得到。G1～G4同时包含了2种网格类型、2种拓扑结构、3种网格密度和3种节点空间分布规律对RANS模拟的影响，在比较后即可确定一套相对最优的命令流脚本，从而既能保证在应用该命令流程式化操作时的时间效率，又能保证其计算精度。

为避免单通道域周期界面网格节点信息插值处理对计算精度的影响，采用全通道域进行计算。选用SST湍流模型来模拟湍流运动，计算时使用修正壁面函数来求解近壁面流动，对流项采用高阶精度格式进行离散。计算得到螺旋桨敞水性能如表1所示。可知，G4网格计算精度最高，其拓扑结构和节点分布的控制步骤即保存为该类型桨叶的命令流脚本文件。

图10 程式化生成DTMB4119螺旋桨桨叶面网格Fig.10 Blade surface mesh of DTMB4119 propeller generated by procedural operation

表1 DTMB4119螺旋桨敞水性能RANS模拟网格因素分析Tab.1 Mesh sensitive analysis for open water characteristics prediction of DTMB4119 propeller

续表1

图11 螺旋桨正、倒车敞水性能预报与校验Fig.11 Validation of ahead and backing open water characteristics of NSRDC4381 and 4383 propellers

图12 NSRDC4383螺旋桨不同负载下0.7R截面压力分布Fig.12 Pressure coefficient distribution versus loading around 0.7R section of NSRDC4383 propeller

其次以100%侧斜NSRDC4383桨和无侧斜4381桨为例，进一步说明程式化操作的实用性。计算得到两个桨正车和倒车敞水性能曲线如图11所示，均与实验值吻合较好。计算得到NSRDC4383桨0.7R截面在不同负载下的压力系数分布如图12所示，可知能很好地再现螺旋桨复杂粘性流场特征。图11中，vs为来流速度，n为转速，D为直径。再者，为分析程式化生成网格在空化模拟中的适用性，采用改进后的Sauer空化模型对4381和4383桨的片空化性能进行模拟。空化模拟数值模型详见作者文献[9-10]。计算得到2个桨的片空化形态与实验值比较如图13所示，可知均能较好地满足计算精度需求。图中J为进速系数，σ为空化数，Cp为压力系数，片空化范围对应为水蒸气体积分数0.5。详细结果分析请参见作者文献[10-12]。

6 结 论

以DTMB4383大侧斜螺旋桨为例，立足于桨叶叶梢点真实几何，详细阐述了螺旋桨水动力性能CFD分析过程中的桨叶几何建模、单通道计算域几何建模及其六面体结构化网格划分的程式化实现过程，并以DTMB4119、NSRDC4381 和 4383 桨 为 例 进行了无空化和空化水动力性能预报的应用分析，证明了程式化生成网格的适用性，可大大缩短计算的前处理时间，可进一步深化其工程应用。

图13 NSRDC4381和4383桨空化形态模拟与校验Fig.13 Validation of sheet cavity pattern of NSRDC4381 and 4383 propellers

[1]Berchiche N,Janson C E.Grid influence on the propeller oper-water performance and flow field[J].Ship Technology Research,2008,55:87-96.

[2]孙铭泽.螺旋桨敞水性能RANS模拟中网格因素的影响分析[D].武汉:海军工程大学,2009.

[3]Hsiao C T,Chahine G L.Scaling of tip vortex cavitation inception for a marine open propeller[C]//27th Symposium on Naval Hydrodynamics.Seoul,Korea,2008.

[4]叶金铭.推进器水动力性能及空泡预报的数值方法和模型试验研究[D].武汉:海军工程大学,2008.

[5]Bosschers J,Vaz G,Starke A R,et al.Computational analysis of propeller sheet cavitation and propeller-ship interaction[C]//MARINE CFD2008,RINA.Southampton,UK,2008.

[6]Vaz G.Modelling of sheet cavitation on hydrofoils and marine propellers using boundary element methods[D].PhD Dissertation,Mechanical Engineering at Universidade T’ecnica de Lisboa,2005.

[7]杨琼方,王永生,张志宏,等.叶片数对喷水推进器性能影响的计算流体动力学分析[J].机械工程学报,2009,45(6):222-228.

[8]杨琼方,王永生,黄 斌,等.融合升力线理论和雷诺时均模拟在螺旋桨设计和水动力性能预报中的应用[J].上海交通大学学报,2011,45(4):486-493.

[9]杨琼方,王永生,张志宏.螺旋桨叶截面空化模拟数值模型的改进与评估[J].北京理工大学学报（自然科学版）,2011,31(12):1401-1407.

[10]杨琼方,王永生,张志宏.螺旋桨空化崩溃性能图谱的多相流模拟[J].华中科技大学学报（自然科学版）,2012,40(2):719-723.

[11]杨琼方,王永生,张志宏.非均匀进流对螺旋桨空化水动力性能的影响[J].水动力学研究与进展（A辑）,2011,26(5):538-550.

[12]杨琼方,王永生,张志宏.螺旋桨梢涡空化初生和空化斗的数值求解[C].第二十三届全国水动力学研讨会暨第十届全国水动力学学术会议文集.北京:海军出版社,2011:467-474.