喷水推进器进流面获取方法及其应用

郭 军， 陈作钢， 戴原星， 陈建平

（1.上海交通大学 海洋工程国家重点实验室，上海200240；2.上海交通大学 高新船舶与深海开发装备协同创新中心，上海200240；3.上海交通大学 船舶海洋与建筑工程学院，上海200240；4.中国船舶及海洋工程设计研究院，上海200011；5.喷水推进技术重点实验室，上海200011）

喷水推进器是利用喷出水流的反作用力来推动船舶前进的一种特殊推进器，具有高航速时推进效率高、操纵性好、抗空化、低噪声等优势，在高性能船舶上得到越来越多的应用.随着计算流体动力学（CFD）的飞速发展，CFD技术在船舶领域发挥着越来越重要的作用.在喷水推进器方面，CFD技术广泛运用于性能研究和优化设计中，现已取得相当进展.在国外，Bulten[1]对喷水推进器进行了数值研究，分析了定常多重参考系（MRF）模型和瞬时刚体运动模型对喷水推进器性能的影响.Takai等[2]对喷水推进船流道形状的优化进行了研究，采用体积力法来模拟泵的效果，研究表明对流道上部曲线和唇部形状进行优化可以减小压力损失，对流道进口进行优化可以提高入口效率.Peri等[3]采用优化软件对高速喷水推进双体船进行了优化设计，以减小入口损失和减小阻力为目标，分别对喷水推进流道形状和船体型线进行了优化，最终优化后阻力减小了6%.Guo等[4]对喷水推进器的进口流动进行了研究，分析了不同进速比下进口的压力分布，研究表明进速比越大进口压力损失越大.Eslamdoost等[5－6]提出了一种压力跳跃的方法来代替泵的作用，对喷水推进器与船体的相互作用进行了研究，分析了喷水推进负推力减额的原因.Lu等[7]对喷水推进泵进行了数值计算，研究了喷水推进泵的水动力性能.在国内，孙存楼等[8]进行了喷水推进船负推力减额机理的研究，认为中高速时黏性阻力降低是推力减额的主要原因.龚杰等[9]进行了喷水推进船的自航数值模拟，比较了虚拟盘法和重叠网格法对船内外流场的影响，研究表明两种方法的外流场计算结果基本相同，内流场采用重叠网格法更接近真实叶轮情况.易文彬等[10]采用试验和数值模拟相结合的方法分析了浸没式喷水推进的水动力特点，试验值和CFD预报值吻合较好.

喷水推进器推力性能是喷水推进系统性能的一项重要内容.试验时喷水推进器推力的测量有直接测量法和动量流量法，动量流量法通过进流面和喷口的动量差来计算推力.由于直接测量法所需的测量设备非常复杂，技术门槛较高，目前仅有国外少数几家单位能够直接测量；动量流量法不但可以选取满足流量要求的任意泵，而且船体与喷水推进之间无须复杂的水密结构，已成为国内外研究机构进行试验和数值计算时采用的主要方法.从第21至第24届国际拖曳水池会议（ITTC），经过喷水推进专家委员会的不断完善，最终推出了基于动量流量法的喷水推进模型试验的推荐指南[11]，形成了比较严谨的理论体系.

根据动量流量法，在CFD数值计算喷水推进器推力的过程中，流量与喷口的速度很容易获得，但进流面的的速度不易获得，主要是进流面的形状无法确定，因此如何准确获得进流面的形状成为准确计算喷水推进器推力的关键.Ding等[12]提出了自定义标量函数的方法来确定喷水推进系统的进流面，但此种方法需要用户在商业软件中编写后处理程序进行二次开发，其要求较高，不适用于大多数研究人员.本文提出一种基于流线法的进流面获取方法，该方法简单实用，适用于大多数的研究人员.使用MATLAB语言开发成程序，使过程大大简化.在此基础上，采用CFD软件STAR－CCM＋对喷水推进器的性能进行数值计算.首先进行了网格无关性验证，之后选取合适的网格，系统研究了喷水推进器在不同转数下的性能，为喷水推进器系统性能的提高和改善奠定了基础.

1 理论模型

1.1 基本方程

流体力学问题要遵守基本守恒定律，包括质量守恒、动量守恒和能量守恒.在船舶流体力学问题中，通常不考虑能量交换，故忽略能量守恒，则不可压缩流体对应的连续性方程和动量守恒方程分别为

式中：ui和uj为速度时均量；xi和xj为坐标；t为时间；ρ为流体密度；p为压力时均量；ν为流体运动黏性系数为雷诺应力.

由于式（2）是不封闭的，需引入湍流模型，所以本文计算中采用SSTk－ω湍流模型.该模型结合k－ω和k－ε模型，在预报具有逆压梯度的流动和流动分离等方面具有优势，其表达形式为

式（3）和（4）中各参数的具体定义以及含义可参见文献[13].

1.2 动量流量法

根据ITTC规程，动量流量法求解喷水推进器推力的流场控制体如图1所示.

图1 流场控制体Fig.1 Control volume of waterjet

进流面A1与来流方向垂直，位于流道斜坡与船体的切点A′之前，距离通常取1倍的流道直径，以避免流道入口的影响.通常认为收缩断面A7与喷口面积A6的直径相同，将喷口面积A6作为控制体的出口面.流场控制体在x方向的总推力为

流体作用在喷水推进器壁面的力，直接传递在船体上，称之为净推力.在计算过程中可以很方便地监测喷水推进器各壁面上的压力和黏性力，通过壁面积分求解出喷水推进器的净推力，即喷水推进器在x方向的合力为

式中：S为喷水推进器的所有壁面面积；px为壁面上的压应力；σx为壁面上的剪切应力.

本文基于流线法获取进流面的基本思想如下：

（1）泵进口面A3设置流线起点，沿上游计算流线，流线分布区域的半径应略大于流道半径，适当设置流线种子的数量以使流线基本布满控制体，这样计算的形状更加接近于进流面的形状.

（2）在进流面的位置A1处建立与来流方向相垂直的平面截面，平面截面的对象为流线，如图2所示.

（3）对进流面上的流线截面数据（图3中蓝点）进行平面点集的凸包计算，从而得到凸包点集最外层的边界点，然后采用孔为平等[14]提出的4次多项式对边界点进行拟合，最终得到进流面的形状曲线（图3中红线）：

式中：c0，c1，c2，c3，c4分别为4次多项式拟合系数.

采用MATLAB语言开发的进流面获取程序如图3所示，进流面获取的原则是进流面和喷口截面的流量相等.程序有6个参数：进流面位置船尾的斜升角；进流面形状中心的y坐标和z坐标；边界曲线的延伸系数，以确保进流面曲线的左右边界超过控制体；缩放系数，由于流线法确定的进流面完全在控制体内部，进流面的流量会略小于喷口流量，而缩放系数可以确保进流面的流量与喷口的流量相等.

图2 流线法Fig.2 Streamline method

图3 进流面获取程序Fig.3 Program of capture area

当从泵进口面A3处计算的流线足够多时，流线所组成的包络体即为流场控制体，流线在位置A1的截面即为理论上的进流面.进流面形状对流线数量（m）的依赖性见图4，以流道直径D对进流面形状进行无量纲处理，Y＝y／D表示无量纲的y坐标，Z＝z／D表示无量纲的z坐标.可见：当m＝500时，得到的进流面形状有微小偏差；当m＞1 000后，得到的进流面形状不再变化，说明当流线数大于1 000后，本文方法可准确得到进流面的形状，也可说明本文方法的合理性.

图4 进流面对流线的依赖性Fig.4 Dependence of capture area on streamlines

2 计算模型

计算域长度为30D，流道入口前为25D，流道入口后为5D，宽度为10D，深度为8D.

网格使用STAR－CCM＋软件生成，网格类型为Trim非结构网格.在叶轮和导叶的导边和随边处生成特征线，进行局部加密.在流道入口处流场变化剧烈，设置网格控制体进行局部加密.计算域总的网格数约为291万个，叶轮旋转区域约为94万个，其余静止区域约为197万个.计算域的网格见图5.

图6所示为喷水推进器进口流道和叶轮表面第一层网格节点与壁面的无量纲距离y＋的分布，喷水推进器流道的y＋为90～160，叶轮和导叶的y＋在90左右，整体上喷水推进器的y＋大部分在30～300.第一层网格节点布置在对数律成立的范围内，以满足壁面函数的要求.

图5 计算域的网格Fig.5 Grids of computational domain

图6 y＋分布云图Fig.6 y＋contour

计算域入口、两侧和底部设为速度进口；计算域出口和推进器喷口为压力出口；计算域顶部、流道壁面和导叶设为无滑移壁面；将计算域分为旋转区域和静止区域，旋转区域采用稳态多参考系（MFR）方法模拟泵的旋转，在旋转区域的速度根据旋转角速度给出，叶轮和导流帽设为静止壁面；旋转区域和静止区域的交界设为交界面.

考虑船体边界层的影响，入口的速度分布采用平板边界层速度分布：式中：v0为设计航速；z为距计算域顶部边界的距离；δ为船模边界层厚度；L为计算域入口至船首的距离；Re为雷诺数.

采用基于分离流的黏性求解器，压力与速度的耦合采用SIMPLE算法计算，对流项离散格式为二阶迎风格式.

3 计算结果及分析

对喷水推进器中的各水动力系数定义如下：

式中：CTg和CTn分别为喷水推进器的总推力系数和净推力系数；KQJ、KH、KQ和η分别为喷水推进器的流量系数、扬程系数、转矩系数和效率；QJ、H和Qshaft分别为流量、扬程和转轴上的转矩；n为转速；RIV和RNV分别为喷水推进器的进速比和喷速比；v3为泵前进口x方向的平均速度；v6为喷口x方向的平均速度；α和β为分别为喷水推进器进流面的动量影响系数和能量影响系数[15]；vx为x方向的速度.

3.1 网格无关性分析

喷水推进器叶轮、导叶和流道以及加密区域的网格大小都与基本尺寸相关，由基本尺寸的变化生成由密到疏的5套网格，分别编号为1～5，来进行网格无关性的分析，5套网格的数目见表1.

在设计转速n0下，进行网格无关性分析，将最密网格（1号网格）的推力和转矩计算结果分别记为CTn1和Qshaft1，其他网格下推力和转矩分别以该值为基准进行归一化处理，网格无关性计算结果见图7.

表1 网格的数目与尺寸Tab.1 Grid sizes and numbers of five grids

图7 网格数目依赖性Fig.7 Dependence of grid number

由图7可以看出，转矩对网格的依赖性不明显，推力对网格的依赖性较明显.当网格数增加至2.91×106时，推力和转矩都基本接近各自的参考基准值.综合考虑计算精度和计算资源，最终采用2号网格方案来进行喷水推进器性能的数值计算.

3.2 喷水推进器性能分析

在设计航速v0下，研究喷水推进器在不同转速（n／n0＝0.4，0.6，0.8，1.0，1.2，1.4，1.6，1.8，2.0）时的水动力性能.

喷水推进器的进流面动量影响系数α、进流面能量影响系数β、进速比RIV和喷速比RNV是影响喷水推进器性能的重要参数；α、RIV和RNV分别为表征A1、A3和A6平均速度的无量纲量.在不同转速下，α、β、RIV和RNV的结果见表2.

表2 影响喷水推进器性能的参数Tab.2 Parameters affecting the performance of waterjet

根据本文开发的程序可以很方便地得到进流面边界的曲线形状，图8所示为不同转速下进流面边界的曲线形状.由于进流面形状左右对称，参数的奇数项系数为0，表3所示为进流面形状在不同转速下的相关参数化表达.其中：w为进流面宽度；c0的绝对值为无量纲的进流面高度.在低转速时，进流面边界的曲线形状不是椭圆形，中间稍微内凹，说明在经验设计中，如果采用椭圆方程计算的进流面和流量会产生较大偏差；随着转速的增加，中间形状由内凹向直线再向外凸过渡，进流面的边界形状逐渐趋向椭圆形；随着转速的增加，进流面的高度和宽度都增加，高度的变化大于宽度的变化.在进行喷水推进试验时，本文计算的进流面形状参数可为获取进流面提供依据.

图8 进流面形状Fig.8 Shapes of capture area

表3 进流面的形状参数Tab.3 Shape parameters of capture area

图9所示为不同转速下进流面A1的速度分布，v为流体的速度.在低于设计转速时，进流几乎全部取自边界层，随着转速的增加流量也相应增加，进流的高度逐渐大于边界层厚度，进流有相当一部分取自边界层以外的流场.

图10所示为不同转速下泵前进口处A3面上速度分布云图.可以看出：泵前进口处速度分布并不均匀，下部速度大，上部速度小；在低转速时，下部速度较均匀，上部速度梯度变化大；在转速较高时，则上部速度较均匀，下部速度梯度变化大；受到泵旋转的影响，速度向旋转方向偏移.

图11所示为不同转速下喷口处A6面上速度分布云图.可以看出：喷口处速度受到泵旋转的影响，分布并不均匀；在喷口中心处有速度较低的区域，低速区的直径大约为喷口直径的20%，低速区的速度大约为喷口最大速度的50%，在低转速时，低速区稍向下偏离喷口的几何中心.

喷水推进器唇部的流场变化剧烈，会产生低压区，容易发生空化现象，图12所示为不同转速下唇部表面压力分布云图和流线图.定义无量纲的压力系数为

式中：p为绝对压力；p0为参考压力（取1个标准大气压）.

由图12可以发现，在唇部发生了流动分流，而发生流动分流的位置随着转速的增加向下移动，这是因为随着转速的增加，流量增加，需要进口从更深的位置抽吸水流.低压区域与流动分流点的变化呈相反趋势，低转速时，低压区域在唇部的下缘；高转速时，低压区域在唇部的上缘.

图9 进流面的速度分布云图Fig.9 Velocity contour of capture area

图10 泵前进口的速度分布云图Fig.10 Velocity contour of pump intake

图11 泵前进口的速度分布云图Fig.11 Velocity contour of pump intake

图13所示为推力系数曲线，总推力系数曲线和净推力系数曲线基本重合，说明ITTC推荐的动量流量法和壁面积分法都可以进行喷水推进器的推力预报，也间接说明本文对进流面获取的准确性，动量流量计算的推力精度较高.在喷水推进计算中，如果采用压力跳越法[5]或其他方法代替泵的作用，壁面积分法就不再适用，只能采用动量流量法来计算推力，因此动量流量法的适用性更广，此时采用本文开发的进流面获取程序将会大大提高效率.

图14所示为喷水推进器的性能曲线，喷水推进当中的流量－扬程曲线在物理含义上相当于敞水螺旋桨中的进速－推力曲线.可以看出：随着流量系数KQJ的增加，扬程系数KH和转矩系数KQ不断下降；效率η在流量系数KQJ＝0.87时达到最大，在最高效率之前效率变化缓慢，此时转速比n／n0＝0.8～2.0，说明喷水推进器在较大的转速范围内都有较高的效率；在最高效率之后效率η下降较快，主要原因是低转速时扬程系数KH下降较快.

图12 唇部的压力系数云图与流线Fig.12 Pressure coefficient contour and streamline in lip

图13 总推力和净推力曲线Fig.13 Gross thrust curve and net thrust curve

图14 喷水推进器性能曲线Fig.14 Performance curves of waterjet

4 结语

本文对喷水推进器进流面的获取方法进行研究，采用基于流线法的进流面获取方法可以准确得到进流面，该方法操作简单方便.开发的进流面获取程序使后处理过程大大简化，提高了效率.在此基础上，对喷水推进器在不同转速下的性能进行了数值预报和分析，给出了进流面形状的参数化表达，为喷水推进试验提供指导.在低转速时，进流面形状不是椭圆形，随着转速增加，进流面形状逐渐趋向椭圆形.进流面形状、泵前进口面速度分布以及唇部流动分离点位置都随转速变化而变化.采用动量流量法计算喷水推进器的推力，不仅精度高，而且应用更广.

致谢 本研究得到了喷水推进技术重点实验室的支持，得益于喷水推进技术学术研讨会（WORKSHOP2018）多位与会专家的技术交流，在此深表谢意.