基于强迫振动和自由振动的风力机翼型涡激振动数值研究

李 羽, 柳 灿, 连 波, 竺晓程

(上海交通大学 机械与动力工程学院, 上海 200240)

随着风力机技术的蓬勃发展,风力机叶片大型化已成趋势,这将导致叶片固有频率降低,加剧了结构对风的敏感性。由于风力机处于自然条件下,需要考虑各种可能工况的影响,如停机状态下的叶片变桨失效,风力机叶片可能会处于非常大的攻角下,引起严重的失速和流动分离现象,可能会发生2种不同形式的大幅振动:失速颤振(Stall induced vibration,SIV)和涡激振动(Vortex-induced vibration,VIV),从而会导致叶片损伤或疲劳失效[[1-5]。现代风力机叶片的涡脱落频率接近叶片的低阶固有频率,因此存在涡激振动的可能性比失速颤振大[6-7]。

受虎门大桥等公共涡振事件的影响,涡激振动现象引起了广泛关注。之前的学者已经研究了圆柱[8-9]、桥梁[10-12]等结构的涡激振动现象,其中Sarpkaya[13]和Williamson等[14]对有关圆柱涡激振动给出了详细的综述。在风力机叶片方面,Skrzypiński等[7]对风力机翼型在大攻角下的涡激振动进行了数值研究,确定了翼型的涡振区间。Zou等[15-16]利用双尾迹涡模型预测了风力机翼型的涡激振动响应。Pellegrino等[17]和Meskell等[18]发现S809翼型在攻角为40°和140°时存在流固频率锁定问题。Heinz等[8]采用数值计算方法模拟了停机状态下风力机叶片的涡振现象,发现来流倾角对涡振的发生具有很重要的影响。Horcas等[19-20]分别采用尾缘襟翼和改变叶尖几何来抑制涡激振动。王渊博等[21]利用双向强流固耦合方法对NREL 5 MW风力机进行研究,发现在变桨故障下该风力机出现了大幅振动现象。总之,以上工作表明大型风力机叶片需要考虑涡振的危害性。

本课题组Hu等[22]采用强迫振动的数值方法,通过改变结构刚度,确定了Du96-W-180风力机翼型的涡振区间,并采用流动模态分解方法研究了尾流流动结构。Lian等[23]则采用自由振动的数值方法,确定了风力机翼型涡振的来流速度范围,分析了翼型涡激振动的发生、发展和演化过程。上述工作分别改变的是结构参数和来流速度。对于实际结构,确定涡振的来流速度区间具有实际参考意义,而在实验和数值研究中,通过改变振动频率的强迫振动方法更容易实现。为了进一步考察在相同频率比下来流速度变化对流场和振动响应的影响,笔者首先采用强迫振动方法,研究了来流速度变化对翼型涡振区间的影响,并与文献[21]的结果进行对比;然后结合自由振动方法,研究不同结构阻尼比对翼型振动响应的影响。

1 计算模型及验证

1.1 研究对象

以代尔夫特理工大学研发的Du96-W-180风力机(以下简称Du96)翼型为研究对象[23]。图1给出了翼型的几何形状和计算域示意图,其中来流速度U为29.89 m/s,弦长l为1 m。本研究中Du96翼型弦长设置为1 m,攻角为90°。进口边界距离翼型10倍弦长,给定速度边界条件。为了使尾迹充分发展,出口距离翼型40倍弦长,并设为压力出口。两侧边界距离翼型各10倍弦长,为自由滑移壁面,翼型表面为无滑移壁面。

图1 二维翼型Du96计算域及边界条件示意图

图2给出了翼型表面的网格示意图。对于翼型固定的情况,文献[22]对网格进行了无关性验证,确定总网格单元数为96 042,并保证翼型表面y+均小于2。采用商业软件ANSYS CFX进行非定常计算,选用SST湍流模型。

图2 翼型表面二维CFD网格示意图

1.2 强迫振动和自由振动计算方法

在翼型强迫振动计算中,利用动网格技术实现翼型沿弦向的周期性正弦运动,非定常时间步设置为振动周期的1/100。

在强迫振动的数值计算中,基于能量法确定气动阻尼比ξa:

(1)

(2)

式中:t0为起始时间;T为振动周期;p为翼型表面压力;v为翼型振动速度矢量;n为翼型表面法向量;dS为翼型表面微元面积;对于每个计算,N不小于10。

若ξa>0,表明翼型一直将振动的能量耗散到流体中,翼型处于稳定状态;反之,若ξa< 0,则说明翼型从流体中吸收能量,翼型振幅增大,即处于气弹不稳定状态。

在翼型自由振动的计算中,根据结构动力学方程,确定翼型的弦向自由振动,可通过动网格变形技术实现。

(3)

式中:y为刚性翼型的瞬时位移;c为线性结构阻尼系数;k为线性刚度系数;L为弹簧振幅。

计算中分别求解流动方程和结构方程。在每个实时步骤中,作用在翼型上的气动载荷都是从流场解中获得的。然后使用Newmark积分法求解结构动力学方程,更新翼型的位移,并相应地改变流场中翼型周围的网格,用于求解下一时间步流体的控制方程。

2 计算结果

2.1 强迫振动

频率比r定义如下:

(4)

式中:fc为翼型振动频率;l为弦长;St为静止翼型的斯特劳哈尔数;fv0为静止翼型的涡脱落频率。

从频率比的定义可以发现,通过改变振动频率可以实现不同频率比,也可以通过改变来流速度,以改变静止翼型的涡脱落频率,从而实现不同频率比。下文将分析2种方法对涡振区间及涡振特征的影响。

在来流雷诺数Re=2.1×106(对应的静止翼型脱落频率为3.738 Hz,St=0.12)时,图3给出了翼型在强迫振动频率为3.560 Hz,10%弦长振幅下的非定常升力系数时程变化曲线和频谱图,由于流动分离和翼型振动的影响,升力系数随时间存在明显的波动。从频谱分析可以发现,在振动稳定后升力系数具有明显的离散峰值频率,升力系数的峰值频率为3.560 Hz,静止翼型下的涡脱落频率不同,而与给定的强迫振动频率一致,显示了“锁频”特征。图3中还存在峰值频率的倍频,说明气动力具有非线性特征。

(a) 升力系数时程变化曲线

基于不同振幅的强迫振动计算,可通过分析升力系数的频谱来确定锁频区间。图4给出了不同来流速度下翼型在不同弦向振幅下的频率锁定区间,得到了具有V型结构的锁频区域,图中的横坐标是基于St得到的来流速度所对应的频率比。由图4可知,所得锁频区间与文献[22]中保持来流速度不变、改变翼型结构固有频率下得到的结果基本一致。

图5给出了2种方法在10%弦长振幅下气动阻尼比随频率比的变化曲线,2条竖虚线中间部分为锁频区域,实心曲线为文献[21]计算的不同结构频率的情况,可以观察到两者的变化趋势类似。其中锁频区域的气动阻尼由负到正,经历了迅速变化,在0.93≤r≤1.07时,翼型处于锁频状态,并且ξa都为负,在r=0.952处达到最小值,此时翼型处于气弹不稳定状态,若解除翼型强迫振动限制,振幅增大;在1.07

图5 10%弦长振幅下气动阻尼比随频率比的变化

表1给出了10%弦长振幅下r=0.93和r=1.11时的气动阻尼比ξa,可以发现,r=0.93时气动阻尼为负,而r=1.11时气动阻尼为正,气动阻尼比随着频率比的增大经历了由负到正的过程。

表1 10%弦长振幅下不同频率比的气动阻尼比

图6分别给出了改变来流速度和改变结构频率下r=0.93和r=1.11时1个振动周期内升力系数的变化曲线以及翼型振动位移曲线。从图6可以发现,对于r=0.93,位移和升力系数之间的相位差约为90°;而对于r=1.11,位移和升力系数基本处于同一相位。

(a) 改变来流速度

对比r=0.93和r=1.11时气动阻尼以及气动力变化曲线可以发现,改变来流速度和改变结构刚度对涡振具有相同的影响趋势,可以认为在本文的研究中锁频区域对来流速度不敏感,可以通过改变振动频率来研究翼型在不同频率比下的涡振响应。

针对频率比为0.93和1.11时强迫振动获得的气动阻尼差异,下面将通过自由振动计算方法来检验2种频率比下的振动响应,并与强迫振动计算的结果进行对比。

2.2 自由振动

在自由振动计算中,将考虑不同结构阻尼对翼型涡振的影响,图7给出了结构阻尼比ξs为2%和6%的自由振动涡振区间,竖虚线代表强迫振动振幅为10%弦长时的涡振边界,横点划线代表强迫振动规定的10%振幅(0.1 m)。可以发现2%结构阻尼比时翼型的振幅明显大于6%结构阻尼比,自由振动计算获得的涡振区间与图5中强迫振动计算的负气动阻尼锁频区间基本吻合。

图7 不同结构阻尼下自由振动翼型涡振区间

图8给出了2%和6%结构阻尼比时的自由振动翼型升力系数频率变化曲线,可以发现2种结构阻尼比下升力系数频率均会在锁频区域锁定到翼型固有频率,而且2%结构阻尼比的锁频范围大于6%结构阻尼比,差异主要集中在大来流风速段。

图8 不同结构阻尼下自由振动翼型升力系数频率图

表2给出了2%、3%和6%结构阻尼比下改变来流速度时r=0.93和r=1.11的升力系数频率,其中翼型固有频率fn=3.560 Hz,不同来流速度下静止翼型涡脱落频率fv0分别为3.827 Hz和3.207 Hz。

表2 r=0.93和r=1.11时不同结构阻尼的翼型升力系数频率

对比不同来流速度下的计算结果可知,当ξs=2%时,升力系数频率均接近翼型固有频率,基本满足涡激振动时锁定到固有频率这一现象;但在ξs=6%的情况下,r=0.93时仍然为涡振锁频振动,而r=1.11时频率接近涡脱落频率,表现为强迫振动。从表2可以发现,不同结构阻尼下锁频的频率与翼型固有频率在数值上存在一定偏差,Sarpkaya[13]认为这是由流固耦合的虚拟质量力引起的。为了进一步确定原因,下面给出不同计算条件下的振动响应时程变化曲线和频谱图。

以静止翼型充分发展的非定常计算结果作为初始流场,开展不同频率比和结构阻尼比的自由振动计算,图9给出了ξs=2%、r=0.93时翼型的振动响应时程曲线及频谱图。从图9可以看出,由于流动分离和脱落的影响,振幅具有明显的波动,而且表现出拍振现象。对比频谱可以发现,频谱中存在最大峰值3.543 Hz和第二峰值3.815 Hz的振动频率,其中最大峰值接近翼型的固有频率,结合振幅幅值可以认为此工况条件下具有锁频特征。

(a) 振动位移时程变化曲线

图10给出了ξs=6%、r=0.93时的振动响应时程曲线及频谱图。当ξs=6%、r=0.93时,由于强迫振动计算中10%弦长振幅下r=0.93时的气动阻尼比(-8.398%)和结构阻尼比之和小于0,振动开始发散,振幅逐渐增大并超过10%弦长,翼型达到最大振幅0.155 m后幅值趋于稳定,峰值频谱为3.602 Hz,也表现出锁频特征。

(a) 振动位移时程变化曲线

图11给出了ξs=2%、r=1.11时的翼型振动响应曲线和频谱图,此时翼型振幅逐渐增大到0.269 m,振动峰值频率为3.497 Hz,也具有锁频特征,但与固有频率略有差异。

(a) 振动位移时程变化曲线

而当ξs=6%、r=1.11时,振幅达到最大值0.093 m之后逐渐减小,最终稳定在0.072 m,且振动频率为3.202 Hz,与翼型结构固有频率有明显差异,而且接近静止翼型的涡脱落频率,结合振幅和频率,可以认为此时振动响应属于强迫振动,如图12所示。

(a) 振动位移时程变化曲线

3 结论

利用动网格技术实现了不同来流速度和不同结构频率下翼型的强迫振动和自由振动计算,对90°攻角下的风力机翼型振动现象进行了数值研究,在本计算条件下,通过改变来流速度与改变结构振动频率,得到相近的V型锁频区域,获得了以下主要结论:

(1) 对于10%弦长振幅的强迫振动数值计算,频率比在0.93～1.11内时,涡脱落频率会锁定到固有频率上,而气动阻尼则经历了从负到正的改变,而且位移和升力的相位差发生变化。

(2) 自由振动数值计算验证了强迫振动方法得到的涡振区间,对于r=0.93,由于10%弦长下强迫振动计算得到大的负气动阻尼,在结构阻尼比为6%时仍然保持涡振状态,而且不同结构阻尼下锁频中的升力系数频率与翼型固有频率存在一定的差异。而对于r=1.11,由于10%弦长下强迫振动计算得到正气动阻尼,在2%结构阻尼比下表现为涡激振动,而在6%结构阻尼比下表现为强迫振动,采用强迫振动方法确定涡振区间边界时需要注意不同结构阻尼的影响。