螺旋桨梢涡空化数值模拟的网格策略研究

张 彬，张志荣

（中国船舶科学研究中心，江苏 无锡 214082）

0 引 言

螺旋桨梢涡是流体在压差作用下从桨叶压力面流向吸力面，经过梢部时发生流动分离而产生的旋涡。一般认为当螺旋桨梢涡涡心压力低于临界压力时将出现梢涡空化，严重影响舰船振动、噪声等性能，因此工程上对于梢涡空化的准确预报有着迫切需求。

研究者们通过PIV和LDV技术对螺旋桨梢涡空化进行了大量试验观测和流场测量,为梢涡空化研究提供了丰富的试验数据[1-3]。随着CFD 的发展和计算能力的提高，数值模拟在梢涡空化研究中发挥着越来越重要的作用。但是由于梢涡涡心的尺寸小且压力梯度大，目前较难精确地模拟梢涡空化。提高螺旋桨梢涡空化数值模拟的精度有许多方法，例如在涡心中设置较多的网格节点，采用合适的湍流处理方法和空化模型等。

Chow 等[4]通过计算验证认为，梢涡数值模拟中至少需要15～20 个网格点覆盖梢涡涡心。Tuomas等[5]分别使用标准k-ε湍流模型和k-ωSST湍流模型，基于三套网格对PPTC螺旋桨进行了梢涡空化数值模拟，发现网格分辨率对梢涡流场数值模拟精度的影响显著大于湍流模型的影响；刘芳远等[6]使用RNGk-ε湍流模型和Zwart-Gerber-Belamri 空化模型对PPTC 桨进行了梢涡数值模拟，通过局部加密梢涡脱泄区域的网格，成功模拟了梢涡空化，推力、扭矩系数与试验接近；胡建等[7]采用大涡模拟（large eddy simulations,LES）方法，基于螺旋加密网格对E779A螺旋桨进行了梢涡空化数值模拟，展示了精细的梢涡流场，梢涡空化形态与试验结果接近，但梢涡空化区域较小；Lloyd等[8]分别使用RANS方法和DDES方法对E779A螺旋桨进行了梢涡空化数值模拟，对其中一个桨叶基于无量纲Q准则判据进行了识别并对梢涡区域进行了加密，发现结合延迟分离涡模拟DDES方法和自适应网格加密方法能够明显降低梢涡涡心的最低压力系数，梢涡空化模拟结果与试验更吻合；Yilmaz等[9]使用STAR-CCM+软件，采用LES方法和Schnerr-Sauer空化模型对E779A螺旋桨进行了空化数值模拟，并在螺旋加密初始网格的基础上，通过自适应网格加密方法对绝对压力接近饱和蒸气压的网格进行了加密，计算所得梢涡空化形态与试验结果比较接近。分析比较文献[7～9]可以得出这一结论，即不同的湍流模拟方法（RANS、DDES和LES）对螺旋桨梢涡空化数值模拟精度有较大影响，但小于网格分辨率的影响。

本文在国内外研究的基础上，以E779A 桨为研究对象，基于开源CFD 软件OpenFOAM 的inter-PhaseChangeFoam 求解器，通过局部网格加密和自适应网格加密方法，精确加密梢涡所处的网格，以相对较少的计算量准确模拟近尾流场的梢涡空化。

1 数值模拟方法

1.1 控制方程

对于空化流求解，基于Volume of Fluid（VOF）方法将两相流可看作以汽液两相按一定比例混合而成的单相流模型。混合相的密度ρm与动力粘性系数μm可以通过如下公式得到：

式中，αl为液相体积分数，αv为汽相体积分数，ρl为液相的密度，ρv为汽相的密度，μl为液相的动力粘性系数，μv为汽相的动力粘性系数。

假定各项流体不可压缩，不考虑流体的热传递效应和体积力，则连续性方程和动量方程为

由于直接数值模拟计算量过大，引入RANS方法得到时均流动的连续性方程和动量方程为

式中，ui和uj为未经雷诺平均的速度矢量，和为雷诺平均后的速度矢量，p为未经雷诺平均的压力，pˉ为雷诺平均后的压力。式（7）中最后一项为雷诺应力项，通过引入k-ωSST湍流模型求解。

为了求解体积分数，引入Schnerr-Sauer 空化模型[10],由式（3）可知只需要求解其中一项的体积分数。汽相体积分数输运方程为

联立式（1）、（6）和（8）可得

Schnerr-Sauer空化模型基于R-P方程和球形气泡半径公式推导得到∂αv/∂t，最终冷凝率和汽化率的表达式为

式中，Cc和Cv均为经验系数，本文都设置为默认值1。pv为液体的饱和蒸气压,本文设置为2350 Pa。气泡半径R满足：

式中，n0为单位液体内的气泡数，设置为1.6×1013。

1.2 几何模型

本文研究对象为INSEAN E779A模型桨，拥有丰富的螺旋桨试验研究数据。E779A模型桨为四叶桨，直径D为0.227 27 m。主要参数如表1所示。坐标系原点设置在桨中心点，各轴方向如图1所示。

表1 E779A模型桨主要参数Tab.1 Particulars of the propeller

1.3 计算域与工况设置

计算域为2.942D的圆柱体，计算域左侧为速度入口边界，到桨中心的距离为1.25D；右侧为压力出口边界，到桨中心的距离为4D。在螺旋桨附近区域设置圆柱体旋转域，在旋转域内使用滑移网格模拟螺旋桨的旋转，其它区域网格静止，最后在区域交界面AMI上进行流场数据交互。旋转域直径为1.32D，其左端面和右端面到桨中心的距离分别为0.25D和0.65D。计算域示意图如图2所示。

为了与试验结果[11]比较，数值模拟在进速系数J=0.71 和空化数σ=1.763 工况下进行，来流速度Uinf=5.808 m/s，转速n=36r/s，时间步长设置为Δt=1.93×10-5s，对应螺旋桨旋转角度0.25°。空化数σ定义如下：

式中，p0为静压，pv为液体的饱和蒸汽压，ρl为液体密度。

2 网格策略

2.1 网格局部加密

2.1.1 网格划分

首先采用六面体网格划分非结构网格，其体网格变化率为2。然后对桨叶表面、导随边和梢部的网格均匀加密，其中导随边和梢部附近的网格尺寸为桨叶表面网格尺寸的1/4。接着在桨叶和桨毂表面添加边界层网格，其高度变化率为1.2，第一层边界层网格高度约为桨叶表面网格尺寸的1/10。最后对螺旋桨附近区域的全部网格进行加密，网格尺寸为桨叶表面网格尺寸的2 倍。桨叶表面网格和螺旋桨附近网格如图3所示。

若对近尾流场全部网格进一步精细加密，将导致网格数量过大。为了在精确模拟梢涡空化的同时减少不必要的计算量，对梢涡涡管所在圆环柱内的两个区域进行局部加密，如图4（a）所示。

图4（a）中蓝色的网格加密区域记为A，每个桨叶附近的红色区域为以螺旋桨中心点为圆心，扇环内径r1=0.45D，扇环外径R1=0.51D，旋转角为36°，轴向从x=0拉伸到x=0.06D的扇环区域。在A区域中加密网格的原因是桨叶附近网格正交性较差而不适用于自适应网格加密，因此需要提前在桨叶附近的梢涡区域布置足够密的网格。

图4（b）中红色的网格加密区域记为B，其为x轴为对称轴，圆环内径r2=0.42D，圆环外径R2=0.51D，轴向从x=0 拉伸到x=0.25D的圆环柱区域。本文为了减少计算量仅加密至轴向x=0.25D，若计算资源充足可尽量将该区域向后延伸。

为了比较说明A和B区域的网格分辨率对梢涡空化数值模拟的影响，根据上述网格划分策略共划分了四套网格，网格信息如表2所示，加密后的G3网格如图4（b）所示。

表2 网格信息Tab.2 Information of different grids

2.1.2 计算结果与分析

基于上述四套网格对螺旋桨梢涡空化进行数值模拟。表3列出了各套网格数值模拟结果的推力系数KT和扭矩系数KQ，结果显示均与试验结果非常接近，证明本文采取的网格划分方法和数值模拟方法比较可靠。

表3 各套网格计算结果的推力系数和扭矩系数Tab.3 Thrust coefficients and torque coefficients of calculation results of each grid

采用汽体体积分数αv=0.2 的等值面表示梢涡空化，结果显示数值模拟的梢涡空化形态均与试验比较吻合，如图5所示。对比四套网格的计算结果，可以观察到梢涡空化区域随着网格的增加而不断增大，其中G3与G4的计算结果差距相对较小。但是数值模拟的梢涡空化长度与试验的差距均很大，仍需要进一步提高梢涡区域的网格分辨率。比较G2与G3的数值模拟结果可知，对A区域的局部加密有利于提高梢涡空化数值模拟精度。

分别截取包含A和B区域在内的两个截面，其无量纲轴向涡量分布如图6 所示。图6（a）显示，随着网格分辨率的提高，梢涡上方的轴向涡量明显增大，而下方的涡量变化很小。这表明梢部端点附近的流动状况更加复杂，速度压力变化较大，需要更高的网格分辨率减少数值耗散。图6（b）中只有G1的轴向涡量相对较小，而其它网格的计算结果差距很小，这是因为G1 网格的整体网格分辨率以及其它三套网格的B区域网格分辨率均较低。比较两个截面的无量纲轴向涡量分布可以发现，随着梢涡从A区域发展到B区域，其轴向涡量衰减十分迅速，这表明需要在B区域中布置足够密的网格以提高数值模拟精度。

从图5与图6可以观察到，G3与G4的计算结果差距较小，为了节约计算资源以及提高计算效率，选择G3网格作为下一步自适应网格加密的初始网格。

2.2 自适应网格加密

2.2.1 网格划分

若对B区域的所有网格进行加密，将大幅增加网格数量。本文采用基于Q准则的自适应网格加密方法，精确提高B区域中梢涡网格的网格分辨率，从而减少不必要的计算量。Q为速度梯度张量的第二不变量，计算公式如下：

式中，Ω为涡张量，S为应变率张量。

自适应网格加密流程如图7所示。首先在一定工况下对局部加密后的网格进行螺旋桨梢涡空化数值模拟，再根据Q方法识别计算结果中的梢涡结构，然后使用topoSet工具标记Q取值大于一定值的网格，最后利用refineMesh 工具对标记的网格进行八叉树加密。可重复自适应加密步骤，直到梢涡区域网格分辨率满足精细数值模拟的需求。螺旋桨近尾流场涡结构复杂，为了避免加密近尾流场中的无关涡结构，在读取并标记网格时Q的最小阈值Qlimit应当设置得较大。

基于G3的梢涡空化数值模拟结果，标记并加密了计算域内轴向x/D=0.06到x/D=0.25的区域内Q＞2.0×106的网格，加密后的网格记作G5，网格总数为2.13×107，自适应加密后的网格如图8所示。

除Q准则之外，自适应网格加密方法也可以通过Δ 准则、λ2准则和Ω方法等识别并加密梢涡区域的网格，或者通过气体体积分数或者绝对压力值识别并加密梢涡空化区域的网格。

本文所用的自适应网格加密并非严格意义上的自适应网格加密，不是在每个时间步后都根据某一个物理量进行网格加密。这是因为没有一个物理量能仅仅表示梢涡结构，例如各种涡判定方法等识别的区域都会包括桨叶表面附近的网格，但靠近桨叶的网格正交性较差而不适用于OpenFOAM 的八叉树加密。因此舍弃了基于dynamicMesh 工具实现的严格意义上的自适应网格加密，换而使用基于topoSet 工具和refineMesh工具实现的“静态”的自适应网格加密。

2.2.2 计算结果与分析

基于自适应网格加密后的G5 网格进行了梢涡空化数值模拟。表4 列出了计算结果的推力系数KT和扭矩系数KQ，与试验结果基本符合。由于自适应网格加密只针对梢涡区域的网格进行加密，对螺旋桨推力和扭矩计算结果的影响非常有限。

表4 G5计算结果的推力系数和扭矩系数Tab.4 Thrust coefficient and torque coefficient of the calculation results with G5

计算结果的梢涡空化形态由汽体体积分数αv的等值面表示，不同的αv取值对梢涡空化形态的影响较大。在此比较了两套网格计算结果中的αv=0.2、0.5 和0.8 三个等值面，如图9 所示。在三种气体体积分数的取值下，自适应加密后的网格的计算结果均优于未加密的网格，并且随着气体体积分数的增大，梢涡空化形态差距也越来越大。结果表明自适应网格加密方法通过准确加密梢涡空化区域的网格，提高了梢涡空化数值模拟精度。

将G5计算结果的αv=0.2的等值面与试验进行比较，如图10所示。数值模拟的梢涡空化形态与试验非常吻合，可以模拟出梢涡空化较为明显的卷起过程和梢涡空化规律性收缩的节点（红色虚线框），说明本文建立的网格策略能够较好地模拟梢涡空化。但是梢涡空化仅仅在轴向上发展至x=0.16D，在长度上与试验差距较大。为了得到更符合试验的计算结果，需要进一步提高网格分辨率或采取更合适的空化模型和湍流模拟方法。此外，数值模拟过高预报了导边附近的片空化区域。

3 结 语

本文结合网格局部加密和自适应网格加密方法，对E779A 螺旋桨进行了螺旋桨梢涡空化数值模拟。数值模拟结果显示，在近尾流场内的梢涡空化形态与试验非常接近，成功模拟出梢涡空化的卷起现象和规律性收缩的节点。若需提高更远区域的梢涡数值模拟精度，则应进一步提高网格分辨率或者采用更合适的湍流处理方法和空化模型。研究建立了基于OpenFOAM 的螺旋桨梢涡空化数值模拟的网格策略，为将来提升螺旋桨梢涡空化预报准确性奠定了基础。