基于网格自适应技术的PPTC桨梢涡空泡数值模拟

2022-04-21陈影，余龙

船舶力学 2022年4期

陈影，余龙

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

0 引言

螺旋桨空化问题一直是螺旋桨水动力研究中的热点。目前研究螺旋桨空泡的手段主要有实验测量和数值模拟两种。由于空化模拟测试成本较高，大部分的研究都是基于数值模拟方法。螺旋桨的空泡主要分为片空泡、梢涡空泡和毂涡空泡等，目前许多国内外学者已经解决了桨叶表面片空泡的数值预测，但是关于梢涡空泡和毂涡空泡的研究尚不充分。

本文采用的螺旋桨为德国波茨坦水池（SVA）的模型桨VP1304，该桨作为基准螺旋桨于第二届和第四届国际船舶推进器专题研讨会（smp’11和smp’15）被数十家单位计算，有大量的试验数据和数值模拟结果[1-2]。Morgut[3]使用Zwart、FCM和Kunz的三个修正之后的传输模型对PPTC螺旋桨进行了CFD模拟，包括均匀流条件和倾斜轴条件；Gaggero[4]通过RANS方法基于OpenFOAM 开源平台计算了PPTC螺旋桨在斜流中的无空化和空化非稳态性能；陈凯杰[5]基于OpenFOAM 开源平台对PPTC 螺旋桨在全湿流和空化流条件下采用RANS和DES两种方法进行了数值模拟。尽管这些文献成功地验证了PPTC螺旋桨的性能参数以及桨叶上片空泡的分布特性，但是无法模拟出梢涡空泡，尤其是脱出至尾流区域的梢涡空泡的模拟。

由于螺旋桨梢涡空泡尺度非常小，故梢涡空泡的数值模拟对网格分辨率要求很高，目前主流的梢涡空泡计算方法是根据初步计算确定梢涡空泡的大致范围，然后确定控制体的螺距和半径，接下来加密控制体内的网格单元。刘登成[6]采用方块截面的控制体，胡健[7]和Yilmaz[8]采用的是螺旋管几何体加密。虽然控制体加密的方法可以适当延长梢涡空泡，但是捕获到的梢涡空泡均呈现出断裂的现象，有一定的缺陷。故已有研究者试图将自适应网格技术应用于梢涡空泡的模拟，以期得到更好的效果。

应用网格自适应技术的关键在于加密标准和加密区域单元尺寸的设置。加密标准直接影响到网格自适应的区域，应尽量保证加密区域与需要提高物理分辨率的区域重合，而在物理量变化较为平缓的区域保持不变或者适当地粗化网格，实现精准加密和提高计算精度，同时控制计算成本。Yilmaz[9]在模拟INSEAN E779AA 空化流时采用绝对压强作为加密标准，当绝对压力小于10 000 Pa（比饱和蒸气压高的压力值）时，进行网格单元加密，取得了较好的效果，但网格增加了近4 倍。Eskilsson 等[10]将自适应网格技术应用于NACA 0015二维水翼的空化流模拟，并使用了4种不同的加密准则：速度u、涡量Q、体积分数α和压强P，结果表明使用体积分数α最能精准地加密空化区域，而使用Q准则不仅能加密空泡区域，还加密了前缘和后缘区域。关于单元尺寸的控制，Kuiper（1981）[11]对V螺旋桨在J=0.3,0.4,0.5时梢涡空化的测量与研究中总结了空泡数和空泡半径之间的经验关系，发现空泡初生时，每个气泡的最小半径始终约为0.25 mm。故加密后的单元尺寸应该不高于0.25 mm，否则不能较好地捕获到梢涡空泡。

本文采用了自适应网格技术研究梢涡空泡，使用Schnerr-Sauer空泡模型和SSTK-Ω湍流模型，采用体积分数α作为判据，提出一种高效率的网格划分策略，实现空化区域的精准加密，选取了PPTC 桨作为计算对象，更精细地模拟了梢涡空泡。

1 数学模型与研究方法

1.1 控制方程

在数值模拟中，将空化流处理为包含液体和蒸气的两相流，采用修正之后的RANS 方程来求解。修正之后的连续性方程和动量方程如下：

式中：下标m表示多相流，ρm为多相流的密度，其余符号与常规RANS方程一致。

1.2 湍流模型

式中的经验封闭系数和其余辅助方程可见文献[12]，此处不再赘述。式（4）中的F1表示混合函数，其在近壁处边界层上等于1，在远离壁面处等于0，将K-Ω模型和K-ε模型结合在一起。

1.3 空泡模型

上述方程中的ρm和μm由气相体积分数α确定，定义如下：

同样地，需要补充一个关于α的输运方程以使方程组封闭，即空泡模型，本文采用Schnerr-Sauer空泡模型[13]。Schnerr-Sauer 空泡模型在Rayleigh-Plesset 空泡模型基础上，忽略了气泡生长加速效应、粘性效应和表面张力效应。关于α的输运方程为

1.4 自适应网格技术

自适应网格技术（Adaptive Mesh Refinement,AMR）是指基于计算结果根据自适应网格标准在计算过程中不断细化网格单元，从而为某些物理值变化特别剧烈的区域提供足够高的网格分辨率，在提高计算精度的同时又保证计算效率。对六面体网格来说，一个加密级别意味着1个父单元细化成8个子单元。

本文采用体积分数α作为控制标准，当α=0或1时，单元位于气相内或者液相内，不进行加密；当α处于0和1之间时，相界面穿过该单元，进行加密。这样就可以有效地自动加密空泡区域，在提高精度的同时又不会使网格数量过大。

2 研究对象与计算域设置

计算中采用SVA 提供的PPTC 螺旋桨几何模型，该桨螺距可调，故在叶片和桨毂之间有很小的缝隙，模拟中予以忽略，其试验模型和几何参数分别见图1 和表1。本文对研讨会smp’11 上发布的空化案例Case2.3[14-15]进行数值模拟，以校验无空化流场及空化流场的计算方法。

表1 PPTC模型桨的几何参数Tab.1 Main particulars of PPTC propeller

图1 PPTC桨试验模型Fig.1 Model of PPTC propeller

关于螺旋桨的进速系数、推力系数、扭矩系数、敞水效率及空泡数的定义为

3 无空化流场计算与网格无关性分析

3.1 计算域设置

计算域（图2）分为旋转域和静止域，采用速度入口和压力出口条件，静止域采用了smp’11[1]提供的空化流场试验段，旋转域直径为1.2D。无空化流的计算采用多重参考系方法计算其稳态性能，空化计算采用滑移网格法计算非定常性能。

图2 计算域及边界条件设置（蓝色区域为旋转域）Fig.2 Computational domain and boundary condition settings for PPTC

3.2 网格划分

本研究采用切割体网格划分技术来划分网格。在划分网格时，整个计算域网格单元基准尺寸设为100 mm，对旋转区域网格加密，加密区尺寸为基准值的8%。对叶片边缘线网格以及叶片面网格进行局部加密，叶片边缘附近最小网格尺寸为基准值的0.25%，桨叶表面最小网格尺寸为基准值的1%。螺旋桨叶片采用5层棱柱层，棱柱层总厚度为基准值的0.2%。网格（图3）的数量为114万。观察计算完成后的Y+分布直方图（图3）可知，大部分网格单元的Y+小于1，基本上所有网格单元的Y+均小于5，这说明近壁单元基本上都位于粘性子层内。为了验证网格无关性，本文设置了三种网格密度：64 万、114万和250万。

图3 114万计算网格示意图与Y+分布直方图Fig.3 Computational mesh on the PPTC propeller and wall Y+distribution histogram

3.3 网格无关性分析

首先通过对三个进速系数（0.6928，1.2621，1.4944）下的无空化流场计算来验证网格无关性，转速均为15 r/s。表2 呈现了粗糙、中等、精细网格在三个进速系数下的KT和10KQ数值模拟结果以及网格增加时KT和10KQ的变化率。从该表中可以看出，变化率最高值为3.54%；而且随着网格单元数量的增加，大部分计算结果呈现出降低的趋势，只有J=0.6928 情况下KT的变化率有少许增加，但该变化率很小，未超过1%。总体而言，采用114 万网格已经满足计算结果对网格的无关性，故将此网格作为后文无空化流和空化流模拟的计算基准网格（G1）。

表2 推力系数和扭矩系数的网格无关性分析Tab.2 Grid independency analysis for the thrust coefficient and the torque coefficient

3.4 水动力性能参数分析

计算条件及计算结果如表3 所示，采用基准网格G1 进行其余工况的计算，表中ε表示计算结果和试验结果之间的相对误差。大部分工况下KT和10KQ的误差小于1%，所有工况的误差小于3%。图4表明，计算值与试验值吻合，说明了数值模拟的可靠性，可为进一步的空化模拟提供参考。

图4 PPTC螺旋桨敞水特征曲线Fig.4 Propeller performance diagrams

表3 PPTC螺旋桨敞水性能计算结果及与试验结果对比（CFD：数值模拟；EFD：模型试验）Tab.3 Comparison of the open water propeller performance coefficients between measured and computed

4 空化流场模拟

4.1 计算参数设置

空化流场模拟时依赖初始条件的设置，故在进行空化流场模拟之前，先计算该空化流场对应的敞水情况，计算收敛后导出速度场和压力场数据，然后将其导入至空化模拟中作为初始场，使计算快速收敛。将模型改为隐式非定常条件，时间步长设为一步旋转1°的时间，时间离散格式为二阶，计算参数与试验[14]（见表4）保持一致，采用基准网格（G1）进行空化流场的模拟。

表4 PPTC桨空化流场模拟参数设置Tab.4 Experimental and computational conditions for cavitation flow

4.2 空化流场与敞水流场推力系数对比

表5 为三个案例中CFD 结果与EFD 结果之间推力系数的对比。整体而言，数值模拟的推力系数和试验值较为接近，略低于试验值；Case2.3.1 的无空化流和空化流的模拟效果最好，误差较小；在Case2.3.2中，空化流的误差偏高，达到5%左右；在Case2.3.3中，无空化流的误差略大。

此外，从表5 还可以看出无空化流和空化流之间推力系数的差异，对每一个案例，空化的存在都会使螺旋桨的推力系数减小。特别是在Case2.3.2 和Case2.3.3 中，这一推力损失现象很突出，损失值达到了14%以上。

表5 无空化流和有空化流推力系数对比Tab.5 Comparison of the thrust coefficients between cavitation flow and no-cavitation flow

4.3 空泡分布情况

由图5 可见，总体而言，三个案例的空泡分布区域与试验较为接近，但细节上有所差别。在Case2.3.1 中，r＞0.95R范围内的片空泡与试验保持一致，接近桨毂的叶根区域空泡分布也与试验很接近；但在导边附近的0.5R至0.9R范围内，数值模拟结果出现了试验结果没有的片空泡。在Case2.3.2中，数值模拟在导边附近捕获到部分空泡，与试验不符；叶根区域的片空泡分布范围略大于试验。在Case2.3.3 中，空泡从吸力面转移至压力面上，数值模拟没有捕捉到叶根区域的片空泡，且压力面上导边附近的空泡范围略小于试验。另外，三个案例的数值模拟均未捕捉到梢涡空泡，这是因为网格分辨率不足造成的。

图5 空泡分布形态对比（数值模拟中用α=0.2等值面表征空泡）Fig.5 Comparison between the computed cavitation and the experimental cavitation

5 梢涡空泡

5.1 网格自适应加密

为了捕获尾流中的梢涡空泡，需提高该处的网格分辨率。本文采用三套网格（图6）来模拟梢涡空泡。第一套网格就是前文所述的基准网格（G1），网格单元数量为114万。第二套网格（G2）在G1基础上应用螺旋管几何加密，螺旋管几何如图6 所示，单元数量为415 万；其中螺旋管的直径为10 mm，其螺距和半径减少量从G1 计算结果中获得，加密基本尺寸为0.5 mm。第三套网格（G3）是以G2 作为初始网格在计算中应用本文采用的网格自适应技术后生成的最终网格，数量为638万；当气体体积分数α在0至1之间时自动加密单元，加密级别取1，加密后单元基本尺寸变为0.25 mm，网格优化贯穿整个计算过程，每5个时间步应用一次网格自适应。若不应用网格自适应方法，直接将螺旋管内单元尺寸设为0.25 mm，则网格量高达2192 万，故本文中的网格划分策略可显著减少网格单元数量，提高计算效率。文中计算梢涡空泡的具体流程图如图7所示。

图6 螺旋管及计算网格Fig.6 Spiral tube geometry and three kinds of computational mesh

图7 梢涡空泡计算流程图Fig.7 Flow chart of tip vortex cavitation simulation

5.2 推力系数对比

表6 为三种网格划分策略的推力系数计算结果的对比，其中KTG1表示采用网格G1 计算得到的推力系数，KTG2、KTG3同理。三套网格计算的结果仅有0.6%以内的微小差异，几乎可忽略。这三套网格之间最大的不同在于对梢涡脱出区域的处理不同，对推力系数的影响很小，在正常的波动范围之内。

表6 三种网格划分策略推力系数对比Tab.6 Comparison between CFD and EFD for three mesh methods

5.3 空泡分布形态对比

图8展示了三种网格划分策略的空泡分布形态，桨叶上的片空泡基本没有变化，但采用螺旋管与网格自适应技术结合的网格划分策略可以明显改善梢涡空泡和毂涡空泡的模拟，特别是尾流区域中的梢涡空泡延长效果显著。

图8 梢涡空泡的延长效果对比图Fig.8 Improvement of tip vortex cavitation extension

5.4 梢涡空泡细节对比

由图9（a）可见，梢涡空泡呈现出卷起的现象，这一现象是由于涡流强度降低或压力增加，导致片空泡和梢涡空泡产生相互作用而引起的，卷起时产生一定规律的节点，节点处的空化涡管直径减小。仅采用螺旋管加密的G2网格模拟出的梢涡（图8）在节点处由于网格分辨率不足，呈现出断裂的现象；但是观察图9（b）可知，应用网格自适应技术后成功地模拟了涡管卷起现象，更精确地获得了流场的空化模式。

图9 Case2.3.1梢涡空泡放大对比图Fig.9 Comparison of the tip vortex cavitation roll-up phenomena between EFD and CFD

5.5 漩涡结构图

采用Q准则来表达漩涡，当Q＞0时表示旋转在流动中占主要地位，其值越大表示涡的强度越大。图10第二列展示了G3计算结果，用Q为10 000的等值面来表示漩涡。与G1计算结果（图10第一列）相比，捕获的梢涡长度要长得多，而且呈现出更为光滑的管状结构，涡管的直径也更加均匀，但是到后面未加密区域时，涡管直径变大，并且很快耗散，若想延长捕获的涡管，可通过增加螺旋管的长度来实现。

图10 梢涡与梢涡空泡对比图Fig.10 Comparison between tip vortex and tip vortex cavitation

观察Q=10 000的等值面与α=0.2的梢涡空泡等值面（图10第三列）可以发现，两者形态十分接近，但是可以观察到Q的等值面范围远远大于α的范围，特别是桨叶表面上，此处可证明采用Q准则作为网格自适应的判断标准时会进行许多不必要的单元加密。若仅想探究螺旋桨的梢涡空泡，从加密精准度来看，气体体积分数α是一个更好的选择。

图11展示了两个案例x=0.3R截面处的涡量分布图，左半边为采用G1计算的结果，右半边为采用G3 计算的结果。可以看到，两种网格的涡量分布图均为周期性分布，周期为5，对应于螺旋桨叶数。但是采用G3网格计算可以获得更精细的涡量分布，由于G1网格在梢涡处分辨率不够高，所以捕捉不到涡核处的最大涡量，而应用螺旋管精细加密和网格自适应的梢涡则显示出更大的涡强。此外，结合图12 可以发现涡量最大值并非出现在梢涡的中心，在向梢涡中心靠近时，涡量强度呈现出先增加后减少的变化趋势。这是由于空泡的作用，当空化流中的压力低于饱和蒸气压时，就会产生空泡，使压力不会再降低，因此在空泡内部压力相对比较均匀，故涡量最大值并非出现在中心。