0 引言
风力机多工作在高寒、沿海、风沙频繁等环境恶劣的地区，其叶片表面经常会附着灰尘、昆虫、雨雪等污垢。污垢增大了叶片表面的粗糙度，并导致风力机的年发电量大大降低[1-3]。对粗糙度敏感性低的翼型的设计，是减轻或消除附着污垢对翼型气动性能影响的有效解决方案，进而可保证风力机在叶片表面粗糙条件下的高效运行。
翼型直接优化设计能很好地解决难以给定恰当目标压力和速率分布的问题，但需要对翼型进行型线参数化表达[4]。文献[5-9]采用Bézier函数、Hicks-Henne函数构造翼型参数化表达式，利用遗传算法并结合XFOIL软件以及人工神经网络模型等进行翼型的优化设计。文献[10]探讨了各种翼型的前缘粗糙度敏感特性，并基于翼型广义泛函集成表征形式优化出一种新翼型。研究发现，广义泛函集成表达方法相比于其他翼型型线表征形式，更易于优化和扩展形成新的翼型型线，但不能很好控制中等和大厚度翼型尾缘处的型线[11-12],因此，文献[12]提出翼型广义泛函集成表达与B样条曲线相结合的方法，以此进行风力机翼型型线的优化设计。
钝尾缘改型设计可以有效改善表面粗糙翼型的气动性能。BAKER等[13]实验研究了对称加厚的不同翼型，发现适度增加尾缘厚度可增大升阻比并降低前缘粗糙敏感度。杨瑞等[14]采用CFD方法模拟薄、钝尾缘翼型的气动性能，结果表明钝尾缘翼型增大了最大升力系数，并减小了前缘污染对升力特性的影响。
鉴于优化设计和钝尾缘改型均能提高表面粗糙翼型的气动性能，本文结合广义泛函集成表达和B样条曲线来形成钝尾缘翼型型线参数化表达式，利用粒子群算法耦合XFOIL软件进行翼型的优化设计，平移优化后钝尾缘翼型吸力面距前缘0.1c(c为弦长)处添加一凸台，来研究表面粗糙翼型钝尾缘优化前后的气动性能和流场特性。
1 钝尾缘翼型型线的表达
采用广义泛函集成表达方法描述上翼面距前缘0.4c和下翼面距前缘0.5c之前的型线，则翼型的坐标分量为[12]
(1)
(2)
式中：a为1/4弦长；θ为幅角；ρ(θ)为翼型形状函数；Ci为多项式系数，C0=1。
为保证翼型型线广义泛函集成表达与B样条曲线在结合点处的连续、光滑，采用3次B样条曲线表示上翼面距前缘0.4c和下翼面距前缘0.5c之后的型线，则型线的坐标分别为
(3)
t∈[0, 1]
式中，P0～P3为上翼面距前缘0.4c的型线的控制变量;
为下翼面距前缘0.5c之后型线的控制变量。
上翼面0.4c、下翼面0.5c之后型线
的点P0,3(0)和
分别经过型线广义泛函集成表达形成的上下翼面型线的最后一个坐标点，P0,3(1)和
经过上下翼面型线终点(1, hj)和(1, h(1-j)) (h为钝尾缘厚度，j为上翼面与钝尾缘的厚度比)。利用式(3)可反算出P0、P3和
则上翼面0.4c和下翼面0.5c之后型线的控制变量只有P1、P2和
型线的坐标具体表示为
(4)
2 表面粗糙钝尾缘翼型的优化设计
基于上述钝尾缘翼型型线参数化表达方法，利用MATLAB编写粒子群优化算法程序，结合XFOIL软件计算翼型气动性能，进行钝尾缘翼型优化设计。平移优化后翼型指定位置坐标添加一凸台来获得表面粗糙的钝尾缘翼型。
2.1 设计变量和约束条件
为较好控制钝尾缘翼型的型线，选取形状函数控制方程的前11项系数C1～C11、B样条曲线控制参数以及尾缘厚度和其分配比为优化设计变量X：
(5)
在钝尾缘翼型优化过程中，需限定控制变量的范围，以避免形成的型线不具有翼型的形状特征，则优化变量边界约束条件为
xmin≤x≤xmax
(6)
式中，x为设计变量X中的任一分量；xmin、xmax分别为x的下限和上限。
大部分风力机叶片在其主要功率产生区采用最大厚度为0.12～0.25的翼型，又由于翼型的相对厚度对叶片结构特性的影响不容忽视，故翼型厚度tr的约束条件为
0.12≤tr≤0.25
(7)
2.2 目标函数
升阻力系数、阻力系数是研究风力机翼型气动性能的两个重要指标。为使翼型在风力机正常运行工况下具有良好的气动性能，选择翼型的最大升阻比为目标函数：
f(X)=max(CL/CD)
(8)
式中，CL、CD分别为翼型的升力系数和阻力系数。
2.3 粗糙翼型型线优化设计
粒子群算法具有实现容易、精度高、收敛快等优点，解决实际问题的效率很高。XFOIL软件基于源面方法和边界层理论计算翼型的气动性能，计算的过程较简单且精度较高[15]。因此，基于建立的钝尾缘翼型型线优化模型，采用粒子群算法耦合XFOIL软件，并通过MATLAB编写优化程序，对S812翼型进行型线优化设计。S812翼型用于风力机叶片的主要功率产生区，最大相对厚度为21%并位于0.39c处，且最大相对弯度为1.57%。
取雷诺数Re=1×106，马赫数Ma=0.04，种群规模为20，最大进化代数为300，学习因子S1、S2为0.5，变量维数为21。为使算法程序既有较强的搜索能力，又有较好的收敛性，惯性权重w=wmax-m(wmax-wmin)/mmax[16]，其中，wmax、wmin分别为惯性权重的最大值和最小值，依据经验取1.2和0.4；m和mmax分别为当前的进化代数和最大的进化代数。
本文平移优化后翼型指定位置的坐标，在翼型吸力面距前缘0.1c处添加一高0.015c、宽0.04c的凸台，以模拟翼型表面粗糙情况，从而得到表面粗糙S812翼型(取名为S812-R翼型，R表示粗糙)及其钝尾缘改型S812-R-BT(BT表示钝尾缘)的形状，如图1所示。S812-R-BT翼型的尾缘厚度为0.039 8c，上下翼面尾缘厚度比为1∶13.16。
图1 优化前后表面粗糙翼型型线
Fig.1 Profiles of rough airfoil before and after optimization
3 优化前后的气动性能和流场特性
利用基于CFD方法的FLUENT软件，计算S812翼型的升力系数、阻力系数，并比较数值计算结果与实验数据，研究表面粗糙翼型钝尾缘优化前后的升力系数、阻力系数、升阻比、压力系数和流场特性，分析钝尾缘优化设计对粗糙翼型气动性能的提升效果。
3.1 数值计算方法
直径50c的半圆形和长50c、宽25c的矩形构成计算域，翼型位于半圆中心。计算域左侧以及上下两侧采用速度进流边界条件，速度由Re=1×106确定；计算域右侧采用压力出流边界条件，表压力给定0；翼型表面满足固壁绝热无滑移条件。利用C型结构化网格划分，并对翼型前后缘以及翼型附近的网格进行局部加密。控制方程组采用连续性方程和二维不可压缩Navier-Stokes方程，湍流模型选用较适合模拟稳态绕流情况的k-ω SST模型[17]。各方程离散采用二阶迎风格式，压力和速度耦合采用SIMPLE算法。速度和连续性方程的残差值分别为10-6和10-4，j和ω的残差值均为10-5。
S812-R翼型的计算域、计算方法与S812翼型相同，但S812-R翼型存在凸台，在网格划分上采用分块结构化网格，即从凸台侧边引出线段分割半圆形计算域，并用映射网格进行划分，得到结构化C型网格。S812翼型的升力系数、阻力系数数值结果与文献[18]提供的实验数据如图2所示，可以看出，升力系数在攻角小于8.2°时与实验值非常接近，在攻角大于8.2°之后大于实验值；阻力系数与实验值相差不大；升力系数、阻力系数的计算值和实验值的总体变化趋势基本一致。因此，可以采用上述数值计算方法分析表面粗糙翼型钝尾缘优化前后的气动性能和流场特性。
图2 S812翼型计算与实验的结果
Fig.2 Calculated and experimental data for S812 airfoil
3.2 表面粗糙翼型钝尾缘优化前后的气动性能
钝尾缘优化前后，表面粗糙翼型的升力系数、阻力系数和升阻比随攻角变化的曲线见图3。由图3a、图3b可以看出，S812-R-BT翼型的升力系数相较于S812R翼型明显增大；S812-R-BT翼型的阻力系数在17.2°攻角之前与S812-R翼型的阻力系数非常接近，17.2°攻角之后大于S812-R翼型的阻力系数。由图3c可以看出，S812-R-BT翼型的升阻比在攻角小于17.2°时明显大于S812R翼型的阻力比，17.2°攻角之后与S812-R翼型十分接近；S812-R-BT翼型的最大升阻比较S812-R翼型的最大升阻比有大幅提高。
图3 S812-R和S812-R-BT翼型的升力系数、阻力系数和升阻比
Fig.3 Lift, drag coefficients and lift-drag ratios of S812-R and S812-R-BT airfoils
3.3 表面粗糙翼型钝尾缘优化前后的压力系数
S812-R翼型及其优化后钝尾缘翼型S812-R-BT的表面压力分布如图4所示。由图4可以看出，17.2°攻角之前，S812-R-BT翼型吸力面的压力系数减小，压力面的压力系数增大；18.3°攻角时，吸力面的压力系数稍微增大，压力面的压力系数也是如此。所以，表面粗糙翼型钝尾缘优化后，上下翼面的压差增加，即S812-R-BT翼型的升力系数高于S812-R翼型，这与图3的规律一致。
图4 S812-R和S812-R-BT翼型的表面压力分布
Fig.4 Pressure distributions around S812-R and S812-R-BT airfoils
3.4 表面粗糙翼型钝尾缘优化前后的流场特性
由图5可以看出，2.1°攻角时，S812-R翼型上翼面前缘凸台后面形成较小的漩涡，尖尾缘处未发生流动分离；S812-R-BT翼型凸台后面没有漩涡，钝尾缘处却存在一对方向相反交替脱落的漩涡，这是因为钝尾缘对气流形成阻挡作用。同时，钝尾缘处的漩涡使吸力面的气流发生下洗作用，因而S812-R-BT翼型压力面的压力增大，吸力面的压力减小，进而引起升力系数增大。
图5 S812-R和S812-R-BT翼型的流线图(α=2.1°)
Fig.5 Streamlines around S812-R and S812-R-BT airfoils(α=2.1°)
由图6可以看出，6.3°攻角时，S812-R翼型凸台后面仍存在小的漩涡， 吸力面尾缘处的漩涡交替脱落，且已经发生流动分离，这对应于图3中升力系数曲线斜率稍有下降的现象；而S812-R-BT翼型依然仅在钝尾缘后面存在一对交替脱落的漩涡，且其尺寸大于2.1°攻角时的漩涡。
图6 S812-R和S812-R-BT翼型的流线图(α=6.3°)
Fig.6 Streamlines around S812-R and S812-R-BT airfoils(α=6.3°)
由图7～图9可以看出，当攻角增大到17.2°时，S812-R-BT翼型凸台后面才产生漩涡。
由图5～图9还可以看出，随着攻角的增加，S812-R翼型及其钝尾缘改型S812-R-BT的凸台后面漩涡逐渐增长，尾缘处的漩涡也持续增长并向吸力面前缘扩张，且S812-R翼型尾缘处的漩涡更靠近前缘；S812-R-BT翼型钝尾缘后面的漩涡使吸力面气流的下洗作用持续发生，因而S812-R-BT翼型的升力系数一直大于S812-R翼型的升力系数。
图7 S812-R和S812-R-BT翼型的流线图(α=13.1°)
Fig.7 Streamlines around S812-R and S812-R-BT airfoils(α=13.1°)
图8 S812-R和S812-R-BT翼型的流线图(α=17.2°)
Fig.8 Streamlines around S812-R and S812-R-BT airfoils(α=17.2°)
图9 S812-R和S812-R-BT翼型的流线图(α=18.3°)
Fig.9 Streamlines around S812-R and S812-R-BT airfoils (α=18.3°)
4 结论
(1)以翼型型线控制方程组的形状函数系数、B样条曲线控制点以及尾缘厚度及其在上翼面分配比为设计变量，利用粒子群算法并结合XFOIL软件进行了S812翼型钝尾缘优化；平移优化后，在翼型吸力面0.1c处坐标添加了一高0.015c、宽0.04c的凸台，得到尾缘厚度为0.039 8c且上下翼面尾缘厚度比为1∶13.16的表面粗糙S812翼型钝尾缘改型。
(2)表面粗糙翼型钝尾缘优化后，升力系数大幅增大；阻力系数在17.2°攻角之前与原型相差不多，在17.2°攻角后阻力系数明显增大；升阻比在17.2°攻角之前显著增大，17.2°攻角后与原型基本一致，最大升阻比明显增大。
(3)钝尾缘翼型尾缘处的漩涡对吸力面的气流造成下洗作用，使得吸力面和压力面的压差变大，推迟了翼型的失速,因而，表面粗糙翼型钝尾缘优化后，气动性能得到明显提高。
[1] SAGOL E, REGGIO M, ILINCA A. Issues Concerning Roughness on Wind Turbine Blades[J]. Renewable and Sustainable Energy Reviews, 2013, 23: 514-525.
[2] CORTEN G P, VELDKAMP H F. Insects Can Halve Wind-turbine Power [J]. Nature, 2001, 412(6842): 41-42.
[3] 焦灵燕, 汪建文, 贺玲丽. 粗糙度对风力机翼型气动性能影响的模拟研究[J]. 可再生能源, 2014, 32(12): 1816-1820.
JIAO Lingyan, WANG Jianwen, HE Lingli. Simulation Study on Effect of Surface Roughness on Aerodynamic Performance of Wind Turbine Airfoil [J]. Renewable Energy Resources, 2014, 32(12): 1816-1820.
[4] 陈婉春, 孙刚. 基于气动特性翼型参数化方法的适用性研究[J]. 力学季刊, 2015, 36(4): 678-689.
CHEN Wanchun, SUN Gang. Applicability of Airfoil Parameterization Methods Based on Aerodynamic Performance [J]. Chinese Quarterly of Mechanics, 2015, 36(4): 678-689.
[5] RIBEIRO A F P, AWRUCH A M, GOMES H M. An Airfoil Optimization Technique for Wind Turbines[J]. Applied Mathematical Modelling, 2012, 36(10): 4898-4907.
[6] DJAVARESHKIAN M H, LATIFI A. Optimization of Wind Turbine Airfoil with Good Stall Characteristics by Genetic Algorithm Using CFD and Neural Network [C]//ASME 2013 International Mechanical Engineering Congress and Exposition. San Diego, 2013: 15-21.
[7] LIU Y X, YANG C, SONG X C. AnAirfoil Parameterization Method for the Representation and Optimization of Wind Turbine Special Airfoil [J]. Journal of Thermal Science, 2015, 24(2): 99-108.
[8] SANAYE S, HASSANZADEH A. Multi-objective Optimization of Airfoil Shape for Efficiency Improvement and Noise Reduction in Small Wind Turbines [J]. Journal of Renewable and Sustainable Energy, 2014, 6(5): 5-15.
[9] 琚亚平, 张楚华. 基于人工神经网络与遗传算法的风力机翼型优化设计方法[J]. 中国电机工程学报, 2009, 29(20): 106-111.
JU Yaping, ZHANG Chuhua. Optimal Design Method for Wind Turbine Airfoil Based on Artificial Neural Network Model and Genetic Algorithm [J]. Proceedings of the CSEE, 2009, 29(20): 106-111.
[10] 陈进, 张石强, 王旭东, 等. 基于粗糙度敏感性研究的风力机专用翼型设计[J]. 空气动力学学报, 2011, 29(2): 142-149.
CHEN Jin, ZHANG Shiqiang, WANG Xudong, et al. Dedicated Wind Turbine Airfoil Design Based on the Roughness Sensitivity Considerations[J]. Acta Aerodynamica Sinica, 2011, 29(2): 142-149.
[11] 王文剑, 黄俊, 宋磊. 基于实用参数控制的翼型参数化方法[J]. 飞机设计, 2013, 33(3): 1-4.
WANG Wenjian, HUANG Jun, SONG Lei. A Geometric Representation of Airfoil Based on Applied Parameters [J]. Aircraft Design, 2013, 33(3): 1-4.
[12] 陈进, 汪泉, 李松林, 等. 翼型集成理论与B样条结合的风力机翼型优化设计方法研究[J]. 太阳能学报, 2014, 35(10): 1930-1935.
CHEN Jin, WANG Quan, LI Songlin, et al. Study of Optimization Design Method for Wind Turbine Airfoil Combining Airfoil Integrated Theory and Bspine [J]. Acta Energiae Solaris Sinica, 2014, 35(10): 1930-1935.
[13] BAKER J P, MAYDA E A, van DAM C P. Experimental Analysis of Thick Blunt Trailing-edge Wind Turbine Airfoils [J]. Journal of Solar Energy Engineering, 2006, 128(4): 422-431.
[14] 杨瑞, 李仁年, 张士昂, 等. 钝尾缘风力机翼型气动性能计算分析[J]. 机械工程学报, 2010, 46(2): 106-110.
YANG Rui, LI Rennian, ZHANG Shiang, et al. Computational Analyses on Aerodynamic Characteristics of Flatback Wind Turbine Airfoils [J]. Journal of Mechanical Engineering, 2010, 46(2): 106-110.
[15] 张旭, 李伟. 非对称钝尾缘厚度对风力机翼型气动性能的影响[J]. 中国电机工程学报, 2013, 33(14): 107-113.
ZHANG Xu, LI Wei. Effect of Asymmetric Blunt Trailing-edge Thickness on Aerodynamic Performance of Wind Turbine Airfoils [J]. Proceedings of the CSEE, 2013, 33(14): 107-113.
[16] SHI Y H, EBERHART R. A Modified Particle Swarm Optimizer [C]// IEEE Conference on Evolutionary Computation. Anchorage, 1998: 69-73.
[17] 任年鑫, 欧进萍. 大型风力机二维翼型气动性能数值模拟[J]. 太阳能学报, 2009, 30(8): 1087-1091.
REN Nianxin, OU Jinping. Numerical Simulation for Pneumatic Characteristics for Two-dimensional Airfoils Large Wind Turbine [J]. Acta Energiae Solaris Sinica, 2009, 30(8): 1087-1091.
[18] RAMSAY R R, GREGOREK G M.Effects of Grit Roughness and Pitch Oscillations on the S812 Airfoil [R]. Golden: National Renewable Energy Laboratory, 1995: 1-149.
