双馈风电机组传动链低电压穿越响应分析

尹尧杰, 褚景春, 姜培学

(1.国电联合动力技术有限公司, 北京 100039; 2.国家能源集团新能源技术研究院有限公司, 北京 102209; 3.清华大学热能工程系, 北京 100084)

在国家碳达峰和碳中和战略的背景下,风电产业迎来了历史性的发展机遇期,中国风电经过30多年的发展,逐渐成为国家电力供应的主力电源,截至2022年底风电机组累计装机容量约3.5×108kW[1]。大规模的风电并入电网,使得风电机组与电网的机电耦合性增强,为了保证整个电网的安全稳定运行,电网对风电机组并网提出了很多的要求,其中一个重要方面就是要求风电机组具备低电压穿越的能力即当电网故障或者扰动引起并网点电压跌落时,在一定电压跌落范围内,风电机组不允许脱网运行[2]。双馈风电机组发电机转子经过变流器和变压器与电网相连,而发电机定子则直接用过变压器连接到电网。电网电压的变化会直接反映到发电机定子电压,引起发电机电磁转矩的变化,因此双馈型风电机组对电网故障比较敏感[3]。风电机组并网点发生电压跌落时,发电机电磁转矩会迅速减小,引起风轮气动功率与发电机电功率的不平衡,产生的振荡载荷通过风电机组传动链进行传递。

中外学者对风电机组低电压穿越的影响进行了很多的研究,如文献[4]提出了一种针对双馈风电机组不对称短路电流的实用计算方法,计及不对称故障后的低电压穿越策略。文献[5]为了改善双馈风电机组故障穿越后机组的稳定性,提出了通过励磁的方法实现低电压穿越的计算机控制策略。文献[6]提出了模型预测转子电流控制与动态电压恢复器的协调控制策略,提高双馈风电机组不脱网运行的能力。文献[7]针对兆瓦级风电机组建立了Bladed仿真模型,研究了低电压穿越过程中机组运行特性以及关键部位载荷特性,结果表明和低电压故障对机组塔架载荷有重要影响。文献[8]建立Bladed-MATLAB联合仿真模型,研究了电网对称与不对称故障对传动链轴系及塔架载荷特性的影响。

从以上分析可以看出大多数学者的研究都集中在风电机组低电压穿越过程控制策略与电气部件影响,仿真采用的机械模型也相对简单,机械模型中并未对齿轮与轴承等部件载荷特性进行深入的分析。

为了准确评估与分析风电机组低电压穿越过程对传动链齿轮与轴承等部件的影响,现分别采用转子有限元法、集中参数建模法以及齿轮啮合理论对双馈风电机组传动链主要部件进行动力学建模。将各个部件的动力学模型进行组装建立双馈风电机组传动链整体的动力学模型,结合风电机组气动载荷计算模型与控制模型可以对传动链的动态响应进行仿真计算。基于某型双馈风电机组低压电穿越过程的测试曲线,对比分析不同形式与类型的电压跌落对风电机组传动链齿轮与轴承载荷的影响。

1 双馈风电机组传动链部件动力学模型

1.1 转子有限元法[9]

风电机组传动链中的盘、弹性轴、轴承以及联轴器等转子系统基本部件均可以采用转子有限元法进行动力学建模。转子有限元法建模的步骤为:将转子系统离散化为有限多个单元节点,以单元节点为广义坐标,通过拉格朗日方程建立单元的运动微分方程。

1.1.1 刚性盘模型

刚性盘模型可以对风电机组传动链中啮合的齿轮副以及发电机转子盘进行建模。刚性盘动能表达式为

(1)

将刚性盘的动能表达式代入拉格朗日方程可以求得刚性盘的运动微分方程为

(2)

1.1.2 弹性轴模型

转子有限元法中通常采用考虑剪切变形的铁木辛柯梁对弹性轴进行建模。风电机组传动链中的弹性轴包括叶片、主轴、齿轮箱内部轴以及发电机轴。梁单元包含2个节点,每个节点包含6个自由度。

利用梁单元节点与相应的插值函数可以推导得出梁单元微段的动能与势能表达式,沿梁单元长度方向进行积分即可得到梁单元的动能与势能表达式。同样的将梁单元的动能与势能表达式代入拉格朗日方程就可以求得梁单元的运动微分方程。梁单元的运动微分方程为

(3)

1.1.3 支承轴承

在转子有限元法中通常采用弹簧-阻尼模型对轴承进行建模,忽略轴承的质量,轴承仅仅提供刚度与阻尼。风电机组传动链中的轴承包含主轴承、齿轮箱内部轴承以及发电机端的轴承。轴承的运动方程为

(4)

式(4)中:Cb为轴承阻尼矩阵;Kb为轴承刚度矩阵;Qb为轴承外力矩阵。

1.1.4 联轴器[10]

风电机组传动链中齿轮箱输出轴通过联轴器与发电机轴相连。转子有限元法中,采用各向同性的弹性部件对联轴器进行建模。联轴器刚度矩阵的表达式为

(5)

式(5)中:Kcoupling为联轴器刚度矩阵;K11为联轴器刚度分块矩阵。

1.2 平行齿轮动力学模型

平行齿轮动力学模型可以采用集中参数建模的方法并结合齿轮啮合理论建立。齿轮啮合模型采用质量-弹簧模型,弹簧沿齿轮啮合方向,弹簧刚度即为齿轮之间的啮合刚度。齿轮啮合模型如图1所示。

oxyz为齿轮啮合总体坐标系;x′y′z′为齿轮啮合局部坐标系;局部坐标系的原点位于齿轮啮合线与齿轮中心连线的交点P,即齿轮发生啮合的点;φ角为齿轮的方位角;d、h分别为啮合齿轮中心在水平方向与垂直方向的距离

齿轮的啮合力在齿轮啮合局部坐标系分解为

(6)

式(6)中:Fn为齿轮啮合力;F′x、F′y、F′z为齿轮啮合力在局部坐标系的分量;φx、φy、φz为啮合力与局部坐标系坐标轴的夹角。

局部坐标系的力可转换成主坐标系中的力,转换关系为

Fxi=F′xsinφ+F′ycosφ

(7)

Fyi=-F′xcosφ+F′ysinφ

(8)

Fzi=F′z

(9)

为了建立齿轮的啮合矩阵需要将主坐标系中每个方向的位移分量向啮合方向分解,得到啮合方向的总位移即可求得齿轮之间的啮合力。

齿轮在啮合方向的平动位移为

uxyz=(sxj-sxi)+(syj-syi)+(szj-szi)

(10)

式(10)中:sxi、sxj为x轴平动位移分量;syi、syj为y轴平动位移分量;szi、szj为z轴平动位移分量。

齿轮在啮合方向的角位移为

uθ=(-sθxj-sθxi)+(-sθyj-sθyi)+(-sθzj-sθzi)

(11)

式(11)中:sθxi、sθxj为x轴方向角位移分量;sθyi、sθyj为y轴方向角位移分量;sθzi、sθzj为z轴方向角位移分量。

齿轮啮合方向的总位移为

umesh=uxyz+uθ

(12)

齿轮的啮合力为

Fn=kgumesh

(13)

式(13)中:kg为齿轮啮合刚度。

将啮合力在主坐标系中进行分解,即可求得啮合力在主坐标系各个方向的力与力矩。然后将力与力矩按照自由度对应顺序放入对应位置即可得到齿轮的啮合矩阵。

1.3 行星齿轮动力学模型

行星齿轮主要部件包括行星架、内齿圈、太阳轮以及行星轮。建立如图2所示的行星齿轮动力学模型。

x、y、z、θx、θy、θz分别为行星齿轮各个部件在各自坐标系的平动位移与转角位移,不同的部件采用不同的下标;r为内齿圈;c为行星架;s为太阳轮;n为行星轮(n=1,2,…,N);N为行星轮个数

行星齿轮动力学模型的建立需要首先分析其外啮合与内啮合位移关系,然后通过对行星齿轮各部件进行运动受力分析建立其运动微分方程,将行星齿轮各部件模型进行组装可以得到行星齿轮的动力学模型。

1.3.1 外啮合位移关系

太阳轮与行星轮的外啮合如图3所示。

图3 行星齿轮外啮合Fig.3 External mesh of planetary gear

由图3可推导得到太阳轮与行星轮之间的啮合位移关系为

δsn=xssinψsncosβb-yscosψsncosβb-zssinβb-
rmsθsxsinψnsinβb+rmsθsycosψnsinβb-
rsθszcosβb-xnsinα′tscosβb+
yncosα′tscosβb+znsinβb+
rnθnysinβb/cosα′ts-rnθnzcosβb

(14)

式(14)中:xs、ys、zs为太阳轮的平动位移;θsx、θsy、θsz为太阳轮的转角位移;rms为太阳轮的节圆半径;rs为太阳轮的基圆半径;xn、yn、zn为行星轮的平动位移;θnx、θny、θnz为行星轮的转角位移;rn为行星轮的半径;ψn为第n个行星轮的方位角;α′ts为太阳轮与行星轮啮合的端面压力角;α′ts为太阳轮与行星轮啮合的端面压力角;βb为螺旋角。

1.3.2 内啮合位移关系

行星轮与内齿圈的内啮合如图4所示。

图4 行星齿轮内啮合Fig.4 Internal mesh of planetary gear

内齿圈与行星轮之间的啮合位移为

δrn=-xrsinψrncosβb-yrcosψrncosβb+
zrsinβb+rmrθrxsinψnsinβb-
rmrθrycosψnsinβb-rrθrzcosβb+
xnsinα′trcosβb+yncosα′trcosβb-
znsinβb+rnθnysinβb/cosα′tr+rnθnzcosβb

(15)

式(15)中:xr、yr、zr为内齿圈的平动位移;θrx、θry、θrz为内齿圈的转角位移;rmr为内齿圈的节圆半径;rr为内齿圈的基圆半径;α′tr为内齿圈与行星轮啮合的端面压力角。

1.3.3 行星齿轮动力学模型

求得行星齿轮内、外啮合位移后,分别对行星齿轮各个部件进行受力分析得到部件的运动微分方程。将行齿轮各个部件的运动微分方程按照对应自由度顺序组装,可以得到行星齿轮整体的运动微分方程为

T(t)+F(t)+c

(16)

式(16)中:M为行星齿轮的质量矩阵;Ωc为行星架的转动角速度;G为行星齿轮的回转矩阵;Kb为行星齿轮的支承刚度矩阵;Km为行星齿轮的啮合刚度矩阵;KΩ为行星齿轮的向心刚度矩阵;T为行星齿轮的外力矩阵;F为行星齿轮的内部激励;c为行星齿轮离心加速度矩阵。

行星齿轮特性矩阵的具体表达式见文献[8]。

2 双馈风电机组传动链弯曲-扭转-轴向耦合模型

将双馈风电机组传动链部件动力学模型按照它们在传动链的位置关系进行组装,即可得到双馈风电机组传动链弯曲-扭转-轴向耦合模型。本文研究以某型兆瓦级风力机传动链为例,建立了其传动链耦合模型,如图5所示。该风力机风轮直径70 m,运行转速为10～20.46 r/min,采用了一级行星齿轮与两级平行齿轮的组合,齿轮箱整体的传动比为87.95。依据转子有限元法将传动链划分为31个节点。

图5 风力机传动链耦合模型Fig.5 Dynamic model of drivetrain

3 双馈风电机组气动与控制计算模型

3.1 双馈风电机组气动载荷计算模型

风电机组气动载荷计算方法有传统动量-叶素理论即BEM(blade element momentum)方法、广义动态尾流方法[11]。BEM方法由于其物流意义明确,易于修正,因此本文研究以传统动量-叶素理论的方法为基础,考虑了叶尖与叶根损失修正模型、尾流修正模型以及Leishman-Beddoes动态失速模型建立了风电机组气动载荷计算程序,其计算流程如图6所示。

3.2 双馈风电机组控制模型

风电机组控制主要分为额定以下的转矩控制与额定以上的变桨控制。当风电机组工作在额定风速以下时,此时变桨角保持在最优桨角,风电机组通过调节发电机转矩保持风轮以最优叶尖速比运行,达到最大的功率系数。当风电机组达到额定转速后,通过变桨控制保持功率恒定,风电机组变桨控制器可以采用PID控制器实现。风电机组转矩-转速控制曲线如图7所示。

图7 风力机转矩转速曲线[12]Fig.7 Torque-versus-speed curve[12]

将建立的双馈风电机组传动链弯曲-扭转-轴向动力学模型与气动载荷计算程序、转速变桨控制结合便可以对双馈风电机组传动链动态响应进行仿真计算。

4 双馈风电机组传动链稳态工况仿真

利用建立的传动链动态响应仿真程序对传动链稳态工况(均匀风速12 m/s,无湍流)的动态响应进行仿真。此时机组工作在满发状态,发电机转速为额定转速1 800 r/min,变桨控制系统通过调整变桨角保持风电机组功率恒定,此时的变桨角为10°。

计算得到风电机组稳态工况各级齿轮的动态啮合力如图8所示。从图8可以看出,传动链齿轮啮合力为单边变化,表明在风电机组稳定运行过程中,代表啮合模型的弹簧处于受压状态,产生了啮合力。行星齿轮内啮合力变化趋势与外啮合力相似。外啮合力的振动幅值略大于内啮合力的幅值变化。平行齿轮的啮合力包含了齿轮啮合的高频振动,一级平行齿轮啮合力幅值要比二级平行齿轮啮合力大,相对行星齿轮的啮合力,平行齿轮的啮合力要小一个数量级,即风电机组传动链载荷沿传动链传递过程中呈现出载荷递减的趋势。

图8 稳态工况传动链各级齿轮啮合力Fig.8 Gear meshing force of drivetrain under steady condition

风电机组传动链中包含了很多的轴承,风电机组运行过程中产生的振动会增大传动链中轴承的动载荷,影响轴承的使用寿命。根据GB/T6391—2010标准[13]中的方法对风电机组传动链所有轴承的当量动载荷进行了计算,计算结果如图9所示。

图9 传动链轴承载荷Fig.9 Bearing load of drivetrain

由图9可见,风电机组传动链所有轴承中,行星轮轴承的当量动载荷最大,时域主要表现为低频振动。过大的轴承动载荷会增大轴承的疲劳载荷,引起轴承磨损、剥落,影响风电机组的安全稳定运行,图10为某双馈风电机组传动链发生的轴承故障与齿轮故障。

图10 传动链轴承与齿轮故障Fig.10 Bearing and gear failure of drivetrain

5 双馈风电机组传动链低电压穿越响应分析

国家电网电力科学研究院在其低电压穿越能力测试规程中规定了需要测试风电机组在大功率即功率在0.9倍额定功率以上以及小功率即功率在0.1～0.3倍额定功率之间发生电压跌落的响应特性。测试内容既包括了并网点电压与电流等电学参数,同时也包含了风速曲线、发电机转速曲线以及变桨曲线等运行参数[14]。

将该型风电机组低电压穿越过程中测试得到风速、转速以及变桨曲线加载到风电机组传动链动态响应程序中,即可分析低电压穿越过程对传动链齿轮与轴承载荷的影响。由于机组振动变化的规律性,本文选取机组的两个极端电压跌落工况即最大跌落和最小跌落进行对比分析,如表1所示。

表1 低电压穿越测试工况[15]Table 1 Test conditions of LVRT[15]

5.1 工况1与工况3的对比

工况1与工况3电压跌落幅值均为跌落至90%,跌落时间也相同,均为2 000 ms,只有跌落时的机组工况不同。工况1在发生电压跌落时,机组工作在满发状态,工况3机组工作在偏载状态,工况1与工况3的风速曲线与电压跌落曲线如图11所示。

图11 工况1与工况3的风速曲线与电压跌落曲线Fig.11 Wind speed and voltage sag curve under test 1 and 3

当电压在4 s发生电压跌落时,工况1由于机组处于满发状态,机组控制系统通过变速变桨调节机组功率平衡,而工况3由于工作在偏载状态仅通过转速调节维持平衡,工况1与工况3的转速与变桨调节曲线如图12所示。

图12 工况1与工况3风电机组转速曲线与变桨曲线Fig.12 Curve of rotate speed and pich under test 1 and 3

将工况1与工况3发生跌落的风速曲线、转速与变桨曲线加载到传动链动态响应仿真程序中,计算得到两种工况行星轮外啮合力与行星轮轴承力的时域响应曲线如图13所示。由图13可见,工况1产生的冲击载荷明显高于工况2,其中行星齿轮外啮合力振动的峰值为工况3的5倍左右,行星轮轴承力的峰值工况1为工况3的6倍左右。

图13 工况1与工况3传动链齿轮与轴承动态响应Fig.13 Dynamic response of gear and bearing under test 1 and 3

5.2 工况1与工况2的对比

工况1与工况2机组均工作在满发状态,但是跌落的幅值不同,工况1电压跌落至90%,工况2电压跌落至20%。工况1与工况2的风速曲线以及工况2的电压跌落曲线如图14所示。

图14 工况1与工况2的风速曲线与电压跌落曲线Fig.14 Wind speed and voltage sag curve under test 1 and 2

工况1与工况2机组都工作在满发状态,因此当机组发生跌落后,机组都是通过变速与变桨同时调节,工况2由于电压跌落的幅值比较大,因此转速与变桨调节幅度相比工况1也大。工况1与工况2的转速和变桨曲线如图15所示。

图15 工况1与工况2风电机组转速曲线与变桨曲线Fig.15 Curve of rotate speed and pich under test 1 and 2

计算得到工况1与工况2行星齿轮外啮合力与行星轮轴承力的动态响应如图16所示。从图16可以看出工况1与工况2都产生了明显的冲击效应,两者动态响应的趋势接近。对两者的动态响应进行统计,工况2行星齿轮外啮合力的峰值是工况1的1.38倍,而行星轮轴承力振动的峰值工况2是工况1的1.18倍。

图16 工况1与工况2传动链齿轮与轴承动态响应Fig.16 Dynamic response of gear and bearing under test 1 and 2

5.3 工况2与工况5的对比

工况2和工况5机组都是满发状态,电压都跌落至20%,跌落时间均为625 ms,只有电压跌落的形式不同,工况2为电压三相对称跌落,工况5为不对称跌落。工况2与工况5的风速与电压跌落曲线如图17所示。

图17 工况2与工况5的风速曲线与电压跌落曲线Fig.17 Wind speed and voltage sag curve under test 2 and 5

工况2与工况5机组都是通过转速与变桨同时调节应对电压跌落引起的功率不平衡,两者的转速变桨曲线如图18所示。由于只是电压跌落的形式不同,工况2与工况5的调节曲线非常接近。

图18 工况2与工况5风电机组转速曲线与变桨曲线Fig.18 Curve of rotate speed and pich under test 2 and 5

将工况2与工况5的风速曲线与转速变桨曲线加载到传动链动态响应仿真程序中,对两组工况动态响应进行仿真,同样提取行星齿轮外啮合力与行星轮轴承力作为对比,对比结果如图19所示。由图19可见,工况2与工况5行星轮外啮合力与行星轮轴承力响应很接近,对振动数据进行分析可以得出行星齿轮外啮合力响应的峰值工况2是工况5的1.10倍,工况2行星轮轴承力响应的峰值是工况5的1.01倍。

图19 工况2与工况5传动链齿轮与轴承动态响应Fig.19 Dynamic response of gear and bearing under test 2 and 5

对以上仿真结果进行小结可知,在风电机组发生电压跌落时,导致发电机电磁转矩变小,风电机组会通过转速与变桨控制调节功率平衡。在此过程中会对风电机组传动链载荷产生冲击。针对行星齿轮外啮合力与行星轮轴承力,这种冲击与发电跌落的机组工况、电压跌落的幅值以及电压跌落的形式相关。具体表现为,跌落时机组载荷越大,电压跌落的幅值越大,产生的冲击载荷越大。另外电压三相对称跌落相比不对称跌落产生的冲击载荷大。因此,所有工况中工况2的跌落工况冲击产生的载荷最大。

6 结论

(1)根据双馈风电机组传动链结构特点,分别采用转子有限元法与集中建模方法对双馈风电机组传动链中的盘、轴、轴承、联轴器以及平行齿轮与行星齿轮等结构进行了动力学建模,并将这些部件动力学模型按照实际的装配关系进行了组装得到传动链整体的弯曲-扭转-轴向耦合的动力学模型。该模型能够对风电机组中的关键部件动力学特性与载荷特性进行准确仿真。

(2)将建立的双馈风电机组传动链耦合动力学模型结合风电机组气动载荷计算模型与转速变桨控制模型形成传动链动态响应仿真程序,并对传动链稳态工况进行了计算。稳态工况中,行星齿轮与平行齿轮均为受载状态,其中行星齿轮外啮合力为齿轮系中啮合力最大部件。对传动链中所有轴承当量动载荷进行了计算,行星轮轴承为所有轴承动载荷最大部件。表明行星齿轮的设计对整个传动链的承载特性至关重要。

(3)将某型兆瓦级风电机组低电压穿越测试的风速、转速以及变桨运行曲线加载到传动链动态响应仿真程序中,对低电压穿越各类工况的动态响应进行了对比分析。结果表明,低电压穿越过程中会对风电机组传动链的齿轮与轴承产生冲击,冲击载荷的大小与机组跌落的工况、跌落的幅值大小以及跌落的形式相关。具体为机组工况载荷越大、电压跌落幅值越大,产生的冲击载荷越大,三相对称跌落产的冲击载荷要比不对称跌落大,并且这三个因素的影响是依次递减的。因此,在机组设计过程中应特别关注机组满发,电压发生三相对称最大跌落值的载荷工况。