大型风电机组是由叶片、轮毂、齿轮传动系统、发电机、机舱、塔架等构件组成的复杂多体系统。风电齿轮箱是在无规律变向载荷的风力作用下运行，振动剧烈、故障频发且维修困难^[1]，因此风电齿轮箱的振动噪声问题受到国内外学者广泛关注。关于齿轮箱内部振动激励及振动导致的辐射噪声，近年来国内外学者已开展了较多研究工作。针对齿轮传动系统的振动，Eritenel等^[2]通过求解行星齿轮系统振动微分方程，得到齿轮系统的非线性振动特性；Liu等^[3]分析了单级齿轮传动摩擦对传动系统动响应的影响；文献[4-6]将制造误差、时变啮合刚度、啮合阻尼、齿侧间隙等因素考虑在齿轮传动系统集中参数法模型中，运用数值方法进行求解；文献[7-8]运用有限元法求解在一定工况下行星齿轮传动系统的动态响应。在齿轮箱辐射噪声方面，文献[9-11]考虑齿轮传动内部激励建立齿轮箱稳态动响应分析模型，采用边界元法进行噪声预估；文献[12-14]采用声振耦合方法分析单级齿轮箱考虑齿轮变刚度激励、传动系统与箱体结合部参数等因素对辐射噪声的影响；李丽君等^[15]基于结构有限元法与声学边界元法研究了振动噪声领域声振强耦合法与弱耦合法问题。

以上文献的研究成果具有重要价值，但是对于多级齿轮箱的辐射振动噪声计算方面，未考虑多级传动系统中各级齿轮副动态传动误差及误差耦合因素。因此，笔者针对风电齿轮箱多级齿轮传动系统的结构特点，综合考虑多种影响因素，建立传动系统的31自由度集中参数模型，采用Runge-Kutta法求解系统微分方程，得到各齿轮副动态传动误差；基于LMS Virtual.Lab建立风电齿轮箱声振耦合模型，研究多级齿轮耦合传动下各齿轮副动态传动误差对齿轮箱振动噪声的影响，并与试验值进行对比验证。

1 传动系统动力学模型分析

文中所研究的齿轮传动系统各级齿轮副基本参数如表 1所示。

针对表 1中的一级行星齿轮及两级平行轴斜齿轮组成的多级齿轮传动系统，采用集中参数法，建立如图 1所示的弯-扭-轴耦合振动分析模型。为简化计算，不考虑齿面摩擦，对于行星直齿轮传动和无轴向力，则仅考虑行星轮、太阳轮、行星架的扭转及径向变形。

齿轮传动系统中，内齿圈固定，其余为活动构件。针对行星架及各齿轮建立独立坐标系，太阳轮、行星架、斜齿轮1~4的坐标原点为各自中心处，坐标方向如图 1中坐标方向所示；行星轮自转坐标系则以行星轮圆心为坐标原点固定在行星架上，考虑切向与径向位移。在定义角位移时，以输入转矩作用下各齿轮转动方向为正，各啮合线上的相对位移方向规定齿面受压方向为正。

1.1 传动系统动力学模型

对于第一级行星齿轮传动，内齿圈固定，行星架输入，太阳轮输出。

设太阳轮与行星轮、行星轮与内齿圈沿啮合点法线方向相对位移分别为x_pis、x_pir(i=1, 2, 3)，则有

$\left\{ \begin{array}{l} {x_{{\rm{p}}i{\rm{s}}}} = {r_{{\rm{bp}}i}}{\theta _{{\rm{p}}i}}-{r_{{\rm{bs}}}}{\theta _{\rm{s}}} + {r_{\rm{c}}}{\theta _{\rm{c}}}\cos \alpha + \left( {{x_{\rm{c}}}-{x_{\rm{s}}}} \right)\cos \left( {{\varphi _i} + \alpha } \right)-\\ \left( {{y_{\rm{c}}} - {y_{\rm{s}}}} \right)\sin \left( {{\varphi _i} + \alpha } \right) + {\eta _{{\rm{p}}i}}\cos \alpha - {\xi _{{\rm{p}}i}}\sin \alpha - {e_{{\rm{p}}i{\rm{s}}}}\left( t \right), \\ {x_{{\rm{p}}i{\rm{r}}}} = {r_{{\rm{bp}}i}}{\theta _{{\rm{p}}i}} - {r_{\rm{c}}}{\theta _{\rm{c}}}\cos \theta - {x_{\rm{c}}}\cos \left( {\alpha + {\varphi _i}} \right) + {y_{\rm{c}}}\sin \left( {\alpha + {\varphi _i}} \right) -\\ {\eta _{{\rm{p}}i}}\cos \alpha + {\xi _{{\rm{p}}i}}\sin \alpha - {e_{{\rm{p}}i{\rm{r}}}}\left( t \right), \end{array} \right.$

(1)

式中：e_pis(t)、e_pir(t) 为法向静态传递误差；x_c、y_c、θ_c为行星架的横向、纵向位移和扭转角；x_s、y_s、θ_s为太阳轮的横向、纵向位移和扭转角；η_pi、ζ_pi、θ_pi为行星轮的切向、径向位移和扭转角；r_bs、r_bpi分别为太阳轮、行星轮的基圆半径，r_c为行星架半径 (即行星轮分布半径)；α为行星级齿轮的压力角；φ_i为第i个行星轮的位置角 (φ_i=2π(i-1)/3，i=1, 2, 3)。

第二级主、从动斜齿面啮合点间因振动和误差产生的沿啮合点法线方向的相对位移x_12n为

$\begin{array}{c} {x_{12{\rm{n}}}} = \left( {{r_{{\rm{b1}}}}{\theta _1}-{r_{{\rm{b2}}}}{\theta _2}} \right)\cos {\beta _{12{\rm{b}}}} + \left( {{z_1}-{z_2}} \right)\sin {\beta _{12{\rm{b}}}} + \\ \left( {{y_1}-{y_2}} \right)\sin {\alpha _{12{\rm{n}}}} + \left( {{x_1} - {x_2}} \right)\cos {\alpha _{12{\rm{t}}}}\cos {\beta _{12{\rm{b}}}} - {e_{12}}\left( t \right), \end{array}$

(2)

式中：β_12b为斜齿轮1、2的基圆螺旋角，α_12n为斜齿轮1、2法面压力角，α_12t为斜齿轮1、2端面压力角；r_b1、r_b2分别为斜齿轮1、2的基圆半径；x₁、y₁、z₁、θ₁为斜齿轮1的横向、纵向、轴向位移和扭转角；x₂、y₂、z₂、θ₂为斜齿轮2的横向、纵向、轴向位移和扭转角；e₁₂(t) 为齿轮副的法向静态传递误差。

第三级主、从动斜齿面啮合点间因振动和误差产生的沿啮合点法线方向的相对位移x_34n为

$\begin{array}{c} {x_{34{\rm{n}}}} = \left( {{r_{{\rm{b3}}}}{\theta _3}-{r_{{\rm{b4}}}}{\theta _4}} \right)\cos {\beta _{34{\rm{b}}}} + \left( {{z_3}-{z_4}} \right)\sin {\beta _{34{\rm{b}}}} + \\ \left( {{y_3}-{y_4}} \right)\sin {\alpha _{34{\rm{n}}}} + \left( {{x_3} - {x_4}} \right)\cos {\alpha _{34{\rm{t}}}}\cos {\beta _{34{\rm{b}}}} - {e_{34}}\left( t \right), \end{array}$

(3)

式中：β_34b为斜齿轮3、4的基圆螺旋角，α_34n为斜齿轮3、4法面压力角，α_34t为斜齿轮3、4端面压力角；r_b3、r_b4分别为斜齿轮3、4的基圆半径；x₃、y₃、z₃、θ₃为斜齿轮3的横向、纵向、轴向位移和扭转角；x₄、y₄、z₄、θ₄为斜齿轮4的横向、纵向、轴向位移和扭转角；e₃₄(t) 为齿轮副的法向静态传递误差。

采用拉格朗日能量法建立风电齿轮传动系统31自由度弯-扭-轴耦合集中参数振动微分方程组，第一级行星传动如式 (4)~式 (6)，第二级斜齿轮传动如式 (7)，第三级斜齿轮传动如式 (8)。

行星架动力学方程组为

$\left\{ \begin{array}{l} {I_{{\rm{cp}}}}{{\ddot \theta }_{\rm{c}}} + {K_{{\rm{c}}\theta }}{\theta _{\rm{c}}} + \sum {{K_{{\rm{sp}}}}\left( t \right){r_{\rm{c}}}\cos \alpha-\sum {{K_{{\rm{pr}}}}\left( t \right){x_{{\rm{p}}i{\rm{r}}}}{r_{\rm{c}}}\cos \alpha + {C_{{\rm{c}}\theta }}{{\dot \theta }_{\rm{c}}}} } \\ \sum {{C_{{\rm{sp}}}}{{\dot x}_{{\rm{p}}i{\rm{s}}}}{r_{\rm{c}}}\cos \alpha-\sum {{C_{{\rm{pr}}}}{{\dot x}_{{\rm{p}}i{\rm{r}}}}{r_{\rm{c}}}\cos \alpha = {T_{{\rm{in}}}}, } } \\ {m_{{\rm{cp}}}}{{\ddot x}_{\rm{c}}} + {K_{{\rm{cr}}}}{x_{\rm{c}}} + \sum {{K_{{\rm{sp}}}}\left( t \right){x_{{\rm{p}}i{\rm{s}}}}\cos \left( {{\varphi _i} + \alpha } \right)}-\sum {{K_{{\rm{pr}}}}\left( t \right){x_{{\rm{p}}i{\rm{r}}}}\cos \left( {{\varphi _i} + \alpha } \right) + {C_{{\rm{c}}x}}{{\dot x}_{\rm{c}}} + } \\ \sum {{C_{{\rm{so}}}}{{\dot x}_{{\rm{p}}i{\rm{s}}}}\cos \left( {{\varphi _i} + \alpha } \right) - \sum {{C_{{\rm{pr}}}}{{\dot x}_{{\rm{p}}i{\rm{r}}}}\cos \left( {{\varphi _i} + \alpha } \right) = 0}, } \\ {m_{{\rm{cp}}}}{{\ddot y}_{\rm{c}}} + {K_{{\rm{c}}y}}{y_{\rm{c}}} - \sum {{K_{{\rm{sp}}}}\left( t \right){x_{{\rm{p}}i{\rm{s}}}}\sin \left( {{\varphi _i} + \alpha } \right) + } \sum {{K_{{\rm{pr}}}}\left( t \right){x_{{\rm{p}}i{\rm{r}}}}\sin \left( {{\varphi _i} + \alpha } \right) + {C_{{\rm{c}}y}}{{\dot y}_{\rm{c}}} - } \\ \sum {{C_{{\rm{sp}}}}{{\dot x}_{{\rm{p}}i{\rm{s}}}}\sin \left( {{\varphi _i} + \alpha } \right) + } \sum {{C_{{\rm{pr}}}}{{\dot x}_{{\rm{p}}i{\rm{r}}}}\sin \left( {{\varphi _i} + \alpha } \right) = 0, } \end{array} \right.$

(4)

式中：I_cp为行星架与行星轮等效转动惯量；m_cp为行星架与行星轮等效质量；K_cj、C_cj(j=x, y) 分别为行星架支撑刚度、支撑阻尼；K_cθ、C_cθ分别为行星架扭转刚度、扭转阻尼；K_sp(t)、K_pr(t) 分别为太阳轮-行星轮、行星轮-齿圈齿轮副时变啮合刚度；C_sp、C_pr分别为太阳轮-行星轮、行星轮-齿圈齿轮副啮合阻尼；T_in为行星架输入扭矩。

太阳轮动力学方程组为

$\left\{ \begin{array}{l} {I_{\rm{s}}}{{\ddot \theta }_{\rm{s}}} + {K_{\theta {\rm{s1}}}}\left( {{\theta _{\rm{s}}}-{\theta _1}} \right)-\sum {{K_{{\rm{sp}}}}\left( t \right){x_{{\rm{p}}i{\rm{s}}}}{r_{{\rm{bs}}}} + {C_{\theta {\rm{s1}}}}\left( {{{\dot \theta }_{\rm{s}}}-{{\dot \theta }_1}} \right) - } \sum {C{_{{\rm{sp}}}}{{\dot x}_{{\rm{p}}i{\rm{s}}}}{r_{{\rm{bs}}}}} = 0\\ {m_{\rm{s}}}{{\ddot x}_{\rm{s}}} + {K_{{\rm{s}}x}}{x_{\rm{s}}} - \sum {{K_{{\rm{sp}}}}\left( t \right){x_{{\rm{p}}i{\rm{s}}}}\cos \left( {{\varphi _i} + \alpha } \right) + {C_{{\rm{s}}x}}{{\dot x}_{\rm{s}}} - \sum {{C_{{\rm{sp}}}}{{\dot x}_{{\rm{p}}i{\rm{s}}}}\cos \left( {{\varphi _i} + \alpha } \right) = 0, } } \\ {m_{\rm{s}}}{{\ddot y}_{\rm{s}}} + {K_{{\rm{s}}y}}{y_{\rm{s}}} + \sum {{K_{{\rm{sp}}}}\left( t \right){x_{{\rm{p}}i{\rm{s}}}}\sin \left( {{\varphi _i} + \alpha } \right) + {C_{{\rm{cy}}}}{{\dot y}_{\rm{s}}} + } \sum {{C_{{\rm{sp}}}}{{\dot x}_{{\rm{p}}i{\rm{s}}}}\sin \left( {{\varphi _i} + \alpha } \right) = 0, } \end{array} \right.$

(5)

式中：I_s为太阳轮转动惯量，m_s为太阳轮质量，K_θs1、C_θs1分别为轴1扭转刚度、阻尼；K_sj、C_sj(j=x, y) 分别为太阳轮横向、纵向支撑刚度、支撑阻尼。

行星轮动力学方程组为

$\left\{ \begin{array}{l} {I_{{\rm{p}}i}}{{\ddot \theta }_{{\rm{p}}i}} + {K_{{\rm{sp}}}}\left( t \right){x_{{\rm{p}}i{\rm{s}}}}\left( {{r_{{\rm{bp}}i}}} \right) + {K_{{\rm{pr}}}}\left( t \right){x_{{\rm{p}}i{\rm{r}}}}\left( {{r_{{\rm{bp}}i}}} \right) + {C_{{\rm{sp}}}}{{\dot x}_{{\rm{p}}i{\rm{s}}}}\left( {{r_{{\rm{bp}}i}}} \right) + {C_{{\rm{sp}}}}{{\dot x}_{{\rm{p}}i{\rm{r}}}}\left( {{r_{{\rm{bp}}i}}} \right) = 0, \\ {m_{{\rm{p}}i}}{{\ddot \xi }_{{\rm{p}}i}}-{K_{{\rm{sp}}}}\left( t \right){x_{{\rm{p}}i{\rm{s}}}}\sin \alpha + {K_{{\rm{pr}}}}\left( t \right){x_{{\rm{p}}i{\rm{r}}}}\sin \alpha + \\ {K_{{\rm{p}}i\xi }}{\xi _i}-{C_{{\rm{sp}}}}{{\dot x}_{{\rm{p}}i{\rm{s}}}}\sin \alpha + {C_{{\rm{pr}}}}{{\dot x}_{{\rm{p}}i{\rm{r}}}}\sin \alpha + {C_{{\rm{p}}i\xi }}{{\dot \xi }_i} = 0, \\ {m_{{\rm{p}}i}}{{\ddot \eta }_{{\rm{p}}i}} + {K_{{\rm{sp}}}}\left( t \right){x_{{\rm{p}}i{\rm{s}}}}\cos \alpha-{K_{{\rm{pr}}}}\left( t \right){x_{{\rm{p}}i{\rm{r}}}}\cos \alpha + \\ {K_{{\rm{p}}i\eta }}{\eta _i} + {C_{{\rm{sp}}}}{{\dot x}_{{\rm{p}}i{\rm{s}}}}\cos \alpha - {C_{{\rm{pr}}}}\left( t \right){x_{{\rm{p}}i{\rm{r}}}}\cos \alpha + {C_{{\rm{p}}i\eta }}{{\dot \eta }_i} - 0, \end{array} \right.$

(6)

式中：I_pi(i=1, 2, 3) 为行星轮转动惯量，m_pi(i=1, 2, 3) 为行星轮质量，K_pij、C_pij(i=1, 2, 3，j=η, ξ) 分别为行星轮支撑刚度、支撑阻尼。

第二级斜齿轮动力学方程组为

$\left\{ \begin{array}{l} {I_1}{{\ddot \theta }_1}-{K_{\theta {\rm{s1}}}}\left( {{\theta _{\rm{s}}}-{\theta _1}} \right) + {K_{12}}\left( t \right){x_{12{\rm{n}}}}{r_{{\rm{b1}}}}\cos {\beta _{12{\rm{b}}}}-{C_{\theta {\rm{s1}}}}\left( {{{\dot \theta }_{\rm{s}}} - {{\dot \theta }_1}} \right) +\\ {C_{12}}{{\dot x}_{12{\rm{n}}}}{r_{{\rm{b1}}}}\cos {\beta _{12{\rm{b}}}} = 0, \\ {m_1}{{\ddot x}_1} + {K_{1x}}{x_1} + {K_{12}}\left( t \right){x_{12{\rm{n}}}}\cos {\alpha _{12{\rm{t}}}}\cos {\beta _{12{\rm{b}}}} + {C_{1x}}{{\dot x}_1} + \\{C_{12}}{{\dot x}_{12{\rm{n}}}}\cos {\alpha _{12{\rm{t}}}}\cos {\beta _{12{\rm{b}}}} = 0, \\ {m_1}{{\ddot y}_1} + {K_{1y}}{y_1} + {K_{12}}\left( t \right){x_{12{\rm{n}}}}\sin {\alpha _{12{\rm{n}}}} + {C_{1y}}{{\dot y}_1} + {C_{12}}{{\dot x}_{12{\rm{n}}}}\sin {\alpha _{12{\rm{n}}}} = 0, \\ {m_1}{{\ddot z}_1} + {K_{1z}}{z_1} + {K_{12}}\left( t \right){x_{12{\rm{n}}}}\sin {\beta _{12{\rm{b}}}} + {C_{1z}}{{\dot z}_1} + {C_{12}}{{\dot x}_{12{\rm{n}}}}\sin {\beta _{12{\rm{b}}}} = 0, \\ {I_2}{{\ddot \theta }_2} - {K_{12}}{x_{12{\rm{n}}}}{r_{{\rm{b2}}}}\cos {\beta _{12{\rm{b}}}} + {K_{\theta 23}}\left( {{\theta _2} - {\theta _3}} \right) - {C_{12}}{{\dot x}_{12{\rm{n}}}}{r_{{\rm{b2}}}}\cos {\beta _{12{\rm{b}}}} + \\{C_{\theta 23}}\left( {{{\dot \theta }_2} - {{\dot \theta }_3}} \right) = 0, \\ {m_2}{{\ddot x}_2} + {K_{2x}}{x_2} - {K_{12}}\left( t \right){x_{12{\rm{n}}}}\cos {\alpha _{12{\rm{t}}}}\cos {\beta _{12{\rm{b}}}} + {C_{2x}}{{\dot x}_2} - \\{C_{12}}{{\dot x}_{12{\rm{n}}}}\cos {\alpha _{12{\rm{t}}}}\cos {\beta _{12{\rm{b}}}} = 0, \\ {m_2}{{\ddot y}_2} + {K_{2y}}{y_2} - {K_{12}}\left( t \right){x_{12{\rm{n}}}}\sin {\alpha _{12{\rm{n}}}} + {C_{2y}}{{\dot y}_2} - {C_{12}}{{\dot x}_{12{\rm{n}}}}\sin {\alpha _{12{\rm{n}}}} = 0, \\ {m_2}{{\ddot z}_2} + {K_{2z}}{z_2} - {K_{12}}\left( t \right){x_{{\kern 1pt} 2{\rm{n}}}}\sin {\beta _{12{\rm{n}}}} + {C_{2z}}{{\dot z}_2} - {C_{12}}{{\dot x}_{12{\rm{n}}}}\sin {\beta _{12{\rm{b}}}} = 0, \end{array} \right.$

(7)

式中：I_i(i=1, 2) 为斜齿轮1、2的转动惯量，m_i(i=1, 2) 为斜齿轮1、2的质量；K_ij、C_ij(i=1, 2；j=x, y, z) 分别为各齿轮支撑刚度、支撑阻尼；K_θ23、C_θ23分别轴2扭转刚度、扭转阻尼；K₁₂(t)、C₁₂分为齿轮副时变啮合刚度、啮合阻尼。

第三级斜齿轮动力学方程组为

$\left\{ \begin{array}{l} {I_3}{{\ddot \theta }_3}-{K_{\theta 23}}\left( {{\theta _2}-{\theta _3}} \right) + {K_{34}}\left( t \right){x_{34{\rm{n}}}}{r_{{\rm{b3}}}}\cos {\beta _{{\rm{34b}}}}-{C_{\theta {\rm{23}}}}\left( {{{\dot \theta }_2} - {{\dot \theta }_3}} \right) + \\{C_{34}}{{\dot x}_{{\rm{34n}}}}{r_{{\rm{b3}}}}\cos {\beta _{34{\rm{b}}}} = 0, \\ {m_3}{{\ddot x}_3} + {K_{3x}}{x_3} + {K_{34}}\left( t \right){x_{34{\rm{n}}}}\cos {\alpha _{34{\rm{t}}}}\cos {\beta _{34{\rm{b}}}} + {C_{3x}}{{\dot x}_3} + \\{C_{34}}{{\dot x}_{34{\rm{n}}}}\cos {\alpha _{34{\rm{t}}}}\cos {\beta _{34{\rm{b}}}} = 0, \\ {m_3}{{\ddot y}_3} + {K_{3y}}{y_3} + {K_{34}}\left( t \right){x_{34{\rm{n}}}}\sin {\alpha _{34{\rm{n}}}} + {C_{3y}}{{\dot y}_3} + {C_{34}}{{\dot x}_{34{\rm{n}}}}\sin {\alpha _{34{\rm{n}}}} = 0, \\ {m_3}{{\ddot z}_3} + {K_{3z}}{z_3} + {K_{34}}\left( t \right){x_{34{\rm{n}}}}\sin {\beta _{34{\rm{b}}}} + {C_{3z}}{{\dot z}_3} + {C_{34}}{{\dot x}_{34{\rm{n}}}}\sin {\beta _{34{\rm{b}}}} = 0, \\ {I_4}{{\ddot \theta }_4} - {K_{34}}{x_{34{\rm{n}}}}{r_{{\rm{b4}}}}\cos {\beta _{34{\rm{b}}}} + {K_{\theta 43}}\left( {{\theta _2} - {\theta _3}} \right) - {C_{34}}{{\dot x}_{34{\rm{n}}}}{r_{{\rm{b4}}}}\cos {\beta _{34{\rm{b}}}} = {T_{{\rm{out}}}}, \\ {m_4}{{\ddot x}_4} + {K_{4x}}{x_4} - {K_{34}}\left( t \right){x_{34{\rm{n}}}}\cos {\alpha _{34{\rm{t}}}}\cos {\beta _{34{\rm{b}}}} + {C_{4x}}{{\dot x}_4} - \\{C_{34}}{{\dot x}_{34{\rm{n}}}}\cos {\alpha _{34{\rm{t}}}}\cos {\beta _{34{\rm{b}}}} = 0, \\ \begin{array}{*{20}{l}} {{m_4}{{\ddot y}_4} + {K_{4y}}{y_4} - {K_{34}}\left( t \right){x_{34{\rm{n}}}}\sin {\alpha _{34{\rm{n}}}} + {C_{4y}}{{\dot y}_4} - {C_{34}}{{\dot x}_{34{\rm{n}}}}\sin {\alpha _{34{\rm{n}}}} = 0, }\\ {{m_4}{{\ddot z}_4} + {K_{4z}}{z_4} - {K_{34}}\left( t \right){x_{34{\rm{n}}}}\sin {\beta _{34{\rm{n}}}} + {C_{4z}}{{\dot z}_4} - {C_{34}}{{\dot x}_{34{\rm{n}}}}\sin {\beta _{34{\rm{b}}}} = 0, } \end{array} \end{array} \right.$

(8)

式中：I_i(i=3, 4) 为斜齿轮3、4的转动惯量；m_i(i=3, 4) 为斜齿轮3、4的质量；K_ij、C_ij(i=3, 4；j=x, y, z) 分别齿轮支撑刚度、支撑阻尼；K₃₄(t)、C₃₄为齿轮副时变啮合刚度、啮合阻尼；T_out为输出扭矩。

1.2 传动系统动力学模型求解

在求解风电齿轮箱传动系统振动微分方程之前，首先在Romax软件中，设置齿轮及传动轴材料属性、轴承安装间隙、齿侧间隙以及齿轮修形参数，建立风电齿轮箱传动系统静力学模型并进行力学分析，得出各轴承支撑刚度、齿轮副时变啮合刚度曲线以及各齿轮副的静态传动误差曲线；再采用matlab软件将所得刚度曲线和静态传动误差曲线按傅里叶级数展开。各级齿轮副的时变啮合刚度、法向静态传递误差分别表示为

$\begin{array}{c} K\left( t \right) = {a_{k0}} + {a_{k1}}\cos \left( {{\omega _m}t} \right) + {b_{k1}}\sin \left( {{\omega _m}t} \right) + \cdots + {a_{ki}}\cos \left( {i{\omega _m}t} \right) + \\ {b_{ki}}\sin \left( {i{\omega _m}t} \right) + \cdots + {a_{kn}}\cos \left( {n{\omega _m}t} \right) + {b_{kn}}\sin \left( {n{\omega _m}t} \right), \end{array}$

(9)

$\begin{array}{c} {e_t} = {a_{e0}} + {a_{e1}}\cos \left( {{\omega _m}t} \right) + {b_{e1}}\sin \left( {{\omega _m}t} \right) + \cdots + {s_{ei}}\cos \left( {i{\omega _m}t} \right) + \\ {b_{ei}}\sin \left( {i{\omega _m}t} \right) + \cdots + {a_{em}}\cos \left( {n{\omega _m}t} \right) + {b_{en}}\sin \left( {n{\omega _m}t} \right), \end{array}$

(10)

式中：a_k0、a_k1、b_k1、…、a_ki、b_ki、…、a_kn、b_kn，a_e0、a_e1、b_e1、…、a_ei、b_ei、…、a_en、b_en为傅里叶级数展开项的系数；ω_m齿轮副啮合频率。

将时变啮合刚度、静态传递误差等所需参数代入式 (1)~式 (8)，采用4-5阶变步长Runge-Kutta法求解风电齿轮箱传动系统振动微分方程组，计算时间4 s，可得各传动构件的振动位移、振动速度、振动加速度。图 2给出了斜齿轮4纵向振动速度ẏ₄、振动加速度ÿ₄的频域曲线。

从频域图中可以看出，振动速度和振动加速度的最大值出现在660 Hz (第三级啮合频率) 附近，同时在其倍频处也出现峰值。斜齿轮4为第三级被动轮，其振动速度、振动加速度受第三级啮合频率影响较大。

将所求解出的齿轮转动角位移代入式 (11) 中，计算得出各齿轮副动态传动误差DTE (t)。

${\rm{DTE}}\left( t \right) = {r_{{\rm{bg}}}}\left[{{\theta _{\rm{g}}}\left( t \right)-\frac{{{z_{\rm{p}}}}}{{{z_{\rm{g}}}}}{\theta _{\rm{p}}}\left( t \right)} \right] = {r_{{\rm{bg}}}}{\theta _{\rm{g}}}\left( t \right) -{r_{{\rm{ bp}}}}{\theta _{\rm{p}}}\left( t \right), $

(11)

式中：θ_p、θ_g分别为主动轮、被动轮在各自坐标系下实际扭转微角度；r_bp、r_bg分别为主动轮、被动轮的基圆半径；z_p、z_g分别为主、被动轮齿数。

表 2给出各级传动误差时域分析结果，图 3给出斜齿轮3-4动态传动误差频域曲线。

从图 3可以看出：斜齿轮副3-4动态传动误差曲线在163 Hz (第二级啮合频率)、660 Hz及其倍频处出现峰值，而最大值出现在660 Hz。对于其他各级齿轮副，太阳轮-行星轮动态传动误差峰值主要表现在29 Hz (第一级啮合频率)、163 Hz及其倍频处，由于太阳轮轴与斜齿轮1同轴，故传动误差受第二级传动影响较大；行星轮-内齿圈动态传动误差峰值主要表现在29 Hz处，其次表现在660 Hz处；斜齿轮副1-2动态传动误差峰值主要表现在29、660 Hz及其倍频处，第二级齿轮传动作为中间传动齿轮副受第一级行星级和第三级斜齿轮传动影响较大。风电齿轮箱三级齿轮耦合传动，各级齿轮副动态传动误差表现出齿轮各级传动之间的耦合效应，传动误差均体现出每级齿轮副的啮合频率及其倍频。

2 风电齿轮箱辐射噪声预估 2.1 风电齿轮箱刚柔耦合模型建立及求解

风电齿轮箱辐射噪声模型的边界条件为传动系统对齿轮箱的激励力，传动系统通过轴承与齿轮箱连接，因此文中通过LMS Virtual.Lab对整个齿轮箱进行刚柔耦合计算，得出箱体处各轴承支反力的频域历程数据，将其作为计算声场的边界条件。

采用LMS Virtual.Lab软件分析减速器动力学特性时，主要步骤如下：

1) 在Ansys软件中对齿轮箱进行柔性化，并计算其约束模态；

2) 将柔性化模型和模态结果导入LMS Virtual.Lab软件中，结合其中所建立的齿轮箱传动系统，构建如图 4所示的风电齿轮箱刚柔耦合计算模型，9个轴承为箱体与转动轴的结合部零件；

3) 在轴承座圆心建立刚性耦合单元，以实现柔性箱体与传动系统之间的数据传递；

4) 根据上节给出的齿轮副动态传递误差及时变啮合刚度，在4对相互啮合的齿轮副间设置齿轮接触力 (gear contact)；

5) 在安装位置处施加固定约束，同时在输入轴处施加输入扭矩与转速，在输出端施加阻力矩；

6) 综合考虑轮齿时变啮合刚度、啮合阻尼、齿侧间隙、轴承支撑刚度与阻尼及由输入输出波动引发的外部激励，采用变步长向后差分法 (BDF) 对齿轮装置进行多体动力学仿真，求解总时间为4 s，得出箱体处各个轴承处支反力的时域和频域历程数据。

表 3给出了各轴承支反力计算结果。其中1~2号轴承连接箱体与行星架，行星架低速输入，轴承力并不大；太阳轮和第二级斜齿轮1同在3~4号轴承所在的轴上，带动两级齿轮传动，且太阳轮处于两轴承外的轴上，因此3号轴承受力最大；第二级斜齿轮2和第三级斜齿轮3同在5~6号轴承所在的轴上，由于斜齿轮受轴向力的原因，5号轴承放松，6号轴承压紧导致轴承力较大；7~9号轴承处于第三级斜齿轮4所在的轴上，同上7号处于预紧状态，8号轴承处于放松状态，9号轴承仅仅辅助轴支撑。

图 5给出了9号轴承支反力频域0~4 000 Hz数据图。9号轴承处于输出轴上，从频域图中可以看出，轴承受第三级输出齿轮啮合频率660 Hz影响较大，频率峰值集中在660 Hz及其倍频处，其频率成分体现出斜齿轮副3-4动态传动误差频率。其他轴承主要频率峰值集中在29、660 Hz及其倍频处，表现频率成分与所在的各级动态传动误差激励频率相近。

2.2 风电齿轮箱声振耦合模型建立及噪声计算

对于声振耦合计算方法，结构振动和声场分布是在一个耦合的环境中同时计算得出，在耦合边界处，结构法线方向的振动速度与流体的振动速度相同，此时声音对结构施加法线方向声压载荷，结构振动速度对声音施加速度输入，声场与结构场互相耦合。耦合环境中结构模型的动力学方程为

$\left( {{\boldsymbol{K}_{\rm{s}}} + j\omega {\boldsymbol{C}_{\rm{s}}}-{\omega ^2}{\boldsymbol{M}_{\rm{s}}}} \right) \cdot \left\{ {{\boldsymbol{u}_i}} \right\} + {\boldsymbol{K}_{\rm{c}}}\left\{ {{\boldsymbol{p}_i}} \right\} = \left\{ {{\boldsymbol{F}_{{\rm{s}}i}}} \right\}, $

(12)

式中：K_s、M_s、C_s分别为结构网格上没有受到约束部分的刚度矩阵、质量矩阵和阻尼矩阵，u_i为位移，K_c为耦合刚度矩阵，p_i为声压，F_si为结构的约束部分、已知载荷和外部声压的激励贡献量。

声学方程为

$\left( {{\boldsymbol{K}_{\rm{a}}} + {\rm{j}}\omega {\boldsymbol{C}_{\rm{a}}} + {\omega ^2}{\boldsymbol{M}_{\rm{a}}}} \right) \cdot \left\{ {{\boldsymbol{p}_i}} \right\}-{\omega ^2}{\boldsymbol{M}_{\rm{c}}}\left\{ {{\boldsymbol{u}_i}} \right\} = \left\{ {{\boldsymbol{F}_{{\rm{a}}i}}} \right\}, $

(13)

式中：K_a、M_a、C_a分别为声学的刚度矩阵、质量矩阵和阻尼矩阵，M_c为耦合质量矩阵，F_ai为已知声压和振动速度边界的激励贡献量^[16]。

将风电齿轮箱箱体结构网格、声学网格和场点网格导入软件的声学模块Acoutics，建立声振耦合计算模型，如图 6所示。

将计算求得的轴承支反力频域历程数据作为声振耦合的边界条件，在9个轴承孔处的刚性耦合点添加频域载荷。采用声振耦合有限元法计算风电齿轮箱辐射噪声，结构网格与声学网格间的数据交换需通过插值计算，声学网格则需要完全包络结构网格，并在交界面建立包络单元，因此在Ansys中建立一个声学球体，球体内部空腔表面完全与齿轮箱外表面相对应，其网格划分时保证箱体与空气交界面处一个波长内包含6个单元；同时在距箱体外表面1 m处建立场点网格。

在声学模块Acoutics设置外界空气属性，空气密度为1.255 kg/m³，空气传播速度为340 m/s，参考声压为2×10^-5 Pa，再设置求解参数，设定计算频段、求解精度，最后采用声振耦合有限元法进行求解齿轮箱表面振动以及场点辐射噪声。

图 7给出部分频段的齿轮箱外部场点辐射噪声声压云图。场点声压最大值出现在500 Hz频率处，此时表面声压为84.2 dB，出现在1和4号轴承座附近，位于输入行星架处的轴承和1号斜齿轮所在的中间轴轴承处。63、125 Hz频段场点声压最大值出现在齿轮箱上端部，250、1 000 Hz频段场点声压最大值出现在齿轮箱左右两侧，2 000 Hz之后各处声压值普遍开始下降。

3 振动噪声试验测试与仿真对比

搭建试验台架并安装风电齿轮箱，进行了振动噪声测试，测试环境温度为8 ℃，湿度66%。测试运行工况为输入转速16.756 r/min，功率1 660 kW。振动噪声测试仪器见表 4。

测试现场及台架如图 8所示，左侧为被测齿轮箱，右侧为陪试减速箱。

3.1 风电齿轮箱振动测点测试

根据测试需求，在齿轮箱上各布置4个测点，测点编号为1~4，其中测点1~2为轴承座位置，测点3和4为支撑位置，分别对4个测点的竖直方向进行了振动测试。图 9为振动及噪声测点布置图。

在声振耦合模型中提取对应测点的节点振动加速度频率历程数据，计算其有效值，将其与试验测点振动加速分析结果对比，如表 5所示，图 10给出了测点1振动加速度试验与仿真频谱值对比图。

对比分析图 10所示的测点1振动加速度频域曲线，二者皆在29、163、326、489、660、1 320 Hz处出现峰值，并且在326 Hz达到最大值，此处轴承座处于第二级主动轮所在的轴上，表现出受第二级齿轮传动影响较大，其余各测点试验与仿真的加速度有效值以及频率分布也比较接近。

3.2 风电齿轮箱辐射噪声测试

在距离齿轮箱1 m远的位置各布置4个测点，测点编号为Ⅰ~Ⅳ，分别对4个测点的空气噪声进行测试。并在声振耦合模型中提取对应测点场点辐射噪声频率历程数据，对数据进行倍频程处理，与试验测点空气噪声结果相对比，如图 11所示。

由图 11可知，测点Ⅰ~Ⅳ的A计权声压级辐射噪声在500 Hz频段处达到最大值，其值分别为81.48、84.72、83.94及81.94 dB。对比分析仿真与测试结果，两者在250~2 000 Hz频段辐射噪声吻合良好，但在63、125、4 000 Hz，仿真比测试结果低。

4 结论

1) 建立了包含时变啮合刚度、啮合阻尼、传递误差等多种因素的多级齿轮传动31自由度弯-扭-轴耦合集中参数动力学模型。

2) 在LMS Virtual.Lab中建立刚柔耦合运动学模型，考虑各级齿轮副动态传动误差对齿轮箱振动的影响；得出齿轮箱轴承座支反力，支反力频率成分与动态传动误差频率谱相对应。

3) 采用声学有限元法，对齿轮箱进行声振耦合振动噪声计算，将轴承支反力作为传动系统对齿轮箱的激励力，得出齿轮箱表面振动和场点辐射噪声。整体结果表明声振耦合仿真值与试验值吻合良好。