CN102708266A - 一种水平轴风力机叶片的极限载荷预测计算方法 - Google Patents
一种水平轴风力机叶片的极限载荷预测计算方法 Download PDFInfo
- Publication number
- CN102708266A CN102708266A CN2012101938161A CN201210193816A CN102708266A CN 102708266 A CN102708266 A CN 102708266A CN 2012101938161 A CN2012101938161 A CN 2012101938161A CN 201210193816 A CN201210193816 A CN 201210193816A CN 102708266 A CN102708266 A CN 102708266A
- Authority
- CN
- China
- Prior art keywords
- foline
- wind
- blade
- section
- angle
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Images
Classifications
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
- Y02E10/00—Energy generation through renewable energy sources
- Y02E10/70—Wind energy
- Y02E10/72—Wind turbines with rotation axis in wind direction
Landscapes
- Wind Motors (AREA)
Abstract
本发明涉及一种水平轴风力机叶片的极限载荷预测计算方法,采用智能优化算法进行极限载荷求解,选取风力机的转速、桨矩角、来流风速、偏航角和方位角为自由变量,建立叶片各截面各方向的载荷同上述自由变量的关系,根据风场类型和设计需要,对各自由变量进行约束,以叶片截面上的载荷为目标函数,使用智能优化算法(比如PSO算法)来求解极限载荷。本发明的方法具有快速而又准确的特点,且计算结果可以方便地应用于叶片铺层优化设计等问题中。
Description
技术领域
本发明涉及风力机叶片的载荷预测计算,特别是水平轴风力机极限载荷的预测计算方法。
技术背景
风力机的载荷源主要有以下几种:空气动力载荷、重力载荷、惯性载荷(包括离心和回转效应)、由于控制系统作用引起的运行载荷(例如,刹车、偏航、变桨距控制和发电机脱网等)。在叶片设计中,主要是考虑由这些载荷源引起的极限载荷和疲劳载荷。到目前为止,大部分风力机主要是因为各种极限状况的出现而失效的,严重的甚至无法修复。因此,极限载荷是叶片铺层设计时的重点。另外,国内的低风况区域分布广泛,随着我国风电产业的发展,设计出能适应低风况的叶片将是大势所趋。但此类叶片的设计与传统的叶片设计相比具有如下特点,叶片本身将变得更长,叶片也更柔,而与整机的接口参数如叶根的最大弯矩等又不允许增加,因此如何控制极限载荷以及开展极限载荷条件下的叶片设计已经成为急需解决的技术难题。而传统的叶片设计理念和方法,大部分均以叶片的最大气动效率,如最大能量输出、最小能量成本,最小质量等经济性指标为设计目标,进行叶片气动、结构设计及优化。但该设计思路无法解决长叶片带来的高载荷问题。所以建立叶片极限载荷的预测计算模型对解决以上问题具有重要意义。
极限载荷预测计算一直是风力机优化设计过程中的一大难点,也是近几年来国外一直研究的热点。目前来说,主要有两类方法来计算预测叶片的极限载 荷,一类是工程上经常使用的方法,也是相对来说最准确的方法,它主要是依据GL或IEC标准定义的工况进行计算分析,对每一个工况都进行气弹模拟,然后根据计算得到的结果进行统计得到。而由于风力机的工作风速范围广,运行工况复杂多变,采用这种方法预测计算出叶片所受到的极限载荷将非常复杂耗时且不易用于叶片铺层的优化设计中。另一类方法主要是基于有限的载荷数据(测量或计算得到的)采用统计学分析方法和外插方法得到叶片的极限载荷,但它们存在以下问题,即如何为已知数据找到一个合适的概率分布函数;如何由短期的载荷分布得到长期的载荷分布;如何定义预测结果的不确定度。所以有必要建立一种新的极限载荷预测模型,它将具有快速而又准确的特点,为叶片的铺层设计甚至是气动设计提供帮助。
发明内容
针对现有技术的缺点和不足,本发明的目的在于提供一种水平轴风力机叶片的极限载荷预测计算方法,该方法具有快速而又准确的特点,且计算结果可以方便地应用于叶片铺层优化设计等问题中。
本发明为解决其技术问题所采取的技术方案为:一种水平轴风力机叶片的极限载荷预测计算方法,采用智能优化算法进行极限载荷求解,其特征在于,所述极限载荷预测计算方法包括如下步骤:
1)选取风力机的转速Ω、桨矩角β2、来流风速V1、偏航角γ和方位角ψ为自由变量;
2)建立叶片各截面各方向的载荷同上述自由变量的关系:
(a)首先建立叶片的气动力Fa和上述自由变量的关系,将叶片分成多个互不相关的叶素,设每个叶素中各截面翼型、来流速度V1、攻角α相同,叶素处的合成气流速度V0作用在长度为dr的叶素上的气动力dFa可分解为法向力dFn和切向力dFt,dFn和dFt可分别表示为,
其中,ρ为空气密度,c为叶素剖面弦长,Cn、Ct分别表示叶素的法向力系数和切向力系数,其表达式为:
其中,Cl、Cd分别为叶素的升力系数和阻力系数,φ为叶素处的入流角;叶素处的入流角φ和合成气流速度V0可分别表示为:
其中,V1为叶素处的来流风速,r为叶素的旋转半径,a为叶素轴向诱导因子,b为叶素切向诱导因子;
3)将离心载荷Fc和重力载荷Fg沿风轮坐标系进行分解,同气动载荷Fa一起作用在叶片上,计算叶片各截面各方向上的载荷;
4)根据风场类型和设计需要,对各自由变量进行如下约束,
其中,Ωupper、Ωlower为转速Ω的上、下限,β2_upper、β2_lower为桨矩角β2的上、下限,Vlower、Vupper为来流风速V1的上、下限,γupper、γlower为偏航角γ的上、下限,ψupper、ψlower为方位角ψ的上、下限,这些参数由风力机实际运行的风场类型、整机参数和设计要求确定;
5)以叶片某一截面的某一方向的载荷、某一截面的不同方向的载荷、叶片多个截面的某一方向的载荷、和/或多个截面的多个方向的载荷为目标函数,以步骤4)设定的约束条件,使用智能优化算法来求解相应载荷的极限值,得到极限载荷。
进一步地,由于控制规律的作用,对不同风速和转速条件下所对应的最小桨距角β2_min按下式来取值:
其中,β2_min_v为不同风速条件下对应的最小桨距角,β2_min_rs为不同转速条件下对应的最小桨距角,Vout为风力机的切出风速;β2_min_v、β2_min_rs按照如下方式进行取值:
(1)当风速Vin<V1<V′时,不同风速条件下对应的最小桨距角β2_min_v为β2_min_v=β2_lower,其中β2_lower为叶片初始安装角,Vin为风力机的切入风速,V′为风力机开始变桨的风速;
(2)当风速V′<V1<Vout时,β2_min_v与风速的对应关系采用以下二次多项式表示:
β2_min_v=B1(V1)2+B2V1+B3-Δβ2,其中Δβ2为常数,且50<Δβ2<150,B1、B2、B3为二次多项式的系数,B1、B2、B3为二次多项式的系数,B1、B2、B3通过对风力机正常运行时不同风速时对应的桨矩角进行二次多项式拟合的方式得到。在实际操作时,可以采用如下方式确定:通过采用FOCUS软件、Bladed等其它气弹软件对风力机进行气弹模拟,我们可以得到风力机正常运行时不同风速时对应的桨矩角,将它们用(ai,bi)表示,i必须大于等于3,ai表示风速值,bi表示对应的桨矩角。ai应取V′<V1<Vout中所有模拟的风速点。根据这些点我们可以在EXCEL表格中作图,并采用回归分 析对这些数据进行二次多项式的拟合,当得到的二次多项式方程在V′<V1<Vout范围内,其R平方值大于等于0.99时,我们取此时拟合得到的二次多项式方程中的二次方系数为B1的值,一次方系数为B2的值,常数项为B3的值。如果采用回归分析得到的二次多项式,其在V′<V1<Vout范围内,R平方值不能大于0.99,则将靠近风速V′和Vout附近的数据点逐个删除,直到满足R平方值大于等于0.99为止。
如果计算得到的β2_min_v≤β2_lower时,取β2_min_v=β2_lower;
(3)当风速Vout<V1<Vupper时,β2_min_v与风速的对应关系采用线性关系式表示如下:
β2_min_v=D1(V1-Vout)+D2-Δβ2,其中D1、D2为线性关系式的系数,采用如下方式确定D1、D2:该直线通过点(Vupper,90-Δβ2)和点(Vout,β2_out),其中β2_out满足β2_out=B1(Vout)2+B2Vout+B3-Δβ2。
(4)当Ω≥Ω′时,Ω′为风力机额定转速附近一设定转速,β2_min_rs与转速Ω的对应关系采用线性关系式表示为:
β2_min_rs=C1Ω+C2,其中C1、C2为线性关系式的系数,通过对风力机运行在额定风速时的一年一遇的极端操作阵风工况进行气弹模拟,得到风力机运行在不同转速时对应的桨矩角,对转速和桨矩角进行线性拟合得到C1、C2。在实际操作时,采用如下方式确定:通过采用FOCUS软件、Bladed等其它气弹软件对风力机运行在额定风速时的一年一遇的极端操作阵风工况(即EOG1)进行气弹模拟,我们可以得到风力机运行在不同转速时对应的桨矩角,将它们用(mi,ni)表示,i必须大于等于2,mi表示转速值,ni表示对应的桨矩角。Ω′取为模拟过程中取得最小桨矩角值时对应的转速。设模拟过程中取得的最大转速为Ω″。则mi应取Ω′≤Ω≤Ω″之 间所有模拟的转速点。根据这些点我们可以在EXCEL表格中作图,并采用回归分析对这些数据进行线性关系式的拟合,当得到的线性关系式在Ω≥Ω′范围内,其R平方值大于等于0.99时,我们取此时拟合得到的线性关系式中的一次方系数为C1的值,常数项为C2的值。如果采用回归分析得到的线性关系式,其在Ω≥Ω′范围内,R平方值不能大于0.99,则将靠近转速Ω′和Ω″附近的数据点逐个删除,直到满足R平方值大于等于0.99为止。
(5)当转速Ω<Ω′时,β2_min_rs=β2_lower。
优选的,采用如下方法计算叶素诱导因子a、b,包括如下计算步骤:
i)设定诱导因子a、b初值:a=a0,b=b0;
j)计算叶素的切向速度Vy0和法向速度Vx0:Vx0=V1(1-a),Vy0=Ωr(1+b);
k)计算叶素处的入流角φ和攻角α:
其中,β1是叶素截面的扭角,β2是叶素的桨矩角;
l)计算叶片损失F,F=FtipFhub,叶尖损失Ftip、轮毂损失Fhub分别表示为:
其中,B是叶片数,Rhub为轮毂半径,R是风轮的旋转半径;
n)求解新的切向诱导因子a:
如果CT≥0.96F,则该叶素载荷过高,新的轴向诱导因子a由下式求解:
其中,B2=1/0.18-4F,B1=-(0.8/0.18-4F),B0=0.16/0.18;
否则,新的轴向诱导因子a为
o)求解新的切向诱导因子b,
p)如果前后两次计算出来的诱导因子a、b变化大于某一容许偏差,则返回
步骤b)循环运算,否则完成计算。
优选的,采用如下方法计算叶素的初始诱导因子a0,b0,包括如下计算步骤:
a)计算除叶根处的圆截面叶素外,其他所有叶素的入流角φ和合成气流速度V0,
其中,
θ=β1+β2,为叶素截面弦线与风轮平面的夹角,
λr为叶素局部尖速比,σ为风轮实度,kl为升力线斜率,Cl0为零攻角对应的升力系数,β1是叶素截面的扭角,β2是叶素的桨矩角;
b)根据以上得到的入流角φ和合成气流速度V0,计算除叶根处的圆截面叶素外的初始诱导因子a0,b0,
c)计算叶根处的圆截面叶素的入流角φ和合成气流速度V0,
φ=arctan(1/λr),
其中,a=σ/(4sinφ+σ);
d)根据以上得到的入流角φ和合成气流速度V0,计算叶根处的圆截面叶素的初始诱导因子a0,b0,
a0=σ/(4sinφ+σ),b0=-a0
优选的,在计算叶素上的气动力dFa时,也可采用如下方法计算叶片各截面各叶素截面的入流角φ和合成气流速度V0,包括如下步骤:
(1)当局部尖速比λr≥4时,根据权利要求4所述的步骤a)、c),针对不同形状的叶素截面直接计算叶素截面的入流角φ和合成气流速度V0,
(2)当局部尖速比λr<4时,根据权利要求3计算叶素截面的入流角φ和合成气流速度V0。
优选的,采用如下方法计算叶片各截面各方向上的载荷,包括如下步骤:
(a)将离心载荷Fc和重力载荷Fg沿风轮坐标系进行分解,同气动载荷Fa一起作用在叶片上,将叶片分成N个叶素,N为大于2的整数,每个叶素对应一个节点,假设每一个叶素上的载荷是均匀分布的,各叶素上的均布载荷pk=Fc+Fg+Fa,k为整数,1≤k≤N,设x方向为风轮旋转方向,y方向为垂直于 风轮旋转平面方向,z方向为旋转中心指向叶尖为正(不考虑预弯);1号节点和N号节点对应的叶素长度dz1、dzN为零:dzN=dz1=0;
(b)计算各叶素截面x方向上的力Tx k和由此产生的力矩My k:
对于节点1:Tx 1=Tx 2+0.5p2dz2,My 1=My 2+0.5Tx 2dz2+p2(dz2)2/8;
对于节点N:Tx N=0,My N=0;
对于节点N-1:Tx N-1=Tx N+0.5pN-1dzN-1,My N-1=My N+0.5Tx NdzN-1+pN-1(dzN-1)2/8;
对于节点N-2到2:
Tx k=Tx k+1+0.5pk+1dzk+1+0.5pkdzk,
My k=My k+1+0.5Tx k+1(dzk+1+dzk)+0.5pk+1dzk+1(0.5dzk+0.25dzk+1)+pk(dzk)2/8;
(c)按照和步骤(b)类似的方法计算各叶素截面y方向上的力Ty k和由此产生的力矩Mx k。
优选的,所述智能优化算法为PSO算法(粒子群优化算法)。
优选的,所述PSO算法为改进的PSO算法,在所述改进的PSO算法中:
(a)惯性权值w按对数规律单调递减,其表达式为:
其中,P为当前迭代次数,gen为最大迭代次数,wmax为最大惯性权值因子,wmin为最小惯性权值因子;
(b)将每一代最优的前n个粒子直接复制到下一代,使每一代得到的最优解得以保存而不被破坏,其余的粒子正常进化后进入下一代;
(c)将每一代中所有不满足约束条件的粒子参与进化生成下一代,但令这些粒子的目标函数值同为一个很小的常数,其余的粒子正常进化后进入下一代。
优选的,n为粒子总数目的约10%。
本发明的水平轴风力机叶片的极限载荷预测计算方法,主要特点在于:
1)将叶片的极限载荷计算看成是求极值问题,采用智能优化算法,如PSO算法(粒子群算法)进行极限载荷求解;
2)依据不同的叶片气动外形参数、风场类型、控制特性等参数进行载荷预测计算;
3)能准确快速的得到风力机叶片各截面的极限载荷值或是同时得到多个载荷分量的极限值,为叶片的铺层设计提供依据。
附图说明
图1为叶素上的气流速度三角形和空气动力分量示意图。
图2为叶片质量计算示意图。
图3为叶片某一方向的载荷分布图。
图4为叶片截面载荷计算流程图。
图5为极限载荷预测计算流程图。
具体实施方式
下面结合附图和具体实施例,进一步阐明本发明,应理解这些实施例仅用于说明本发明而不用于限制本发明的范围,在阅读了本发明之后,本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。
由于风力机的极限载荷预测也是一种求取极值问题,可以采用优化算法进行求解。所以我们将改进的PSO算法以及改进的多目标PSO算法用于风力机极限载荷的求解,建立了一种新的极限载荷预测模型。以下是整个极限载荷预测模型的建模过程。
(一)自由变量选取
对于风力机叶片所受的气动力F和气动力矩M,可以用以下函数表示。
(F,M)=f(γ,χ,t,β2,Ω,β1,c,V1,ψ,ρ,Cl,Cd)
式中,γ为偏航角,χ为锥角,t为倾角,β2是桨矩角,Ω为风轮转速,β1为叶片的截面扭角,c为叶片的截面弦长,V1为来流风速(包括其大小和方向),ψ为方位角,ρ为空气密度,Cl和Cd分别为该截面翼型的升、阻力系数。
由于叶片的气动载荷主要同风力机的转速Ω、桨矩角β2、来流风速V1、偏航角γ、方位角ψ等相关,所以选取转速、桨矩角、来流风速、偏航角和方位角为自由变量。
(二)载荷计算方法
1.气动载荷计算
首先建立叶片的气动力Fa和上述自由变量的关系,将叶片分成多个互不相关的叶素,设每个叶素中各截面翼型、来流速度V1、攻角α相同,图1为叶素上的气流速度三角形和空气动力分量示意图。叶素处的合成气流速度V0作用在长度为dr的叶素上的气动力dFa可分解为法向力dFn和切向力dFt,dFn和dFt可分别表示为,
其中,ρ为空气密度,c为叶素剖面弦长,Cn、Ct分别表示叶素的法向力系数和切向力系数,其表达式为:
其中,Cl、Cd分别为叶素的升力系数和阻力系数,φ为叶素处的入流角;叶素处的入流角φ和合成气流速度V0可分别表示为:
其中,V1为叶素处的来流风速,r为叶素的旋转半径,a为叶素轴向诱导因子,b为叶素切向诱导因子。
采用如下方法计算叶素诱导因子a、b,包括如下计算步骤:
q)设定诱导因子a、b初值:a=a0,b=b0;
r)计算叶素的切向速度Vy0和法向速度Vx0:Vx0=V1(1-a),Vy0=Ωr(1+b);
s)计算叶素处的入流角φ和攻角α:
其中,β1是叶素截面的扭角,β2是叶素的桨矩角;
t)计算叶片损失F,F=FtipFhub,叶尖损失Ftip、轮毂损失Fhub分别表示为:
其中,B是叶片数,Rhub为轮毂半径,R是风轮的旋转半径;
v)求解新的切向诱导因子a:
如果CT≥0.96F,则该叶素载荷过高,新的轴向诱导因子a由下式求解:
其中,B2=1/0.18-4F,B1=-(0.8/0.18-4F),B0=0.16/0.18;
否则,新的轴向诱导因子a为
w)求解新的切向诱导因子b,
x)如果前后两次计算出来的诱导因子a、b变化大于某一容许偏差,则返回步骤b)循环运算,否则完成计算。
采用如下方法计算初始诱导因子a0,b0,包括如下计算步骤:
a)计算除叶根处的圆截面叶素外,其他所有叶素的入流角φ和合成气流速度V0,
其中,
θ=β1+β2,为叶素截面弦线与风轮平面的夹角,
λr为叶素局部尖速比,σ为风轮实度,kl为升力线斜率,Cl0为零攻角对应的升力系数,β1是叶素截面的扭角,β2是叶素的桨矩角;
b)根据以上得到的入流角φ和合成气流速度V0,计算除叶根处的圆截面叶素外的初始诱导因子a0,b0,
c)计算叶根处的圆截面叶素的入流角φ和合成气流速度V0,
φ=arctan(1/λr),
其中,a=σ/(4sinφ+σ);
d)根据以上得到的入流角φ和合成气流速度V0,计算叶根处的圆截面叶素 的初始诱导因子a0,b0,
a0=σ/(4sinφ+σ),b0=-a0
2.离心载荷和重力载荷的计算
将离心载荷和重力载荷沿风轮坐标系进行分解,同气动载荷一起作用在叶片上。假设每一个叶素上的载荷是均匀分布的,如图3所示,x方向为风轮旋转方向,y方向为垂直于风轮旋转平面方向,它是垂直于纸面向里,z方向为旋转中心指向叶尖为正(不考虑预弯)。在图3中,均布载荷pk为pk=Fc+Fg+Fa。1号节点和N号节点的对应的叶素长度为零,即dzN=dz1=0。然后在此基础上计算各截面的力Tx k和由此产生的力矩My k:
对于节点1:Tx 1=Tx 2+0.5p2dz2,My 1=My 2+0.5Tx 2dz2+p2(dz2)2/8;
对于节点N:Tx N=0,My N=0;
对于节点N-1:Tx N-1=Tx N+0.5pN-1dzN-1,MyN-1=My N+0.5Tx NdzN-1+pN-1(dzN-1)2/8;
对于节点N-2到2:
Tx k=Tx k+1+0.5pk+1dzk+1+0.5pkdzk,
My k=My k+1+0.5Tx k+1(dzk+1+dzk)+0.5pk+1dzk+1(0.5dzk+0.25dzk+1)+pk(dzk)2/8。
对于y方向,其计算方法类似。
综合上面有关气动载荷、离心载荷和重力载荷的计算方法,得到叶片截面载荷计算流程如图4所示。
(三)约束条件处理
以某一1.5MW风力机为计算对象,根据风场类型、整机参数和设计需要,对各自由变量进行如下约束,
另外,由于控制规律的作用,对不同风速和转速条件下所对应的桨距角取值范围进行如下约束:
β2≥β2_min,
其中β2_min是不同风速和转速条件下所对应的最小桨距角,它采用如下步骤确定:
(1)当风速3≤V1≤11m/s时,最小桨矩角与风速的对应关系采用下式表示
β2_min_v=β2_lower=-1,
(2)当风速11<V1≤25m/s时,最小桨矩角与风速的对应关系采用以下二次多项式表示。
β2_min_v=-0.0579(V1)2+3.5789V1-30.611-Δβ2,
其中Δβ2=10.86
由于桨矩角始终不能小于β2_lower,所以有,
β2_min_v=β2_lower(β2_min_v<β2_lower),
(3)当风速25<V1≤38.26m/s时,将最小桨矩角与风速的对应关系表示如下,
β2_min_v=5.0908(V1-25)+22.674-Δβ2,
此外,桨矩角总是根据转矩来控制的,而转矩与转速直接相关,所以当风轮的转速增加时,其对应的最小桨矩角也会增加。
(4)当Ω≥17.4时,最小桨矩角与转速的对应关系表示如下,
β2_min_rs=1.9969Ω-35.681,
本模型中,根据额定风速时的一年一遇的极端操作阵风(EOG)工况确定不同转速时,转速与其对应的最小桨矩角之间的关系并拟合得到上式的线性表达式。
(5)当Ω<Ω′时,最小桨矩角与转速的对应关系表示如下,
β2_min_rs=β2_lower=-1
(6)综合以上各式,得到不同风速和转速条件下所对应的最小桨矩角β2_min为,
需要注意的是,由于风力机以及所适用的风场类型不一样,约束条件的具体数值会有所不同,但都采用以上的方法进行确定。
(四)极限载荷计算流程
由于叶片各截面各方向的极限载荷不一定同时出现,所以需要对某一截面的某一方向的极限载荷进行单独求解。如在计算风力机的载荷包络线时,常选取某一截面的弦向载荷分量进行求解。有时为了实际需要,也要求对不同截面的不同方向的极限载荷进行同时求解。如在校核分析截面的极限强度是否满足设计要求时,需要对截面内两个互相垂直方向的载荷分量进行同时求解。可以采用各种智能优化算法来进行叶片极限载荷的求解,比如PSO算法,也可以对PSO算法进行改进,采用改进的PSO算法。
这里介绍一种改进的PSO算法,在保留传统PSO算法的主要特点外,对其进行如下改进:
(a)惯性权值w按对数规律单调递减,其表达式为:
其中,P为当前迭代次数,gen为最大迭代次数,wmax为最大惯性权值因子,wmin为最小惯性权值因子;
(b)将每一代最优的前n个粒子直接复制到下一代,使每一代得到的最优解得以保存而不被破坏,其余的粒子正常进化后进入下一代;
(c)将每一代中所有不满足约束条件的粒子参与进化生成下一代,但令这些粒子的目标函数值同为一个很小的常数,其余的粒子正常进化后进入下一代。
最后将以上模型同改进的PSO算法一起结合,建立一种新的风力机叶片极限载荷预测模型。其计算流程如图5所示。对于多个载荷分量同时求解的情况,与此类似,只是在得到目标函数值以后增加了pareto解集的筛选过程。
以上所述仅是本发明的优选实施方式,应当指出,对于本技术领域的普通技术人员来说,在不脱离本发明原理的前提下,还可以做出若干改进和润饰,这些改进和润饰也应视为本发明的保护范围。
Claims (10)
1.一种水平轴风力机叶片的极限载荷预测计算方法,采用智能优化算法进行极限载荷求解,其特征在于,所述极限载荷预测计算方法包括如下步骤:
1)选取风力机的转速Ω、桨矩角β2、来流风速V1、偏航角γ和方位角ψ为自由变量;
2)建立叶片各截面各方向的载荷同上述多个自由变量的关系:
(a)建立叶片的气动力Fa和上述多个自由变量的关系,将叶片分成多个互不相关的叶素,设每个叶素中各截面翼型、来流速度V1、攻角α相同,叶素处的合成气流速度V0作用在长度为dr的叶素上的气动力dFa可分解为法向力dFn和切向力dFt,dFn和dFt可分别表示为,
其中,ρ为空气密度,c为叶素剖面弦长,Cn、Ct分别表示叶素的法向力系数和切向力系数,Cn、Ct的表达式为:
其中,Cl、Cd分别为叶素的升力系数和阻力系数,φ为叶素截面处的入流角;叶素截面处的入流角φ和合成气流速度V0可分别表示为:
其中,V1为叶素截面处的来流风速,r为叶素截面的旋转半径,a为叶素轴向诱导因子,b为叶素切向诱导因子;
3)将离心载荷Fc和重力载荷Fg沿风轮坐标系进行分解,同气动载荷Fa一起作用在叶片上,计算叶片各截面各方向上的载荷;
4)根据风场类型和设计需要,对各自由变量进行如下约束,
其中,Ωupper、Ωlower为转速Ω的上、下限,β2_upper为桨矩角β2的上限、β2_min为最小桨矩角,Vlower、Vupper为来流风速V1的上、下限,γupper、γlower为偏航角γ的上、下限,ψupper、ψlower为方位角ψ的上、下限,这些参数由风力机实际运行的风场类型、整机参数和设计要求确定;
5)以叶片某一截面的某一方向的载荷、某一截面的不同方向的载荷、叶片多个截面的某一方向的载荷、和/或多个截面的多个方向的载荷为目标函数,以步骤4)设定的约束条件,使用智能优化算法来求解相应载荷的极限值,得到极限载荷。
2.根据权利要求1所述的水平轴风力机叶片的极限载荷预测计算方法,其特征在于,在上述步骤4)中最小桨距角β2_min按下式来取值:
其中,β2_min_v为不同风速条件下对应的最小桨距角,β2_min_rs为不同转速条件下对应的最小桨距角,Vout为风力机的切出风速;β2_min_v、β2_min_rs按照如下方式进行取值:(1)当风速Vin<V1<V′时,不同风速条件下对应的最小桨距角β2_min_v为β2_min_v=β2_lower,其中β2_lower为叶片初始安装角,Vin为风力机的切入风速,V′为风力机开始变桨的风速;
(2)当风速V′<V1<Vout时,β2_min_v与风速V1的对应关系采用以下二次多项式表示:
β2_min_v=B1(V1)2+B2V1+B3-Δβ2,其中Δβ2为常数,且50<Δβ2<150,B1、B2、B3为二次多项式的系数,B1、B2、B3通过对风力机正常运行时不同风速对应的桨矩角进行二次多项式拟合的方式得到;
如果计算得到的β2_min_v≤β2_lower时,取β2_min_v=β2_lower;
(3)当风速Vout≤V1<Vupper时,β2_min_v与风速V1的对应关系采用线性关系式表示如下:
β2_min_v=D1(V1-Vout)+D2-Δβ2,其中D1、D2为线性关系式的系数,D1、D2采用如下方式确定:该直线通过点(Vupper,90-Δβ2)和点(Vout,β2_out),其中β2_out满足β2_out=B1(Vout)2+B2Vout+B3-Δβ2;
(4)当Ω≥Ω′时,Ω′为风力机额定转速附近一设定转速,β2_min_rs与转速Ω的对应关系采用线性关系式表示为:
β2_min_rs=C1Ω+C2,其中为线性关系式的系数,通过对风力机运行在额定风速时的一年一遇的极端操作阵风工况进行气弹模拟,得到风力机运行在不同转速时对应的桨矩角,对转速和桨矩角进行线性拟合得到C1、C2;
(5)当转速Ω<Ω′时,β2_min_rs=β2_lower。
3.根据权利要求1所述的水平轴风力机叶片的极限载荷预测计算方法,其特征在于,采用如下方法计算叶素诱导因子a、b:
a)设定诱导因子a、b初值:a=a0,b=b0;
b)计算叶素的切向速度Vy0和法向速度Vx0:Vx0=V1(1-a),Vy0=Ωr(1+b);
c)计算叶素截面处的入流角φ和攻角α:
其中,β1是叶素截面的扭角,β2是叶素的桨矩角;
d)计算叶片损失F,F=FtipFhub,叶尖损失Ftip、轮毂损失Fhub分别表示为:
其中,B是叶片数,Rhub为轮毂半径,R是风轮的旋转半径;
f)求解新的切向诱导因子a:
如果CT≥0.96F,则该叶素载荷过高,新的轴向诱导因子a由下式求解:
其中,B2=1/0.18-4F,B1=-(0.8/0.18-4F),B0=0.16/0.18;
h)如果前后两次计算出来的诱导因子a、b变化大于某一容许偏差,则返回步骤b)循环运算,否则完成计算。
4.根据权利要求3所述的水平轴风力机叶片的极限载荷预测计算方法,其特征在于,采用如下方法计算叶素的初始诱导因子a0,b0,包括如下计算步骤:
a)计算除叶根处的圆截面叶素外,其他所有叶素截面的入流角φ和合成气流速度V0,
其中,
为叶尖损失,
θ=β1+β2,为叶素截面弦线与风轮平面的夹角,
λr为叶素局部尖速比,σ为风轮实度,kl为升力线斜率,Cl0为零攻角对应的升力系数,β1是叶素截面的扭角,β2是叶素的桨矩角;
b)根据步骤a)得到的入流角φ和合成气流速度V0,计算除叶根处的圆截面叶素外的初始诱导因子a0,b0,
c)计算叶根处的圆截面叶素的入流角φ和合成气流速度V0,
其中,a=σ/(4sinφ+σ);
d)根据步骤c)得到的入流角φ和合成气流速度V0,计算叶根处的圆截面叶素的初始诱导因子a0,b0,
a0=σ/(4sinφ+σ),b0=-a0
5.根据权利要求4所述的水平轴风力机叶片的极限载荷预测计算方法,其特征在于,当局部尖速比λr≥4时,根据权利要求4所述的步骤a)、c)计算叶素截面的入流角φ和合成气流速度V0。
7.根据上述任一项权利要求所述的水平轴风力机叶片的极限载荷预测计算方法,其特征在于,采用如下方法计算叶片各截面各方向上的载荷,包括如下步骤:
(a)将离心载荷Fc和重力载荷Fg沿风轮坐标系进行分解,同气动载荷Fa一起作用在叶片上,将叶片分成N个叶素,N为大于2的整数,每个叶素对应一个节点,假设每一个叶素上的载荷是均匀分布的,各叶素上的均布载荷pk=Fc+Fg+Fa,k为整数,1≤k≤N,设x方向为风轮旋转方向,y方向为垂直于风轮旋转平面方向,z方向为旋转中心指向叶尖为正;1号节点和N号节点对应的叶素长度dz1、dzN为零:dzN=dz1=0;
(b)计算各叶素截面x方向上的力Tx k和由此产生的力矩My k:
对于节点1:Tx 1=Tx 2+0.5p2dz2,My 1=My 2+0.5Tx 2dz2+p2(dz2)2/8;
对于节点N:Tx N=0,My N=0;
对于节点N-1:Tx N-1=Tx N+0.5pN-1dzN-1,My N-1=My N+0.5Tx NdzN-1+pN-1(dzN-1)2/8;
对于节点N-2到2:
Tx k=Tx k+1+0.5pk+1dzk+1+0.5pkdzk,
My k=My k+1+0.5Tx k+1(dzk+1+dzk)+0.5pk+1dzk+1(0.5dzk+0.25dzk+1)+pk(dzk)2/8;
(c)按照和步骤(b)类似的方法计算各叶素截面y方向上的力Ty k和由此产生的力矩Mx k。
8.根据上述任一项权利要求所述的水平轴风力机叶片的极限载荷预测计算方法,其特征在于,所述智能优化算法为PSO算法。
10.根据权利要求9所述的水平轴风力机叶片的极限载荷预测计算方法,其特征在于,n为粒子总数目的约10%。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201210193816.1A CN102708266B (zh) | 2012-06-12 | 2012-06-12 | 一种水平轴风力机叶片的极限载荷预测计算方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201210193816.1A CN102708266B (zh) | 2012-06-12 | 2012-06-12 | 一种水平轴风力机叶片的极限载荷预测计算方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN102708266A true CN102708266A (zh) | 2012-10-03 |
CN102708266B CN102708266B (zh) | 2014-01-01 |
Family
ID=46901029
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201210193816.1A Active CN102708266B (zh) | 2012-06-12 | 2012-06-12 | 一种水平轴风力机叶片的极限载荷预测计算方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN102708266B (zh) |
Cited By (17)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103953503A (zh) * | 2014-04-18 | 2014-07-30 | 南车株洲电力机车研究所有限公司 | 风力发电机组偏航制动力矩控制装置及其方法 |
CN105320785A (zh) * | 2014-07-09 | 2016-02-10 | 南车株洲电力机车研究所有限公司 | 一种用于风电机组载荷计算的工况全自动生成方法及装置 |
CN105422391A (zh) * | 2015-12-22 | 2016-03-23 | 南车株洲电力机车研究所有限公司 | 一种风力发电机组极限载荷辨识方法 |
CN106768917A (zh) * | 2016-11-23 | 2017-05-31 | 中国科学院工程热物理研究所 | 一种风力机叶片现场载荷测试与评估方法 |
CN107194028A (zh) * | 2017-04-19 | 2017-09-22 | 中国航空工业集团公司金城南京机电液压工程研究中心 | 一种rat用叶片设计方法 |
CN107203364A (zh) * | 2017-05-26 | 2017-09-26 | 哈尔滨工程大学 | 一种用于压气机全工况特性的预测和辨识方法 |
CN107762728A (zh) * | 2016-08-19 | 2018-03-06 | 北京天诚同创电气有限公司 | 偏航和变桨控制方法、控制系统及风力发电机组 |
CN109598030A (zh) * | 2018-11-14 | 2019-04-09 | 南京航空航天大学 | 一种风力机叶尖损失修正计算方法 |
CN109715937A (zh) * | 2016-09-15 | 2019-05-03 | 乌本产权有限公司 | 用于确定运行载荷和设计塔建筑物的方法、塔建筑物和风能设备 |
CN110023621A (zh) * | 2016-10-17 | 2019-07-16 | 诺迈士科技有限公司 | 确定风力涡轮上的载荷 |
WO2019165753A1 (zh) * | 2018-02-28 | 2019-09-06 | 北京金风科创风电设备有限公司 | 风力发电机组的载荷预测方法和装置 |
CN110501127A (zh) * | 2019-08-28 | 2019-11-26 | 湘潭大学 | 一种基于损伤状态倾角斜率的等截面梁损伤识别方法 |
CN111651841A (zh) * | 2020-05-30 | 2020-09-11 | 扬州大学 | 基于圆周割线改进型粒子群算法的叶片临界颤振系统参数辨识方法 |
US11261846B2 (en) | 2019-11-01 | 2022-03-01 | General Electric Company | System and method for designing and operating a wind turbine power system based on statistical analysis of operational and/or grid data thereof |
CN115045859A (zh) * | 2022-05-30 | 2022-09-13 | 西安交通大学 | 一种离心鼓风机复合叶轮设计方法 |
CN115982897A (zh) * | 2023-03-21 | 2023-04-18 | 浙江华东测绘与工程安全技术有限公司 | 一种海上风机叶片的空气动力载荷等效构建方法及装置 |
CN117744409A (zh) * | 2024-02-19 | 2024-03-22 | 南京航空航天大学 | 海上浮式风机叶片变形与叶轮轮毂载荷预测方法及系统 |
Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20100259046A1 (en) * | 2007-11-06 | 2010-10-14 | Sridhar Kota | Active control surfaces for wind turbine blades |
CN102332043A (zh) * | 2011-09-16 | 2012-01-25 | 中国科学院工程热物理研究所 | 一种基于结构尺寸参数优化的风力机叶片优化设计方法 |
-
2012
- 2012-06-12 CN CN201210193816.1A patent/CN102708266B/zh active Active
Patent Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20100259046A1 (en) * | 2007-11-06 | 2010-10-14 | Sridhar Kota | Active control surfaces for wind turbine blades |
CN102332043A (zh) * | 2011-09-16 | 2012-01-25 | 中国科学院工程热物理研究所 | 一种基于结构尺寸参数优化的风力机叶片优化设计方法 |
Non-Patent Citations (2)
Title |
---|
《工程热物理学报》 20080531 廖猜猜等 "基于一种改进的PSO算法在风力机叶片优化中的应用" 第773-776页 1-10 第29卷, 第5期 * |
廖猜猜等: ""基于一种改进的PSO算法在风力机叶片优化中的应用"", 《工程热物理学报》 * |
Cited By (28)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103953503A (zh) * | 2014-04-18 | 2014-07-30 | 南车株洲电力机车研究所有限公司 | 风力发电机组偏航制动力矩控制装置及其方法 |
CN105320785A (zh) * | 2014-07-09 | 2016-02-10 | 南车株洲电力机车研究所有限公司 | 一种用于风电机组载荷计算的工况全自动生成方法及装置 |
CN105320785B (zh) * | 2014-07-09 | 2018-05-29 | 南车株洲电力机车研究所有限公司 | 一种用于风电机组载荷计算的工况全自动生成方法及装置 |
CN105422391B (zh) * | 2015-12-22 | 2018-02-02 | 南车株洲电力机车研究所有限公司 | 一种风力发电机组极限载荷辨识方法 |
CN105422391A (zh) * | 2015-12-22 | 2016-03-23 | 南车株洲电力机车研究所有限公司 | 一种风力发电机组极限载荷辨识方法 |
CN107762728B (zh) * | 2016-08-19 | 2019-08-16 | 北京天诚同创电气有限公司 | 偏航控制方法、控制系统及风力发电机组 |
CN107762728A (zh) * | 2016-08-19 | 2018-03-06 | 北京天诚同创电气有限公司 | 偏航和变桨控制方法、控制系统及风力发电机组 |
CN109715937A (zh) * | 2016-09-15 | 2019-05-03 | 乌本产权有限公司 | 用于确定运行载荷和设计塔建筑物的方法、塔建筑物和风能设备 |
CN110023621B (zh) * | 2016-10-17 | 2024-01-02 | 诺迈士科技有限公司 | 确定风力涡轮上的载荷 |
CN110023621A (zh) * | 2016-10-17 | 2019-07-16 | 诺迈士科技有限公司 | 确定风力涡轮上的载荷 |
CN106768917A (zh) * | 2016-11-23 | 2017-05-31 | 中国科学院工程热物理研究所 | 一种风力机叶片现场载荷测试与评估方法 |
CN107194028A (zh) * | 2017-04-19 | 2017-09-22 | 中国航空工业集团公司金城南京机电液压工程研究中心 | 一种rat用叶片设计方法 |
CN107194028B (zh) * | 2017-04-19 | 2020-10-27 | 中国航空工业集团公司金城南京机电液压工程研究中心 | 一种rat用叶片设计方法 |
CN107203364A (zh) * | 2017-05-26 | 2017-09-26 | 哈尔滨工程大学 | 一种用于压气机全工况特性的预测和辨识方法 |
CN107203364B (zh) * | 2017-05-26 | 2020-12-22 | 哈尔滨工程大学 | 一种用于压气机全工况特性的预测和辨识方法 |
WO2019165753A1 (zh) * | 2018-02-28 | 2019-09-06 | 北京金风科创风电设备有限公司 | 风力发电机组的载荷预测方法和装置 |
CN109598030B (zh) * | 2018-11-14 | 2019-09-10 | 南京航空航天大学 | 一种风力机叶尖损失修正计算方法 |
CN109598030A (zh) * | 2018-11-14 | 2019-04-09 | 南京航空航天大学 | 一种风力机叶尖损失修正计算方法 |
CN110501127A (zh) * | 2019-08-28 | 2019-11-26 | 湘潭大学 | 一种基于损伤状态倾角斜率的等截面梁损伤识别方法 |
CN110501127B (zh) * | 2019-08-28 | 2021-01-22 | 湘潭大学 | 一种基于损伤状态倾角斜率的等截面梁损伤识别方法 |
CN110501127B9 (zh) * | 2019-08-28 | 2021-04-09 | 湘潭大学 | 一种基于损伤状态倾角斜率的等截面梁损伤识别方法 |
US11261846B2 (en) | 2019-11-01 | 2022-03-01 | General Electric Company | System and method for designing and operating a wind turbine power system based on statistical analysis of operational and/or grid data thereof |
CN111651841A (zh) * | 2020-05-30 | 2020-09-11 | 扬州大学 | 基于圆周割线改进型粒子群算法的叶片临界颤振系统参数辨识方法 |
CN111651841B (zh) * | 2020-05-30 | 2024-01-26 | 扬州大学 | 基于圆周割线改进型粒子群算法的叶片临界颤振系统参数辨识方法 |
CN115045859A (zh) * | 2022-05-30 | 2022-09-13 | 西安交通大学 | 一种离心鼓风机复合叶轮设计方法 |
CN115982897A (zh) * | 2023-03-21 | 2023-04-18 | 浙江华东测绘与工程安全技术有限公司 | 一种海上风机叶片的空气动力载荷等效构建方法及装置 |
CN115982897B (zh) * | 2023-03-21 | 2023-08-15 | 浙江华东测绘与工程安全技术有限公司 | 一种海上风机叶片的空气动力载荷等效构建方法及装置 |
CN117744409A (zh) * | 2024-02-19 | 2024-03-22 | 南京航空航天大学 | 海上浮式风机叶片变形与叶轮轮毂载荷预测方法及系统 |
Also Published As
Publication number | Publication date |
---|---|
CN102708266B (zh) | 2014-01-01 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN102708266B (zh) | 一种水平轴风力机叶片的极限载荷预测计算方法 | |
Xudong et al. | Shape optimization of wind turbine blades | |
Pourrajabian et al. | Aero-structural design and optimization of a small wind turbine blade | |
Wu et al. | Modeling turbine wakes and power losses within a wind farm using LES: An application to the Horns Rev offshore wind farm | |
Chen et al. | Numerical analysis of unsteady aerodynamic performance of floating offshore wind turbine under platform surge and pitch motions | |
Delafin et al. | Effect of the number of blades and solidity on the performance of a vertical axis wind turbine | |
De Tavernier et al. | Airfoil optimisation for vertical‐axis wind turbines with variable pitch | |
Sebastian et al. | A comparison of first-order aerodynamic analysis methods for floating wind turbines | |
Rogowski et al. | Steady and unsteady analysis of NACA 0018 airfoil in vertical-axis wind turbine | |
Chaudhary et al. | Modeling and optimal design of small HAWT blades for analyzing the starting torque behavior | |
Branlard et al. | A multi-purpose lifting-line flow solver for arbitrary wind energy concepts | |
CN111734585A (zh) | 风力发电机的极限载荷的确定方法、装置及可读存储介质 | |
Sakib et al. | Parked and operating loads analysis in the aerodynamic design of multi-megawatt-scale floating vertical axis wind turbines | |
Rogowski et al. | Numerical analysis of a small-size vertical-axis wind turbine performance and averaged flow parameters around the rotor | |
Brown et al. | Rapidly recovering wind turbine wakes with dynamic pitch and rotor speed control | |
John et al. | Helical vortex theory and blade element analysis of multi‐bladed windmills | |
MURATOĞLU et al. | Numerical analyses of a straight bladed vertical axis Darrieus wind turbine: Verification of DMS algorithm and Qblade code | |
Hasan | Design and performance analysis of small scale horizontal axis wind turbine for nano grid application | |
CN109611268A (zh) | 一种双叶轮水平轴风力机设计优化方法 | |
Martin et al. | Site specific optimization of rotor/generator sizing of wind turbines | |
Hosseinkhani et al. | Performance Prediction of a SANDIA 17-m Vertical Axis Wind Turbine Using Improved Double Multiple Streamtube | |
Hong et al. | The design and testing of a small-scale wind turbine fitted to the ventilation fan for a livestock building | |
Koulatsou et al. | Artificial Time Histories of Wind ActionsFor Structural Analysis of Wind Turbines | |
Kesikbaş | Investigations of upscaling effects for aerodynamic design of large wind turbine rotors by using bem theory and optimization | |
Singh et al. | Design and Optimisation of a 20kW Horizontal Axis Wind Turbine using HARP_Opt |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
C14 | Grant of patent or utility model | ||
GR01 | Patent grant |