螺旋桨敞水性能数值预报方法研究

孔为平 王金宝 冯学梅

（中国船舶及海洋工程设计研究院 上海200011）

0 引 言

螺旋桨敞水性能是船舶快速性能的一个基本组成部分，与船舶阻力性能预报一样，其准确的定量预报是迄今为止船舶计算流体力学（CFD）发展比较受关注的一个方面。

螺旋桨敞水性能、尾流场的数值计算方法的开发与验证，构成了目前螺旋桨水动力性能研究内容的主要部分，受到历届ITTC、欧盟框架计划项目（FP5-LEADING EDGE项目、FP6-VIRTUE项目）以及许多计算流体力学工作组（CFD Workshop）的持续关注。

本部分计算对象为E779A型螺旋桨（E779A桨）和30.8万吨超级油轮螺旋桨（30.8万吨VLCC桨）,进行了敞水性能计算,并考察了网格、计算域、压力离散格式等对敞水性能预报结果的影响。其中，E779A桨为意大利船模水池中心INSEAN的标桨，其敞水性能测量数据也从该水池获得。INSEAN围绕E779A桨敞水性能和空化特性进行了一系列完整细致的试验研究，欧盟第六框架计划（EU-FP6，即6th Framework Program of the European Union）的项目 VIRTUE（The Virtual Tank Utility in Europe）下的第四工作包 “数值空泡水池”（WP4，The Numerical Cavitation Tank）也选择E779A桨作为研究对象，组织多个单位，采用多种求解器，对处于均匀来流中和船后尾流场中螺旋桨的空泡现象进行计算比较，并分别在2007年和2008年针对该桨组织了两次相关的CFD Workshop［1］。 30.8万吨VLCC螺旋桨的试验数据则由中国船舶及海洋工程设计研究院水池提供。

冯学梅等曾针对E779A及PPTC桨，成功地进行了桨的敞水、尾流、空泡计算［2-3］。计算结果与试验结果都相当一致，成为与会19份计算结果中最为突出者。在此基础上，本文对其中一些数值计算作了分析和调整，以期使相关的算法更为合理，并形成更趋成熟、快捷的流程。

1 E779A螺旋桨敞水性能计算

1.1 计算对象和工况

欧盟框架计划项目将E779A桨的测量数据资料视为机密，但随时日推延，E779A桨水动力性能数值和水声测量数据已有部分公开，并被广泛用作基准校验的对象，进行数值计算的比较研究。该桨带有纵斜，为固定螺距的右旋四叶桨。桨模及其主参数见图1和表1。

图1 E779A桨模

表1 E779A桨模主要参数

表2 E779A桨敞水性能计算工况

1.2 计算方案

鉴于来流均匀和螺旋桨几何上的周期性，为进行敞水计算，只需取单个桨叶所在的单通道（对于4叶桨而言，周向跨90°）作为计算域即可进行计算分析，由此可节省计算时间,提高计算效率。计算域的进口面、出口面均为90°的扇面，与桨轴垂直，彼此沿周向有所错位；以桨毂和轴为内边界面；外边界面取在圆柱体表面上，其直径为螺旋桨直径的数倍；周向侧面由进、出口面附近的子午面和连贯其间的螺纹面组成。本部分敞水计算中，取进口在上游2.11D处，出口在下游3.08D处，外边界所在圆柱面半径为1.54D，桨轴为从进口一直延伸至出口的圆柱，其直径为0.123 5D。来流方向与z轴方向相反，具体参照图2。

图2 计算区域及坐标

这里的螺旋桨敞水计算中，整个计算区域处于一个以转速n绕桨轴作旋转运动的参考坐标系中，桨叶相对此参考坐标系静止不动。因此可选用FLUENT软件提供的运动参考坐标系模型（即MRF模型）［4］。

采用基于压力耦合的粘性求解器。压力与速度的耦合引用SIMPLE方法，湍流模型为SSTκ-ω模型。对压力项，考察采用了标准（Standard）和二阶差分（Second Order）格式离散。动量方程其余项、湍流模型方程湍流动能（Turbulence Kinetic Energy）和湍流耗散率（Specific Dissipation Rate）项均采用二阶迎风格式（Second Order Upwind）作离散。

在进口边界处设置为速度进口条件，给定均匀来流的各速度分量；出口边界给定表压为0，即与参考点静压相等；外边界同样设为速度进口边界；而两个周向侧面置为旋转周期性边界。

1.3 计算结果及分析

以下讨论网格选择、计算域的考察、压力标准离散格式和压力二阶离散格式对计算收敛情况影响的比较以及收敛过程、时间的考察等工作，最后给出敞水计算的相对偏差和相关曲线。因E779A桨的敞水性能测量数据从意大利INSEAN水池获得，并未授权公开，故以下计算所得敞水性能数据与试验结果的比较都以相对偏差给出。

1.3.1 网格选择

为考察网格对敞水计算结果的影响，对E779A桨生成3套网格进行计算比较。在CAD软件Gambit中生成约84万个四面体非结构化网格（G84），在Gambit和TGrid中生成约103万个和125万个带边界层网格的混合网格（G103和G125），其中第一层边界层网格离壁距离为0.000 022D。网格特征信息见表3；网格G84、网格G103与网格G125之比较见图3。

表3 网格特征

图3 网格比较

针对上述3套网格，选择进速系数J=0.397进行 计算考察（见表4），其中压力项离散采用二阶格式。

表4 J=0.397时，不同网格计算结果与试验值相对偏差比较

表4中：

Δ表示相对偏差；

对叶面压力分布进行了考察，可以看出：压力量级达到1e4时，导边及外缘压力梯度非常大，参见图4。所以，为更好捕捉真实的流动特性，桨叶边缘应加密处理。相比较而言，从网格分布和边界层布置来看，G125都更适应该敞水计算。

图4 E7779A桨叶面压力分布

由表4可以看出：G84没用边界层，Kq相差较大；G103计算出的Kt、Kq相对偏差都超过了 3%；G125网格的叶面轮廓进行了加密，计算结果都在2%以内。下面将对G125在其他进速系数条件下进行详细考察。

1.3.2 计算域影响

前面提到本次敞水计算区域。冯学梅等人针对该E779A桨的敞水计算在计算域上采用以下方法：取进口在上游2D处，出口在下游6D处，外边界所在圆柱面半径为3D，桨轴为从进口一直延伸至出口的圆柱，其直径为0.123 5D。计算域相应位置明显大于本计算域，为考察计算的可信度，在该计算域中提取出本计算域边界相应位置上的压力分布等相关量。

由图5和图6可以看出，在z=-0.7 m处压力速度等值线都几乎是周向均匀分布，并且沿径向只在靠近壁面的小范围内变化，压力量级基本为1e2，与桨叶上的压力量级1e4相比，若取其为0，其影响应可忽略不计。速度分布也都在1.1 m/s左右，这与在J=0.397下所给定的速度入口条件基本吻合。可以认为在z=-0.7 m处已经达到了无穷远条件。

图5 z=-0.7 m时，桨叶面上的压力分布（左：等值线；右：沿径向分布）

图6 z=-0.7 m时，桨叶面上的速度分布（左：等值线；右：沿径向分布）

同样，由图7可以看出，在侧边界（r=0.35 m）上，无论是压力分布还是速度分布，都和所给定的侧面边界条件基本吻合。因此，本文所取计算域能满足计算要求。

图7 r=0.35 m时，压力和速度分布

1.3.3 压力离散格式

对于压力二阶离散格式和标准离散格式的考察。标准离散格式是用计算出的单元体中心压力和速度来插值到单元面上。二阶压力格式直接计算出单元面上压力和速度值。二阶格式一般适用于高速旋转、高雷诺数、复杂几何曲面的流动计算。而标准离散格式，则有较好的数值性能，比如收敛稳定。

由图8和图9可以看出，当垂轴量度幅值相同时（一般在迭代500次左右），残差、推力和扭矩就开始进入稳定状态；迭代1 000次时的推力和扭矩则非常稳定。压力标准离散格式收敛更好，残差量级都降到1e-3以下，而且推力、扭矩都收敛于一条直线。相比之下，压力二阶离散格式收敛不太理想，连续方程和Omega残差都在1e-3以上，推力和扭矩呈波动式收敛，数值波动约为±0.68%。从数值上来看，压力标准格式推力、扭矩绝对值略低于压力二阶格式；与试验数据相比较，压力标准格式更接近试验值。综合看来，对于本敞水计算雷诺数1.9e6，转速11.79 rad/s，使用标准压力满足计算精度要求，并且能收到更好的收敛效果。

图8 J=0.397时的压力标准离散格式收敛曲线

图9 J=0.397时的压力二阶离散格式收敛曲线

1.3.4 收敛时间

本文对该敞水计算的收敛时间进行了考察分析，选择了几个进速系数，给出了相关计算收敛历史曲线如图 10～14。

图10 J=0.397时的收敛曲线

图11 J=0.498时的收敛曲线

图12 J=0.695时的收敛曲线

图13 J=0.845时的收敛曲线

图14 J=0.946时的收敛曲线

表5 收敛结果整体对比（500～1 200步或1 500步）

表6 相关迭代步和试验值结果对比

综合上述图表可知，本敞水计算经500步后，无论推力T还是转矩Q都已渐趋稳定，其后的变化均达到千分之一左右甚至更小的相对偏差。使用6核，500步计算时间大约为45 min，即最少45 min可以完成一个进速系数点计算。

1.3.5 结果比较和分析

由表7和图15可见（其中cfd为计算值，efd为试验值），除去进速系数最大的J=0.995和J=1.094以外，G125的推力系数、扭矩系数以及敞水效率的相对偏差基本保持在2.5%以内。其中，边界层网格G125的敞水性能曲线几乎与试验曲线重叠。

表7 敞水性能计算结果相对偏差

图15 敞水曲线比较

由图15可知，敞水效率在J=0.95时最大；而计算结果相对试验结果偏差过大的情况，则出现在J=0.995和J=1.094。此处，因推力、扭矩本身很小，故无论计算还是试验的难度均较大。

2 30.8万吨VLCC桨敞水性能计算

30.8 万吨VLCC桨为备用桨，该桨没有纵斜，为可变螺距的右旋四叶桨（见图16）。

图16 30.8万吨VLCC桨模

其主要参数为：

桨模直径D=0.2 m，螺距比P0.7/D=0.667，盘面比Ae/Ao=0.54，毂径比 d/D=0.146，纵斜角=0°，叶数=4，右旋向。

其敞水性能计算工况为：

进速系数 J=0.0～0.8，转速 n=25.947 r/s，水温为10.5℃，水运动粘性系数v=1.295 3×10-6m2/s，水密度 ρ=1 000.6 kg/m3。

根据E779A桨的计算经历，对30.8万吨VLCC桨的计算不再考察网格的影响，而是直接采用带边界层网格的网格布置。设置叶片上最小网格尺寸为0.001D，第一层边界层网格离壁大小为0.000 02D，在Gambit和TGrid软件中生成约80万个带边界层的四面体非结构化网格。计算方案与边界条件的设置与E779A桨相同。

由表8和图17可以看出计算所得敞水性能与试验结果的比较。可见，敞水性能曲线与相应的试验曲线相当吻合，当进速系数J=0.3～0.6时，无论推力、扭矩还是效率偏差都在1%以内。再次验证了这套敞水性能数值模拟方案的可行性和可靠性。

表8 30.8万吨VLCC桨敞水性能结果比较

图17 30.8万吨VLCC桨敞水性能曲线

3 结 论

在敞水计算中，为得到精确结果，除了边界条件和计算域合理外，网格布置也很关键，应该根据流动特性，合理密布关键区域尤其是桨叶上的网格，以及采用边界层网格技术，而标准压力插值格式的运用，也使计算结果精确且效率更高。

4 致 谢

本文在写作过程中，得到了蔡荣泉研究员的悉心指导，也得益于于海、吴琼等同事的热心帮助，在此表示衷心感谢。

［1］MARNET-CFD. ［2012-03-10］.Final report and review of the state-of the-art in the application of CFD in the Maritime and Offshore Industries［EB/OL］.https：//pronet.wsatkins.co.uk/marnet/.

［2］冯学梅.水面舰船水动力快速性能数值研究［D］.上海：上海交通大学，2001.

［3］蔡荣泉，冯学梅.浅谈舰船计算流体力学（CFD）实用化［J］.船舶，2012，23（2）：75-84.

［4］Fluent 6.1 Tutorial Guide［Z］.2003.