风力机叶片动力学响应分析

势瓦＋丁鲁＋＝瓦Ｐｙ２一一硝Ｍ虿雾
截面上的应力与弯矩的关系为：
Ｉ Ｙｌ仃ｌｄＡｌ＋ｌ Ｙｌ盯２ｄＡ２＋…＝Ｍ。
ｏ ４ｌ
ｏ＾２
可得：
等．ｆ．。，，ｊｄＡ。＋等』屯ｙ２ｄＡｚ＋…＝Ｍ。
由此得挠度与弯矩的关系：
肘。＝（Ｅｌ，，ｌ＋如，，２＋…）罟｛
代入振动方程得：
习习ａ（２（Ｅ州ｔｄ，＋Ｅ２ｄ如＋＋如．．＋．．）．＋础。型３）２＋＝ｐ础ｙ－ａ型ｔ２ 等
模型，但Ｂｅｄｄｏｅｓ—Ｌｅｉｓｈｍａｎ模型考虑非定常效应时使 用二阶微分方程来描述，且需通过叶型试验要确定 的参数较多，因此比较复杂。这里介绍Ｋｒｅｎｋ模型，
Ｋｒｅｎｋ模型是在考虑风力机叶片实际状况的基础上
万方数据
太阳能学报
３１卷
对Ｂｅｄｄｏｅｓ．Ｌｅｉｓｈｍａｎ模型作了适当的简化。 在Ｋｒｅｎｋ模型中用ｃ，。表示定常条件下不考虑
本文采用Ｍａｒｃ有限元软件实现风力机叶片动 力学响应计算。选用ＦＯＲＣＤＴ载荷子程序施加气动 载倚。根据Ｍａｒｃ子程序的调用特点，在计算的每一 个增量步中，对于每一个将载荷定义成ＦＯＲＣＤＴ子 程序的节点，该子程序都要被调用一次。在静力分 析中，位移向量被传递到子程序中；在动态分析中， 速度和加速度量也被传递过来［２ Ｊ。由此可见，只要 将计算叶片气动载荷的方法编成ＦＯＲＣＤＴ子程序， 就可实现在计算中每个增量步互相交换边界条件的 目的。
肼。＝（Ｇｌ‘＋Ｇ２缸＋…）关
由此得振动方程：
７１，订ａ２０＝未（Ｇ。［ｐｌ＋Ｇ２”…）箬饥
式中，Ｇ；——截面各铺层材料的面内剪切模量； ，。；——截面不同材料的极惯性矩。且有：
，，＝∑，ｐｉ
根据以上的推导，梁单元使用各向同性材料， Ｅ、Ｇ分别由下式确定：
，；：∑，“ＥＶ：Ｉ竿Ｆ
Ｔ可＿ ∑，，ｉＥ
∑，，ｉＧ；
１ 叶片动力学响应计算基本思路
在前述的三类方法中，第一类方法由于自身假 设的限制，应用不多，第三类方法由于理论上的难 度，亟待进一步的发展，而第二类方法得到了广泛应 用，该方法在每个时问步交换流体域和固体域的信 息，这样，一个域求解的结果就被当作边界条件施加 到另一个域中，即在每一步将计算得到的流体力施 加到叶片上，得到叶片的新佗置，而叶片的新位置又 被作为新的流固边界去计算下一步的流体力。
是将叶片沿展向分成许多微段，这些微段称为叶素。 假设每个叶素上的流动相互之间没有干扰，即将叶 素看成二维叶型。
风通过风轮后产生漩涡，因而有轴向诱导速度 和环向诱导速度，分别用轴向诱导因子口、环向诱导 因子口’来表示，轴向诱导速度与来流速度相反，环 向诱导速度与叶片旋转方向相反，由此可以得到考 虑诱导速度后气流相对于叶素的速度的大小。由相 对气流速度、准静态升力、阻力和力矩曲线即可求得 叶素上的力和力矩，将作用在每个叶素上的力和力 矩沿叶片展向积分，就可以求得作用在叶片上的力 和力矩。
关键词：风力机叶片；气动弹性；气动阻尼系数
中图分类号：ＴＫ８３
文献标识码：Ａ
０引 言
随着风力机尺寸的增加，一方面叶片的叶尖切 向速度小断提高，另一方面考虑到重鼍问题，叶片将 越来越细长，在这种趋势下，研究风力机叶片在风载 荷下的动力学响应，从而判断叶片的运动稳定性就 显得越来越霞要。由于叶片处于气流中，流场的流 动对叶片产生的作用力会使叶片产生变形，而叶片 的变形会影响到流场流动，若要较准确的得到叶片 的动力学响应，就必须考虑叶片与流场的这种相互 影响作用，即必须考虑叶片的气动弹性响应问题。
的大小，选用修正的动量叶素理论和梁单元作为流体和结构模型，并将动态失速模型与动虽叶素理论结合计算气
动力；以Ｍａｒｃ通用有限元计算程序为基础平台，将气动力计算程序写入子程序中，并根据子程序的调用特点，实现了 在每个增量步中气动和结构边界条件的相互交换。最后给出计算算例，得到的计算结果与实际使用情况相符合。
分离时的升力曲线；Ｃ。．表示定常条件下考虑分离时 的升力曲线；下标加ｄ表示对应的非定常条件下的 数值；口为攻角。定义分离系数．厂，物理意义如图３ 所示，它表示分离点的位置。
３计算结果及分析
以上面提到的某型７５０ｋＷ风力机叶片为例，计 算了叶片在风速为１０～２５ｍ／ｓ时的响应，图４～图６ 仅绘出２５ｍ／ｓ风速下叶尖节点３个方向的响应历 程，另外几种情况下叶尖节点３个方向的响应历程 图形与此类似。
性基本一致。 ２．３耦合计算过程
ＦＯＲＣＤＴ子程序计算流程如下：在每一个增量 步，Ｍａｒｃ主程序将节点号及弯曲、摆振和扭转方向位 移传给ＦＯＲＣＤＴ载荷子程序，子程序根据节点号得
到截面翼型类型以及该点的初始展向位置，初始安 装角及弦长，根据挥舞、摆振和扭转方向位移修正展 向位置及安装角，并计算攻角，由攻角插值获得升力 系数、阻力系数和力矩系数，根据修正的动量叶素理
万方数据
９１８
太 阳 能 学报
３ｌ卷
这是因为： １）风力机叶片由多层玻璃布铺敷而成，铺层十
分复杂，首先是层数多，最多处有３０多层，其次每一 层的铺敷范围极不规律，只有极个别的层是全型面 铺设的，图ｌ为部分铺层的示意图【３】。如果选用三 维模型，为描述这种复杂的铺层，就需要尽量多的结 构单元，而一般流体域网格要远大于固体域，因此整 个模型规模很大，导致计算时间不可接受；
由于修正的动量叶素理论是比较常用的方法， 限于篇幅这里不再详细叙述，具体过程和用到的公 式可参见文献［４］。
２．２结构模型［５’６】 与动量叶素理论相配套，叶片结构采用梁单元
模拟。由图ｌ可知叶片铺层层数多且形式复杂，无 论沿弦向还是沿展向材料都不存在均匀性，将叶片 简化成梁后，每个实际截面变成了一个点，梁单元的 材料属性如何给出，才能保证梁单元模型的动力学 特性与真实叶片的动力学特性一致，是需要考虑的 问题。在分析叶片的振动时，根据受力情况，叶片主 要有３个方向的运动，即挥舞振动（Ｚ向，坐标参考 图１），摆振振动（Ｙ向）和扭转振动（绕ｘ轴）。Ｙ向 弯曲振动微分方程为：
编号为ｌ、２、３……１６，其中１、２节点之间是圆柱段， 之后各段选用不同的翼型。该叶片的额定工作风速 为１４ＩＩｌ／ｓ。
根据稳态升阻力系数曲线，使用上述模型迭代
计算得到叶片各截面在额定工作风速附近不同风速
下的稳态攻角，如图２所示。
６０
３
４
５０
ｆ４０
＼
援３０ 铎
２０
０１２３４
１０
节摹节符宵节 点赛点 点 ５６７８９１ｌＩｌ ５６
目前求解气弹响应的方法，根据出现的先后顺 序可以分为三大类：第一类是将流体和固体区域解 耦，将流体对结构的作用折算到结构动力学方程中 去，或者反之。这类方法为了能实现将一个区域的 作用效果用另一个区域的参量表示出来，做了许多 假设和限制。第二类方法将流体区域和结构区域的 方程单独求解，但是在每个时间步交换边界条件。 这类方法中一个区域对于另一个区域的响应是滞后 的，只是时间步越小这种滞后的影响越小。第三类 方法用一个统一的方程描述流体域和结构域，结构 和流体运动在每个时间步用同一个求解器同时求 解。由于流体与结构的刚度相差很大，因此这类方 法实现起来比较困难，有待进一步的发展…。
１０ １２ １４ １６ １８ ２０ ２２ ２４ ２６ 风速／ｍ·８－‘
图２各截面攻角随风速变化 ｎｇ．２ Ａｔｔａｃｋ锄ｇｌｅ ｃｈ粕ｇ嗍ｗｉｔｈ ｗｉｎｄ ｓｐｏｅｄ 表１中列出了各截面翼型升力系数直线段结束 对应的攻角以及最大升力系数位置对应的攻角，由 图２与表１可以看出，风速在ｌＯｍ／ｓ时，叶片各截面 还处于升力系数曲线的线性段，随后，各截面逐渐出 现失速，可见在叶片工作过程中，失速现象时有发 生。当叶型进行俯仰振荡运动时，其失速攻角比叶 型静止时的失速攻角要大，另外，叶型空气动力特性 随攻角变化的曲线出现迟滞现象，在攻角增大的过 程中，最大升力系数大于静态值；而在攻角减小的过 程中，则小于静态值，这种现象称之为动态失速。在 对风力机叶片做气动弹性分析时动态失速是一个很 重要的现象，需要重点考虑。 表１升力曲线直线段结束及最大值位置对应攻角 Ｔａｂｌｅ １ ａｔｔａｃｋ跗ｌｇｌｅ ｏｆ ｓｔｒａｉｇｈｔ ｌｉｎｅ ｓｅｇｍｅｎｔ
图３分离函数厂的物理意义 Ｆｉｇ．３ Ｔｈｅ ｍｅａ｜ｌｉｌｌｇ ｏｆ ｓｅｐａ＿ｍｔｉｏｎ ｆｕｎｃｔｉｏｎ ｆ 则有如下关系：
ｃ。＝（掣）ｃｍ
在非定常流中升力曲线的时变性来自于两个方 面，一方面是升力线斜率的时变，另一方面足分离点 位置的时变。
考虑在无分离的条件下升力线斜率的时变性， 当攻角发生小的扰动时，叶型上会暂时性的产生小 的分离，然后气流会重新附体，这一滞后过程导致了 环量的变化，从而导致升力的变化，Ｋｒｅｎｋ模型用脉 冲响应函数来描述这个过程，最终得到升力线的滞 后响应表示如下：
论计算微段的气动升力、气动阻力及气动力矩，并将 其折算到节点位置上，由于Ｍａｒｃ主程序要求 ＦＯＲＣＤＴ载荷子程序返回的是载荷增量，因此这里 还需要取出上一增量步的载荷值，从而获得增量值。
Ｍａｒｃ主程序根据返回的载荷增量，求解结构的响应， 如此实现反复迭带。
以某型７５０ｋＷ风力机叶片为例，根据以上的思 路，将叶片分成１５段，每一段由一个梁单元来表示， 这样每一个节点代表一个截面，由叶根至叶尖分别
叶片为薄壁结构，可得扭转振动微分方程：
７１ｐ ３。２。０＝券＋ｍＩ
截面上的剪应力与扭矩的关系为：
可得：
Ｉ ｒＩ ｒｌｄＡｌ＋Ｉ ｒ２ ｒ２ｄＡ２＋…＝Ｍ。
ｏ ４ｌ
ｏ＾２
鬟（Ｇ。，＾Ｉ ｒ２ｄＡ。＋Ｇ２Ｊ．＾２ ｒ；ｄＡ：力学响应分析
９１９
由此得到角位移和扭矩的关系：
ｆｄｌ（ｔ）＋ｃｃＪＩ Ｃｌ（ｔ）＝Ａｌ ｅｍ（ｔ）
｛ｅ２（ｆ）＋∞２ ｃ２（￡）＝Ａ２ｅ【ｏ（ｔ）
【ＣＩｏ．ｄ（ｔ）＝Ｃｍ（口）一ｃｌ（ｔ）一ｃ２（￡） 考虑分离点位置的时变性，在分离状态下，当攻 角改变时分离点位置也会改变，而位置的改变滞后 与攻角的改变，用一阶微分方程表示这种滞后关系， 得：