近年来，随着海上风力机等风电设备的持续发展^[1]，风力机越来越趋向于大型化，而额定功率大的风力机使用的叶片也更大、更长。叶片的气动外形直接决定着风力机的气动效率及叶片所承受的载荷情况。更大、更长的叶片在为风力机带来高功率的同时也使叶片本身所受的气动载荷提高。过大的载荷可能使叶片更易发生形变，从而影响风力机的气动性能与运行稳定性^[2]。因此，在对大型风电叶片进行优化时，不仅要考虑气动效率的提升，还要尽可能降低叶片所受的气动载荷。

国内外学者在对风力机气动外形进行优化设计时，叶片设计变量一般选取组成叶片的各截面翼型形线，或叶片的弦长与扭角沿展向的分布情况。王威等^[3]基于熵产理论，对翼型进行了考虑最大升阻比与最小熵产率的多目标优化。司景喆等^[4]基于粒子群算法，针对特定的风场条件对叶片的弦长与安装角分布进行了气动优化。Wen等^[5]采用GABP（genetic algorithm back propagation）人工神经网络对S809翼型进行了改良，得到了具有更大升力系数与最大升阻比的翼型。汪泉等^[6]以风轮功率系数与噪声声压级为优化目标，得到了使叶片气动性能更好、噪声更低的气动外形。陈刚等^[7]以气动性能对攻角的敏感程度为改进方向，得到了性能更好的风力机翼型。

本文将二维翼型与三维叶片的几何优化相结合，综合考虑各截面翼型形线与叶片扭角沿展向的分布，同时对高于额定风速的叶片桨距角进行改良。基于翼型参数化方法、CFD（computational fluid dynamics）计算、动量叶素理论（blade element momentum，BEM）及多岛遗传算法搭建自动优化平台，以美国NREL实验室的5 MW大型风电叶片为例，完成以最大气动效率与最小气动载荷为目标的多目标优化设计。

1 优化方法及验证 1.1 气动性能计算方法

在对叶片进行优化的过程中，需反复计算各截面翼型形线 、扭角产生变化后新叶片的气动效率与载荷，计算量较大。考虑计算量与计算准确性，分别求得各叶素处翼型气动参数，再基于动量叶素理论进行叠加，从而求得叶片整体气动性能。

Xfoil气动分析程序采用无黏的线性涡面元法，通过使用以层流线性稳定理论为基础的e^N函数控制法来准确识别流动转捩。由Xfoil气动分析计算所得的风力机翼型升阻力数据与实验结果较为符合^[8]，具有一定的可靠性，因此，本文以其作为气动参数计算的方法。

得到各截面翼型气动参数后，由动量叶素理论^[9]将叶片分为多个叶素，能够快速、准确地得到叶片的气动性能^[10]。由动量定理，定义轴向、切向诱导速度V_a，V_b与轴、切向诱导因子a，b，并以此得到叶素位置处的相对速度V_rel、入流角$\phi$与当地攻角α，其相关表达式为

$ \left\{ \begin{gathered} {V_a} = {v_0}(1 - a) \\ {V_b} = wr(1 + b) \\ \end{gathered} \right. $

(1)

$ {V_{{\text{rel}}}} = \sqrt {{v_0^2}{{(1 - a)}^2} + {w^2}{r^2}{{(1 + b)}^2}} $

(2)

$ \left\{ \begin{array}{*{20}{l}} \alpha = \phi - \theta \\ \phi = \arctan \left[ {\dfrac{{{v_0}(1 - a)}}{{wr(1 + b)}}} \right] \\ \end{array} \right. $

(3)

式中：v₀表示来流速度，m/s；w表示旋转角速度；r表示当地叶素位置处的半径；θ表示当地叶素位置处的安装角。

图1为翼型截面的速度三角形和受力情况。

图中，L，D分别表示升力、阻力。

基于动量理论和叶素理论，分别得到在叶素dr上的推力dT与扭矩dM的表达式，并引入叶片损失系数进行修正，即

$ \left\{ {\begin{array}{*{20}{l}} {{{\rm{d}}} T = 4{\text{π}}r\rho {v_0}a(1 - a){{\rm{d}}} r = \dfrac{1}{2}Brc{v_0}^2{C_{\text{n}}}{{\rm{d}}} r} \\ {{{\rm{d}}} M = 4{\text{π}}{r^3}\rho {v_0}w(1 - a)b{{\rm{d}}} r = \dfrac{1}{2}Brc{v_0}^2{C_{\text{t}}}r{{\rm{d}}} r} \end{array}} \right. $

(4)

$ \left\{ {\begin{array}{*{20}{l}} {{C_{\text{n}}} = {C_{\text{l}}}\cos\,\phi + {C_{\text{d}}}\sin \,\phi } \\ {{C_{\text{t}}} = {C_{\text{l}}}\sin \,\phi - {C_{\text{d}}}\cos\, \phi } \end{array}} \right. $

(5)

$ F = {F_{\text{t}}}{F_{\text{h}}} $

(6)

$ \left\{ {\begin{array}{*{20}{l}} {{F_{\text{t}}} = \dfrac{2}{{\text{π}}}\arccos ({{\text{e}}^{ - \frac{B}{2}\frac{{R - r}}{{r\sin \,\phi }}}})} \\ {{F_{\text{h}}} = \dfrac{2}{{\text{π}}}\arccos ({{\text{e}}^{ - \frac{B}{2}\frac{{r - {R_{{\text{hub}}}}}}{{{R_{{\text{hub}}}}\sin\, \phi }}}})} \end{array}} \right. $

(7)

式中：ρ表示来流密度；B表示风轮叶片数；C_n，C_t表示法、切向力系数；C_l，C_d表示升、阻力系数；F表示损失系数；F_t，F_h表示叶尖、轮毂损失；c表示弦长；R为叶轮半径；R_hub表示轮毂半径。

联立式（4）~（7），可得轴、切向诱导因子a，b为

$ \left\{ \begin{gathered} a = \dfrac{1}{{\dfrac{{4F{{\sin }^2}\phi }}{{\sigma {C_{\text{n}}}}} + 1}} \\ b = \dfrac{1}{{\dfrac{{4F\sin\, \phi \cos \,\phi }}{{\sigma {C_{\text{t}}}}} + 1}} \\ \end{gathered} \right. $

(8)

$ \sigma = \frac{{Bc}}{{2{\text{π}}r}} $

(9)

式中，σ表示弦长位置处的实度，即相应半径位置处叶片弦长之和与叶轮周长之比。

根据式（8）进行迭代计算，并将此计算过程用于所有叶素，并沿叶片展向积分，求得叶片的气动性能值。风力机的风能利用率C_p的表达式为

$ \left\{ \begin{array}{*{20}{l}} M = \displaystyle\int_0^R {4{\text{π}}\rho \lambda {v_0}^2{r^3}b} (1 - a){{\rm{d}}} r \\ {C_{\text{p}}} = \dfrac{P}{{\dfrac{1}{2}\rho {v_0}^3S}} = \dfrac{{Mw}}{{\dfrac{1}{2}\rho {v_0}^3{\text{π}}{R^2}}} \end{array} \right. $

(10)

式中：M表示风轮扭矩；λ表示尖速比；P表示输出功率；S表示扫风面积。

叶根弯矩M_t的表达式为

$ \left\{ \begin{array}{*{20}{l}} {M_{{i}}} = \displaystyle\int_{{r_i}}^{{r_{i + 1}}} {\dfrac{1}{2}\rho {v_0}^2c{C_{\rm{n}}}(r - {r_0})} {{\rm{d}}} r \\ {M_{\text{t}}} = \displaystyle\sum\limits_{i = 1}^n {{M_{{i}}}} \end{array} \right. $

(11)

式中，M_i表示叶素处所受挥舞弯矩。

1.2 计算方法验证

本文以美国NREL实验室的5 MW大型风电叶片^[11]为优化目标，其相关参数见表1。

沿叶片展向选取7个截面、6种翼型参与优化，其几何参数与翼型种类见表2。

本文的气动性能计算方法为基于前述理论的自编程序。为验证计算方法的可靠性，将不同风速下的计算结果与设计值及CFD方法计算结果进行对比。CFD方法通过建模与网格划分，基于k-ω SST湍流模型，采用定常、SIMPLE算法。网格采用四面体网格进行划分，总的网格量为820万。选取了包括额定风速在内的7种风速进行对比，对比结果如图2所示。研究结果表明，本文计算方法和CFD方法所得结果与叶片设计值相比，输出功率均较为接近，因此，本文方法具有一定的可靠性。

1.3 优化方法

为避免优化结果出现局部最优解，本文优化的设计变量将二维翼型与三维叶片的优化相结合，沿叶片展向选取表2给出的7个翼型截面位置的翼型形线及相应扭角为变量。采用Hicks-Henne^[12]型函数扰动法对各翼型形线进行参数化^[13]。以翼型的前缘点和尾缘点为界将翼型分为上下两条形线，在翼型的原始形状上叠加控制参数与型函数以实现翼型几何形状的变化。变化后的翼型形状表达通用式为

$ y(x) = {y_0}(x) + \sum\limits_{i = 1}^n {{c_i}{f_i}(x)} $

(12)

式中：x，y分别为翼型数据的横坐标和纵坐标；y₀为原始翼型的纵坐标；c_i为控制参数；${f_i}(x)$代表型函数的基函数。

$ {f}_{i}(x)=\left\{\begin{array}{*{20}{l}} {{x}^{0.25}(1-x){\text{e}}^{-20x},}&{i=1}\\ {{\mathrm{sin}}^{3}(\text{π}{x}^{e(i)}),}&{1 < i < n}\\ {{x}^{0.5}{(1-x)}^{0.1}{\text{e}}^{-20x},}&{i=n} \end{array} \right.$

(13)

$ e(i) = {\text{ln}}\;0.5/{\text{ln}}\;{x_i},\quad 0 \leqslant {x_i} \leqslant 1 $

(14)

由上述参数化方法，通过控制参数c_i的数值变化，结合扭角即可实现每步优化过程中翼型与叶片形状的改变。本文上下两条形线各选取5个控制点，令x_i分别为$0.1，0.3，0.5，0.7，0.9$ 。

以高风能利用率与低叶根弯矩为优化目标，进行叶片的多目标优化。考虑到多目标设计的复杂性，引入加权法将其转化为单目标设计^[14]。其目标函数可表示为

$ f(x) = \max ({u_1}P/{s_1} - {u_2}{M_{\rm{t}}}/{s_2}) $

(15)

式中：u表示权重因子；s表示比例因子，分别根据不同物理量的数量级选取。

本文对设计过程中所需要考虑到的几何参数、变量及目标进行约束，以确保优化结果的可信性与合理性，加快最优点搜寻速度。

为得到合理的翼型形线与叶片几何形状，对型函数控制参数c_i及扭角β的取值范围设定约束，其表达式为

$ \left\{ \begin{array}{*{20}{l}} {c_i} \leqslant 0.02 \\ {\beta _{\min }} \leqslant \cdots \leqslant {\beta _{i - 1}} \leqslant {\beta _i} \leqslant {\beta _{i + 1}} \leqslant \cdots \leqslant {\beta _{\max }} \end{array} \right. $

(16)

为防止翼型几何参数改变过度对叶片产生结构、强度方面的影响，约束各叶素处翼型的最大相对厚度T_m，使其与原始叶片厚度$T_{{\rm{m}}_0}$基本相当，其表达式为

$ \left| {{T_{\text{m}}} - {T_{{{\text{m}}_{\text{0}}}}}} \right| \leqslant 0.05 $

(17)

对输出功率与叶根弯矩同样给出约束，确保优化叶片的输出功率P不小于优化叶片的输出功率P₀；优化叶片的叶根弯矩M_t小于原始叶片的叶根弯矩$M_{{\rm{t}}_0} $。约束表达式为

$ \left\{ {\begin{array}{*{20}{l}} {P > {P_0}} \\ {{M_{\text{t}}} < {M_{{{\text{t}}_{\text{0}}}}}} \end{array}} \right. $

(18)

多岛遗传算法是一种具有较强可靠性的优化方法，曾被运用于叶片、翼型的气动优化^[15]。于上述基础上，采用多岛遗传算法对原始叶片进行优化。在得到效果最好的结果后，对超出额定风速工况下叶片的桨距角也进行优化，使其达到风力机的额定功率。图3给出了整体优化流程。

2 优化结果与分析 2.1 叶片几何

经前述优化流程，完成迭代计算后，得到优化叶片沿展向的扭角分布与7个叶素截面的翼型形线。表3给出了优化后各叶素截面翼型形线的几何形状与原始翼型的对比。表中：T_m与x_tm分别表示翼型的最大相对厚度与位置；C_am与x_cam分别表示翼型的最大弯度与位置。对比优化前后翼型几何情况可以看到，在经过遗传算法的寻优后，叶根处DU40，DU35钝尾缘翼型最大相对位置前移，叶中处DU30，DU25钝尾缘翼型与叶尖处NACA64尖尾缘翼型最大相对位置均后移。从叶根到叶尖处的钝、尖尾缘翼型的最大弯度均有所降低。其中，尖尾缘翼型与靠近叶尖处的部分钝尾缘翼型最大弯度位置前移，其余钝尾缘翼型最大弯度均后移。

图4给出了优化后的各叶素截面翼型形线与原始翼型形线的对比图。图中，r为叶片展向长度，弦长和厚度均为无量纲值。

图5给出了优化前后叶片沿展向安装角分布情况的对比。由图5可知，优化后的叶片安装角相较于原使叶片，在使用钝尾缘翼型的靠近叶根处呈减小趋势，在使用尖尾缘翼型的靠近叶尖处呈增大趋势。

图6给出了优化叶片在超出额定风速时的叶片桨距角与原始叶片的对比。由图6可知，优化叶片的桨距角相对原始叶片有所减小，最大减小角度约为1°。

2.2 优化性能分析

对优化叶片与原始叶片的气动性能进行了比较。由于翼型形状与叶片安装角分布发生了变化，在额定工况下，优化叶片基本保持了额定功率的输出，同时其叶根弯矩降低了约5%。

超出额定风速后，叶片将进行变桨以保证额定功率。因此对叶片变桨前的风能利用率、叶片变桨前与变桨后的叶根弯矩进行变风况性能的分析。

由图7可以看出，优化叶片的风能利用率相较于原始叶片在低风速下均有所提高，平均提高了约1%。叶片风能利用率的提高幅度随着风速增大而降低，在叶片设计工况下基本与原始叶片持平。

由图8可见，原始叶片与优化叶片的叶根弯矩均在额定工况下达到最大。其中，设计工况下的叶根弯矩相较于原始叶片降低了约5%。变桨前低风速下，相较于原始叶片，叶根弯矩平均降低了约5.8%；变桨后高风速下，相较于原始叶片，叶根弯矩平均降低了约4.4%。优化叶片于变桨前后，叶根弯矩均得到了一定程度的降低。

3 结　论

针对海上风力机等大型风电设备叶片功率高、负载大的特点，对NREL的5 MW风力机叶片的各截面翼型形线、扭角及变桨后的桨距角进行了多目标优化，得到了以下结果：

a. 相较于原始叶片，优化后的叶片在设计工况下输出功率基本与原始叶片保持一致，其叶根弯矩降低了5%。优化后的各截面翼型最大厚度位置、弯度及弯度位置改变，叶片扭角呈先增大后减小的趋势，变桨角减小。

b. 对优化叶片变桨前的气动效率与变桨前后的气动载荷进行变风况性能验证。相较于原始叶片，低风速下的优化叶片风能利用率平均提升了1%；低、高风速下的叶片弯矩分别平均降低了5.8%与4.4%。