CN102750410B - 一种水平轴风力机叶片铺层的优化设计方法 - Google Patents

Abstract

本发明涉及一种极限载荷条件下叶片铺层的优化设计方法，包括：对叶片的初始铺层进行参数化建模；得到叶片展向各截面的铺层信息；计算叶片的质量线密度分布和刚度分布；根据所述质量线密度分布和刚度分布计算叶片的固有频率；根据叶片各截面的极限载荷分布及叶片的结构特性参数求取叶片的最小极限安全因子以及最大叶尖偏移；判断叶片的固有频率、最小极限安全因子以及最大叶尖偏移是否满足约束条件；采用智能优化算法对叶片铺层进行优化，得到最优的叶片铺层方案。本发明的具有：能对不同的叶片铺层进行优化设计，可根据不同的设计目标和要求来改变所需优化设计的铺层参数，设计速度快，优化后的铺层可以直接用于实际的叶片生产等众多优点。

Description

一种水平轴风力机叶片铺层的优化设计方法

技术领域

本发明涉及水平轴风力机叶片铺层的设计技术领域，特别是水平轴风力机极限载荷的预测计算方法以及极限载荷条件下的风力机叶片铺层优化设计方法。

技术背景

现代大型风力机主要以水平轴的结构型式为主。叶片作为其关键部件之一，它的性能直接关系到整个风电机组的风能转换效率和制造成本，而现在的水平轴风力机叶片主要是由复合材料铺层制成。所以开展风力机叶片铺层的优化设计对降低风电机组成本及提高叶片的可靠性具有重要意义。

风力机叶片的铺层设计需要考虑的设计要求很多，如叶片的质量、叶片的固有频率、叶根的极限载荷、叶片的最大叶尖偏移、叶片的最小极限安全因子等。其中叶根的极限载荷、叶片的最大叶尖偏移和叶片的最小极限安全因子与叶片受到的载荷直接相关，尤其是与叶片受到的极限载荷相关。到目前为止,大部分风力机主要是因为各种极限状况的出现而失效的,严重的甚至无法修复。因此，极限载荷是叶片铺层设计时的重点。另外，国内的低风况区域分布广泛，随着我国风电产业的发展，设计出能适应低风况的叶片将是大势所趋。但此类叶片的设计与传统的叶片设计相比具有如下特点，叶片本身将变得更长，叶片也更柔，而与整机的接口参数如叶根的最大弯矩等又不允许增加，因此如何控制极限载荷以及开展极限载荷条件下的叶片设计已经成为急需解决的技术难题。而传统的叶片设计理念和方法，大部分均以叶片的最大气动效率，如最大能量输出、最小能量成本,最小质量等经济性指标为设计目标，进行叶片气动、结构设计及优化。

在得到风力机的极限载荷后，由于叶尖偏移与载荷成正比，安全因子与载荷成反比，可以将叶片各截面的极限载荷同时加载到叶片上，进行最大叶尖偏移计算和最小极限安全因子计算评估。而工程设计上使用的方法主要还是采用多工况的气弹模拟计算，然后统计得到叶片的最大叶尖偏移和最小极限安全因子。这样做虽然比较准确但是非常耗时。而且如果叶片的铺层不满足其中一项要求，则需要重新修改铺层，重新进行各工况的气弹模拟计算，这将大大增加计算时间，而事实上有一些计算是可以不用做的。同时也不易用于叶片铺层的优化设计中。

目前，工程上优化设计叶片铺层所使用的方法主要依据设计者的经验手动的调整铺层、然后计算该叶片铺层是否满足以上的设计要求，不满足再返回来手动的调整铺层，再计算，再调整…，直到满足设计要求为止。所以其设计时间长，而且往往可以继续优化，得到性能更好的叶片。为了减少设计时间又能得到性能更好的叶片，也有一些学者将优化算法用于叶片结构的优化设计中，但他们在进行结构分析计算时，通常选取有限的几个工况进行分析，由于极限载荷、最大叶尖偏移、最小安全因子不一定出现在选取的工况中，所以无法保证优化后的叶片运行在所有工况下都满足设计要求。需要进行全工况的校核分析，如果不满足设计要求还需要重新优化、重新校核…，直到满足设计要求为止。虽然能在一定程度上减少设计时间同时获得较好结果，但还需要对以上优化设计方法进一步的改进来更快更好的设计叶片，以使设计出的叶片具有更高的市场竞争力。

发明内容

针对目前使用的设计方法存在的设计时间长且不易用于优化设计的问题，本发明的目的在于提供一种水平轴风力机叶片铺层的优化设计方法。

本发明为解决其技术问题所采取的技术方案是：

一种极限载荷条件下叶片铺层的优化设计方法，其特征在于，所述优化设计方法包括如下步骤：

步骤S1、获得风力机叶片的初始铺层方案，对叶片的初始铺层进行参数化建模；

步骤S2、根据参数化建模后的叶片铺层，得到叶片展向各截面的铺层信息；

步骤S3、根据叶片展向各截面的铺层信息，采用经典层合板理论并假设周向应力为零，计算叶片的质量线密度分布和刚度分布；

步骤S4、根据叶片的质量线密度分布和刚度分布计算叶片的固有频率；

步骤S5、根据叶片各截面的极限载荷分布及叶片的结构特性参数(主要是叶片截面的挥舞刚度、摆振刚度、扭转刚度、抗拉刚度、质量线密度、质量中心和弹性中心参数)求取叶片的最小极限安全因子以及最大叶尖偏移；

步骤S6、判断叶片的固有频率、最小极限安全因子以及最大叶尖偏移是否满足约束条件，如果满足则沿叶片分段积分得到叶片的总质量；否则，直接令叶片质量为一很小的常数值，如-99999。

步骤S7、采用智能优化算法对叶片铺层进行优化，改变叶片的铺层，继续步骤S2至步骤S6，直至得到最优的叶片铺层方案。

优选地，步骤S3中，采用经典层合板理论并假设周向应力为零计算叶片各截面的刚度特性、质量线密度，得到叶片沿展向的质量线密度分布和刚度分布。

优选地，步骤S4中，将叶片简化成叶根固定在轮毂的悬臂梁，根据叶片模态计算方法得到质量矩阵和刚度矩阵

${\stackrel{\‾}{M}}_{\mathrm{ij}}={\∫}_0^L\stackrel{\·}{m}\left(x\right){\mathrm{\ψ}}_i\left(x\right){\mathrm{\ψ}}_j\left(x\right)\mathrm{dx},$

${\stackrel{\‾}{K}}_{\mathrm{ij}}={\∫}_0^L{\mathrm{EI}}_F\left(x\right){\mathrm{\ψ}}_i^{\″}\left(x\right){\mathrm{\ψ}}_j^{\″}\left(x\right)\mathrm{dx}+{\∫}_0^L{\mathrm{\ψ}}_i^{\′}\left(x\right){\mathrm{\ψ}}_j^{\′}\left(x\right){\∫}_x^L\left(\stackrel{\·}{m}\left(x\right){\mathrm{\Ω}}^{2}r\right)\mathrm{drdx},$

其中，其中，为截面质量线密度；EI_F(x)为弹性主轴坐标系下的截面挥舞刚度；L为叶片单元段的长度；Ω为风轮转速；r为截面的旋转半径；ψ_i(x)和ψ_j(x)为形状函数，分别取为xⁱ⁺¹和x^j+1，式中1≤i,j≤n,n为大于等于5的整数，一般取为5。

计算由质量矩阵和刚度矩阵构成的常系数方程的特征值，得到叶片的挥舞频率。

优选地，步骤S5中，计算叶片各截面的最大载荷，将各截面的最大载荷同时加载到叶片上，通过静力学分析得到叶片的最大叶尖偏移及截面各材料承受的最大应力σ_max(x_E,y_E)；设叶片截面上(x_E,y_E)位置处，材料的许用拉应力为[σ_t]、许用压应力为[σ_c]，材料的安全系数为γ_m，则相应的最小极限拉应力安全因子SF_{t_min}、最小极限压应力安全因子SF_{c_min}分别为

${\mathrm{SF}}_{t\_\min }=\frac{\left[{\mathrm{\σ}}_t\right]}{{\mathrm{\γ}}_m\mathrm{\σ}(x_E,y_E)},{\mathrm{SF}}_{c\_\min }=\frac{\left[{\mathrm{\σ}}_c\right]}{{\mathrm{\γ}}_m\mathrm{\σ}(x_E,y_E)}.$

优选地，步骤S7中，采用智能优化算法对叶片铺层进行优化时，目标函数根据不同的叶片设计目标和要求进行设定，可以为叶片质量、叶片一阶挥舞频率、叶片最大叶尖偏移和最小极限安全因子等。

优选地，步骤S7中，采用智能优化算法对叶片铺层进行优化时，

a)选取主、副梁帽铺层的展向位置x和铺层厚度t做为设计变量，不考虑梁帽铺层的宽度变化；

b)以叶片质量最小为目标函数；

c)约束条件：

c1、叶片的最大叶尖偏移d_{tip_max}满足：d_{tip_max}≤d_{tip_set}，

其中，d_{tip_set}为某一设定值，该值满足d_{tip_set}≤d_{tip_max0}+[d]-d_{tip_act}，

其中，[d]为叶片的实际最大叶尖偏移的允许值，d_{tip_act}为初始叶片铺层时叶片实际的最大叶尖偏移，d_{tip_max0}为初始叶片铺层时根据叶片各截面的极限载荷分布计算得到的最大叶尖偏移；

c2、叶片的一阶挥舞频率Flap_1(x,t)满足：

|Flap_1(x,t)-(Flap_1)_set|≤Δ

其中，(Flap_1)_set是设定的一阶挥舞频率，取为大于额定转速三倍频的10％的值；Δ是容许度，取为额定转速三倍频的5％的值。

c3、叶片的最小极限安全因子SF_min(为最小拉应力安全因子SF_{t_min}或最小压应力安全因子SF_{c_min})满足：

|SF_min|≥SF_set≥1

其中，SF_set为某一设定的安全系数值，一般取为初始叶片最小极限安全因子与其最小疲劳安全因子的比值。如果没有可以参考的初始叶片值，则直接取为1.1；

c4、对梁帽铺层的展向位置x和铺层厚度t进行几何约束：

$\left\{\begin{array}{c}x\≤x\≤x_{\mathrm{up}}\\ t_{\mathrm{down}}\≤t\≤t_{\mathrm{up}}\end{array}\right.$

其中，梁帽的铺层厚度t还必须为单层材料厚度t₀的整数倍。

优选地，所述智能优化算法为PSO算法或遗传算法。

优选地，所述PSO算法为改进的PSO算法，在所述改进的PSO算法中：(a)惯性权值w按对数规律单调递减，其表达式为：

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

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

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

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

进一步地，所述改进的PSO算法中，n为粒子总数目的约10％。

同现有技术相比，本发明的极限载荷条件下叶片铺层的优化设计方法具有以下显著的优点：

1)能对不同的水平轴风力机叶片铺层进行优化设计，且可以根据不同的设计目标和要求来改变所需优化设计的铺层参数；

2)设计速度快，能够考虑大部分的叶片铺层设计要求，如质量、频率、最大叶尖偏移、叶片的极限强度要求；

3)优化后的铺层可以直接用于实际的叶片生产制造。

附图说明

图1为叶片截面铺层信息示意图。

图2材料主坐标系与参考坐标系。

图3叶片截面铺层几何示意图。

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

图5为载荷变换坐标系。

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

图7为几何变换坐标系。

图8为梁帽厚度沿叶片展向分布示意图。

图9为减重优化设计程序流程图。

具体实施方式

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

本发明的极限载荷条件下叶片铺层的优化设计方法，主要包括7个步骤：步骤S1、获得风力机叶片的初始铺层方案，对叶片的初始铺层进行参数化建模；步骤S2、根据参数化建模后的叶片铺层，得到叶片展向各截面的铺层信息；步骤S3、根据叶片展向各截面的铺层信息，采用经典层合板理论并假设周向应力为零计算叶片的质量线密度分布和刚度分布；步骤S4、根据叶片的质量线密度分布和刚度分布计算叶片的固有频率；步骤S5、根据叶片各截面的极限载荷分布及叶片的结构特性参数(主要是叶片截面的挥舞刚度、摆振刚度、扭转刚度、抗拉刚度、质量线密度、质量中心和弹性中心参数)求取叶片的最小极限安全因子以及最大叶尖偏移；步骤S6、判断叶片的固有频率、最小极限安全因子以及最大叶尖偏移是否满足约束条件，如果满足则沿叶片分段积分得到叶片的总质量；否则，直接令叶片质量为一很小的常数值，如-99999；步骤S7、采用智能优化算法对叶片铺层进行优化，改变叶片的铺层，继续步骤S2至S6，直至得到最优的叶片铺层方案。

对于步骤S1，关于对叶片的铺层进行参数化建模，现有技术中具有多种成熟的做法，这里简单介绍一种铺层参数化建模的方法。如图1所示，将叶片截面被分成很多sector，每个sector包含有相同的铺层。根据叶片的实际铺层情况，可以参考商用软件FOCUS的输入文件样式，编写参数化的文本文件。它主要由三个文件组成：铺层定义文件，线定义文件和材料特性文件。铺层定义文件给出每一种铺层在展向的位置，厚度分布，所采用的材料以及铺层在周向的边线。线定义文件主要是用来定义所有铺层的边线，通过它能找到铺层在任意展向截面的周向位置。材料特性文件就是给出叶片所用材料的特性参数，如密度、模量等。

进行步骤S2时，需根据步骤S1参数化建模后的叶片铺层，得到叶片展向各截面的铺层信息，即根据步骤S1中产生的三个文件进行数据处理以得到我们关心的一些截面的铺层信息。它有一个假设，那就是假设铺层厚度在两已知截面内线性分布，以方便未知截面的数据处理。

进行步骤S3时，根据叶片展向各截面的铺层信息，采用经典层合板理论并假设周向应力为零并假设周向应力为零，计算叶片的质量线密度分布和刚度分布。

对于单位宽度的单层板，当纤维主轴同xoy坐标轴成θ角时，以逆时针为正(如图2所示)。则单层折算刚度系数Q_ij(i,j＝1,2,6)可以表示成：

$\left\{\begin{array}{c}{Q}_{11}=Q_{11}^{\′}{\cos }^4\mathrm{\θ}++2(Q_{12}^{\′}+{2Q}_{66}^{\′}){\sin }^2\mathrm{\θ}{\cos }^2\mathrm{\θ}+Q_{22}^{\′}{\sin }^4\mathrm{\θ}\\ Q_{12}=(Q_{11}^{\′}+Q_{22}^{\′}-{4Q}_{66}^{\′}){\sin }^2\mathrm{\θ}{\cos }^2\mathrm{\θ}+Q_{12}^{\′}({\cos }^4\mathrm{\θ}+{\sin }^4\mathrm{\θ})\\ Q_{22}=Q_{11}^{\′}{\sin }^4\mathrm{\θ}++2(Q_{12}^{\′}+{2Q}_{66}^{\′}){\sin }^2\mathrm{\θ}{\cos }^2\mathrm{\θ}+Q_{22}^{\′}{\cos }^4\mathrm{\θ}\\ Q_{16}=(Q_{11}^{\′}-Q_{12}^{\′}-{2Q}_{66}^{\′})\mathrm{sim\θ}{\cos }^3\mathrm{\θ}+(Q_{11}^{\′}-Q_{12}^{\′}+{2Q}_{66}^{\′}){\sin }^3\mathrm{\θ}\cos \mathrm{\θ}\\ Q_{26}=(Q_{11}^{\′}-Q_{12}^{\′}-{2Q}_{66}^{\′}){\sin }^3\mathrm{\θ}\cos \mathrm{\θ}+(Q_{11}^{\′}-Q_{12}^{\′}+{2Q}_{66}^{\′})\sin \mathrm{\θ}{\cos }^3\mathrm{\θ}\\ Q_{66}=(Q_{11}^{\′}+Q_{22}^{\′}-{2Q}_{12}^{\′}-{2Q}_{66}^{\′}){\sin }^2\mathrm{\θ}{\cos }^2\mathrm{\θ}+Q_{66}^{\′}({\cos }^4\mathrm{\θ}+{\sin }^4\mathrm{\θ})\end{array}\right.$

式中，

$\left\{\begin{array}{c}{Q}_{11}^{\′}=E_1/(1-v_{12}{v}_{21})\\ Q_{12}^{\′}=v_{12}{E}_2/(1-v_{12}{v}_{21})\\ Q_{22}^{\′}=E_2/(1-v_{12}{v}_{21})\\ Q_{22}^{\′}=G_{12}\end{array}\right.$

所以，由以上的四个独立的弹性常数，1方向的弹性模量E₁、2方向的弹性模量E₂、平面内的泊松比ν₁₂、面内剪切模量G₁₂和正交各向异性材料的几何参数可计算出单层折算刚度系数。

在此基础上，计算单位宽度多层板的刚度表达式如下(如图3所示)。

$A_{\mathrm{ij}}=\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{\left({\stackrel{\‾}{Q}}_{\mathrm{ij}}\right)}_k(z_k-z_{k-1})=\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{\left({\stackrel{\‾}{Q}}_{\mathrm{ij}}\right)}_k{t}_k$

$B_{\mathrm{ij}}=\frac{1}{2}\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{\left({\stackrel{\‾}{Q}}_{\mathrm{ij}}\right)}_k({z_k}^2-{z_{k-1}}^2)=\frac{1}{2}\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{\left({\stackrel{\‾}{Q}}_{\mathrm{ij}}\right)}_k{t}_k{\stackrel{\‾}{z}}_k$

$D_{\mathrm{ij}}=\frac{1}{3}\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{\left({\stackrel{\‾}{Q}}_{\mathrm{ij}}\right)}_k({z_k}^3-{z_{k-1}}^3)=\frac{1}{3}\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{\left({\stackrel{\‾}{Q}}_{\mathrm{ij}}\right)}_k{t}_k({{\stackrel{\‾}{z}}_k}^2+\frac{{t_k}^2}{12})$

其中，A_ij、B_ij、D_ij(i＝1,2,6)为拉伸刚度、耦合刚度及弯曲刚度。t_k为第k层板的厚度，第k层中心点与中面的垂直距离， ${\stackrel{\‾}{Q}}_{16}=Q_{16}-\frac{Q_{12}{Q}_{26}}{Q_{22}};{\stackrel{\‾}{Q}}_{61}=Q_{61}-\frac{Q_{62}{Q}_{21}}{Q_{22}};{\stackrel{\‾}{Q}}_{66}=Q_{66}-\frac{Q_{62}{Q}_{26}}{Q_{22}}.$

计算单位宽度多层板的质量线密度的表达式如下(如图3所示)。

${\stackrel{\‾}{m}}_k={\mathrm{\ρ}}_k{t}_k$

其中，ρ_k为第k层材料的密度。

将单位宽度多层板的刚度及质量线密度沿叶片截面的周向积分即可得到整个截面的刚度和质量线密度。

进行步骤S4时，根据叶片的质量线密度分布和刚度分布计算叶片的固有频率，即通过对叶片进行模态分析得到叶片的固有频率。到目前为止，叶片模态计算的方法已经比较成熟，如直接在ANSYS中铺层建模计算、将叶片简化为直梁模型、还有锥梁模型计算等。将叶片简化成悬臂梁模型，然后采用梁单元求解既能大大减少计算时间，而且其计算精度能满足工程设计需要。本发明中，优选将叶片简化成叶根固定在轮毂的悬臂梁。在只考虑旋转平面外的振动时，叶片旋转的动能、势能为：

$T=\frac{1}{2}{\∫}_0^L\stackrel{\·}{m}\left(x\right){\left({\stackrel{\·}{w}}_F\right)}^2\mathrm{dx},U=\frac{1}{2}{\∫}_0^L{\mathrm{EI}}_F\left(x\right){\left(w_F^{\″}\right)}^2\mathrm{dx}+\frac{1}{2}{\∫}_0^L{\left(w_F^{\′}\right)}^2{\∫}_x^L\left(\stackrel{\·}{m}\left(x\right){\mathrm{\Ω}}^{2}r\right)\mathrm{drdx}.$

在不考虑外载和阻尼影响时，系统机械能守恒：令叶片位移得叶片振动的微分方程为： $M\stackrel{\·\·}{q}+\mathrm{Kq}=0,$ 其中，

q(t)仅是时间t的函数。为一仅与x相关而与时间无关的函数。它需要满足位移边界条件。为了降低质量矩阵和刚度矩阵的大小，采用瑞利-里兹法进行频率计算，将取为n+1次多项式，且满足叶根固定连接条件。表达式如下这样，质量矩阵M和刚度矩阵K可进一步表示为n维的质量矩阵和刚度矩阵且由新的刚度矩阵和质量矩阵得到的频率与原特征值相等。与中各元素的表达式为，

${\stackrel{\‾}{M}}_{\mathrm{ij}}={\∫}_0^L\stackrel{\·}{m}\left(x\right){\mathrm{\ψ}}_i\left(x\right){\mathrm{\ψ}}_j\left(x\right)\mathrm{dx},{\stackrel{\‾}{K}}_{\mathrm{ij}}={\∫}_0^L{\mathrm{EI}}_F\left(x\right){\mathrm{\ψ}}_i^{\″}\left(x\right){\mathrm{\ψ}}_j^{\″}\left(x\right)\mathrm{dx}+{\∫}_0^L{\mathrm{\ψ}}_i^{\′}\left(x\right){\mathrm{\ψ}}_j^{\′}\left(x\right){\∫}_x^L\left(\stackrel{\·}{m}\left(x\right){\mathrm{\Ω}}^{2}r\right)\mathrm{drdx}$

在得出新的质量矩阵和刚度矩阵后，叶片的挥舞频率的求解就是计算由质量矩阵和刚度矩阵构成的常系数方程的特征值。而求解特征值的方法很多。

进行步骤S5时，根据叶片各截面的极限载荷分布及叶片的结构特性参数(主要是叶片截面的挥舞刚度、摆振刚度、扭转刚度、抗拉刚度、质量线密度、质量中心和弹性中心参数)求取叶片的最小极限安全因子以及最大叶尖偏移。

由于叶尖偏移与叶片所受的载荷成正比，载荷越大，其叶尖偏移也越大。因此，要得到叶片的最大叶尖偏移，只需得到叶片各截面的最大载荷即可。计算叶片各截面的最大载荷，然后将各截面的最大载荷同时加载到叶片上，最后通过静力学分析得到叶片的最大叶尖偏移。下面是叶尖变形的具体计算方法：

(1)如图5所示，将叶片坐标系下的弯矩转换到截面的弹性主轴坐标系下。则有，

$\left\{\begin{array}{c}{M}_{\mathrm{YE}}=M_{\mathrm{YB}}\cos ({\mathrm{\β}}_1+\mathrm{\υ})+M_{\mathrm{XB}}\sin ({\mathrm{\β}}_1+\mathrm{\υ})\\ M_{\mathrm{XE}}=-M_{\mathrm{YB}}\sin ({\mathrm{\β}}_1+\mathrm{\υ})+M_{\mathrm{XB}}\cos ({\mathrm{\β}}_1+\mathrm{\υ})\end{array}\right.$

其中M_XB、M_YB为叶片坐标系下的弯矩；M_XE、M_YE为截面弹性主轴坐标系下的弯矩。β₁为截面扭角，υ为弹性主轴坐标系与截面弦向坐标系的夹角。

(2)材料力学中的梁理论可以得到关于弹性主轴的曲率为，

${\mathrm{\κ}}_{\mathrm{YE}}=\frac{M_{\mathrm{YE}}}{{\mathrm{EI}}_E},{\mathrm{\κ}}_{\mathrm{XE}}=\frac{M_{\mathrm{XE}}}{{\mathrm{EI}}_F}.$

(3)将这些曲率转换回叶片坐标系，得到关于叶片坐标系下的曲率为，

$\left\{\begin{array}{c}{\mathrm{\κ}}_{\mathrm{XB}}={\mathrm{\κ}}_{\mathrm{XE}}\cos ({\mathrm{\β}}_1+\mathrm{\υ})+{\mathrm{\κ}}_{\mathrm{YE}}\sin ({\mathrm{\β}}_1+\mathrm{\υ})\\ {\mathrm{\κ}}_{\mathrm{YB}}=-{\mathrm{\κ}}_{\mathrm{XE}}\sin ({\mathrm{\β}}_1+\mathrm{\υ})+{\mathrm{\κ}}_{\mathrm{YE}}\cos ({\mathrm{\β}}_1+\mathrm{\υ})\end{array}\right.,$

此外，角位移、位移与曲率之间的关系如下，

$\frac{{\mathrm{d\θ}}_{\mathrm{XB}}}{\mathrm{dz}}={\mathrm{\κ}}_{\mathrm{XB}},\frac{{\mathrm{d\θ}}_{\mathrm{YB}}}{\mathrm{dz}}={\mathrm{\κ}}_{\mathrm{YB}},\frac{{\mathrm{du}}_{\mathrm{YB}}}{\mathrm{dz}}={\mathrm{\θ}}_{\mathrm{XB}},\frac{{\mathrm{du}}_{\mathrm{XB}}}{\mathrm{dz}}={-\mathrm{\θ}}_{\mathrm{YB}}.$

(4)假设角位移和位移在相邻两节点之间线性变化，根据图6及以上各式，则各节点的位移为，

对于Y_B方向的位移u_YB ^k：

当节点号k＝1时，其相应的转角与位移为：θ_XB ¹＝0,u_YB ¹＝0；

当节点号k＝2至N时，其相应的转角与位移为，

θ_XB ^k＝θ_XB ^k-1+0.5(κ_XB ^k+κ_XB ^k-1)(z^k-z^k-1)

${u_{\mathrm{YB}}}^k={u_{\mathrm{YB}}}^{k-1}+{{\mathrm{\θ}}_{\mathrm{XB}}}^{k-1}(z^k-z^{k-1})+(\frac{1}{6}{{\mathrm{\κ}}_{\mathrm{XB}}}^k+\frac{1}{3}{{\mathrm{\κ}}_{\mathrm{XB}}}^{k-1}){(z^k-z^{k-1})}^2$

当M_XB、M_YB为极限弯矩时，则由步骤(1)至(4)计算得到的叶尖处位移，也就是u_YB ^N即为叶片在Y_B方向的最大叶尖偏移。

同理，对于X_B方向的最大叶尖偏移位移u_XB ^N可仿照类似方式进行计算。

关于叶片的极限强度及最小极限安全因子的计算

叶片的极限强度主要由截面所受的载荷以及叶片的材料铺层决定。在实际的设计过程中，为了确定叶片是否满足极限强度要求，根据标准，需要进行叶片的静态加载实验。在实验开始前，首先通过对叶片的全工况分析得到叶片各截面的极限载荷，也就是通常说的载荷包络线。然后将某一方向上得到的载荷包络线等效加载到叶片的对应方向上，看此时的叶片是否会发生破坏。如果叶片没被破坏，则说明该叶片的结构是安全可靠的，否则需要改进。由于挥舞方向的载荷最大，刚度较小，所以主要计算叶片在极限挥舞弯矩作用下的材料应变，应力。但同时也将考虑取得最大挥舞弯矩时，叶片截面所受的摆振弯矩，由于截面拉力对材料所受的拉、压应力影响小，在此不予考虑。下面是计算过程：

(a)将叶片坐标系下的弯矩转换到截面的弹性主轴坐标系下，设载荷的安全系数为γ_F，则有，

$\left\{\begin{array}{c}{\mathrm{M\κ}}_{\mathrm{YE}}={\mathrm{\γ}}_F{M}_{\mathrm{YB}}\cos ({\mathrm{\β}}_1+\mathrm{\υ})+{\mathrm{\γ}}_F{M}_{\mathrm{XB}}\sin ({\mathrm{\β}}_1+\mathrm{\υ})\\ M{}_{\mathrm{XE}}=-{\mathrm{\γ}}_F{M}_{\mathrm{YB}}\sin ({\mathrm{\β}}_1+\mathrm{\υ})+{\mathrm{\γ}}_F{M}_{\mathrm{XB}}\cos ({\mathrm{\β}}_1+\mathrm{\υ})\end{array}\right.,$

(b)如图7所示，设弹性主轴坐标系的原点在弦向坐标系中的坐标值为(x_EC,y_EC)，然后将截面的各铺层坐标转换到弹性主轴坐标系下：

$\left\{\begin{array}{c}{x}_E=(x_c-x_{\mathrm{EC}})\cos \left(\mathrm{\υ}\right)-(y_c-y_{\mathrm{EC}})\sin \left(\mathrm{\υ}\right)\\ y_E=(x_c-x_{\mathrm{EC}})\sin \left(\mathrm{\υ}\right)+(y_c-y_{\mathrm{EC}})\cos \left(\mathrm{\υ}\right)\end{array}\right.,$

(c)求得截面各位置处的铺层对应的应变(以拉为正)为，

$\mathrm{\ϵ}(x_E,y_E)=\frac{M_{\mathrm{XE}}}{{\mathrm{EI}}_F}{y}_E+\frac{M_{\mathrm{YE}}}{{\mathrm{EI}}_E}{x}_E$

其中EI_F、EI_E为截面在弹性主轴坐标系x_Eoy_E下的挥舞刚度、摆振刚度；

(d)由此可求得截面各位置处相应的应力(以拉为正)为，

$\mathrm{\σ}(x_E,y_E)={\stackrel{\‾}{Q}}_{11}\mathrm{\ϵ}(x_E,y_E),$

(e)在此基础上，通过进一步计算还可以得到材料的极限安全因子。设截面(x_E,y_E)位置处，材料的许用拉应力为[σ_t]、许用压应力为[σ_c]，材料的安全系数为γ_m，则相应的拉应力安全因子SF_t、压应力安全因子SF_c为

${\mathrm{SF}}_t=\frac{\left[{\mathrm{\σ}}_t\right]}{{\mathrm{\γ}}_m\mathrm{\σ}(x_E,y_E)},{\mathrm{SF}}_c=\frac{\left[{\mathrm{\σ}}_c\right]}{{\mathrm{\γ}}_m\mathrm{\σ}(x_E,y_E)},$

本文在计算各材料的安全因子时，拉应力安全因子均为正，压应力安全因子均为负。所以，SF_t≥1且SF_c≤-1则说明该叶片是安全可靠地，否则是不安全的。这样我们就能更直观的判断出材料的极限强度要求是否满足。

当M_XB、M_YB为极限弯矩时，则由步骤(a)至(d)可计算得到截面各材料的最大应力σ_max(x_E,y_E)；设叶片截面上(x_E,y_E)位置处，材料的许用拉应力为[σ_t]、许用压应力为[σ_c]，材料的安全系数为γ_m，则相应的最小拉应力安全因子SF_{t_min}、最小压应力安全因子SF_{c_min}分别为

${\mathrm{SF}}_{t\_\min }=\frac{\left[{\mathrm{\σ}}_t\right]}{{\mathrm{\γ}}_m\mathrm{\σ}(x_E,y_E)},{\mathrm{SF}}_{c\_\min }=\frac{\left[{\mathrm{\σ}}_c\right]}{{\mathrm{\γ}}_m\mathrm{\σ}(x_E,y_E)}.$

所以，只要SF_{t_min}≥1且SF_{c_min}≤-1则说明该叶片是安全可靠地，否则是不安全的。

进行步骤S6时，沿叶片分段积分得到叶片的总质量。如图4所示，在得到叶片的质量线密度分布后，将叶片划分成多段，通过叶片的分段积分，就可以得到叶片的总质量为其中为截面的质量线密度；m_add为附加质量，如人孔板、楔形条、电缆等附件的质量。但由于这些附加质量在厚度上不均匀、宽度太窄、不易定位等原因，参数化建模复杂，对其一般不考虑它对截面的刚度贡献，而只考虑它的质量。而且由于在叶片的铺层优化过程中，这些附加质量不变，所以铺层质量的变化同样能反映出了整个叶片的质量变化，可以在铺层优化过程中不考虑这些附加质量。

进行步骤S7时，采用如下方法进行优化设计建模：

(1)设计变量

由于梁帽是叶片的主要承力部分且在叶片的总质量中占比较大，所以一般选取梁帽的偶曾参数为设计变量，它包括梁帽铺层的位置、厚度、宽度等参数。在本优化设计模型中，选取梁帽铺层的展向位置和铺层厚度做为设计变量，不考虑梁帽铺层的宽度变化。只是在该模型中，将主梁帽和副梁帽的铺层同时优化，如图8所示，以增加自由变量数，拓宽寻优范围。一般来说，主副梁帽的铺层在压力面和吸力面是一样的，该模型在主梁帽上选取4个位置点，在副梁帽上也选取4个点进行优化建模，由于每个位置点包含有一个铺层位置变量x_i和铺层厚度变量t_i，所以一共有16个设计变量。由于梁帽铺层中，各层材料是一样的，则有t_i＝n_i×t₀,t₀是单层厚度，n_i为铺层层数。由于梁帽各位置处的厚度是各铺层厚度的叠加，所以n_i为正整数。

(2)目标函数

通常来讲，叶片的质量越大，制造所需材料越多，成本越高，同时还会对塔架，旋转轴等其它整机部件的要求提高。所以要求叶片质量越小越好。因此，本模型以叶片质量最小为设计目标。它是设计变量的函数，表示为：Min m(x_i,t_i)(i＝1,…,8)。

(3)约束条件

风力机叶片铺层的优化设计问题是一个复杂的约束优化问题，有许多设计要求需要考虑。本模型中，对以下的设计要求进行约束处理：

第一个是最大叶尖偏移约束条件。根据设计要求，叶片的最大叶尖偏移d_{tip_max}应该低于允许值[d]，以免叶片与塔架发生碰撞。

将各截面极限载荷同时加载到叶片上得到叶片的最大叶尖偏移d_{tip_max}，会导致其大于叶片的实际最大叶尖偏移甚至是允许值[d]，但它并不是一个真实的值。如果对于初始叶片其实际的最大叶尖偏移为d_{tip_act}，而采用上述方法计算得到的最大叶尖偏移为d_{tip_max0}，则当叶片铺层减少时，由于最大叶尖偏移与载荷成正比，采用上述方法计算得到的最大叶尖偏移增量Δd_{tip_max}将大于实际的最大叶尖偏移Δd_{tip_act}，当叶片铺层减少时，如果令Δd_{tip_max}＝[d]-d_{tip_act}，则对应的实际最大叶尖偏移dt_{ip_act}有d_{tip_act}+Δd_{tip_act}≤[d]。所以，本模型在设置最大叶尖偏移约束条件时，一般让d_{tip_max}小于某一设定值d_{tip_set}，该值满足d_{tip_set}≤d_{tip_max0}+[d]-d_{tip_act}。

最后，对采用上述方法计算得到的最大叶尖偏移进行如下约束，

d_{tip_max}≤d_{tip_set}。

第二个是对叶片的一阶挥舞频率进行约束，以避免叶片发生共振。表示如下Flap_1(x_i,t_i)-(Flap_1)_set|≤Δ，上式中，Flap_1(x_i,t_i)是一阶挥舞频率，(Flap_1)_set是设定的一阶挥舞频率，Δ是容许度。

第三个是极限应力条件。为了使叶片安全可靠，材料实际所承受应力σ不能超过材料的许用应力[σ]。σ是考虑了载荷安全系数和材料安全系数的结果。[σ]是材料的许用拉应力或许用压应力。用安全因子表示为，|SF_min|≥1，SF_min为最小极限拉应力安全因子或最小极限压应力安全因子。

考虑到叶片除了受极限应力外，还受疲劳应力的作用。为了保证叶片同时满足疲劳设计要求，叶片的刚度不能降低太多，所以本模型通过提高对极限安全因子的要求，即让它不小于某一设定的安全系数值SF_set，从使优化后的叶片满足疲劳设计要求。用不等式表示为，|SF_min|≥SF_set≥1。

最后，我们还需要对设计变量进行几何约束。它们表示如下，

$\left\{\begin{array}{cc}{x}_{i\_\mathrm{down}}\≤x_i\≤x_{i\_\mathrm{up}}& \\ t_{i\_\mathrm{down}}\≤t_i\≤t_{i\_\mathrm{up}}& (i=1,...,8)\\ t_2=t_3\mathrm{and}{t}_6=t_7& \end{array}\right.$

为了与实际的叶片铺层更相符，梁帽的厚度t_i还必须为单层材料厚度的整数倍，即t_i＝nt₀。

(4)计算流程

由以上的分析可以看到，风力机叶片的铺层优化设计是一个非常复杂的约束优化设计过程。该模型继续采用改进的PSO算法来进行搜索寻优。具体流程图如图9所示。在所述改进的PSO算法中：

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

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

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

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

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

n为粒子总数目的约10％。

由于在该模型中，自由变量、目标函数和约束条件可以根据不同的设计目标和要求进行改变；且该模型能够考虑大部分的设计要求，如质量、频率、最大叶尖偏移、叶片的极限强度要求。可以将其应用于不同叶片的铺层设计。此外，由于在进行最大叶尖偏移、叶片的最小极限安全因子计算时采用了简化方法计算，大大减少了整个优化设计过程需要的时间。

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

Claims (9)

1.一种极限载荷条件下叶片铺层的优化设计方法，其特征在于，所述优化设计方法包括如下步骤：

步骤S1、获得风力机叶片的初始铺层方案，对叶片的初始铺层进行参数化建模；

步骤S2、根据参数化建模后的叶片铺层，得到叶片展向各截面的铺层信息；

步骤S3、根据叶片展向各截面的铺层信息，采用经典层合板理论并假设周向应力为零，计算叶片的质量线密度分布和刚度分布；

步骤S4、根据叶片的质量线密度分布和刚度分布计算叶片的固有频率；

步骤S5、根据叶片各截面的极限载荷分布及叶片的结构特性参数求取叶片的最小极限安全因子以及最大叶尖偏移；所述结构特性参数包括叶片截面的挥舞刚度、摆振刚度、扭转刚度、抗拉刚度、质量线密度、质量中心和弹性中心；

步骤S6、判断叶片的固有频率、最小极限安全因子以及最大叶尖偏移是否满足约束条件，如果满足则沿叶片分段积分得到叶片的总质量；否则，直接令叶片总质量为一很小的常数值；

步骤S7、采用智能优化算法对叶片铺层进行优化，改变叶片的铺层，继续步骤S2至步骤S6，直至得到最优的叶片铺层方案。

2.根据权利要求1所述的极限载荷条件下叶片铺层的优化设计方法，其特征在于，步骤S3中，采用经典层合板理论并假设周向应力为零计算叶片各截面的刚度特性、质量线密度，得到叶片沿展向的质量线密度分布和刚度分布。

3.根据权利要求1所述的极限载荷条件下叶片铺层的优化设计方法，其特征在于，步骤S4中，将叶片简化成叶根固定在轮毂的悬臂梁，根据叶片截面的质量线密度和刚度值得到质量矩阵和刚度矩阵

${\stackrel{\_}{M}}_{\mathrm{ij}}={\∫}_0^L\stackrel{.}{m}\left(x\right){\mathrm{\ψ}}_i\left(x\right){\mathrm{\ψ}}_j\left(x\right)\mathrm{dx},$

${\stackrel{\_}{K}}_{\mathrm{ij}}={\∫}_0^L{\mathrm{EI}}_F\left(x\right){\mathrm{\ψ}}_i^{\′\′}\left(x\right){\mathrm{\ψ}}_j^{\′\′}\left(x\right)\mathrm{dx}+{\∫}_0^L{\mathrm{\ψ}}_i^{\′}\left(x\right){\mathrm{\ψ}}_j^{\′}\left(x\right){\∫}_x^L\left(\stackrel{.}{m}\left(x\right){\mathrm{\Ω}}^{2}r\right)\mathrm{drdx},$

其中，为截面质量线密度；EI_F(x)为弹性主轴坐标系下的截面挥舞刚度；L为叶片单元段的长度；Ω为风轮转速；r为截面的旋转半径；ψ_i(x)和ψ_j(x)为形状函数，分别取为xⁱ⁺¹和x^j+1，式中1≤i,j≤n,n为大于等于5的整数；

计算由质量矩阵和刚度矩阵构成的常系数方程的特征值，得到叶片的固有频率。

4.根据权利要求1所述的极限载荷条件下叶片铺层的优化设计方法，其特征在于，步骤S5中，计算叶片各截面的最大载荷，将各截面的最大载荷同时加载到叶片上，通过静力学分析得到叶片的最大叶尖偏移及截面各材料承受的最大应力σ_max(x_E,y_E)；设叶片截面上(x_E,y_E)位置处，材料的许用拉应力为[σ_t]、许用压应力为[σ_c]，材料的安全系数为γ_m，则相应的最小拉应力安全因子SF_{t_min}、最小压应力安全因子SF_{c_min}分别为 ${\mathrm{SF}}_{t\_\min }=\frac{\left[{\mathrm{\σ}}_t\right]}{{\mathrm{\γ}}_m\mathrm{\σ}(x_E,y_E)},$ ${\mathrm{SF}}_{c\_\min }=\frac{\left[{\mathrm{\σ}}_c\right]}{{\mathrm{\γ}}_m\mathrm{\σ}(x_E,y_E)}.$

5.根据权利要求1所述的极限载荷条件下叶片铺层的优化设计方法，其特征在于，步骤S7中，采用智能优化算法对叶片铺层进行优化时，目标函数根据不同的叶片设计目标和要求进行设定，为叶片质量、叶片一阶挥舞频率、叶片最大叶尖偏移或最小极限安全因子。

6.根据权利要求1所述的极限载荷条件下叶片铺层的优化设计方法，其特征在于，所述智能优化算法为PSO算法或遗传算法。

7.根据权利要求6所述的极限载荷条件下叶片铺层的优化设计方法，其特征在于，所述PSO算法为改进的PSO算法，在所述改进的PSO算法中：

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

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

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

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

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

8.根据权利要求7所述的极限载荷条件下叶片铺层的优化设计方法，其特征在于，n为粒子总数目的约10％。

9.根据权利要求1所述的极限载荷条件下叶片铺层的优化设计方法，其特征在于，步骤S6中，所述很小的常数值为-99999。