CN105117539A - 风力机叶片模态频率及其双峰高斯分布拟合方法 - Google Patents

Abstract

本发明涉及一种风力机叶片模态频率及其双峰高斯分布拟合方法，包括以下步骤：1)以风力机叶片的实体模型为研究对象，根据瞬时最小势能原理建立结构运动微分方程，利用Block？Lanczos方法求解特征方程；2)通过有限元分析软件ANSYS数值计算不同转速下风力机叶片的前十阶模态频率，得到模态频率随转速变化的两种不同的变化规律，给出相应的变化机理分析；3)对风力机叶片第一阶和第十阶模态频率分别随转速的变化曲线均进行双峰高斯拟合，并开展对比验证。与现有技术相比，本发明具有精确度高等优点。

Description

风力机叶片模态频率及其双峰高斯分布拟合方法

技术领域

本发明涉及一种风力机叶片模拟技术，尤其是涉及一种风力机叶片模态频率及其双峰高斯分布拟合方法。

背景技术

目前风力发电受到各国政府和企业的高度关注，风电装机量稳步上升，而风力机的安全运行是保证风力机稳定发电的基础，叶片作为风力机的关键部件，其性能好坏直接影响风力发电的工作效率、生产成本与安全性。由于其展向长、弦向短、柔性较好以及所受载荷复杂多变的特点，叶片容易发生强烈振动甚至共振，因此叶片振动模态的计算分析十分重要。目前国内外学者利用有限元和实验方法对结构的固有频率和振型等模态参数进行了研究，但对不同转速下叶片模态频率的变化规律尚无深入研究。

发明内容

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种风力机叶片模态频率及其双峰高斯分布拟合方法。

本发明的目的可以通过以下技术方案来实现：

一种风力机叶片模态频率及其双峰高斯分布拟合方法，其特征在于，包括以下步骤：

1)以风力机叶片的实体模型为研究对象，根据瞬时最小势能原理建立结构运动微分方程，利用BlockLanczos方法求解特征方程；

2)通过有限元分析软件ANSYS数值计算不同转速下风力机叶片的前十阶模态频率，得到模态频率随转速变化的两种不同的变化规律，给出相应的变化机理分析；

3)对风力机叶片第一阶和第十阶模态频率分别随转速的变化曲线均进行双峰高斯拟合，并开展对比验证，验证结果表明该拟合方法精确度高。

所述的步骤1)具体为：

101)叶片由有限元离散化处理后，运用瞬时最小势能原理导出以下结构运动微分方程：

$\[M\]\left\{\stackrel{\·\·}{u}\right\}+\[C\]\left\{\stackrel{\·}{u}\right\}+\[K\]\left\{u\right\}=\left\{F\right\}---\left(1\right)$

式中，{u}、及{F}分别表示有限元节点位移、速度、加速度和作用在叶片上外力的矢量；[M]、[C]分别表示质量和阻尼矩阵，刚度矩阵[K]＝[K₀]+[K_a]，[K₀]为叶片小变形的有限元刚度矩阵，[K_a]为动刚度矩阵；

102)在计算旋转叶片的振动模态时，动刚度矩阵表示为

[K_a]＝[K_s]+[K_g](2)

式中，[K_s]为叶片刚体运动与弹性变形耦合而产生的应力刚度矩阵，[K_s]为由离心惯性力导致面内变形所引起的几何刚度矩阵；

103)在式(1)中，令{F}＝{0}，这时叶片处于自由振动状态，在求解叶片结构的固有频率时，不考虑阻尼作用，即[C]＝[0]，因此，式(1)可写为

$\[M\]\left\{\stackrel{\·\·}{u}\right\}+\[K\]\left\{u\right\}=\left\{0\right\}---\left(3\right)$

设代入式(3)后得到叶片结构动力问题的特征值方程，即有

式中，ω_i和分别为第i阶的特征值和特征向量，分别对应叶片结构的第i阶模态频率和模态振型；ω_i＝2πf_i，f_i为自振频率。

所述的步骤2)具体为：

201)对实体模型进行精细网格划分，通过网格无关性验证后，得到了有限元模型，改变叶片的转速，计算出叶片在不同转速下的前十阶模态频率；

202)随着转速的增加，旋转叶片的第一阶到第六阶模态频率均呈现非线性递减，其中第一阶模态的减少趋势具有代表性；

203)第七阶至第十阶模态频率随转速增加而增加，其中第十阶模态的变化趋势最为明显。

随着转速增加，叶片刚体运动与变形的耦合作用会导致模态频率的减少，而离心惯性力作用则会导致模态频率的增加。

所述的步骤3)具体为：

301)设叶片第i阶的模态频率f_i随转速ω变化满足双峰高斯分布

$f_i=a_1\exp (-{\left(\frac{\mathrm{\ω}-b_1}{c_1}\right)}^2)+a_2\exp (-{\left(\frac{\mathrm{\ω}-b_2}{c_2}\right)}^2)---\left(5\right)$

式中，a₁、b₁、c₁、a₂、b₂、c₂均为双峰高斯公式的待估参数，

302)根据实测数据得到叶片第一阶、第十阶模态频率f₁、f₁₀随转速ω变化曲线的拟合公式，分别写为

$f_1=2.63\exp \left(-{\left(\frac{\mathrm{\ω}+75.11}{55.09}\right)}^2\right)-2.014\×{10}^{-2}\exp \left(-{\left(\frac{\mathrm{\ω}-11.49}{6.462}\right)}^2\right)---\left(6\right)$

$f_{10}=-4.896\exp (-{\left(\frac{\mathrm{\ω}-23.5}{6.106}\right)}^2)+4.254\×{10}^7\exp (-{\left(\frac{\mathrm{\ω}-456.0}{116.4}\right)}^2)---\left(7\right).$

与现有技术相比，本发明结合国内外风力机叶片相关理论研究和实际应用工况，建立了叶片的实体模型，计算了不同转速下叶片的前十阶模态频率，揭示了模态频率随转速的变化规律，给出了相应的理论分析，可为风力机叶片的优化设计和进一步的振动特性分析提供了技术参考，对风机叶片的安全运行具有重要的参考价值。

附图说明

图1为叶片的实体模型及有限元网格剖分示意图；

图2为叶片的约束方式示意图；

图3为叶片第一阶模态频率随转速的变化曲线图；

图4为叶片第十阶模态频率随转速的变化曲线图；

图5为叶片第一阶模态频率随转速变化的拟合曲线及对比验证图；

图6为叶片第十阶模态频率随转速变化的拟合曲线及对比验证图。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。

实施例

一种风力机叶片模态频率及其双峰高斯分布拟合方法，包括以下步骤：

1)以风力机叶片的实体模型为研究对象，根据瞬时最小势能原理建立结构运动微分方程，利用BlockLanczos方法求解特征方程；

2)通过有限元分析软件ANSYS数值计算不同转速下风力机叶片的前十阶模态频率，得到模态频率随转速变化的两种不同的变化规律，给出相应的变化机理分析；

3)对风力机叶片第一阶和第十阶模态频率分别随转速的变化曲线均进行双峰高斯拟合，并开展对比验证，验证结果表明该拟合方法精确度高。

1、风力机叶片的动力特征方程

叶片由有限元离散化处理后，运用瞬时最小势能原理可导出以下结构运动微分方程

$\[M\]\left\{\stackrel{\·\·}{u}\right\}+\[C\]\left\{\stackrel{\·}{u}\right\}+\[K\]\left\{u\right\}=\left\{F\right\}---\left(1\right)$

式中，{u}、及{F}分别表示有限元节点位移、速度、加速度和作用在叶片上外力的矢量；[M]、[C]分别表示质量和阻尼矩阵，刚度矩阵[K]＝[K₀]+[K_a]，[K₀]为叶片小变形的有限元刚度矩阵，[K_a]为动刚度矩阵。

风力机叶片绕转轴以一定角速度大范围转动时，离心惯性力影响了叶片的变形，此时叶片的转动与变形之间相互耦合而导致叶片刚度变化，出现“动力刚化”现象。在计算旋转叶片的振动模态时，动刚度矩阵可表示为

[K_a]＝[K_s]+[K_g](2)

式中，[K_s]为叶片刚体运动与弹性变形耦合而产生的应力刚度矩阵，[K_s]为由离心惯性力导致面内变形所引起的几何刚度矩阵。

在式(1)中，令{F}＝{0}，这时叶片处于自由振动状态。此外，在求解叶片结构的固有频率时，一般不考虑阻尼作用，即[C]＝[0]。因此，式(1)可写为

$\[M\]\left\{\stackrel{\·\·}{u}\right\}+\[K\]\left\{u\right\}=\left\{0\right\}---\left(3\right)$

设代入式(3)后可得到叶片结构动力问题的特征值方程，即有

式中，ω_i和分别为第i阶的特征值和特征向量，分别对应叶片结构的第i阶模态频率和模态振型；ω_i＝2πf_i，f_i为自振频率。

2、应用软件的求解过程

运用ANSYS软件对风力机叶片的模态进行求解，分三步进行。

1)实体建模

采用目前装机容量为5MW的风力机叶片作为研究对象,将翼型坐标进行变换，运用UG建立叶素和曲面，并通过曲线组功能建立了叶片实体模型，如图1所示。

2)设置实体模型的约束方式

将叶片模型导入Geometry模块中，并在Model模块对实际运行工况进行设置。对轮毂施加CylindricalSupport约束，使得叶片只有RZ的自由度，如图2所示。

3)不同转速下的前十阶模态频率

在Inertial中选择RotationalVelocity，以轮毂中心为转轴，设置叶片的转速，在AnalysisSettings中设置模态阶数为10，见图2。对实体模型进行精细网格划分，通过网格无关性验证后，得到了有限元模型，如图1所示。改变叶片的转速，由此可计算出叶片在不同转速下的前十阶模态频率。

3、不同转速下叶片模态频率的计算与分析

设叶片在展向和弦向的弹性模量相同，材料为各向同性。玻璃钢叶片的材料参数见表1。

表1

密度	泊松比	弹性模量
			1700kg/m3	0.17	17.6GPa

根据ANSYS软件的求解方法，经相关的网格划分和条件设置，转速分别为3.5rpm、6.9rpm、9.5rpm、12.1rpm、15.0rpm、20.0rpm，计算得到了六种不同转速下叶片的前十阶模态频率(Hz)，如表2所示。

表2

分析表2后不难得知：随着转速的增加，旋转叶片的第一阶到第六阶模态频率均呈现非线性递减，其中第一阶模态的减少趋势具有一定代表性，如图3所示，这一变化主要来源于叶片刚体运动与弹性变形耦合所带来的刚度弱化作用，此时应力刚度的减少占主导地位。相反，第七阶至第十阶模态频率随转速增加而增加，其中第十阶模态的变化趋势最为明显，如图4所示，这是由于离心惯性力所引起的面内变形使叶片刚度增大的缘故，此时几何刚度的增加占主导地位，从而导致了模态频率呈缓慢上升的趋势。

由此可见，随着转速增加，叶片刚体运动与变形的耦合作用会导致模态频率的减少，而离心惯性力作用则会导致模态频率的增加。因此，.要考虑在不同工况下两者的联合效应，分析哪一因素起主导作用，这样才能准确反映旋转叶片的振动特性

4、利用高斯函数法进行转速-模态频率曲线拟合

下面结合有限元计算的数据点给出转速-模态频率曲线的拟合方法，利用实例并对拟合曲线的精准度进行验证。

设叶片第i阶的模态频率f_i随转速ω变化满足双峰高斯分布

$f_i=a_1\exp (-{\left(\frac{\mathrm{\ω}-b_1}{c_1}\right)}^2)+a_2\exp (-{\left(\frac{\mathrm{\ω}-b_2}{c_2}\right)}^2)---\left(5\right)$

式中，a₁、b₁、c₁、a₂、b₂、c₂均为双峰高斯公式的待估参数。分别结合表2中六种不同叶片转速下的第一阶、第十阶模态频率，即分别为图3、4中的实心点数据，利用上述方法，可得到叶片第一阶、第十阶模态频率f₁、f₁₀随转速ω变化曲线的拟合公式，分别写为

$f_1=2.63\exp \left(-{\left(\frac{\mathrm{\ω}+75.11}{55.09}\right)}^2\right)-2.014\×{10}^{-2}\exp \left(-{\left(\frac{\mathrm{\ω}-11.49}{6.462}\right)}^2\right)---\left(6\right)$

$f_{10}=-4.896\exp (-{\left(\frac{\mathrm{\ω}-23.5}{6.106}\right)}^2)+4.254\×{10}^7\exp (-{\left(\frac{\mathrm{\ω}-456.0}{116.4}\right)}^2)---\left(7\right)$

图5和图6分别给出了第一阶、第十阶有限元计算的数据点及其双峰高斯拟合曲线，可以看出，数据点基本分布在拟合曲线上，吻合较好。

实例：

当叶片转速ω＝17.5rpm时，通过本发明的计算步骤可得到叶片第一阶和十阶模态频率的有限元计算值(分别对应于图5和图6中的空心点)，分别为f₁ ^*＝0.1458Hz和为例。将ω＝17.5rpm分别代入式(6)和(7)，得到叶片的第一阶、第十阶模态频率的拟合值分别为f₁＝0.1473Hz和f₁₀＝27.3390Hz。以有限元计算值为基准，可得到第一阶、第十阶的双峰高斯分布拟合的相对误差δ_f1和δ_f10，分别为

${\mathrm{\δ}}_{f_1}=\left|\frac{f_1-{f_1}^{*}}{{f_1}^{*}}\right|=\left|\frac{0.1473-0.1458}{0.1458}\right|\≈1.029\%---\left(8\right)$

${\mathrm{\δ}}_{f_{10}}=\left|\frac{f_{10}-{f_{10}}^{*}}{{f_{10}}^{*}}\right|=\left|\frac{27.3390-27.3110}{27.3110}\right|\≈0.103\%---\left(9\right)$

由此可见，本发明的风机叶片模态频率计算方法以及双峰高斯分布拟合公式的确定方法切实可行，对于风力机模态频率f随转速ω变化曲线，不管是单调上升还是单调下降，双峰高斯分布拟合均具有很高的拟合精度。

Claims (5)

1.一种风力机叶片模态频率及其双峰高斯分布拟合方法，其特征在于，包括以下步骤：

1)以风力机叶片的实体模型为研究对象，根据瞬时最小势能原理建立结构运动微分方程，利用BlockLanczos方法求解特征方程；

2)通过有限元分析软件ANSYS数值计算不同转速下风力机叶片的前十阶模态频率，得到模态频率随转速变化的两种不同的变化规律，给出相应的变化机理分析；

3)对风力机叶片第一阶和第十阶模态频率分别随转速的变化曲线均进行双峰高斯拟合，并开展对比验证。

2.根据权利要求1所述的一种风力机叶片模态频率及其双峰高斯分布拟合方法，其特征在于，所述的步骤1)具体为：

101)叶片由有限元离散化处理后，运用瞬时最小势能原理导出以下结构运动微分方程：

$\[M\]\left\{\stackrel{\·\·}{u}\right\}+\[C\]\left\{\stackrel{\·}{u}\right\}+\[K\]\left\{u\right\}=\left\{F\right\}---\left(1\right)$

式中，{u}、及{F}分别表示有限元节点位移、速度、加速度和作用在叶片上外力的矢量；[M]、[C]分别表示质量和阻尼矩阵，刚度矩阵[K]＝[K₀]+[K_a]，[K₀]为叶片小变形的有限元刚度矩阵，[K_a]为动刚度矩阵；

102)在计算旋转叶片的振动模态时，动刚度矩阵表示为

[K_a]＝[K_s]+[K_g](2)

式中，[K_s]为叶片刚体运动与弹性变形耦合而产生的应力刚度矩阵，[K_s]为由离心惯性力导致面内变形所引起的几何刚度矩阵；

103)在式(1)中，令{F}＝{0}，这时叶片处于自由振动状态，在求解叶片结构的固有频率时，不考虑阻尼作用，即[C]＝[0]，因此，式(1)可写为

$\[M\]\left\{\stackrel{\·\·}{u}\right\}+\[K\]\left\{u\right\}=\left\{0\right\}---\left(3\right)$

设代入式(3)后得到叶片结构动力问题的特征值方程，即有

式中，ω_i和分别为第i阶的特征值和特征向量，分别对应叶片结构的第i阶模态频率和模态振型；ω_i＝2πf_i，f_i为自振频率。

3.根据权利要求1所述的一种风力机叶片模态频率及其双峰高斯分布拟合方法，其特征在于，所述的步骤2)具体为：

201)对实体模型进行精细网格划分，通过网格无关性验证后，得到了有限元模型，改变叶片的转速，计算出叶片在不同转速下的前十阶模态频率；

202)随着转速的增加，旋转叶片的第一阶到第六阶模态频率均呈现非线性递减，其中第一阶模态的减少趋势具有代表性；

203)第七阶至第十阶模态频率随转速增加而增加，其中第十阶模态的变化趋势最为明显。

4.根据权利要求3所述的一种风力机叶片模态频率及其双峰高斯分布拟合方法，其特征在于，随着转速增加，叶片刚体运动与变形的耦合作用会导致模态频率的减少，而离心惯性力作用则会导致模态频率的增加。

5.根据权利要求1所述的一种风力机叶片模态频率及其双峰高斯分布拟合方法，其特征在于，所述的步骤3)具体为：

301)设叶片第i阶的模态频率f_i随转速ω变化满足双峰高斯分布

$f_i=a_1\exp (-{\left(\frac{\mathrm{\ω}-b_1}{c_1}\right)}^2)+a_2\exp (-{\left(\frac{\mathrm{\ω}-b_2}{c_2}\right)}^2)---\left(5\right)$

式中，a₁、b₁、c₁、a₂、b₂、c₂均为双峰高斯公式的待估参数，

302)根据实测数据得到叶片第一阶、第十阶模态频率f₁、f₁₀随转速ω变化曲线的拟合公式，分别写为

$f_1=2.63\exp \left(-{\left(\frac{\mathrm{\ω}+75.11}{55.09}\right)}^2\right)-2.014\×{10}^{-2}\exp \left(-{\left(\frac{\mathrm{\ω}-11.49}{6.462}\right)}^2\right)---\left(6\right)$

$f_{10}=-4.896\exp (-{\left(\frac{\mathrm{\ω}-23.5}{6.106}\right)}^2)+4.254\×{10}^7\exp (-{\left(\frac{\mathrm{\ω}-456.0}{116.4}\right)}^2)---\left(7\right).$