三维水翼梢涡气核运动特性的数值模拟

郭春雨， 薛嵘， 胡健， 王恋舟

(1.哈尔滨工程大学 船舶工程学院，黑龙江 哈尔滨 150001； 2.西南交通大学 力学与工程学院，四川 成都 611756)

梢涡在有限展长的翼型绕流现象中普遍存在。当流动流过有攻角的有限长水翼的端部，水翼上下表面的压力差会驱使流动从下表面的高压侧流向上表面的低压侧，从而形成三维的旋涡，我们把它称为梢涡[1]。梢涡具有强烈的旋转效应，使得涡核中心的压力值显著降低。在水动力学领域，梢涡引起的压降会驱使水中的气核卷入涡核内部，产生空化。梢涡空化会降低水翼结构的水动力性能，并会伴随出现强烈的噪声。为了避免梢涡空化带来的危害，需要对它的形成机理进行分析。

学者们进行了大量的试验来研究梢涡的细微结构。Giuni等[2-3]通过烟迹试验，得出梢涡是由起始于压力面的主涡和2个分别起始于吸力面和压力面的次级涡演化形成的。徐良浩等[4]使用2D-PIV和3D-PIV对梢涡的二维和三维速度场进行了测量，获得了不同工况下的涡核位置、尺寸及速度分布，得出了涡核在三维空间内的摆动特性。Maine等[5]通过对3种不同翼型剖面的实验观测，研究了水翼的纵倾和侧斜对梢涡形成的影响。梢涡试验往往伴随着空化初生现象。Arndt等[6]通过LDV试验指出，梢涡空化发生的必要条件是分布在水中的气核进入涡核中心并发展为足够大的气泡。在这个基础上，Arndt等[7]进一步研究了压力变化和水中含气量对梢涡空化的影响。他们得出在含气量多的水中，梢涡空化多发生于水翼梢部；而在含气量较少的水中，梢涡空化的起始位置在梢部流动后端。

随着计算机运算能力的逐步发展，采用各种湍流模型来模拟梢涡逐步成为可能。相比各种试验手段，数值模拟可以观察到更加细微的涡结构，从而更好地去观察梢涡的演化，进而了解梢涡空化初生的规律[8]。尽管如此，由于梢涡尺寸很小，压力梯度很大，模拟梢涡是一件具有挑战性的工作[9]。Decaix等[10]采用RANS模型对截面形状为NACA 0009的水翼梢涡进行了数值模拟，数值结果表明RANS方法可以准确模拟出梢涡的涡核位置、平均速度场等的时均特性。蒲汲君等[11]使用k-w、DES和LES模型分别对水翼的梢涡流场进行了数值模拟，并使用气泡静力平衡方程计算了空化初生现象。结果表明LES模型可以很好地模拟梢涡处流场，此外应用气泡静力方程可以准确计算初始梢涡空泡数。当涉及空化初生问题时，可以耦合球形空泡模型来模拟气泡运动[12]。Rayleigh[13]建立了球形气泡模型，用于描述球形气泡在变化的压力场中的半径变化。许多研究致力于分析气泡及颗粒在水中的运动，Maxey等[14]给出了一种描述刚性球体的运动方程，可以很好地描述球形空泡的运动规律。Park等[15]利用他的计算方程，计算了气泡在水翼梢部的运动。

本文以剖面为NACA 0015的三维水翼为研究对象，采用大涡模拟(large eddy simulation, LES)方法来模拟水翼梢涡。将计算结果与试验值[16]进行对比，满足精度要求之后，在这个基础上开展基于离散元 (discrete element method, DEM)模型的气核在水中运动研究，分析气核在梢涡作用下的演化规律。整个计算工作在数值模拟软件STAR-CCM+下完成。

1 几何模型与数值模型及网格划分

1.1 几何模型

几何模型与文献[17]中试验所用翼型保持一致。如图1所示，翼型截面形状为NACA 0015，弦长C=0.15 m，展弦比6.6，攻角12°，端部装有圆形导流罩。计算域布置以及边界条件如图2所示。为了节省计算资源，将水翼的展长选为3.3倍弦长，远离梢部的一端贴近壁面，并将壁面设置为对称边界条件，从而节省一半的网格数量。由此，该计算域的宽为4.5C，高为4.5C，水翼前端距速度入口距离为3.3C，后端距压力出口距离9C，其余壁面均设置为无滑移壁面条件。计算环境为25 ℃的水，此时密度998 kg/m3，动力粘度0.000 889 9 kg/(m·s)，流速10 m/s，从而雷诺数Re=1.5×106。

1.2 湍流模型

本研究采用大涡模拟(large eddy simulation, LES)求解流体域中任何位置的空间湍流结构。网格尺寸超出选用尺寸的大尺度旋涡进行直接求解，小于选用网格尺寸的细微旋涡结构采用亚格子进行模拟。从而在复杂的流动模拟中得到细微旋涡结构的流动图像。

图1 水翼模型及导流罩Fig.1 Hydrofoil model with round end cap

图2 计算域布置条件Fig.2 Compute region conditions

(1)

式中φ表示速度分量、压力、能量或组分浓度。

将分解的求解变量插入纳维-斯托克斯方程可得到已滤波物理量的方程。滤波后的质量和动量的传输方程可以写为：

(2)

(3)

湍流应力张量现在表示为亚网格尺度应力。这些应力产生自较大的已求解涡流与较小的未解析涡流之间的相互作用，使用Boussinesq近似进行建模，为：

(4)

(5)

亚网格尺度湍流粘度μt必须由能够解释小涡流对已求解流体的作用的亚网格尺度模型来描述。采用WALE亚网格尺度模型，使用混合长度假设来对亚网格尺度应力建模：

μt=ρΔ2Sw

(6)

式中：Δ为长度尺度或网格滤波器宽度，使用网格单元体积定义如下：

Δ=min(κd,CwV1/3)

(7)

式中：Cw为模型系数，典型值介于0.5(对于均匀各向同性衰变湍流)～0.325(对于通道流)之间;κ为冯卡门常数；Sw的定义见文献[17]。

1.3 DEM模型

本研究采用DEM模型对气核的运动特性进行模拟。离散元方法(DEM)是一种基于拉格朗日方法的数值方法，通常用来模拟大量相互作用的离散对象(通常是颗粒)如何运动。根据格雷斯等[18]的研究，气核直径通常在30～300 μm，在这一尺度下，气核在运动过程中可以维持较为稳定的球形外形，从而在本研究中，将气核简化为不可形变的球形颗粒，应用DEM方法对气核进行建模。

DEM方法以牛顿第二定律作为控制方程：

(8)

(9)

式中：mi、vi和ωi分别表示颗粒i的质量、速度和角速度；Ii为颗粒i的转动惯量，对于球形粒子，Ii=2miRi2/5，Ri为颗粒i的半径；Fg=mig为颗粒i的重力；,Ffluid为流体对颗粒i的作用力(包括阻力、升力、附加质量力和浮力等)；Fij为颗粒i与颗粒j或墙壁之间的碰撞力以及其他作用在粒子上的非接触力；Tij为接触力矩，表示除颗粒重力以外的接触力在粒子上产生的扭矩。

STAR-CCM+中的DEM模型在拉格朗日架构内实现，并使用拉格朗日相模型对粒子进行定义。本文使用DEM方法对气核进行模拟，并与CFD结合使用，从而求解由颗粒组成的离散相的气核在流体介质中的运动。其中，流场的求解由大涡模拟完成；气核的运动基于拉格朗日方法,由式(8)、(9)进行求解。STAR-CCM+将对DEM粒子和流体之间的相互作用进行建模，并将LES计算得到的流场中各网格处的力代入到式(8)、(9)中，对气核的运动进行求解。

1.4 网格划分

采用Star-CCM+对水翼模型进行面网格重构，形成表面三角化良好的高质量面网格，以面网格为基础生成带有边界层和切割体网格的非结构体网格。边界层共划分12层棱柱层网格，并将y+值控制在1以下，从而满足LES模拟的计算要求。对水翼表面和梢涡区域的体网格进行加密处理，最终确认计算域网格总数约为780万，网格细节如图3所示。

图3 全局网格、加密区和边界层Fig.3 Overall mesh condition, refining region and boundary layer mesh

2 计算结果分析

2.1 数值验证

为验证所采用的数值模型是否满足计算要求，以数值模拟结果与文献[16]中的相关试验工况的结果进行对比，来验证模拟结果的准确性。监测水翼表面受到的力在x和y方向上的分解值，即阻力和升力，并将升力换算为力系数并与试验值进行对比。可以得出本次数值模拟得到的升力系数为1.18，对比文献[16]中对应的试验结果的升力系数为1.24，相对误差在4.8%以内，可以认为本文的计算结果与文献中的结果较为吻合，从而本文采用的数值模型是可靠的。

图4为梢涡附近流场无因次垂向速度沿弦长方向分布结果与试验结果的对比。在水翼随边后端0.2倍弦长距离的位置布置一排监测点，用于监测梢涡附近的垂向速度值，如图5所示。本研究采用的大涡模拟方法是一种瞬态的模拟方法，因此各监测点需要将监测一段时间后的速度值进行平均化，再将平均化后的垂向速度值UZ除以来流速度U∞进行无因次化，由于梢涡轨迹近似于X轴方向，从而认为UZ/U∞的值作为梢涡的无因次切向速度。从图5中可以看出，数值模拟对涡核位置附近的切向速度值捕捉与试验值较为相符，除了在峰值和谷值处有一定的偏差。这是因为无量纲化速度的极值位于梢涡涡核位置，由于涡核高强度的旋转带来了切向速度的极值。试验中使用的涡量计布置位置为一系列等间距点，在峰值附近，布置点较少而漏掉了极值；而模拟中在涡核位置附近布置了更多的监测点，从而监测到涡核的极值。通过对数值模拟中速度场的验证，可以认为本研究对于流场的模拟具有较高的精确度。

图4 水翼后端0.2倍弦长处梢涡附近垂向速度对比Fig.4 Comparison of vertical velocity across the tip vortex at r=0.2C

图5 水翼后端0.2倍弦长处垂向速度监测位置Fig.5 The location of vertical velocity monitoring at r=0.2C

2.2 结果分析

2.2.1 流场结构

本节从梢涡结构方面来分析梢涡流场数值模拟的结果。使用Q判据来可视化涡结构，Q判据的计算公式为：

(10)

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

图6为Q=10 000 s-2的瞬时梢涡等值面，并采用涡量幅值进行着色。为了展现整个梢涡位置涡量场的细节，分别从水翼侧面和后方视角展示了涡量场的整体与局部。由图6(a)展现的整体涡结构可知，本文采用的LES模型结合水翼后端位置网格加密，可以捕捉到精细的梢涡结构，且梢涡主涡的强度值在1 600 s-1附近可以保持在5倍弦长之上的距离，这同Bailey等[19]的试验结果相符。从等值面涡量场可得，本研究对于梢涡区域的网格细化程度满足分析要求。

图6(b)展示了近梢部端涡量场的涡结构。从图中可以看出，梢部涡系由2部分组成，序号为1的涡来源于水翼压力面梢部的位置边界层发生分离，进而产生涡强很高的梢涡；序号为2的涡来源于水翼吸力面靠近梢部的位置，这部分的流动在从水翼随边分离后，在梢涡的强旋转效应下形成了包裹在梢涡外边的涡系。

图7为不同截面上的涡量云图，截面的选取参考图8，选取水翼随边为起始位置，每隔0.5倍弦长距离共作6个截面，形成图7的截面。从涡量云图可以观察到，梢涡位置具有多重涡结构组成的涡系。压力面形成一个主涡和一个二次涡，在向下游迁移到1.5倍弦长附近发生融合；吸力面发生流动分离并形成多个较弱的涡心，在向下游发展的过程中涡量逐渐耗散。

图6 Q=10 000 s-2等值面的涡量场及细节展示Fig.6 Vortex field and detail display of iso-surface with Q=10 000 s-2

2.2.2 气核运动特性分析

本节在得到稳定发展的涡量场的基础上，向流场中加入DEM粒子所表示的气核，来研究水翼梢部气核的空间演化过程。

拉格朗日相的DEM粒子需要通过喷射器加入到欧拉相的流场中进行求解。如图9为DEM粒子的喷射点位置，在水翼的端部，沿弦长方向每隔1/10弦长设置一排喷射点，这些喷射点沿圆形导流罩长度方向间距为1 mm。随后在拉格朗日相下设置气核的物理属性，在本研究中，气核的半径为100 μm，密度为气体密度。如此设置后在流场中喷射点的位置各出现1个气核，随后气核将在DEM方程的控制下，基于受到流场中的各种力随梢涡向水翼后端迁移。

图7 水翼随边后端涡量云图截面分布Fig.7 The distribution of vorticity contour section behind hydrofoil′s trailing edge

图8 不同截面上涡量云图Fig.8 Vorticity contour in typical sections

图9 DEM粒子喷射点位置Fig.9 DEM injector positions distribution

如图10所示，选取喷射点位于端部中轴线并距端部垂直距离为1 mm处的气核流迹进行分析。图中用速度场对气核流迹进行着色，深色线条为和气核初始位置相同的流线。从流线轨迹可知，初始流动位置越靠近水翼前缘，在向下游迁移的过程中会出现远离梢涡核心的趋势；而初始位置接近水翼后缘的流动会被梢涡捕获。对比经过各气核初始位置的流线，可以对这一现象进行解释。由于梢涡在水翼压力面后端开始产生，因此水翼端部后半端的流动受旋涡影响较大，从而这些位置的流线紧紧跟随梢涡运动；水翼端部前半段产生的流线表明，经过这部分的流动会沿水翼展向迁移，在随边与水翼分离后，受到图6(b)中涡系2的影响，随着包裹在涡核外部的涡系旋转运动。气核具有一定的流线跟随性，从而伴随着2种不同的流线运动状态，气核也具有2种不同的运动状态。

图10 初始位置距端部1 mm处气核流迹及流线对比Fig.10 Comparison between nuclei track and streamline with initial location from round end=1 mm

气核的流迹虽然表现出一定的流线跟随性，但是与流线还是有一些区别。从图6(b)细节可以看出，气核的迁移轨迹极易受到涡强的影响。在水翼压力面梢涡的主涡刚形成的位置产生的气核，会迅速迁移到旋涡中心，跟随涡核向下流旋转运动；在水翼前端部产生的气核，在随流线迁移到与水翼分离后，迁移轨迹位于梢涡外部包裹的涡系内。可以得出，气核先随流线向后迁移，当轨迹附近出现强涡量区域后，气核会向涡强较强的位置进行迁移，从而不再跟随流线。

接下来选取部分气核，对它们流迹的特征进行分析。如图11分别为1、2、5、7、9号的喷射器发射的气核流迹，在向下游发展的过程中垂向速度的变化曲线。对迁移距离除以弦长得到L/C，垂向速度除以来流速度得到V/U进行无量纲化处理。由于水平长度7C以后的网格位于梢涡加密区外，因此不再取之后的速度进行分析。从曲线中可以看出，为1号和2号的气核由于流迹位于旋涡核心外部，在迁移的过程中速度振荡频率较小。而其余气核流迹位于涡核内部，可以按照振荡特性将轨迹分为3个阶段。如图11(c)所示，第1阶段是从气核刚进入流场到迁移到水翼尾端位置，这一阶段气核开始进入刚形成的梢涡内部，垂向速度的变化表现出较大的随机性。这是由于这一阶段旋涡刚开始形成，而气核不再跟随流线运动并开始进入涡核内部，造成了无规律的垂向速度变化。第2阶段是气核从到水翼尾端位置迁移到距尾端水平距离为4倍弦长处。这一阶段涡核具有相对较大的涡强，气核随涡核向下游段迁移过程中具有较为稳定的振荡频率和幅值。第3阶段位于水平距离4倍弦长之后，在这一阶段，涡核的强度逐步下降，而气核的垂向速度变化频率也随之变缓慢。

图11 部分气核流迹上的垂向速度分布Fig.11 Vertical velocity distribution in typical nuclei track

图12中，(a)～(e)分别为1、2、5、7、9号的喷射器发射出的气核运动过程中相对压力的变化曲线。对迁移距离除以弦长得到L/C，对计算得到的相对压力转化为相对压力系数Cp，完成无量纲化处理。压力变化曲线与速度变化曲线有一定的相似之处，位于梢涡核心位置处和梢涡核心以外区域的流迹在压力变化上具有2种明显不同的趋势。1号和2号的气核在迁移过程中，压力首先经历了一段明显的上升，在迁移到水翼后端1倍弦长的位置处相对压力达到峰值；随后下降到-0.3左右。而5、7和9号的气核在迁移到水翼后端1倍弦长的位置处压力先急剧下降，随后再急剧上升，最后平稳上升。从气核的迁移轨迹来看，1号和2号的气核在刚开始迁移的过程中随流线向展长的方向进行迁移，因而经历了一次压力上升；随后被梢涡外端涡系捕获后，压力再次下降并稳定。而5、7、9号的气核刚开始迁移时被涡核所捕获，故而压力急剧下降；随后在随涡核向后迁移的过程中，随着涡强降低，压力也随之上升；而在涡强下降到一定程度后，下降趋势逐步变慢，此时压力上升趋势也开始放缓。

图12 部分气核流迹上的相对压力分布Fig.12 Relative pressure distribution in typical nuclei track

3 结论

1)采用LES模型对带导流罩的三维水翼受到的升阻力、无因次速度分布进行计算，与文献中数据符合较好，计算结果具有良好的可信度。

2)带导流罩的圆形导边水翼的梢涡流场有着独特的结构。在压力面出现主涡和二次涡，在吸力面出现一系列小涡。主涡和二次涡在发展的过程中合并，形成梢涡核心，而吸力面形成的涡系在向后迁移的过程中，包裹在核心外部，形成总体的涡系。

3)气核在迁移过程中，在水翼首部出现的气核的迁移轨迹位于梢涡核心外，而水翼首部后端的气核迁移轨迹位于梢涡核心内。此外，涡核会对气核的迁移产生很大的影响，使其偏离流线的方向；气核被涡核捕获后，垂向速度会有较大的振荡值，受到的压力也会产生显著的变化。