CN103530482B - 一种非均匀入流中螺旋桨噪声数值预报方法 - Google Patents

Abstract

本发明公开了一种非均匀入流中螺旋桨噪声数值预报方法，步骤包括，首先，划分CFD计算网格，进行稳态迭代计算淌水性能参数和入流口速度验证模型的准确性；接下来，进行CFD非稳态计算，记录螺旋桨片空化周期形态以及片空化面积变化；然后，通过CFD流场数据，进行声学边界元数值计算，计算螺旋桨负载噪声；最后，通过螺旋桨片空化单极子辐射噪声以及负载噪声计算螺旋桨噪声。

Description

一种非均匀入流中螺旋桨噪声数值预报方法

技术领域

本发明涉及一种基于非均匀入流中螺旋桨噪声数值预报方法，属于螺旋桨辐射噪声数值预报技术领域。

背景技术

螺旋桨噪声包括空化噪声和非空化噪声，非空化状态时即螺旋负载噪声，空化状态时包括空化噪声和负载噪声两部分，其中螺旋桨空化的黏性数值模拟受多相流模型、湍流模型、空化模型和相变临界压力的影响，国内外对片空化进行成功模拟的很少，也没有较优的空化模型和湍流模型，成功模拟片空化的都是针对均匀入流；螺旋桨空化辐射噪声的研究仅限于采用面元法和单个脉冲球形空泡辐射噪声理论想结合，或者采用FW-H方程，将空化的脉动体积转化为噪声，这些方法都没有很好地跟实际数据对比验证，因此非均匀入流中螺旋桨空化单极子辐射噪声是一个研究难题，目前成功模拟螺旋桨噪声仅限于探讨性研究，没有与实际数据对比验证，因此准确数值预报螺旋桨辐射噪声也是一个研究难题。

发明内容

发明目的：针对现有技术中存在的问题与不足，本发明提供一种非均匀入流中螺旋桨噪声数值预报方法，该方法通过CFD数值模拟仿真非均匀入流中螺旋桨片空化的周期形态和空化面积变化，通过将螺旋桨的片空化辐射噪声等效为一个大的球形空化辐射，将螺旋桨的片空化面积等效为空化长度，通过空化长度依据单空泡空化辐射噪声原理求解螺旋桨片空化单极子辐射噪声，并通过声学数值仿真螺旋桨负载噪声，最后得到螺旋桨噪声。

技术方案：一种非均匀入流中螺旋桨噪声数值预报方法，包括以下步骤：第一步：划分网格螺旋桨计算域网格、检查网格划分质量并定义边界条件；第二步：用一阶差分插值计算入流速度；第三步：在CFD软件中导入计算网格并进行相关网格操作；第四步：在CFD软件中选择计算模型（定义求解器、选择湍流方程及设置对应参数、选择多相流模型及设置对应参数、设置操作环境、设定边界条件）；第五步：在CFD软件中设置求解参数、进行初始化设置以及收敛设置；第六步：在CFD软件中通过设置不同的操作环境进行稳态计算计算螺旋桨的淌水水动力参数值，并与实际值进行对比验证，通过后处理对比入流口速度与实际速度值，验证模型设置的准确性；第七步：将稳态计算结果作为非稳态计算的初始值，在CFD软件中进行非稳态迭代计算；第八步：在CFD软件中进行后处理，显示螺旋桨片空化的周期形态以及记录螺旋桨空化面积变化；第九步：将CFD计算的每个时间步的流场数据导入到声学数值仿真软件中，进行流场数据转换；第十步：划分声学网格和场点网格，进行声学边界元法的相关设置；第十一步：进行声学边界元数值计算，导出场点声压频率响应计算螺旋桨负载噪声；第十二步：计算螺旋桨噪声，并验证准确性。该数值预报方法通过数值预报螺旋桨片空化的面积，根据球形空泡辐射噪声理论，将螺旋桨空化辐射噪声看成是单个大空泡辐射噪声，从而计算螺旋桨片空化单极子辐射噪声，并通过CFD流场数据数值计算螺旋桨负载噪声，最后通过螺旋桨空化单极子辐射噪声和负载噪声计算螺旋桨噪声。

本发明的原理是利用单空泡空化辐射噪声原理求取螺旋桨片空化单极子辐射噪声，通过声学数值仿真螺旋桨负载噪声，最后得到螺旋桨噪声。

本发明方法的具体实现方式如下：

第一步，对螺旋桨计算域进行网格划分：

设定坐标原点在螺旋桨中心点处，螺旋桨的旋转轴为x轴，螺旋桨的直径为D，假定x正方向为下游，x负方向为上游；建立计算动域和计算静止域。进行网格划分，检查网格质量；定义流体计算域和边界条件；

第二步，采用一阶差分插值计算非均匀入流速度面的速度值：

实测入非均匀入流速度面的三向速度值（轴向、径向和切向三个方向），测试点的位置分布在r/R＝0.2,0.4,0.6,0.8,1.0…0.2M_test，其中R为螺旋桨半径，r为测试点的半径，M_test为测试半径的个数；

入流速度面的面网格划分采用结构网格划分，径向的节点数为M，一般情况下，M_test＜M，在径向方向进行一阶插值计算，径向插值的点数为M_r；

第三步，导入计算网格并进行网格相关操作：

将网格文件读入到CFD软件中（Fluent）进行网格检查，确保最小网格体积大于零，否则重新划分网格；进行计算域尺寸的调整，调整网格的尺寸比例，使最后尺寸符合实际模型的大小；将网格进行交换和光滑处理，依次按交换系数从小到大进行光滑和交换操作，直到每个交换系数的交换网格数目为零；读入非均匀入流速度文本文件；

第四步，计算模型设置：

在CFD软件中进行计算模型设置，定义求解器，按照默认的设置；选择K-epsilon模型，在K-epsilon模型下保存默认选项Stadard，在Near-Wall Treatment选项下选择标准壁面函数；流体介质的选取，选取water-liquid和water-vapor，根据实际值设置流体介质的密度和粘性系数；设置操作环境，在Operating Pressure中根据实际值写入环境压力；将螺旋桨的旋转速度的单位选为rpm；设置边界条件，包括选择流体动域的坐标系，稳态计算中选取MRF，非稳态计算中选取Moving Mesh，定义动域旋转轴为x轴，根据实际值设定选准速度，速度方向遵循右手定则，默认流体静止域为静止坐标系，非均匀入流速度面中速度采用笛卡尔坐标系，分别选择profile读入文本文件中的v_x，v_y，v_z，静止域外围入流速度面中速度采用笛卡尔坐标系，x方向的速度值根据实际的入流速度进行设定，压力出口选择静压力为零，定义螺旋桨的桨叶和桨毂壁面速度无滑移，粗糙系数为0。定义交界面；

第五步，设置求解参数、初始化以及收敛条件：

设置松弛因子，其中vapor，Turbulent Kinetic Energy，Specifiec DissipationRate，Turbulent Viscosity的松弛因子为λ，其它保持不变，定义差分方程形式，其中Pressure为Standard，其它的均为First Order Upwind形式；将入流口速度面1的速度作为初始值；设置各参数的残差；

第六步，在CFD软件中通过设置不同的操作环境进行稳态计算计算螺旋桨的淌水水动力参数值(推力系数和力矩系数)，并与实际值进行对比验证，通过后处理对比入流口速度与实际速度值，验证模型设置的准确性：

推力系数和力矩系数的计算公式：

$\left\{\begin{array}{c}{k}_T=\frac{T}{\mathrm{\ρ}{n}^2{D}^4}\\ k_Q=\frac{Q}{\mathrm{\ρ}{n}^2{D}^5}\end{array}\right.$ 式（1）

其中式（1）中 $T=\sqrt{T_x^2+T_y^2+T_z^2},Q=\sqrt{Q_x^2+Q_y^2+Q_z^2},$ T_x，T_y，T_z，Q_x，Q_y，Q_z从后处理Report Forces中读取，Q_x，Q_y，Q_z的力矩参考点为原点，ρ＝1000kg/m³，n为螺旋桨转速，D为螺旋桨直径；

通过Display Contours显示入流口轴向、切向和径向速度的等值面图，与实际速度进行对比，验证模型的准确性；

第七步，进行非稳态迭代计算：

将计算模型更改为非稳态计算，选择多相流模型为混合二相模型，设置空化类型为Singhal空化模型，并设置相应的参数，如饱和压力，单位体积的气泡数目等，设置混合相的第一相为water-liquid，第二相位water-vapor；将第六步中的稳态迭代结果作为初始条件，进行非稳态迭代计算，迭代时间步长为Δt，每个迭代步的迭代次数为m（根据具体例子选取不同的迭代次数）；

第八步，在CFD软件中进行后处理，显示螺旋桨片空化的周期形态以及记录螺旋桨空化面积变化：

在动域螺旋桨叶片上设置汽化体积分数为a的等值线图，通过Display Coutour显示螺旋桨吸力面上汽化体积分数α_v＝a的等值线图，对比相同条件下的实际片空化图，验证空化模型的准确性；记录螺旋桨每个迭代步的空化面积其中N为总迭代步数；

第九步，将CFD计算的每个时间步的流场数据导入到声学数值仿真软件中，并将时域流场数据进行数据转换：

在CFD软件中导出声源数据文件，其中声源是螺旋桨壁面包括螺旋桨桨叶面和桨毂面，这里声学仿真软件采用的Virtual Lab.，则声源数据文件为CGNS格式，将CGNS格式文件导入到声学仿真软件中，将CFD计算得到的文件由单元中心转换到单元节点上，保存为CFDData文件；

第十步，划分声学网格和场点网格，进行声学边界元法的相关设置：

将螺旋桨桨毂面进行封闭，按照每个波长至少6个单元的原则划分声学网格，建立以螺旋桨中心为中心，半径为5D的球形场点网格，以及生成p₁（0，0.5m，0），p₂（0，1.0m，0），p₃（0，1.5m，0）和p₄（0，2.0m，0）四个场点；

导入CFDData文件，将CFD节点流场数据进行快速傅里叶变换转移到声学网格上；设置声学边界条件，选取偶极子声压作为速度边界条件，采用不可压流动计算；

第十一步，进行声学响应计算声压频率响应函数计算，导出场点p₁，p₂，p₃，p₄的声压响应曲线，设置参考声压为1upa计算声压级频率响应，即螺旋桨负载噪声的声压级计算；

第十二步，计算螺旋桨噪声：

如果螺旋桨不产生空化，则在CFD仿真计算中省去第七步和第八步，则第十一步的计算结果即为螺旋桨噪声；如果螺旋桨产生空化，采用下面的公式进行计算，

根据式（2）计算螺旋桨空化面积的一阶差分，

${D1}_i=\frac{S_{c_{i+1}}-S_{c_i}}{\mathrm{\Δt}}$ 式（2）

其中i＝1,2,…N-1，根据式（3）计算空化面积的二阶差分，

${D2}_j=\frac{{D1}_{j+1}-{D1}_j}{\mathrm{\Δt}}$ 式（3）

其中j＝1,2,…N-2，根据式（4）计算等效空化长度，

${\mathrm{DV}}_k=6S_{c_k}{D1}_k^2+3S_{c_k}^2{D2}_k$ 式（4）

其中k＝1,2,…N-2，然后根据公式（5）计算螺旋桨空化单极子辐射噪声声压，

$p_k=\frac{{\mathrm{\ρ}}_{\mathrm{water}\_\mathrm{liquid}}}{4\mathrm{\πr}}{\mathrm{DV}}_k,k=\mathrm{1,2},...N-2$ 式（5）

根据公式（6）进行离散傅里叶变换，

$P\left(l\right)={\left|\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}{p}_k\exp \right\{-j2\mathrm{\πnl}/N\left\}\right|}^2$ 式（6）

在式（6）中，l＝0,1,…N/2+1，||表示绝对值运算，j表示虚数单位，即 $j=\sqrt{-1};$

第十一步中计算得到的负载噪声声压级为p_loading(l)，其中l＝0,1,…N/2+1，则螺旋桨噪声按照式（7）计算，

p(l)＝P(l)+p_loading(l)式 （7）

其中l＝0,1,…N/2+1。

有益效果：与现有技术相比，本发明具有以下有益效果：

1.改进空化模型，选择相对应的湍流模型，成功模拟非均匀入流条件下螺旋桨片空化周期形态，对比国际数据的淌水性能参数以及周期空化图片，有效验证空化模型的准确性，为下一步预报非均匀入流下螺旋桨空化单极子声源辐射噪声提供保证。

2.通过预报的周期片空化面积，将周期片空化面积转换为空化特征长度，通过空化特征长度以及脉动球体体积声源辐射噪声原理求解螺旋桨空化单极子声源辐射噪声，通过对比实际数据的叶频信息与叶倍频信息，验证预报的准确性。

3.进行CFD仿真，通过声学数值计算将流体信息转换为螺旋桨负载噪声，并通过国际数据对比，验证螺旋桨负载噪声预报的准确性，为进一步预报螺旋桨噪声提供保证。

4.通过螺旋桨的负载噪声和螺旋桨片空化单极子声源辐射噪声的频域相加得到空化状态下的螺旋桨噪声，并与实际的数据对比验证数值预报方法的正确性，开辟了新的数值预报螺旋桨辐射噪声的方法。

附图说明

图1为本发明实施例的方法流程图；

图2为网格计算模型图；

图3为三周期非均匀周向入流速度分布图（在旋转区域圆内）；

图4为测试点（0，2.0m，0）处，DTMB4119空化非空化下负载噪声声压级；

图5为测试点（0，2.0m，0）处，DTMB4119空化和非空化下噪声声压级；

图6为DTMB4119螺旋桨空化和非空化状态下半径为5D的球形面的声压级分布图（参考声压为1uPa），其中图（a）为DTMB4119无空化，图（b）为DTMB4119空化；

图7为民船桨周向片空化模拟，黑色图片为实验图，白色图片为汽相体积分数α_v＝0.1的等值面，其中图(a)J＝0.6546,U₀＝4.062m/s,σ＝2.8054，图(b)J＝0.6,U₀＝3.752m/s,σ＝2.7；

图8为两种不同转速民船桨的辐射噪声解调谱图，其中图(a)n＝25rps，图(b)n＝28rps；

图9为两种不同转速条件下民船桨片空化单极子声源辐射声压分布图；

图10为两种不同转速条件下民船桨片空化单极子辐射噪声声压频谱图；

图11为民船桨测试点（0，2.0m，0）处仿真数据的声压级，其中图(a)n＝15rps中，图(b)n＝20rps。

具体实施方式

下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

如图1所示，非均匀入流中螺旋桨噪声数值预报方法，包括以下步骤：

第一步，对螺旋桨计算域进行网格划分：

设定坐标原点在螺旋桨中心点处，螺旋桨的旋转轴为x轴，螺旋桨的直径为D，假定x正方向为下游，x负方向为上游；建立中心在原点的圆柱体1，圆柱体1的轴向方向为x轴，直径为1.06D，长为0.8D，上游距离0.32D，下游距离0.48D；建立中心在原点的圆柱体2，圆柱体2的轴向方向为x轴，直径为1.471D，长为4.25D，上游距离0.5D，下游距离3.75D；圆柱体1减去螺旋桨实体区域为计算动域，圆柱体2减去圆柱体1的圆环柱区域为计算静域；计算静止域再细划分为五个区域，如图2所示①②③④⑤，其中将区域⑤等分为四个扇形区域从而划分结构网格；动域采用非结构网格划分，静止区域采用结构网格划分，其中区域①⑤的网格划分较区域②③④的网格密，划分网格遵循自下而上的原则，首先每个区域的线网格，然后是面网格，最后是体网格；划分好网格后，检查网格质量，对于螺旋桨计算区域，一般网格数量大概在两百万到三百万之间；定义区域①②③④⑤为静止域，包围螺旋桨的圆柱体区域为动域；如图2所示定义速度入口面1和2，压力出口面，交界面即为各个区域的交界处（设置交界面的条件是两个区域的交界面不共面或者划分的网格不一致），壁面为螺旋桨桨叶压力面、吸力面和桨毂面（以上操作在软件Gambit中进行）；

第二步，采用一阶差分插值计算入流速度面1的速度值：

实测入流速度面1的三向速度值（轴向、径向和切向三个方向），测试点的位置分布在r/R＝0.2,0.4,0.6,0.8,1.0…0.2M_test，其中R为螺旋桨半径，r为测试点的半径，M_test为测试半径的个数；

入流速度面1的面网格划分采用结构网格划分，径向的节点数为M，一般情况下，M_test＜M，在径向方向进行一阶插值计算，径向插值的点数为M_r；

插值计算公式为

$v_r=\left\{\begin{array}{c}\\ v_{0.2R}+\frac{v_{0.4R}-v_{0.2R}}{0.2R}\&CenterDot;(r-0.2R),\\ r>0.2R\\ v_{0.2m_{1}R}+\frac{{v_{0.2(m_1+1)R}-v}_{0.2m_{1}R}}{0.2R}\&CenterDot;(r-0.2m_{1}R),\\ 0.2m_{1}R<r<0.2(m_1+1)R,1\&le;m_1<M_{\mathrm{test}},0.2R<r<0.2M_{\mathrm{test}}R\\ v_{0.2M_{\mathrm{test}}R}+\frac{v_{0.2M_{\mathrm{test}}R}-v_{0.2(M_{\mathrm{test}}-1)R}}{0.2R}\&CenterDot;(r-0.2M_{\mathrm{test}}R),\\ r>0.2M_{\mathrm{test}}R\end{array}\right.$ 式（1）

其中v_r为半径r的插值速度，θ_i＝iθ₀,i＝0,1,2,…N,N＝π/θ₀+1，θ₀为测试点之间的间隔角度，θ_i为从上游看去顺时针方向与y轴的夹角，包括周向、切向和径向三个方向的速度，分别为r≤1.417D，v_r径向插值个数为M_r，M_r+M_test≤M；

对v_r进行坐标转换，坐标转换公式为

$\left\{\begin{array}{c}{v}_{x_i}=v_{r,{\mathrm{\&theta;}}_i}^a,x_i=D_1\\ v_{y_i}=+v_{r,{\mathrm{\&theta;}}_i}^t\cos \left({\mathrm{\&theta;}}_i\right)-v_{r,{\mathrm{\&theta;}}_i}^r\sin \left({\mathrm{\&theta;}}_i\right),y_i=r\cos \left({\mathrm{\&theta;}}_i\right)\\ v_{z_i}=v_{r,{\mathrm{\&theta;}}_i}^t\sin \left({\mathrm{\&theta;}}_i\right)+v_{r,{\mathrm{\&theta;}}_i}^r\cos \left({\mathrm{\&theta;}}_i\right),z_i=r\sin \left({\mathrm{\&theta;}}_i\right)\end{array}\right.$ 式（2）

测试点的个数为N(M_test+M_r)，进行上述坐标转换后拼接为矩阵

$\left(\begin{array}{cccccc}{x}_1& y_1& z_1& v_{x_1}& v_{y_1}& v_{z_1}\\ x_2& y_2& z_2& v_{x_2}& v_{y_2}& v_{z_2}\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ x_{N(M_{\mathrm{test}}+M_r)}& y_{N(M_{\mathrm{test}}+M_r)}& z_{N(M_{\mathrm{test}}+M_r)}& v_{x_{N(M_{\mathrm{test}}+M_r)}}& v_{y_{N(M_{\mathrm{test}}+M_r)}}& v_{z_{N(M_{\mathrm{test}}+M_r)}}\end{array}\right)$

将上述矩阵利用边界分布剖面函数 $\begin{array}{c}(\left(\begin{array}{ccc}\mathrm{velocity}& \mathrm{point}& N(M_{\mathrm{test}}+M_r)\end{array}\right)\\ \left(\begin{array}{ccc}x& x_1& x_2......x_{N(M_{\mathrm{test}}+M_r)}\end{array}\right)\\ \left(\begin{array}{ccc}y& y_1& y_2......y_{N(M_{\mathrm{test}}+M_r)}\end{array}\right)\\ \left(\begin{array}{ccc}z& z_1& z_2......z_{N(M_{\mathrm{test}}+M_r)}\end{array}\right)\\ \left(\begin{array}{ccc}{v}_x& v_{x_1}& v_{x_2}......v_{x_{N(M_{\mathrm{test}}+M_r)}}\end{array}\right)\\ \left(\begin{array}{ccc}{v}_y& v_{y_1}& v_{y_2}......v_{y_{N(M_{\mathrm{test}}+M_r)}}\end{array}\right)\\ \left(\begin{array}{ccc}{v}_z& v_{z_1}& {v_{z_2}......v}_{z_{N(M_{\mathrm{test}}+M_r)}}\end{array}\right)\end{array}$ 写成速度文本文件；

第三步，导入计算网格并进行网格相关操作：

将网格文件读入到CFD软件中（Fluent）进行网格检查，确保最小网格体积大于零，否则重新划分网格；进行计算域尺寸的调整，调整网格的尺寸比例，使最后尺寸符合实际模型的大小；将网格进行交换和光滑处理，依次按交换系数从小到大进行光滑和交换操作，直到每个交换系数的交换网格数目为零；读入第二步中profile文本文件；

第四步，计算模型设置：

在CFD软件中进行计算模型设置，定义求解器，按照默认的设置；选择K-epsilon模型，在K-epsilon模型下保存默认选项Stadard，在Near-Wall Treatment选项下选择标准壁面函数；流体介质的选取，选取water-liquid和water-vapor，根据实际值设置流体介质的密度和粘性系数；设置操作环境，在Operating Pressure中根据实际值写入环境压力；将螺旋桨的旋转速度的单位选为rpm；设置边界条件，包括选择流体动域的坐标系（稳态计算中选取MRF，非稳态计算中选取Moving Mesh），定义动域旋转轴为x轴，根据实际值设定选准速度，速度方向遵循右手定则，默认流体静止域为静止坐标系，入流速度面1中速度采用笛卡尔坐标系，分别选择profile读入文本文件中的v_x，v_y，v_z，入流速度面2中速度采用笛卡尔坐标系，x方向的速度值根据实际的入流速度进行设定，压力出口选择静压力为零，定义螺旋桨的桨叶和桨毂壁面速度无滑移，粗糙系数为0。定义交界面；

第五步，设置求解参数、初始化以及收敛条件：

设置松弛因子，其中vapor，Turbulent Kinetic Energy，Specifiec DissipationRate，Turbulent Viscosity的松弛因子为0.5，其它保持不变，定义差分方程形式，其中Pressure为Standard，其它的均为First Order Upwind形式；将入流口速度面1的速度作为初始值；设置各参数的残差；

第六步，在CFD软件中通过设置不同的操作环境进行稳态计算计算螺旋桨的淌水水动力参数值(推力系数和力矩系数)，并与实际值进行对比验证，通过后处理对比入流口速度与实际速度值，验证模型设置的准确性：

推力系数和力矩系数的计算公式：

$\left\{\begin{array}{c}{k}_T=\frac{T}{\mathrm{\&rho;}{n}^2{D}^4}\\ k_Q=\frac{Q}{\mathrm{\&rho;}{n}^2{D}^5}\end{array}\right.$ 式（3）

其中式（3）中 $T=\sqrt{T_x^2+T_y^2+T_z^2},Q=\sqrt{Q_x^2+Q_y^2+Q_z^2},$ T_x，T_y，T_z，Q_x，Q_y，Q_z从后处理Report Forces中读取，Q_x，Q_y，Q_z的力矩参考点为原点，ρ＝1000kg/m³，n为螺旋桨转速，D为螺旋桨直径；

通过Display Contours显示入流口轴向、切向和径向速度的等值面图，与实际速度进行对比，验证模型的准确性；

第七步，进行非稳态迭代计算：

将计算模型更改为非稳态计算，选择多相流模型为混合二相模型，设置空化类型为Singhal空化模型，并设置相应的参数，如饱和压力，单位体积的气泡数目等，设置混合相的第一相为water-liquid，第二相位water-vapor；将第六步中的稳态迭代结果作为初始条件，进行非稳态迭代计算，迭代时间步长为Δt，每个迭代步的迭代次数为m（根据具体例子选取不同的迭代次数）；

第八步，在CFD软件中进行后处理，显示螺旋桨片空化的周期形态以及记录螺旋桨空化面积变化：

在动域螺旋桨叶片上设置等值面α_v＝0.1，通过Display Coutour显示螺旋桨吸力面上汽化体积分数α_v＝0.1的等值线图，对比相同条件下的实际片空化图，验证空化模型的准确性；记录螺旋桨每个迭代步的空化面积其中N为总迭代步数；

螺旋桨的转速为nrps，在迭代步时间Δt内转到的角度为360nΔt，螺旋桨转一圈所需的迭代步数为1/(nΔt)，记录螺旋桨片空化面积，保证N≥2/(nΔt)，即至少转动两圈的数据；

第九步，将CFD计算的每个时间步的流场数据导入到声学数值仿真软件中，并将时域流场数据进行数据转换：

在CFD软件中导出声源数据文件，其中声源是螺旋桨壁面包括螺旋桨桨叶面和桨毂面，这里声学仿真软件采用的Virtual Lab.，则声源数据文件为CGNS格式，将CGNS格式文件导入到声学仿真软件中，将CFD计算得到的文件由单元中心转换到单元节点上，保存为CFDData文件；在声学数值仿真软件Virtual Lab.中选择Data Transfer Analysis Case，将Meshing Mapping选为Creat a New Mesh Mapping，Source选为Node and Eelement下的Export_Q_000001->centroid3d（CFD保存文件时的默认名称），Target选为Export_Q_000001->node and elements，最后在Data Transfer Solution Set.1中Update。

第十步，划分声学网格和场点网格，进行声学边界元法的相关设置：

将螺旋桨桨毂面进行封闭，按照每个波长至少6个单元的原则划分声学网格，建立以螺旋桨中心为中心，半径为5D的球形场点网格，以及生成p₁（0，0.5m，0），p₂（0，1.0m，0），p₃（0，1.5m，0）和p₄（0，2.0m，0）四个场点；

导入CFDData文件，将CFD节点流场数据进行快速傅里叶变换转移到声学网格上；通过Check Id Conflicts进行网格节点及单元编号冲突检查，定义流体材料为Water，设定材料属性（密度和速度），定义流体属性，选取偶极子声压作为速度边界条件，采用不可压流动计算；

第十一步，通过Acoustic Response Analysis Case进行声学计算，添加偶极子声源作为边界条件；通过Vector to Function Conversion Case计算声压频率响应函数，设置参考声压为1upa计算声压级频率响应，即螺旋桨负载噪声的声压级计算；

第十二步，计算螺旋桨噪声：

如果螺旋桨不产生空化，则在CFD仿真计算中省去第七步和第八步，则第十一步的计算结果即为螺旋桨噪声；如果螺旋桨产生空化，采用下面的公式进行计算，

根据式（2）计算螺旋桨空化面积的一阶差分，

${D1}_i=\frac{S_{c_{i+1}}-S_{c_i}}{\mathrm{\&Delta;t}}$ 式（2）

其中i＝1,2,…N-1，根据式（3）计算空化面积的二阶差分，

${D2}_j=\frac{{D1}_{j+1}-{D1}_j}{\mathrm{\&Delta;t}}$ 式（3）

其中j＝1,2,…N-2，根据式（4）计算等效空化长度，

${\mathrm{DV}}_k=6S_{c_k}{D1}_k^2+3S_{c_k}^2{D2}_k$ 式（4）

其中k＝1,2,…N-2，然后根据公式（5）计算螺旋桨空化单极子辐射噪声声压，

$p_k=\frac{{\mathrm{\&rho;}}_{\mathrm{water}\_\mathrm{liquid}}}{4\mathrm{\&pi;r}}{\mathrm{DV}}_k,k=\mathrm{1,2},...N-2$ 式（5）

根据公式（6）进行离散傅里叶变换，

$P\left(l\right)={\left|\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{p}_k\exp \right\{-j2\mathrm{\&pi;nl}/N\left\}\right|}^2$ 式（6）

在式（6）中，l＝0,1,…N/2+1，||表示绝对值运算，j表示虚数单位，即 $j=\sqrt{-1};$

第十一步中计算得到的负载噪声声压级为p_loading(l)，其中l＝0,1,…N/2+1，则螺旋桨噪声按照式（7）计算，

p(l)＝P(l)+p_loading(l) 式（7）

其中l＝0,1,…N/2+1。

具体来说：

实施例1

验证的桨模为DTMB4119桨，无侧斜无纵斜分布的三叶螺旋桨，螺旋桨直径为0.305m，盘面比0.6，叶切面为NACA-66，毂径比为0.2。CFD计算区域模型如下：上游的距离为0.7D，下游的距离为3.75D，旋转区域的半径为0.16165m，长为0.14m，静止区域的半径为1.471D，采用Gambit划分网格，旋转区域采用非结构网格划分，网格数目为572382，静止区域采用结构网格划分，网格数目为1762748。CFD计算的参数设置如下：静止压力为101325Pa，进口入流速度为5.808m/s，螺旋桨转速为20rps，采用滑移网格技术，湍流模型采用k-εRNG模型，计算得到该桨敞水性能与实验值比较如表1所示。

表1DTMB4119桨敞水性能测试验证

J	计算推力系数kT	实验值kT	计算力矩系数10kQ	实验值10kQ
					0.5	0.285	0.2826	0.477	0.4851
0.7	0.2	0.2059	0.36	0.3787
					0.833	0.146	0.1471	0.28	0.2943
0.9	0.12	0.1152	0.239	0.2474

非均匀入流采用三周期非均匀轴向入流，轴向入流速度的计算公式如下：

$v_{\mathrm{axis}}=v_{\mathrm{ship}}[A_0+\underset{n-1}{\overset{4}{\mathrm{\&Sigma;}}}(A_n\cos n{\mathrm{\&theta;}}_0+B_n\sin n{\mathrm{\&theta;}}_0)]$

其中螺旋桨的轴向为x方向，θ₀是从下游方向看去，与y方向的夹角，v_n为船速，即在试验中为均匀入流速度。三周期非均匀轴向入流的参数是现有的测试数据，例如MIT PhDdissertation1979年《Prediction of steady and unsteady performance of marinepropeller with or without cavitation by numerical lifting surface theory》中公开的内容。

CFD计算的参数设置如下：静止压力为101325Pa，入流速度中间旋转圆盘的采用三周期非均匀轴向入流速度，从0到360度，间隔为5度，其它入流面的轴向速度为5.808m/s，在空化状态下设置空化指数为1.7，相应的静止压力做调整。图3为三周期非均匀入流轴向速度分布图。

图4为螺旋桨在空化和非空化状态的负载噪声的声压级曲线比较，从图中可以看出，空化状态下的声压级较非空化状态下较小。图6螺旋桨空化和非空化状态下半径为5D的球形面的声压级分布，从图中可以得知，非空化状态下有偶极子对称性，而在空化状态下，偶极子对称性不明显。图5为螺旋桨空化和非空化状态下螺旋桨噪声声压级曲线。

实施例2

某民船桨的几何参数如下：叶片数为4，螺旋桨直径为0.2482m，盘面比为0.55，侧斜角为32度。

图7为民船桨周向片空化模拟对比验证图，在轴向均匀入流速度U₀为3.25m/s，环境压力113000Pa，螺旋桨转速n分别为25rps和28rps，每一步迭代时间为0.0005秒，两种不同转速对应的进速系数J分别为0.5238、0.4677，其空化指数σ分别为5.747、4.58，计算该桨片空化单极子声源辐射声压，其中r＝1。图9为螺旋桨转动一周0.04秒内的声压分布图。两种不同转速的片空化辐射噪声声压频谱图，如图10。每个迭代步的时间为0.0005秒，则频率分辨率为2000Hz，依据采样定理，有效分析频率为1000Hz。

从图10可以看出，采用文中的方法可以模拟螺旋桨空化单极子声源辐射噪声，功率谱密度图清晰显示了螺旋桨的叶频信息，两种不同转速的叶频频率值分别是100和112，与螺旋桨转速吻合，叶频幅值分别是0.04367，0.955，叶频幅值随进速系数和空化指数的减小而增大。对比图8，可以看出，预报的空化噪声的缓变分量部分与实际水听器测量的声压功率谱很好的吻合，并准确预报了螺旋桨的叶频信息。

图11是民船桨测试点（0，2.0m，0）处测得的声压级频率响应图，仿真的最大频率分辨率为1000Hz，根据总声压级计算公式：

$\mathrm{SPL}=20{\log }_{10}[\left(\underset{i}{\mathrm{\&Sigma;}}p\left(i\right)\right)/p_0]$

其中p(i)为对应频率值的频率响应值，p₀为参考声压，值为1upa。根据公式计算1000Hz以内的总声压级，当转速为15rps时，在测试点（0，0.5m，0）、（0，1.0m，0）、（0，1.5m，0）和（0，2.0m，0）分别计算的结果为200.2233、199.5429、199.074和198.6839，而实际的总声压级为198.1594，有较好的吻合；当转速为20rps时，在上述测试点分别计算的结果为199.1840、198.4708、198.0929和198.0929，而实际的总声压级为195.2954，也有较好的吻合。总声压级的计算精度与分辨率相关，分辨率越小越精确，在文中的仿真分析中，转速为15rps时的分辨率为7.5Hz，转速为20rps时的分辨率为10Hz，因此在计算实际数据的总声压级时也保持实际数据和仿真数据的分辨率一致。仿真很好体现偶极子声源分布特征，原因有两点，要出现很明显的偶极子分布特征必须保证非均匀入流时单峰值周向非均匀入流，此外偶极子明显出现在叶频和叶倍频处。表1为民船桨在叶频、二倍叶频以及1000Hz以内仿真数据和实际数据的声压级比较，数据分析显示，仿真的噪声声压级在叶频和倍叶频以及低频内的声压级有较好的一致性。

表1民船桨实际数据和仿真数据的声压级比较表

Claims (1)

1.一种非均匀入流中螺旋桨噪声数值预报方法，其特征在于，步骤包括：

首先，划分CFD计算网格，进行稳态迭代计算淌水性能参数和入流口速度验证模型的准确性；CFD是Computational Fluid Dynamics的简称，即计算流体动力学；

接下来，进行CFD非稳态计算，记录螺旋桨片空化周期形态以及片空化面积变化；

然后，通过CFD流场数据，进行声学边界元数值计算，计算螺旋桨负载噪声；

最后，通过螺旋桨片空化单极子辐射噪声以及负载噪声计算螺旋桨噪声；

划分CFD计算网格，进行稳态迭代计算淌水性能参数和入流口速度验证模型的准确性具体包括以下步骤：

第一步，对螺旋桨计算域进行网格划分：

设定坐标原点在螺旋桨中心点处，螺旋桨的旋转轴为x轴，螺旋桨的直径为D，假定x正方向为下游，x负方向为上游；建立计算动域和计算静止域；进行网格划分，检查网格质量；定义流体计算域和边界条件；

第二步，采用一阶差分插值计算非均匀入流速度面的速度值：

实测非均匀入流速度面的三向速度值，测试点的位置分布在r/R＝0.2,0.4,0.6,0.8,1.0…0.2M_test，其中R为螺旋桨半径，r为测试点的半径，M_test为测试半径的个数；

入流速度面1的面网格划分采用结构网格划分，径向的节点数为M，M_test＜M，在径向方向进行一阶插值计算，径向插值的点数为M_r；

第三步，导入计算网格并进行网格相关操作：

将网格文件读入到CFD软件中进行网格检查，确保最小网格体积大于零，否则重新划分网格；进行计算域尺寸的调整，调整网格的尺寸比例，使最后尺寸符合实际模型的大小；将网格进行交换和光滑处理，依次从交换系数从小都大进行光滑和交换操作，知道每个交换系数的交换网格数目为零；读入非均匀入流速度面的文本文件；

第四步，计算模型设置：

在CFD软件中进行计算模型设置，定义求解器，按照默认的设置；选择K-epsilon模型，在K-epsilon模型下保存默认选项Stadard，在Near-Wall Treatment选项下选择标准壁面函数；流体介质的选取，选取water-liquid和water-vapor，根据实际值设置流体介质的密度和粘性系数；设置操作环境，在Operating Pressure中根据实际值写入环境压力；将螺旋桨的旋转速度的单位选为rpm；设置边界条件，包括选择流体动域的坐标系，稳态计算中选取MRF，非稳态计算中选取Moving Mesh，定义动域旋转轴为x轴，根据实际值设定选准速度，速度方向遵循右手定则，默认流体静止域为静止坐标系，非均匀入流速度面中速度采用笛卡尔坐标系，分别选择profile读入文本文件中的v_x，v_y，v_z，静止域外围入流速度面中速度采用笛卡尔坐标系，x方向的速度值根据实际的入流速度进行设定，压力出口选择静压力为零，定义螺旋桨的桨叶和桨毂壁面速度无滑移，粗糙系数为0；定义交界面；

第五步，设置求解参数、初始化以及收敛条件：

设置松弛因子，其中vapor，Turbulent Kinetic Energy，Specifiec DissipationRate，Turbulent Viscosity的松弛因子为λ，其它保持不变，定义差分方程形式，其中Pressure为Standard，其它的均为First Order Upwind形式；将入流口速度面1的速度作为初始值；设置各参数的残差；

第六步，在CFD软件中通过设置不同的操作环境进行稳态计算计算螺旋桨的淌水水动力参数值，即推力系数和力矩系数，并与实际值进行对比验证，通过后处理对比入流口速度与实际速度值，验证模型设置的准确性：

推力系数和力矩系数的计算公式：

其中式(1)中 T_x，T_y，T_z，Q_x，Q_y，Q_z从后处理ReportForces中读取，Q_x，Q_y，Q_z的力矩参考点为原点，ρ＝1000kg/m³，n为螺旋桨转速，D为螺旋桨直径；

通过Display Contours显示入流口轴向、切向和径向速度的等值面图，与实际速度进行对比，验证模型的准确性；

进行CFD非稳态计算，记录螺旋桨片空化周期形态以及片空化面积变化具体包括以下步骤：

第七步，进行非稳态迭代计算：

将计算模型更改为非稳态计算，选择多相流模型为混合二相模型，设置空化类型为Singhal空化模型，并设置相应的参数，参数包括饱和压力、单位体积的气泡数目，设置混合相的第一相为water-liquid，第二相位water-vapor；将第六步中的稳态迭代结果作为初始条件，进行非稳态迭代计算，迭代时间步长为Δt，每个迭代步的迭代次数为m；

第八步，在CFD软件中进行后处理，显示螺旋桨片空化的周期形态以及记录螺旋桨空化面积变化：

在动域螺旋桨叶片上设置汽化体积分数为a的等值线图，通过Display Coutour显示螺旋桨吸力面上汽化体积分数α_v＝a的等值线图，对比相同条件下的实际片空化图，验证空化模型的准确性；记录螺旋桨每个迭代步的空化面积其中N为总迭代步数；

通过CFD流场数据，进行声学边界元数值计算，计算螺旋桨负载噪声具体包括以下步骤：

第九步，将CFD计算的每个时间步的流场数据导入到声学数值仿真软件中，并将时域流场数据进行数据转换：

在CFD软件中导出声源数据文件，其中声源是螺旋桨壁面包括螺旋桨桨叶面和桨毂面，这里声学仿真软件采用的Virtual Lab.，则声源数据文件为CGNS格式，将CGNS格式文件导入到声学仿真软件中，将CFD计算得到的文件由单元中心转换到单元节点上，保存为CFDData文件；

第十步，划分声学网格和场点网格，进行声学边界元法的相关设置：

将螺旋桨桨毂面进行封闭，按照每个波长至少6个单元的原则划分声学网格，建立以螺旋桨中心为中心，半径为5D的球形场点网格，以及生成p₁(0，0.5m，0)，p₂(0，1.0m，0)，p₃(0，1.5m，0)和p₄(0，2.0m，0)四个场点；

导入CFDData文件，将CFD节点流场数据进行快速傅里叶变换转移到声学网格上；设置声学边界条件，选取偶极子声压作为速度边界条件，采用不可压流动计算；

通过螺旋桨片空化单极子辐射噪声以及负载噪声计算螺旋桨噪声具体包括以下步骤：

第十一步，进行声学响应计算声压频率响应函数计算，导出场点p₁，p₂，p₃，p₄的声压响应曲线，设置参考声压为1upa计算声压级频率响应，即螺旋桨负载噪声的声压级计算；

第十二步，计算螺旋桨噪声：如果螺旋桨不产生空化，则在CFD仿真计算中省去第七步和第八步，则第十一步的计算结果即为螺旋桨噪声；如果螺旋桨产生空化，采用下面的公式进行计算，

根据式(2)计算螺旋桨空化面积的一阶差分，

其中i＝1,2,…N-1，根据式(3)计算空化面积的二阶差分，

其中j＝1,2,…N-2，根据式(4)计算等效空化长度，

其中k＝1,2,…N-2，然后根据公式(5)计算螺旋桨空化单极子辐射噪声声压，

根据公式(6)进行离散傅里叶变换，

在式(6)中，l＝0,1,…N/2+1，| |表示绝对值运算，j表示虚数单位，即

第十一步中计算得到的负载噪声声压级为p_loading(l)，其中l＝0,1,…N/2+1，则螺旋桨噪声按照式(7)计算，

p(l)＝P(l)+p_loading(l) 式(7)

其中l＝0,1,…N/2+1。