水翼三维非稳态空化流动特性数值研究

邓 飞， 张 松， 梁 兴， 李伟平

（1.南昌工程学院南昌市低品位能源转化与系统节能技术重点实验室，南昌 330099；2.渝能（集团）有限责任公司，重庆 401121）

0 引言

水翼空化具有强烈的非定常特性，易引起流场空泡团产生，形成周期性脉动造成水翼水动力载荷变化，从而出现水翼振动和流场噪声等问题［1］，影响流体机械的安全可靠运行。一直以来，研究空化发生的机理、特性和抑制技术都是流体机械研究的重点和难点［2］。研究表明，水翼安放角度［3］、雷诺数、绕流速度［4］等水力因素决定了水翼的绕流特征。谢春梅等［5］针对高雷诺数下三维水翼的梢涡流场，采用基于浸没边界法的壁模型大涡模拟方法，对NACA16020 水翼进行了数值研究，系统观察了流场的各类统计特性。Wei 等［6］对翼型在6°攻角下的水动力及空化性能进行数值研究，当空化发生时水动力特性将剧烈恶化。刘胡涛等［7］对高雷诺数的NACA0012 水翼绕流开展了数值计算，探讨了初始攻角条件、绕流速度以及刚心位置三个因素对水翼振动的影响，水翼振动状态依赖于系统的初始条件，当绕流速度不断增大时水翼将发生颤振现象。田北晨等［8］基于建立的多尺度算法对绕NACA66 水翼空化流动进行模拟，在回射流发展阶段，离散空泡逐渐增加并分布在回射流扰动区。Yang等［9］采用MRT-LBM方法对空泡溃灭进行热力学分析并验证了MRT-LBM空泡溃灭的有效性。

目前，针对水翼空化现象的机理研究主要集中在空泡初生阶段［10-11］，对不同绕流速度下的水翼空化特性及其准周期性演变过程研究较少。本文针对不同绕流速度下的NACA6408 水翼非稳态绕流过程，数值研究各流动因素的变化规律以及其对水翼空化过程的影响，旨在获取水翼发生空化的临界绕流速度和压力等重要参数，以及空化过程的准周期性特征等。

1 数值模型与网格划分

本文采用计算流体动力学软件Fluent 16.0 进行数值模拟研究，如图1 所示为NACA6408 三维水翼计算域模型，依据该水翼的相对弯度和相对厚度，对其及流体域进行建模。整个流体域长度为10b，其中b为水翼弦长，进口距水翼前缘距离为3b，出口距水翼尾缘为6b，计算域的进口面、出口面尺寸为150 mm×150 mm。

图1 NACA6408三维水翼计算域模型

本文利用ICEM 网格划分软件对三维NACA6408水翼计算域进行结构化网格划分，针对水翼的头部、尾部及上下壁面等部位进行局部网格加密。经过网格无关性计算验证，当网格尺寸小于0.07 mm 时，水翼的升力系数趋于稳定［12］，为合理利用计算机资源，提高数值模拟的准确性、经济性，采用边界层为0.07 mm的网格尺寸，网格总数为5501430。采用k-ε 湍流模型［13］和Zwart-Gerber-Belarnri（ZGB）空化模型［14］进行模拟计算，与Schnerr-Sauer（SS）空化模型比较，本文采用的Zwart-Gerber-Belarnri（ZGB）空化模型使用空泡密度表示气相的体积分数［15］，能够更好地模拟出非定常空化的准周期性，捕捉到更强烈的云状空化。

2 计算结果与讨论

计算模型的b＝150 mm，计算域长度则为1.5 m，在绕流速度υin＝1、5、10、15 和20 m／s 5 种工况下分别进行数值模拟计算。在沿翼型主流方向等间距截取11 个截面i（i＝1，2，…，11）。由于翼型结构特点，同一个截面将被翼型分成上下两部分（见图1），令水翼吸力面以上的截面为B 部分，即Bi，压力面以下的截面为A部分，即Ai。

通过对绕流攻角为4°时的5 种不同绕流速度工况依次计算，数值结果表明翼型在υin＝1、5 和10 m／s工况下不会发生空化，在υin＝15 和20 m／s 时产生空化现象。在前3 种工况下，截面A 部分的压力在不同位置处变化较小，主要是由于NACA6408 翼型结构的压力面曲率很小；截面B部分的平均压力pav先减小后增加，即pav从截面B1减小至截面B4后增大，直至截面B11，原因是截面面积的变化使速度呈现出先增大后减小的规律，pav则是先减小后增大。此时同一截面的A与B部分的pav和湍动能κ 变化平稳，不存在空化现象。

如表1 所列为不同υin工况下以截面1、4、6、8、11为代表的B部分截面pav值。从表中可以看出，在同一υin下，不同时间（t＝0.05 和0.1 s）与不同截面间的pav略有不同，但是维持在同一个数量级。υin＜10 m／s时没有发生空化，截面B1～B4间的pav随着速度的增加而逐渐减小，相同速度下截面B6～B11间的pav逐渐增大；υin＝5 m／s时，pav≥89.528 kPa；当υin＝10 m／s时，pav≥52. 192 kPa；而当υin＝15 m／s 时，pav≥11.890 kPa，当υin＝20 m／s 时，pav≥9.763 kPa。因此发生空化时，pav≤52.192 kPa。

表1 不同工况下截面B部分的pav变化kPa

如图2 所示为υin＝15 m／s下翼展中间截面体积分数及压力云图。由图2（a）可知，在此工况下翼型的空化呈现一定的周期性，且周期T＝75 ms。由于此工况为非定常流动，无法对每个时刻进行描述，故将一个周期均分为8 等分以方便分析。在T／8 时，空化开始出现于沿翼型吸力面主流方向b／4 ～b／2 处且附着于吸力面，随着时间推移空化过程越来越明显；在2T／8时，翼型吸力面b／4 ～b／2 处附着型空泡开始向钩状空泡转变；在3T／8 时，附着型空泡长度达到最大，其在钩状空泡向后延伸至翼型吸力面4b／5 处开始脱落；在4T／8 ～5T／8 时，空泡在翼型吸力面4b／5 处部分脱落，脱落后的空泡流向下游；在6T／8 时，剩余空泡继续附着于翼型吸力面b／3 ～5b／6 处且空泡长度随着时间推移而逐渐变小；在7T／8 时，空泡长度缩短至翼型吸力面b／3 ～2b／3 处，最终在8T／8 时逐渐减小为零并进入下一个周期，空化过程呈现出准周期性规律。由图2（b）中可知，在6T／8 时，低压区的覆盖范围最大，随后低压区的覆盖范围开始减小；在8T／8 时，低压区的覆盖范围达到最小，且与初始时刻的压力p相同，从而形成一个完整的周期。结合图2（a）可知，此时翼型吸力面低压区的p低于空化发生的临界压力。

图2 在υin ＝15 m／s时，翼展中间截面的体积分数及压力云图

2.1 水翼各截面压力p变化分析

针对υin＝15 m／s工况水翼表面发生空化时A、B两部分各截面处的压力p、湍动能κ、速度分量u、v和w等参数变化进行分析。图3 所示υin＝15 m／s时，翼型主流方向i位置的A和B 部分在一个周期内最大、最小和平均压力折线图。由图可知，A 部分截面的最大压力pmax和最小压力pmin随着截面呈现逐渐减小的趋势，即在A2～A4基本趋于平稳；A5～A11呈现逐渐下降。结合NACA6408 翼型结构以及在υin＝1、5 和10 m／s工况下无空化现象，可以判定A1～A4压力下降不是截面面积变化造成的，而是受空化发展的影响。B部分截面的pmax和pmin在B1～B4间呈断崖式下降，原因是翼型结构的影响。B4～B11呈现上升趋势，上升过程波动幅度存在较大的差异。

图3 在υin ＝15 m／s时，11个不同截面位置的压力变化折线图

在0 ～2T／8 段，pav呈现大幅度上升，整体平均增加了12 kPa 左右，原因是因为空泡附着于翼型吸力面，且空泡占比很小，对流场影响较小，从而使pav逐渐上升；在3T／8 ～6T／8 段，B4～B11平均压力pav在降至最低后几乎不再变化，中间位置的截面压力维持在110 kPa左右，这是由于此时附着型空泡开始向钩状空泡转变，并且钩状空泡脱离翼型吸力面形成独立的空泡，此时为空化的核心阶段，压力整体处于低压状态；在7T／8 ～8T／8 时刻B1～B11的pav开始逐渐升高，并恢复至1T／8 的压力水平，这是由于脱落的空泡随主流流向下游并逐渐溃灭的原因。

综合分析图3 可知，A 部分截面的pmax、pmin和pav的波动幅值都明显大于B 部分，因此水翼翼型结构对压力脉动的幅值有着重要影响，空化显著的区域其压力脉动更加明显，可以利用压力脉动变化来探究空化现象特性。

2.2 各截面κ分布特性

图4 所示υin＝15 m／s下A、B部分各截面κ折线图。由图可知，在一个周期内，A部分截面的κ几乎不存在变化，。此时流场A 部分截面相对稳定，可以排除翼型结构对κ的影响。B部分截面中，在B1～B5区段内湍动能随时间变化较小，流场相对稳定；而在B6～B11内湍动能变化较为明显，在B11处湍动能最大值增加至3.6 J／kg，说明B6～B11内部空泡的产生使得流场紊乱。空化发生区域在0 时刻湍动能变化明显，为空化初始阶段，此时流场相对紊乱，能量耗散大；在2T／8 和3T／8 时产生附着型空泡，流场相对稳定；4T／8至8T／8 时空化由附着型空泡向钩状空泡转化，在脱离翼型吸力面后最终溃灭，流场相对紊乱，能量耗散大，如图4 的B部分截面所示。

图4 在υin ＝15 m／s时，11个不同截面位置的κ变化折线图

因此，可看出湍动能代表流场的紊乱程度，空化的发生与湍动能大小有关，同时与翼型结构和回射流等有关。湍动能又随速度增加而增大，因此流场速度也是判断流场紊乱的因素之一。

2.3 水翼各截面速度分布

为了分析各截面的速度对流场脉动特性的影响规律，如图5 所示υin＝15 m／s 下A、B 部分截面在同1个周期内X、Y、Z3 个方向的平均分速度u、v和w折线图。

图5 在υin ＝15 m／s时，X、Y、Z方向11个不同截面位置的平均分速度变化

（1）A部分截面。X方向上A1～A3的主流方向速度u基本平稳；A4～A11间速度u有小幅度上升，原因是翼型截面积变化造成的影响。Y方向上速度v在A1～A4速度为正且变化缓慢，而在A5～A11间速度v方向为负且逐渐增大，原因是水流以较大的υin冲击翼型头部结构从而产生反方向的速度。Z方向的速度w基本维持在0 附近，整体比较稳定。

（2）B部分截面。X方向上B1～B4区间主流方向速度u逐渐增大；在B5～B11区间u逐渐减小且变化幅度较大，原因是翼型截面积变化以及空化现象的产生。Y方向上在B1～B4速度v为正且逐渐减小；在B5～B11速度v方向为负，原因是空化的产生形成复杂旋涡；随着时间推移，Y方向上速度零点的出现时间沿截面逐渐后移，说明产生的旋涡逐渐往后移动，证明了空泡的脱离并且向下游移动，B9～B11速度v方向为负且逐渐增大，主要是由于空泡随主流流向下游或者发生溃灭，导致Y方向速度v最后逐渐减小。Z方向上速度w维持在0 附近且有正有负，说明空化过程产生的旋涡并不是平面型，而是复杂的立体型旋涡。

综上所述，X主流方向速度u主要受υin的影响，Y和Z方向的速度v和w主要受空化过程的影响，速度方向发生改变说明有空泡的脱离、溃灭和立体型旋涡的产生。针对水翼空化过程的X、Y、Z3 个方向分速度分析是探讨空化现象特性的重要研究思路之一。

3 结语

本文基于计算流体动力学技术数值研究了NACA6408 水翼在绕流攻角4°时不同绕流速度工况的流场特性，通过数值分析υin＝15 m／s下水翼空化发生过程，获得流场p、κ以及各方向分速度等参数的变化规律。结果表明：①随着进口速度的增加，各截面上的平均湍动能以及速度的波动幅值明显增加，压力、湍动能以及流场速度的变化与空化过程密切相关。②水翼空化发生区域压力、湍动能以及流场速度变化明显，水翼压力面以下区域的参数变化相对缓慢，水翼吸力面以上区域变化相对明显。压力面变化曲率较大导致进口水流流经翼型吸力面时流场相对紊乱。③在υin＝1、5 和10 m／s工况时水翼表面不会产生空化现象，在υin＝15 和20 m／s时发生空化，同时空化基本在B5～B11发生。④空化发生时非主流方向（Y和Z方向）的速度存在方向变换，形成了复杂的立体型旋涡。本文的数值计算结果对深入认识水翼空化演变机理，提高水翼抗空化性能具有重要的参考意义。