0 引　言

全球气候变暖的问题加重，使得清洁能源的开发利用成为当前能源行业发展的关注热点^[1]，风能发电已然成为能源行业重要的战略部署。叶片作为风力机能量转化的重要部件，其结构特性是保证风机稳定运行和提高风机能量转换率的关键，所以在设计初期对相应的叶片结构参数进行校核必不可少。刚度对叶片挠度安全性起着决定性作用，其大小将直接影响风机运行的可靠性^[2]，由于叶片结构与材质的不规则性，刚度的获取多以后期叶片实验为主^[3]。因此，研究叶片的刚度计算方法对叶片设计初期的安全校核以及叶片优化前后的刚度评估有着重要意义^[4]。

目前国内外对于叶片刚度的研究现状，将刚度与叶片的气动弹性相结合，探究不同程度的刚度对叶片运行的影响^[5]，其中在叶片刚度计算方面，多以简化叶片的数学模型或是简化叶片结构模型进行计算，多数为叶片截面刚度分布计算。例如通过研究叶片结构发现主梁帽刚度占据总刚度的90%以上，从而将叶片简化为工字梁结构进行刚度计算^[6]。在叶片数学模型建模上，利用有限差分法和有限元法的数学方法，将叶片等效为变截面悬臂梁，再利用壳单元，简化叶片结构建立有限元模型，从而获得叶片刚度分布，为叶片刚度分析提供参考^[7]。因为刚度直接影响挠度，所以以挠度作为中间参考值，用挠度参数代替刚度参数体现叶片刚度大小^[8]。在截面刚度计算方面，利用截面刚度沿叶片展向分布曲线校核刚度^[9]，使用 Matlab优化工具箱计算出叶片各截面的参数，再结合材料力学理论计算挠度沿叶高的分布曲线，从而完成刚度校核^[10]。

通过上述研究可知，叶片的刚度计算多以截面刚度分布曲线或是以挠度作为中间值进而分析叶片整体刚度，但刚度参考参数不够宏观直接。针对此现状问题，参考材料力学中变截面构件挠度计算的等效刚度法^[11-12]，根据风机叶片安装的边界条件以及叶片本身变截面的几何特征，将其简化成变截面悬臂梁，用等效刚度法求得的等效刚度值体现叶片刚度，相对刚度分布曲线而言，单一的刚度参数值，可以更加宏观直接地体现叶片刚度。但原始等效刚度法的线性近似用于叶片存在着误差较大的问题。针对该问题，以叶片的截面数据特征为基础对原始的等效刚度法进行优化，通过模拟实验得到叶尖位移值验证优化的可靠性^[13]，用优化后的算法计算得到的等效刚度值代表叶片整体刚度，为叶片设计初期的刚度提供参考值。

1 变截面悬臂梁的等效刚度法优化 1.1 变截面悬臂梁等效刚度法的基础原理

等效刚度法是通过等效系统法将变刚度的梁用等刚度的梁来代替，其中等效系统的每一个挠曲线等于原来变刚度梁的挠曲线。若由一已知的变刚度梁得到了一个等效系统，就可用等效系统对原结构进行计算。

变截面构件的弹性曲线一般微分方程为

积分二次得到变截面挠度曲线：

则变刚度可表示为某一截面刚度与函数的乘积，如下式所示：

等刚度挠度曲线：

由于式（2）、式（4）相同，且边界条件相同，所以

则有

即

通过上述讨论说明，一个变刚度的梁总是存在一个等刚度的等效系统，但这里的等刚度并不等于原结构的刚度，因为是通过改变弯矩从而达到刚度等效，也就是说，若$ f(x) = 1 $，则$ {E_x}{I_x} = {E_1}{I_1} $，即等效系统的刚度等于原结构的刚度。

1.2 基于叶片等效刚度法的计算优化

根据等效系统原理，变刚度与等效刚度的悬臂梁的叶尖挠度相同，可建立公式：

其中$l$为固定端到自由端的长度，${I_x}$为距离坐标系原点（固定端）$x$时的截面惯性矩，$x$为在$x$轴上与原点的距离，${I_{{\rm{eq}}}}$为等效刚度系统中的等效惯性矩。

由式（7）即可推出等效刚度：

在经典的等效刚度法计算中，主要用于变刚度构件截面惯性矩形成的函数积分求解，数学计算坐标示意图，如图1所示。

原始的计算方法是通过取变截面悬臂梁的固定端和自由端的惯性矩值，将其近似成简单的线性函数，如式（9）所示，从而完成积分求值。然而叶片沿展向方向截面惯性矩变化程度大，若仍取叶片的固定端和自由端的惯性矩值进行线性近似会产生较大误差，并且不能反映叶片中间段的特殊翼型对叶片整体带来的刚度影响。因此利用叶片建模截面数据的特点以及原始等效积分法存在线性近似积分法误差过大的缺点，通过这两个方面进行基于叶片等效刚度法的计算优化。

式中：$ {I}_{自由端} $——变截面悬臂梁自由端的截面惯性矩值，m⁴；

$ {I}_{固定端} $——变截面悬臂梁固定端的截面惯性矩值，m⁴；

$I(x)$——截面惯性矩函数；

$x$——距自由端距离。

根据叶片截面数据具有实验样本的数据特征并且呈现非线性，在积分法上可以选取梯形积分和辛普森积分，其次在区间分割相同的条件下，辛普森积分公式比梯形积分公式的计算精度高^[13]，但辛普森积分法仅适用于等分成偶数个小区间的区间上，所以根据叶片建模存在部分截面连续等间距和部分截面非等间距的特点，采用辛普森积分法和梯形积分法分别对等间距截面区间和非等间距截面进行分段积分求和，提高计算精度。

假设对函数$ {y_i} $在区间$(0,l)$积分，其中$(0,p)$为单独的非等分区间，$(p,l)$为$ 2n $个等分区间，每个区间长度为$ h $，结合上述区间特征对等分区间采用辛普森公式，非等分区间采用梯形积分，最后进行分段积分求和，其中$ {y_i} = \dfrac{{{x^3}}}{{{I_x}}} $推出积分求解公式，可得到优化后的积分求解公式如式（11）所示，再将式（11）代入式（8）即可优化原始等效刚度法的计算进而求解。

2 计算案例

以小型风机的某S翼型叶片的摆振刚度计算为例，简述叶片建模初期的模型数据和特征以及如何应用优化后的等效刚度法对叶片取值的刚度计算。

2.1 叶片模型及材料数据

以某S翼型叶片为例，叶片结构总长700 mm，由10个特征翼型面组成，每个翼型面之间通过放样将其连接，其中0号和1号翼型面间距为35 mm，1号至9号翼型面相邻的面为等间距70 mm，叶根部位长度105 mm，叶片模型如图2所示，叶片材料参数如表1所示。

2.2 基于叶片的等效刚度法计算

由于叶片安装方式以及叶片主体的特点，将叶片近似成变截面悬臂梁模型，基于叶片的安装方式近似将10号截面作为固定端截面，0号截面作为自由端截面，以固定端为$ x $轴坐标原点，$ x $轴为沿叶片展向距离固定端的截面序号，$ y $轴为截面惯性矩可建立坐标如图3（a）所示。

由于等效刚度法的计算是基于挠度系统能量守恒推理得出，但在等效系统中采用共轭梁法计算挠度，原始的叶片坐标系会因为共轭梁系统建立从而发生对调，导致原来的固定端变成自由端，自由端变成固定端，变化后的坐标原点首尾互换，如图3（b）所示，导致截面的排列是从叶尖到叶根进行$x$的取值，通过Solidworks的截面属性功能计算得到叶片建模时各个截面的惯性矩。叶片的截面数据如表2所示。

根据叶片建模的截面间距的分布特点，0号至1号截面间距为35 mm和9号至10号截面间距为64.5 mm采用梯形积分法，1号至9号截面为8等分间距为70 mm采用辛普森积分法，加以分段积分求和，基于式（11）代入式（8）得到等效刚度的优化后叶片刚度的求解公式，如下式所示：

式中：$ E $——叶片的杨氏模量，N/mm²；

${I_{{\rm{eq}}}}$——等效惯性矩，mm^４；

$ H $——固定端到自由端的长度，mm；

$ {x_i} $——$ i $号截面距离坐标系原点（固定端）的长度，mm；

$ {I_i} $——$ i $号截面的截面惯性矩，mm^４。

计算求得该S翼型叶片的等效刚度156.19 $ {\text{N/}}{{\text{m}}^{\text{2}}} $。

3 模拟试验 3.1 模型建立

叶片模型为某S翼型小型水平轴风力机叶片，长700 mm，通过10个特征翼型坐标在CAD中绘制翼型截面，将各翼型截面在Solidworks中结合截面间隔和法兰连接盘放样生成单只叶片，将三维模型导入ANSYS中，以$ x $轴和$ {\textit{z}} $轴形成$ xo{\textit{z}} $为平面平放叶片，如图4所示。

3.2 网格划分

网格采用ANSYS Workbench中的mesh模块进行划分，由于风力机叶片结构的不规则和叶片表面曲度高，结构化网格对复杂外形的贴体网格生成较难，非结构化网格对于不规则的结构具有更好的适应性，所以采用非结构化网格对叶片实体的网格进行划分，非结构化的网格单元尺寸选取5 mm，网格数量与质量如表3所示，网格划分密度如图5所示，网格质量较为理想。

3.3 计算类型选取和边界条件设置

叶片为三维实体结构，具有$ xy{\textit{z}} $ 3个方向的自由度，并具有一定弹性与位移，模拟主要计算叶片的结构位移所以选择Solid186单元类型计算。

叶片等效成悬臂梁，所以边界条件的选取与叶片的安装方式相结合。由于法兰连接，所以选取圆盘处作为固定端，如图6所示。添加重力加速度，垂直于叶片平面施加均布载荷，在均布载荷的取值上，根据叶片的常规风速范围，取6 m/s，8 m/s，10 m/s，12 m/s这4个风速在20 ℃所对应的22.5 Pa，40 Pa，62.5 Pa，90 Pa风压值，风压计算公式如式（12）所示，计算获得的风压值作为均布载荷的压强值。

式中：$ P $——风压值，Pa；

$ \rho $——空气密度，1.205 kg/m³（20 ℃时）；

$ v $——风速，m/s。

3.4 模拟试验结果

用模拟试验的摆振方向位移值作为中间参数，验证摆振刚度计算值的可靠性，所以选取方向性位移计算板块，对叶片载荷方向下进行位移求解，可得叶片位移云图，如图7所示。

悬臂梁的自由端挠度为最大挠度，在叶片上对应的最大挠度就是叶尖位移，模拟结果可得到22.5、40、62.5、90 Pa下均布载荷的叶尖位移值如表4所示。

3.5 试验值与理论值的结果分析

以模拟试验的叶尖位移为中间值作为参考标准，验证优化后等效刚度法计算结果是否达到优化目的。在不同风压值下，用原等效刚度法得到的刚度值和优化后计算得到的等效刚度值代入式（7），即可获得优化前后相应的叶尖位移值如表5所示。将优化前后的等效刚度值与模拟试验的叶尖位移值对比可知，优化后的等效刚度法减小了该方法应用于叶片计算的误差，实现了量级式的误差减小，证明等效刚度法优化后应用于叶片的适用性。

优化后的等效刚度法计算的位移值与模拟试验的误差值如表5所示。由表中误差值可知，随着载荷值的增加误差值随之增大，误差最大为43.64%，在低载荷值22.5 Pa下的误差值达到最小为9.68%。

在计算方法和数据正确且不考虑数学计算本身存在的误差情况下，对误差的分析回归到叶片的数学建模问题，在对叶片采取梁的数学计算模型中包含相应的梁的假设，由于叶片在叶片的模拟试验时模型并未按照“平截面假设”后的计算模型进行计算，存在基于叶片的展向以及纵向两个方向的剪切变形量，由于纵向剪切变形的存在有效减弱了外部载荷的作用效果，在梁的计算中并未考虑剪切形变，并且剪切变形量与力呈线性关系，所以外加载荷增大剪切变形量越大的情况下，计算误差会随载荷增加而增大。由于梁作为特殊的平板模型，将叶片简化成梁在不考虑转角自由度的剪切变形时计算弯曲挠度是可行的。

所以优化后的等效刚度法估算叶片的弯曲刚度在不考虑叶片的剪切变形的情况下，对叶片的刚度评估是可靠的，可为叶片设计或优化提供新的刚度参数参考，为叶片的刚度计算提供新方法，方法简单，便于工程应用。

4 结束语

1）提出了叶片新的刚度计算方法。基于材料力学的等效刚度法，根据叶片截面数据特征，改善原始等效刚度法的数学计算方法，使其在叶片的计算上减小误差，适用于叶片刚度计算。

2）验证了优化等效刚度法的可行性。通过均布载荷的模拟试验，将模拟实验计算的叶尖位移作为中间参考参数，在22.5 Pa的载荷下误差值达到9.68%，说明在减小剪切变形的影响下，误差在工程误差的允许范围内，即该方法计算得到的刚度值可为评估叶片刚度提供参考。