0 引言
为从风中吸取更多的能量,降低风力发电成本,风力机尺寸日益增大,风力机叶片也变得更细长,叶片结构柔度增加,气动载荷、重力、惯性载荷及操纵载荷同时随着叶片长度的增加急剧上升[1],其带来的气弹问题也越发突出。气动弹性产生的振动将会显著影响风力机的运行寿命。因此,需要进行气动结构耦合分析来预测非稳态气动载荷与叶片弹性变形的复杂相互作用。
叶素动量（blade-element momentum, BEM）理论、涡尾迹方法和计算流体力学（computational fluid dynamics, CFD）是目前较常用的风力机气动分析方法[2]。BEM由于其简单易用、计算量小等优点被广泛使用。但由于传统的二维BEM模型忽略了气流沿叶片的展向流动,不能反映叶素单元间的相互影响,导致其不能很好地预测风轮的状态,需要人为引入修正模型,这在一定程度上降低了该理论的普适性和可靠性,尤其是对于气动流场存在三维展向流动以及非定常特性的情况[3]。CFD是一种高计算精度的叶片气动分析方法,通过求解N-S控制方程获得叶片详细的流场信息和气动载荷,可以准确地计算湍流、非定常等复杂流场结构对风力机气动性能的影响,但CFD不仅需要较高的计算成本,而且数值计算时需要与结构界面做特殊处理,这也使其在气弹分析方面应用较少[4]。与BEM和CFD这两种方法对比,涡尾迹方法从构造叶片表面的涡量出发,用不同涡量分布形式来模拟叶片与三维流场之间的相互作用。陈进格等[5]采用升力面配合自由尾迹预测法,计算尾流场对叶素的诱导作用,适用于大型柔性叶片三维气动分析。李德源等[6]采用升力线理论建立叶片的气动模型,计算5 MW叶片的气动性能并与引入修正模型的BEM对比,结果表明升力线气动模型具有较高的计算稳定性。GRECO等[7]采用一种自由涡尾迹和面元法全耦合的气动特性计算模型计算叶片在不发生严重的流量分离时的气动性能,不需BEM那样引入叶尖、叶根修正模型,不需CFD那样高的计算成本,便能得到精确的结果。
在气弹耦合分析中,赵鹰等[8]采用修正的BEM理论,结合基于有限元复合材料的叶片结构,建立了叶片的气弹耦合模型,结果表明叶片气动弹性效应会造成叶片功率输出下降,导致达不到设计功率。SAYED等[9]采用时间精确的CFD-计算结构力学（computational structural dynamics, CSD）弱耦合方法,研究了DTU10MW风力机的气弹响应,发现考虑叶片变形时风轮功率与推力提高了约1%和0.3%。钱晓航等[10]运用几何精确梁方法构建NREL5MW和IEA15MW风力机叶片,结合叶素动量理论建立了非线性叶片气弹耦合模型,采用线性分析与非线性分析方法计算了湍流风速时两个风力机组的动态响应,结果表明对于大型柔性叶片的风力机,应考虑几何非线性对风力机叶片气动弹性响应的影响,以准确评估风力机设计的安全性和稳定性。然而,目前综合考虑大型风力机整机气弹耦合特性的研究相对较少。
1 风力机多体动力学建模
1.1 整机机组离散模型
大型风力机的叶片为细长的柔性体,翼型截面形状与大小沿展向变化,非线性变形明显。为此,采用6个“超级单元”模型将叶片离散为由刚体、弹簧、阻尼器和运动副组成叶片的多刚体系统[13,14],其中超级单元模型方法是通过将需要分析的构件分成若干单元,每个单元由4个刚体组成,相邻刚体之间通过运动副连接,并用弹簧和阻尼器约束刚体间的相对运动,采用相关计算方法获取能准确模拟构建运动的弹簧刚度系数。同样用3个超级单元离散塔架,用一个刚体简化机舱,塔顶与机舱用万向节连接,用来模拟机舱的偏航与俯仰运动,将主轴划分为两个刚体,中间用转动铰连接,其中一个刚体与机舱用转动铰连接,用来模拟主轴的转动,另一个刚体与轮毂固定连接,用一个刚体简化轮毂,三个叶片的第一个刚体分别与轮毂以转动铰连接,用来模拟叶片的变桨,从而将整机离散为71个刚体、116个自由度,由弹簧阻尼约束的刚柔混合多体系统,如图1所示。经过单元相关性数值分析,风力机系统可以用相对较少的自由度,准确地描述整机和叶片的结构动态特性,满足时域气弹响应分析要求。
Fig. 1 The hybrid multibody system model of DTU 10MW wind turbine图1 DTU10MW整机混合多体模型 |
1.2 风力机多体动力学
取各转动铰的旋转自由度为广义坐标$q(t)$,则整机多体系统的拓扑构型的广义坐标阵为：
$q(t)={{\left[ {{q}_{1}}\text{(}t\text{),}\cdot \cdot \cdot \text{,}{{q}_{s}}(t) \right]}^{\text{T}}}$ (1)
式中：${{q}_{\text{s}}}$的下标s代表风力机多体系统的自由度;上标T代表转置。根据虚功形式的动力学普遍方程,考虑整机受到各种几何或运动约束,如塔底的固定约束、主轴的转速约束,系统的动力学方程具体形式如下：
$\mathbf{Z}\ddot{q}\mathbf{+}\Phi _{\mathbf{q}}^{\text{T}}\lambda \mathbf{=}z$ (2)
其中：$\mathbf{Z}$为广义质量阵;${{\Phi }_{q}}$为约束方程的雅可比矩阵;$\lambda $为拉格朗日乘子阵;$\Phi _{\mathbf{q}}^{\text{T}}\lambda $为刚体间所有理想约束力的贡献;$z$为广义力阵。
采用双循环隐式积分方法[15]求解风力机整机的动力学微分代数方程。该方法采用LU分解法将系统的自由度分为独立自由度${{q}_{\text{I}}}$与非独立自由度${{q}_{\text{D}}}$,则可将动力学方程写为
$\left[ \begin{matrix} {{Z}_{\text{DD}}} & {{Z}_{\text{DI}}} \\ {{Z}_{\text{ID}}} & {{Z}_{\text{II}}} \\ \end{matrix} \right]\left\{ \begin{matrix} {{{\ddot{q}}}_{\text{D}}} \\ {{{\ddot{q}}}_{\text{I}}} \\ \end{matrix} \right\}+{{\left[ \begin{matrix} {{\Phi }_{{{q}_{\text{D}}}}} & {{\Phi }_{{{q}_{\text{I}}}}} \\ \end{matrix} \right]}^{\text{T}}}\lambda =\left\{ \begin{matrix} {{z}_{\text{D}}} \\ {{z}_{\text{I}}} \\ \end{matrix} \right\}$ (3)
约束方程对时间求二阶导数,得到
${{\Phi }_{{{q}_{\text{I}}}}}{{\ddot{q}}_{\text{I}}}\text{+}{{\Phi }_{{{q}_{\text{D}}}}}{{\ddot{q}}_{\text{I}}}-\gamma \text{=}0$ (4)
联立式（3）、式（4）可将动力学方程缩并至仅包含独立自由度：
$\tilde{Z}{{\ddot{q}}_{\text{I}}}=\tilde{z}$ (5)
式中：$\tilde{Z}={{Z}_{\text{II}}}-{{Z}_{\text{ID}}}\Phi _{{{q}_{\text{D}}}}^{-\text{1}}{{\Phi }_{{{q}_{\text{I}}}}}\text{+}\Phi _{{{q}_{\text{I}}}}^{\text{T}}\Phi _{{{q}_{\text{D}}}}^{-\text{T}}{{Z}_{\text{DD}}}\Phi _{{{q}_{\text{D}}}}^{-\text{1}}{{\Phi }_{{{q}_{\text{I}}}}}-\Phi _{{{q}_{\text{I}}}}^{\text{T}}\Phi _{{{q}_{\text{D}}}}^{-\text{T}}{{Z}_{\text{DI}}}$;$\tilde{z}={{z}_{\text{I}}}-{{Z}_{\text{ID}}}\Phi _{{{q}_{\text{D}}}}^{-\text{T}}\gamma -\Phi _{{{q}_{\text{I}}}}^{\text{T}}\Phi _{{{q}_{\text{D}}}}^{-\text{T}}{{z}_{\text{D}}}\text{+}\Phi _{{{q}_{\text{I}}}}^{\text{T}}\Phi _{{{q}_{\text{D}}}}^{-\text{T}}{{Z}_{\text{DD}}}\Phi _{{{q}_{\text{D}}}}^{-\text{T}}\gamma $。
求解上述方程的基本思想是先利用式（3）、式（4）中已按独立非独立自由划分的矩阵组集出缩并后的质量阵$\widetilde{Z}$、力阵$\widetilde{z}$,并由式（5）求解出独立自由度的广义加速度${{\ddot{q}}_{\text{I}}}$,然后与非独立自由度的广义加速度${{\ddot{q}}_{\text{D}}}$重新组集$\ddot{q}$,最后由$\ddot{q}$利用数值积分方法求解速度$\dot{q}$和位移$q$。
2 基于升力线叶片的气动载荷计算
Fig. 2 Schematic diagram of the vortex method model图2 涡方法模型示意图 |
Fig. 3 Schematic diagram of spiral tail vortex图3 螺旋尾涡示意图 |
线涡对空间中的目标点位置处产生的诱导速度可通过Biot-Savart公式描述：
$\text{d}{{W}_{i}}=\frac{\Gamma }{4\text{ }\!\!\pi\!\!\text{ }}\frac{l\times \text{d}s}{{{\left| l \right|}^{3}}}$ (6)
$s=r\cos (\theta +{{\theta }_{k}})i+r\sin (\theta +{{\theta }_{k}})j+r\theta \tan (\phi )k$ (7)
式中：r为线涡脱落的叶片展向位置;$\theta $为螺旋尾涡的方位角;${{\theta }_{k}}$为目标点所在的叶片与尾涡起始的叶片之间的夹角;$\phi $为螺旋尾涡的俯仰角,与来流风速U、叶片转速$\Omega $有关,$\phi \text{=}{{\tan }^{-1}}\left( \frac{U}{r\Omega +{{w}_{\text{y}}}} \right)$。
诱导速度${{W}_{n}}$垂直于未受扰动相对风速$W'$,并给出附着涡任意位置${r}'$[已知惯性坐标(${x}',{y}',{z}'$)]处由所有尾涡引起的诱导速度${{W}_{\text{n}}}$为[17]：
${{W}_{n}}\left( {{\xi }'} \right)=\underset{{{\xi }_{\text{hub}}}}{\overset{1}{\mathop \int }}\,\frac{\frac{\text{d}\Gamma }{\text{d}\xi }}{4\text{ }\!\!\pi\!\!\text{ }{{R}_{\text{L}}}}\underset{K=1}{\overset{N}{\mathop \sum }}\,\underset{0}{\overset{\infty }{\mathop \int }}\,\frac{{{N}_{1}}{\xi }'+{{N}_{2}}{{\lambda }_{0}}}{{{D}_{1}}^{3/2}{{D}_{2}}^{1/2}}\text{d}\theta \text{d}\xi $ (8)
${{R}_{\text{L}}}=R\cos \psi {{N}_{1}}={{\xi }^{2}}-\xi {\xi }'\cos \left( \theta +{{\theta }_{k}} \right)$
${{N}_{2}}=\frac{{{V}_{0}}}{\Omega R}\left\{ \xi \left[ \theta \sin \left( \theta +{{\theta }_{k}} \right)+\cos \left( \theta +{{\theta }_{k}} \right) \right]-{\xi }' \right\}$
${{D}_{1}}={{\xi }^{2}}+{{{\xi }'}^{2}}-2\xi {\xi }'\cos \left( \theta +{{\theta }_{k}} \right)+{{\left( \frac{{{V}_{0}}\theta }{\Omega R} \right)}^{2}}$
${{D}_{2}}={{\lambda }_{0}}^{2}+{{{\xi }'}^{2}}$,$~{{\theta }_{k}}=\frac{2\pi \left( N-K \right)}{N}$
$W'=\sqrt{{{({{V}_{0}}+{{V}_{\text{h}}})}^{2}}+{{\left( \Omega r+{{V}_{\text{b}}} \right)}^{2}}}$
${{\alpha }_{i}}=\text{ta}{{\text{n}}^{-\text{1}}}\left( \frac{{{W}_{n}}}{W'} \right)$,$W={{W}_{n}}/\cos ({{\alpha }_{i}})$
式中：R为风轮半径,$\xi =x/R$,${\xi }'={x}'/R$为叶片展向无量纲位置参数;${{\lambda }_{0}}={{V}_{0}}/(R\Omega )$;$\psi $为风轮锥角;$N$为叶片数;$K$为目标点所在叶片序号;${{\alpha }_{i}}$为诱导攻角;$W$为相对风速;${{V}_{\text{h}}}$、${{V}_{\text{b}}}$分别为叶片挥舞、摆振方向的振动速度。尾涡远端对目标点的诱导速度逐渐减弱,计算适当圈数尾涡的诱导速度贡献即可,本文计算14圈尾涡对目标点诱导速度的贡献。
由Kutta-Joukowski定理可知,环量与升力系数的关系为
$\Gamma =0.5cW{{C}_{\text{L}}}$ (9)
式中：c为翼型弦长;${{C}_{L}}$为升力系数。
由于诱导速度表达式（6）和式（8）中环量项$\Gamma \text{(}\xi \text{)}$未知。设附着涡环量为Fourier级数,并且满足边界条件$\Gamma \text{(}\xi \text{)}=\Gamma \text{(}1\text{)=0}$,由于叶片靠近轮毂中心的截面为圆筒状,其截面升力系数在任意攻角下恒为0,即环量恒为0,则设环量函数为：
$\Gamma (\xi )=\underset{m=1}{\overset{\infty }{\mathop \sum }}\,{{A}_{m}}\sin (m\text{ }\!\!\pi\!\!\text{ }\xi )$ (10)
式中：无穷级数求和可用有限项求和来近似计算,项数m越多,环量$\Gamma (\xi )$的计算精度越高,但所需计算时间也越长,一般m不宜取过大。综合考虑计算速度与计算精度,本文m取13。
在小攻角范围内,升力系数与有效攻角之间可近似为线性关系：
${{C}_{\text{L}}}={{a}_{0}}{{\alpha }_{\text{e}}}+{{b}_{0}}$ (11)
式中：${{a}_{0}}$、${{b}_{0}}$分别为斜率与截距。
综合式（6）~ 式（10）
$\left[\begin{array}{l}\sum_{m=1}^{M} m A_{m} \int_{\xi_{\mathrm{hub}}}^{1} \frac{\cos (m \pi \xi)}{4 R_{\mathrm{L}}^{2} \Omega \sqrt{\lambda_{0}^{2}+\xi^{\prime 2}}} \frac{I}{\xi-\xi^{\prime}} \mathrm{d} \xi+ \\ \tan ^{-1}\left(\frac{\lambda_{0}}{\xi^{\prime}}\right)-\beta\end{array}\right] a_{0}+b_{0} \\ =\frac{2}{c R_{\mathrm{L}} \Omega \sqrt{\lambda_{0}^{2}+\xi^{\prime 2}}} \sum_{m=1}^{M} A_{m} \sin (m \pi \xi)$ (12)
通过m元一次方程（12）的求解,可获得${{A}_{m}}$的初始解,进而可计算出初始环量$\Gamma (\xi )$与有效攻角${{\alpha }_{\text{e}}}$,结合初值联立式（9）和式（10）建立新的m元一次方程并求解得出新的${{A}_{m}}$,利用迭代法,可以得到更高数值精度的有效攻角${{\alpha }_{\text{e}}}$和诱导攻角${{\alpha }_{i}}$。通过翼型升阻力系数数据表插值得到${{C}_{L}}$和${{C}_{\text{D}}}$后,通过沿着叶片展向积分得到轴向推力F、转动力矩Q和功率P。
$F=\underset{{{\xi }_{\text{hub}}}}{\overset{1}{\mathop \int }}\,\left( \frac{1}{2}\rho Nc{{W}^{2}} \right)\left( {{C}_{\text{L}}}\cos {\phi }'+{{C}_{\text{D}}}\sin {\phi }' \right)R\cos \psi \text{d}\xi $(13)
$Q=\underset{{{\xi }_{\text{hub}}}}{\overset{1}{\mathop \int }}\,\left( \frac{1}{2}\rho Nc{{W}^{2}} \right)\left( {{C}_{\text{L}}}\sin {\phi }'-{{C}_{\text{D}}}\cos {\phi }' \right){{R}^{2}}\cos \psi \xi \text{d}\xi $(14)
$P=\underset{{{\xi }_{\text{hub}}}}{\overset{1}{\mathop \int }}\,\left( \frac{1}{2}\rho Nc{{W}^{2}} \right)\left( {{C}_{\text{L}}}\sin {\phi }'-{{C}_{\text{D}}}\cos {\phi }' \right){{R}^{2}}\Omega \text{cos}\psi \xi \text{d}\xi $(15)
${\phi }'=\text{ta}{{\text{n}}^{-\text{1}}}\frac{{{V}_{\text{0}}}+{{V}_{\text{h}}}}{\Omega {x}'+{{V}_{\text{b}}}}-{{\alpha }_{i}}~ \text{ }\!\!~\!\!\text{ }\!\!~\!\!\text{ }~~~$ (16)
3 气弹耦合程序设计
整机气弹耦合方程是一个典型的非线性代数微分方程,求解方程组可以采用数值积分方法[14]。整机的气弹耦合程序的非线性代数微分方程的建立以及数值积分流程如图4所示,首先程序将会读取系统的基本数据和仿真数据,如系统的拓扑构型参数、刚体的惯性参数、铰的参数等,接着计算不随时间变化的中间常值数据,如通路阵等。完成数据的输入与前处理后,程序进入数值积分过程,将上一时刻气动模块计算的气动载荷以外力的形式加入整机的多体动力学模块中,通过求解动力学方程得到这个时刻的广义坐标、广义速度,进而计算出各刚体气动中心的速度、扭角等,并与风速一起作为这个时刻升力线模块的输入参数,通过升力线模块计算有效攻角、入流角等参数,结合升阻力系数表,计算出各刚体上的气动载荷,并作为下一时刻整机的多体动力学模块的外力,如此循环直至仿真结束时间tend。本文采用MATALB中的ode15s求解器来求解系统的常微分方程组,积分精度设为1 × 10-5,输出步长为0.01 s。
Fig. 4 Flow chart of aeroelastic coupling analysis图4 气弹耦合分析流程图 |
4 计算结果与分析
4.1 稳态风工况
如图5所示,当风速低于12 m/s时,升力线模型计算的风轮功率略高于HAWC的结果,而推力基本一致,仅在风速11 m/s下相差较大,当风速大于12 m/s时,升力线模型计算的风轮功率小于HAWC的结果,推力基本一致,仅在高风速下略小于HAWC的推力。
Fig. 5 Power and thrust of each wind speed图5 各风速功率、推力 |
综上可知,叶素动量理论通过引入诱导速度修正和叶尖损失修正来考虑尾流场对叶素的诱导影响,而升力线理论未考虑任何修正,直接计算尾流场对叶素的诱导速度,两组模型有效攻角和气动载荷沿叶展方向的趋势基本一致,表明升力线理论结合超级单元结构的气弹耦合模型是有效的、可行的。
Fig. 6 Distribution of effective angle of attack and relative wind speed along spanwise direction图6 有效攻角、相对风速沿展向分布 |
Fig. 7 Distribution of lift and drag coefficients along spanwise direction图7 升力、阻力系数沿展向分布 |
Fig. 8 Distribution of tangential force and axial force along spanwise direction图8 切向力、轴向力沿展向分布 |
4.2 非稳态风况气弹响应
Fig. 9 The curves of wind speed with time in the hub图9 轮毂处风速随时间变化曲线 |
考虑DTU10MW风力机风轮直径达到178 m,轮毂高度达119 m,叶片在转动过程中,叶素所处高度变化大,因此引入风速随高度变化的指数规律分布加以修正：
$~\frac{V}{{{V}_{0}}}={{\left( \frac{H}{{{H}_{0}}} \right)}^{\alpha }}$ (17)
式中：$V$为距离地面高度H处的平均风速;${{V}_{0}}$为距离地面高度${{H}_{0}}$处的平均风速;α为切变系数,切变系数主要与地表粗糙度和大气稳定度有关,本文取0.12。
图10为考虑随机风时风力机功率、推力的时域响应曲线。在考虑随机风时,功率与推力的变化趋势与随机风的变化趋势一致,功率最大可达11.53 MW,超过额定功率10 MW,最小4.17 MW,变化幅度为7.36 MW,功率推力最大可达1.49 MN,最小0.97 MN,变化幅度为0.52 MN,在设计大型风力机时,需要考虑功率的大幅度变化对电网电能质量的影响以及推力的大幅度变化对风力机安全运行的影响。
Fig. 10 Curve of power and thrust with time图10 功率、推力随时间变化曲线 |
图11为考虑随机风时叶片叶尖相对于轮毂中心的挥舞、摆振位移时域响应曲线。由图可知,考虑随机风时,叶尖挥舞位移最大可达7.50 m,最小4.75 m,变化幅度为2.75 m,叶尖摆振位移最大为0.441 m,最小为0.161 m,变化幅值为0.28 m,总的变化趋势与随机风一致。在设计大型风力机叶片时,需要考虑叶片叶尖的大位移可能会导致叶片与塔架的碰撞,如设计预弯叶片增大叶尖与塔架的间距,或设计自适应降载的后掠叶片减小叶尖位移。相较于叶尖挥舞位移,叶尖摆振位移随风速的快速变化振动得更剧烈。
Fig. 11 Curve of tip flapwise and edgewise displacement with time图11 叶尖挥舞、摆振位移随时间变化曲线 |
图12为考虑随机风时塔顶的横向、纵向位移的时域响应曲线。可以看出,塔顶横向、纵向位移的变化趋势与随机风的变化趋势一致,随机风对塔顶位移的影响很大,在大型风力机的设计过程中,需要校核随机风况下塔架的强度,并需要开发更能反映实际风场特性的湍流风模型。
Fig. 12 Curve of transverse and longitudinal displacement of tower top with time图12 塔顶横向、纵向位移随时间变化曲线 |
5 结论
基于升力线理论与超级单元模型耦合建立风力机整机的气弹耦合模型,并与HAWC模型的稳态风况对比验证,然后用建立的气弹耦合模型进行随机风的数值仿真,结果表明：
（1）基于升力线理论,可以在不考虑修正的情况下计算出与修正后的叶素动量理论相近的风力机气动性能和响应。
（2）所建立的气弹模型能够很好地完成气动和结构的耦合仿真,可以很好地模拟风力机在随机风况下的气动载荷与结构位移的动态响应,该结果对风力机叶片设计具有一定的参考价值。
（3）相对于需要引入修正模型的叶素动量理论而言,升力线理论考虑了风力机三维尾流场对叶片的诱导速度,更适合分析积叠线弯曲、三维流动显著的叶片,如后掠、预弯叶片。
