基于多体动力学的双馈式风机复合主轴传动链建模与仿真分析

韩新月，张林中，王小虎，李英昌，刘嘉琳

（1.国电联合动力技术有限公司，北京 100039；2.北京市风电设备可靠性工程技术研究中心，北京 100039）

0 引言

随着对风电机组相关技术研究的日益深入[1]，越来越多的厂商发现，基于平台计算载荷，通过行业商业软件对关键零部件强度进行校核设计的方法存在较大缺陷。主要原因是，软件中建模为理想化模型且过于简化，与实际运行机组之间存在制造、安装误差，无法对各部件之间的相互作用关系进行精准模拟，造成设计裕度过大。除此之外，通用软件无法输出一些关键部位，特别是传动链内部结构，如齿轮、齿轮箱内轴承等的载荷，为关键部件的设计带来困难。

为解决以上问题，本文就风电机组传动链模型计算方法、模型搭建、仿真、结果对比等方面展开阐述，建立近实物的传动链模型，以期降低机组重量。利用多体动力学模型代替物理样机对其候选设计的各种特性进行测试和评价，以减少研发成本，缩短研发时间。

1 风电机组传动链动力学模型

风电机组多柔体动力学模型的载荷计算是一个非常复杂的过程[2]～[5]，对载荷的精确估计需要采用精确的计算机模拟，计算机模拟的基础是基于数学模型仿真，因此数学模型搭建是设计核心。

1.1 传动链数学模型

利用叠加原理，推导出在惯性坐标系下轮毂的角速度EωH为

式中：EωB为塔顶在惯性坐标系中的角速度；BωN为机舱相对于塔顶的角速度；NωH为轮毂（风轮）相对于机舱的相对角速度。

EωB与塔架的变形相关。

式中：q˙7，q˙8分别为塔筒一阶模态塔顶纵向和横向位移；q˙9，q˙10分别为塔筒二阶模态塔顶纵向和横向位移；H为塔筒高度；φ1T，φ2T分别为塔筒的一、二阶模态；h为对应塔筒段的高度，其值为0～H。

BωN与偏航速率和倾角的速率相关。

式中：q˙5为机舱倾角；q˙6为机舱与风向夹角；d2，d3分别为偏航坐标系中的纵向位移和垂直方向位移。

NωH与低速轴的角速度相关。

式中：q˙4为风电机组风轮的方位角；e1为低速轴方位坐标系中的轴，通常指叶片1在0方位角的方向。

轮毂的角速度可以改写为

如果忽略塔架轴向的变形，弹性体塔架任意点的线速度EVT为

式中：h为塔架上任意点到塔底的高度；a1，a3分别为塔顶坐标系下的横向位移和纵向位移。

机舱质心D在惯性系中的线速度EVD为

式中：EVO为塔顶在惯性坐标系下的线速度；EωD为机舱相对塔筒的角速度；rOD为机舱质心相对塔筒的距离。

低速轴端点P在惯性系中的线速度EVP为

式中：EωN为机舱相对惯性坐标系的角速度；rDP为机舱质心相对低速轴的距离。

轮毂中心点Q在惯性系中的线速度EVQ为

式中：EωH为轮毂相对惯性坐标系的角速度；rPQ为低速轴至叶根的距离。

叶片上任意点S在惯性系中的线速度EVS为

式中：HVS为叶片某段截面相对轮毂的线速度；rQS为叶片某段至叶根的距离。

通过整理，角速度和线速度可以分解为

风力发电机组结构上任意点X的加速度EaX可以在式（12）的基础上直接求导获得。

动力学方程采用Kane方法建立，直接采用牛顿定律，完整系统的动力学Kane方程为

式中：Fr为广义主动力；Fr*为广义惯性力。

将前面部分推导的方程带入，得到的耦合动力学方程可以描述整个风力发电机组结构在风载作用下的力学行为。风力发电机组动力学方程可以转化为

写成矩阵形式为

式中：Crs为广义加速度项的系数；fr为低阶项；qi为广义与柔性体有关的坐标。

风力发电机传动链的基本分析模型是把各个部件都看作惯性元件，它们都具有弹性和阻尼。

系统的转矩和惯量的简单模型为

式中：Jm为齿轮箱低速级的部件转动惯量，由于齿轮箱内部各轴的转动惯量远小于风轮及主轴的转动惯量，所以可以近似看作是风轮转动惯量；Je为高速轴及发电机转子的转动惯量；G为传动比；B为阻尼；ω为角速度；Ta，Te分别为气动转矩和电磁转矩，忽略变流器的动态特性。

一般认为传动机构属于刚性器件，一阶惯性环节即可表示该机构的特性。传动机构的运动方程为

式中：Ta为传动机构输入转矩；Te为传动机构输出转矩；tk为轮毂惯性时间常数。

在简化模型中可将传动轴的惯量等效到发电机转子中，齿轮箱为理想的刚性齿轮组。

1.2 传动链结构及拓扑关系

在传动链模型中，包含了叶片、轮毂、主轴、主轴上的轴承、齿轮箱、电机和机架等风电机组的大部分部件（图1）。

图1 传动链结构Fig.1 Structure of drive train

传动链拓扑作为双馈风电机组的重要部件，是整个风电机组可靠性的重要一环，齿轮箱与主机架、高速轴、低速轴连接，整个传动链的拓扑关系如图2所示。在该模型中，坐标系的x轴平行于主轴，y轴水平，z轴向上。z轴与重力方向成5°角。在图中，接头标记“α”是关于x的旋转方向。

齿轮箱采用2级行星1级平行结构，齿轮箱在子结构中建模，采用力元242（FE242：花键联轴器）和力元5（FE5：弹簧阻尼平行凸轮轴位置）对齿轮箱花键联轴器进行了仿真。测力元件5由高刚度弹簧构成，在x轴向平移方向上的间隙为±0.5 mm。

本文采用有限元模型对齿轮箱中的部分部件进行了柔性体仿真。齿轮箱中的轴承在动力学模型中采用具有刚度矩阵的力单元41（FE41：弹簧-阻尼矩阵Cmp）进行模拟。

1.3 新型主轴优化设计

风力发电机组主轴是传动系统的重要组成部分。主轴前端通过螺栓与风轮刚性连接，中部与轴承连接，后部与齿轮箱低速轴（或电机）连接。主轴承受力大且复杂，受力形式主要为轴向力、径向力、弯矩、转矩和剪切力。目前，机舱重量越来越重，主轴作为机舱的一个重要部件，主轴的重量和性能也是重要设计参数。本文主轴设计采用碳纤维和金属材料的复合主轴，形式如图3所示。图中，d1为主轴复合材料体内径，d2为复合材料外径，d3为主轴金属体外径。复合主轴中间采用复合材料、外侧采用金属34CrNiMo6，内腔采用缠绕式预浸料碳纤维，外侧采用金属。

图3 复合主轴剖面图Fig.3 Compound mainshaft sectional drawing

基于图4的优化流程对主轴重量成本进行优化。在保证主轴满足极限强度和疲劳强度的前提下，优化主轴内径、外径和缠绕角度等参数。本程序假设传动链承受的极限和疲劳载荷是不变的，基于此假设，在Isight优化软件中进行主轴的参数优化，搭建了优化平台（图5）。

图4 复合主轴优化设计流程图Fig.4 Flow chart of composite mainshaft optimization design

图5 复合主轴优化平台Fig.5 Composite mainshaft optimization optimization platform

优化目标：主轴重量和成本。

优化变量：d1，d2，d3和碳纤维缠绕线左右角度。

约束条件：全生命周期内的极限强度和疲劳强度满足IEC61400规范和GL2010规范要求。

d1，d2，d3和碳纤维缠绕线左右角度共计5个参数设为优化变量，以单位成本重量最低为优化目标，约束为疲劳强度校核、极限强度角度等安全性评估，根据图4优化出一款复合材料主轴。

优化程序说明：

①循环次数设为i，自初始i=0开始计算，每循环一次i+1；

②第i次循环对应的尺寸参数为d1i，d2i，d3i，角度参数为γ1i，γ2i，当i=0时，上述参数均为初始值；

③根据极限载荷和疲劳载荷，结合上述的尺寸和材料数据便可进行安全性评估；

④判断是否满足标准设计要求，是的话满足主轴成本和重量；

⑤否则i+1，重复上述步骤②～④；

⑥计算单位成本重量；

⑦yi是否最小，判断依据为|yi+1-yi|/yi＜5%；

⑧若最小，退出循环；

⑨得到主轴的尺寸参数d1i，d2i，d3i，角度参数γ1i，γ2i；

⑩否则继续循环，直至找到满足条件。

在图5中，基于复合材料主轴，对部件进行极限载荷和疲劳载荷校核，两者并行求解计算。其中极限强度校核调用了有限元求解软件ANSYS，疲劳强度校核调用了Ncode软件，针对疲劳和极限的安全系数依据IEC61400规范和GL2010规范要求判断是否满足要求，即为安全裕度的计算，若安全裕度满足要求，则进行重量和成本评估，采用遗传算法进行优化计算。基于上述优化程序，对某一机型主轴进行优化设计，原主轴重量为20 209 kg，优化后为17 512 kg，重量降低了13.35%，成本降低了7.15%。

1.4 传动链仿真模型

利用模态缩减求解传动链各部件中的柔性体[6]，[7]。齿轮箱内部行星架、太阳轮轴、高速轴等模型如图6所示。将这些部件按照拓扑关系装配定位到系统中，分别建出齿轮箱子部件，传动链系统模型如图7所示。

图6 齿轮箱定位及装配关系Fig.6 Gearbox positioning and assembly relationship

图7 传动链系统模型Fig.7 Drive train models

2 仿真结果与分析

2.1 激励频率

风力发电机组作为旋转机械的一类，从风轮至电机整个传动链的旋转便是风电机组的激励来源，故旋转频率作为激励频率。若激励频率与固有频率接近，则极易发生共振，必须根据风轮旋转速度推算传动链各个速度等级的激励频率。风轮转子的工作转速为6.05～12.49 r/min，额定转速为10.864 r/min。根据齿轮箱各级齿数（表1）和风轮转速，推出轮毂转速的1p，2p，3p；根据各级齿轮齿数计算各级转频以及啮频，表2列出了相应的激励频率和啮合频率。

表1 各级齿数Table 1 Number of teeth at all levels

表2 激振频率Table 2 Excitation frequency

由表2可知，第三级齿轮副（二次谐波）啮合力的最高频率为2 438.847 Hz。因此，在共振分析中，所有高于2 438.847 Hz的固有频率均可以忽略。

2.2 固有频率

振动频率与初始条件无关，仅与系统的固有频率有关。固有频率与系统的固有特性有关，固有频率大小会影响到产品稳定性，因此研究子部件和系统的固有频率是产品设计的一个重要步骤。在计算传动链的固有频率之前，已进行时域计算，以保持模型处于额定工作状态。轮毂上的输入负载为2 409.523 kN·m的恒定扭矩。由于平衡，发电机的转速为184.307 rad/s（约1 760 r/min）。表3为不同转速下的固有频率。

表3 不同转速下的固有频率Table 3 Natural frequencies at different speeds

由表3可知，切入速度和切出速度之间的频率变化不明显，相同顺序频率的能量分布非常相似。因此，额定转速下的固有频率可以用来表示传动系模型的动态特性。根据激振频率和传动链频率，分成4个频率范围绘制坎贝尔曲线（图8）。

图8 坎贝尔图Fig.8 Campbell diagram

在传动过程中可以看出齿轮箱动态应变和传动链内部结构动态响应，图9为齿轮箱动态应力及传动链过程动态响应云图。

图9 齿轮箱动态应力及传动链过程动态响应云图Fig.9 Dynamic stress of gearbox and dynamic response of dynamic stress of drive train

2.3 结果对比

根据上述设计复合材料主轴、传动链动力学建模及仿真，以某风况下的机组传动链的位移、速度、加速度和载荷4个方面时序进行对比（表4）。

表4 振动特性对比Table 4 Comparison of vibration characteristics

由表4可知：主轴优化前后与测试曲线趋势一致，且载荷均值、幅值比较一致，说明仿真结果准确性高；优化后的位移、速度、加速度和载荷的极限值均有不同程度的降低，分别降低了1.0%，42.7%，13.8%和2.1%，表明了系统更加稳定。

表5为仿真与实测频率对比。

表5 仿真与实测频率对比Table 5 Comparison of simulated and measured frequencies

由表5可知，多体动力学仿真平台频域仿真结果与实际测试的动态频率平均相差1.51%，体现了仿真模型与实测模型的一致性。

3 结论

风电机组的多体动力学模型在一定程度上具有与物理样机相当的功能真实度。本文通过对风电机组静态及运行时动态的模拟，分析并且显示出实际风电机组传动链及零部件的运动和载荷情况。

①设计了新型复合主轴，编写了复合主轴优化设计程序，有效降低了传动链重量，提升了主轴性能。

②对复合材料主轴进行优化，优化后的传动链位移、速度、加速度和载荷的极限值均有不同程度的降低，分别降低了1.0%，42.7%，13.8%和2.1%，表明了系统更加稳定。

③对仿真结果进行后处理，多体动力学仿真平台频域仿真结果与实际测试的动态频率平均相差1.51%，体现了仿真模型与实测模型的一致性。