一种水平轴风力机叶片的极限载荷预测计算方法

摘要

本发明涉及一种水平轴风力机叶片的极限载荷预测计算方法，采用智能优化算法进行极限载荷求解，选取风力机的转速、桨矩角、来流风速、偏航角和方位角为自由变量，建立叶片各截面各方向的载荷同上述自由变量的关系，根据风场类型和设计需要，对各自由变量进行约束，以叶片截面上的载荷为目标函数，使用智能优化算法（比如PSO算法）来求解极限载荷。本发明的方法具有快速而又准确的特点，且计算结果可以方便地应用于叶片铺层优化设计等问题中。

说明书

技术领域

本发明涉及风力机叶片的载荷预测计算，特别是水平轴风力机极限载荷的预测计算方法。 

技术背景

风力机的载荷源主要有以下几种：空气动力载荷、重力载荷、惯性载荷（包括离心和回转效应）、由于控制系统作用引起的运行载荷（例如，刹车、偏航、变桨距控制和发电机脱网等）。在叶片设计中，主要是考虑由这些载荷源引起的极限载荷和疲劳载荷。到目前为止,大部分风力机主要是因为各种极限状况的出现而失效的,严重的甚至无法修复。因此，极限载荷是叶片铺层设计时的重点。另外，国内的低风况区域分布广泛，随着我国风电产业的发展，设计出能适应低风况的叶片将是大势所趋。但此类叶片的设计与传统的叶片设计相比具有如下特点，叶片本身将变得更长，叶片也更柔，而与整机的接口参数如叶根的最大弯矩等又不允许增加，因此如何控制极限载荷以及开展极限载荷条件下的叶片设计已经成为急需解决的技术难题。而传统的叶片设计理念和方法，大部分均以叶片的最大气动效率，如最大能量输出、最小能量成本,最小质量等经济性指标为设计目标，进行叶片气动、结构设计及优化。但该设计思路无法解决长叶片带来的高载荷问题。所以建立叶片极限载荷的预测计算模型对解决以上问题具有重要意义。 

极限载荷预测计算一直是风力机优化设计过程中的一大难点，也是近几年来国外一直研究的热点。目前来说，主要有两类方法来计算预测叶片的极限载 荷，一类是工程上经常使用的方法，也是相对来说最准确的方法，它主要是依据GL或IEC标准定义的工况进行计算分析，对每一个工况都进行气弹模拟，然后根据计算得到的结果进行统计得到。而由于风力机的工作风速范围广，运行工况复杂多变，采用这种方法预测计算出叶片所受到的极限载荷将非常复杂耗时且不易用于叶片铺层的优化设计中。另一类方法主要是基于有限的载荷数据（测量或计算得到的）采用统计学分析方法和外插方法得到叶片的极限载荷，但它们存在以下问题，即如何为已知数据找到一个合适的概率分布函数；如何由短期的载荷分布得到长期的载荷分布；如何定义预测结果的不确定度。所以有必要建立一种新的极限载荷预测模型，它将具有快速而又准确的特点，为叶片的铺层设计甚至是气动设计提供帮助。 

发明内容

针对现有技术的缺点和不足，本发明的目的在于提供一种水平轴风力机叶片的极限载荷预测计算方法，该方法具有快速而又准确的特点，且计算结果可以方便地应用于叶片铺层优化设计等问题中。 

本发明为解决其技术问题所采取的技术方案为：一种水平轴风力机叶片的极限载荷预测计算方法，采用智能优化算法进行极限载荷求解，其特征在于，所述极限载荷预测计算方法包括如下步骤： 

1）选取风力机的转速Ω、桨矩角β₂、来流风速V₁、偏航角γ和方位角ψ为自由变量； 

2）建立叶片各截面各方向的载荷同上述自由变量的关系： 

（a）首先建立叶片的气动力F_a和上述自由变量的关系，将叶片分成多个互不相关的叶素，设每个叶素中各截面翼型、来流速度V₁、攻角α相同，叶素处的合成气流速度V₀作用在长度为dr的叶素上的气动力dF_a可分解为法向力dF_n和切向力dF_t，dF_n和dF_t可分别表示为， 

$\left(\begin{array}{c}{\mathrm{dF}}_n=\frac{1}{2}{{\mathrm{\rho V}}_0}^2{\mathrm{cC}}_n\mathrm{dr}\\ {\mathrm{dF}}_t=\frac{1}{2}{{\mathrm{\rho V}}_0}^2{\mathrm{cC}}_t\mathrm{dr}\end{array}\right),$

其中，ρ为空气密度，c为叶素剖面弦长，C_n、C_t分别表示叶素的法向力系数和切向力系数，其表达式为： 

$\left(\begin{array}{c}{C}_n=C_l\cos \phi +C_d\sin \phi \\ C_t=C_l\sin \phi -C_d\cos \phi \end{array}\right),$

其中，C_l、C_d分别为叶素的升力系数和阻力系数，φ为叶素处的入流角；叶素处的入流角φ和合成气流速度V₀可分别表示为： 

$\phi =\arctan \frac{(1-a)V_1}{(1+b)\mathrm{\Omega r}},$$V_0=\sqrt{{(1-a)}^2{V_1}^2+{(1+b)}^2{\left(\mathrm{\Omega r}\right)}^2},$

其中，V₁为叶素处的来流风速，r为叶素的旋转半径，a为叶素轴向诱导因子，b为叶素切向诱导因子； 

（b）叶素单位长度上的离心载荷 重力载荷 其中 为叶素单位长度上的质量，包括集中质量在该截面的分布； 

3）将离心载荷F_c和重力载荷F_g沿风轮坐标系进行分解，同气动载荷F_a一起作用在叶片上，计算叶片各截面各方向上的载荷； 

4）根据风场类型和设计需要，对各自由变量进行如下约束， 

$\left(\begin{array}{c}{\Omega }_{\mathrm{lower}}＜\Omega ＜{\Omega }_{\mathrm{upper}}\\ V_{\mathrm{lower}}＜V_1＜V_{\mathrm{upper}}\\ {\beta }_{2\_\mathrm{lower}}＜{\beta }_2＜{\beta }_{2\_\mathrm{upper}}\\ {\gamma }_{\mathrm{lower}}＜\gamma ＜{\gamma }_{\mathrm{upper}}\\ {\psi }_{\mathrm{lower}}＜\psi ＜{\psi }_{\mathrm{upper}}\end{array}\right),$

其中，Ω_upper、Ω_lower为转速Ω的上、下限，β_{2_upper}、β_{2_lower}为桨矩角β₂的上、下限，V_lower、V_upper为来流风速V₁的上、下限，γ_upper、γ_lower为偏航角γ的上、下限，ψ_upper、ψ_lower为方位角ψ的上、下限，这些参数由风力机实际运行的风场类型、整机参数和设计要求确定； 

5）以叶片某一截面的某一方向的载荷、某一截面的不同方向的载荷、叶片多个截面的某一方向的载荷、和/或多个截面的多个方向的载荷为目标函数，以步骤4）设定的约束条件，使用智能优化算法来求解相应载荷的极限值，得到极限载荷。 

进一步地，由于控制规律的作用，对不同风速和转速条件下所对应的最小桨距角β_{2_min}按下式来取值： 

$\left(\begin{array}{cc}{\beta }_{2\_\min }=\max ({\beta }_{2\_\min \_\mathrm{rs}},{\beta }_{2\_\min \_v})& (V_1＜V_{\mathrm{out}})\\ {\beta }_{2\_\min }={\beta }_{2\_\min \_v}& (V_1>V_{\mathrm{out}})\end{array}\right),$

其中，β_{2_min_v}为不同风速条件下对应的最小桨距角，β_{2_min_rs}为不同转速条件下对应的最小桨距角，V_out为风力机的切出风速；β_{2_min_v}、β_{2_min_rs}按照如下方式进行取值： 

（1）当风速V_in<V₁<V′时，不同风速条件下对应的最小桨距角β_{2_min_v}为β_{2_min_v}=β_{2_lower}，其中β_{2_lower}为叶片初始安装角，V_in为风力机的切入风速，V′为风力机开始变桨的风速； 

（2）当风速V′<V₁<V_out时，β_{2_min_v}与风速的对应关系采用以下二次多项式表示： 

β_{2_min_v}=B₁(V₁)²+B₂V₁+B₃-Δβ₂，其中Δβ₂为常数，且50<Δβ₂<150，B₁、B₂、B₃为二次多项式的系数，B₁、B₂、B₃为二次多项式的系数，B₁、B₂、B₃通过对风力机正常运行时不同风速时对应的桨矩角进行二次多项式拟合的方式得到。在实际操作时，可以采用如下方式确定：通过采用FOCUS软件、Bladed等其它气弹软件对风力机进行气弹模拟，我们可以得到风力机正常运行时不同风速时对应的桨矩角，将它们用(a_i,b_i)表示，i必须大于等于3，a_i表示风速值，b_i表示对应的桨矩角。a_i应取V′<V₁<V_out中所有模拟的风速点。根据这些点我们可以在EXCEL表格中作图，并采用回归分 析对这些数据进行二次多项式的拟合，当得到的二次多项式方程在V′<V₁<V_out范围内，其R平方值大于等于0.99时，我们取此时拟合得到的二次多项式方程中的二次方系数为B₁的值，一次方系数为B₂的值，常数项为B₃的值。如果采用回归分析得到的二次多项式，其在V′<V₁<V_out范围内，R平方值不能大于0.99，则将靠近风速V′和V_out附近的数据点逐个删除，直到满足R平方值大于等于0.99为止。 

如果计算得到的β_{2_min_v}≤β_{2_lower}时，取β_{2_min_v}=β_{2_lower}； 

（3）当风速V_out<V₁<V_upper时，β_{2_min_v}与风速的对应关系采用线性关系式表示如下： 

β_{2_min_v}=D₁(V₁-V_out)+D₂-Δβ₂，其中D₁、D₂为线性关系式的系数，采用如下方式确定D₁、D₂：该直线通过点（V_upper，90-Δβ₂）和点（V_out，β_{2_out}），其中β_{2_out}满足β_{2_out}=B₁(V_out)²+B₂V_out+B₃-Δβ₂。 

（4）当Ω≥Ω′时，Ω′为风力机额定转速附近一设定转速，β_{2_min_rs}与转速Ω的对应关系采用线性关系式表示为： 

β_{2_min_rs}=C₁Ω+C₂，其中C₁、C₂为线性关系式的系数，通过对风力机运行在额定风速时的一年一遇的极端操作阵风工况进行气弹模拟，得到风力机运行在不同转速时对应的桨矩角，对转速和桨矩角进行线性拟合得到C₁、C₂。在实际操作时，采用如下方式确定：通过采用FOCUS软件、Bladed等其它气弹软件对风力机运行在额定风速时的一年一遇的极端操作阵风工况（即EOG₁）进行气弹模拟，我们可以得到风力机运行在不同转速时对应的桨矩角，将它们用(m_i,n_i)表示，i必须大于等于2，m_i表示转速值，n_i表示对应的桨矩角。Ω′取为模拟过程中取得最小桨矩角值时对应的转速。设模拟过程中取得的最大转速为Ω″。则m_i应取Ω′≤Ω≤Ω″之 间所有模拟的转速点。根据这些点我们可以在EXCEL表格中作图，并采用回归分析对这些数据进行线性关系式的拟合，当得到的线性关系式在Ω≥Ω′范围内，其R平方值大于等于0.99时，我们取此时拟合得到的线性关系式中的一次方系数为C₁的值，常数项为C₂的值。如果采用回归分析得到的线性关系式，其在Ω≥Ω′范围内，R平方值不能大于0.99，则将靠近转速Ω′和Ω″附近的数据点逐个删除，直到满足R平方值大于等于0.99为止。 

（5）当转速Ω<Ω′时，β_{2_min_rs}=β_{2_lower}。 

优选的，采用如下方法计算叶素诱导因子a、b，包括如下计算步骤： 

i)设定诱导因子a、b初值：a=a₀，b=b₀； 

j)计算叶素的切向速度V_y0和法向速度V_x0：V_x0=V₁(1-a)，V_y0=Ωr(1+b)； 

k)计算叶素处的入流角φ和攻角α： 

$\phi =\arctan \frac{(1-a)V_1}{(1+b)\mathrm{\Omega r}},$α=φ-β₁-β₂， 

其中，β₁是叶素截面的扭角，β₂是叶素的桨矩角； 

l)计算叶片损失F，F=F_tipF_hub，叶尖损失F_tip、轮毂损失F_hub分别表示为： 

$F_{\mathrm{tip}}=\frac{2}{\pi }{\cos }^{-1}{e}^{-\left(\frac{B(R-r)}{2r\sin \phi }\right)},$$F_{\mathrm{hub}}=\frac{2}{\pi }{\cos }^{-1}{e}^{-\left(\frac{B(r-R_{\mathrm{hub}})}{2R_{\mathrm{hub}}\sin \phi }\right)},$

其中，B是叶片数，R_hub为轮毂半径，R是风轮的旋转半径； 

m)计算风轮的推力系数C_T， 其中 为风轮实度； 

n)求解新的切向诱导因子a： 

如果C_T≥0.96F，则该叶素载荷过高，新的轴向诱导因子a由下式求解： 

$a=\frac{1-B_1+\sqrt{{B_1}^2-4B_2(B_0-C_T)}}{2B_2},$

其中，B₂=1/0.18-4F，B₁=-(0.8/0.18-4F)，B₀=0.16/0.18； 

否则，新的轴向诱导因子a为$a={[1+\frac{4F{\sin }^2\phi }{{\mathrm{\sigma C}}_n}]}^{-1};$

o)求解新的切向诱导因子b，$b={[\frac{4F\sin \phi \cos \phi }{{\mathrm{\sigma C}}_t}-1]}^{-1};$

p)如果前后两次计算出来的诱导因子a、b变化大于某一容许偏差，则返回 

步骤b)循环运算，否则完成计算。 

优选的，采用如下方法计算叶素的初始诱导因子a₀，b₀，包括如下计算步骤： 

a）计算除叶根处的圆截面叶素外，其他所有叶素的入流角φ和合成气流速度V₀， 

$V_0=\mathrm{\Omega r}\sqrt{1+\frac{1}{{\lambda }_r^2}},({\lambda }_r>0),$

$\left(\begin{array}{cc}\phi =\frac{4F_{\mathrm{tip}}-{\mathrm{\sigma k}}_l{\lambda }_r+\sqrt{|{(4F_{\mathrm{tip}}-{\mathrm{\sigma k}}_l{\lambda }_r)}^2+16{\lambda }_r^2{\mathrm{\sigma F}}_{\mathrm{tip}}(k_l\theta -C_{l0})|}}{8{\lambda }_r{F}_{\mathrm{tip}}}& ({\lambda }_r\ge 1.5)\\ \phi =\arctan (1/{\lambda }_r)& ({\lambda }_r＜1.5)\end{array}\right),$

其中， 

$F_{\mathrm{tip}}=\frac{2}{\pi }\arccos \left(e^{-f}\right),$为叶尖损失， 

$f=\frac{B}{2}\frac{R-r}{r\sin {\phi }^{\prime }},$

${\phi }^{\prime }=\frac{4-{\mathrm{\sigma k}}_l{\lambda }_r+\sqrt{|{(1-{\mathrm{\sigma k}}_l{\lambda }_r)}^2+16{\lambda }_r^2\sigma (k_l\theta -C_{l0})|}}{8{\lambda }_r},$

θ=β₁+β₂，为叶素截面弦线与风轮平面的夹角， 

λ_r为叶素局部尖速比，σ为风轮实度，k_l为升力线斜率，C_l0为零攻角对应的升力系数，β₁是叶素截面的扭角，β₂是叶素的桨矩角； 

b）根据以上得到的入流角φ和合成气流速度V₀，计算除叶根处的圆截面叶素外的初始诱导因子a₀，b₀, 

$a_0=\frac{V_1-V_0\sin \phi }{V_1},$$b_0=\frac{V_0\cos \phi -\mathrm{\Omega r}}{\mathrm{\Omega r}}$

c）计算叶根处的圆截面叶素的入流角φ和合成气流速度V₀， 

φ＝arctan(1/λ_r)，$V_0=(1-a)\mathrm{\Omega r}\sqrt{1+\frac{1}{{\lambda }_r^2}},$

其中，a=σ/(4sinφ+σ)； 

d）根据以上得到的入流角φ和合成气流速度V₀，计算叶根处的圆截面叶素的初始诱导因子a₀，b₀, 

a₀=σ/(4sinφ+σ)，b₀=-a₀

在上述步骤a）、c）中：非偏航时，λ_t=Ωr/V₀；偏航时， V_{1_inplane}、V_{1_outplane}分别为来流风速V₁分解在风轮平面内的速度和垂直于风轮平面的速度。 

优选的，在计算叶素上的气动力dF_a时，也可采用如下方法计算叶片各截面各叶素截面的入流角φ和合成气流速度V₀，包括如下步骤： 

（1）当局部尖速比λ_r≥4时，根据权利要求4所述的步骤a）、c），针对不同形状的叶素截面直接计算叶素截面的入流角φ和合成气流速度V₀， 

（2）当局部尖速比λ_r<4时，根据权利要求3计算叶素截面的入流角φ和合成气流速度V₀。 

优选的，采用如下方法计算叶片各截面各方向上的载荷，包括如下步骤： 

（a）将离心载荷F_c和重力载荷F_g沿风轮坐标系进行分解，同气动载荷F_a一起作用在叶片上，将叶片分成N个叶素，N为大于2的整数，每个叶素对应一个节点，假设每一个叶素上的载荷是均匀分布的，各叶素上的均布载荷p^k=F_c+F_g+F_a，k为整数，1≤k≤N，设x方向为风轮旋转方向，y方向为垂直于 风轮旋转平面方向，z方向为旋转中心指向叶尖为正（不考虑预弯）；1号节点和N号节点对应的叶素长度dz¹、dz^N为零：dz^N=dz¹=0； 

（b）计算各叶素截面x方向上的力T_x^k和由此产生的力矩M_y^k： 

对于节点1：T_x¹=T_x²+0.5p²dz²，M_y¹=M_y²+0.5T_x²dz²+p²(dz²)²/8； 

对于节点N：T_x^N=0，M_y^N=0； 

对于节点N-1：T_x^N-1=T_x^N+0.5p^N-1dz^N-1，M_y^N-1=M_y^N+0.5T_x^Ndz^N-1+p^N-1(dz^N-1)²/8； 

对于节点N-2到2： 

T_x^k=T_x^k+1+0.5p^k+1dz^k+1+0.5p^kdz^k， 

M_y^k=M_y^k+1+0.5T_x^k+1(dz^k+1+dz^k）+0.5p^k+1dz^k+1(0.5dz^k+0.25dz^k+1)+p^k(dz^k)²/8； 

（c）按照和步骤（b）类似的方法计算各叶素截面y方向上的力T_y^k和由此产生的力矩M_x^k。 

优选的，所述智能优化算法为PSO算法（粒子群优化算法）。 

优选的，所述PSO算法为改进的PSO算法，在所述改进的PSO算法中： 

（a）惯性权值w按对数规律单调递减，其表达式为： 

$w=w_{\max }-\frac{\log (1+P)\times (w_{\max }-w_{\min })}{\log (1+\mathrm{gen})},$

其中，P为当前迭代次数，gen为最大迭代次数，w_max为最大惯性权值因子，w_min为最小惯性权值因子； 

（b）将每一代最优的前n个粒子直接复制到下一代，使每一代得到的最优解得以保存而不被破坏，其余的粒子正常进化后进入下一代； 

（c）将每一代中所有不满足约束条件的粒子参与进化生成下一代，但令这些粒子的目标函数值同为一个很小的常数，其余的粒子正常进化后进入下一代。 

优选的，n为粒子总数目的约10%。 

本发明的水平轴风力机叶片的极限载荷预测计算方法，主要特点在于： 

1）将叶片的极限载荷计算看成是求极值问题，采用智能优化算法，如PSO算法（粒子群算法）进行极限载荷求解； 

2）依据不同的叶片气动外形参数、风场类型、控制特性等参数进行载荷预测计算；

3）能准确快速的得到风力机叶片各截面的极限载荷值或是同时得到多个载荷分量的极限值，为叶片的铺层设计提供依据。 

附图说明

图1为叶素上的气流速度三角形和空气动力分量示意图。 

图2为叶片质量计算示意图。 

图3为叶片某一方向的载荷分布图。 

图4为叶片截面载荷计算流程图。 

图5为极限载荷预测计算流程图。 

具体实施方式

下面结合附图和具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。 

由于风力机的极限载荷预测也是一种求取极值问题，可以采用优化算法进行求解。所以我们将改进的PSO算法以及改进的多目标PSO算法用于风力机极限载荷的求解，建立了一种新的极限载荷预测模型。以下是整个极限载荷预测模型的建模过程。 

（一）自由变量选取 

对于风力机叶片所受的气动力F和气动力矩M，可以用以下函数表示。 

(F,M)＝f(γ,χ,t,β₂,Ω,β₁,c,V₁,ψ,ρ,C_l,C_d) 

式中，γ为偏航角，χ为锥角，t为倾角，β₂是桨矩角，Ω为风轮转速，β₁为叶片的截面扭角，c为叶片的截面弦长，V₁为来流风速（包括其大小和方向），ψ为方位角，ρ为空气密度，C_l和C_d分别为该截面翼型的升、阻力系数。 

由于叶片的气动载荷主要同风力机的转速Ω、桨矩角β₂、来流风速V₁、偏航角γ、方位角ψ等相关，所以选取转速、桨矩角、来流风速、偏航角和方位角为自由变量。 

（二）载荷计算方法 

1.气动载荷计算 

首先建立叶片的气动力F_a和上述自由变量的关系，将叶片分成多个互不相关的叶素，设每个叶素中各截面翼型、来流速度V₁、攻角α相同，图1为叶素上的气流速度三角形和空气动力分量示意图。叶素处的合成气流速度V₀作用在长度为dr的叶素上的气动力dF_a可分解为法向力dF_n和切向力dF_t，dF_n和dF_t可分别表示为， 

$\left(\begin{array}{c}{\mathrm{dF}}_n=\frac{1}{2}{{\mathrm{\rho V}}_0}^2{\mathrm{cC}}_n\mathrm{dr}\\ {\mathrm{dF}}_t=\frac{1}{2}{{\mathrm{\rho V}}_0}^2{\mathrm{cC}}_t\mathrm{dr}\end{array}\right),$

其中，ρ为空气密度，c为叶素剖面弦长，C_n、C_t分别表示叶素的法向力系数和切向力系数，其表达式为： 

$\left(\begin{array}{c}{C}_n=C_l\cos \phi +C_d\sin \phi \\ C_t=C_l\sin \phi -C_d\cos \phi \end{array}\right),$

其中，C_l、C_d分别为叶素的升力系数和阻力系数，φ为叶素处的入流角；叶素处的入流角φ和合成气流速度V₀可分别表示为： 

$\phi =\arctan \frac{(1-a)V_1}{(1+b)\mathrm{\Omega r}},$$V_0=\sqrt{{(1-a)}^2{V_1}^2+{(1+b)}^2{\left(\mathrm{\Omega r}\right)}^2},$

其中，V₁为叶素处的来流风速，r为叶素的旋转半径，a为叶素轴向诱导因子，b为叶素切向诱导因子。 

采用如下方法计算叶素诱导因子a、b，包括如下计算步骤： 

q)设定诱导因子a、b初值：a=a₀，b=b₀； 

r)计算叶素的切向速度V_y0和法向速度V_x0：V_x0=V₁(1-a)，V_y0=Ωr(1+b)； 

s)计算叶素处的入流角φ和攻角α： 

$\phi =\arctan \frac{(1-a)V_1}{(1+b)\mathrm{\Omega r}},$α=φ-β₁-β₂， 

其中，β₁是叶素截面的扭角，β₂是叶素的桨矩角； 

t)计算叶片损失F，F=F_tipF_hub，叶尖损失F_tip、轮毂损失F_hub分别表示为： 

$F_{\mathrm{tip}}=\frac{2}{\pi }{\cos }^{-1}{e}^{-\left(\frac{B(R-r)}{2r\sin \phi }\right)},$$F_{\mathrm{hub}}=\frac{2}{\pi }{\cos }^{-1}{e}^{-\left(\frac{B(r-R_{\mathrm{hub}})}{2R_{\mathrm{hub}}\sin \phi }\right)},$

其中，B是叶片数，R_hub为轮毂半径，R是风轮的旋转半径； 

u)计算风轮的推力系数C_T， 其中 为风轮实度； 

v)求解新的切向诱导因子a： 

如果C_T≥0.96F，则该叶素载荷过高，新的轴向诱导因子a由下式求解： 

$a=\frac{1-B_1+\sqrt{{B_1}^2-4B_2(B_0-C_T)}}{2B_2},$

其中，B₂=1/0.18-4F，B₁=-(0.8/0.18-4F)，B₀=0.16/0.18； 

否则，新的轴向诱导因子a为                                                  $a={[1+\frac{4F{\sin }^2\phi }{{\mathrm{\sigma C}}_n}]}^{-1};$

w)求解新的切向诱导因子b，$b={[\frac{4F\sin \phi \cos \phi }{{\mathrm{\sigma C}}_t}-1]}^{-1};$

x)如果前后两次计算出来的诱导因子a、b变化大于某一容许偏差，则返回步骤b)循环运算，否则完成计算。 

采用如下方法计算初始诱导因子a₀，b₀，包括如下计算步骤： 

a）计算除叶根处的圆截面叶素外，其他所有叶素的入流角φ和合成气流速度V₀， 

$V_0=\mathrm{\Omega r}\sqrt{1+\frac{1}{{\lambda }_r^2}},({\lambda }_r>0),$

$\left(\begin{array}{cc}\phi =\frac{4F_{\mathrm{tip}}-{\mathrm{\sigma k}}_l{\lambda }_r+\sqrt{|{(4F_{\mathrm{tip}}-{\mathrm{\sigma k}}_l{\lambda }_r)}^2+16{\lambda }_r^2{\mathrm{\sigma F}}_{\mathrm{tip}}(k_l\theta -C_{l0})|}}{8{\lambda }_r{F}_{\mathrm{tip}}}& ({\lambda }_r\ge 1.5)\\ \phi =\arctan (1/{\lambda }_r)& ({\lambda }_r＜1.5)\end{array}\right),$

其中， 

$F_{\mathrm{tip}}=\frac{2}{\pi }\arccos \left(e^{-f}\right),$为叶尖损失， 

$f=\frac{B}{2}\frac{R-r}{r\sin {\phi }^{\prime }},$

${\phi }^{\prime }=\frac{4-{\mathrm{\sigma k}}_l{\lambda }_r+\sqrt{|{(1-{\mathrm{\sigma k}}_l{\lambda }_r)}^2+16{\lambda }_r^2\sigma (k_l\theta -C_{l0})|}}{8{\lambda }_r},$

θ=β₁+β₂，为叶素截面弦线与风轮平面的夹角， 

λ_r为叶素局部尖速比，σ为风轮实度，k_l为升力线斜率，C_l0为零攻角对应的升力系数，β₁是叶素截面的扭角，β₂是叶素的桨矩角； 

b）根据以上得到的入流角φ和合成气流速度V₀，计算除叶根处的圆截面叶素外的初始诱导因子a₀，b₀, 

$a_0=\frac{V_1-V_0\sin \phi }{V_1},$$b_0=\frac{V_0\cos \phi -\mathrm{\Omega r}}{\mathrm{\Omega r}}$

c）计算叶根处的圆截面叶素的入流角φ和合成气流速度V₀， 

φ＝arctan(1/λ_r)，$V_0=(1-a)\mathrm{\Omega r}\sqrt{1+\frac{1}{{\lambda }_r^2}},$

其中，a=σ/(4sinφ+σ)； 

d）根据以上得到的入流角φ和合成气流速度V₀，计算叶根处的圆截面叶素 的初始诱导因子a₀，b₀, 

a₀=σ/(4sinφ+σ)，b₀=-a₀

在上述步骤a）、c）中：非偏航时，λ_t=Ωr/V₀；偏航时， V_{1_inplane}、V_{1_outplane}分别为来流风速V₁分解在风轮平面内的速度和垂直于风轮平面的速度。 

2.离心载荷和重力载荷的计算 

图2为叶片质量计算示意图。根据图2，单位长度上的离心载荷为， 其中 为单位长度上的质量，包括集中质量在该截面的分布。单位长度上的重力载荷为， 

将离心载荷和重力载荷沿风轮坐标系进行分解，同气动载荷一起作用在叶片上。假设每一个叶素上的载荷是均匀分布的，如图3所示，x方向为风轮旋转方向，y方向为垂直于风轮旋转平面方向，它是垂直于纸面向里，z方向为旋转中心指向叶尖为正（不考虑预弯）。在图3中，均布载荷p^k为p^k=F_c+F_g+F_a。1号节点和N号节点的对应的叶素长度为零，即dz^N=dz¹=0。然后在此基础上计算各截面的力T_x^k和由此产生的力矩M_y^k： 

对于节点1：T_x¹=T_x²+0.5p²dz²，M_y¹=M_y²+0.5T_x²dz²+p²(dz²)²/8； 

对于节点N：T_x^N=0，M_y^N=0； 

对于节点N-1：T_x^N-1=T_x^N+0.5p^N-1dz^N-1，My^N-1=M_y^N+0.5T_x^Ndz^N-1+p^N-1(dz^N-1)²/8； 

对于节点N-2到2： 

T_x^k=T_x^k+1+0.5p^k+1dz^k+1+0.5p^kdz^k， 

M_y^k=M_y^k+1+0.5T_x^k+1(dz^k+1+dz^k）+0.5p^k+1dz^k+1(0.5dz^k+0.25dz^k+1)+p^k(dz^k)²/8。 

对于y方向，其计算方法类似。 

综合上面有关气动载荷、离心载荷和重力载荷的计算方法，得到叶片截面载荷计算流程如图4所示。 

（三）约束条件处理 

以某一1.5MW风力机为计算对象，根据风场类型、整机参数和设计需要，对各自由变量进行如下约束， 

$\left(\begin{array}{c}9.00＜\Omega ＜27.14\\ 3.00＜V_1＜38.26\\ -1.0＜{\beta }_2＜90.0\\ -8.0＜\gamma ＜8.0\\ 0.0＜\psi ＜120.0\end{array}\right)$

另外，由于控制规律的作用，对不同风速和转速条件下所对应的桨距角取值范围进行如下约束： 

β₂≥β_{2_min}， 

其中β_{2_min}是不同风速和转速条件下所对应的最小桨距角，它采用如下步骤确定： 

（1）当风速3≤V₁≤11m/s时，最小桨矩角与风速的对应关系采用下式表示 

β_{2_min_v}=β_{2_lower}=-1， 

（2）当风速11<V₁≤25m/s时，最小桨矩角与风速的对应关系采用以下二次多项式表示。 

β_{2_min_v}=-0.0579(V₁)²+3.5789V₁-30.611-Δβ₂， 

其中Δβ₂=10.86 

由于桨矩角始终不能小于β_{2_lower}，所以有， 

β_{2_min_v}=β_{2_lower}(β_{2_min_v}<β_{2_lower})， 

（3）当风速25<V₁≤38.26m/s时，将最小桨矩角与风速的对应关系表示如下， 

β_{2_min_v}=5.0908(V₁-25)+22.674-Δβ₂， 

此外，桨矩角总是根据转矩来控制的，而转矩与转速直接相关，所以当风轮的转速增加时，其对应的最小桨矩角也会增加。 

（4）当Ω≥17.4时，最小桨矩角与转速的对应关系表示如下， 

β_{2_min_rs}=1.9969Ω-35.681， 

本模型中，根据额定风速时的一年一遇的极端操作阵风（EOG）工况确定不同转速时，转速与其对应的最小桨矩角之间的关系并拟合得到上式的线性表达式。 

（5）当Ω<Ω′时，最小桨矩角与转速的对应关系表示如下， 

β_{2_min_rs}=β_{2_lower}=-1 

（6）综合以上各式，得到不同风速和转速条件下所对应的最小桨矩角β_{2_min}为， 

$\left(\begin{array}{cc}{\beta }_{2\_\min }=\max ({\beta }_{2\_\min \_\mathrm{rs}},{\beta }_{2\_\min \_v})& (V_1＜25)\\ {\beta }_{2\_\min }={\beta }_{2\_\min \_v}& (V_1>25)\end{array}\right)$

需要注意的是，由于风力机以及所适用的风场类型不一样，约束条件的具体数值会有所不同，但都采用以上的方法进行确定。 

（四）极限载荷计算流程 

由于叶片各截面各方向的极限载荷不一定同时出现，所以需要对某一截面的某一方向的极限载荷进行单独求解。如在计算风力机的载荷包络线时，常选取某一截面的弦向载荷分量进行求解。有时为了实际需要，也要求对不同截面的不同方向的极限载荷进行同时求解。如在校核分析截面的极限强度是否满足设计要求时，需要对截面内两个互相垂直方向的载荷分量进行同时求解。可以采用各种智能优化算法来进行叶片极限载荷的求解，比如PSO算法，也可以对PSO算法进行改进，采用改进的PSO算法。 

这里介绍一种改进的PSO算法，在保留传统PSO算法的主要特点外，对其进行如下改进： 

（a）惯性权值w按对数规律单调递减，其表达式为： 

$w=w_{\max }-\frac{\log (1+P)\times (w_{\max }-w_{\min })}{\log (1+\mathrm{gen})},$

其中，P为当前迭代次数，gen为最大迭代次数，w_max为最大惯性权值因子，w_min为最小惯性权值因子； 

（b）将每一代最优的前n个粒子直接复制到下一代，使每一代得到的最优解得以保存而不被破坏，其余的粒子正常进化后进入下一代； 

（c）将每一代中所有不满足约束条件的粒子参与进化生成下一代，但令这些粒子的目标函数值同为一个很小的常数，其余的粒子正常进化后进入下一代。 

最后将以上模型同改进的PSO算法一起结合，建立一种新的风力机叶片极限载荷预测模型。其计算流程如图5所示。对于多个载荷分量同时求解的情况，与此类似，只是在得到目标函数值以后增加了pareto解集的筛选过程。 

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。