一种风力机叶片抗颤振翼型簇设计方法与流程

1.本发明涉及属于风力机叶片翼型设计技术领域，具体涉及一种风力机叶片抗颤振翼型簇设计方法。


背景技术：

2.叶片主要是由复合材料制成的薄壁空心结构，随着叶片的长度增加，其刚度降低，柔性特征突出，易形成气动弹性不稳定性，即颤振。颤振会使叶片发生疲劳破坏，产生裂纹，严重时会导致叶片直接断裂。为了保证大型风力发电设备的正常工作，必须重视风力机叶片的抗颤振性能。
3.风力机叶片的气动外形作为捕获风能的关键因素，通过对风机叶片翼型的优化设计，提高其颤振临界速度，是改善抗颤振性能的有效方法；风能作为一种清洁、取之不尽的能源在世界中被广泛应用，伴随着风力发电技术不断升级，提高风力机叶片的抗颤振性能成为叶片设计的重点。为此，提出一种风力机叶片抗颤振翼型簇设计方法。


技术实现要素：

4.本发明所要解决的技术问题在于：如何通过对风机叶片翼型的优化设计，提高其颤振临界速度，改善抗颤振性能，提供了一种风力机叶片抗颤振翼型簇设计方法。
5.本发明是通过以下技术方案解决上述技术问题的，本发明包括以下步骤：
6.s1：翼型表达
7.对原有翼型添加形状扰动函数来表达新的翼型；
8.s2：颤振临界速度影响因素分析
9.对风机叶片进行敏感性分析，得到颤振临界速度影响因素，提高颤振临界影响速度因素中的扭转刚度以提高抗颤振性能；
10.s3：选择约束条件
11.确定翼型的约束条件为厚度约束条件与升阻比约束条件；
12.s4：确定翼型目标函数
13.将翼型在攻角α0下的升阻比与翼型的极惯性矩i作为目标函数，通过采用多目标遗传算法对风力机叶片翼型进行优化设计。
14.更进一步地，在所述步骤s1中，所述原有翼型为荷兰du翼型系列中的du35。
15.更进一步地，在所述步骤s1中，在原有翼型基础上添加形状扰动函数，保证x轴坐标不变，得到新翼型的型值点坐标：
16.上翼面翼型的型值点坐标：
17.y
u
(i)＝y
ou
(i) δy
u
(i)
18.x
u
(i)＝x
ou
(i)
19.下翼面翼型的型值点坐标：
20.y
l
(i)＝y
ol
(i) δy
l
(i)
21.x
l
(i)＝x
ol
(i)
22.其中，y
u
(i)代表新翼型上翼面第i个点的纵坐标，y
ou
(i)代表原翼型上翼面第i个点的纵坐标，δy
u
(i)代表上翼面第i个点的形状扰动函数，y
l
(i)代表新翼型下翼面第i点的纵坐标，y
ol
(i)代表原始翼型上翼面第i个点的横坐标，δy
l
(i)代表下翼面第i个点的形状扰动函数；x
u
(i)代表新翼型上翼面第i个点的横坐标，x
ou
(i)代表原始翼型上翼面第i个点的横坐标，x
ol
(i)代表原始翼型下翼面第i个点的横坐标。
23.更进一步地，上、下翼面形状扰动函数如下：
24.上翼面形状扰动函数：
[0025][0026]
下翼面形状扰动函数：
[0027][0028]
其中，下标u和l分别代表上翼面和下翼面，i代表翼型型值点索引，k与f(k)分别是控制翼型厚度与弯度分布关键点变量的个数与系数，g为给点矢量，f(k)代表hicks
‑
henne型函数，ξ和η为指数因子。
[0029]
更进一步地，在所述步骤s2中，翼型的扭转刚度表达式如下：
[0030][0031]
其中，g为剪切模量，r为翼型所在界面距离叶根的半径，i为翼型极惯性矩。
[0032]
更进一步地，所述极惯性矩i的计算过程如下：
[0033]
s21：计算出翼型的面积和质心坐标点：
[0034]
翼型的面积：
[0035][0036]
其中，s为翼型的面积；b为翼型弦长，单位标准翼型值弦长取1；u代表横坐标，v代表纵坐标；θ1(u)表示翼型下翼面曲线函数；θ2(u)为上翼面曲线函数；
[0037]
翼型质心坐标：
[0038][0039][0040]
其中，ρ(u，v)代表不同坐标点处的密度，假设叶片是均匀的，则ρ(u，v)为常数；
[0041]
则上式可进一步推导：
[0042]
[0043][0044]
s22：坐标原点位于翼型的前缘点，求出翼型对x，y方向的惯性矩：
[0045][0046][0047]
s23：平移坐标系，使原点与翼型质心重合，通过平行轴定理，求出翼型相对质心轴系的惯性矩：
[0048]
i1＝i
x
‑
ay2[0049]
i2＝i
y
‑
ax2[0050]
其中：i1代表y轴相对于质心的惯性矩；
[0051]
a代表翼型的面积；
[0052]
y代表y轴与质心轴的垂直距离；
[0053]
i2代表x轴相对于质心的惯性矩；
[0054]
x代表x轴与质心的垂直距离；
[0055]
则翼型极惯性矩为：
[0056]
i＝i1 i2。
[0057]
更进一步地，在所述步骤s3中，厚度约束条件为t1≤t≤t2，式中t为翼型厚度，t1、t2为设定值。
[0058]
更进一步地，在所述步骤s3中，优化后的翼型在不同攻角下的升阻比不低于原翼型的升阻比，选择[α1，α2，
…
，α
i
]攻角下的升阻比作为约束条件，升阻比约束条件为：
[0059]
c
l
′
/c
d
′
≥c
l
/c
d
[0060]
式中c
l
'为优化后的升力系数，c
d
'为优化后的阻力系数，c
l
为原有翼型的升力系数，c
d
为原有翼型的阻力系数。
[0061]
更进一步地，在所述步骤s4中，目标函数为：
[0062]
maximize＝[f1(x)，f2(x)]
[0063]
f1(x)＝i
[0064]
f2(x)＝c
l
/c
d
。
[0065]
更进一步地，在所述步骤s4中，多目标遗传算法为nsga
‑
ii。
[0066]
本发明相比现有技术具有以下优点：
[0067]
1)通过采用三次样条拟合翼型的型值点，得到了上下翼面的曲线函数，将翼型分割成为50等分，得到51个点的横坐标值，带入到上下翼面曲线函数中，得到51个点横坐标对应的纵坐标，其中上下翼面各51个点，确保了重新拟合的翼型曲线与原曲线一致，采用分段三次样条来拟合这51个点，从而得到50个子域；利用此方法可以得到极惯性矩的具体数值，并且实现优化寻优过程中的自动计算。
[0068]
2)基于isight优化设计平台实现了在迭代进程中自动更新翼型几何外形、网格模型、自动求解目标及约束函数值，将翼型几何外形更新程序、翼型流场ug几何生成模块、icem网格划分模块、cfx模块集成到一起，进行数据传递，从而实现了优化数学模型的迭代
进程。
[0069]
3)提出了风力机叶片抗颤振翼型多目标优化设计方法，建立了翼型截面极惯性矩与翼型几何外形的关系，通过多目标遗传算法得到了优化后的翼型簇集，其不仅提高了抗颤振稳定性能并且增大了其气动性能(升阻比更大)。
[0070]
4)通过在原翼型上添加形状扰动函数，不断在原有翼型基础上更新，直到寻找到帕累托最优解集。
[0071]
5)优化后的翼型集合与原翼型相比较，在大攻角条件下具有良好的抑制层流分离的能力，实现转捩控制，使得翼型失速延期。
附图说明
[0072]
图1是本发明实施例一中的原始翼型图(du35)；
[0073]
图2是本发明实施例一中翼型表面压力及剪力分布示意图；
[0074]
图3是本发明实施例一中优化设计的帕累托最优解集；
[0075]
图4a、b、c是本发明实施例一中的优化后三种翼型与原始翼型的比较图；
[0076]
图5是本发明实施例一中优化后翼型与原始翼型的升阻比比较图；
[0077]
图6是本发明实施例一中的isight平台流程图；
[0078]
图7是本发明实施例二中优化设计中的翼型流场几何域；
[0079]
图8是本发明实施例二中优化设计中的网络拓扑图；
[0080]
图9a是本发明实施例二中优化设计的翼型全局计算网格；
[0081]
图9b是本发明实施例二中优化设计的翼型局部计算网格；
[0082]
图10是本发明实施例二中优化设计中的初始翼型和优化翼型绕流流线图。
具体实施方式
[0083]
下面对本发明的实施例作详细说明，本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。
[0084]
实施例一
[0085]
本实施例以du35为设计原始翼型，但该方法不限于此翼型，本发明从翼型的形状拟合、极惯性矩的计算方法、数值计算的方法选取等方面进行了创新设计。
[0086]
翼型的拟合曲线
[0087]
本实施例是选取函数线性叠加的方法来表达新的翼型，通过在原有翼型基础上添加形状扰动函数，以幅值参数和指数因子为设计变量，对原有翼型进行更新。
[0088]
极惯性矩的计算
[0089]
在对极惯性矩计算过程中，首先建立极惯性矩与几何外形的关系，然后计算出翼型的面积和质心坐标；平移坐标系，使坐标原点和质心重合，根据平行轴原理，可求出翼型相对于质心轴系的惯性矩，最终得到翼型的极惯性矩。
[0090]
数值计算方法
[0091]
前面虽然给出了极惯性矩计算的具体方法，但要得到具体数值，并实现优化过程中寻找最优解的自动计算，需要采用数值计算的方法，采用高斯
‑
勒让德积分对求解域进行
积分计算。
[0092]
(1)借助ug生成计算域的确定几何模型，流场区域采用半椭圆型流场，这样做的目的是方便引入不同的入流攻角。利用ugupdate，此为ug和ansys的集成计算算程序，将ug的prt文件先储存成x_t文件格式，再导入ansys进行计算，既可以更新prt文件，又可以存储成x_t文件。
[0093]
(2)利用cad软件ug导出的parasolid格式的文件，并将其导入到icem中进行网格划分。应该注意网格划分不能过密或者过疏，否则会造成计算精度过差或者计算资源浪费。
[0094]
(3)在cfxpre中应该合理的创建计算域和边界条件，有利于控制数值计算的精度。
[0095]
本实施例的风力机叶片抗颤振翼型多目标优化设计方法，具体包括以下步骤：
[0096]
步骤一：翼型的表达方式是进行优化的先决条件，本发明的翼型的几何形状是在原有翼型基础上添加形状扰动函数，本发明选取解析函数线性叠加法来表达新的翼型。既继承了原翼型的总体特征，又可以节省时间。原始翼型如图1所示。
[0097]
在原有翼型基础上添加形状扰动函数，保证x轴坐标不变，得到新翼型的型值点坐标。
[0098]
上翼面翼型的型值点坐标：
[0099]
y
u
(i)＝y
ou
(i) δy
u
(i)
[0100]
x
u
(i)＝x
ou
(i)
[0101]
下翼面翼型的型值点坐标：
[0102]
y
l
(i)＝y
ol
(i) δy
l
(i)
[0103]
x
l
(i)＝x
ol
(i)
[0104]
其中，y
u
(i)代表新翼型上翼面第i个点的纵坐标，y
ou
(i)代表原翼型上翼面第i个点的纵坐标，δy
u
(i)代表上翼面第i个点的形状扰动函数，y
l
(i)代表新翼型下翼面第i点的纵坐标，y
ol
(i)代表原始翼型上翼面第i个点的横坐标，δy
l
(i)代表下翼面第i个点的形状扰动函数；x
u
(i)代表新翼型上翼面第i个点的横坐标，x
ou
(i)代表原始翼型上翼面第i个点的横坐标，x
ol
(i)代表原始翼型下翼面第i个点的横坐标；
[0105]
其中，上翼面形状扰动函数：
[0106][0107]
下翼面形状扰动函数：
[0108][0109]
式中下标u和l分别代表上翼面和下翼面，i代表翼型型值点索引，k代表形状模式索引，g为给点矢量，f(k)代表hicks
‑
henne型函数，ξ和η指数因子。
[0110]
取g(k)＝[0.1，0.2，0.3，0.5，1，3，5，8]；k与f(k)分别是控制翼型厚度与弯度分布关键点变量的个数与系数。采用hicks
‑
henne型函数，对翼型上下表面各关键点形状的影响平滑均衡，既可以保证翼型关键点的形状调整，又可以保证关键点与翼型形状平滑。
[0111]
所述步骤一的原翼型选择为荷兰du翼型系列中的du35。
[0112]
步骤二：翼型极惯性矩计算：
[0113]
对于风机叶片进行敏感性分析，得到扭转刚度、重心位置、空气密度及弹性轴位置
这四个因素对于颤振临界速度影响较大。将重心位置，弹性轴位置向前缘移动或者将叶片放在低密度区域均可以提高抗颤振性能，但这些都是不易实现的。我们对翼型的扭转刚度进行分析，发现提高翼型的扭转刚度可以提高抗颤振性能，翼型的扭转刚度与三个因素有关；
[0114][0115]
其中，为扭转刚度；g为剪切模量，与材料的选择有关，剪切模量高的材料其扭转刚度好，则抗颤振性能越高；
[0116]
r为翼型所在界面距离叶根的半径，随着半径的增大，其扭转刚度变小，则抗颤振性能减小；
[0117]
i为翼型极惯性矩，与翼型外形密切相关，翼型抗颤振性能随极惯性矩增大而增加；
[0118]
极惯性矩i取决于翼型几何外形，但极惯性矩i为几何外形的隐式表达，需要建立几何外形与极惯性矩之间的关系。首先，计算出翼型的面积和质心坐标点。
[0119]
求翼型的面积：
[0120][0121]
其中，s为翼型的面积；b为翼型弦长，单位标准翼型值弦长取1；u代表横坐标，v代表纵坐标；θ1(u)表示翼型下翼面曲线函数；θ2(u)为上翼面曲线函数。
[0122]
翼型质心坐标：
[0123][0124][0125]
其中，ρ(u，v)代表不同坐标点处的密度，故假设叶片是均匀的，则ρ(u，v)为常数。
[0126]
则上式可进一步推导：
[0127][0128][0129]
坐标原点位于翼型的前缘点，可求出翼型对x，y方向的惯性矩：
[0130][0131]
[0132]
平移坐标系，使原点与翼型质心重合，通过平行轴定理，求出翼型相对质心轴系的惯性矩
[0133]
i1＝i
x
‑
ay2[0134]
i2＝i
y
‑
ax2[0135]
其中：i1代表y轴相对于质心的惯性矩；
[0136]
a代表翼型的面积；
[0137]
y代表y轴与质心轴的垂直距离；
[0138]
i2代表x轴相对于质心的惯性矩；
[0139]
x代表x轴与质心的垂直距离；
[0140]
翼型极惯性矩为：
[0141]
i＝i1 i2[0142]
步骤三：约束条件的选择，翼型沿着风力机叶片展向分布时，其在不同叶片处存在着厚度要求，所以要对翼型厚度t进行约束，同时要保证优化后的翼型在不同攻角下的升阻比不低于原翼型的升阻比，选择[α1，α2，
…
，α
i
]攻角下的升阻比作为约束条件；
[0143]
厚度约束条件t1≤t≤t2[0144]
升阻比约束条件c
l
′
/c
d
′
≥c
l
/c
d
[0145]
式中c
l
'为优化后的升力系数，c
d
'为优化后的阻力系数，c
l
为原翼型的升力系数，c
d
为原翼型的阻力系数。
[0146]
本实施例选择du35翼型为案例，其厚度约束条件：
[0147]
0.3≤t≤0.4
[0148]
在追求极惯性矩的同时，肯定会影响翼型的气动性能。升阻比作为衡量翼型气动性能的重要指标，因为du35翼型在8
°
攻角附近其升阻比最大(在步骤四中将8
°
攻角下的升阻比作为一个目标函数)，所以选择[0
°
，2
°
，4
°
，6
°
，10
°
，12
°
，14
°
]攻角下的升阻比作为约束条件。
[0149]
升阻比约束条件：
[0150]
c
l
′
/c
d
′
≥c
l
/c
d
[0151]
其中，c
l
′
为优化后的升力系数，c
d
′
为优化后的阻力系数，c
l
为原翼型的升力系数，c
d
为原翼型的阻力系数；
[0152]
采用计算流体力学数值计算方法计算翼型上下翼面的压力及剪切力，从而得到翼型上下表面的升阻力系数，并进一步得到翼型的升阻比，优化后翼型升阻比不低于原翼型。翼型表面压力和剪力分布示意图如图2所示。
[0153]
步骤四：
[0154]
通过对步骤2进行分析，随着叶片增长，其半径r必然增加，则叶片扭转刚度减小，增加了叶片的柔性，使得叶片更容易发生颤振。翼型极惯性矩越大，其扭转刚度越高，则抗颤振性能越强，故将翼型的极惯性矩i作为一个设计目标；不同的翼型在不同的攻角下具有可变的升阻比，这里，先定义一个攻角α0，在此攻角下翼型具有最大的初始升阻比，将翼型在α0下的升阻比作为另一个目标函数。
[0155]
在isight平台中通过采用多目标遗传算法nsga
‑
ii对风力机叶片翼型进行优化设计。利用isight平台对matlab、ug、icem、cfx等软件进行集成运算，建立仿真分析流程，设定
和修改设计变量以及约束条件和目标函数，自动进行多次分析循环。可以自动的进行翼型坐标点的更新，自动进行几何建模、网格划分以及cfx的自动化计算。
[0156]
在本实施例中，将极惯性矩和8
°
攻角下的升阻比作为目标函数：
[0157]
maximize＝[f1(x)，f2(x)]
[0158]
f1(x)＝i
[0159]
f2(x)＝c
l
/c
d
[0160]
在isight平台中可以对cad/cae软件进行集成，用户在isight平台中自行修改设计变量、约束条件和升阻比，采用多目标遗传因子算法nsga
‑
ii对仿真优化流程自动进行循环。多目标优化的帕累托最优解集如图3所示。
[0161]
在帕累托最优解集中选择三组翼型与原有翼型比较，如图4a，4b，4c所示；
[0162]
选择的三组翼型与原有翼型的升阻比进行比较，如图5所示。
[0163]
步骤五：通过步骤一到步骤四，实现了对翼型抗颤振多目标的设计。此优化流程设计过程中，重点在于如何实现迭代进程中自动更新翼型的几何模型，网格划分，自动求解目标以及求解约束函数。本案例中将翼型几何外形更新程序(matlab文件)、翼型流场几何生成模块(ug文件)、网格划分模块(icem)、数值计算模块(cfx)共同集成到isight平台中，实现数据传递并进行优化迭代进程。同时isight平台可以将设计变量，约束条件，目标函数值的变化历程以可视化的方式直观的显示出来，还可以提供多种数据分析和后处理工具。在isight平台的流程图如图6所示。
[0164]
实施例二
[0165]
本实施例以du35为原始翼型，该翼型最大相对厚度为0.35，位于27.50％翼弦处，该翼型弧度为
‑
1.96％，位于30％的翼弦处。
[0166]
计算域：
[0167]
借助ug生成确定计算域的几何模型，流场区域采用半椭圆几何形状，完整的椭圆长轴为100b短轴长度为50b，采用布尔运算从半椭圆流场中减去翼型几何曲线构成的几何体，便可得到翼型绕流流场的计算域的几何模型。如图7所示。
[0168]
计算网格：
[0169]
从ug导出parasolid格式文件，然后导入到网格划分软件icem中。在翼型边界层区域，考虑到近壁面粘滞作用，分别采用内o
‑
block和外o
‑
block拓扑生成边界层网格，首先，使用内o
‑
block确定翼型几何拓扑，并将内部block删除，其次，使用外o
‑
block确定翼周边界层拓扑区域，以生成边界层网格。近壁层第一层网格高度为0.2mm，确保粘性层第一层至少有一层网格。对网格数目进行划分，网格数为82852。网格拓扑图如图8所示，网格图如图9a、b所示。
[0170]
计算方法：
[0171]
在cfx
‑
pre中对流体分析进行前处理，设置参数。设置流体类型为理想气体，收敛精度为1e
‑
6。cfx提供了丰富的后处理选项，通过后处理可以直观地显示出流场计算结果，部分流线图如图10所示。
[0172]
综上所述，上述实施例的风力机叶片抗颤振翼型簇设计方法，通过在原有的翼型基础上添加翼型扰动函数来修改翼型的几何外形，提高其抗颤振稳定性能并增强其气动性能(升阻比更大)；通过经典抗颤振理论的研究，得出翼型的极惯性矩对其抗颤振性能具有
突出的影响，其值可通过高斯
‑
勒让德积分计算得出；并将最大升阻比和极惯性矩作为优化目标，在isight平台中融合翼型扰动函数、高斯
‑
勒让德积分计算、网格划分、流体力学计算、多目标遗传算法来获取优化后的新翼型簇,其抗颤振性能及气动性能均得到了显著的提高，值得被推广使用。
[0173]
尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。