绕水翼非定常空化流动的涡动力学分析

赵 宇，王国玉，黄 彪，胡常莉，陈广豪，吴 钦（北京理工大学 机械与车辆工程学院，北京 100081）

绕水翼非定常空化流动的涡动力学分析

赵 宇，王国玉，黄 彪，胡常莉，陈广豪，吴 钦
（北京理工大学 机械与车辆工程学院，北京 100081）

文章采用数值方法对绕水翼的非定常空化流动进行了流动计算，并采用涡动力学方法对非定常空化流场特性进行了分析。计算中，采用基于当地涡旋修正的湍流模型封闭控制方程，通过实验对该计算结果进行验证；采用涡量法分析了空化的发生和发展对涡量输运过程的影响。结果表明：基于当地涡旋修正的湍流模型可以准确地模拟水翼周围非定常空化流场结构；水翼空化流场具有强烈的非定常特性，空化的发展和流场中的涡结构有着紧密的联系；反向射流的作用会导致附着空穴尾缘附近速度梯度的非定常变化，是导致涡量产生和消失的重要因素；空化区域内部水气两相的质量交换会导致流场内体积变化率的变化，使得空化区域内部流体的非正压性显著增强，上述两者导致空化区域内部的涡量分布呈现强烈的非定常特性；空化现象引起的水气两相之间的相互转换以及对湍流的作用都是影响涡结构空间分布的重要因素。

空化；湍流模型；非定常；涡动力学分析

0 引 言

作为高速水动力学研究领域的核心问题，空化是指在一定条件下，液相介质中发生的一种相变现象。在水流中，当其中的压力降低到水的饱和蒸汽压时，水流中将产生充满气体和蒸汽的空穴，这时，水流中发生了空化。这些空穴是一些单个的球形空泡，多个空泡形成的空泡群或外形不规则的空腔。空化是包含相变、非定常、多维湍流等多种复杂影响因素的流动现象，由空化非定常特性引起的一定尺度的空泡团的脉动、脱落是引起导弹出水过程水动力载荷、螺旋桨空泡噪音、振动和水力机械空蚀的主要原因[1-3]。涡运动是流体中普遍存在的一种运动形式，是湍流的一种基本结构[4]。涡运动因其固有的非定常特性和非线性机制，又因其对自然界的探索和工程应用都具有巨大的意义，一直是流体力学的研究前沿[5]。空化现象的非定常发展过程与涡运动息息相关[6]，空泡的产生和发展会产生复杂的涡结构[7]，而涡的形成可强化流场的湍流强度，使得湍流和空泡之间的交互作用变得复杂[8]。

近年来，国内外学者开始关注空化流场涡结构的演变过程与空化发展的相互作用，Gopalan和Katz[9]采用粒子测速系统和高速录像系统对喷口结构的片状空化流场进行研究，实验结果表明涡的产生是空化流动的重要特性。Wang等[10]采用实验方法研究了水翼周围的附着空化湍流场的动态特性，结果表明云状空化下大尺度涡结构和反向射流存在密切的关系。张博等[11]采用数值与实验的方法研究了水翼周围的云状空化流动，发现空穴尾部大尺度的涡运动引起了近壁处的逆压梯度，从而在壁面和附着空穴之间产生了液相的反向射流，反向射流逐渐向水翼前缘运动，从而导致了附着空穴的断裂。Huang等[12]采用实验和数值相结合的方法对水翼非定常云状空化流动进行分析，发现涡结构的脱落是导致云状空化流场强烈非定常性的主要原因。由此可见，非定常空化现象与流场中的旋结构紧密相关，了解分析空化发展过程与流场涡的产生、演化的相互作用是更加深入理解非定常空化流动机理的重要途径。涡动力分析方法是研究空化流场涡结构的重要方法，涡量输运方程是分析空化和涡结构相互作用的有力工具。Nagendra等[13]采用大涡模拟的方法对文丘里管内部典型空化流场进行计算，并采用涡量输运方程对文丘里管空化流动进行研究，结果表明空化现象的发生降低了涡动力方程中的拉伸项。黄彪[14]对水翼周围云状空化流场进行研究，和无空化的流场结构对比发现，空化现象的产生可以有效地增加涡量输运方程中的斜压矩产生项。空化流动的涡动力特性逐渐成为空化研究的热点方法，但上述研究缺少从涡量输运的角度分析典型空化流动的非定常特性。

本文采用数值方法对绕CLARK-Y水翼的非定常空化流场涡结构进行研究，对水翼周围非定常空化现象与涡结构的相互关系进行分析；并引入涡动力方程，分析空化现象对涡旋结构发展的影响。

1 数值方法

1.1 控制方程

计算采用单流体模型，适用于牛顿流体的笛卡尔坐标系下忽略体积力和热传导守恒形式的N-S方程如下：

式中：ρm是水汽混合相密度，ρl是液相密度，ρv是汽相密度，αv是汽相体积分数，αl是液相体积分数，u是速度，p是压力，μm是混合相的粘性系数，μl和μv分别是液相和汽相的粘性系数，μT是湍流粘性系数。下标（i,j,k）表示笛卡尔坐标系下的坐标轴方向。和-分别表示凝结和蒸发速率。本文中，采用Kubota空化模型[15]：

式中：αnuc是空化核子的体积分数，RB是空泡直径，pv是饱和蒸汽压，p是当地压力。Cdest和Cprod是凝结和蒸发系数。本文中，上述参数取值如下：αnuc=5×10-4，RB=1×10-6m，Cdest=50，Cprod=0.01。

1.2 湍流模型

本文采用修正标准k-ε湍流模型来封闭上述方程组，该模型在标准k-ε湍流模型基础上引入当地涡旋运动系数，考虑了当地涡旋运动对湍流的影响。在标准k-ε湍流模型中，湍动能方程和耗散方程如下：

式中：Pk表示湍动能产生项，湍流粘性μt定义如下：

本文将当地涡旋运动效应修正引入到标准k-ε湍流模型中。在k方程中的湍动能产生项上乘以修正系数Pk→Pk·fr。修正因子的表达式如下[16]：

图1 计算模型及边界条件设置Fig.1 Geometry and Boundary Conditions

1.3 几何及边界条件设置

计算采用上述计算模型，流动介质参数与25℃的水和水蒸气保持一致。水相密度和运动粘度分别设置为ρl=999.19 kg/m3和μl=ρlvl=1.139×10-3Pa·s；汽相设置为ρv=0.023 08 kg/m3和μv=9.862 6×10-6Pa·s。汽化压力设置为pv=3 169 Pa。几何模型和边界条件如图1所示，CLARK-Y水翼弦长0.07 m，最大相对厚度为11.7%，翼展0.07 m。水翼攻角设置为8°，入口处设置为法向速度入口，速度大小为10 m/s，从而保证雷诺数Re=U∞c/vl=7×105。出口处采用压力边界条件，设定环境压力为43 169 Pa，从而保证空化数为（σ=（p-p∞）/（0.5ρv2））0.8。水翼表面和流场顶部和底部均设置为无滑移壁面边界条件，前后边界设置为对称面。在计算中，时间步长Δt=1*10-4s，保证计算平均科朗特数CFL=U∞Δt/Δx=1。

1.4 模型验证

为了验证数值方法的可行性，将数值模拟结果和实验结果进行对比。实验在一闭式空化水洞[9]内进行。该空化水洞试验段尺寸为0.70×0.07×0.19 m，流速范围为2～15 m/s，可达到的最小空化数为0.3。在实验中，采用高速摄像系统来观察空泡形态的演变历程，采用流体动力测量系统测量水翼受到的升力。

从高速摄像系统得到的系列空穴形态图中可以发现水翼吸力面空泡表现为准周期性的发展、断裂过程，呈现出显著的非定常特性。根据空穴形态选取特征时刻，图2给出了各特征时刻的实验所得空穴形态图和数值计算结果的对比情况。从图中可以看出，数值方法预测的空化发展过程与实验结果相吻合，能准确地捕捉空泡的初生、发展以及断裂溃灭过程。

图2 特征时刻空穴形态对比图Fig.2 Instantaneous outline of cavities of experiment and simulation

定义升力系数为：

式中：U∞为距实验段上游入口210 mm处参考断面上的平均速度，A为水翼的有效升力面积，Fy为水翼所受到的阻力和升力。为了研究数值模型对预测翼型空化水动力特性的影响，图3给出了计算得到的翼型升力随时间的波动情况，并与实验结果进行对比。从图中可以看出，实验和数值计算得到的水翼升力系数均呈现准周期性波动情况且波动的幅值基本一致，实验测得的平均升力系数约为0.74，数值计算得到的平均升力系数为0.71，二者基本吻合。

图3 翼型升力系数曲线对比Fig.3 Lift coefficient of experiment and simulation

2 涡量分析方法

为了分析水翼周围的流场结构，定义特征涡量如下：

式中：δ≈0.1c=0.007 m为实验测得的无空化条件下水翼尾缘的湍流边界层厚度。

为了更清晰地反映流场内的涡旋结构，引入Q判据[17]：

Q是在伽利略变换下的速度梯度张量第二不变量，式中，Si j表示形变率张量，ωi j表示涡量，分别定义为：，其中i和j代表了不同的方向。对Q＞0的区域定义为涡流所在，因为从上式可看出，Q可说是涡度和轴向的形变率的差值，当差值为正，表示旋转的趋势大于轴向的形变，也就是此区域受涡流所主宰。

引入涡量输运方程：

式中第一项（I）表示由于流场的速度梯度引起的涡线的伸缩和弯曲，从而导致涡量的大小和方向都发生变化，在绕水翼的非定常空化流动中，空化现象引起流场内的速度梯度发生明显的变化，从而导致涡量的绝对值发生变化，并且该速度梯度分布具有强烈的非定常性，从而导致速度梯度产生项具有强烈的非定常性。第二项（II）表示流体微团的体积变化引起的涡量大小变化。对于不可压缩流体，该项绝对值为0，但由于空化现象的发生，使得流场内部空化区域的体积变化率发生变化，从而引起涡量的绝对值发生变化。第三项（III）表示由于不平行的压力梯度和密度的梯度导致的斜压矩对涡量的作用效果。对于正压流体，其密度仅是压力的函数，则密度和压力项具有相同的变化梯度，该项为0，但在空化流场中，压力和密度的梯度并不总是平行的，特别是在空穴界面及闭合区域内，从而会引起涡量的变化。第四项（Ⅳ）表示涡量的粘性扩散效应。定义特征涡量输运项如下：

3 结果分析

图4分别给出了特征时刻下的沿翼展方向的涡量（根据特征涡量进行量纲一化处理）等值面分布情况，以及对应时刻的Q分布情况。综合图2和图4可以看出在t1=0.25T0时刻，水翼前缘位置开始出现少量片状附着空穴，水翼吸力面尾缘附近仍能观察到上个周期脱落的空泡团。水翼前缘的附着空穴逐渐发展，脱落空泡团逐渐脱离水翼表面，向下游运动。从涡量分布图中可以看出，分别从水翼前缘和尾缘开始，产生方向相反的涡量聚集带状区域并逐渐向下游延伸。前缘涡（-）和尾缘涡（+）相互作用，在水翼吸力面附近汇合，并逐渐脱落，汇合区域对应于水翼吸力面的空化区域。

到t2=0.5T0时刻，附着空穴长度增加到最大，在该时刻，空穴边界开始出现不稳定的波动，特别是在附着空穴的末端附近；此时前缘涡带逐渐向下游发展，并较为稳定地附着在水翼吸力面上，附着空穴位于前缘涡带中，尾缘涡控制范围较小。从Q分布中可以看出，附着空穴内部Q绝对值较小，说明在附着空穴内部即前缘涡带核心位置的流动较为稳定，旋转效应和变形效应相对平衡，流动没有较强的涡旋特性。而在附着空穴的末端附近，Q绝对值较大，且正负区域间隔分布，这说明在该区域内部，旋转效应和变形效应处于非平衡状态，两者相互转化，相互制衡；该区域对应于前缘涡带和尾缘涡的交汇处，方向相反的涡结构相互作用，导致该区域内部出现强烈的不稳定性。

图4 特征时刻涡量分布和Q分布图Fig.4 Instantaneous contours of dimensionless vorticity and Q

随着空化现象的发展，水翼尾部附近的指向前缘涡带核心区域的反向压力梯度逐渐发展增大到足以克服由于流动产生的正方向动量，在近壁面区域产生反向射流，在反向射流前端形成局部高压，并推动近壁面的液相流体向水翼前缘运动；从图中可以看出在前缘涡带（-）和水翼壁面之间存在具有一定厚度的正向涡带（+），这是由于前缘涡和尾缘涡的相互作用导致尾缘涡结构的断裂和脱落，一部分带有正向涡量的流体脱离水翼尾缘，另一部分带有正向涡量的流体在反向射流的作用下沿壁面向水翼前缘运动。到t3=0.75T0时刻，反向射流前端逐渐向水翼前缘移动，到达附着空穴前端附近，附着空穴在高速水流的作用下发生断裂，反向射流区域内流体带有正向的涡量（+）。在断裂位置处和水翼尾缘附近，存在方向相反的涡结构相互作用，旋转效应和变形效应相互制约，相互转换。

随着空化的进一步发展，由于惯性的作用，反向射流继续向水翼前缘运动，延缓了水翼前缘附着空穴的迅速增长，另一方面，附着空穴断裂的另一部分在主流的作用下向下游运动，从而形成了大尺度的脱落空泡团。到t4=0.9T0时刻，附着空穴逐渐发展，开始逐渐覆盖水翼吸力面，大尺度的脱落空泡团逐渐脱离水翼吸力面。从涡量云图中可以看出，附着空穴处于前缘涡带的位置处，脱落的空泡内部涡结构较为复杂，具有相反方向涡量的流体存在强烈的相互作用，这是由于具有正向涡量的反向射流冲击具有相反方向的前缘涡带，使之发生断裂，从而形成脱落的涡旋结构，脱落的涡旋结构内部同时保留着反向射流的正向涡量以及前缘涡带的反向涡量，所以导致脱落空泡内部存在较为复杂的涡旋结构。

图5中给出了不同时刻下沿翼展方向的涡量各输运项等值面云图分布情况，图中各输运项根据特征涡量输运项进行量纲一化处理。从图中可以看出，I项发展规律具有强烈的非定常特性，在t1= 0.25T0时刻，绝对值较大的位置对应于在附着空穴尾缘和水翼尾缘附近，在t2=0.5T0时刻，附着空穴逐渐发展到覆盖整个吸力面，绝对值较大区域出现在水翼尾缘，从t3=0.75T0时刻到t4=0.9T0时刻，绝对值较大位置对应于反向射流的前端，从水翼尾缘逐渐向前缘发展。这是由于在附着空穴尾缘，在反向射流的作用下出现速度旋涡，该区域内部的速度梯度分布较为复杂，因而导致该区域内部的I项的绝对值较大。

图5 特征时刻涡量输运项分布图Fig.5 Instantaneous contours of dimensionless vorticity transport term

绝对值较大的II项分布逐渐覆盖整个吸力面，并且具有一定的非定常特性：在t1=0.25T0时刻，在脱落空泡团的核心位置的II项绝对值较大，从t2=0.5T0时刻到t4=0.9T0时刻，附着空穴逐渐发展，附着空穴核心位置对应于II项绝对值较大的区域。这是由于II项表征的是流体微团的可压缩性，在附着空穴和脱落空泡核心位置，水蒸气含量较大，水气之间存在强烈的质量交换，从而导致强烈的可压缩效果。

绝对值较大的III项区域分布规律具有强烈的非定常特性。在t1=0.25T0时刻，绝对值较大的区域主要分布在水翼尾缘附近，水翼前缘出现绝对值较大区域但体积较小，在t2=0.5T0时刻，绝对值较大区域相互交错，覆盖整个水翼吸力面，从t2=0.5T0时刻到t4=0.9T0时刻，水翼吸力面的绝对值较大的III项区域出现断裂，并且断裂位置逐渐向水翼前缘运动。综上所述，绝对值较大的III项区域和空化区域相一致，这是由于空化现象的发生，使得空化区域内部的流场结构处于非平衡状态，水蒸气含量不断变化，并且速度旋涡，流体的非正压性得到显著加剧，从而导致流场内部的密度梯度与压力梯度方向不一致。

绝对值较大的Ⅳ项分布逐渐覆盖整个吸力面，并且随着空化现象的发展呈现出一定的非定常特性。这是由于空化现象的发生发展导致水翼吸力面处于气液混合状态，水气之间存在的相互转换导致该区域内部流体的粘度发生显著的非定常变化，而同时，空化现象的发生也导致水翼吸力面的湍流强度产生强烈的非定常变化，上述两种因素的作用，导致水翼吸力面附近的Ⅳ项取值较大。

4 结 论

本文采用数值方法对绕水翼的非定常空化涡动力流场特性进行研究，分析了水翼周围非定常空化现象与涡结构的相互关系，并引入涡动力方程，分析空化现象对涡结构发展的影响，得到主要结论如下：

（1）水翼空化的发展和流场中的涡结构有着紧密的联系，基于当地涡旋修正的湍流模型可以准确地模拟水翼周围非定常空化流场结构。

（2）反向射流的作用会导致附着空穴尾缘附近速度梯度的非定常变化，是导致涡量产生和消失的重要因素。

（3）空化区域内部水气两相的质量交换会导致流场内体积变化率的变化，空化区域内部流体的非正压性显著增强，上述两者导致空化区域内部的涡量分布呈现强烈的非定常特性。

（4）空化现象引起的水气两相之间的相互转换以及对湍流的作用都是影响涡结构空间分布的重要因素。

[1]刘筠乔,鲁传敬,李 杰,等.导弹垂直发射出筒过程中通气空泡流研究[J].水动力学研究与进展A辑,2007,22(5): 549-554. Liu Yunqiao,Lu Chuanjing,Li Jie,Cao Jiayi.An investigation of ventilated cavitating flow in vertical launching of a missile[J].Journal of Hydrodynamics(Ser.A),2007,22(5):549-554.

[2]Cervone A,Bramanti C,Rapposelli E.Experimental characterization of cavitation instabilities in a two-bladed axial indeucer[J].Journal of Propulsion and Power,2006,22(6):1389-1395.

[3]顾 巍,何友声,胡天群.轴对称体空泡流的噪声特性与空泡界面瞬态特征[J].上海交通大学学报,2000,34(8):1026-1030. Gu Wei,He Yousheng,Hu Tianqun.Noise feature of cavitation flows on axisymmetric bodies and transient characteristics of cavity interface[J].Journal of Shanghai Jiaotong University,2000,34(8):1026-1030.

[4]童秉纲,尹协远,朱克勤,等.涡运动理论[M].合肥:中国科学技术大学出版社,2009. Tong Binggang,Yin Xieyuan,Zhu Keqin.Vortex dynamics theroy[M].Anhui:Press of University of Science and Technology of China,2009.

[5]Wu J Z,Ma H Y and Zhou M D.Vorticity and vortex dynamics[M].Berlin:Springer-Verlag,2006.

[6]Li C Y,Ceccio S L.Interaction of single travelling bubbles with the boundary layer and attached cavitation[J].Journal of Fluid Mechanics,1996,322:329-353.

[7]Belahadji B,Franc P,Michel M.Cavitation in the rotational structures of a turbulent wake[J].J Fluid Mech.,1995,287: 383-403.

[8]Arakeri V H,Acosta A J.Viscous effects in the inception of cavitation on axisymmetric bodies[J].Journal of Fluids Engineering,1973,95(4):519-527.

[9]Gopalan S,Kate J.Flow structure and modeling issues in the closure region of attached cavitation[J].Physics of Fluids, 2000,12(4):3414-3431.

[10]Wang G,Senocak I,Shyy W.Dynamics of attached turbulent cavitating Flows[J].Progress in Aerospace Sciences,2001, 37:551-581.

[11]张 博,王国玉,黄 彪,等.云状空化非定常脱落机理的数值与实验研究[J].力学学报,2009,41(5):651-659. Zhang Bo,Wang Guoyu,Huang Biao,et al.Numerical and experimental studies on unsteady shedding mechanisms of cloud cavitation[J].Chinese Journal of Theoretical and Applied Mechanics,2009,41(5):651-659.

[12]Huang B,Wang G.Experimental and numerical investigation of unsteady cavitating flows through a 2D hydrofoil[J].Science China Technological Science,2011,54(7):1801-1812.

[13]Nagendra Dittakavi,Aditya Chunekar,Steven Frankel.Large eddy simulation of turbulent-cavitation interactions in a venturi nozzle[J].Journal of Fluids Engineering,2010,132(12):

[14]黄 彪.非定常空化流动机理及数值计算模型研究[D].北京:北京理工大学,2012. Huang Biao.Physical and numerical investigation of unsteady cavitating flows[D].Beijing:Beijing Institute of Technology, 2012.

[15]Kubota A,Kato H,Yamaguchi H.A new modeling of cavitating flows:A numerical study of unsteady cavitation on a hydrofoil section[J].Journal of Fluid Mechanics,1992,240(1):59-96.

[16]Smirnov P E,Menter F R.Sensitization of the SST turbulence model to rotation and curvature by applying the Spalart-Shur correction term[C]//ASME Paper GT 2008-50480.Berlin,Germany,2008.

[17]Haller G.An objective definition of a vortex[J].J Fluid Mech,2005,525:1-26.

Vortex dynamic analysis of unsteady cavitating flows around a hydrofoil

ZHAO Yu,WANG Guo-yu,HUANG Biao,HU Chang-li,CHEN Guang-hao,WU Qin
(School of Mechanical and Vehicular Engineering,Beijing Institute of Technology,Beijing 100081,China)

Unsteady cavitating flows around a hydrofoil are calculated.Vortex dynamics are used to analyze the flow structure.Local swirling correction turbulence model is adopted to close the governing equations.Validation of the LSC turbulence model can be reached comparing with experimental results.Vorticity transport equations are conducted to analyze the influence of cavitation on vorticity transport process. The results show that:The LSC turbulence model can reflect the influence of the local swirling on cavitating flows well.It can get a good compromise with experimental results.Unsteady cavitation process contains developing of complex vortex structures.Re-entrant jet near trailing of attached cavity leads to distinct changes of velocity gradient,which has great influence on production and dissipation of vorticity.Intensive mass transfer between liquid and vapor phase may induce the volumetric expansion or contraction and baroclinic torque.Two factors mentioned above keep changing and result in unsteady distributions of vorticity.Cavitation may lead to unsteady mass and momentum transfers between liquid and vapor phase,which will change turbulence and vorticity distributions distinctly.

cavitation;turbulence model;unsteady;vortex dynamics analysis

O35 TV131.2

A

10.3969/j.issn.1007-7294.2015.08.002

1007-7294（2015）08-0896-09

2015-04-23

国家自然科学基金重点项目（51239005，51479002）

赵 宇（1989-），男，博士研究生；

王国玉（1961-），男 教授，博士生导师，E-mail:wangguoyu@bit.edu.com。