一种关于水下螺旋桨梢涡空化的数值预报方法

摘要

本发明公开了一种关于水下螺旋桨梢涡空化的数值预报方法，属于螺旋桨优化设计领域。一种关于水下螺旋桨梢涡空化的数值预报方法，包括以下步骤：(a)基本网格确定；(b)精细网格确定；(c)最优网格确定；(d)利用最优网格预报所需工况条件下的梢涡空化。本发明针对螺旋桨梢涡空化进行数值预报；通过与相关实验结果对比，本发明较为有效地预报了E779A型螺旋桨在几种不同工况下的梢涡空化，因此，本发明对螺旋桨设计中空化性能的预测与评估具有重要作用，可有效减少设计成本和设计周期，其应用前景非常好。

说明书

技术领域

本发明涉及螺旋桨优化设计领域，具体地说，一种关于水下螺旋桨梢涡空化 的数值预报方法。

背景技术

螺旋桨空化不仅会降低螺旋桨性能,产生空化剥蚀,导致船体振动,而且 会产生足以暴露自己的空化噪声。随着现代船舶对负荷和船速的要求越来越高， 螺旋桨空化现象难以避免。而梢涡空化常常是螺旋桨最早出现的空化，其现象的 准确预报可以对空化初生的判断提供重要依据。梢涡空化的数值预报也是最难实 现的空化数值预报之一。数值捕获梢涡空泡，并获取其周围流场的压力、速度和 流场结构在内的全面空化流场信息可以进一步增强对噪声产生机理的认识。梢涡 空化一直是众多相关学者的的研究热点和难点，也是推进器设计制造乃至造船行 业所关注的热点之一。

目前，对空化流场的研究主要有实验和数值分析两种方法。完全依赖实验存 在耗资费时、螺旋桨模型修改不方便以及传感器在流场中的安装会对流场造成影 响等不足。近年来，随着计算机技术的快速发展，计算机运算能力的大大提高， 使得利用现代计算流体力学对梢涡空化进行数值预报成为可能，数值预报成为空 化预报的重要手段。现在对螺旋桨空化数值模拟较为成功是片空化，而梢涡空化 (tip vortex cavitation)数值预报难度较大。目前常用的预报方法为，通过确立合 理的计算域空间，建立高质量的网格，选用合适的湍流模型和空化模型，在基于 混合多相流模型的基础上求解雷诺应力RANS方程得到螺旋桨周围流场的数值 解，已经能基本实现螺旋桨空化的数值预报，但是梢涡空化的预报效果不明显。

目前，国外对梢涡空化研究主要侧重于水翼。美国Dynaflow Inc.公司将修正 的球形空泡动力学模型(Rayleigh–Plesset)和非球形空泡动力学模型相结合，并内 置到UnRANS方程求解器中进行求解，对梢涡空化初生现象进行数值预报，他 们考虑到流场中空化核的尺寸和空间分布对空化初生的影响。韩国首尔大学研究 人员采用Eulerian–Lagrangian方法，并同样考虑空泡核的尺寸分布特征对梢涡空 化及其空化噪声进行了数值分析。挪威研究人员利用商业软件Fluent软件对 NACA水翼的梢涡空化进行了数值模拟，其结果与实验结果较为一致。此外，瑞 典研究人员利用大涡模拟方法对螺旋桨空化进行了数值预报，其螺旋桨梢涡空化 预报结果与实验现象较为一致，但片空化预报区域明显比实验显示结果大。

所以，急需一种能够有效的预报螺旋桨梢涡空化的数值预报方法。

发明内容

螺旋桨梢涡空化主要特征有：

1.在螺旋桨各种空化形式中，梢涡空化常常最早出现，并伴随着空化初生现 象；

2.在尾流场中螺旋桨叶梢附近产生，沿轴向呈现螺旋线形状分布，且痕迹较 长；

3.在梢涡中心线区域有强烈的涡旋运动，并产生强烈的低压；

4.梢涡最初在叶梢附近产生，在脱离桨叶后，梢涡的螺旋线半径会进行收缩；

5.梢涡会随着螺旋桨尾流的涡面结构向尾流方向进行滑溜；

6.梢涡空化发出的噪声一般分布在中高频段，其频谱呈现连续特征。

本发明的原理就是依据上述的梢涡空化的主要特征，综合使用空化模型、湍 流模型和混合两相流模型来数值求解纳维—斯托克斯(N-S)方程组，得到螺旋 桨流体环境中的压力、速度和反映空化的汽相体积分数等各项物理参数值，本发 明提供一种关于水下螺旋桨梢涡空化的数值预报方法，可对螺旋桨设计中空化性 能进行预测与评估。

本发明采用如下技术方案：

一种关于水下螺旋桨梢涡空化的数值预报方法，包括以下步骤：

(a)基本网格确定；

(b)精细网格确定；

(c)最优网格确定；

(d)利用最优网格预报所需工况条件下的梢涡空化。

本发明方法的具体实现方式包括以下步骤：

步骤一，利用建模软件建立螺旋桨三维几何模型，并将其导入网格划分软件。

步骤二，在网格划分软件中建立三种备选网格，并将其导入计算程序：

三种备选网格的计算域相同，计算域为圆柱形，其速度入流边界距离螺旋桨 中心为1.5D，D为螺旋桨直径，下游压力出口边界距离为5D，螺旋桨中心至 侧面外围距离为2.5D，三种备选网格的网格数量逐渐增加，分别约为200万、 300万和400万。

步骤三，空化模型和湍流模型设定：

采用全空化模型和重整化群湍流模型(RNG k-ε湍流模型)，并对其重要参数 进行修正，对空化模型中相变率参数的修正和湍流模型中湍流黏度系数的修正采 用C语言编写，再利用宏调用(DEFINE_TURBULENT_VISCOSITY等)形式嵌 入计算程序。

步骤四，数值计算参数设定：

数值参数设定包括工况条件、边界条件和数值算法的相关参数设定；工况条 件主要设定螺旋桨旋转速度，环境压力和入流速度值，从而确定螺旋桨无量纲参 数，即进速系数(J)和空化数(σ_n)；对于边界条件设定，速度入口边界采用入 流速度值，远场边界条件同样采用入流速度值设定，下游压力出口界面的出口压 力设置为静压力；对于数值算法，纳维—斯托克斯(N-S)方程中对流项采用二 阶迎风格式离散，扩散项采用二阶中心差分格式离散，速度压力耦合采用适合非 结构网格的SIMPLE算法，使用逐点Gauss-Seidel迭代求解离散方程，利用代数 多重网格加速计算收敛，对于非定常计算采用滑动网格计算技术，提高计算的准 确性。由于多相流模型、空化模型和滑动网格计算对计算机资源消耗较大，采用 并行计算技术来缩短计算时间；

质量守恒连续性(continuity)残差收敛标准为三阶，方程中其它物理量残差 收敛标准为四阶，利用单相流的收敛解作为多项流求解的初始值，并将稳态解作 为非稳态计算的初始值，为了确保计算收敛，螺旋桨旋转速度逐步加大到预定值， 并适当缩小松弛因子。

Navier-Stokes方程，也称纳维—斯托克斯方程，简称N-S方程。下面表达式 为张量形式的N-S方程，其中下标k与j,i一样表示坐标轴。

$\left(\begin{array}{c}\frac{\partial }{\partial t}\left(\rho u_i\right)+\frac{\partial }{\partial x_i}\left(\rho u_i{u}_j\right)=\rho f_i-\frac{\partial p}{\partial x_i}+\frac{\partial }{\partial x_j}\mu [(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i})-\frac{2}{3}{\delta }_{\mathrm{ij}}\frac{\partial u_k}{\partial x_k}\\ =\rho f_i-\frac{\partial p}{\partial x_i}+\frac{\partial }{\partial x_i}\left(\mu \frac{\partial u_i}{\partial x_j}\right)+\frac{1}{3}\frac{\partial }{\partial x_j}\left(\mu \frac{\partial u_k}{\partial x_k}\right)\end{array}\right)$

式中，μ为流体的动力黏度，δ_ij为“Kronecker Delta”张量。

步骤五，进行数值计算：

由于空化模型加入N-S方程后，计算的稳定性降低，容易出现奇异现象。 因此，为了能使数值计算平稳进行，采用逐级分步骤的计算过程。具体来说，在 螺旋桨工况参数中，环境压力和入流速度可以直接设定到工况值，而螺旋桨转速 采用分级增加，直到增加到预定工况值。另外，先计算无空化模型流场分布，等 到计算稳定后再打开空化模型。先对压力、密度、动量和汽相分数等参数进行一 阶精度离散格式计算，计算稳定后，再将离散精度提高到二阶或QUCIK等。为 了保证二阶计算的稳定性，将亚松弛因子适当降低。压力、动量、汽相分数、湍 流动能、湍流耗散率和湍流黏性等参数的亚松弛因子分别设定为：0.25、0.6、0.2、 0.7、0.7、0.9。

步骤六，根据水动力数值结果确立基本网格：

利用这三种备选网格对螺旋桨水动力进行数值计算，并对计算结果进行分 析；当水动力参数(推力系数和转矩系数)的计算结果随着网格数增加后趋于平 稳后，则认为数值计算随着网格增加而趋于稳定，则选定计算稳定后网格数量最 少的网格作为基本网格。

步骤七，将基本网格的文件导入网格划分软件。

步骤八，在基本网格的基础上，根据梢涡形状模型，利用网格划分软件建立 梢涡区域网格：

首先，根据等距螺旋线模型在专业建模工具中建立一个与模型参数相符的螺 旋线形状的管道区域几何模型，这个区域是假定的梢涡空化发生区域，再将螺旋 线几何形状模型导入基本网格中，并进行网格划分，等距螺旋线形状区域网格从 临近桨叶叶梢开始，等距螺旋线半径预设为0.82R，螺旋线区域的单元网格尺寸 约为0.0001D，原基本网格中每条边的单元尺寸不变。

步骤九，重复步骤三的空化模型和湍流模型设定，步骤四的数值计算参数设 定和步骤五的进行数值计算。

步骤十，根据梢涡空化数值预报结果确立精细网格：

这一步骤主要为确立等距螺旋线模型中的重要参数r。螺旋线半径r的长度 约为(0.81-0.83)R，根据梢涡空化数值结果来调整并最终确定这一参数，具体方 法为，当步骤九的计算结果中梢涡空化基本位于这一假定区域(即为步骤八中设 定的r所建立的螺旋线管道区域就是梢涡空化发生区域)，则精细网格确立，可 以进行下一步骤；反之，当梢涡空化数值预报基本不在这一假定区域时，则根据 位置误差对参数r进行调整，并返回步骤八和步骤九重新计算，直到结果基本相 符，则确立精细网格。

步骤十一，将精细网格文件导入网格划分软件。

步骤十二，建立尾流距离分别为7D和9D的另外两种计算域网格：

在精细网格的基础上，建立另外两种计算域，其尾流距离分别为7D和9D， 原来精细网格中每条边的单元尺寸不变，划分增加的尾流部分区域网格，其网格 单元尺寸为原精细网格在5D处网格单元尺寸。

步骤十三，重复步骤三的空化模型和湍流模型设定，步骤四的数值计算参数 设定和步骤五的进行数值计算。

步骤十四，与尾流距离为5D的结果进行比较确定最优网格。

对这三种计算域网格进行数值计算后进行分析，当随着尾流距离增加，梢涡 空化数值预报精度不再明显增加时，则确立最佳尾流距离及最优网格。

步骤十五，利用最优网格和数值计算参数对所需工况进行数值预报。

更进一步地，所述的全空化模型及其参数修正为：

当p<p_v时，蒸汽产生率为：

$R_e=\mathrm{Ce}\frac{k}{\gamma }{\rho }_l{\rho }_v\sqrt{\frac{2}{3}\frac{p_v-p}{{\rho }_l}}(1-f_v)$

当p>p_v时，汽相变液相，蒸汽凝固率R_c为：

$R_c=\mathrm{Cc}\frac{k}{\gamma }{\rho }_l{\rho }_v\sqrt{\frac{2}{3}\frac{p-p_v}{{\rho }_l}}{f}_v$

其中，f_v＝α_vρ_v/ρ_m为汽相质量分数，汽化系数Ce＝0.02和凝结系数Cc＝0.01 为经验参数。

根据量纲分析在相变率(Re和Rc)表达式中采用k而不是。在FLUENT 软件环境下可利用自定义函数UDF对空化模型中参数相变率Re进行修正，修正 采用C语言编写后调入计算程序。

更进一步地，所述的重整化群湍流模型及其参数修正为：

重整化群湍流模型RNG k-ε是对瞬态N-S方程用重整化群(Renormalization  Group，简称RNG)的数学方法推导出来的模型。它通过在大尺度运动项和修正 黏度项中体现小尺度的影响，而使这些小尺度运动系统地从控制方程中除去。其 k方程和ε方程分别为：

$\frac{\partial }{\partial t}\left({\rho }_{m}k\right)+\frac{\partial }{\partial x_j}\left({\rho }_{m}k{u}_{\mathrm{mj}}\right)=\frac{\partial }{\partial x_j}\left[\left({\alpha }_k\mu \right)\frac{\partial k}{\partial x_j}\right]+G-{\rho }_mϵ$

$\frac{\partial }{\partial t}\left({\rho }_mϵ\right)+\frac{\partial }{\partial x_j}\left({\rho }_mϵu_{\mathrm{mj}}\right)=\frac{\partial }{\partial x_j}\left[\left({\alpha }_{ϵ}\mu \right)\frac{\partial ϵ}{\partial x_j}\right]-R+C_{1ϵ}\frac{ϵ}{k}G-C_{2ϵ}{\rho }_m\frac{{ϵ}^2}{k}$

式中，湍流动能耗散率(Turbulent Dissipation Rate)湍流动能k和耗散率ε的有效湍流普朗特数的倒数a_k＝a_ε＝1.39；模型参数C _1ε＝1.47，C_2ε＝1.68；黏性系数为μ＝μ_t+μ_m，μ_m为混合流黏度系数；修改了湍流黏 度系数μ_t＝[ρ_v+α_l¹⁰(ρ_l-ρ_v)]C_μk²/ε，C_μ＝0.085，使其更适用于高雷诺数的非定常两 相流模拟，从而可以更好地模拟螺旋桨空化。

更进一步地，所述的等距螺旋线数学模型及其参数如下：

$\left(\begin{array}{c}x=r\cos \left(\frac{2\pi }{k}z\right)\\ y=r\sin \left(\frac{2\pi }{k}z\right)\end{array}\right)$

其中，x、y、z分别为笛卡尔坐标系的三个坐标轴，常数r为螺旋线半径，r 的长度约为(0.81-0.83)R，R为螺旋桨半径；常数k为螺旋线的每旋转一周在x轴 上前进的距离，k＝U_∞/n，U_∞为入流速度，n为螺旋桨转速。

更进一步地，所述的网格划分的方法为：采用分区域混合网格划分方法划分 网格：螺旋桨周围流场区域采用非结构网格方法划分，网格由桨毂到叶梢逐渐减 小，叶梢处面网格为三角形，大小约为0.001D，桨榖处为0.015D，D为螺旋桨 直径；在桨叶表面采用边界层网格，边界层网格共有4层，其相邻两层高度比为 1.1，第一层网格单元高度约为0.0008D，使得无量纲参数20<y⁺<300，采用结构 网格划分螺旋桨外围规则形状的计算域。

本发明的方法可以在一般通用的CFD流体计算软件(CFX，FLUENT等) 中实现，网格划分可以采用GAMBIT等软件实现。首先将初步设计的螺旋桨数 字模型导入网格划分软件，并按照本发明中的方法进行网格划分。网格模型导入 计算平台，按照设计工况进行数值计算，并将数值结果返回螺旋桨设计人员，进 行空化性能评估，再确定所设计的螺旋桨是否达到设计要求。此外，本发明还利 用批处理文件在操作系统平台进行并行数值计算。

有益效果：

(1)梢涡空化预报中影响梢涡空化数值计算结果的主要技术环节有：模型 (湍流模型和空化模型)、数值算法和网格，由于模型和数值算法对不同对象(桨 模和工况)的数值计算有较强的适应性和鲁棒性，而网格对不同对象的数值计算 有较大差异性。本发明通过将模型，算法和网格等重要条件进行组合，找到一个 有效的技术途径来更好地预报梢涡空化，本发明对螺旋桨设计中空化性能的改善 具有重要应用价值和应用前景。

(2)本发明基于黏性流理论，综合使用空化模型、湍流模型和混合两相流 模型来数值求解纳维—斯托克斯(N-S)方程组，得到螺旋桨流体环境中的压力、 速度和反映空化的汽相体积分数等各项物理参数值，本发明技术是基于近来发展 迅速的现代计算流体力学和空泡动力学，将网格生成技术、数值计算方法和并行 计算技术引入螺旋桨梢涡空化的数值预报，体现了多学科与多领域技术的融合 性。

(3)由于水翼在流体做直线平移运动，而螺旋桨则是在流体中做旋转运动， 其周围流场情况比水翼复杂。本发明针对螺旋桨梢涡空化进行数值预报；通过与 相关实验结果对比，本发明较为有效地预报了E779A型螺旋桨在几种不同工况下 的梢涡空化，因此，本发明对螺旋桨设计中空化性能的预测与评估具有重要作用， 可有效减少设计成本和设计周期，其应用前景非常好。

附图说明

图1为本发明的全流道螺旋桨计算域的示意图；

图2(a)为桨叶表面示意图；

图2(b)为图2(a)中A处的边界层网格示意图；

图3为叶梢附近网格分布的示意图；

图4为E779A螺旋桨几何模型；

图5(a)为本发明的基本网格确定流程图；

图5(b)为本发明的精细网格确定流程图；

图5(c)为本发明的最优网格确立流程图；

图5(d)为本发明的主流程图；

图6为桨叶尾流中梢涡在叶梢附近似呈现螺旋线形状分布示意图

图7为本发明的梢涡空化数值预报与实验结果的对比。

图8(a)为0°位置E779A桨模非均匀入流条件下在1/8周期内的梢涡空化 数值模拟结果；

图8(b)为18°位置E779A桨模非均匀入流条件下在1/8周期内的梢涡空化 数值模拟结果；

图8(c)为36°位置E779A桨模非均匀入流条件下在1/8周期内的梢涡空化 数值模拟结果；

图8(d)为54°位置E779A桨模非均匀入流条件下在1/8周期内的梢涡空化 数值模拟结果；

图8(e)为72°位置E779A桨模非均匀入流条件下在1/8周期内的梢涡空化 数值模拟结果；

图8(f)为81°位置E779A桨模非均匀入流条件下在1/8周期内的梢涡空化 数值模拟结果。

图中：1、速度入流边界，2、远场边界，3、压力出口边界。

具体实施方式

下面结合附图和实施例对本发明作进一步详细的说明。

实施例

如图5(a)、图5(b)、图5(c)和图5(d)所示，一种关于水下螺旋桨梢涡空化 的数值预报方法，具体实现方式包括以下步骤：

步骤一，利用专业建模软件建立螺旋桨三维几何模型，并将其导入网格划分 软件。

步骤二，在网格划分软件中建立三种备选网格，并将其导入计算程序。

三种备选网格的计算域相同，根据螺旋桨周围流场特征，建立圆柱形计算域， 如图1所示，其速度入流边界距离螺旋桨中心为1.5D，D为螺旋桨直径，下游 压力出口边界距离为5D，螺旋桨中心至侧面外围距离为2.5D，三种备选网格的 网格数量逐渐增加，分别约为200万、300万和400万；网格划分的方法为：采 用分区域混合网格划分方法划分网格：螺旋桨周围流场区域采用非结构网格方法 划分，网格由桨毂到叶梢逐渐减小，叶梢处面网格(如图3所示)为三角形，大 小约为0.001D，桨榖处为0.015D，D为螺旋桨直径；由于空化主要分布在叶面， 因此这一区域网格质量要求较高，为了更好地适应壁面函数，在桨叶表面采用边 界层网格，采用边界层网格以提高对梢涡空化的预报精度，如图2(a)所示， 桨叶表面示意图，图2(b)为图2(a)中A处的边界层网格示意图，桨叶表面 边界层网格共有4层，其相邻两层高度比为1.1，第一层网格单元高度约为 0.0008D，使得无量纲参数20<y⁺<300，采用结构网格划分螺旋桨外围规则形状 的计算域。

步骤三，空化模型和湍流模型设定：

采用全空化模型和重整化群湍流模型(RNG k-ε湍流模型)，并对其重要参数 进行修正，对空化模型中相变率参数的修正和湍流模型中湍流黏度系数的修正采 用C语言编写，再利用宏调用(DEFINE_TURBULENT_VISCOSITY等)形式嵌 入计算程序；

全空化模型及其参数修正为：

当p<p_v时，蒸汽产生率为：

$R_e=\mathrm{Ce}\frac{k}{\gamma }{\rho }_l{\rho }_v\sqrt{\frac{2}{3}\frac{p_v-p}{{\rho }_l}}(1-f_v)$

当p>pv时，汽相变液相，蒸汽凝固率Rc为：

$R_c=\mathrm{Cc}\frac{k}{\gamma }{\rho }_l{\rho }_v\sqrt{\frac{2}{3}\frac{p-p_v}{{\rho }_l}}{f}_v$

其中，f_v＝α_vρ_v/ρ_m为汽相质量分数，汽化系数Ce＝0.02和凝结系数Cc＝0.01 为经验参数。

根据量纲分析在相变率(Re和Rc)表达式中采用k而不是在FLUENT 软件环境下可利用自定义函数UDF对空化模型中参数相变率Re进行修正，修正 采用C语言编写后调入计算程序；

重整化群湍流模型RNG k-ε是对瞬态N-S方程用重整化群(Renormalization  Group，简称RNG)的数学方法推导出来的模型。它通过在大尺度运动项和修正 黏度项中体现小尺度的影响，而使这些小尺度运动系统地从控制方程中除去。其 k方程和ε方程分别为：

$\frac{\partial }{\partial t}\left({\rho }_{m}k\right)+\frac{\partial }{\partial x_j}\left({\rho }_{m}k{u}_{\mathrm{mj}}\right)=\frac{\partial }{\partial x_j}\left[\left({\alpha }_k\mu \right)\frac{\partial k}{\partial x_j}\right]+G-{\rho }_mϵ$

$\frac{\partial }{\partial t}\left({\rho }_mϵ\right)+\frac{\partial }{\partial x_j}\left({\rho }_mϵu_{\mathrm{mj}}\right)=\frac{\partial }{\partial x_j}\left[\left({\alpha }_{ϵ}\mu \right)\frac{\partial ϵ}{\partial x_j}\right]-R+C_{1ϵ}\frac{ϵ}{k}G-C_{2ϵ}{\rho }_m\frac{{ϵ}^2}{k}$

式中，湍流动能耗散率(Turbulent Dissipation Rate)湍流动能k和耗散率ε的有效湍流普朗特数的倒数a_k＝a_ε＝1.39；模型参数C _1ε＝1.47，C_2ε＝1.68；黏性系数为μ＝μ_t+μ_m，μ_m为混合流黏度系数；修改了湍流黏 度系数μ_t＝[ρ_v+α_l¹⁰(ρ_l-ρ_v)]C_μk²/ε，C_μ＝0.085，使其更适用于高雷诺数的非定常两 相流模拟，从而可以更好地模拟螺旋桨空化。

步骤四，数值计算参数设定：

数值参数设定包括工况条件、边界条件和数值算法的相关参数设定；工况条 件主要设定螺旋桨旋转速度，环境压力和入流速度值，从而确定螺旋桨无量纲参 数，即进速系数(J)和空化数(σ_n)；对于边界条件设定，速度入口边界采用入 流速度值，远场边界条件同样采用入流速度设定，下游压力出口界面的出口压力 设置为静压力；对于数值算法，纳维—斯托克斯(N-S)中对流项采用二阶迎风 格式离散，扩散项采用二阶中心差分格式离散，速度压力耦合采用适合非结构网 格的SIMPLE算法，使用逐点Gauss-Seidel迭代求解离散方程，利用代数多重网 格加速计算收敛，对于非定常计算采用滑动网格计算技术，提高计算的准确性。 由于多相流模型、空化模型和滑动网格计算对计算机资源消耗较大，采用并行计 算技术来缩短计算时间；

质量守恒连续性(continuity)残差收敛标准为三阶，方程中其它物理量残差 收敛标准为四阶，利用单相流的收敛解作为多项流求解的初始值，并将稳态解作 为非稳态计算的初始值，为了确保计算收敛，螺旋桨旋转速度逐步加大到预定值， 并适当缩小松弛因子；

步骤五，进行数值计算：

由于空化模型加入N-S方程后，计算的稳定性降低，容易出现奇异现象。 因此，为了能使数值计算平稳进行，采用逐级分步骤的计算过程。具体来说，在 螺旋桨工况参数中，环境压力和入流速度可以直接设定到工况值，而螺旋桨转速 采用分级增加，直到增加到预定工况值。另外，先计算无空化模型流场分布，等 到计算稳定后再打开空化模型。先对压力、密度、动量和汽相分数等参数进行一 阶精度离散格式计算，计算稳定后，再将离散精度提高到二阶或QUCIK等。为 了保证二阶计算的稳定性，将亚松弛因子适当降低。压力、动量、汽相分数、湍 流动能、湍流耗散率和湍流黏性等参数的亚松弛因子分别设定为：0.25、0.6、0.2、 0.7、0.7、0.9。

步骤六，根据水动力数值结果确立基本网格：

利用这三种备选网格对螺旋桨水动力进行数值计算，并对计算结果进行分 析；当水动力参数(推力系数和转矩系数)的计算结果随着网格数增加后趋于平 稳后，则认为数值计算随着网格增加而趋于稳定，则选定计算稳定后网格数量最 少的网格作为基本网格。

步骤七，将基本网格的文件导入网格划分软件。

步骤八，在基本网格的基础上，根据梢涡形状模型，利用网格划分软件建立 梢涡区域网格：

这里关键技术问题是如何合理地建立一个反映螺旋桨梢涡空化区域的几何 模型，这也是本发明方法的重要环节之一。具体来说，根据梢涡空化在尾流中分 布近似地呈现螺旋线形状特征，在螺旋桨尾流区域建立等距螺旋线形状的计算 域，即以螺旋线为中心划分一个狭长的计算域空间。等距螺旋线形状区域网格从 临近桨叶叶梢开始，并向尾流方向延伸。等距螺旋线的参数与螺旋桨几何形状参 数和工况参数相关。具体的等距螺旋线数学模型及其参数如下：

$\left(\begin{array}{c}x=r\cos \left(\frac{2\pi }{k}z\right)\\ y=r\sin \left(\frac{2\pi }{k}z\right)\end{array}\right)$

其中，x、y、z分别为笛卡尔坐标系的三个坐标轴，常数r为螺旋线半径，r 的长度约为(0.81-0.83)R，R为螺旋桨半径；常数k为螺旋线的每旋转一周在x轴 上前进的距离，k＝U_∞/n，U_∞为入流速度，n为螺旋桨转速；图6为国外某参考 文献中提供的桨叶尾流中梢涡在叶梢附近似呈现螺旋线形状分布示意图。

梢涡空化区域网格制作具体方法为：首先，根据上述等距螺旋线模型在专业 建模工具中建立一个与模型参数相符的螺旋线形状的管道区域几何模型，这个区 域是假定的梢涡空化发生区域，再将螺旋线几何形状模型导入基本网格中，并进 行网格划分，等距螺旋线形状区域网格从临近桨叶叶梢开始，等距螺旋线半径预 设为0.82R，螺旋线区域的单元网格尺寸约为0.0001D，原基本网格中每条边的 单元尺寸不变。

步骤九，重复步骤三的空化模型和湍流模型设定，步骤四的数值计算参数设 定和步骤五的进行数值计算。

步骤十，根据梢涡空化数值预报结果确立精细网格：

精细网格指在基本网格的基础上，对梢涡空化发生区域进行精细网格划分， 主要是提高这一区域的网格密度。

在这一步骤中主要确立等距螺旋线模型中的重要参数r。螺旋线半径r的长 度约为(0.81-0.83)R，这里并没有很好方法一次性确立这一参数，只有根据梢涡 空化数值结果来调整并最终确定这一参数，具体来说先采用步骤八中所建立网格 (按照r为0.82R来建立的一个假定的梢涡空化区域网格)，并经过步骤九的数 值计算。当计算结果中梢涡空化基本位于这一假定区域(即r为0.82R梢涡区域)， 则精细网格确立，可以进行下一步骤。反之，当梢涡空化数值预报基本不在这一 假定区域时，则根据位置误差对参数r进行调整，并重新计算，直到结果基本相 符。

步骤十一，将精细网格文件导入网格划分软件。

步骤十二，建立尾流距离分别为7D和9D的另外两种计算域网格：

在精细网格的基础上，建立另外两种计算域，其尾流距离分别为7D和9D， 原来精细网格中每条边的单元尺寸不变，划分增加的尾流部分区域网格，其网格 单元尺寸为原精细网格在5D处网格单元尺寸；

步骤十三，重复步骤三的空化模型和湍流模型设定，步骤四的数值计算参数 设定和步骤五的进行数值计算。

步骤十四，与尾流距离为5D的结果进行比较确定最优网格：

螺旋桨中心至计算域的压力出口距离对梢涡空化的准确预报非常重要。距离 大了浪费资源，距离小了会引起数值计算过程中尾流耗散增加，从而降低数值预 报中的梢涡强度，使得梢涡空化在数值模拟中不能很好预报。而且这一参数会随 着螺旋桨形状和工况参数的不同而发生变化。因此，在本发明方法中需要对这一 参数进行测定。这里分别建立三种计算域，其尾流距离分别为5D、7D和9D。 其中第一种计算域(5D)在前面步骤中已经完成，后面两种计算域只需增加尾 流部分区域的网格。对这三种计算域网格进行数值计算后进行分析。当随着尾流 距离增加，梢涡空化数值预报精度不再明显增加时，则确立最佳尾流距离及最优 网格。

步骤十五，利用最优网格和数值计算参数对其它工况进行数值预报。

本发明较为有效地预报了E779A型螺旋桨(模型见图4)在几种不同工况下的 梢涡空化。图7显示了两种工况条件下，E779A螺旋桨梢涡空化数值模拟结果与 其它文献中报道的实验结果。图中采用汽相体积分数av＝0.2的等值面(图中黑 色区域)作为空化数值结果。图7数值结果中梢涡空化尾迹非常明显，并与实验 结果吻合。说明本发明方法对梢涡空化数值预报效果显著。图8为非均匀入流条 件下近1/4周期内的梢涡空化数值模拟结果。图8(a)(b)(c)(d)(e)(f)中 黑色区域分别为逆时针旋转的E779A桨模的1号桨叶在0°，18°，36°，54°， 72°和81°等位置时刻的空化预报。图中显示梢涡空化尾迹非常明显，其长度 和粗细随着桨叶旋转位置变化而变化，表明本发明中使用方法对非均匀入流条件 下梢涡空化的数值模拟同样有显著效果。