基于面元-涡粒子法的螺旋桨气动特性及噪声研究

刘乾, 刘汉儒, 李家辉, 尚珣, 陈南树, 王掩刚

1.西北工业大学 动力与能源学院, 陕西 西安 710072;2.中国空气动力研究发展中心 气动噪声控制重点实验室, 四川 绵阳 621000

随着航空业的发展,对飞行器的高效、节能、低噪声要求越来越高。涡扇发动机对燃料的利用效率约40%,而分布式电推进系统对电能的利用率能够超过70%[1-2],小型螺旋桨在低速设计工况下能长期保持较高气动效率。要使分布式推进系统的噪声满足标准,需要使用非定常气动及声学仿真手段进行预测,并采取降噪措施处理使之满足适航规定。螺旋桨气动噪声主要来源于空间中的非定常压力脉动和雷诺应力[3],为了实现叶轮机械噪声性能的预测、评估并提出改进的方向,必须对叶轮机械开展非定常数值计算。

面元法是目前计算精度最高的势流方法[4-6],Baltazar等[7]使用面元法对导管螺旋桨性能进行计算,发现尾涡计算准确性对性能预测极为重要。为了准确模拟尾迹作用,涡粒子法被研究人员用作尾迹模型。涡粒子法采用离散的拉格朗日粒子求解涡量输运方程[8-9],数值耗散更低,能够准确模拟尾迹的发展情况[10],将面元法和涡粒子法结合可以实现计算速度与精度的统一。Willis等[11]使用面元-涡粒子法对翼身振动问题进行了求解,并且获得了准确的结果。胡昊等[12]使用面元-涡粒子法对风力机气动特性进行计算,发现面元-涡粒子法能够给出风力机的其他流动细节。

在上述研究中,将面元法与涡粒子法结合进行叶轮机械气动性能模拟被认为是具有研究价值的工作,但将该方法用于气动噪声预测的研究有限。面元法高效的非定常数值模拟性能,使之在声学预测中具备独特优势,把面元-涡粒子法与气动声学模型结合,对叶轮机械气动及声学性能的预测有较高的工程实用价值,适于叶轮机械的多工况、多学科设计。

本文将分别介绍面元法、涡粒子法、Lowson频域声学模型以及面元-涡粒子法的理论,结合有限体积法验证面元-涡粒子法的仿真精度,并深入对比分析2种方法的数值计算结果,进一步揭示螺旋桨气动特性和声学特性。

1 理论方法

1.1 面元法

面元法以流体连续性为理论基础,求解三维流场中势函数:

2φ=0

(1)

使用Green函数求解(1)式,并引入源、汇、偶极等作为势流单元和不可穿透固体表面S,则在S上的P,Q两点间存在:

(2)

根据固体壁面、尾迹面元不可穿透边界条件,无限远处势函数诱导速度为0假设

(3)

可使用偶极和源表示速度势,从(3)式得到:

本文中使用涡线段等效面元计算诱导速度。如果涡线段与场点较近会产生数值计算的奇异性,因此引入涡环半径r减弱奇异性

式中:R1是涡线段起点和场点间的空间向量;R2是涡线段终点和场点间的空间向量;R0是涡线段起点和终点间的空间向量;h是涡线段和场点间的垂直距离。从(6)式能看出,当h过小时,r的存在将会减弱诱导速度数值奇异性;计算较远处诱导速度时,r/h接近0,对计算精度基本无影响。

为了保证面元能够准确捕捉前缘、尾缘处的流动情况,应该有效加密面元。为了更好地满足Kutta条件,把圆尾缘处理为尖尾缘,并使用3次样条插值方法生成离散坐标点。

本文使用余弦方法对螺旋桨的径向进行面元划分,控制吸、压力面的周向和径向离散节点数量一致,使用余弦方法离散可以保证前、尾缘处的面元密度更大,能准确描述叶片曲率。螺旋桨吸力面或压力面的径向、弦向面元数分别为NR和NC,所以单个螺旋桨叶片可被划分为2×NR×NC个面元。以NACA0012二维标准翼型示意,使用余弦方法划分面元[13]。

图1 NACA0012的余弦离散

1.2 涡粒子法

涡粒子法使用离散的拉格朗日粒子求解涡量输运方程。对不可压Navier-Stokes方程取散度后可以得到涡量输运方程

(7)

使用涡粒子离散空间中的涡量,设第i个涡粒子的强度和位置分别为αi和xi,则涡量场离散为

(8)

式中:涡粒子强度αi是涡量和体积的乘积;ξσ是使涡量光滑分布的磨光算子(光滑核函数),普遍使用高斯形式的核函数(也称截断函数)[14]

(9)

使用涡粒子模拟流场的发展时,粒子强度因为旋涡拉伸和黏性耗散发生变化,因此涡粒子法的求解过程分为位置和涡强两部分之间的迭代

(10)

(11)

(11)式右边由旋涡拉伸项和黏性扩散项构成,分别使用核函数和涡强交换方法(PSE)对2项求解,得到

使用四阶Runge-Kutta算法[15]求解(12)～(13)式。(12)式中ρ=|x-xj|/σ是无量纲长度参数,K(ρ)是用于计算诱导速度的核函数。

光滑核函数直接决定空间中涡量的分布情况,所以核函数的选择会显著影响数值计算精度。根据矩特性判断光滑核函数的精度,对r阶精度的光滑核函数有如下约束[16]:

(16)

给出几组不同精度的光滑函数和与之相对应的Green函数，如表1和图2所示。

表1 常用的光滑函数和对应的Green函数

图2 3种不同精度的Green函数

在表1常用的光滑函数和对应的Green函数中有

(17)

如图2所示,随着阶数的增加,涡量分布越集中,数值计算的误差也越小;Fun2和Fun3能把涡量集中在相同范围内,但Fun2比Fun3阶数更低,计算量更小,所以使用二阶精度的Gaussian分布函数计算诱导速度。

1.3 Lowson声学模型

叶轮机械中压力脉动有较强的周期性,使用时域方法计算声压要求输入多个时间步的压力脉动,计算成本较高;而频域方法仅需要计算叶片特征频率谐频上的声压,所以使用频域下的Lowson声压表达式[17]

cn=an+ibn=

(18)

将时间微分项替换,引入周向和轴向力分量,分部积分后得到Lowson公式

(19)

Fi=(-T,-Dsinθ,Dcosθ)

(xi-yi)=(x,y-Rcosθ,-Rsinθ)

为了提高适用性,将片条理论假设下的声压公式变换到三维坐标系[18]

对第n阶旋转偶极子声源的自由远场声压表达式有

(22)

1.4 面元-涡粒子法

使用简化的Hess等效原则,实现面元到涡粒子的转化。为了满足尾缘的Kutta条件,添加一列缓冲尾迹面元作为面元和涡量的过渡,并将其转化的涡粒子靠近尾缘边等效为涡线(图3尾迹面元转化为涡粒子中红色线段),其余3条边转化为涡粒子:

图3 尾迹面元转化为涡粒子

在求解过程中,涡粒子产生的诱导速度作为自由来流速度的一部分,用于求解面元上的源汇分布;面元对各个涡粒子产生诱导速度和速度梯度,作为自由来流一部分参与运算。具体迭代过程见图4。

CFD计算结果与声学模型耦合需要使用非定常数值计算方法获得时域下叶表压力脉动,对时域压力脉动做时频变换,将频域声压脉动作为声源项输入声类比模型,流程如图5所示。

图4 面元法和涡粒子法的耦合过程

图5 CFD和声学模型耦合

2 数值计算结果

2.1 前处理参数设置

在轴向来流速度15 m/s、转速为5 000 r/min的工况下(进距比为4.643),对螺旋桨进行仿真。

2.1.1 有限体积法网格尺度

流场分为旋转域和外域,旋转域半径取螺旋桨半径的1.5倍,外域取螺旋桨半径的10倍。使用URANS和声类比理论做叶轮机械气动噪声分析时,选择3BPF作为分析频率的上限,采样频率要高于6BPF使之满足奈奎斯特采样定理。

3BPF对应的声波波长是0.243 m,为了能够严格满足点声源假设,网格尺度需要小于1/6声波波长,所以网格尺度应该小于0.040 5 m。以表2中的网格尺度参数作为基准。

表2 网格大小设置参数 m

划分多面体网格如图6所示,为了验证该算例的准确性,需要进行网格无关性验证:

1) 边界层网格高度

在表2所示算例中,控制边界层厚度为7×10-4m,网格层数为15层,保证近壁面处网格高度小于10-5m。改变近壁面网格高度后,求得推力相对误差小于0.1%,说明边界层网格高度合适。

2) 尾迹区网格密度

尾迹区域中,在原网格基础上叠加5%基础尺寸网格,求得推力为11.474 N,与基准算例11.437 N相对误差为0.32%,认为尾迹区网格密度满足要求。

图6 旋转域网格示意图

分别对边界层、叶表、尾迹区进行网格加密并求解,证实使用该算例网格可以获得准确结果。

2.1.2 面元离散尺度

结合图1所示的面元余弦离散方法,对螺旋桨叶片的吸、压力面分别进行弦向和展向50×20,50×30,50×40,50×50共4种方式离散,并探究不同面元离散密度对数值仿真的影响。

表3 不同面元离散密度和有限体积法求得推力

图7 螺旋桨叶片面元离散示意图

如图8所示,面元离散密度与计算精度有如下关系:

1) 随着展向面元离散密度增加,面元-涡粒子法求出的推力逐渐接近有限体积法计算结果;

2) 当面元密度过高时,面元-涡粒子法的推力偏离有限体积法计算结果;

3) 面元长宽比接近1时,计算结果更准确。

上述现象说明,为了提高计算精度,需要把面元离散密度控制在合适的范围内,弦向50×展向40的离散密度能够保证计算准确。

图8 不同面元离散密度下的推力

2.2 气动性能

2.2.1 推力对比

选择轴向来流风速为15 m/s,转速为4 000,4 500,5 000,5 500,6 000 r/min的工况进行仿真,可获得面元-涡粒子法和有限体积法推力。

如图9所示,使用面元-涡粒子法与有限体积方法获得的推力趋势一致;同时在4 000 r/min时误差为8%,在4 500 r/min以上时推力误差降低。

图9 不同转速下螺旋桨推力大小

2.2.2 叶表压力分布对比

选择转速为5 000 r/min,轴向前飞速度15 m/s的工况。在25%,50%,75%叶高处,对比面元-涡粒子法和有限体积法解出的压力分布。

结合图10,对比不同叶高处的静压分布发现:

1) 在前缘滞止点处和吸、压力面大部分区域内,2种方法能获得基本一致的静压。

2) 在距离压力面前缘0.1倍弦长位置处,使用面元-涡粒子法算得静压偏低。在无黏假设下,面元法求出的气体加速作用更明显,而真实气体存在黏性,在叶背处加速作用稍弱,产生静压计算误差。

图10 25%叶高处螺旋桨压力分布对比

2.2.3 尾迹流场

从图11a)中可以看出,螺旋桨尾迹区域有明显的周期性叶尖涡存在。赵寅宇[19]证实在叶尖涡粒子对诱导速度计算有主要作用,而在图11a)中也能看到叶尖涡强度比桨盘其他展向位置的涡量要更高。图11b)是在尾迹区内周向速度大小,从图中可以看出,每经过高涡量区域会产生轴向加速作用。

图11 垂直方向尾迹区的涡量和速度分布

为了进一步揭示流场信息,下面对z轴坐标为-0.4 m处的切面上涡量和诱导速度进行分析。

为了提高数值计算收敛性,在使用有限体积法计算时考虑桨毂的影响。从图12可以看出,螺旋桨叶尖半径区域的尾迹区涡量仍然较高,呈现对称分布,与面元-涡粒子法中的1-1旋涡对应,在桨毂半径附近区域形成了较小的旋涡,也即是3-3旋涡。

图12 出口水平切面上的尾迹涡量

为了进一步验证2种方法解出的尾迹速度分布一致性,提取对角线(图12中黑色虚线)上速度分布，如图13所示。

从图13中可以看出,面元-涡粒子法和有限体积法解出的速度分布均呈现出“M”形,使用上述方法能够获得一致的尾迹速度分布。

图13 距出口轴向0.4 m处横截面轴向速度分布

2.3 声学性能对比

在本节中将结合各倍频下的声压级,对比面元-涡粒子法和有限体积法分别与Lowson声学模型结合求出的声学结果,分析螺旋桨声学特性以及不同方法之间的差异。取垂直于螺旋桨桨盘平面布置观测点,观测点半径为螺旋桨半径的5倍。

图14 不同叶片通过频率下的总声压级

结合图14可以看出,使用2种不同的方法解出的总声压级分布形式具有如下特征:

1) 使用2种方法解出的指向性一致,均呈现出以螺旋桨桨盘面为对称面,“8”形的分布形式,具有典型的偶极子声源分布特征[20];

2) 随着倍频的不断增加,2种方法预测出的声压级均有明显降低,但声压级指向性特征未发生改变,仍然轴向声压级最高,而在桨盘平面上最低。

上述现象说明,面元-涡粒子法能够较好捕捉到声压级指向性特征,有效预测螺旋桨噪声。

2.4 计算效率对比

为了评估面元-涡粒子法的计算效率,对比分别使用有限体积法与面元-涡粒子法计算相同时间步的耗时情况。数值计算在同一台服务器上完成,服务器的CPU型号是英特尔锐龙Gold 6132,该CPU运算主频是2.6 GHz,运算性能参数对比如表4所示。

表4 计算耗时对比

对比单核运算耗时可以发现,面元-涡粒子法计算效率明显高于有限体积法。面元-涡粒子法能有效提高螺旋桨的非定常气动和声学仿真效率,节省计算资源,在设计阶段具有较高的应用价值。

3 结 论

本文发展了基于面元-涡粒子法的叶轮机械非定常气动及噪声的快速预测方法。使用该方法与二阶有限体积法预测某型螺旋桨非定常气动性能,发现在推力、叶表压力以及尾迹速度分布等气动方面和声学方面具有较好的预测精度,研究结论如下:

1) 对比不同面元离散密度的计算结果,发现需要将面元尺寸控制在合适的范围内,使之同时满足能准确描述叶表几何尺寸和诱导作用准确的需求。使用面元-涡粒子法可以准确获得螺旋桨的推力、叶表压力以及尾迹速度分布,在不同叶高截面处压力分布最大相对误差为5%以内;

2) 涡粒子法在处理尾迹涡强发展时具有低的数值耗散优势。与二阶有限体积法相比,使用面元-涡粒子法可以获得明显的叶尖涡诱导作用,并且由于不存在转静交界面的数值耗散,尾迹涡强连续性也随之提高;

3) 使用面元-涡粒子法可以获得与有限体积法指向性一致的声压级分布,1BPF下前向辐射角60°的声压级误差为5%,该特征角度下,2BPF下误差为8%。除此之外,面元-涡粒子法的计算效率更高,将该方法与声学模型混合适用于螺旋桨非定常气动噪声性能的快速评估和进一步优化设计。