一种考虑翼型附着流和后缘分离流的动态失速计算方法

1.本发明涉及风力机叶片技术领域。具体地说是一种考虑翼型附着流和后缘分离流的动态失速计算方法。


背景技术：

2.叶片动态失速是指当叶片沿变桨轴俯仰振荡运动时，其动态失速攻角比静态失速攻角要大，且翼型动态升力系数随攻角变化出现迟滞环的现象。其主要特征是叶片周围流场存在复杂的非定常分离涡和大尺度漩涡结构，涡状紊流从翼型前缘或后缘区域分离，引发气动升力急剧下降、阻力迅速增大的失速问题，严重破坏流场稳定性。
3.大气湍流、偏航和塔影效应、风剪切等多种客观因素导致风力机叶片长时间运行在动态失速区，而由二维静态升阻力气动数据得到的翼型气动载荷并不能真实反应其动态失速特性。因此，建立准确通用的翼型动态失速模型就显得尤为重要。


技术实现要素：

4.为此，本发明所要解决的技术问题在于提供一种精准通用的考虑翼型附着流和后缘分离流的动态失速计算方法。
5.为解决上述技术问题，本发明提供如下技术方案：
6.一种考虑翼型附着流和后缘分离流的动态失速计算方法，包括如下步骤：
7.步骤s1，输入翼型运行参数，确定翼型的几何攻角变化函数；
8.步骤s2，根据翼型的静态升力系数变化曲线，确定开始分离攻角、完全分离攻角、零升力攻角、线性升力系数斜率；
9.步骤s3，建立考虑翼型附着流和后缘分离流的动态失速模型；
10.步骤s4，采用四阶龙格库塔法求解动态失速模型；
11.步骤s5，计算得到翼型的非定常动态升力系数和动态分离过程。
12.上述一种考虑翼型附着流和后缘分离流的动态失速计算方法，在步骤s1中，翼型运行参数包括翼型的平均攻角、翼型的攻角幅值、翼型的来流风速、翼型的弦长、翼型的衰减频率、马赫数、模型的压力滞后时间常数、模型的边界层滞后时间常数。
13.上述一种考虑翼型附着流和后缘分离流的动态失速计算方法，翼型的几何攻角α为：
14.α＝a b
×
sin(2ukt/c)
15.式中，α为翼型的几何攻角，t为翼型的运行时间，a为翼型的平均攻角，b为翼型的攻角幅值，u为翼型的来流风速，k为翼型的衰减频率，c为翼型的弦长。
16.上述一种考虑翼型附着流和后缘分离流的动态失速计算方法，在步骤s2中，
17.翼型的两个完全分离攻角α
±
fs
的静态升力系数和两个开始分离攻角α
±
fz
的静态升力系数的通用表达式为：
[0018][0019][0020]
式中，α
±
fs
为翼型的两个完全分离攻角，α
±
fz
为翼型的两个开始分离攻角，为翼型的静态升力系数，为翼型的几何攻角等于α
±
fs
时的静态升力系数，为翼型的几何攻角等于α
±
fz
时的静态升力系数，c
l,a
为翼型的附着流线性区升力系数斜率，α0为翼型的零升力攻角；
[0021]
翼型的附着流线性区升力系数斜率c
l,a
为：
[0022][0023]
式中，c
l,a
为翼型的附着流线性区升力系数斜率，为翼型的几何攻角等于α时的静态升力系数，α为翼型的几何攻角；
[0024]
和根据翼型的静态升力系数与攻角变化曲线得到。
[0025]
上述一种考虑翼型附着流和后缘分离流的动态失速计算方法，在步骤s3中，考虑翼型附着流和后缘分离流的动态失速模型为：
[0026][0027][0028][0029][0030]
式中，u为翼型的来流风速，为翼型的来流风速变化速率，c为翼型的弦长，a1,a2,b1,b2为模型的四个经验常数，α为翼型的几何攻角，t
p
是模型的压力滞后时间常数，tf是模型的边界层滞后时间常数，为翼型的附着流非稳态升力系数，f
′
为翼型的等效准静态分离点，x1、x2、x3、x4为模型的四个状态变量，为模型的四个状态变量随时间的变化率；
[0031]
其中，翼型的附着流非稳态升力系数为：
[0032][0033]
式中，为翼型的附着流非稳态升力系数，αe为翼型的有效攻角，α0为翼型的零升力攻角；
[0034]
翼型的有效攻角αe为：
[0035]
αe＝α(1-a
1-a2) x1 x2。
[0036]
上述一种考虑翼型附着流和后缘分离流的动态失速计算方法，翼型的等效准静态
分离点f
′
为：
[0037][0038]
式中，f
′
为翼型的等效准静态分离点，αf为翼型的准静态攻角，为翼型的几何攻角等于αf时的静态升力系数，
[0039]
翼型的准静态攻角αf为：
[0040]
αf＝x3/c
l,α
α0[0041]
式中，αf为翼型的准静态攻角，x3为模型的状态变量，c
l,a
为翼型的附着流线性区升力系数斜率，α
±
fs
为翼型的两个完全分离攻角，α
±
fz
为翼型的两个开始分离攻角。
[0042]
上述一种考虑翼型附着流和后缘分离流的动态失速计算方法，在步骤s4中，采用四阶龙格库塔法求解动态失速模的四个一阶微分方程组：
[0043]
初始条件：t＝0时刻的四个状态变量x1、x2、x3、x4的值为[0；0；0；0]，相对误差为1e-5
，绝对误差为[1e-5
,1e-5
,1e-5
,1e-5
]；运行时间区间t＝[0,40]；最终求解得到x1、x2、x3、x4为四个列向量。
[0044]
上述一种考虑翼型附着流和后缘分离流的动态失速计算方法，在步骤s5中，翼型的非定常动态升力系数为：
[0045][0046]
式中，为翼型的非定常动态升力系数，c
l,a
为翼型附着流线性区升力系数斜率，为翼型在有效攻角αe时的完全分离流升力系数，x4为模型的状态变量，αe为翼型的有效攻角，α0为翼型的零升力攻角。
[0047]
上述一种考虑翼型附着流和后缘分离流的动态失速计算方法，
[0048]
翼型的完全分离流升力系数为：
[0049][0050]
式中，为翼型的完全分离流升力系数，为翼型的几何攻角等于αe时的静态升力系数，αe为翼型的有效攻角，α
±
fs
为翼型的两个完全分离攻角，α
±
fz
为翼型的两个开始分离攻角。
[0051]
本发明的技术方案取得了如下有益的技术效果：
[0052]
根据本发明一种考虑翼型附着流和后缘分离流的动态失速计算方法，考虑翼型附着流和后缘分离流分离过程，通过四个一阶微分方程建立了能够反应翼型动态失速过程的动态失速模型，利用四阶龙格库塔法求解得到不同几何攻角变化下的翼型非定常动态升力系数及翼型气流动态分离过程。该方法运算速度快，只需要输入风速、静态升力系数、时间常数等翼型运行参数，即可掌握不同几何攻角变化下的翼型动态分离过程，得到精度较高的非定常动态升力系数迟滞环曲线，为分析翼型气流分离过程、评估翼型动态失速性能以及进一步探究叶片动态失速影响因素、揭示叶片动态失速流动机理及等方面提供新选择，具有极其重要的科学意义和工程实用价值。
附图说明
[0053]
图1本发明实施例中一种考虑翼型附着流和后缘分离流的动态失速计算方法的总体流程图；
[0054]
图2本发明实施例中翼型的静态升力系数变化曲线图；
[0055]
图3本发明实施例中动态失速模型状态变量变化曲线图；
[0056]
图4本发明实施例中翼型的动态分离点随几何攻角变化曲线图；
[0057]
图5本发明实施例中翼型的升力系数随几何攻角变化曲线图。
具体实施方式
[0058]
本实施例中，以美国国家可再生能源实验室5mw海上风力机叶片du21翼型为例，对考虑翼型附着流和后缘分离流的动态失速计算方法进行具体说明。
[0059]
图1是本发明实施例中考虑翼型附着流和后缘分离流的动态失速计算方法的总体
流程图。
[0060]
具体计算方法如下。
[0061]
步骤s1，如表1所示，输入翼型du21的平均攻角、翼型的攻角幅值、翼型的来流风速、翼型的弦长、衰减频率、马赫数、压力滞后时间常数、边界层滞后时间常数等运行参数，确定翼型的几何攻角变化函数，然后进入步骤s2。
[0062]
表1翼型du21运行参数表
[0063][0064]
翼型几何攻角α为：
[0065]
α＝a b
×
sin(2ukt/c)
[0066]
式中，α为翼型的几何攻角，t为翼型的运行时间，a为翼型的平均攻角，b为翼型的攻角幅值，u为翼型的来流风速，k为翼型的衰减频率，c为翼型的弦长。
[0067]
由表1可知，翼型的运行时间t＝40秒，翼型的平均攻角a＝9
°
，翼型的攻角幅值b＝8
°
，翼型的来流风速u＝11.4m/s，翼型的衰减频率k＝0.075，翼型的弦长c＝3m。
[0068]
步骤s2，根据图2给出的翼型du21的静态升力系数变化曲线，确定开始分离攻角、完全分离攻角、零升力攻角、线性升力系数斜率，然后进入步骤s3。
[0069]
翼型的两个完全分离攻角α
±
fs
的静态升力系数和两个开始分离攻角α
±
fz
的静态升力系数的通用表达式为：
[0070][0071][0072]
式中，α
±
fs
为翼型的两个完全分离攻角，即流动完全分离时的攻角，包括α
fs
与α-fs
，α
±
fz
为翼型的两个开始分离攻角，即开始流动分离时的攻角，包括α
fz
与α-fz
，为翼型du21的静态升力系数；
[0073]
为翼型的几何攻角等于α
±
fs
时的静态升力系数，为翼型的几何攻角等于α
±
fz
时的静态升力系数，开始分离攻角的静态升力系数和完全分离攻角时的静态升力系数可通过图2插值得到。
[0074]cl,a
为翼型的附着流线性区升力系数斜率，α0为翼型的零升力攻角，α0＝-4.2
°
；
[0075]
翼型的附着流线性区升力系数斜率c
l,a
为：
[0076][0077]
式中，c
l,a
为翼型的附着流线性区升力系数斜率，可以由图2得到；为翼型
的几何攻角等于α时的静态升力系数，α为翼型的几何攻角，有步骤s1得到；
[0078]
步骤s3，建立考虑翼型附着流和后缘分离流的动态失速模型，然后进入步骤s4。
[0079]
考虑翼型附着流和后缘分离流的动态失速模型为：
[0080][0081][0082][0083][0084]
式中，u为翼型的来流风速，为翼型的来流风速变化速率，c为翼型的弦长，a1,a2,b1,b2为模型的四个经验常数(参见表1)，α为翼型的几何攻角，t
p
是模型的压力滞后时间常数，tf是模型的边界层滞后时间常数，为翼型的附着流非稳态升力系数，f
′
为翼型的等效准静态分离点，x1、x2、x3、x4为模型的四个状态变量，为模型的四个状态变量随时间的变化率；
[0085]
其中，翼型的附着流非稳态升力系数为：
[0086][0087]
式中，为翼型的附着流非稳态升力系数，αe为翼型的有效攻角，α0为翼型的零升力攻角(α0＝-4.2
°
)；
[0088]
翼型的有效攻角αe为：
[0089]
αe＝α(1-a
1-a2) x1 x2。
[0090]
其中，翼型的等效准静态分离点f
′
为：
[0091][0092]
式中，f
′
为翼型的等效准静态分离点，αf为翼型的准静态攻角，为翼型的几何攻角等于αf时的静态升力系数，
[0093]
翼型的准静态攻角αf为：
[0094]
αf＝x3/c
l,α
α0[0095]
式中，αf为翼型的准静态攻角，x3为模型的第三个状态变量，c
l,a
为翼型的附着流线性区升力系数斜率，α
±
fs
为翼型的两个完全分离攻角，α
±
fz
为翼型的两个开始分离攻角。
[0096]
步骤s4，采用四阶龙格库塔法求解动态失速模型，然后进入步骤s5。
[0097]
采用四阶龙格库塔法求解动态失速模的四个一阶微分方程组。初始条件：t＝0时刻的四个状态变量x1、x2、x3、x4的值为[0；0；0；0]，相对误差为1e-5
，绝对误差为[1e-5
,1e-5
,1e-5
,1e-5
]。运行时间区间t＝[0,40]，求解时长仅用了7.712秒。最终求解得到x1、x2、x3、x4为四个557
×
1的列向量，图3给出了本发明实施例中动态失速模型的四个状态变量随时间的变化曲线。
[0098]
步骤s5，计算得到翼型的非定常动态升力系数和动态分离过程。
[0099]
图4是本发明实施例中翼型的动态分离点随几何攻角变化曲线图。
[0100]
图5是本发明实施例中翼型的升力系数随几何攻角变化曲线图。
[0101]
翼型的非定常动态升力系数为：
[0102][0103]
式中，为翼型的非定常动态升力系数，c
l,a
为翼型附着流线性区升力系数斜率，由步骤s2得到，为翼型在有效攻角αe时的完全分离流升力系数，x4为模型的第四个状态变量，即动态分离点，αe为翼型的有效攻角，α0为翼型的零升力攻角。
[0104]
翼型的完全分离流升力系数为：
[0105][0106]
式中，为翼型的完全分离流升力系数，为翼型几何攻角等于有效攻角αe时的静态升力系数，可通过有效攻角αe插值图2给出的翼型du21的静态升力系数曲线得到，α
±
fs
为翼型的两个完全分离攻角，α
±
fz
为翼型的两个开始分离攻角。
[0107]
图4给出了翼型du21的动态分离点和随几何攻角变化曲线。动态分离点随着几何
攻角的增大，翼型动态失速现象逐渐明显，动态分离点位置由1减小到0.3，边界层分离由后缘逐渐向前缘移动。之后几何攻角减小，动态分离点位置由0.3增大到1，直到翼型不发生动态失速。
[0108]
图5给出的翼型du21的动态升力系数变化曲线和静态升力现象曲线偏差很大，偏离了静态值的同时还出现非常明显的迟滞环，相对于静态失速攻角9.3
°
时的极大升力系数1.403，动态升力系数在失速攻角17.85
°
时达到最大值2.317。失速攻角滞后8.55
°
，而最大动态升力系数相比最大静态升力系数增加65.2％。
[0109]
本实施例结果与被广大科研学者一致认可的风力机专用软件gh bladed计算结果完全一致，证实本实施例结果准确可靠。
[0110]
实施例作用与效果
[0111]
根据本发明一种考虑翼型附着流和后缘分离流的动态失速计算方法，考虑翼型附着流和后缘分离流分离过程，通过四个一阶微分方程建立了能够反应翼型动态失速过程的动态失速模型，利用四阶龙格库塔法求解得到不同几何攻角变化下的翼型非定常动态升力系数及翼型气流动态分离过程。
[0112]
该方法运算速度快，只需要输入风速、静态升力系数、时间常数等翼型运行参数，即可掌握不同几何攻角变化下的翼型动态分离过程，得到精度较高的非定常动态升力系数迟滞环曲线，为分析翼型气流分离过程、评估翼型动态失速性能以及进一步探究叶片动态失速影响因素、揭示叶片动态失速流动机理及等方面提供新选择，具有极其重要的科学意义和工程实用价值。
[0113]
显然，上述实施例仅仅是为清楚地说明所作的举例，而并非对实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。而由此所引伸出的显而易见的变化或变动仍处于本专利申请权利要求的保护范围之中。