一种抗气动弹性变形的风力机翼型设计方法

摘要

本发明公开了一种抗气动弹性变形的风力机翼型设计方法，结合翼型参数化表征方法与动量叶素理论，建立了翼型优化设计数学模型，考虑气弹变形对翼型的影响，设计出了一种抗气动弹性变形的风力机翼型WQ?D180，并将新翼型与传统翼型进行了气弹变形及气动性能分析，表明该翼型具有良好的气动性能及抗气弹变形性能。由于该翼型具有抗气弹变形的能力，即叶片实际运行时风载荷对叶片翼型截面变形影响小，使得叶片实际性能与理论性能相差不大，从而提高了叶片实际运行性能。

说明书

技术领域

本发明属于风力机技术领域，涉及一种抗气动弹性变形的风力机翼型设计方法， 尤其涉及一种利用有限元法与翼型集成理论相结合的抗气动弹性变形的风力机翼型设计 新方法。

背景技术

风力机叶片截面气动弹性变形主要是受到气动力、惯性力及弹性力的耦合作用 下，使叶片翼型本身发生变形。而叶片截面主要由气动外形及复合材料内部结构组成，以往 研究大都集中在如何设计高性能的风力机翼型，而忽略了翼型耐气弹变形的能力。由于风 载荷的作用，使得叶片翼型本身发生变形，实际叶片翼型气动力将偏离原设计气动载荷，使 得叶片实际运行时翼型气动性能不如未变形翼型。

发明内容

为了解决上述技术问题，本发明提供了一种利用有限元法与翼型集成理论相结合 的抗气动弹性变形的风力机翼型设计新方法。

本发明所采用的技术方案是：一种抗气动弹性变形的风力机翼型设计方法，其特 征在于，包括以下步骤：

步骤1：基于翼型参数化集成表达式，确定目标函数1和目标函数2，设计变量及约 束条件；

步骤2：采用粒子群算法进行抗气动弹性翼型优化，初始化步骤1中的变量；

步骤3：将初始化的变量导入翼型参数化集成表达式中，形成初始翼型集，并采用 约束条件过滤掉不符合翼型特性的几何廓线；

步骤4：判断初始翼型集中的元素是否为翼型；

若是，则执行下述步骤5；

若否，则回转执行上述步骤2；

步骤5：通过步骤1中得到的目标函数1和目标函数2，通过目标函数1和目标函数2 表达式计算翼型适应度值，该适应度值要求具有高的气动特性及抗气动弹性变形性能；

步骤6：根据适应度值更新初始翼型集中个体最优和全局最优解；

步骤7：判断是否满足终止条件；

若否，则进行粒子群算参数自适应调整，并回转执行上述步骤3；

若是，则输出新翼型。

作为优选，步骤1中是采用泛函集成理论对翼型进行参数化表征，以气弹变形最小 及最大升阻比为目标函数，设计变量及约束条件；

其中翼型参数化集成表达式为：

$\left(\begin{array}{c}x=(r+a^2/r)cos\theta \\ y=(r-a^2/r)sin\theta \end{array}\right);$

式中，x为翼型横坐标，y为翼型纵坐标；r为翼型在平面ζ中的矢径，a_k、b_k为系数，θθ 为幅角，a为1/4翼型弦长；为翼型控制方程函数；

其中以光滑条件下升阻比最大为目标函数1：

f₁(x)＝max(C_L/C_D)；

其中C_L为升力系数，C_D为阻力系数；

以叶片截面最小气弹变形为目标函数2：

f₂(x)＝min(L_i-L_i-1)；

其中L_i为第i次气弹变形翼型尾缘处的位移，L_i-1为第i-1次翼型尾缘处的位移；

选取翼型控制方程函数的第1到第8项系数作为优化设计的变量，确定设计变 量为：

X＝(a₁,b₁,a₂,b₂,a₃,b₃,a₄,b₄)；

对翼型控制方程前8项系数进行约束，其约束条件如表1：

表1设计变量范围

对翼型最大相对厚度进行约束，选取叶尖附近翼型，最大相对厚度约为18％。

作为优选，步骤5中所述计算目标函数1的抗气动弹性翼型适应度值，是计算一定 攻角条件下的气动参数，所述气动参数包括升力系数、阻力系数及升阻比。

作为优选，步骤5中所述计算目标函数2的抗气动弹性翼型适应度值，其具体实现 包括以下子步骤：

步骤5.1：利用目标函数2优化翼型；

步骤5.2：利用叶片截面气动弹性耦合有限元模型计算并提取翼型变形前后升力 大小；

步骤5.3：将初始翼型导入到叶片截面气动弹性耦合有限元模型中并施加气动载 荷；

步骤5.4：采用有限元法计算并提取变形后的翼型数据；

步骤5.5：将变形后翼型的升力大小与变形前翼型的升力大小进行比较；

若比较结果大于预设阀值，则回转执行所述步骤5.1；

若比较结果小于或等于预设阀值，则表明该设计翼型具有抗气动弹性特性，求解 结束输出到适应度值中。

作为优选，步骤5.2中所述叶片截面气动弹性耦合有限元模型，其构建过程包括以 下子步骤：

(1)基于复杂叶片表面形状函数拟合的方法，编制MATLAB与APDL语言之间的数据 传递机制，建立叶片二元翼段参数化模型；

(2)采用shell181单元，建立单向层合板、双向层合板及三向层合板，生成叶片气 动外形；然后定义各层的材料密度、厚度及铺设角；

(3)将叶片截面简化成变截面的悬臂梁，在此处附加施加约束对其施加气动载荷， 对其进行气动弹性变形分析。

(4)通过叶片气动力计算，编制压力分布数值拟合程序，并建立气动力多项式函数 与叶片二元翼段结构的数据传递机制，实现气动力与弹性体的耦合模型。

作为优选，步骤5.3中所述将初始翼型施加气动载荷，作用在叶片气动载荷即压力 分布为：

$p=\frac{1}{2}{\mathrm{\rho v}}_{rel}^2{C}_p+p_{\infty };$

其中ρ为空气密度取1.205kg/m³；C_p为压力系数，p_∞为标准大气压强；v_rel为作用在 叶片上的相对速度；

$v_{rel}=\sqrt{v_x^2+v_y^2}=\sqrt{{(1-a)}^2{v}_0^2+{(1+b)}^2{\left(\omega r\right)}^2};$

其中，v₀为风速；ω为风轮旋转角速度(rad/s)；r为叶素沿叶片展向位置；a、b为轴 向诱导因子和周向诱导因子；v_x、v_y为来流风速在x、y方向的速度分量；

$a=\frac{2+Y_1-\sqrt{4Y_1(1-F)+Y_1^2}}{2(1+{\mathrm{FY}}_1)};$

$b=\frac{1}{(1-aF)Y_2/(1-a)-1};$

式中Y₂＝4Fsinφcosφ/(σC_tF₁)；C_t＝C_lsinφ-C_dcos φ；σ＝Bc/(2πr)，B为叶片数；F₁为法向力系数和切向力系数修正因子；

式中，R为风轮半径；g＝exp[-0.125(Bλ-21)]+0.1；为叶尖速比，Ω为风轮 转速。

作为优选，步骤5.5中所述预设阀值为0.02。

作为优选，步骤7中所述终止条件是迭代数达到400。

本发明专利的有益效果：

结合翼型参数化表征方法与动量叶素理论，建立了翼型优化设计数学模型，考虑 气弹变形对翼型的影响，设计出了一种抗气动弹性变形的风力机翼型WQ-D180，并将新翼型 与传统翼型进行了气弹变形及气动性能分析，表明该翼型具有良好的气动性能及抗气弹变 形性能。由于该翼型具有抗气弹变形的能力，即叶片实际运行时风载荷对叶片翼型截面变 形影响小，使得叶片实际性能与理论性能相差不大，从而提高了叶片实际运行性能。

附图说明

图1：本发明实施例的流程图；

图2：本发明实施例中计算目标函数2的抗气动弹性翼型适应度值的流程图；

图3：本发明实施例的叶片内部结构及材料示意图；

图4：本发明实施例的叶片截面气动弹性耦合有限元模型示意图；

图5：本发明实施例的典型叶片截面的拟合压力分布与离散压力分布对比图；

图6：本发明实施例的WQ-D180新翼型变形前后翼型示意图；

图7：本发明实施例的DU96W180变形前后翼型示意图；

图8(a)：本发明实施例的DU96W180翼型变形前后升力系数示意图；

图8(b)：本发明实施例的DU96W180翼型变化前后升阻比示意图；

图9(a)：本发明实施例的WQ-D180翼型变形前后升力系数示意图；

图9(b)：本发明实施例的WQ-D180翼型变化前后升阻比示意图。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发 明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不 用于限定本发明。

本发明基于叶素动量理论与复合材料铺层理论，推导出作用于叶片翼段气动载荷 计算数学表达式，针对某2MW实际叶片截面，考虑实际复合材料叶片铺层情况，采用有限元 法建立了叶片翼段气弹耦合设计模型；针对叶片变形较大的叶尖附近翼段，约束翼型最大 相对厚度为18％，在气弹载荷作用下设计抗气动弹性变形的风力机新型翼型。

请见图1，本发明提供的一种抗气动弹性变形的风力机翼型设计方法，包括以下步 骤：

步骤1：基于翼型参数化集成表达式，确定目标函数1和目标函数2，设计变量及约 束条件；

其中翼型参数化集成表达式为：

$\left(\begin{array}{c}x=(r+a^2/r)cos\theta \\ y=(r-a^2/r)sin\theta \end{array}\right)---\left(1\right);$

式中，x为翼型横坐标，y为翼型纵坐标；r为翼型在平面ζ中的矢径，a_k、b_k为系数，θ 为幅角，a为1/4翼型弦长；为翼型控制方程函数；

其中以光滑条件下升阻比最大为目标函数1：

f₁(x)＝max(C_L/C_D)；

其中C_L为升力系数，C_D为阻力系数；

以叶片截面最小气弹变形为目标函数2：

f₂(x)＝min(L_i-L_i-1)；

其中L_i为第i次气弹变形翼型尾缘处的位移，L_i-1为第i-1次翼型尾缘处的位移；

选取翼型控制方程函数的第1到第8项系数作为优化设计的变量，确定设计变 量为：

X＝(a₁,b₁,a₂,b₂,a₃,b₃,a₄,b₄)；

对翼型控制方程前8项系数进行约束，其约束条件如表1：

表1设计变量范围

对翼型最大相对厚度进行约束，选取叶尖附近翼型，最大相对厚度约为18％。

步骤2：采用粒子群算法进行抗气动弹性翼型优化，初始化步骤1中的变量；

步骤3：将初始化的变量导入翼型参数化集成表达式中，形成初始翼型集，并采用 约束条件过滤掉不符合翼型特性的几何廓线；

步骤4：判断初始翼型集中的元素是否为翼型；

若是，则执行下述步骤5；

若否，则回转执行上述步骤2；

步骤5：通过步骤1中得到的目标函数1和目标函数2，通过目标函数1和目标函数2 表达式计算翼型适应度值，该适应度值要求具有高的气动特性及抗气动弹性变形性能；

目标函数1采用RFOIL气动计算软件计算一定攻角条件下的升力系数、阻力系数及 升阻比等气动参数。而目标函数2则采用本专利提出的叶片截面气动弹性耦合有限元模型 计算并提取翼型变形前后升力大小。将初始翼型导入到叶片截面气动弹性耦合有限元模型 中去并施加气动载荷；采用有限元法计算并提取变形后的翼型数据，将变形后的翼型与变 形前翼型进行气动性能比较(升力大小比较)，如果相差很小(0.02)，则表明该设计翼型具 有抗气动弹性特性，求解结束输出到适应度值中。其具体过程流程图2所示(虚线箭头表示 变形后翼型气动变形特性)：其具体实现包括以下子步骤：

步骤5.1：利用目标函数2优化翼型；

步骤5.2：利用叶片截面气动弹性耦合有限元模型计算并提取翼型变形前后升力 大小；

针对某2MW实际风力机叶片延展向0.75～0.85R处翼段，该叶片翼段截面如图3所 示。叶片主梁主要承受挥舞方向的载荷，而腹板则主要承受拍打方向的载荷。

由于叶片沿展向变形较大的区域在靠近叶尖附近，而靠近叶根较厚叶片变形较 小。因此，本发明选择叶片相对厚度为18％翼型的二元翼段进行抗气弹特性的翼型设计。鉴 于MATLAB能够实现RFOIL、ANSYS之间的参数化建模与调用，采用ANSYS自带的APDL语言建立 叶片二元翼段参数化模型；RFOIL计算翼段的气动力；利用MATLAB编制气动插值程序将两者 结合起来设计抗气动弹性变形的风力机翼型。建立叶片截面气动弹性耦合有限元模型，其 构建过程包括以下子步骤：

(1)基于复杂叶片表面形状函数拟合的方法，编制MATLAB与APDL语言之间的数据 传递机制，建立叶片二元翼段参数化模型。

(2)根据图2，采用shell181单元，建立单向层合板、双向层合板及三向层合板，生 成叶片气动外形；然后定义各层的材料密度、厚度及铺设角。

(3)由于叶片气动力中心为延弦向0.25倍弦长处，相关资料将叶片截面简化成变 截面的悬臂梁，在此处附加施加约束对其进行弹性力学分析。鉴于此，本发明亦将此处作为 约束条件，施加气动载荷，对其进行气动弹性变形分析。

(4)通过叶片气动力计算，编制压力分布数值拟合程序，并建立气动力多项式函数 与叶片二元翼段结构的数据传递机制，实现叶片截面气动弹性耦合有限元模型，如图4所 示。

步骤5.3：将初始翼型导入到叶片截面气动弹性耦合有限元模型中并施加气动载 荷；

由于风力机叶片受到气动载荷的作用，使得叶片本身发生弹性变形，而这种变形 反过来又影响叶片气动特性。作用在叶片气动载荷即压力分布可表示为：

$p=\frac{1}{2}{\mathrm{\rho v}}_{rel}^2{C}_p+p_{\infty }---\left(3\right)$

其中ρ为空气密度取1.205kg/m³；C_p为压力系数，对于每个翼型截面，可采用RFOIL 软件比较准确的计算翼型在失速前的气动特性；p∞为标准大气压强；v_rel为作用在叶片上 的相对速度。在(3)式中，相对速度v_rel是叶片转速与风速的合成速度，基于修正动量叶素理 论即可求得相对速度。

$v_{rel}=\sqrt{v_x^2+v_y^2}=\sqrt{{(1-a)}^2{v}_0^2+{(1+b)}^2{\left(\omega r\right)}^2}---\left(4\right)$

其中，v₀为风速；ω为风轮旋转角速度(rad/s)；r为叶素沿叶片展向位置；a、b为轴 向诱导因子和周向诱导因子；vx、v_y为来流风速在x、y方向的速度分量；

轴向诱导因子a和周向诱导因子b为

$a=\frac{2+Y_1-\sqrt{4Y_1(1-F)+Y_1^2}}{2(1+{\mathrm{FY}}_1)}---\left(5\right)$

$b=\frac{1}{(1-aF)Y_2/(1-a)-1}---\left(6\right)$

式中Y2＝4Fsinφcosφ/(σC_tF₁)；C_t＝C_lsinφ-C_dcosφ；σ＝Bc/(2πr)，B为叶片数，取3。

F₁为法向力系数和切向力系数修正因子

式中，R为风轮半径；g＝exp[-0.125(Bλ-21)]+0.1；为叶尖速比，Ω为风轮 转速。

通过式(5)、(6)不断迭代求解a及b，代入式(4)求出其相对速度，再代入式(3)即可 计算出每个叶素上的真实压力分布。为了简化气动力计算，本发明选取沿叶片展向0.75～ 0.85翼段，对此翼段施加气动载荷。对计算出来的气动载荷以多项式函数的形式加载到叶 片二元翼段模型中。图5为典型叶片截面的拟合压力分布与离散压力分布对比图(风速 12.5m/s，额定转速为20r/min)。由图可以看出，采用多项式函数拟合，可以准确地表示叶片 截面离散的压力分布点。

步骤5.4：采用有限元法计算并提取变形后的翼型数据；

步骤5.5：将变形后翼型的升力大小与变形前翼型的升力大小进行比较；

若比较结果大于预设阀值0.02，则回转执行所述步骤5.1；

若比较结果小于或等于预设阀值0.02，则表明该设计翼型具有抗气动弹性特性， 求解结束输出到适应度值中。

步骤6：根据适应度值更新初始翼型集中个体最优和全局最优解；

步骤7：判断是否满足终止条件(迭代数为400)；

若否，则进行粒子群算参数自适应调整，并回转执行上述步骤3；

若是，则输出新翼型。

本实施例采用粒子群算法对耐气弹风力机翼型进行优化，通过不断迭代最终优化 设计出气弹变形前后气动性能影响小、翼型廓线变形不明显的风力机专用翼型。如图6所 示，为优化翼型气动弹性变形前后翼型廓线，取名为WQ-D180翼型。由图可见，该设计新翼型 变形前后廓线变化很小，这对翼型气动性能影响较小。

为了说明新翼型抗气动弹性性能优越，针对某2MW实际风力机叶片在叶片，沿展向 0.80R处相关参数为：叶片截面弦长为1.2m，扭角2.2°，局部攻角为8.8°，最大相对厚度约为 18％。选取荷兰翼型DU96-W-180与WQ-D180翼型进行弹性变形前后气动性能对比分析。如图 7为DU96-W-180翼型变形前后翼型廓线。由于气动载荷的作用，使得翼型尾缘附近的弹性变 形较大，而翼型前缘部分厚度相对较大，几乎没有变形。为了定量分析翼型几何参数的变 化，表2给出了这两种翼段气弹变形前后的关键参数。由表2可知，通过风力机翼段气弹耦合 变形，这两种翼型最大的变化在于最大弯度比变形前要大，这势必会对翼型的气动特性有 影响，而其他几何参数，如：最大厚度及其弦向位置变化不大。

表2两种翼型变形前后几何参数

为了说明气动弹性变形对翼型气动性能的影响，利用RFOIL软件模拟计算在粗糙 和光滑条件下这两种翼型气弹变形前后的气动数据。如图8(a)～9(b)所示，分别为两种翼 型变形前后的升力系数及升阻比对比。总体上来看，相比变形前，这两种翼型变形后的升力 系数无论是自由转捩(光滑条件)还是固定转捩(粗糙条件)，均有不同程度增大，如图8(a) 及9(a)所示；而对于升阻比，表3给出这两种翼型变形前后关键气动参数对比情况。结合图 表可知，对比这两种翼型变形前后气动特性情况，升阻比变化较大的为DU96-W-180翼型，在 自由转捩条件下，变形后最大升阻比降低了约3.5％。而WQ-D180新翼型变形前后，最大升力 系数、最大升阻比并无明显变化，表明WQ-D180翼型抗气动弹性变形能力优于DU96-W-180翼 型，主要原因是WQ-D180翼型最大弯度所处弦向位置比较靠后，翼型气动弹性变形前后其几 何轮廓线变化不大，如图6所示，因此翼型变形前后其气动性能变化并无明显变化。

表3翼型气动性能参数对比

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本 发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权 利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发 明的请求保护范围应以所附权利要求为准。