螺旋桨梢涡及梢涡空泡数值模拟

2019-05-14刘芳远傅慧萍

船舶力学 2019年4期

刘芳远，傅慧萍，李杰

（1.上海交通大学船舶海洋与建筑工程学院，上海200240；2.高新船舶与深海开发装备协同创新中心，上海200240）

0 引言

螺旋桨空泡是一种汽化空泡，即水因降压到饱和蒸汽压力导致汽化，水汽通过界面，进入气核并使之膨胀。螺旋桨桨叶上出现的空泡有涡空泡、泡状空泡、片空泡和云雾空泡4 种类型。螺旋桨空泡会对螺旋桨的水动力性能和噪声性能产生重要影响。准确的数值模拟方法可以对螺旋桨性能进行预报，并为螺旋桨的设计提供理论依据。

近廿年来，国内外对空化流提出的数值模拟方法有很多，大致可以分为势流方法和Navier-Stokes方程求解及欧拉—拉格朗日混合方法三种。尽管作为一个水动力学热点问题已经被研究多年，但囿于目前的数值模拟准确性，仍然有众多学者在空化模型及与空化相匹配的湍流模型，以及数值计算方法等方面开展相关的研究。尤其是对于螺旋桨梢涡及梢涡空泡的模拟，目前仍然是一个难点。德国波兹坦水池（SVA）曾于2011 年在汉堡主办了第二届船舶推进器国际研讨会（SMP’11），该研讨会采用盲算的方式进行比较分析，其中发布的第3 个空化案例（Case 2.3.3）在物理模型试验中出现了显而易见的梢涡空泡，数十家与会单位竟无一能算出。

国内，韩宝玉等（2011）[1]针对梢涡空化问题，应用雷诺平均的N-S 方程并结合改进的VOF（volume fraction）模型，对椭圆水翼梢涡空化特性进行了研究，采用的湍流模型为代数雷诺应力模型。为了准确模拟梢涡涡核内系统旋转和流线曲率影响，将旋转和曲率修正方法应用于湍流模型中。但是对水翼的水动力没有进行定量的校核，只是给出了所谓消失空泡数。Shi 等（2015）[2]应用分离涡模型（DES）对螺旋桨尾流中梢涡流场分布进行了数值研究，但还没有算过空化。Ji 等（2013）[3]采用大涡模拟（LES）和基于质量输运方程的空化模型预测了一个3 维扭曲水翼的非定常空化流。辛公正（2014）[4]采用数值RANS 方法模拟了水翼和螺旋桨的梢涡流动，并不求解空泡，而是寻找最小压力系数与起始空泡数之间的关系。Liu 等（2016）[5]基于OpenFOAM 平台，比较了线性与非线性湍流模型对螺旋桨梢涡空化的数值模拟能力，发现非线性k-ε 模型能够更加准确地捕捉到梢涡空化现象。但没有同时对螺旋桨的推力系数和扭矩系数进行校核。将螺旋桨空泡水动力性能计算做得比较好的有冯学梅等（2012）[6]，他们基于FLUENT 软件，采用多块结构化网格对业内选作标桨的E779A 桨和PPTC 桨进行了均匀来流下的空泡数值模拟。

目前主流的空化计算方法是结合空化模型，求解两相流的RANS 方程。欧拉—拉格朗日混合方法主要用于研究空化初生[7]。Hsiao 等（2005）[8]以水翼为研究对象，采用特定区域发放气核的研究方法，考虑气核尺寸分布进行梢涡空化初生噪声尺度效应研究。真正与螺旋桨相关是在2008 年第27 届海军水动力学会议上，他们发表了敞水螺旋桨的梢涡空化初生尺度效应研究[9]。最近的研究成果是关于螺旋桨及气体扩散作用对气核分布的影响[10]。国内熊鹰等（2013）[11]也采用这种方法对螺旋桨空化初生进行尺度效应研究，但均不直接模拟空泡。

本文基于螺旋桨空泡数值模拟的现状，以PPTC 桨为研究对象，通过一种新颖的网格划分方式，求解均匀来流下的螺旋桨空泡及其水动力。目标是同时求解出与试验相符的空泡形状，以及螺旋桨的推力和扭矩系数。

1 控制方程

采用多相流模型中的均质混合流模型，此时气液两相流被当作包含微气泡的单相连续介质，通过引进气相体积分数α 来定义混合相的密度ρm，仅求解一套控制方程。基于混合密度的均质混合流的连续方程为

动量方程为

式中：下标l、v 分别指液相和气相；μm为混合流的动力粘性系数，定义与混合密度类似；μt,m为湍流引起的涡粘性系数，需要对混合流采用湍流模型进行求解。本文采用RNG k-ε 两方程模型求解涡粘系数，该模型在求解旋转流动和剪切流时有优势，壁面第一层网格节点的无量纲距离y+要求在30～100 之间。

关于蒸汽质量分数f 的控制方程为

式中：Re和Rc分别为由汽化和凝结引起的相变率；f 按下式定义：

根据对（5）式右端源项的建模方式，可以区分出不同的空化模型。常用的有Zwart-Gerber-Belamri模型和Schnerr-Sauer 模型。本文采用Zwart 模型：

当p≤pv时，

当p＞pv时，

式中：pv为水的饱和蒸气压，与温度相关；气核半径RB默认取值1 μm；αnuc为水中所含气核的体积分数，默认取值5×10-4；Ce和Cc分别为汽化系数和凝结系数，默认取值50 和0.01。

2 研究对象与网格

PPTC 桨是指德国波兹坦水池（SVA）设计的可调螺距桨VP1304，也是第二届和第四届船舶推进器国际研讨会（SMP’11 和SMP’15）的算例桨，其在均流和斜流情形下的敞水性能、伴流分布、空化性能及压力脉动等试验数据通过互联网为全世界的水动力学研究人员所共享。本文即以SMP’11 的算例桨PPTC’11 为研究对象，对研讨会上发布的空化案例Case 2.3 进行数值模拟，以校验空化模型及相应的空化流场计算方法。

图1 研究对象与计算域Fig.1 Computed object and computational field

PPTC 为直径D=250 mm 的5 叶右旋可调螺距桨，这使得桨叶在导边和随边处与桨毂之间存在0.3 mm 的间隙，数值模拟时予以忽略。SMP’11 的空泡测试在SVA 的空泡筒K15A 中完成，轴无倾角。分别在15、20、25 r/s 三个转速下对敞水性能进行了测试，以探讨雷诺数对敞水性能的影响；敞水试验在高压时进行，以避免空泡。图1 为SVA 提供的PPTC 桨模与空泡筒试验段几何模型。为了更好地模拟均匀定常来流下桨的水动力性能，按SMP’11 的要求，计算时将试验段变换成等截面的圆柱段，其它如桨模安装位置和试验段长度保持不变。

依据流动最大雷诺数及湍流模型的适用范围估算了壁面第一层网格节点的距离，约为2.4D‰，据此在近壁区进行边界层网格加密。尤其是对桨叶梢部的梢涡脱出区域进行了区域划分及域内网格加密，以保证涡核区域内径向15 个网格节点的密度（如图2 所示）。生成全流道计算域网格，共8 560 925 个混合型单元。

图2 梢部网格加密Fig.2 Local grid refinement of the tip vicinity

3 无空化流场计算

PPTC 桨被设计专门用于产生梢涡，并展开广泛的模型试验，用以校验螺旋桨的各种势流或粘流分析程序。在无空化情形下，SVA 在空泡筒对3 个转速下的PPTC 敞水性能进行了测试。本章选择3条敞水性能曲线中进速系数等于1 附近的3 个试验点进行数值模拟。计算条件同试验条件（见表1），表中pA为环境压力。旨在获得可靠的空化流场初始解，并对网格划分及计算方法进行无空化计算校验。

表1 无空化计算条件Tab.1 Computational conditions w/o CAV

计算在单相流框架下进行。湍流模型采用旋转流适用的RNG k-ε 模型，近壁区的处理采用标准壁面函数。压力速度耦合采用PISO（Pressure-Implicit with Splitting of Operators）算法，即隐式算子分割算法，该算法适用于高倾斜度网格。空间离散均采用二阶以上高精度离散方式：梯度扩散项采用基于单元体的格林—高斯离散；压力采用高旋流适用的PRESTO 格式；动量方程和湍流方程均采用适用于四面体网格的二阶迎风格式。表2 给出了推力系数KT和扭矩系数KQ的计算结果。由表可见：KT和KQ的计算误差在3.0%～5.4%以内；随着转速提高（n=15～25 r/s），雷诺数增大（Re0.7=1.0×106～1.6×106，以0.7R 半径处的合速度及桨叶弦长为特征量），计算精度略有降低。

表2 无空化计算结果与误差分析Tab.2 Results and error analysis w/o CAV

为了准确定义并描述旋转的流场，需要对漩涡结构进行定义。关于漩涡的判别法则有许多，比较流行的是涡量准则、Q 准则、Δ 准则和λ2准则。以上这些判别方式，都定义了一个可以逐点计算数值的函数，然后判定每一点位于漩涡内或者漩涡外。尽管判定方式不少，但是对涡进行准确而严格的定义却是较为困难的。通常，漩涡被认为是涡量较高的区域，但是对该区域的准确阈值却没有严格的定义，随着流动性质的不同，阈值的判定将会严重影响结果。而且，在没有漩涡的强剪切流动中涡量也可能较大，所以，利用涡量准则来判定漩涡是不够准确的。利用流线来定义涡并进行描述也是不可靠的，因为即使在简单的匀速平移运动中，流线的拓扑结构也容易被改变。而后三种方法（Q 准则、Δ 准则和λ2准则）则在漩涡判断的方式上具有等价性。本文采用Q 准则对漩涡进行定义和直观的表达。

Q 准则由Hunt 等在1988 年提出[12]，它定义流场中速度梯度张量▽V 的第二不变量Q 具有正值的区域为漩涡。另外，它要求漩涡区域的压强低于周围的压强。Q 定义为：

式中：Ω 为涡量，S 为变形率，皆为二阶张量。符号“||||”表示张量的二范数。S 和Ω 分别代表了流场中一点的变形和旋转。Q 准则反映了流场中一个流体微团旋转和变形之间的一种平衡，Q＞0 则反映了旋转在流动中占据统治地位，同时也体现了该旋转区域的压力为附近区域的最小值（极小值）。图3 给出了以转速为特征量的无量纲等值面图。由图可见，采取本文所述网格生成方法得到的梢涡流场较之采用传统的不对梢涡进行区域划分的网格计算结果要显著得多，所捕获的梢涡长度要长得多：梢涡从导边后缘拖出，并向螺旋桨后方移动并逐渐耗散；当值增大，梢涡区域向涡核中心线收缩。

图3 等值面Fig.3 ISO surfaces of

4 空化流场计算

螺旋桨空化性能预报对于分析螺旋桨的设计工况和非设计工况都意义重大。SMP’11 对工作状态下的螺旋桨空化性能进行了观测，均在转速n=25 r/s 下进行，按不同的推力系数和空化数分为3 个案例，分别标注为Case 2.3.1、Case 2.3.2 和Case 2.3.3，试验条件见表3。其中，Case 2.3.1 和Case 2.3.3 为非设计工况，Case 2.3.2 为设计工况。本章将对Case 2.3.1 和Case 2.3.2 进行数值模拟。计算条件同试验条件。

表3 空化计算条件Tab.3 Computational conditions with CAV

表3 中σn为以旋转线速度为特征速度的空化数，定义如下：由此可以换算出环境压力pA。

先计算无空化情形，此时环境压力可以取作一个大气压。校验推力系数（扭矩系数缺失），并作为空化流场的初始解。然后将计算切换至均质混合流框架下，改变环境压力后，开启空化模型并设置空化参数。此时，与无空化计算相比，增添了一个气相体积分数的控制方程。该方程的空间离散选用QUICK（Quadratic upwind interpolation）格式。此格式用于六面体网格具有三阶精度，用于四面体网格或混合网格时只具有二阶精度。其他在无空化计算时已有的设置则保持不变。

4.1 算例1（Case 2.3.1）

该算例属于非设计工况，载荷偏高情形。先进行无空化计算（此时pA=101.325 kPa），将计算得到的推力系数与试验值进行比较。由表4 可见，推力系数的计算误差为2.9%。重新设置环境压力（此时pA=42.207 kPa），开启空化计算。由表4 可见，空化将引起推力和扭矩系数下降；KT和KQ的计算误差均在3.0%以内。

表4 计算结果与误差分析（Case 2.3.1）Tab.4 Results and error analysis for Case 2.3.1

空泡形态一般采用汽相体积分数琢的等值面来表示（0<琢≤1），因此空泡的大小长短与琢的取值相关。图4 给出了空泡形态的计算结果（琢=0.2）与试验结果（EFD）对比。其中，“SS”表示吸力面；“PS”表示压力面；“SVSS”表示吸力面侧视图。图4（a）的试验结果显示，此工况下，除了发生显著的梢涡空化外，吸力面叶根处也发生了片空化。图4（b）～（d）的3 个视图全方位地展示了梢涡空泡、吸力面叶根处的片空泡及毂涡空泡，计算与试验符合较好。

图4 空泡形态计算与试验比较Fig.4 Observed and calculated cavitation(Case 2.3.1)

图5 等值面Fig.5 ISO surfaces of (Case 2.3.1)

4.2 算例2（Case 2.3.2）

此工况为设计工况，中等载荷情形。先进行无空化计算（此时pA=101.325 kPa），将计算得到的推力系数与试验值进行比较。由表5 可见，推力系数的计算误差为4.5%。重新设置环境压力（此时pA=30.528 kPa），开启空化计算。由表5 可见：空化将引起推力和扭矩系数下降；KT的计算误差为2.4%，KQ为4.5%。

表5 计算结果与误差分析（Case 2.3.2）Tab.5 Results and error analysis for Case 2.3.2

图6 给出了空泡形态的计算结果与试验结果对比。图6（a）的试验结果显示，此工况下的梢涡空泡在靠近随边处发生并拖出，吸力面叶根处发生了较大面积的片空化，并伴随有泡状空泡。图6（b）～（d）分3 个视图给出了数值模拟结果（琢=0.2）。由图可见，除了梢涡空泡、吸力面叶根处的片空泡，压力面叶根处也发生了较小的片空泡，计算与试验符合较好。此工况下，没有出现毂涡空泡。

图6 空泡形态计算与试验比较Fig.6 Observed and calculated cavitation(Case 2.3.2)

图7 等值面Fig.7 ISO surfaces of (Case 2.3.2)

4.3 水中含气率对空化计算的影响

数值实验表明，当水中气核体积分数琢nuc在一定范围内取不同值时，对空泡形态的影响并不显著，却会影响宏观量—推力和扭矩系数的计算。表6 给出了算例Case 2.3.2 中琢nuc对KT和KQ的影响。由表可见，琢nuc越大，KT和KQ有降低的趋势，从而更接近试验值，但此时的计算容易发散。如何让琢nuc的取值与模型试验中的水中含气量等环境变量相关联，或者说，进一步探究或改良现有空化模型是下一步研究需要解决的问题。