CN113569360A - 一种风力机叶片抗颤振翼型簇设计方法 - Google Patents

一种风力机叶片抗颤振翼型簇设计方法 Download PDF

Info

Publication number
CN113569360A
CN113569360A CN202110962098.9A CN202110962098A CN113569360A CN 113569360 A CN113569360 A CN 113569360A CN 202110962098 A CN202110962098 A CN 202110962098A CN 113569360 A CN113569360 A CN 113569360A
Authority
CN
China
Prior art keywords
airfoil
flutter
profile
resistant
cluster
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
Application number
CN202110962098.9A
Other languages
English (en)
Other versions
CN113569360B (zh
Inventor
孟瑞
赵如新
高强
卢宗远
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Anhui University of Technology AHUT
Original Assignee
Anhui University of Technology AHUT
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Anhui University of Technology AHUT filed Critical Anhui University of Technology AHUT
Priority to CN202110962098.9A priority Critical patent/CN113569360B/zh
Publication of CN113569360A publication Critical patent/CN113569360A/zh
Application granted granted Critical
Publication of CN113569360B publication Critical patent/CN113569360B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/17—Mechanical parametric or variational design
    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F03—MACHINES OR ENGINES FOR LIQUIDS; WIND, SPRING, OR WEIGHT MOTORS; PRODUCING MECHANICAL POWER OR A REACTIVE PROPULSIVE THRUST, NOT OTHERWISE PROVIDED FOR
    • F03D—WIND MOTORS
    • F03D1/00—Wind motors with rotation axis substantially parallel to the air flow entering the rotor 
    • F03D1/06—Rotors
    • F03D1/065—Rotors characterised by their construction elements
    • F03D1/0675—Rotors characterised by their construction elements of the blades
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/27—Design optimisation, verification or simulation using machine learning, e.g. artificial intelligence, neural networks, support vector machines [SVM] or training a model
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/12—Computing arrangements based on biological models using genetic models
    • G06N3/126—Evolutionary algorithms, e.g. genetic algorithms or genetic programming
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/04—Constraint-based CAD
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • 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

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Biophysics (AREA)
  • Computer Hardware Design (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Artificial Intelligence (AREA)
  • Evolutionary Biology (AREA)
  • Software Systems (AREA)
  • Computing Systems (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Mathematical Physics (AREA)
  • Mathematical Analysis (AREA)
  • Algebra (AREA)
  • Data Mining & Analysis (AREA)
  • Fluid Mechanics (AREA)
  • Medical Informatics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Physiology (AREA)
  • Genetics & Genomics (AREA)
  • Biomedical Technology (AREA)
  • Computational Linguistics (AREA)
  • Computational Mathematics (AREA)
  • General Health & Medical Sciences (AREA)
  • Molecular Biology (AREA)
  • Sustainable Development (AREA)
  • Sustainable Energy (AREA)
  • Chemical & Material Sciences (AREA)
  • Combustion & Propulsion (AREA)
  • Mechanical Engineering (AREA)
  • Structures Of Non-Positive Displacement Pumps (AREA)
  • Wind Motors (AREA)

Abstract

本发明公开了一种风力机叶片抗颤振翼型簇设计方法，属于风力机叶片翼型设计技术领域，包括以下步骤：S1：翼型表达；S2：颤振临界速度因素分析；S3：选择约束条件；S4：确定翼型目标函数。本发明通过在原有的翼型基础上添加翼型扰动函数来修改翼型的几何外形，提高其抗颤振稳定性能并增强其气动性能(升阻比更大)；通过经典抗颤振理论的研究，得出翼型的极惯性矩对其抗颤振性能具有突出的影响，其值可通过高斯‑勒让德积分计算得出；并将最大升阻比和极惯性矩作为优化目标，在isight平台中融合翼型扰动函数、高斯‑勒让德积分计算、网格划分、流体力学计算、多目标遗传算法来获取优化后的新翼型簇,其抗颤振性能及气动性能均得到了显著的提高。

Description

一种风力机叶片抗颤振翼型簇设计方法
技术领域
本发明涉及属于风力机叶片翼型设计技术领域，具体涉及一种风力机叶片抗颤振翼型簇设计方法。
背景技术
叶片主要是由复合材料制成的薄壁空心结构，随着叶片的长度增加，其刚度降低，柔性特征突出，易形成气动弹性不稳定性，即颤振。颤振会使叶片发生疲劳破坏，产生裂纹，严重时会导致叶片直接断裂。为了保证大型风力发电设备的正常工作，必须重视风力机叶片的抗颤振性能。
风力机叶片的气动外形作为捕获风能的关键因素，通过对风机叶片翼型的优化设计，提高其颤振临界速度，是改善抗颤振性能的有效方法；风能作为一种清洁、取之不尽的能源在世界中被广泛应用，伴随着风力发电技术不断升级，提高风力机叶片的抗颤振性能成为叶片设计的重点。为此，提出一种风力机叶片抗颤振翼型簇设计方法。
发明内容
本发明所要解决的技术问题在于：如何通过对风机叶片翼型的优化设计，提高其颤振临界速度，改善抗颤振性能，提供了一种风力机叶片抗颤振翼型簇设计方法。
本发明是通过以下技术方案解决上述技术问题的，本发明包括以下步骤：
S1：翼型表达
对原有翼型添加形状扰动函数来表达新的翼型；
S2：颤振临界速度影响因素分析
对风机叶片进行敏感性分析，得到颤振临界速度影响因素，提高颤振临界影响速度因素中的扭转刚度以提高抗颤振性能；
S3：选择约束条件
确定翼型的约束条件为厚度约束条件与升阻比约束条件；
S4：确定翼型目标函数
将翼型在攻角α0下的升阻比与翼型的极惯性矩I作为目标函数，通过采用多目标遗传算法对风力机叶片翼型进行优化设计。
更进一步地，在所述步骤S1中，所述原有翼型为荷兰DU翼型系列中的DU35。
更进一步地，在所述步骤S1中，在原有翼型基础上添加形状扰动函数，保证X轴坐标不变，得到新翼型的型值点坐标：
上翼面翼型的型值点坐标：
yu(i)＝you(i)+Δyu(i)
xu(i)＝xou(i)
下翼面翼型的型值点坐标：
yl(i)＝yol(i)+Δyl(i)
xl(i)＝xol(i)
其中，yu(i)代表新翼型上翼面第i个点的纵坐标，you(i)代表原翼型上翼面第i个点的纵坐标，Δyu(i)代表上翼面第i个点的形状扰动函数，yl(i)代表新翼型下翼面第i点的纵坐标，yol(i)代表原始翼型上翼面第i个点的横坐标，Δyl(i)代表下翼面第i个点的形状扰动函数；xu(i)代表新翼型上翼面第i个点的横坐标，xou(i)代表原始翼型上翼面第i个点的横坐标，xol(i)代表原始翼型下翼面第i个点的横坐标。
更进一步地，上、下翼面形状扰动函数如下：
上翼面形状扰动函数：
Figure BDA0003222689120000021
下翼面形状扰动函数：
Figure BDA0003222689120000022
其中，下标u和l分别代表上翼面和下翼面，i代表翼型型值点索引，k与f(k)分别是控制翼型厚度与弯度分布关键点变量的个数与系数，g为给点矢量，f(k)代表Hicks-Henne型函数，ξ和η为指数因子。
更进一步地，在所述步骤S2中，翼型的扭转刚度表达式如下：
Figure BDA0003222689120000023
其中，G为剪切模量，R为翼型所在界面距离叶根的半径，I为翼型极惯性矩。
更进一步地，所述极惯性矩I的计算过程如下：
S21：计算出翼型的面积和质心坐标点：
翼型的面积：
Figure BDA0003222689120000024
其中，S为翼型的面积；b为翼型弦长，单位标准翼型值弦长取1；u代表横坐标，v代表纵坐标；θ1(u)表示翼型下翼面曲线函数；θ2(u)为上翼面曲线函数；
翼型质心坐标：
Figure BDA0003222689120000031
Figure BDA0003222689120000032
其中，ρ(u，v)代表不同坐标点处的密度，假设叶片是均匀的，则ρ(u，v)为常数；
则上式可进一步推导：
Figure BDA0003222689120000033
Figure BDA0003222689120000034
S22：坐标原点位于翼型的前缘点，求出翼型对x，y方向的惯性矩：
Figure BDA0003222689120000035
Figure BDA0003222689120000036
S23：平移坐标系，使原点与翼型质心重合，通过平行轴定理，求出翼型相对质心轴系的惯性矩：
I1＝Ix-Ay2
I2＝Iy-Ax2
其中：I1代表Y轴相对于质心的惯性矩；
A代表翼型的面积；
y代表Y轴与质心轴的垂直距离；
I2代表X轴相对于质心的惯性矩；
x代表X轴与质心的垂直距离；
则翼型极惯性矩为：
I＝I1+I2。
更进一步地，在所述步骤S3中，厚度约束条件为t1≤t≤t2，式中t为翼型厚度，t1、t2为设定值。
更进一步地，在所述步骤S3中，优化后的翼型在不同攻角下的升阻比不低于原翼型的升阻比，选择[α1，α2，…，αi]攻角下的升阻比作为约束条件，升阻比约束条件为：
cl′/cd′≥cl/cd
式中cl'为优化后的升力系数，cd'为优化后的阻力系数，cl为原有翼型的升力系数，cd为原有翼型的阻力系数。
更进一步地，在所述步骤S4中，目标函数为：
maximize＝[f1(x)，f2(x)]
f1(x)＝I
f2(x)＝cl/cd。
更进一步地，在所述步骤S4中，多目标遗传算法为NSGA-II。
本发明相比现有技术具有以下优点：
1)通过采用三次样条拟合翼型的型值点，得到了上下翼面的曲线函数，将翼型分割成为50等分，得到51个点的横坐标值，带入到上下翼面曲线函数中，得到51个点横坐标对应的纵坐标，其中上下翼面各51个点，确保了重新拟合的翼型曲线与原曲线一致，采用分段三次样条来拟合这51个点，从而得到50个子域；利用此方法可以得到极惯性矩的具体数值，并且实现优化寻优过程中的自动计算。
2)基于isight优化设计平台实现了在迭代进程中自动更新翼型几何外形、网格模型、自动求解目标及约束函数值，将翼型几何外形更新程序、翼型流场ug几何生成模块、icem网格划分模块、cfx模块集成到一起，进行数据传递，从而实现了优化数学模型的迭代进程。
3)提出了风力机叶片抗颤振翼型多目标优化设计方法，建立了翼型截面极惯性矩与翼型几何外形的关系，通过多目标遗传算法得到了优化后的翼型簇集，其不仅提高了抗颤振稳定性能并且增大了其气动性能(升阻比更大)。
4)通过在原翼型上添加形状扰动函数，不断在原有翼型基础上更新，直到寻找到帕累托最优解集。
5)优化后的翼型集合与原翼型相比较，在大攻角条件下具有良好的抑制层流分离的能力，实现转捩控制，使得翼型失速延期。
附图说明
图1是本发明实施例一中的原始翼型图(DU35)；
图2是本发明实施例一中翼型表面压力及剪力分布示意图；
图3是本发明实施例一中优化设计的帕累托最优解集；
图4a、b、c是本发明实施例一中的优化后三种翼型与原始翼型的比较图；
图5是本发明实施例一中优化后翼型与原始翼型的升阻比比较图；
图6是本发明实施例一中的isight平台流程图；
图7是本发明实施例二中优化设计中的翼型流场几何域；
图8是本发明实施例二中优化设计中的网络拓扑图；
图9a是本发明实施例二中优化设计的翼型全局计算网格；
图9b是本发明实施例二中优化设计的翼型局部计算网格；
图10是本发明实施例二中优化设计中的初始翼型和优化翼型绕流流线图。
具体实施方式
下面对本发明的实施例作详细说明，本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。
实施例一
本实施例以DU35为设计原始翼型，但该方法不限于此翼型，本发明从翼型的形状拟合、极惯性矩的计算方法、数值计算的方法选取等方面进行了创新设计。
翼型的拟合曲线
本实施例是选取函数线性叠加的方法来表达新的翼型，通过在原有翼型基础上添加形状扰动函数，以幅值参数和指数因子为设计变量，对原有翼型进行更新。
极惯性矩的计算
在对极惯性矩计算过程中，首先建立极惯性矩与几何外形的关系，然后计算出翼型的面积和质心坐标；平移坐标系，使坐标原点和质心重合，根据平行轴原理，可求出翼型相对于质心轴系的惯性矩，最终得到翼型的极惯性矩。
数值计算方法
前面虽然给出了极惯性矩计算的具体方法，但要得到具体数值，并实现优化过程中寻找最优解的自动计算，需要采用数值计算的方法，采用高斯-勒让德积分对求解域进行积分计算。
(1)借助ug生成计算域的确定几何模型，流场区域采用半椭圆型流场，这样做的目的是方便引入不同的入流攻角。利用ugupdate，此为ug和ansys的集成计算算程序，将ug的prt文件先储存成x_t文件格式，再导入ansys进行计算，既可以更新prt文件，又可以存储成x_t文件。
(2)利用CAD软件ug导出的parasolid格式的文件，并将其导入到icem中进行网格划分。应该注意网格划分不能过密或者过疏，否则会造成计算精度过差或者计算资源浪费。
(3)在cfxpre中应该合理的创建计算域和边界条件，有利于控制数值计算的精度。
本实施例的风力机叶片抗颤振翼型多目标优化设计方法，具体包括以下步骤：
步骤一：翼型的表达方式是进行优化的先决条件，本发明的翼型的几何形状是在原有翼型基础上添加形状扰动函数，本发明选取解析函数线性叠加法来表达新的翼型。既继承了原翼型的总体特征，又可以节省时间。原始翼型如图1所示。
在原有翼型基础上添加形状扰动函数，保证X轴坐标不变，得到新翼型的型值点坐标。
上翼面翼型的型值点坐标：
yu(i)＝you(i)+Δyu(i)
xu(i)＝xou(i)
下翼面翼型的型值点坐标：
yl(i)＝yol(i)+Δyl(i)
xl(i)＝xol(i)
其中，yu(i)代表新翼型上翼面第i个点的纵坐标，you(i)代表原翼型上翼面第i个点的纵坐标，Δyu(i)代表上翼面第i个点的形状扰动函数，yl(i)代表新翼型下翼面第i点的纵坐标，yol(i)代表原始翼型上翼面第i个点的横坐标，Δyl(i)代表下翼面第i个点的形状扰动函数；xu(i)代表新翼型上翼面第i个点的横坐标，xou(i)代表原始翼型上翼面第i个点的横坐标，xol(i)代表原始翼型下翼面第i个点的横坐标；
其中，上翼面形状扰动函数：
Figure BDA0003222689120000061
下翼面形状扰动函数：
Figure BDA0003222689120000062
式中下标u和l分别代表上翼面和下翼面，i代表翼型型值点索引，k代表形状模式索引，g为给点矢量，f(k)代表Hicks-Henne型函数，ξ和η指数因子。
取g(k)＝[0.1，0.2，0.3，0.5，1，3，5，8]；k与f(k)分别是控制翼型厚度与弯度分布关键点变量的个数与系数。采用Hicks-Henne型函数，对翼型上下表面各关键点形状的影响平滑均衡，既可以保证翼型关键点的形状调整，又可以保证关键点与翼型形状平滑。
所述步骤一的原翼型选择为荷兰DU翼型系列中的DU35。
步骤二：翼型极惯性矩计算：
对于风机叶片进行敏感性分析，得到扭转刚度、重心位置、空气密度及弹性轴位置这四个因素对于颤振临界速度影响较大。将重心位置，弹性轴位置向前缘移动或者将叶片放在低密度区域均可以提高抗颤振性能，但这些都是不易实现的。我们对翼型的扭转刚度进行分析，发现提高翼型的扭转刚度可以提高抗颤振性能，翼型的扭转刚度与三个因素有关；
Figure BDA0003222689120000071
其中，
Figure BDA0003222689120000072
为扭转刚度；G为剪切模量，与材料的选择有关，剪切模量高的材料其扭转刚度好，则抗颤振性能越高；
R为翼型所在界面距离叶根的半径，随着半径的增大，其扭转刚度变小，则抗颤振性能减小；
I为翼型极惯性矩，与翼型外形密切相关，翼型抗颤振性能随极惯性矩增大而增加；
极惯性矩I取决于翼型几何外形，但极惯性矩I为几何外形的隐式表达，需要建立几何外形与极惯性矩之间的关系。首先，计算出翼型的面积和质心坐标点。
求翼型的面积：
Figure BDA0003222689120000073
其中，S为翼型的面积；b为翼型弦长，单位标准翼型值弦长取1；u代表横坐标，v代表纵坐标；θ1(u)表示翼型下翼面曲线函数；θ2(u)为上翼面曲线函数。
翼型质心坐标：
Figure BDA0003222689120000074
Figure BDA0003222689120000075
其中，ρ(u，v)代表不同坐标点处的密度，故假设叶片是均匀的，则ρ(u，v)为常数。
则上式可进一步推导：
Figure BDA0003222689120000076
Figure BDA0003222689120000077
坐标原点位于翼型的前缘点，可求出翼型对x，y方向的惯性矩：
Figure BDA0003222689120000078
Figure BDA0003222689120000079
平移坐标系，使原点与翼型质心重合，通过平行轴定理，求出翼型相对质心轴系的惯性矩
I1＝Ix-Ay2
I2＝Iy-Ax2
其中：I1代表Y轴相对于质心的惯性矩；
A代表翼型的面积；
y代表Y轴与质心轴的垂直距离；
I2代表X轴相对于质心的惯性矩；
x代表X轴与质心的垂直距离；
翼型极惯性矩为：
I＝I1+I2
步骤三：约束条件的选择，翼型沿着风力机叶片展向分布时，其在不同叶片处存在着厚度要求，所以要对翼型厚度t进行约束，同时要保证优化后的翼型在不同攻角下的升阻比不低于原翼型的升阻比，选择[α1，α2，…，αi]攻角下的升阻比作为约束条件；
厚度约束条件t1≤t≤t2
升阻比约束条件cl′/cd′≥cl/cd
式中cl'为优化后的升力系数，cd'为优化后的阻力系数，cl为原翼型的升力系数，cd为原翼型的阻力系数。
本实施例选择DU35翼型为案例，其厚度约束条件：
0.3≤t≤0.4
在追求极惯性矩的同时，肯定会影响翼型的气动性能。升阻比作为衡量翼型气动性能的重要指标，因为DU35翼型在8°攻角附近其升阻比最大(在步骤四中将8°攻角下的升阻比作为一个目标函数)，所以选择[0°，2°，4°，6°，10°，12°，14°]攻角下的升阻比作为约束条件。
升阻比约束条件：
cl′/cd′≥cl/cd
其中，cl′为优化后的升力系数，cd′为优化后的阻力系数，cl为原翼型的升力系数，cd为原翼型的阻力系数；
采用计算流体力学数值计算方法计算翼型上下翼面的压力及剪切力，从而得到翼型上下表面的升阻力系数，并进一步得到翼型的升阻比，优化后翼型升阻比不低于原翼型。翼型表面压力和剪力分布示意图如图2所示。
步骤四：
通过对步骤2进行分析，随着叶片增长，其半径R必然增加，则叶片扭转刚度减小，增加了叶片的柔性，使得叶片更容易发生颤振。翼型极惯性矩越大，其扭转刚度越高，则抗颤振性能越强，故将翼型的极惯性矩I作为一个设计目标；不同的翼型在不同的攻角下具有可变的升阻比，这里，先定义一个攻角α0，在此攻角下翼型具有最大的初始升阻比，将翼型在α0下的升阻比作为另一个目标函数。
在isight平台中通过采用多目标遗传算法NSGA-II对风力机叶片翼型进行优化设计。利用isight平台对matlab、ug、icem、cfx等软件进行集成运算，建立仿真分析流程，设定和修改设计变量以及约束条件和目标函数，自动进行多次分析循环。可以自动的进行翼型坐标点的更新，自动进行几何建模、网格划分以及cfx的自动化计算。
在本实施例中，将极惯性矩和8°攻角下的升阻比作为目标函数：
maximize＝[f1(x)，f2(x)]
f1(x)＝I
f2(x)＝cl/cd
在isight平台中可以对CAD/CAE软件进行集成，用户在isight平台中自行修改设计变量、约束条件和升阻比，采用多目标遗传因子算法NSGA-II对仿真优化流程自动进行循环。多目标优化的帕累托最优解集如图3所示。
在帕累托最优解集中选择三组翼型与原有翼型比较，如图4a，4b，4c所示；
选择的三组翼型与原有翼型的升阻比进行比较，如图5所示。
步骤五：通过步骤一到步骤四，实现了对翼型抗颤振多目标的设计。此优化流程设计过程中，重点在于如何实现迭代进程中自动更新翼型的几何模型，网格划分，自动求解目标以及求解约束函数。本案例中将翼型几何外形更新程序(matlab文件)、翼型流场几何生成模块(ug文件)、网格划分模块(icem)、数值计算模块(cfx)共同集成到isight平台中，实现数据传递并进行优化迭代进程。同时isight平台可以将设计变量，约束条件，目标函数值的变化历程以可视化的方式直观的显示出来，还可以提供多种数据分析和后处理工具。在isight平台的流程图如图6所示。
实施例二
本实施例以DU35为原始翼型，该翼型最大相对厚度为0.35，位于27.50％翼弦处，该翼型弧度为-1.96％，位于30％的翼弦处。
计算域：
借助ug生成确定计算域的几何模型，流场区域采用半椭圆几何形状，完整的椭圆长轴为100b短轴长度为50b，采用布尔运算从半椭圆流场中减去翼型几何曲线构成的几何体，便可得到翼型绕流流场的计算域的几何模型。如图7所示。
计算网格：
从ug导出parasolid格式文件，然后导入到网格划分软件ICEM中。在翼型边界层区域，考虑到近壁面粘滞作用，分别采用内O-block和外O-block拓扑生成边界层网格，首先，使用内O-block确定翼型几何拓扑，并将内部block删除，其次，使用外O-block确定翼周边界层拓扑区域，以生成边界层网格。近壁层第一层网格高度为0.2mm，确保粘性层第一层至少有一层网格。对网格数目进行划分，网格数为82852。网格拓扑图如图8所示，网格图如图9a、b所示。
计算方法：
在cfx-pre中对流体分析进行前处理，设置参数。设置流体类型为理想气体，收敛精度为1e-6。cfx提供了丰富的后处理选项，通过后处理可以直观地显示出流场计算结果，部分流线图如图10所示。
综上所述，上述实施例的风力机叶片抗颤振翼型簇设计方法，通过在原有的翼型基础上添加翼型扰动函数来修改翼型的几何外形，提高其抗颤振稳定性能并增强其气动性能(升阻比更大)；通过经典抗颤振理论的研究，得出翼型的极惯性矩对其抗颤振性能具有突出的影响，其值可通过高斯-勒让德积分计算得出；并将最大升阻比和极惯性矩作为优化目标，在isight平台中融合翼型扰动函数、高斯-勒让德积分计算、网格划分、流体力学计算、多目标遗传算法来获取优化后的新翼型簇,其抗颤振性能及气动性能均得到了显著的提高，值得被推广使用。
尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。

Claims (10)

1.一种风力机叶片抗颤振翼型簇设计方法，其特征在于，包括以下步骤：
S1：翼型表达
对原有翼型添加形状扰动函数来表达新的翼型；
S2：颤振临界速度影响因素分析
对风机叶片进行敏感性分析，得到颤振临界速度影响因素，提高颤振临界速度影响因素中的扭转刚度以提高抗颤振性能；
S3：选择约束条件
确定翼型的约束条件为厚度约束条件与升阻比约束条件；
S4：确定翼型目标函数
将翼型在攻角α0下的升阻比与翼型的极惯性矩I作为目标函数，通过采用多目标遗传算法对风力机叶片翼型进行优化设计。
2.根据权利要求1所述的一种风力机叶片抗颤振翼型簇设计方法，其特征在于：在所述步骤S1中，所述原有翼型为荷兰DU翼型系列中的DU35。
3.根据权利要求1所述的一种风力机叶片抗颤振翼型簇设计方法，其特征在于：在所述步骤S1中，在原有翼型基础上添加形状扰动函数，保证X轴坐标不变，得到新翼型的型值点坐标：
上翼面翼型的型值点坐标：
yu(i)＝you(i)+Δyu(i)
xu(i)＝xou(i)
下翼面翼型的型值点坐标：
yl(i)＝yol(i)+Δyl(i)
xl(i)＝xol(i)
其中，yu(i)代表新翼型上翼面第i个点的纵坐标，you(i)代表原有翼型上翼面第i个点的纵坐标，Δyu(i)代表上翼面第i个点的形状扰动函数，yl(i)代表新翼型下翼面第i点的纵坐标，yol(i)代表原有翼型上翼面第i个点的横坐标，Δyl(i)代表下翼面第i个点的形状扰动函数；xu(i)代表新翼型上翼面第i个点的横坐标，xou(i)代表原有翼型上翼面第i个点的横坐标，xol(i)代表原有翼型下翼面第i个点的横坐标。
4.根据权利要求3所述的一种风力机叶片抗颤振翼型簇设计方法，其特征在于：上、下翼面形状扰动函数如下：
上翼面形状扰动函数：
Figure FDA0003222689110000021
下翼面形状扰动函数：
Figure FDA0003222689110000022
其中，下标u和l分别代表上翼面和下翼面，i代表翼型型值点索引，k与f(k)分别是控制翼型厚度与弯度分布关键点变量的个数与系数，g为给点矢量，g为给点矢量，f(k)代表Hicks-Henne型函数，ξ和η为指数因子。
5.根据权利要求1所述的一种风力机叶片抗颤振翼型簇设计方法，其特征在于：在所述步骤S2中，翼型的扭转刚度表达式如下：
Figure FDA0003222689110000023
其中，G为剪切模量，R为翼型所在界面距离叶根的半径，I为翼型极惯性矩。
6.根据权利要求5所述的一种风力机叶片抗颤振翼型簇设计方法，其特征在于：所述极惯性矩I的计算过程如下：
S21：计算出翼型的面积和质心坐标点：
翼型的面积：
Figure FDA0003222689110000024
其中，S为翼型的面积；b为翼型弦长，单位标准翼型值弦长取1；u代表横坐标，v代表纵坐标；θ1(u)表示翼型下翼面曲线函数；θ2(u)为上翼面曲线函数；
翼型质心坐标：
Figure FDA0003222689110000025
Figure FDA0003222689110000026
其中，ρ(u，v)代表不同坐标点处的密度，假设叶片是均匀的，则ρ(u，v)为常数；
则上式可进一步推导：
Figure FDA0003222689110000027
Figure FDA0003222689110000028
S22：坐标原点位于翼型的前缘点，求出翼型对x，y方向的惯性矩：
Figure FDA0003222689110000031
Figure FDA0003222689110000032
S23：平移坐标系，使原点与翼型质心重合，通过平行轴定理，求出翼型相对质心轴系的惯性矩：
I1＝Ix-Ay2
I2＝Iy-Ax2
其中：I1代表Y轴相对于质心的惯性矩；
A代表翼型的面积；
y代表Y轴与质心轴的垂直距离；
I2代表X轴相对于质心的惯性矩；
x代表X轴与质心的垂直距离；
则翼型极惯性矩为：
I＝I1+I2。
7.根据权利要求6所述的一种风力机叶片抗颤振翼型簇设计方法，其特征在于：在所述步骤S3中，厚度约束条件为t1≤t≤t2，式中t为翼型厚度，t1、t2为设定值。
8.根据权利要求7所述的一种风力机叶片抗颤振翼型簇设计方法，其特征在于：在所述步骤S3中，优化后的翼型在不同攻角下的升阻比不低于原翼型的升阻比，选择[α1，α2，…，αi]攻角下的升阻比作为约束条件，升阻比约束条件为：
cl′/cd′≥cl/cd
式中cl'为优化后的升力系数，cd'为优化后的阻力系数，cl为原有翼型的升力系数，cd为原有翼型的阻力系数。
9.根据权利要求8所述的一种风力机叶片抗颤振翼型簇设计方法，其特征在于：在所述步骤S4中，目标函数为：
maximize＝[f1(x)，f2(x)]
f1(x)＝I
f2(x)＝cl/cd。
10.根据权利要求1所述的一种风力机叶片抗颤振翼型簇设计方法，其特征在于：在所述步骤S4中，多目标遗传算法为NSGA-II。
CN202110962098.9A 2021-08-20 2021-08-20 一种风力机叶片抗颤振翼型簇设计方法 Active CN113569360B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110962098.9A CN113569360B (zh) 2021-08-20 2021-08-20 一种风力机叶片抗颤振翼型簇设计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110962098.9A CN113569360B (zh) 2021-08-20 2021-08-20 一种风力机叶片抗颤振翼型簇设计方法

Publications (2)

Publication Number Publication Date
CN113569360A true CN113569360A (zh) 2021-10-29
CN113569360B CN113569360B (zh) 2024-03-22

Family

ID=78172308

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110962098.9A Active CN113569360B (zh) 2021-08-20 2021-08-20 一种风力机叶片抗颤振翼型簇设计方法

Country Status (1)

Country Link
CN (1) CN113569360B (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114065423A (zh) * 2021-11-12 2022-02-18 西北工业大学 快速评估航空发动机风扇叶片颤振的方法
CN116776600A (zh) * 2023-06-21 2023-09-19 安徽工业大学 基于自适应代理模型的风力机叶片优化设计方法及系统
GB2623958A (en) * 2022-10-31 2024-05-08 Smallspark Space Systems Design of an aerodynamic component

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104612892A (zh) * 2014-12-30 2015-05-13 中国科学院工程热物理研究所 一种风力机翼型的多学科优化设计方法
CN112507471A (zh) * 2020-10-23 2021-03-16 湖北工业大学 一种大攻角范围条件下垂直轴风力机翼型设计方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104612892A (zh) * 2014-12-30 2015-05-13 中国科学院工程热物理研究所 一种风力机翼型的多学科优化设计方法
CN112507471A (zh) * 2020-10-23 2021-03-16 湖北工业大学 一种大攻角范围条件下垂直轴风力机翼型设计方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
张强等: "基于EI加点准则与代理模型的风力机专用翼型气动与结构优化设计", 中国电机工程学报, 11 August 2021 (2021-08-11), pages 1 - 11 *
王伟等: "风力机叶片三维模型颤振问题", 噪声与振动控制, vol. 33, no. 5, 31 October 2013 (2013-10-31), pages 90 - 96 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114065423A (zh) * 2021-11-12 2022-02-18 西北工业大学 快速评估航空发动机风扇叶片颤振的方法
GB2623958A (en) * 2022-10-31 2024-05-08 Smallspark Space Systems Design of an aerodynamic component
CN116776600A (zh) * 2023-06-21 2023-09-19 安徽工业大学 基于自适应代理模型的风力机叶片优化设计方法及系统
CN116776600B (zh) * 2023-06-21 2024-04-12 安徽工业大学 基于自适应代理模型的风力机叶片优化设计方法及系统

Also Published As

Publication number Publication date
CN113569360B (zh) 2024-03-22

Similar Documents

Publication Publication Date Title
CN113569360B (zh) 一种风力机叶片抗颤振翼型簇设计方法
WO2021077855A1 (zh) 一种直升机旋翼翼型确定方法及系统
CN104612892B (zh) 一种风力机翼型的多学科优化设计方法
CN104408260B (zh) 一种潮流能水轮机叶片翼型设计方法
CN107391891B (zh) 一种基于模型融合方法的大展弦比机翼优化设计方法
Tao et al. Robust optimization for a wing at drag divergence Mach number based on an improved PSO algorithm
WO2020134388A1 (zh) 一种基于随机等几何分析的叶片高刚度设计方法
Lee et al. Performance improvement of horizontal axis wind turbines by aerodynamic shape optimization including aeroealstic deformation
CN112507471B (zh) 一种大攻角范围条件下垂直轴风力机翼型设计方法
CN111400834A (zh) 风力发电机组叶片翼型气动优化设计方法、模型及装置
Song et al. Airfoil optimization using a machine learning-based optimization algorithm
CN114154434A (zh) 一种无尾飞翼布局的多约束精细化气动优化设计方法
Zhaoyun et al. Optimization of tandem blade based on improved particle swarm algorithm
CN116822417B (zh) 一种风力机叶片表面风压分析方法
CN117494611A (zh) 一种基于有粘当地流活塞理论的超声速颤振分析方法
CN111563320B (zh) 一种结构与水弹性一体化螺旋桨设计方法
Martinez et al. Aerodynamic analysis of wind turbine rotor blades
CN112199777A (zh) 一种适用于模化仿生前缘流场特征的方法
Pang et al. Intelligence algorithm for optimization design of the low wind speed airfoil for wind turbine
CN114065423B (zh) 快速评估航空发动机风扇叶片颤振的方法
Pinto et al. On Rotor Aeroacoustic Optimization for Urban Air Mobility
CN113486512A (zh) 一种功能梯度变厚度叶片模型的颤振分析方法
CN113626935A (zh) 一种高巡航效率跨音速月形机翼设计方法
Balasubramanian et al. Novel curvature-based airfoil parameterization for wind turbine application and optimization
Abdulkareem et al. Numerical Investigation of the effect of changing the thickness of airfoils used in wind turbines on the lift to drag ratio

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant