调距桨调距过程中叶元体受力的CFD分析方法研究

杨琼方，王永生，刘承江

（海军工程大学机械工程系 舰船新型推进技术研究室,武汉430033）

调距桨调距过程中叶元体受力的CFD分析方法研究

杨琼方，王永生，刘承江

（海军工程大学机械工程系 舰船新型推进技术研究室,武汉430033）

对调距桨作定距桨工况和作调距桨工况时叶元体水动力进行了CFD计算与分析。叶元体采用NACA4412翼型。采用κ-ε和SST两种湍流模型求解叶元体攻角为6.4°时压力系数分布并与面元法计算结果和实验曲线进行比较，校验了CFD分析方法的可信性。在此基础上对叶元体攻角每改变2°进行一次稳态工况计算，分析调距桨在60°转角范围内叶元体水动力随攻角的变化关系。变化曲线与实验曲线吻合良好，且计算精度要高于面元法程序XFOIL。然后采用滑移网格技术对调距桨整个调距过程进行了动态模拟，并将计算结果与作定距桨使用时的稳态工况值时进行了对比。结果表明：在攻角-12°～18°范围内，叶元体所受水动力、力矩在稳态和动态时的值基本相等。动态时叶元体正负失速角均较稳态延迟2°～4°。在大攻角范围内，叶元体稳态时所受阻力明显较动态时大、升力和制动力矩均较动态时小。研究结果表明现行采用的利用稳态敞水实验值来确定调距机构强度的方法是存在风险的；建议取稳态值和动态值两者中大者为强度设计的基本依据。

船舶；调距桨；定距桨；叶元体；水动力；数值模拟

1 引 言

采用求解RANS方程对定距桨（Fixed Pitch Propeller,FPP）敞水性能（KT～J,KQ～J）进行数值分析是目前众多研究机构和大学将计算流体力学（Computational Fluid Dynamics,CFD）应用于螺旋桨设计与性能分析中的一个十分活跃的研究领域。文献[1]通过求解RANS方程对三叶定距桨在设计工况和非设计工况下的敞水性能进行了CFD模拟，计算所得推力系数KT、力矩系数KQ、速度和压力分布以及湍流量与实验数据均有较好的一致性。文献[2]利用RANS求解器对导管螺旋桨的不可压粘性流场进行了模拟，在不同进速系数J下得到KT、KQ和速度分布与实验吻合良好。作者所在研究室开展有关这类工作[3-6]的步骤通常是：（1）采用CAD软件（SolidWorks,UG,Por/E,CATIA等）对定距桨进行几何建模；（2）依据物理过程特点选取流场控制体；（3）将几何模型导入网格划分软件（ICEM,TuoboGrid），采用适当的网格类型、网格拓扑结构和网格节点数对流场进行空间离散；（4）选择合适的湍流模型、空间离散格式、时间积分格式以及数值算法，设置相应的边界条件和初始条件，采用基于有限体积法的CFD软件Fluent或CFX进行数值计算；（5）计算结果的后处理，使流场“可视化”；（6）校核与验证，包括壁面y+值分析、湍流模型的影响、网格无关性分析等，并将计算结果与实验/试验数据进行比对分析；（7）CFD计算结果的分析与研究，从中汲取与研究内容相关的认知结果。

但是，利用CFD方法来分析调距桨（Controllable Pitch Propeller,CPP）调距过程中桨叶或其叶元体水动力和力矩的相关研究尚未见报导。调距桨桨叶绕调距轴匀速转动时桨叶水动力分析是一件富有工程意义的研究工作：（1）调距机构强度设计主要依据是桨叶上的水动力和离心惯性力[7]。桨叶固定在桨毂上某一位置时所受的水动力、力矩和调距过程中瞬时处于同一位置时所受水动力、力矩从理论上来说是不相同的。作者认为，强度设计应从桨叶稳态受力和动态受力中两者择其大者为外来载荷。（2）推进系统紧急加速过程、紧急倒车过程时调距桨水动力性能的CFD准确计算可以大大提高推进系统动态性能数值计算的精度。利用CFD来分析调距桨调距过程水动力性能的方法较现行惯用的利用调距桨敞水稳态性能计算上述过程动态性能的准稳态方法是一个较大进步。

本文对调距桨桨叶叶元体在调距过程中的受力进行了CFD分析。采用κ-ε和SST两种湍流模型计算了NACA4412翼型表面的压力系数分布并与面元法计算结果和实验值进行了比较，检验了本研究方法的可信性。在此基础上分析了桨叶定距时叶元体的升力、阻力和力矩大小及其随叶元体攻角变化的规律，然后采用滑移网格技术对调距过程中叶元体运动进行了动态模拟，并将动态过程中叶元体的升力、阻力和力矩随攻角的变化规律与桨叶在桨毂上相应位置固定时稳态受力（升力、阻力和力矩）进行了分析和比较。

2 理论基础

本文利用基于中心节点控制和有限体积法的计算流体力学程序来求解RANS方程。不可压粘性流体数值求解的控制方程如下：

其中，σk=1.0，Cε1=1.44，Cε2=1.92，σε=1.3。 选择 SST 湍流模型时，结合了 κ-ε 和 κ-ω 两种湍流模式，在自由流动区域使用κ-ε模式，而在近壁面区域（y+＜2.5 ）使用κ-ω模型中的低雷诺数公式，两者之间通过混合函数来过渡，这样可以不需要使用壁面函数，能够较好地模拟粘性底层的流动[1，8]。

3 数值计算

3.1 校核与验证

桨叶叶元体取为NACA4412非对称翼型，如图1所示。建立几何模型时坐标原点取于弦长中点。叶元体流场控制体取为14c×6c，c为弦长。采用结构化网格对流场空间进行离散。叶元体周围采用O形网格，以保证前缘和尾缘处的网格质量。外层采用L形网格。采用基于附加修正技术的多重网格用于加速求解。为了能充分求解边界层的流动，近壁面区域应布置至少10层网格[9-10]，第一层网格厚度应尽可能小，以使所有壁面y+在1～30之间[1]。叶元体网格结构及流体控制体尺寸如图2所示。网格节点数约18万。

3.2 定距时叶元体流体动力计算

本文计算分析针对调距桨目前主要采用的调距方式进行，即正倒车转换时通过零螺距位置，调距角范围约45°～60°（图9）。本研究时取60°转角。定义正车时攻角为正。定距时叶元体攻角每变换2°进行一次稳态工况时的受力计算。所有攻角计算均采用与校验分析时同样的流场控制体尺寸、网格拓扑结构、网格元素和节点数、初始条件、边界条件及其他相应的求解设置。计算得叶元体阻力系数、升力系数、力矩随攻角变化关系如图10所示。图10中同时给出了程序XFOIL的计算结果和NACA4412翼型在两种雷诺数下的实验升力特性曲线。经比较可得：CFD计算叶元体在攻角-12°～14°范围内受力随攻角变化曲线与面元法计算结果及实验曲线均吻合良好，攻角大于 14°时，CFD计算精度要高于面元法程序XFOIL。这也再次校验了上述分析方法的正确性。

叶元体升力为零时的攻角称为零升力角α0。理论和实验均表明，α0主要与翼型拱度有关。对NACA4412翼型，最大相对拱度=4%，即最大拱度为2c/5。依据NACA四位数字翼型零升力角的近似估算公式：α0=-100，即NACA4412翼型零升力角为-4°。图10所示的CFD计算结果中当升力为0时对应攻角正好为-4°，与经验公式估算值吻合。叶元体升力最大时对应攻角为失速角。此后，随攻角增加升力迅速下降，阻力急剧增加。由图10可知，在该雷诺数下，失速角位于 14°～16°。攻角为-2°时阻力系数最小。在攻角-12°～14°范围内，叶元体所受升力随攻角呈直线增加，阻力变化平缓，绕弦线中心俯仰力矩随攻角变化也呈线性关系。此时无明显流动分离，压差阻力很小，主要受粘性阻力作用。稳态计算过程表明：当攻角很大时，数值计算收敛性变差，如图11所示攻角为26°和-22°时升力系数与阻力系数收敛过程，这主要是由于叶元体流动大面积分离引起的，如图12所示攻角为26°和-22°时叶元体周围流体速度流线，分离涡清晰可见。这时湍流模型对流动变化剧烈的流体的适应性也会对计算结果产生一定影响。

3.3 调距时叶元体流体动力计算

调距时叶元体周围流体是非定常流动，对其要进行动态过程计算。对动态调距过程设计如下：恒速调距；从零螺距开始变到最大正车螺距；然后从最大正车螺距变到最大倒车螺距（负螺距）；再从最大倒车螺距（负螺距）变回到零螺距位置。调距过程的动态计算可利用动网格技术或滑移网格技术来实现。动网格技术是通过指定物体壁面网格节点的运动模式、周围流场网格节点通过自动插值实现网格变形从而进行实时模拟。利用动网格技术进行计算时，网格在变形过程中扭曲度变大，质量变差，因而壁面第一层网格难以满足厚度要求，即壁面y+值要求，从而使计算精度降低。叶元体变形过程中某时刻网格拓扑结构如图13所示。滑移网格技术是指通过设置旋转域和静止域之间的交界面来进行数据交换，各个子域的网格节点相对位置并不改变，如图14所示。本文采用滑移网格技术，利用CFX Expression Language（CEL）语言编程控制叶元体运动规律来实现动态模拟。内旋转域半径取为1.8c。嵌入程序模块因篇幅原因在此省略。

用上述方法计算得到叶元体升力系数、阻力系数和力矩随攻角变化关系如图15所示。由图15可知：零升力角同样是-4°，阻力最小时攻角也是-2°。失速角位于16°～18°，较稳态时稍有延迟。在攻角-16°～16°范围内，叶元体所受升力随攻角线性增加，阻力变化平缓，绕弦线中心俯仰力矩随攻角变化也呈线性关系。该范围较稳态工况计算值稍大。在零螺距到正最大螺距区间，正车和倒车过程中同一攻角位置时叶元体升力、阻力和力矩基本吻合，但处于零螺距到负正大螺距区间时，在大攻角-18°～-30°范围内，同一攻角位置正车时叶元体所受阻力较倒车时大，升力小，制动力矩也小。这主要是由于翼型的非对称引起的。攻角为-20°时正车和倒车过程中对应该瞬时叶元体速度流线如图16所示。正车过程中处于该攻角时叶元体流体分离要早，且范围要大，涡流更加明显。

4 对比与分析

将动态调距时叶元体的升力系数、阻力系数和力矩随攻角的变化规律与桨叶在桨毂上相应位置时稳态工况的升力系数、阻力系数和力矩合画在同一个坐标系中并进行比较，如图17所示。由图17可知：攻角-12°～18°范围内，叶元体所受升力、阻力和力矩在定距工况和调距工况时基本相同。攻角2°～18°范围内，叶元体稳态时所受升力较动态时稍大，力矩稍小。从攻角18°到正最大螺距和从攻角-12°到负最大螺距转角范围内，叶元体在稳态时所受阻力明显较动态时大、升力和力矩均较动态时小。从图18可以看到，现行采用的利用稳态受力情况来计算调距机构强度的方法从安全性上来说是存在一定风险的，也不符合实际情况。作者认为，应先对调距桨作定距桨工况时的受力和作调距桨工况时的受力进行计算，然后对其分析，择其大者作为设计强度计算的基本外界载荷。

5 结 论

本文通过求解RANS方程，对调距桨作定距桨工况和作调距桨工况时叶元体升力、阻力和力矩进行了计算与分析。有以下结论：

（1）作定距桨工况即稳态工况时，叶元体攻角为6.4°的CFD计算叶元体压力系数与实验值能很好地吻合。在雷诺数为1.0×106时，稳态计算得翼型零升力角为-4°，失速角为14°～16°，阻力系数最小时攻角为-2°，这些计算结果与经验公式估算值一致。

（2）作调距桨工况时，计算叶元体零升力角和阻力系数最小角与稳态工况时相同，正负失速角均较稳态时延迟2°～4°。在攻角-12°～18°范围内，叶元体所受水动力、力矩与定距工况时基本一致，均是阻力随攻角变化平缓，升力和力矩随攻角增加而呈线性增加。

（3）大攻角时，叶元体作定距桨工况所受阻力明显较作调距桨工况时大、升力和力矩均较作调距桨工况时小。由此说明稳态工况和动态工况叶元体受力是有差别的，这表明：过去惯用的利用稳态工况受力值来计算调距机构强度的方法是不合理的、是存在风险的。

（4）调距机构强度设计应从桨叶稳态受力和动态受力中两者择其大者。本文所做二维叶元体水动力计算正是调距过程中桨叶水动力分析的一个结果。叶元体在调距过程中水动力矩是大于稳态过程相应值的，现在设计部门惯用的将稳态工况受力值加上一定安全系数来进行强度设计会使桨毂直径过大、调距机构变重、使桨效率降低。利用CFD方法准确计算动态过程水动力大小能解决这一问题。

（5）本研究是按国内外调距机构常用速度（30秒完成从最大正螺距变至最大负螺距）来进行CFD计算的。显然，调距角速度的大小对叶元体水动力是有影响的。调距的角速度对叶元体水动力的影响尚有待进一步分析。叶元体大攻角时湍流模型对流体流动的适应性也有待进一步研究。

[1]Lucia Sileo,Aldo Bonfiglioli,Vinicio Magi．RANS simulation of the flow past a marine propeller under design and off-design conditions[C]//4th Annual Conference of the Computational Fluid Dynamics Society of Canada.Canada,2006.

[2]Sanchez-Caja A,Rautanhiemo P,Siikonen T.Simulation of incompressible viscous flow around a ducted propeller using a RANS equation solver[C]//23rd Symposium on Navy Hydrodynamics.Val de Reuil,France,2000:527-539.

[3]杨琼方,王永生,张志宏等．大侧斜螺旋桨敞水性能的RANSEs模拟[J]．湖南科技大学学报（自然科学版）,2008,23(4):66-70.

[4]杨琼方,王永生,张志宏等．喷水推进泵三维造型研究及其流体动力性能CFD分析[J].武汉理工大学学报（交通科学与工程版）,2008,32(5):921-924.

[5]杨琼方,王永生,张志宏等．叶片数对喷水推进器性能影响的计算流体动力学分析[J].机械工程学报,2009,45(6):222-228.

[6]Yang Qiongfang,Wang Yongsheng.Research on the optimum blades number of the mixed flow pump based on CFD[C]//International Conference on Waterjet Propulsion 5.London:RINA,2008.

[7]吴耀祖,邵中立,瞿守恒.舰船原理[M].武汉:海军工程大学,2004:216-218.

[8]Abdel-Maksoud M,Heinke H J.Scale effects on ducted propellers[C]//24th Symposium on Navy Hydrodynamics.Fukuoka,Japan,2002:744-759.

[9]Anders Smaerup Olsen.Investigation of the viscous flow around two model propellers in uniform inflow[R].Finland:Department of Naval Architecture and Offshore Engineering Technical University of Denmark,2001.

[10]Turbulence and Wall Function Theory[M].CFX-Solver Theory,2005:97-108.

[11]XFOIL 6.9 User Primer[M].Mark Drela,MIT Aero＆Astro.Harold Youngren,Aerocraft,Inc.,1989.

Research on CFD method for hydrodynamic forces of a CPP’s blade section

YANG Qiong-fang,WANG Yong-sheng,LIU Cheng-jiang

(Institute of Marine Propulsion Technology,Department of Mechanical Engineering,Naval University of Engineering,Wuhan 430033,China)

The hydrodynamic forces and torque of a CPP’s blade section are calculated and analyzed by CFD for both the FPP’s and CPP’s working conditions.The blade section is NACA4412 airfoil.The procedure of RANSEs simulation used in this research is introduced first.The pressure coefficient distribution of the element when the attack angle is 6.4 degrees is calculated with both κ-ε and SST κ-ω turbulence models and the results correspond well to the panel method rssolves and experimental data,indicating that the CFD method is validated.A steady state is calculated at every 2 degrees,which is used to analyze how the hydrodynamic forces and torque are changed while the attack angle change during the CPP’s whole working region.The transient process of changing the pitch is carried out using the sliding mesh,and the results are compared to the FPP working condition at the same attack angle.The comparison shows that the hydrodynamic forces and torque of the element agree well with each other for both steady and transient results when the attack angle lies in-12～18 degrees,but out of this region,the drag of the element is significantly bigger for the steady process,but the lift and brake torque are smaller;the attack angle when the lift of element attaches to the maximum is delayed about 2 to 4 degrees for the transient process.It is clear that the method used now to calculate the intension of the device for changing pitch based on the steady state force is risky.It is suggested that the larger value from either a steady state or a transient state should be a basic force to calculate the device’s intension.

ship;CPP(controllable pitch propeller);FPP(fixed pitch propeller);blade section;hydrodynamic forces;CFD

U664.34

A

1007-7294（2010）05-0458-08

2009-10-09

杨琼方（1984-），男，海军工程大学机械工程系博士生。

book=465,ebook=52