分布式粗糙前缘对NACA0012翼型失速特性的影响

李 艺， 白俊强， 张彦军， 赵 轲

(1. 西北工业大学 航空学院， 西安 710072； 2. 第一飞机设计研究院，西安 710089)

由于飞机机翼、螺旋桨桨叶和风力机叶片具有多样化的运行环境和长时间、多工况的使用条件，其翼型表面容易受到风沙、灰尘以及冰晶的堆积，或由于腐蚀及污染等现象导致表面粗糙度发生变化，进而影响机翼、桨叶的气动性能.研究粗糙度对翼型气动特性的影响对分析飞机和风力机的气动特性、改善风力机输出功率十分关键，也对固定转捩风洞试验有指导意义.

分布式粗糙带最初用于模拟翼型或平板的强制转捩流动[1].之后，Ramsay等[2-3]通过实验研究前缘粗糙对NREL-S系列翼型气动特性的影响，结果表明，相比光滑翼型，粗糙前缘减小翼型最大升力系数，增大阻力系数，降低零升俯仰力矩系数.Reuss等[4]对NACA4415翼型的粗糙度影响进行实验研究，得到相似的结果.Kerho等[5]对NACA0012翼型0° 迎角时前缘分布式粗糙带对边界层发展和转捩过程进行实验研究，相比光滑翼型，粗糙带使转捩位置前移，其转捩过程中的湍流强度比自然转捩小.包能胜等[6-7]对风力机翼型表面增加粗糙带进行初步理论计算和实验研究，分析其对翼型性能的影响，结果表明压力面后缘增加粗糙度后翼型升力系数和阻力系数均较光滑叶片有所增加，而增加前缘粗糙度则在不同翼型上表现出不同的影响效果.Li等[8]对比了光滑和粗糙DU95-W-180翼型的仿真结果，并对比了不同粗糙高度、粗糙带长度及安装位置对翼型气动特性的影响.研究发现随着粗糙高度逐渐增加，升力系数迅速减小，当粗糙高度达到临界粗糙高度后升力系数基本保持不变.在翼型后缘放置粗糙带时，随着粗糙高度增加，升力系数先减小再增大.李仁年等[9]根据叶片易形成沙粒堆积和风蚀，采用数值方法研究风力机翼型粗糙表面气动性能发现翼型后缘处附加粗糙带可以提高翼型升力系数.焦灵燕等[10]研究不同表面粗糙厚度对多个风力机叶片翼型的影响，数值模拟结果表明翼型升阻比均下降，而S822翼型更适合在低雷诺数风沙环境中使用.Joseph等[11]通过对粗糙表面DU96-W-180翼型进行风洞实验，发现随着粗糙度雷诺数的增加，翼型升力减小、阻力增大，转捩位置前移.粗糙元带来的Tollmien-Schlichting波快速增长导致转捩位置改变，阻力增加几乎完全由转捩区域前移造成，而非来自粗糙元本身的阻力增量.李虹杨等[12-13]对粗糙表面上的转捩流动进行模拟，并分析其对涡轮叶片流动转捩及传热特性的影响.结果表明粗糙度会显著改变湍流边界层壁面温度分布，且粗糙度对吸力面转捩位置影响是非线性的，当粗糙度大于某一临界值时转捩位置发生突变.Zhang[14]通过实验研究了GA(W)-1翼型表面不同尺寸、不同分布模式粗糙度的影响.结果表明，相比分布模式，粗糙度高度对升力和失速特征的影响更大.Kruse等[15]采用包含粗糙度修正的新k-ω(k为湍动能，ω为比耗散率)湍流模型对前缘粗糙的NACA633-418翼型进行数值模拟，并实验进行对比.结果显示计算的升力减小与阻力增大量与试验较为一致，但对锯齿形粗糙带的模拟精度不够，模型存在局限性.Wang等[16]研究了低雷诺数下表面粗糙度对压力机翼型气动性能的影响，重点分析不同粗糙区长度和粗糙高度对前缘分离泡结构及边界层发展的作用.结果表明，在过渡粗糙区域和完全粗糙区，粗糙度变化对翼型性能的影响不同.

1 数值模拟方法

采用基于三维可压缩雷诺平均Navier-Stokes方程的求解器进行流场模拟，控制方程形式为

∬∂ΩG(Q)·ndS

(1)

式中：t为时间；Ω为控制体；Q为守恒变量；V为控制体体积；∂Ω为控制体的边界；n为控制面的外法线单位向量；F(Q)为无黏通量；G(Q)为黏性通量；S为控制体表面积.求解器无黏通量的空间离散格式为二阶迎风Roe格式，黏性通量的空间离散格式为二阶中心差分.时间推进方式采用伪时间推进，算法为近似因子分解隐式推进.湍流模型采用Menter切应力输运模型.

1.1 全湍流动模拟方法

Menter切应力输运模型方程如下[17]：

(2)

式中：ρ为流体密度；ui为速度(i及下文j，k对应物理空间坐标下的3个方向)；xi为坐标轴分量；Pk、Pω、Dk、Dω为方程源项；ω为比耗散率；μl为流体黏性系数；μt为涡黏性系数；Cω、σω、σω2为模型常数；F1为混合函数.涡黏性系数μt由下式确定：

(3)

式中：a1为常数，取0.31；W为涡量绝对值；F2为耦合函数.

方程源项Pk、Pω、Dk、Dω的定义根据爱因斯坦求和约定，可简写为

(4)

(5)

湍流模型控制方程中的常数及函数的设置可参考文献[17]，本文主要关注其边界条件的改动.

在湍流模型中，原始的壁面边界条件为

(6)

式中：v为运动黏度；Δy为壁面距最近网格点间的垂直距离；β为常数，取0.09.

考虑粗糙壁面的影响，本文采用Wilcox[21]提出的湍流壁面边界条件修正.壁面比耗散率ωWall由下式得到：

(7)

(8)

1.2 转捩流动模拟方法

(11)

式中：Pθ t为方程源项；σθ t和cθ t为模型常数[19]；Reθ t为通过经验关系式进行求得的转捩雷诺数；Fθ t为边界层指示因子；tL为时间尺度，用来平衡量纲；U为速度.

间歇因子γ的控制方程为

(12)

方程的产生项和破坏项为

(13)

式中：ca1、ca2、ce1及ce2为常数；ε为应变率；函数Flength控制转捩区长度；函数Fonset判断转捩是否发生；函数Fturb在层流边界层外及黏性底层关闭再层流化源项[22].

由于层流分离泡区域能量耗散较大，在分离泡再附着时需要提高湍流动能的产生项，以使其完成再附着的数值模拟.因此需要对间歇因子进行修正：

γeff=max{γ,γsep}

(14)

(15)

式中:γeff为有效间歇因子；常数s1为2；Reν为涡量雷诺数；RT为黏性比.

转捩模型通过有效间歇因子γeff影响湍动能输运方程的源项来达到模拟流动中湍流增长的目的，利用γeff修改的k方程源项为

(16)

(17)

考虑粗糙度的影响，通过粗糙增长因子输运方程[20]确定转捩动量厚度雷诺数的数值，实现在壁面粗糙条件下对流动转捩及当地摩擦因数准确的数值模拟.

粗糙增长因子Ar的输运方程为

(18)

(19)

Pθ t, mod=

(20)

式中：b为混合函数.粗糙度的影响通过函数FAr体现，而混合函数b保证流动转捩后关闭粗糙度对方程源项的影响，函数形式如下：

(21)

FAr=

(22)

图1 平板网格示意图Fig.1 Grid of flat plate

式中：λθ为Thwaites压力梯度因子；函数F(λθ)采用Langtry等[19]提出的公式进行计算，他们对试验数据进行拟合，表征压力梯度因子λθ和边界层边缘处湍流度Tu的影响.

F(λθ)=

(23)

经过验证，该公式可以有效反映分布式粗糙元在自由转捩中的作用[18, 20].

2 模拟方法验证

由于湍流模型耦合粗糙边界的模拟已经非常成熟[23-25]，本文仅对转捩模型与粗糙增长因子的耦合效果进行验证.选择对零压力梯度平板及NACA0012翼型绕流进行模拟，验证该粗糙表面流动转捩预测方法的有效性.

2.1 零压力梯度平板

Feindt[26]曾对零压力梯度粗糙平板进行转捩流动实验，并给出实验中关于粗糙雷诺数Reks=ρuτks/μ及边界层雷诺数Rex=ρUinX/μl(Uin为平板流场域入口速度；X为沿边界层坐标自绕流物体前缘算起的距离)的关系，Langel等[20,27]均使用该算例对转捩模型进行标定，下文将引入他们的模拟结果进行对比.

Uin为20 m/s，湍流度Tu为1%，基于弦长的雷诺数Re为1.3×106.采用7个不同等效粗糙高度ks进行模拟，分别为0、90、130、210、245、275以及 320 μm，其中0对应光滑平板.

选取长度L为1 m的平板按实验给定的雷诺数进行模拟计算，平板网格示意图如图1所示.网格分布及边界条件布置也在图中展示.流场沿x轴范围为 -0.2～1.5 m，z轴范围为0 ～0.22 m.流场入口距离平板0.2 m，入口速度为20 m/s，且等速度分布.平板上方网格量为301×101，壁面法向网格距离增长率为1.1，壁面第一层网格无量纲高度y+为0.3.

不同粗糙高度下转捩位置和摩擦因数Cf的对比如图2所示，其中Rext为转捩位置处的边界层雷诺数.引入 Dassler 等[27]和Langel等[20]等的数值模拟结果进行对比.从结果看，随着粗糙高度增大，转捩位置逐渐向上游移动，且转捩后全湍区域的摩擦因数稍大于光滑表面转捩后的摩擦因数.当粗糙度较小时，如ks<90 μm时，粗糙高度对流动的转捩位置影响较小，当粗糙高度到达到临界值后，其对转捩的影响迅速扩大.该粗糙转捩模型较好地反映了这一特征，与实验数据吻合.

图2 表面摩擦因数分布及转捩位置Fig.2 Skin friction coefficient distributions and transition locations

2.2 NACA0012翼型

NACA0012翼型粗糙前缘实验由Kerho等[5]在伊利诺伊大学的亚音速风洞中完成.本算例的自由来流马赫数为0.1，基于弦长的雷诺数为1.25×106.翼型前缘在 0.006 112c～0.025 8c(c为翼型弦长)位置处粘贴粗糙带.粗糙颗粒高度为350 μm.模拟迎角为0° ，湍流度Tu为0.1%.

表1 用于收敛性分析的不同网格尺寸Tab.1 Different mesh sizes in mesh convergence analysis

采用粗、中、密3套网格进行数值计算，网格点布置如表1所示，远场到翼型的距离为100倍弦长.密网格沿翼型周向网格点数为830个，沿法向布点160个，翼型法向网格增长率为1.1，第一层网格高度y+约为0.3.而粗网格的网格量约为细网格的1/4.具体网格分布信息如表1所示.图3所示为使用3套网格计算光滑翼型0° 迎角时表面摩擦因数分布对比.可以看到3套网格计算的结果差别较小，预测的转捩位置基本一致，仅在摩擦因数峰值处有细微差别.考虑到计算量及模拟精度，并且为了保证在转捩位置处流向网格间距小于0.003c，选择中等网格进行后续数值计算，其分布如图4所示.

图3 光滑翼型0° 迎角表面摩擦因数分布Fig.3 Distribution of surface friction coefficient at an angle of attack of 0° for smooth airfoil

图4 NACA0012翼型中等网格Fig.4 Medium grid for NACA0012 airfoil

图5所示为翼型上表面4个不同流向站位的速度型分布图，站位分别在x/c=0.2，0.3，0.4，0.5处，图中：h为距翼型表面的法向距离，u为当地流向速度，U∞为来流速度.将光滑表面、粗糙表面翼型模拟结果与实验结果[5]进行对比.对于光滑翼型，4个站位处均保持典型层流速度型分布，流动转捩发生距前缘更远的位置.粗糙带后，转捩发生位置向上游大幅移动，数值模拟结果捕捉到该现象，与实验吻合较好.由于数值模拟预测的转捩位置比实验结果稍向下游移动，因此在x/c=0.2处出现少许偏差.

图5 不同流向站位速度型对比Fig.5 Comparison of velocity profiles at several streamwise locations

上面两个算例表明该粗糙表面转捩预测方法能够比较准确地反映粗糙表面对流动转捩的影响，具有较好的模拟精度，可以将方法应用于下文的算例计算分析.

3 NACA0012翼型数值模拟结果及分析

使用NACA0012翼型作为研究模型，自由来流马赫数为0.086，基于c=1 m的雷诺数为2.0×106.在翼型前缘10%弦长以前的位置上粘贴有粗糙带，颗粒的等效粗糙高度为0.001 m和0.003 m.前缘粗糙翼型的几何如图6所示，其前缘红色部分为粗糙表面.网格则采用上一节中翼型的中等网格.

图6 NACA0012翼型粗糙表面分布Fig.6 Rough surface distribution on NACA0012 airfoil

3.1 翼型绕流全湍模拟分析

对光滑翼型(C1)和两个前缘粗糙翼型进行数值模拟，两个前缘粗糙翼型(分别标记为R1和R2)的等效粗糙高度为0.001 m和0.003 m.计算得到各翼型的升力系数CL和阻力系数CD随迎角α变化如图7所示，各翼型的最大升力系数CLmax及失速迎角α′如表2所示.

从图7可以看出，翼型前缘有粗糙带后，翼型的最大升力系数大幅下降，阻力系数增大，失速迎角显著减小.由表2可知，相比C1翼型的最大升力系数，R1翼型减小30.59%、R2翼型减小33.63%；相比C1翼型的失速迎角，R1翼型和R2翼型均减小5°.

表2 各翼型最大升力系数及失速迎角Tab.2 CLmax and stall angle of airfoils

各个翼型失速迎角附近的流场流线图如图8～10所示.在全湍流动状态下，光滑NACA0012翼型首先在后缘出现分离，之后分离随迎角增大向上游延伸，分离区逐渐扩大，直到翼型因上表面出现大面积分离而失速，这是典型的翼型后缘失速状态.而在前缘附加粗糙带后，翼型的失速过程基本保持一致，仍为后缘失速.但是相比光滑翼型，前缘粗糙翼型在12° 迎角附近即发生后缘分离，分离区出现的迎角减小，这说明边界层维持流动稳定的能力减弱.

图7 不同粗糙度翼型升力系数、阻力系数对比Fig.7 Variation of CL and CD for airfoils at different roughnesses

图8 C1翼型失速迎角附近流线图Fig.8 Streamlines and contours of streamwise velocity for C1 airfoil near stall angle

图9 R1翼型失速迎角附近流线图Fig.9 Streamlines and contours of streamwise velocity for R1 airfoil near stall angle

为进一步分析全湍流动中粗糙前缘对翼型气动特性的影响规律，提取表面压力、摩擦力分布及绕流边界层的参数进行对比.图11所示为10° 迎角翼型的压力系数Cp和表面摩擦因数分布.

从压力系数分布图可以看到，在10° 迎角时，粗糙翼型与光滑翼型的压力分布差别较大.当翼型前缘附加粗糙带后，前缘的负压峰值降低，且粗糙高度越高，负压峰值越低.再对比表面摩擦因数分布，在翼型前缘附加粗糙带使表面摩擦因数迅速增大，而当流体经过粗糙带后，表面摩擦因数迅速下降.

截取翼型流向不同弦长站位处法向流场变量分布进行对比，变量包括当地速度、湍动能和涡黏性，对比结果如图12～14所示.

图11 α=10° 翼型压力系数分布和翼型表面摩擦因数分布Fig.11 Pressure coefficients and surface friction coefficients of airfoils at α=10°

图12 α=10° 翼型不同站位速度型分布对比Fig.12 Comparisons of velocity profiles at several streamwise locations at α=10°

图13 α=10° 翼型不同站位湍动能分布对比Fig.13 Comparisons of turbulence kinetic energy at several streamwise locations at α=10°

由图12可知，附加前缘粗糙带后，翼型边界层厚度显著增大，速度梯度减小，速度型分布变得陡峭.而且随着等效粗糙高度增大，沿壁面法向边界层内速度增长变得缓慢，边界层厚度增大.相比光滑翼型，粗糙翼型在后缘附近维持前凸的速度型分布更加困难，已经逐渐接近分离状态.流场流线图也显示，粗糙翼型在12° 迎角即形成后缘分离区.前缘粗糙带减小翼型开始出现后缘分离的迎角.

在湍动能和涡黏性分布对比图中，同一站位下，粗糙翼型的边界层湍动能和涡黏性更大，且随着等效粗糙高度的增大而增大.这是由于翼型前缘附加粗糙带，使翼型表面摩擦力增大，壁面造成的流体涡量增加，使得边界层中流体的动量交换和湍流输运增强，边界层内的湍动能及涡黏性增大.而这也会增加边界层内能量的消耗，使其克服逆压梯度的能力减弱，分离起始点向上游移动.

图14 α=10° 翼型不同站位涡黏性分布对比Fig.14 Comparisons of eddy viscosity at several streamwise locations at α=10°

3.2 翼型绕流转捩模拟分析

表3 各翼型最大升力系数及失速迎角Tab.3 CLmax and stall angle of airfoils

图15 不同粗糙度翼型升力系数、阻力系数对比Fig.15 Variation of CL and CD for airfoils at different roughnesses

在转捩状态下，翼型前缘粘贴粗糙带后，翼型的最大升力系数大幅提高，阻力系数增大，失速迎角显著增大.相比C2翼型的最大升力系数，R3翼型的最大升力系数增加27.48%；相比C2翼型的失速迎角，R3翼型的失速迎角增加5°.

两个翼型在失速迎角附近的流场流线图如图16和17所示，图中x、y为空间坐标.从图中看，转捩状态下，光滑NACA0012翼型在9° 迎角时前缘出现已分离泡，这是由于翼型头部前缘半径较小而且在翼型头部有很高的负压峰，紧接负压峰后是一个很强的逆压梯度，翼面上该处的层流边界层在逆压梯度的作用下出现层流分离，并且流动很快地转捩为湍流，之后流动再附着于壁面形成分离泡.随迎角增大，头部负压前移，分离点也逐渐前移，分离气泡亦随之扩大.当迎角达到11° 时，分离气泡因不能重新附着壁面而突然破裂，导致翼型上表面突然完全分离，翼型失速.失速时上翼面出现的大范围分离引起翼型升力急剧下降.此时翼型失速类型为前缘失速.

而对于粗糙翼型，在翼型前缘没有层流分离泡存在，翼型上表面的分离是从后缘开始发展的，随迎角增大而分离区逐渐向上游延伸、扩大，最终在背风区形成大范围分离.此时翼型失速类型为后缘失速.由于失速延迟，失速迎角增大，粗糙翼型的最大升力系数得到显著提高.因此在此流动条件下，前缘粗糙带使NACA0012翼型的失速类型由前缘失速变为后缘失速.

图16 C2翼型失速附近前缘分离泡及绕流流线图Fig.16 Streamlines and contours of streamwise velocity for C2 airfoil near stall

图17 R3翼型失速附近流线图Fig.17 Streamlines and contours of streamwise velocity for R3 airfoil near stall

为进一步分析转捩流动中粗糙前缘对翼型气动特性的影响规律，提取压力系数、表面摩擦因数分布及绕流边界层的参数进行对比.10° 迎角时，翼型的压力系数分布和表面摩擦因数分布如图18所示.

从压力系数分布对比图可以看到，在10° 迎角时，粗糙翼型与光滑翼型的压力分布差别较小.而在表面摩擦因数分布对比图中，两个翼型的表面摩擦因数在前缘处的分布有较大差别.图19所示为两个翼型前缘流线图.由图可知，光滑翼型前缘存在层流分离泡，使前缘处表面摩擦因数小于0，之后由于流动发生转捩，表面摩擦因数迅速增大；而粗糙翼型的前缘没有分离泡存在，其表面摩擦因数在流动转捩后快速增大，而且由于前缘粗糙带的作用，表面摩擦因数的上扬幅度也较大，并在粗糙带结束后迅速下降.

截取翼型流向不同弦长站位处法向流场变量分布进行对比，变量包括当地速度、湍动能和涡黏性，对比结果如图20～22所示.

从图20看出，相同站位处，相比光滑翼型，粗糙翼型的边界层厚度更高，边界层速度型更陡峭.光滑翼型在前缘产生层流分离泡，且由分离泡引起流动转捩；而粗糙翼型虽然没有分离泡产生，但是粗糙带促使流动转捩，转捩后翼型表面过大的表面摩擦力加速边界层内的动量交换，湍流脉动增强，边界层能量耗散增大.同时过多的能量消耗会使得粗糙翼型边界层抵抗下游逆压梯度的能力减弱，随着迎角增大，后缘将出现分离并逐渐扩大.

图18 α=10° 翼型压力系数分布对比和表面摩擦因数分布对比Fig.18 Pressure coefficients and surface friction coefficients of airfoils at α=10°

图19 α=10° 翼型前缘附近流线图Fig.19 Streamlines near leading edge for airfoils at α=10°

图20 α=10° 翼型不同站位速度型分布对比Fig.20 Comparisons of velocity profiles at several streamwise locations at α=10°

图21 α=10° 翼型不同站位湍动能分布对比Fig.21 Comparisons of turbulence kinetic energy at several streamwise locations at α=10°

图22 α=10° 翼型不同站位涡黏性分布对比Fig.22 Comparisons of eddy viscosity at several streamwise locations at α=10°

从湍动能和涡黏性分布对比图中也可以看到，由于粗糙翼型过大的表面摩擦力对边界层能量衰减、湍流脉动增加、动量交换增强的巨大作用，使得粗糙翼型的边界层相比光滑翼型的边界层湍动能增大、涡黏性增大.

虽然粗糙翼型边界层能量衰减更大，在逆压梯度中抵抗分离的能力减弱，但是由于其边界层中强烈的掺混作用，粗糙翼型前缘没有出现分离泡，流动分离首先从后缘出现，失速延迟，失速迎角较大.而光滑翼型由于前缘层流分离泡破裂导致失速，失速迎角较小.

4 结论

使用考虑粗糙度影响的湍流模型和转捩模型对NACA0012翼型进行数值模拟，分析粗糙前缘对翼型边界层发展及失速类型的影响，得到以下结论：

(1) 在全湍状态下，光滑翼型的失速类型为后缘失速，而附加前缘粗糙带后，粗糙带使翼型前缘表面摩擦力增大，边界层混掺和湍流输运能力增强，壁面边界对边界层流动能量的损耗增强，由于边界层动能的衰减使边界层抵抗逆压梯度的能力减弱，后缘分离出现的迎角减小.因此该流动状态下，粗糙前缘未改变翼型失速类型，相比光滑翼型，粗糙翼型的失速迎角减小、最大升力系数减小、阻力系数增大.

(2) 在转捩状态下，光滑翼型由于前缘分离泡破裂、分离区突然扩展至整个翼型上表面而失速，失速类型为前缘失速.而附加前缘粗糙带后，由于湍流脉动增强，流动掺混作用增强，在光滑翼型失速迎角时，粗糙翼型前缘没有分离泡存在，翼型失速延迟，失速迎角增大，同时在粗糙翼型后缘产生分离，且分离区随迎角增大而扩大，最后造成翼型上表面大范围分离，其失速类型为后缘失速.因此该流动状态下，粗糙前缘改变翼型失速类型，相比光滑翼型，粗糙翼型的失速迎角增大，最大升力系数增大，阻力系数增大.