螺旋桨旋涡发放数值模拟

张营川 马 骏

大连理工大学船舶工程学院，辽宁大连 116024

1 引 言

螺旋桨振动噪声，尤其是共振引起的异常噪声——“唱音”，是舰船三大噪声源之一。以水下隐蔽性为主要特征的潜艇和声自导鱼雷，为避免过早被敌方发现，又能及早发现敌方，必须具有很低的自噪声，这是提高潜艇、声自导鱼雷隐蔽性和战术性能的重要措施。因此，研究螺旋桨“唱音”具有重要意义。

螺旋桨“唱音”是由湍流激励引起螺旋桨桨叶结构共振而引发的非线性声学现象，这一现象的本质可以归结到螺旋桨的旋涡发放问题中，也就是说，在螺旋桨的旋转过程中，桨叶的随边产生了旋涡，当旋涡发放的频率与桨叶的固有频率接近时，螺旋桨桨叶由于共振便开始发出清脆的鸣音。虽然对于螺旋桨噪声的研究已开展很多［1-2］，但至目前为止，对于螺旋桨“唱音”现象的研究还只停留在依靠模型实验与经验公式阶段［3-5］，而且只能通过分析测量得到的螺旋桨的声学数据来研究避免“唱音”现象的发生。

近年来，螺旋桨推进领域的数值研究有所发展，出现了较多文献［6-9］，但对于螺旋桨“唱音”问题的数值研究还少有涉及。如果能够从计算机辅助工程的角度出发，建立大小合理的计算模型并能得到工程上精确的螺旋桨桨叶在不同来流速度和不同转速下的旋涡发放频率，就可以在螺旋桨的设计阶段通过输入螺旋桨的模型参数及转速和航速来计算出螺旋桨是否会发生“唱音”现象，这是非常具有实际应用意义的，其结果不仅可以减少模型实验的时间和费用，而且还可避免在船舶下水之后因发生“唱音”现象而再来修改桨叶随边形状等一系列耗费人力物力的工作。但在当前的计算能力下，不做任何简化，详细地建立螺旋桨与周围流体模型进行流固耦合计算还存在困难。

因此，本文从数值模拟的角度在螺旋桨 “唱音”问题上做了一些初步研究，即基于LS-DYNA软件以及任意拉格朗日—欧拉（ALE）方法对螺旋桨和周围流体建立了有限元模型，通过适当的简化，观察到了桨叶后方旋涡的运动特性，并给出了对应的解释，最后计算并分析了桨叶的结构响应及相关影响因素。这些结论对于后续开展螺旋桨涡激共振和“唱音”现象方面的数值研究具有一定的参考价值。

2 ALE方法的基本公式

ALE方法是描述流体运动的两种基本方法。在拉格朗日描述下，网格随流体一起运动，网格点的速度与当地流体微团的速度相同。在欧拉描述下，网格的空间位置固定，即网格点的速度为零，流体微团穿过网格单元所构成的控制体。文献［10］的方法是将这两种基本方法统一起来，允许网格以任意速度运动。当网格速度为零时，为欧拉方法，当网格速度等于流体运动速度时，则为拉格朗日方法。

在ALE方法的描述中，引入了拉格朗日和欧拉坐标之外的第3个任意参照坐标。与参照坐标相关的材料微商可采用下式描述：

式中，f为ALE方法中与坐标和时间相关的任意变量；Xi为拉格朗日坐标；xi为欧拉坐标；wi为相对速度。相对速度通过下式表达：

式中，v为物质速度；u为网格速度。

因此，用材料时间导数和参照几何构形的时间导数两者之间的替换关系可以推导出所需的ALE方程。

ALE算法的控制方程可由下列守恒方程给出：

1）质量守恒方程

式中，ρ为流体密度。

2）动量守恒方程

控制固定域上的牛顿流体流动问题的增强形式由控制方程和对应的初始及边界条件组成，控制流体问题的方程是Navier-Stokes方程的ALE描述：

式中，bi为单位质量的体积力分量；σij为Cauchy应力张量。

Cauchy应力张量能够以Stokes本构方程的形式表示，具体表达式如下：

式中，p为压力；δij为kronec ker δ函数；μ为流体的动力粘性系数。

式（5）可以与下列边界条件联立求解：

式（6）、式（7）中，nj为边界外法线单位向量；Ui0为指定初始边界。

并且，对于计算域有下式成立：

式（8）、式（9）中，Γ 为计算域的完整边界；Γ1，Γ2为计算域的部分边界。

假设整个计算域的速度场在t＝0时刻已知，即

式中，E为内能。

3 平面结构旋涡发放现象数值模拟

1）圆柱剖面

螺旋桨在旋转时其桨叶和流体的相互作用是非常复杂的，在模型实验中也难以观察到桨后流场情况，直接模拟这一流固耦合过程难免出现误差。因此，首先应针对最简单的平面圆柱在来流下的旋涡发放现象进行模拟，然后根据结果调整一系列的计算参数，只有这样，才能使后续的螺旋桨桨叶及流体的模拟尽量接近实际情况。

图1所示为平面圆柱和流体的有限元模型示意图以及局部网格放大图。由于采用了ALE方法处理圆柱与流体的流固耦合问题，因此在建模过程中，圆柱和流体可以分别建模和划分网格，相互之间互不干扰，网格可以重叠。在数值求解过程中，各时刻交替地求解流体和固体区域，并且在交替过程中通过罚函数方程在耦合界面进行有关物理量的传递，从而达到不同求解域的相互耦合。

圆柱的直径为B=0.001m，来流速度为v=4.5 m／s，Strouhal数取为 St=0.2，可以得出该圆柱的旋涡发放频率ns约为：

则圆柱的旋涡发放周期约为0.001 s。

图2为数值模拟得到的圆柱在垂直来流方向的受力曲线。流体流经圆柱时，圆柱每发放出一个旋涡，圆柱在垂直来流方向的受力便会改变一次大小和方向，因此，通过圆柱在垂直来流方向的受力曲线，可以计算出圆柱的旋涡发放频率约为0.001 s，这个结果与公式计算得到的结果基本吻合。图3给出了圆柱旋涡发放时，流体单元的应变云图以及整个流体域的速度云图。图4给出了圆柱旋涡发放时，整个流体域的流体速度云图（含速度方向）以及局部放大图像。

尽管如此，数值计算结果还是存在着不足，这主要是由于在建模过程中，考虑到实际的计算能力，为了减小整个有限元模型的网格数量，一方面因是采用单层三维网格来研究厚度很薄的圆柱模型与来流的相互作用，因此没有建立有一定长度的圆柱模型参与模拟；另一方面，本文对流体区域的网格大小进行了多次尝试，在网格较大时，得出的旋涡发放周期和圆柱受力都不准确，而在网格较小时，能得到相对较好的计算结果，但网格的细化使计算时间大幅增加。所以在计算条件允许的前提下，根据圆柱直径为0.001m的情况，流体网格大小应设定为0.00005m。

2）机翼剖面

在上述计算的基础上，本文数值模拟了机翼剖面的旋涡发放现象。选取的翼型为NACA 0012，翼展 0.6m，弦长 0.1 m。 计算模型右侧边界面施加了流入边界条件，来流速度为4.5 m／s，左侧边界面施加了流出边界条件，其他边界面施加的是无反射边界条件，整个水域施加了初始流入速度，可以保证流场的连贯；时间步长Δt受到Courant稳定条件的限制取得非常小，在计算过程中，还采用了LS-DYNA软件中的质量缩放方法来节省计算时间。在计算中，考虑了单层机翼剖面受到流体力作用后自身的应力和应变情况。

图5给出了机翼剖面开始旋涡发放以后的流体速度等值线图（左）和速度矢量图（右），并且在需要放大的位置进行了标记。图6为速度矢量图的放大图，在图中用圆圈标记了旋涡的位置以及旋涡旋转的方向，清晰地展示了单层机翼剖面在来流速度较低时的旋涡发放现象。

从单层机翼剖面的尾部最末端选取一个节点，图7中的曲线表示的就是该节点垂直来流方向的位移随时间的变化情况。该条曲线可以用来读取机翼剖面的旋涡发放频率等重要信息。除流体数据外，数值模拟还得到了机翼剖面在流体力作用下自身的应力分布情况，从图8中可以看出，颜色深的部分为机翼剖面的高应力区域，而颜色较深区域中的结构单元应力比较小。

需要说明的是，由于控制旋涡发放的是机翼尾部最后端的宽度 （与机翼的整体尺寸相比非常小），因此流体网格必须按照此宽度进行设定才能准确地与固体网格耦合，从而导致整个计算模型的单元数量非常多。另外，求解含有结构应力应变的流固耦合问题要比单独求解刚性固体的流固耦合问题复杂得多，需要的计算时间也更长。为了保证能在短时间内模拟出机翼剖面的旋涡发放现象，只能将流体网格适当放大，但由此会导致旋涡发放频率与实际有所差别。要得到比较精确的翼型剖面旋涡发放频率，就需要将本文的模型更加细化。

4 固定螺旋桨桨叶旋涡发放现象的数值模拟

本文选定四叶的MAU4-55螺旋桨作为研究对象。虽然应用本文方法不仅可以求出桨叶后方的流场特性，而且还可计算出桨叶内部应力等结构力学性能，但是数值模拟螺旋桨整桨的流固耦合现象需要相当长的计算时间，为了提高计算速度，本文假设：1）只研究螺旋桨的一个桨叶；2）将桨叶根部固定，桨叶不发生运动，在桨叶导边前方施加平行来流，来流速度为4.5 m／s，从桨叶导边向随边方向流动；3）调整桨叶的位置，使桨叶的叶面尽量与来流方向平行；4）桨叶与流体网格较大，以便能在短时间内模拟出三维的漩涡发放现象。

图9给出了简化后的流固耦合有限元模型，其中深色网格为桨叶固体网格，浅色为流体网格。

通过数值模拟固定桨叶的旋涡发放，给出流体单元的压力云图如下。

1）桨叶后方的单个旋涡

如图10所示，本文的数值模型和计算方法能够模拟出固定的螺旋桨桨叶在一定来流速度下的旋涡生成现象，而且展示了桨叶后部旋涡的拉扯破裂情况。经分析得出，流体在流经桨叶叶面达到桨叶随边时，产生了旋涡发放，生成一个旋涡。由于桨叶根部结构与桨叶顶部结构的偏转方向不同（图9中右图），桨叶根部对应的旋涡下部与桨叶顶部对应的旋涡上部的运动方向相反，整个旋涡在继续运动过程中，使得其因自身不同部分的拉扯而最终导致旋涡断裂。这一点与圆柱形状结构的旋涡发放明显不同。

2）桨叶后方的2个旋涡

随着时间的推移，桨叶后方不断产生新的旋涡。如图11所示，当桨叶后方存在2个旋涡时，一方面，这两个旋涡在运动过程中可能会发生接触，从而产生相互作用，两个旋涡形状也会发生一定的变化；另一方面，由于每个旋涡都有各自的流体力学特性，因此有的旋涡是在上下不同部分的拉扯下发生破裂，有的旋涡则不会破裂，或者破裂之后由于吸力还能够再次组合成一个完整的旋涡。

3）桨叶后方的3个旋涡

如图12左图所示，从图中可以清晰地观察到3个独立的旋涡，这表明本文的模型已经完全可以模拟出固定桨叶后方的旋涡，从开始生成一直到移动出计算模型边界的整个运动过程，给出了桨叶后方流场的实时变化情况。在右图中，3个独立的旋涡在运动过程中发生了接触，进而相互之间的边界不再清晰，而且导致3个旋涡的流体压力整体下降。这可以从能量守恒的角度加以解释，即由于在旋涡的相互接触过程与流场的紊乱过程中都消耗了旋涡的能量，因此从流体的压力云图上便显示出旋涡压力的下降，即从高压力变为低压力。

接下来，将对桨叶在流体作用下的内部结构力学性能进行分析。

1）桨叶叶梢与桨叶根部的响应比较

选取节点进行速度响应比较，数值模拟结果如图13所示 （Z方向速度即为垂直来流方向的速度）。其中，曲线A代表桨叶叶梢节点的数值模拟结果，曲线B代表桨叶根部节点的数值模拟结果。

选取单元进行内部应力响应比较，数值模拟结果如图14所示。其中，曲线A代表桨叶叶梢单元的数值模拟结果，曲线B代表桨叶根部单元的数值模拟结果。

根据数值模拟结果，在本文的算例中，桨叶叶梢节点的速度响应要大于桨叶叶根节点的速度响应，这主要是因为桨叶的叶根部分比桨叶的叶梢部分更靠近约束条件，而且叶根的厚度也比叶梢的厚度大得多，因此，在桨叶旋涡发放产生的涡激力作用下，叶根部分的节点难以发生运动。

对等截面的悬臂梁来说，靠近约束的位置其单元应力较大；而对于变截面的悬臂梁来说其单元应力与梁的厚度有很大关系，厚度大的位置其单元应力较小。本文的桨叶模型近似于变截面的悬臂梁，但厚度大的位置却靠近约束，所以，尽管图10的数值模拟结果显示桨叶叶根部分的单元应力比叶梢部分更大，但这只是针对本文的算例而言，对于具体的桨叶，需具体计算才能得到正确的结论。

空白组大鼠踝关节表面光滑，滑膜细胞无增生，周围未见炎症细胞浸润。模型组大鼠关节间隙变窄，关节软骨退变，局部可见骨质破坏，滑膜表面粗糙，滑膜细胞增生、层次紊乱，滑膜及周围软组织可见大量炎症细胞浸润，以淋巴细胞、巨噬细胞浸润为主，可见丰富的血管翳形成，纤维母细胞和脂肪细胞增生活跃。与模型组比较，各给药组大鼠关节软骨退变、骨破坏及滑膜细胞增生明显减轻，血管翳形成较少，炎症细胞数量减少，仅见少量纤维母细胞及脂肪细胞增生，详见图2。

2）桨叶随边与桨叶叶面中部的响应比较

选取桨叶随边和叶面中部的节点进行速度响应比较，数值模拟结果如图15所示。同时，选取桨叶随边和叶面中部的单元进行有效应力响应比较，数值模拟结果如图16所示。其中，曲线A代表桨叶随边的数值模拟结果，曲线B代表桨叶中部的数值模拟结果，选取的节点和单元至桨轴的距离基本相同。

图中线条表明，由于选取的两个节点至桨叶根部约束的距离基本相同，而桨叶随边的厚度又小于桨叶中部的厚度，所以，桨叶随边的节点速度响应要比桨叶中部的大。另一方面，由于流体流经桨叶叶面时，在桨叶随边产生了旋涡，旋涡脱落产生的涡激力大部分作用在桨叶的随边，而且桨叶随边的厚度要小于桨叶中部的厚度，因此，桨叶随边单元的有效应力响应大于桨叶中部单元的有效应力响应。

5 结 论

综上所述，经过一系列的假定和简化，采用数值模拟的方法研究螺旋桨固定桨叶的旋涡发放现象是完全可行的，并且观察得到了3个主要结论：

1）螺旋桨桨叶根部与顶部结构的偏向不同，在流体流动到桨叶随边产生旋涡发放时，旋涡上部与下部的运动方向也不相同，因此，部分旋涡会在运动过程中发生破裂，但旋涡破裂是否发生、破裂发生后旋涡是否还会重组，这与旋涡本身的流体特性有关。

2）独立的旋涡在逐渐远离桨叶的过程中，容易与其他旋涡发生碰撞而导致旋涡形状的改变和旋涡边界的破坏，以至整个流场的紊乱，这些都消耗了旋涡的能量，最终降低了旋涡的压力。

3）与桨叶叶根部分相比，桨叶叶梢部分的节点速度较大。但对于不同的桨叶，需具体计算才能比较桨叶叶梢与叶根部分的内部应力大小。与桨叶中部相比，桨叶随边部分的节点速度较大，单元有效应力也较大。桨叶内部有效应力的主要影响因素在于选取的单元所在位置的桨叶厚度，以及单元至桨叶根部的距离和单元至旋涡脱落位置的距离。

上述研究只是螺旋桨“唱音”问题的基础性工作，要从数值模拟的角度完整研究螺旋桨“唱音”的机理和各种影响因素，还需从以下几个方面进行努力：

首先，螺旋桨“唱音”现象的实质是桨叶后方的旋涡发放频率接近于桨叶的固有频率，因此，通过数值模拟的方法求出旋转中的螺旋桨桨叶后方旋涡的发放频率是问题的关键。由于受条件所限，本文仅给出了定性的分析。其原因在于，螺旋桨桨叶的随边厚度直接影响着旋涡发放的频率，而通常桨叶随边的厚度相对于整个螺旋桨的大小来说非常小，而在流固耦合数值模拟过程中，又必须要以随边的厚度来决定流体和固体网格的大小。因此，要计算得出较接近于实际情况的桨叶旋涡发放频率，整个数值模型就需要用大量的单元来进行描述。但完成这项工作离不开大型计算设备的支持，同时对研究者的模型简化能力也是一种考验。

其次，本文采用的简化模型是将桨叶固定，在桨叶导边前方施加流速来模拟螺旋桨转动时桨叶和流体的相对运动。这种方法虽然反应出了桨叶涡激共振问题的本质，但相对于完整的螺旋桨旋转过程来说，还不能体现各个桨叶后方所产生的旋涡之间的相互作用现象。

最后，目前螺旋桨“唱音”问题的相关研究都只考虑了螺旋桨转速与“唱音”现象之间的关系，没有加入螺旋桨前方来流的因素。而实际上螺旋桨在静水中旋转和在不同来流下旋转，桨叶随边后方的旋涡发放方向与频率都会发生变化。也就是说，带入旋涡发放频率计算公式中的速度v应该是一种合速度，即螺旋桨旋转引起的桨叶随边与静止流体的相对速度和前方来流速度的合速度。虽然由于螺旋桨桨叶形状不规则，难以从理论上计算出合速度的大小和方向，但在研究螺旋桨“唱音”问题时，在数值模拟过程中，必须考虑前方来流与螺旋桨转速的共同作用，才能接近于真实情况。

［1］SEOl H，JUNG B，SUH J C.Prediction of non-cavitating underwater propeller noise ［J］.Journal of Sound and Vibration，2002，257（1）：131-156.

［2］KEHR Y Z，KAO J H.Numerical prediction of the blade rate noise induced by marine propellers［J］.Journal of Ship Research，2004，48（1）：1-14 .

［3］于大鹏，赵德有.螺旋桨鸣音的混沌动力特性研究［J］.振动与冲击，2009，28（12）：47-52.

［4］李洁雅，欧礼坚，王志勇.对螺旋桨引起的谐鸣振动的探讨与实践［J］.中国修船，2001，15（6）：26-27.

［5］华汉金.螺旋桨鸣音产生机理及防治方法 ［J］.船舶，2002（2）：20-23.

［6］刘丹，陈凤馨.CFD在计算船舶螺旋桨敞水性能中的应用研究［J］.现代制造工程，2010（4）：18-20.

［7］赖华威，刘月琴，吴家鸣.基于CFD方法的螺旋桨性能计算与分析［J］.船海工程，2009，38（4）：131-135.

［8］胡健，黄胜，王培生.螺旋桨水动力性能的数值预报方法［J］.江苏科技大学学报，2008，22（1）：1-6.

［9］蔡荣泉，陈凤明，冯学梅.采用Fluent软件的螺旋桨敞水性能计算分析［J］.船舶力学，2006，10（5）：41-48.

［10］李裕春，时党勇，赵远.ANSYS 10.0／LS-DYNA 基础理论与工程实践［M］.北京：中国水利水电出版社，2008.