动桨射流数值模拟

张玉倩，曾剑锋，张弘强，张星烨，王志纯

(1.天津天科工程监理咨询事务所，天津 300452；2.南开大学 软件学院，天津 300071；3.天津大学 水利工程仿真与安全国家重点实验室，天津 300072；4.天津大学 建筑工程学院，天津 300072；5.法国达飞海运集团公司，马赛 13000；6.交通运输部天津水运工程科学研究所，天津 300456)

随着全球航运行业的不断发展，船舶螺旋桨射流也愈发受到国内外研究人员的关注。无论是分析船桨、船舵、船体之间的相互干扰，或是探讨附加装置的节能效果，都要对螺旋桨射流流场作深入研究。

国内外专家学者关于螺旋桨射流的研究对象，大部分都是固定于水体中的螺旋桨(以下简称为定桨)。其中，Albertson等[1]基于轴向动量理论，采用平面自由射流理论研究螺旋桨射流，得到了平面射流初速度公式，为之后的研究奠定了理论基础。之后，Blaauw与Kaa[2]、Berger与Cederwall[3]、Verhey[4]、Fuehrer与Romisch[5]在Albertson的理论基础上，通过建立螺旋桨物理模型，对Albertson理论作了修正，提出了一系列推求射流流场各位置流速的计算方法。Stewart[6]、Hamill与Johnston[7]、Hashmi[8]、Lam等[9]对螺旋桨射流流场的起始流速及流速衰减情况作了进一步研究，通过物模试验对Albertson等经验公式作了修正，得到了更接近于真实情况的螺旋桨射流流速分布规律。

然而，在实船航行过程中，螺旋桨是推动船体并随之一起运动的，对于某些特定条件下的定桨物模试验，其对实船螺旋桨射流真实情况的反映会有所欠缺。温春鹏等[10]在分析适航资源分布的同时，发现船舶在航行过程中会对浮泥产生扰动，尾流有较为明显的浑浊带。张定军等[11]通过船舶操纵模拟试验，研究了船舶港内作业对港池水域的要求。目前，预报螺旋桨水动力特性的计算方法主要有升力线法、升力面法、面元法、CFD(Computational Fluid Dynamics)方法等。虽然CFD方法起步较晚，但是其适用范围更广、计算精度更高。因此，本文在前人的定桨物模试验基础上，采用CFD软件FLUENT对推进过程中的螺旋桨(以下简称为动桨)射流进行三维数值计算，以得到更符合实际情况的动桨射流流场分布规律。同时，采用自航船模进行动桨射流试验，通过高精度的PIV(Particle Image Velocimetry)系统进行流场测量，以验证数值计算方法的可行性。

常规CFD计算大多涉及的是稳态问题，然而对于旋转、推进同时进行的螺旋桨，其必然存在运动边界的问题。因此，本文采用动网格技术，对存在运动边界的非定常流动进行数值计算。另外，在螺旋桨射流流场中，射流轴向速度占据主要地位[12]，因此本文所提及的射流均表示轴向射流。

1 计算对象

1.1 计算域

图1 计算域Fig.1 Calculation domain

由于实船螺旋桨的外观、尺寸及空间形态极其复杂，且本文研究对象为螺旋桨桨后射流流场，并未涉及螺旋桨推力、扭矩、表面压力分布等水动力特性。因此，在计算域的建立过程中对螺旋桨进行概化，采用边界类型中的Fan模型替代直径(以下简称为桨径，用D表示)0.015 m、转速3 000 rpm的螺旋桨。计算域如图1所示，Gx、Gy、Gz轴的交点为三维坐标原点，螺旋桨中轴线(以下简称为桨轴线)坐标为(x，0.075，-0.015)，坐标轴定义为：x轴为纵向，与桨轴线一致，指向下游；y轴为横向，与桨平面一致；z轴为垂向，服从右手系。

若螺旋桨射流在发展过程中过早碰到水槽边壁，则射流流速将会发生反射，从而影响整个射流流场的计算结果。为减小边壁反射对射流流场的不利影响，并综合计算机的硬件配置及计算精度，计算域采用长0.9 m、宽0.15 m、高0.07 m的长方体区域。其中，z=-0.045～0的下部区域为水体部分(即研究区域)，z=0～0.025的上部区域为空气部分。

1.2 计算网格

图2 计算网格Fig.2 Computational grids

计算网格划分是数值计算过程中最为耗时的环节之一，但也是直接影响计算精度与效率的关键因素之一。计算网格如图2所示，本文采用局部加密方法对网格进行划分，不但可以提高计算精度，而且可以避免流场变化平缓区域的计算资源浪费。

在紧邻螺旋桨的区域内建立一个长0.9 m、宽0.017 m、高0.017 m的长方体，将螺旋桨的运动轨迹完全包裹在内，并采用适应性较好的四面体网格对该长方体及其相邻区域进行加密。在向外平缓过渡的过程中，网格尺寸逐渐增大、结构逐渐稀疏，并采用计算较易收敛的六面体网格。完成计算网格划分，网格总数近90万个。

2 计算方法

2.1 控制方程

对于计算域内形式复杂的三维非定常不可压缩湍流流动，本文采用Reynolds平均法，即对非稳态Navier-Stokes(N-S)方程采用时间平均法，得到时均形式的控制方程。

目前，Reynolds平均法是使用最为广泛的湍流计算方法，其核心思路是不直接求解瞬时N-S方程，而是求解时均化Reynolds方程[13]。实际工程关注的重点是湍流所引起的平均流场变化与流场整体分布，正如本文着重研究桨后射流流场。

2.2 动网格技术

动网格技术是用于处理运动边界所引起的非定常流动的常用方法，该方法在各个研究领域内都能够得到较为广泛的应用[14]。

在数值计算过程中，为了使网格能够适应运动边界移动所引起的变化，需要对计算网格进行修正。本文采用弹性光滑与局部重划相结合的方法，对计算网格进行修正。

2.3 UDF技术

常规螺旋桨射流数值计算的研究对象一般也是定桨，需要通过改变来流大小来控制进速系数(用J表示)，以替代实际情况下的进速系数变化。

为真实模拟动桨射流，实现真正意义上的“进速系数”，本文采用可被动态连接到FLUENT求解器中的UDF(User Defined Function)技术。采用编译器VC++对UDF进行编译，控制螺旋桨沿x轴正方向分别作J=0.2的匀速直线运动。计算前，对整个过程中的动网格进行更新检查，通过调整相关参数使得网格质量满足计算要求。

2.4 湍流模式与多相流模式

大部分实际工程中，流体运动都处于湍流状态，螺旋桨射流也不例外。相对于其它湍流模式，RNGk-ε模型更适合对旋流、射流等进行计算[15]。因此，本文采用RNGk-ε模式作模拟。

处理多相流的计算方法，主要有欧拉-拉格朗日法、欧拉-欧拉法。在FLUENT中，有VOF模式、混合模式、欧拉模式等3种欧拉-欧拉多相流模式可供选择。由于计算内容包含流体与固体表面的相互作用，根据FLUENT对多相流选择的基本原则，采用VOF模式作模拟。

2.5 边界条件设置

在计算域中，入口为Velocity_Inlet，出口为Outflow，螺旋桨与流体的接触面为Fan，水槽底面、边壁均为Wall，动静交接面传递方式采用混合面法[16]。

对流项采用二阶迎风格式进行离散，对扩散项采用具有二阶精度的中心差分方式进行离散。计算初期，分别选取0.002 s、0.01 s两个时间步长进行计算。计算结果显示，以上两个时间步长更新后的网格质量均能满足计算精度的要求，且计算结果较为接近，说明采用0.01 s的时间步长可以得到与时间步长无关的数值解。因此，本文选取0.01 s作为时间步长进行计算，每一时间步长内迭代50次。

图3 纵截面射流流速等值图Fig.3 Jet velocity contours in 4 longitudinal sections

3 计算结果分析

3.1 进速系数为0.2时的流场分析

3.1.1 纵截面流场分析

计算完成之后，对进速系数J=0.2时的射流流场进行分析。从y=0.075的xoz平面(即包含桨轴线的xoz平面，以下简称为中轴面)开始，依次提取y=0.075+1/4D=0.078 75、y=0.075+1/2D=0.082 5、y=0.075+3/4D=0.086 25(即桨侧距离桨轴线1/4、1/2、3/4倍桨径)等3个xoz平面。纵截面射流流速等值图如图3所示，通过这4幅纵截面图，对不同横向位置的二维射流流场进行分析。

由图3-a可得，射流流场以桨轴线为中心呈对称分布，流速分布呈双峰型，断面最大射流流速位于桨叶中部。在远离桨平面的扩散过程中，射流流速不断减小、衰减程度逐渐减小。此外，射流流速在桨后极短距离内激增，并于0.37D位置达到最大流速，该计算结果也与“Hamill发现最大流速值在0.35D区域内无衰减”[17]这一研究成果相吻合。

由图3-b可得，纵截面横向移动1/4D距离，射流流场仍以桨轴线为中心呈对称分布，断面最大射流流速位置已由桨叶中部逐渐偏向桨中心，且双峰型已逐渐向单峰型衰减。由图3-c可得，纵截面横向移动1/2D距离，双峰型已完全衰减为单峰型，断面最大射流流速位于桨中心。由图3-d可得，纵截面横向移动3/4D距离，流速分布呈细长椭圆型，断面最大射流流速位于桨中心。

纵观4幅纵截面射流流速等值图，随着纵截面距中轴面横向距离的增大，射流流速及其衰减程度不断减小，流场影响范围逐渐减小，双峰型特征愈发减弱，断面最大射流流速位置已由桨叶中部逐渐集中于桨中心。

3.1.2 中轴面流速变化分析

图4 中轴面射流流速变化图Fig.4 Jet velocity graph of center axial plane

对于射流流场最具有代表性、流速最为显著的中轴面，有必要对其流速变化进行深入分析。在中轴面上，沿水深方向依次提取4条垂直于z轴的直线，往深处依次为Z1(x，0.075，-D)、Z2(x，0.075，-1.25D)、Z3(x，0.075，-1.5D)、Z4(x，0.075，-1.75D)，由此作出不同水深位置的射流流速曲线。中轴面射流流速沿水深方向的变化曲线图如图4所示，图中Vx表示螺旋桨射流的轴向速度大小(下同)，x/D表示桨后若干倍桨径的纵向位置。

由图4可得，不同水深位置的4条射流流速曲线，在分布规律上均保持一致：从桨平面开始，射流流速激增，并于桨后0.37D的位置达到峰值；峰值过后，射流流速迅速减小，衰减程度较大；逐渐的，衰减程度减小，射流流速缓慢减小。并且，桨后越远位置，射流流速衰减程度越小。

图5 桨后三维射流流线图Fig.5 3D diagram of jet streamline behind propeller

由于断面最大射流流速位于桨叶中部，因此，水深位置为桨叶中部的直线Z2，相对于水深位置为桨轴线的直线Z1，射流流速更大，且为相同x轴位置的流速最大值。并且，随着水深不断增大，射流流速衰减程度逐渐减小。

3.1.3 三维流场分析

为了能够更加直观的对桨后射流流场进行分析，提取桨后及其相邻区域内的三维射流流线图，对其总体分布趋势进行分析，桨后三维射流流线图如图5所示(图例与图3-e相同)。

由图5可得，射流流线以桨轴线为中轴线，围绕其向后呈螺旋式发展。桨平面后的射流流速均高于外流域，这也是螺旋桨向后推水、水体反作用于螺旋桨而产生推力的结果。较大的流速变化均高度集中于桨平面及桨后流域，随着远离桨平面，外流域的变化与影响程度迅速减小。内半径处，射流在螺旋桨作用下向内收缩，直径也逐渐减小，反映了螺旋桨对水体的抽吸作用。此外，该图也与Lam等[15]用CFD方法模拟出的螺旋桨射流流线图相吻合。

3.2 模型验证

本文采用自航船模进行动桨射流试验，并采用高精度的PIV系统进行流场测量，以验证数值仿真方法的可行性。本文采用由美国TSI公司生产的立体PIV流场测量系统，该系统由激光器、同步器、光电编码器、CCD摄像机等精密部件组成。采用的自航船模由国内船模试验研究经验非常丰富的交通运输部天津水运工程科学研究所研制，并在吃水深度、静水航速、零舵角直航操纵性、“Z”型操纵性等方面与实船进行了相似性率定。

图6 物模试验布置Fig.6 Layout of physical model test

在长20.0 m、宽0.4 m、高0.5 m的试验水槽中，水深0.1 m，静水；自航船模螺旋桨桨径0.015 m，转速3 000 rpm，桨中心位于水面以下0.015 m位置，螺旋桨推动船体以进速系数0.2匀速直线通过水槽中部试验段(即PIV系统可拍摄段)，物模试验布置如图6所示。在试验数据中提取中轴面上位于桨后0.015 m、0.030 m、0.060 m、0.150 m(即桨后1、2、4、10倍桨径)等4条垂直于x轴的直线，得到桨后不同距离射流沿水深方向的流速分布。相应的，在数值仿真结果中提取中轴面上位于桨后0.015 m、0.030 m、0.060 m、0.150 m等4条垂直于x轴的直线，与物模试验所得到的桨后不同距离射流沿水深方向的流速分布进行对比验证。

物模试验与数值计算转速(n)比尺为1:1、桨径(d)比尺为1:1，且射流流速V～nd[18]，因此射流流速(V)比尺为1:1。物模试验与数值计算所得到的射流流速分布比较如图7所示，图中z/D表示水面以下若干倍桨径的垂向位置。

7-a 桨后1倍桨径位置射流流速分布7-b 桨后2倍桨径位置射流流速分布7-c 桨后4倍桨径位置射流流速分布7-d 桨后10倍桨径位置射流流速分布图7 射流流速分布比较Fig.7 Comparisons of jet velocity distributions

由图7-a、7-b可得，相对于试验结果，模拟结果中流速沿程衰减稍快、断面流速最大值稍大，但两者总体趋势吻合较好。由图7-c可得，相对于试验结果，模拟结果中流速沿程衰减稍慢、断面流速最大值稍大，两者总体趋势依旧吻合较好。由图7-d可得，相对于试验结果，模拟结果中仅流速沿程衰减稍慢，其余两者均吻合较好。

纵观4幅射流流速分布对比图，达到第一个流速峰值之前，两条曲线的斜率基本一致，流速及其衰减程度均基本相同；经过第一个流速峰值之后大约0.6D的垂向深度范围内，数值计算所得到的流速略大于物模试验结果，但偏差均在5%以内，且二者流速衰减程度基本相同；再往下的垂向范围内，两条曲线的斜率又基本趋于一致，流速及其衰减程度均基本相同。

数值计算与物模试验中的射流流速分布总体趋势保持一致，最大射流流速近似相等，射流流速衰减程度基本相同。至于二者在部分范围内的射流流速、最大射流流速出现位置、轴对称情况等方面还存在些许差异，主要有两个方面的原因：一方面，物模试验过程中会出现偶然性误差，因此会有部分数据偏离整体趋势，甚至不符合实际规律；另一方面，数值计算严格按照选定公式进行迭代计算，因此结果较物模试验更为理想化，数值亦会偏大。

对于水流流向复杂、流速多变的螺旋桨射流流场来说，以上些许差异，均可允许。综上所述，数值计算结果与物模试验结果基本吻合，验证了数值计算方法的可行性。

4 结论

实船航行过程中，由于螺旋桨运动使得桨平面、桨后流场不断发生变化，常规CFD稳态分析方法的计算精度较低，不能真实有效地模拟动桨射流情况。因此，本文采用UDF及动网格技术，对动桨射流进行三维数值计算，得出的主要结论如下：

(1)射流中轴面上，流场以桨轴线为中心呈对称分布，流速分布呈双峰型，且断面最大流速位于桨叶中部。射流流速在桨后极短距离内激增，并于0.37D位置达到最大流速。

(2)偏离桨轴线的横向距离越大，射流流速及其衰减程度越小，流场影响范围越小，双峰型特征愈发减弱，断面最大流速由桨叶中部逐渐集中于桨中心。偏离桨轴线的垂向距离增大，射流流速先增后减，衰减程度逐渐减小。

(3)采用自航船模进行动桨射流试验，并采用高精度的PIV系统进行流场测量，验证了数值计算方法的可行性，说明本文提出的方法适用于动桨射流研究，并能为较复杂工况下的动桨射流研究提供一定的理论依据。

致谢：参与本研究工作的还有长沙理工大学水利工程学院胡旭跃教授、沈小雄教授，在此并致谢意。