基于响应面方法的风力机叶片多目标优化设计研究

邓 磊，乔志德，宋文萍，高永卫

（西北工业大学 翼型、叶栅空气动力学国家级重点实验室，陕西 西安 710072）

0 引 言

随着风电技术的迅猛发展，风电设备的设计能力和制造技术日趋成熟，产品进入商品化阶段，风电在与传统电能的竞争中，经济因素已从某种意义上上升到首要因素，因此必然的趋势就是单机容量的大型化。单机容量的增加则意味着叶片长度增加，由于叶片重量与风轮半径近似成三次方关系，随之而来的是叶片重量随风轮半径的急剧增加。因此，新型风力机专用翼型的设计以及有更高的风能捕捉能力的叶片设计成为热点。为适应风机发展的大型化趋势，研究者提出了高升力、高升阻比风力机翼型的设计思想，通过使用高升力翼型以减小叶片的弦长，并进而减小叶片的重量，提高叶片的风能捕捉能力［1－2］。

随着计算机技术的发展，优化设计技术逐步应用于风力机叶片的设计中，并且由于其灵活、实用等优点成为研究的热点［3－6］。优化方法是通过对几何外形的优化迭代，以达到预期的气动参数或者风力机性能参数；这种方法优点是设计周期短，提高了设计水平。目前叶片的优化设计方法中，设计变量主要是叶片外形，包括径向弦长、扭转角分布和厚度分布等；设计目标包括在设计状态运行时最大化功率系数、年发电量（AEP）或最小化单位发电成本等，也有考虑噪声等的优化问题。

近年来，全局优化方法－响应面方法［7］逐渐引起了研究者的兴趣。响应面方法是试验设计与数理统计相结合的优化方法：当实验结果（实验室实验结果或者数值计算结果）与已知参数间的函数关系式为隐式时，通过响应面方法可以设计实验结果与参数变量间的关系，并在实验测量值或数值分析基础上，对指定的试验点集合通过连续的试验，求得参数变量的系数，最终建立响应与参变量间的函数关系式。当建立了响应面模型，并形成适当的最优化问题后，剩下的是用数学规划来求出最优解。

本文使用西北工业大学翼型研究中心设计的NPU－WA族风力机专用翼型，基于响应面优化设计技术，对某1.5MW水平轴变速变距型风力机进行了多目标、多约束叶片优化设计研究。

1 基于BEM理论的叶片性能预测

近年来随着计算机技术的发展，基于Navier－Stokes（N－S）方程的CFD数值模拟逐渐用于风轮的气动性能预测中［8－10］，其优点是不需要提供径向使用翼型的可靠的气动性能数据，但是由于风轮运行状态非常复杂和计算花费因素，目前尚不能有效地用于实际叶片设计中。而基于动量－叶素理论的风轮气动性能预测方法由于简单快捷的特点，使用非常广泛；而且在提供可靠的翼型气动性能前提下，预测结果有相当的精度，但缺点是对二维翼型气动性能数据的依赖性高，并被认为性能预测的主要误差即来自于翼型的气动性能数据。在本文的设计研究中，需要进行大量的风轮性能的预测，本文选择基于BEM理论的PROPID软件来进行风轮的性能分析。PROPID是在性能预测程序PROP基础上发展而来的水平轴风力机叶片性能分析和反设计程序［11］，本文使用其PROPID51的版本，使用了Prantl叶尖损失、轮毂损失模型和Corrigan Post－stall模型，在常速和变速风轮的运行中，PROPID的性能预测结果均有很好的精度［12］。

2 响应面方法

2.1 响应面方法

响应面方法是采用试验设计理论，对给定的设计点集合进行试验，得到目标函数的响应面模型，来预测非设计点的响应值。由于二阶多项式模型形式比较灵活，对真实响应近似程度较好，且便于求解，因此很多工程问题采用二阶多项式模型。如果采用完全的二阶多项式作为响应面模型，构造模型所需的试验次数与设计变量的平方成正比，限制了在更多设计变量和多目标设计问题中的应用。为减小计算量，可以使用不含二阶交叉项的简化模型：

上式可以看出nrc＝2nv＋1，其中nv为设计变量数，nrc为模型包含的项数；ns为试验次数，一般取为nrc的1.5～3倍［13］；如进行26个设计变量的优化，构造一次完全二阶多项式模型需要试验700次左右，而构造简化模型只需试验80－120次，大大减小了试验次数。使用简化模型在减少了大量的二阶交叉项的同时，并没有改变模型的二阶多项式特性；但是对同一响应空间的拟合精度会有所降低，若是要求达到同样的精度，就需要缩小设计空间，以弥补舍掉交叉项带来的损失［14］。

2.2 设计变量

风力机叶片的优化设计中，设计变量一般选择为叶片径向站位上的弦长和扭转角的分布。本文中，叶片外形的变化扰动是通过基本形函数的线性组合来实现的：

其中，Δy表示弦长或者扭转角的变化量；δk取值（－1，0，1），控制局部扰动方向；γk为设计变量控制因子，其值的大小对优化过程和优化效率有直接关系，对优化结果也有影响；δkγk为控制变量；bk为形函数：

上述形函数在设计变量所在位置xk取得最大值，在远离xk时逐渐减小。由于在叶片优化中，设计变量所在的位置一般不会位于叶根或者叶尖位置，因此在叠加形状函数后，叶片在中段形状变化较大，在叶根和叶尖位置变化量减小。在本文计算中，径向使用11个截面位置，每个截面位置有弦长和扭转角两个设计变量，共有22个设计变量。试验设计满足D－优化准则；响应面模型使用不含二阶交叉项的二阶多项式模型，构造一次响应模型进行80次试验。

2.3 目标函数

本文基于响应面方法进行叶片的多目标、多约束的优化设计研究。

目标函数：

约束条件：

A表示叶片的面积，下标0表示原始叶片的值；AEP表示年发电量，Cp，max表示功率系数最大值。

对多目标优化时，本文通过“统一目标函数”的方法，将多个目标函数统一到一个目标函数中：

其中，f为统一目标函数，q为分目标的个数，ωj为根据重要性给出的第j个目标函数的权值。frj表示通过转化设计目标法得到的在（0～1）之间的无量纲的目标函数。

构造完响应面模型并形成适当的最优化问题之后，就需要用数学规划来求出优化解。本文使用BOX复合形搜索法［15］来进行搜索，这种方法对带约束的优化问题有很好的鲁棒性。

3 某1.5MW风力机叶片优化设计

3.1 问题描述

本文以某1.5MW风力发电机叶片设计为例，进行叶片优化设计研究。径向使用西北工业大学翼型研究中心设计的NPU－WA－族风力机专用翼型，风轮主要设计参数如表1所示。初始叶片长度38m，最大扭转角15°，径向使用翼型 NPU－WA－400／350／300／250／210／180／150，其中数字前两位表示最大相对厚度，最后一位表示设计次数。如400表示40%相对厚度翼型初始设计。布置位置如表3左侧所示（翼型所在位置为使用此翼型的中心位置），厚度的分布近似于文献中的较厚的一类分布［16］，这类分布兼顾了翼型的气动性能，也利于减轻叶片重量。弦长分布为某叶片生产商给出的初始分布。

在PROPID的输入文件中，使用默认的沿径向均匀布置的弦长和扭转角分布，因此需要将表2中左侧径向弦长和厚度分布插值到这些均匀分布的位置上，如表2右侧所示。

表1 某1.5MW风力机设计参数Table 1 Design parameters for 1.5MW HAWT

表2 某1.5MW风力机叶片外形Table 2 Primary geometric data for 1.5MW HAWT blade

3.2 使用翼型CFD计算气动数据的设计

本节使用RANS方程和转捩位置计算耦合求解的方法计算径向翼型的气动性能用于叶片设计，其中转捩位置使用基于线性稳定性分析的eN方法来计算。在叶片70%以内翼型有15mm的后缘厚度，约0.5%～0.6%弦长；在叶片70%径向以外翼型有3mm的后缘厚度，约0.2%～0.3%弦长；由于后缘厚度小，因此在二维翼型性能计算中，使用尖后缘计算。

为研究在统一目标函数中分目标的权值分配对设计结果的影响，本文使用三个权值分配组合，其中定义为Opti.0.2－0.8的叶片表示年发电量目标的权值是0.2，最大功率系数权值是0.8，以此类推。年发电量的计算时，假设风速满足威布尔分布，尺度参数A＝0.892244，形状参数K＝1.7。

由表3可以看出，权值的分配对设计结果有比较大的影响，其中随着最大化功率系数目标权值的增加，叶片面积的减小非常明显。图1比较了设计结果外形，可以看出在功率系数的权值较大时弦长的减小非常明显，同时优化增加了叶片内侧的弦长和扭转角。需要指出的是原始叶片的外侧的扭转角在减小到0°之后又增加，这是工程的需要而不是气动的需要；本算例进行气动设计时，没有考虑到这个变化。

表3 用CFD计算数据的设计结果Table 3 Design results using CFD calculated data

图2为设计结果的主要气动性能比较。由功率系数曲线可以看出，优化后增加了最大功率系数，但是由于叶片弦长的减小，造成了低尖速比时功率系数变差；由功率曲线可以看出，优化后额定风速也有所降低；由径向升力系数的分布可以看出，优化后径向30%向外的升力有明显提高，在径向50%以外的升力系数在1.1～1.2之间，这有利于叶片弦长的减小。由本算例设计结果可以看出，优化后提高了叶片在设计点的气动性能，但是同时可以看出在叶片的多目标设计中，设计结果的性能随目标权值的变化。

3.3 使用翼型风洞测量气动数据的设计

2010年NPU－WA－族风力机翼型在西北工业大学NF－3风洞进行了试验，二元实验段尺寸为1.6m（宽）×3.0m（高），模型尺寸1.595m（展长）×0.8m（弦长），试验雷诺数分别为1.0×106、1.5×106、3.0×106、4.0×106和5.0×106，并进行了自由转捩和固定转捩的试验。本节使用自由转捩试验测量数据进行了优化设计，并使用固定转捩数据分析了设计结果以考察前缘污染对设计结果气动性能的损失的影响。在本算例中，两个设计目标的权值分配为0.5－0.5。

图3比较了优化前后的叶片的外形，可以看出优化叶片减小了弦长分布，叶片面积减小了11.31%，这对叶片重量的减轻有显著的贡献；扭转角的变化也和CFD数据设计结果类似。

图3 使用风洞测量数据设计结果外形比较Fig.3 Geometric properties of design results using measured data

图4（a）比较了功率系数曲线，可以看出优化叶片增加了设计点的最大功率系数，但是在低尖速比，特别是尖速比小于6时功率系数比原始叶片有比较明显的减小，这是叶片弦长减小的结果。同时从图中前缘污染的功率系数可以看出，优化叶片的最大功率系数虽然比初始叶片略大，但是整体性能差于原始叶片，说明只使用最大功率系数作为设计目标还需要进一步的研究。图4（b）比较了功率曲线，可以看出优化叶片在干净前缘时额定风速比原始叶片小，但是在前缘污染时，额定风速有所增加，但是仍在设计要求范围内。图4（c）比较了叶片径向升力系数，可以看出优化叶片增加了径向30%以外的升力系数，且升力系数近似于均匀地分布在1.2左右，这个升力系数对应于5°～6°攻角，是翼型最大升阻比对应的攻角范围。在前缘污染时，径向20%左右升力系数有非常明显的减小，这是由于这个位置上布置是的NPUWA－400翼型，厚翼型在前缘污染时气动性能损失很大。图4（d）比较了优化前后径向雷诺数的比较，可以看出优化叶片减小了雷诺数分布，这是叶片弦长减小的结果；30%径向以外雷诺数在300－450万之间，这个雷诺数还没有达到NPU－WA翼型族的设计雷诺数，翼型族还可以用于更大尺寸的叶片设计中。

图4 优化前后气动性能比较（由风洞测量数据）Fig.4 Comparisons of aerodynamic performance（from measured data）

表4 表面条件对气动性能的影响Table 4 Aerodynamic performance with clean and soiled surface conditions

表4比较了前缘条件对叶片气动性能的影响。可以看出优化叶片的最大功率系数均比原始叶片要大，但是在前缘污染时年发电量的减小却高于原始叶片，这是由于在前缘污染时功率系数的显著减小造成的，因此在此后的研究中，需要对此进行进一步的研究。

4 结 论

通过本文的研究设计，可以得出如下结论：

（1）响应面方法可以成功地应用于风力机叶片的多目标、多约束的优化设计中，但是目标函数的权值分配对设计结果有较大的影响；

（2）最大功率系数的权值增加引起径向弦长分布的减小，并进而减小叶片的面积，将有效地减小叶片的重量，但是将引起低尖速比时叶片性能降低；

（3）使用NPU－WA族风力机翼型进行叶片设计，增加了叶片径向升力系数分布。径向升力系数近似均匀分布，攻角在翼型的最大升阻比对应攻角附近，这和NPU－WA族翼型设计思想吻合；

（4）由于本文进行的是单设计点的设计，因此设计结果在非设计点的性能有降低，特别表现在前缘污染时气动性能下降，在此后的研究中需要考虑多设计点的优化问题。

［1］DAN M.Somer，the S830，S831，and S832airfoils［R］.NREL／SR－500－36339，2005.

［2］DAN M.Somer，the S833，S834，and S835airfoils［R］.NREL／SR－500－36340，2005.

［3］LANZFAME R，MESSINA M.Optimal wind turbine design to maximize energy production［J］.Journal of Power and Energy，2009，223（2）：93－101.

［4］WANG X D，WEN Z S，WEI J Z.Shape optimization of wind turbine blades［J］.Wind Energy，2009，12（8）：781－803.

［5］MENDEZ J，GREINER D.Wind blade chord and twist angle optimization using genetic algorithms［A］.The 5thInternational Conference on Engineering Computational Technology（ECT 2006）［C］.Las Palmas de Gran Cnaria，Spain，2006.

［6］KUSIAK A，ZHENG A，SONG Z.Power optimization of wind turbines with data mining and evolutionary computation［J］.Renewable Energy，2010，35：695－702.

［7］MYERS R H，MONTGOMERY D C.Response surface methodology：process and product optimization using designed experiments［M］.NY：Wiley ＆Sons，2002：1－78，351－401.

［8］GUPTA S G，LEISHMAN J.Performance predictions of the NREL phase VI combined experiment rotor using a free－vortex wake model［R］.AIAA 2006－390.

［9］PARK Y M，CHANG B H.Numerical simulation of wind turbine scale effects by using CFD［R］.AIAA 2007－216.

［10］许建华，宋文萍，韩忠华.基于嵌套网格的风力机叶片扰流数值模拟［J］.太阳能学报，2010，31（1）：91－95.

［11］DAYTON A GRIFFIN.NREL Advanced Research Turbine （ART）aerodynamic design of ARRT－2B rotor blades［R］.NREL／SR－500－28473，2000.

［12］TANGLER J，DAVID KOCUREK J.Wind turbine post－stall airfoil performance characteristics guidelines for blade－element momentum methods［R］.NREL／CP－500－36900，2005.

［13］AHN J，KIM H J，LEE D H，et al.Response surface method for airfoil design in transonic flow［J］.Journal of Aircraft，2001，38（2）：231－238.

［14］熊俊涛，基于响应面方法的气动优化设计［D］.［硕士学位论文］.西北工业大学，2005.

［15］卢险峰.最优化方法应用基础［M］.上海：同济大学出版社，2003：107－137.

［16］TPI Composites.Innovative design approaches for large wind turbine blades－final report［R］.SAND2004－0074，2004.