伴随着国家支持和环境对发展清洁能源的需求，风力发电市场发展迅猛，2013年中国新增装机容量16 088 MW，累计装机容量91 412 MW，新增装机和累积装机两项数据均居世界第一。随着风力发电市场的发展，风机频发的故障成了主要待解决的问题，提高风机设计可靠性成了一个风机的主要研究方向^[1]。

风力发电机组主要由气动系统、传动链系统、电气系统和控制系统组成，在风力发电机发电的过程中，其传动链往往受到时变风载，紧急停车，塔筒的晃动以及冬夏温差等因素的影响，这些因素都影响着传动链的工作性能，导致传动系统故障频发。目前国内外学者对风电机组传动系统进行了一定的研究，胥良考虑极限和额定载荷工况，对兆瓦级风电齿轮箱的关键零部件进行静力学分析，得到齿轮箱箱体、行星架的应力云图以及极限载荷系数^[2]。陈爽利用campell分析、模态能量分布分析和阶次扫频分析等手段对风电齿轮箱传动系统潜在共振进行了分析^[3]。徐向阳利用动力学方法，对柔性销轴式风电齿轮箱的均载、耦合特性与动态特性等展开了研究^[4]。朱才朝等人建立了风电齿轮箱外啮合和内啮合的有限元模型，提出了一种计算齿轮啮合过程中实际接触齿数、齿间载荷分配与齿面载荷分配的方法，并对风电齿轮箱的共振特性展开相关研究^[5-6]；秦大同等人针对风力发电机对齿轮箱的高可靠性要求，提出了一种加速疲劳试验技术方法来验证齿轮箱的疲劳程实时在线测试和评价^[7]。综上所述，现有文献在对风机的理论、风场测试和台架实验进行研究时，均忽略了塔筒的摆动与扭转对风电机组传动链的影响。研究以主流兆瓦级风点机组传动链结构为研究对象，考虑风机塔筒的摆动与风机传动链的扭转振动耦合效应，利用集中参数法对风机传动链动态特性进行研究。

1 兆瓦级风电机组传动链结构与原理

目前风机主要有直驱型风机、半直驱型风机和双馈型风机，其中以双馈型为主。国内各风场以水平式的风机结构为主，如图 1所示为水平式某兆瓦风力发电机组传动链结构简图，机舱位于近100 m的塔筒顶端，齿轮箱位于塔筒中部，通过主轴轴承与扭力臂固定在机舱中，利用叶片将风能转化为机械能。叶片固定在轮毂上并通过轮毂带动主轴转动，主轴通过锁紧盘与齿轮箱行星架相连接，通过齿轮箱升速将主轴上的低转速大扭矩动能转化为高转速低扭矩的动能，带动发电机进行发电。

图 2所示为风电齿轮箱结构简图，主轴通过锁紧盘带动行星架转动，行星轮安装在行星架上，同时与内齿圈和太阳轮啮合，进而带动太阳轮转动，太阳轮与低速轴为一体设计，低速轴通过花键带动中间级主动轮转动，第二级和第三极为大齿轮与小齿轮进行啮合升速，最终将主轴的低转速大扭矩动能转化为适应发电机发电频率的高转速低扭矩的动能。

2 兆瓦级风电机组传动链动力学建模

某三点支撑的兆瓦级风电机组传动链连接拓扑关系如图 3所示，其中叶片根部与轮毂固定连接，轮毂与主轴连接，主轴利用锁紧盘与齿轮箱上的行星架相连接，同时通过主轴轴承支撑在机舱上，齿轮箱通过左右2个扭力臂支撑在机舱上，同时通过齿轮高速轴与联轴器连接，联轴器与发电机转子连接，发电机定子通过轴承支撑机舱上，机舱由塔架支撑。

齿轮箱是传动链的重要部件，其结构和受力较为复杂。为了对传动链的动态特性进行深入分析，将图 3中的齿轮箱内部连接拓扑关系进行细化，如图 4所示为齿轮箱内部结构拓扑关系图。风力发电机拓扑图符合含义如表 1所示，齿轮箱拓扑图符合含义如表 2所示。

由于叶片连接在轮毂上，不考虑变桨时，可以将其考虑为固定连接，而轮毂固定于主轴前端，因此将叶片与轮毂的质量和转动惯量换算到主轴上，根据上述兆瓦级风电机组传动链拓扑图模型和齿轮箱拓扑关系模型，利用拉格朗日法建立系统动力学模型。其中系统动能为

$\begin{array}{c} T = 1/2\left( {{I_{mx}}\dot \theta _{mx}^2 + {I_{my}}\dot \theta _{my}^2 + {I_{mz}}\dot \theta _{mz}^2 + {I_{cx}}\dot \theta _{cx}^2 + {I_{cy}}\dot \theta _{cy}^2 + {I_{cz}}\dot \theta _{cz}^2 + } \right.\\ {I_{sx}}\dot \theta _{sx}^2 + {I_{sy}}\dot \theta _{sy}^2 + {I_{sz}}\dot \theta _{sz}^2 + {I_{rx}}\dot \theta _{rx}^2 + {I_{ry}}\dot \theta _{ry}^2 + {I_{rz}}\dot \theta _{rz}^2 + \\ \sum\limits_{i = 1}^N {{I_{piu}}\dot \theta _{piu}^2 + \sum\limits_{i = 1}^N {{I_{piv}}\dot \theta _{piv}^2 + \sum\limits_{i = 1}^N {{I_{piw}}\dot \theta _{piw}^2 + + {I_{2wx}}\dot \theta _{2wx}^2 + } } } \\ {I_{2wy}}\dot \theta _{2wy}^2 + {I_{2wz}}\dot \theta _{2wz}^2 + {I_{2px}}\dot \theta _{2px}^2 + {I_{2py}}\dot \theta _{2py}^2 + {I_{2pz}}\dot \theta _{2pz}^2 + \\ {I_{3wx}}\dot \theta _{3wx}^2 + {I_{3wy}}\dot \theta _{3wy}^2 + {I_{3wz}}\dot \theta _{3wz}^2 + {I_{3px}}\dot \theta _{3px}^2 + {I_{3py}}\dot \theta _{3py}^2 + \\ {I_{3pz}}\dot \theta _{3pz}^2 + {I_{gx}}\dot \theta _{gx}^2 + {I_{gy}}\dot \theta _{gy}^2 + {I_{gz}}\dot \theta _{gz}^2 + {I_{py}}\dot \theta _{py}^2 + {I_{pz}}\dot \theta _{pz}^2\\ {m_m}\dot x_n^2 + {m_m}\dot y_m^2 + {m_m}\dot z_m^2 + {m_c}\dot x_c^2 + {m_c}\dot y_c^2 + {m_c}\dot z_c^2 + {m_s}\dot x_s^2 + \\ {m_s}\dot y_s^2 + {m_s}\dot z_s^2 + {m_r}\dot x_r^2 + {m_r}\dot y_r^2 + {m_r}\dot z_r^2 + \sum\limits_{i = 1}^N {{m_{pi}}\dot u_{pi}^2 + } \\ \sum\limits_{i = 1}^N {{m_{pi}}\dot v_{pi}^2 + \sum\limits_{i = 1}^N {{m_{pi}}\dot w_{pi}^2 + } } {m_{2w}}\dot x_{2w}^2 + {m_{2w}}\dot y_{2w}^2 + {m_{2w}}z\dot x_{2w}^2 + \\ {m_{2p}}\dot x_{2p}^2 + {m_{2p}}\dot y_{2p}^2 + {m_{2p}}\dot z_{2p}^2 + {m_{3w}}\dot x_{3w}^2 + {m_{3w}}\dot y_{3w}^2\\ {m_{3w}}\dot z_{3w}^2 + {m_{3p}}\dot x_{3p}^2 + {m_{3p}}\dot y_{3p}^2 + {m_{3p}}\dot z_{3p}^2 + {m_g}\dot x_g^2 + \\ \left. {{m_g}\dot y_g^2 + {m_g}\dot z_g^2 + {m_p}\dot x_p^2 + {m_p}\dot z_p^2 + {m_t}\dot x_t^2 + {m_t}\dot z_t^2} \right), \end{array}$

(1)

式中：T为系统动能；下标m表征传动主轴；c表征行星架；s表征太阳齿轮；r表征齿圈；pi表征第i个行星齿轮；2w表征齿轮箱中间级大齿轮；2p表征齿轮箱中间级小齿轮；3w表征齿轮箱高速级大齿轮；3p表征齿轮箱高速级小齿轮；g表征齿轮箱箱体；p表征齿轮箱支撑；t表征塔筒；u、v、w为行星齿轮坐标系，x，y，z为坐标系方向自由度，θ_x，θ_y，θ_z为扭转方向自由度，I为转动惯量，m为构建质量。

系统势能

$U = {U_{{\rm{spi}}}} + {U_{{\rm{rpi}}}} + {U_{\rm{b}}} + {U_{2wp}} + {U_{3wp}},$

(2)

式中：U_spi为太阳行星齿轮啮合势能；U_rpi为齿圈行星齿轮啮合势能；U_b为轴承支撑势能；U_2wp为中间级齿轮啮合势能；U_3wp为高速级齿轮啮合势能。

其中

${U_{{\rm{spi}}}} = 1/2\sum\limits_{i = 1}^N {{K_{{\rm{spi}}}}\delta _{{\rm{spi}}}^2, }$

(3)

${U_{{\rm{rpi}}}} = 1/2\sum\limits_{i = 1}^N {{K_{{\rm{rpi}}}}\delta _{{\rm{rpi}}}^2, }$

(4)

$\begin{array}{c} {U_b} = 1/2\left( {{K_{m\theta x}}\theta _{mx}^2 + {K_{m\theta y}}\theta _{my}^2 + {K_{mx}}x_m^2 + {K_{my}}y_m^2 + {K_{mz}}z_m^2} \right) + \\ 1/2\left( {{K_{c\theta x}}\theta _{cx}^2 + {K_{c\theta y}}\theta _{cy}^2 + {K_{cx}}x_c^2 + {K_{cy}}y_c^2 + {K_{cz}}z_c^2} \right) + \\ 1/2\left( {{K_{s\theta x}}\theta _{sx}^2 + {K_{s\theta y}}\theta _{sy}^2 + {K_{sx}}x_s^2 + {K_{sy}}y_s^2 + {K_{sz}}z_c^2} \right) + \\ 1/2\sum\limits_{i = 1}^N {\left( {{K_{piu}}\theta _{piu}^2 + {K_{piv}}\theta _{piv}^2 + {K_{piu}}\delta _{piu}^2 + {K_{piv}}\delta _{piv}^2 + {K_{piw}}\delta _{piw}^2} \right) + } \\ 1/2\left( {{K_{r\theta x}}\theta _{rx}^2 + {K_{r\theta y}}\theta _{ry}^2 + {K_{r\theta z}}\theta _{rz}^2 + {K_{rx}}x_r^2 + {K_{ry}}y_r^2 + {K_{rz}}z_r^2} \right) + \\ 1/2\left( {{K_{2w\theta x}}\theta _{2wx}^2 + {K_{2w\theta y}}\theta _{2wy}^2 + {K_{2wx}}x_{2w}^2 + {K_{2wy}}y_{2w}^2 + {K_{2wz}}z_{2w}^2} \right) + \\ 1/2\left( {{K_{2p\theta x}}\theta _{2px}^2 + {K_{2p\theta y}}\theta _{2py}^2 + {K_{2px}}x_{2p}^2 + {K_{2py}}y_{2p}^2 + {K_{2pz}}z_{2p}^2} \right) + \\ 1/2\left( {{K_{3w\theta x}}\theta _{3wx}^2 + {K_{3w\theta y}}\theta _{3wy}^2 + {K_{3wx}}x_{3w}^2 + {K_{3wy}}y_{3w}^2 + {K_{3wz}}z_{3w}^2} \right) + \\ 1/2\left( {{K_{3p\theta x}}\theta _{3px}^2 + {K_{3p\theta y}}\theta _{3py}^2 + {K_{3px}}x_{3p}^2 + {K_{3py}}y_{3p}^2 + {K_{3pz}}z_{3p}^2} \right) + \\ 1/2\left( {{K_{g\theta x}}\theta _{gx}^2 + {K_{g\theta y}}\theta _{gy}^2 + {K_{g\theta z}}\theta _{gz}^2 + {K_{gx}}x_g^2 + {K_{gy}}y_g^2 + {K_{gz}}z_g^2} \right) + \\ 1/2\left( {{K_{p\theta y}}\theta _{py}^2 + {K_{px}}x_p^2 + {K_{pz}}z_p^2} \right) + 1/2\left( {{K_{t\theta y}}x_{ty}^2 + {K_{tx}}x_t^2 + {K_{tz}}z_t^2} \right), \end{array}$

(5)

${U_{2wp}} = 1/2 \cdot {K_{2wp}}\delta _{2wp}^2,$

(6)

${U_{3wp}} = 1/2 \cdot {K_{3wp}}\delta _{3wp}^2,$

(7)

式中，K_spi为太阳行星齿轮啮合刚度；δ_spi为太阳行星齿轮啮合向等效振动位移^[8-11]。K_rpi为齿圈行星齿轮啮合刚度；δ_rpi为齿圈行星齿轮啮合向等效振动位移；K_mθx为主轴x轴扭转向支撑刚度；δ_2wp为中间级齿轮啮合向等效振动位移；δ_3wp为高速级齿轮啮合向等效振动位移。

系统耗散能

$R = {R_{{\rm{spi}}}} + {R_{{\rm{rpi}}}} + {R_b} + {R_{2wp}} + {R_{3wp}},$

(8)

式中：R_spi为太阳行星齿轮啮合耗散能；R_rpi为齿圈行星齿轮啮合耗散能；R_b为系统支撑耗散能；R_2wp为中间级齿轮啮合耗散能；R_3wp为高速级齿轮啮合耗散能。

其中

${R_{{\rm{spi}}}} = 1/2\sum\limits_{i = 1}^N {{C_{{\rm{spi}}}}\delta _{spi}^2, }$

(9)

${R_{{\rm{rpi}}}} = 1/2\sum\limits_{i = 1}^N {{C_{{\rm{rpi}}}}\delta _{{\rm{rpi}}}^2, }$

(10)

$\begin{array}{c} {R_b} = 1/2\left( {{C_{m\theta x}}\dot \theta _{mx}^2 + {C_{m\theta y}}\dot \theta _{my}^2 + {C_{mx}}\dot x_m^2 + {C_{my}}\dot y_m^2 + {C_{mz}}\dot z_m^2} \right) + \\ 1/2\left( {{C_{c\theta x}}\dot \theta _{cx}^2 + {C_{c\theta y}}\dot \theta _{cy}^2 + {C_{cx}}\dot x_c^2 + {C_{cy}}\dot y_c^2 + {C_{cz}} + \dot z_c^2} \right) + \\ 1/2\left( {{C_{s\theta x}}\dot \theta _{sx}^2 + {C_{s\theta y}}\dot \theta _{sy}^2 + {C_{sx}}\dot x_s^2 + {C_{sy}}\dot y_s^2 + {C_{sz}} + \dot z_s^2} \right) + \\ 1/2\sum\limits_{i = 1}^N {\left( {{C_{piu}}\dot \theta _{piu}^2 + {C_{piv}}\dot \theta _{piv}^2 + {C_{piu}}\dot \delta _{piu}^2 + {C_{piv}}\dot \delta _{piv}^2 + {C_{piw}}\delta _{piw}^2} \right) + } \\ 1/2\left( {{C_{r\theta x}}\dot \theta _{rx}^2 + {C_{r\theta y}}\dot \theta _{ry}^2 + {C_{r\theta z}}\dot \theta _{rz}^2 + {C_{rx}}\dot x_r^2 + {C_{ry}}\dot y_r^2 + {C_{rz}}\dot z_r^2} \right) + \\ 1/2\left( {{C_{2w\theta x}}\dot \theta _{2wx}^2 + {C_{2w\theta y}}\dot \theta _{2wy}^2 + {C_{2wx}}\dot x_{2w}^2 + {C_{2wy}}\dot y_{2w}^2 + {C_{2wz}}\dot z_{2w}^2} \right) + \\ 1/2\left( {{C_{2p\theta x}}\dot \theta _{2px}^2 + {C_{2p\theta y}}\dot \theta _{2py}^2 + {C_{2px}}\dot x_{2p}^2 + {C_{2py}}\dot y_{2p}^2 + {C_{2pz}}\dot z_{2p}^2} \right) + \\ 1/2\left( {{C_{3w\theta x}}\dot \theta _{3wx}^2 + {C_{3w\theta y}}\dot \theta _{3wy}^2 + {C_{3wx}}\dot x_{3w}^2 + {C_{3wy}}\dot y_{3w}^2 + {C_{3wz}}\dot z_{3w}^2} \right) + \\ 1/2\left( {{C_{3p\theta x}}\dot \theta _{3px}^2 + {C_{3p\theta y}}\dot \theta _{3py}^2 + {C_{3px}}\dot x_{3p}^2 + {C_{3py}}\dot y_{3p}^2 + {C_{3pz}}\dot z_{3p}^2} \right) + \\ 1/2\left( {{C_{g\theta x}}\dot \theta _{gx}^2 + {C_{g\theta y}}\dot \theta _{gy}^2 + {C_{g\theta z}}\dot \theta _{gz}^2 + {C_{gx}}\dot x_g^2 + {C_{gy}}\dot y_g^2 + {C_{gz}}\dot z_g^2} \right) + \\ 1/2\left( {{C_{p\theta y}}\dot \theta _{py}^2 + {C_{px}}\dot x_p^2 + {C_{pz}}\dot z_p^2} \right) + 1/2\left( {{C_{t\theta y}}\dot x_{ty}^2 + {C_{tx}}\dot x_t^2 + {C_{tz}}\dot z_t^2} \right), \end{array}$

(11)

${R_{2pw}} = 1/2 \cdot {C_{2pw}}\dot \delta _{2pw}^2,$

(12)

${R_{3wp}} = 1/2 \cdot {C_{3wp}}\delta _{3wp}^2,$

(13)

式中，C_spi为太阳行星齿轮啮合阻尼；${\dot \delta }$_spi为太阳行星齿轮啮合线方向等效振动速度。C_rpi为齿圈行星齿轮啮合阻尼；${\dot \delta }$_rpi为齿圈行星齿轮啮合线方向等效振动速度。C_mθx为主轴x轴扭转向支撑阻尼。C_2wp为中间级齿轮啮合阻尼；${\dot \delta }$_2wp为中间级齿轮啮合线方向等效振动速度。C_3wp为高速级齿轮啮合阻尼；为高速级齿轮啮合线方向等效振动速度。

根据拉格朗日方程，定义拉格朗日算子L=T-U，对系统各自由度求偏微分

$\left\{ \begin{array}{l} \frac{{\rm{d}}}{{{\rm{d}}t}}\left( {\frac{{\partial L}}{{\partial \dot e}}} \right)-\frac{{\partial L}}{{\partial e}} + \frac{{\partial R}}{{\partial \dot e}} = {F_e}\\ \frac{{\rm{d}}}{{{\rm{d}}t}}\left( {\frac{{\partial L}}{{\partial \dot f}}} \right)-\frac{{\partial L}}{{\partial f}} + \frac{{\partial R}}{{\partial \dot f}} = {F_f}\\ \frac{{\rm{d}}}{{{\rm{d}}t}}\left( {\frac{{\partial L}}{{\partial \dot g}}} \right)-\frac{{\partial L}}{{\partial g}} + \frac{{\partial R}}{{\partial \dot g}} = {F_g}\\ \frac{{\rm{d}}}{{{\rm{d}}t}}\left( {\frac{{\partial L}}{{\partial \dot h}}} \right) - \frac{{\partial L}}{{\partial h}} + \frac{{\partial R}}{{\partial \dot h}} = {F_h}\\ \frac{{\rm{d}}}{{{\rm{d}}t}}\left( {\frac{{\partial L}}{{\partial \dot l}}} \right) - \frac{{\partial L}}{{\partial l}} + \frac{{\partial R}}{{\partial \dot l}} = {F_l} \end{array} \right.,$

(14)

其中，e为传动主轴自由度，为6方向自由度；f为齿轮箱传动系统自由度，共66个自由度；g为齿轮箱箱体自由度，为6方向自由度；h为机舱自由度，为3方向自由度(x，z及θy)；l为塔筒自由度，为3方向自由度(x，z及θy)。整理求导后系统各方向方程，可得基于柔性支撑的传动链耦合系统动力学方程

$\boldsymbol{Mx + Cx + Kx} = \boldsymbol{F}\left( t \right)$

(15)

式中，M为系统质量矩阵，$\boldsymbol{x}\left( {\ddot x, \dot x} \right)$为系统振动位移(速度，加速度)矩阵；C为系统阻尼矩阵；K为系统刚度矩阵；F为系统外部载荷矩阵，其中发电机考虑为系统扭矩负载。

3 兆瓦级风电机组传动链动态特性分析

表 3、4所示为某3点支撑兆瓦级风电机组系统及齿轮箱参数。

忽略阻尼的影响，对上述方程求解，提取耦合系统前16阶固有频率如表 5所示。

由于主轴与齿轮箱各部件多为旋转部件，这些部件的转频或啮合频率均可能激发传动链的共振，对整个风机是潜在的威胁，需要对传动链进行潜在共振分析。表 6为传动链主要激励频率，包括各轴的转频以及各级啮合频率。

利用风电系统传动链固有频率与系统主要激励频率作Campell图，如图 5所示。从图中可以看出，在转速为1 700.9 r/min时存在共振点A，为中间级2倍啮合频率与第八阶系统固有频率266.3 Hz交点，额定转频裕度为5.03%；在转速为1 695 r/min时，存在共振点B，为高速级啮合频率3倍频与第十六阶固有频率1 950.9 Hz的交点，额定转频裕度为5.3%；在转速为1 649.4 r/min时存在共振点C，为高速级2倍啮合频率与第十五阶系统固有频率1 264.4 Hz交点，额定转频裕度为7.9%。其余共振点离额定转频较远。

为了验证分析模型的正确性，考虑塔筒支撑作用，采用Simpack建立了其柔性支撑传动链刚柔耦合多体动力学模型如图 6所示，得到系统前6阶固有频率如表 7所示，其对应振型如图 7所示。可以看出前面集中参数分析结果与有限元仿真结果基本一致，同时从图 7所示振型可以看出其第二阶固有频率为塔筒的摆动，第四阶固有频率为塔筒的摆动与齿轮箱行星架等转动部件的扭转相耦合。从图 5的(a)、(b)中可以看出，传动链第二阶固有频率和第四固有频率与部分转频存在交点，是风机传动链的潜在共振点，因此进一步证实了考虑塔筒扭摆的必要性。

4 结论

考虑风机塔筒的摆动与风机传动链的扭转耦合效应，采用集中参数法建立了基于柔性塔筒支撑的兆瓦级风电机组传动链耦合振动分析模型，并对传动链系统动态固有特性进行求解，利用Campell图对风机传动链潜在共振点进行了分析，得出切入转速与切出转速之间的潜在共振频率主要为：中间级2倍啮合频率与第八阶系统固有频率交点266.3 Hz，高速级啮合频率3倍频与第十六阶固有频率的交点1 950.9 Hz和高速级2倍啮合频率与第十五阶系统固有频率交点1 264.4 Hz，其余共振点离额定转频较远。利用Simpack建立了风机传动链多体系统动力学模型，仿真分析结果表明集中参数模型与有限元仿真结果基本一致，验证了分析方法的正确性。研究为风机传动链设计提供参考，具有重要的理论意义与工程应用价值。