工程中对转子进行动力学分析时，采用梁单元对转子进行建模具有精度较高、计算速度快等优势。然而传统的梁单元无法考虑转子的实际三维形状对弯曲刚度的影响，因此在用梁单元表达带有截面突变的阶梯转子时，往往带来一定的误差^[1]。误差大小与材料性质(弹性模量、Poisson比等)、尺寸以及阶梯处的连接方式等因素有关。弯曲刚度的表达误差是构成转子动力学分析中转子临界转速、振型和不平衡响应等关键参数计算误差的主要来源之一。

为消除上述误差，可以采取2种方式：一种方式是采用三维单元进行转子的建模，但会导致计算量的大幅增加，而且不利于与其他模型联立开展复杂的动力学分析；另一种方式是将阶梯处使用一个等效的梁单元进行表达，使得该段梁单元的弯曲刚度与实际转子的弯曲刚度相近，从而减少误差。

在建立等效梁单元模型时，目前常用的方法是按照45°原则对阶梯处进行等效^[2-4]，即采用45°锥形作为弯曲刚度计算的有效尺寸(以下简称“45°法”)。Nation^[5]研究了45°法在阶梯转子扭转刚度等效中的应用，通过ANSYS仿真验证了45°以外区域对扭转刚度的影响很小。Jweeg等^[6]研究了阶梯处的损伤对轴的动态性能的影响。Li等^[7]提出了一种基于扭转刚度的阶梯转子等效方法并应用在燃气轮机性能计算中，通过ANSYS计算转角变化，从而推算出阶梯段的等效扭转刚度和惯性矩。于军等^[8]提出将阶梯处用一段锥形代替，并将锥形进一步简化为直径等于两端平均直径的圆柱等效段，等效段长度通过轴的最大挠度确定，该方法得到了一定的应用^[9]。张贤彪等^[10]也采用了类似的方法，计算了等效段长度的变化规律。国内外学者针对连续变截面梁的研究已经相对成熟^[11-12]。拉杆转子中的截面变化等效问题也是研究的热点之一，拉杆转子阶梯处通常采用一定的接触刚度对转子进行建模^[13-14]。相比传统梁单元建模方法，以上建模方法虽然能够在一定程度上减少误差，但均未直接针对阶梯转子的弯曲刚度进行分析，因此仍不可避免地存在等效误差。

为进一步提高阶梯转子的建模精度，本文提出一种基于弯曲刚度的等效梁单元建模方法。通过ANSYS三维实体单元计算得到阶梯转子的弯曲刚度，作为梁单元建模的参考值，对比不同梁单元等效建模方法的特点，得到基于弯曲刚度的等效梁单元建模方法，从而为提高阶梯转子梁单元建模精度提供新的方法和依据。转子中的阶梯可以分为轴阶梯和圆盘阶梯，在不考虑普通阶梯轴的截面突变影响时，将轴直接简化为梁单元将造成计算中采用的弯曲刚度高于实际值；而忽略圆盘对轴的弯曲刚度的贡献将造成计算中采用的弯曲刚度低于实际值。本文将分别针对轴阶梯和圆盘阶梯2种阶梯形式进行分析，并将所提出的建模方法应用于转子动力学计算中，以验证方法的有效性。

1 阶梯轴的等效建模方法

将阶梯处进行等效建模的目的是使等效后的梁单元在阶梯处的弯曲刚度与实际一致。

1.1 基于弯曲刚度的等效建模方法

对于梁单元模型，弯曲刚度EI由轴的弹性模量E和截面惯性矩I乘积得到。考虑直径为d的圆截面，EI表示为

$ E I=\frac{\pi E d^4}{64}. $

(1)

对于实体单元模型，截面的等效弯曲刚度(EI)_solid根据弯矩M和中轴线挠度w进行计算：

$ (E I)_{\text {solid }}=\frac{M}{w^{\prime \prime}} . $

(2)

采用ANSYS三维实体单元模型计算了阶梯轴在施加静载荷后的形变，单元类型为SOLID186。轴弹性模量E=2.11×10¹¹ Pa，Poisson比ν=0.3；施加静载荷F=1 000 N，位于轴末端，方向垂直于中轴线；不考虑重力影响。图 1a中，建立坐标系O-xyz，取阶梯轴左端截面圆心为原点；阶梯轴左半部分的长度l₁=100 mm，直径d₁=40 mm；右半部分的长度l₂=100 mm，直径d₂=32 mm；网格基本尺寸为1 mm，网格单元个数为399 152。计算结果如图 1b所示，通过建立路径，获取中轴线上各点的挠度，由式(2)计算得到各截面的(EI)_solid。

以ANSYS计算得到的(EI)_solid为参考值，图 2对比了基于不同梁单元模型计算得到的弯曲刚度。在传统梁单元模型中，梁单元弯曲刚度随着轴径变化而突变，阶梯处的弯曲刚度明显大于参考值。采用45°法等效时，等效梁单元几何形状为锥形，弯曲刚度变化趋势更接近参考值，然而阶梯处的弯曲刚度与参考值仍存在一定差异。

由图 2可知，现有梁单元模型计算得到的弯曲刚度均与参考值存在一定差异，锥形梁单元的弯曲刚度变化趋势接近参考值，因此本文采用锥形梁单元进行弯曲刚度的等效。不同于45°法的固定锥角，本方法采用的锥形梁单元的长度l(或锥角α)由以下原则确定：将截面突变处作为锥形梁单元终点，锥形梁单元起点通过最小二乘法计算得到，满足等效区间内锥形梁单元的的等效弯曲刚度(EI)_eq与参考值数值偏差最小(以下简称“等效刚度法”)，其数学表达为

$ \min S(l)=\sum\limits_{i=1}^{n}{{{\left( {{(EI)}_{\text{eq}}}(l, i)-{{(EI)}_{\text{solid}}}(i) \right)}^{2}}}. $

(3)

其中i为等效区间内各计算截面的序号。

由于(EI)_eq是l的函数，因此取等效区间内(EI)_eq和(EI)_solid差值的平方和S(l)最小，即可得到最优等效梁单元长度l_e。

1.2 等效计算结果

对于图 1中的阶梯转子，采用等效刚度法计算得到l_e=9.8 mm，由此得到等效锥形梁单元各截面的(EI)_eq。图 3为(EI)_eq与参考值的对比，可以发现由等效刚度法计算得到的阶梯处弯曲刚度与参考值非常一致，很好地反映了阶梯处弯曲刚度的变化。

图 3中，除了锥形梁单元等效区间以外，在等效区间外靠近d₁的一侧，存在一部分弯曲刚度加强区，即参考弯曲刚度大于按照梁单元计算得到的弯曲刚度；在等效区间外靠近d₂的一侧，存在一部分弯曲刚度削弱区，即参考弯曲刚度小于按照梁单元计算得到的弯曲刚度。考虑到上述2部分区域在弯曲刚度计算上的抵消作用，且等效处理困难，因此在等效时予以忽略。

在工程计算中，如果按照上述方法对每一个算例都进行一次求解，则过程过于繁琐，不利于工程应用。因此本文给出了l_e的变化规律，通过查询图表或代入公式直接得到l_e。由于d₁、d₂等尺寸参数均可能对l_e产生影响，因此计算了不同尺寸下的l_e及对应的锥角α，并对尺寸参数进行了归一化处理，结果如图 4所示。

由图 4可知，随着阶梯相对直径d₂/d₁的增加，等效梁单元相对长度l_e/d₁及对应α呈近似二次多项式变化，拟合公式为：

$ \frac{l_{\mathrm{e}}}{d_1}=-0.7159\left(\frac{d_2}{d_1}\right)^2+0.7566 \frac{d_2}{d_1}+0.0943, $

(4)

$ \alpha=-58.063^{\circ}\left(\frac{d_2}{d_1}\right)^2+16.429^{\circ} \frac{d_2}{d_1}+46.083^{\circ} . $

(5)

随着d₂/d₁在0.5~0.95变化，α变化范围约为40°~10°，可见，采用固定的45°锥角对弯曲刚度进行等效存在不可避免的偏差，尤其是在d₂/d₁较大时。

2 整体式圆盘的等效建模方法

转子中的圆盘可以视为一种特殊的阶梯轴，与普通阶梯轴不同，圆盘两端有2处相邻的截面突变，可能存在相互干扰，需要重新考虑其等效规律。

2.1 基于弯曲刚度的等效建模方法

转动机械通常存在叶轮、推力盘等结构，在转子动力学分析中，一般将其等效为刚性圆盘，如果不考虑圆盘对轴的刚度的加强作用，将造成计算误差。王正^[4]给出了圆盘组对轴的弯曲刚度的大致影响规律。Riley^[15]通过理论分析与试验相结合的方法，得到了不同尺寸的圆盘对扭转刚度的影响规律。圆盘对轴的弯曲刚度的加强作用仍需进一步研究。

圆盘是叶轮等结构的简化形式，其在轴上的安装方式包括整体锻造和热套装配等，热套装配的过盈量对弯曲刚度的加强作用具有显著的影响，影响规律较为复杂，且难以保证计算精度。本文仅研究通过整体锻造等方式加工得到的整体式圆盘的等效建模方法。

整体式圆盘的等效梁单元建模方法与普通阶梯转子的等效梁单元建模方法相似，通过ANSYS三维实体单元计算带圆盘转子在施加载荷时的挠度，由式(2)计算得到对应的弯曲刚度(EI)_solid，作为参考值。根据(EI)_solid变化规律，在圆盘阶梯的截面突变处用对称的2段锥形进行等效，等效锥形的最大截面直径d_e(或锥角α)同样通过最小二乘法确定。

2.2 等效计算结果

转子轴径d=40 mm，圆盘直径D=60 mm，厚度B=20 mm，其余参数与图 1中算例相同。采用等效刚度法求得d_e=49.5 mm，对应的等效弯曲刚度(EI)_eq与参考值的对比如图 5所示。

从图 5可以看出，圆盘对轴的弯曲刚度具有明显的增强作用，提出的等效方法能够很好地反映出轴在圆盘处弯曲刚度的变化趋势；当采用45°法进行等效时，圆盘处的弯曲刚度明显大于参考值。

为了便于工程计算，给出了不同尺寸下圆盘对轴的弯曲刚度的加强作用规律，尺寸参数经过归一化处理，结果如图 6所示。B/d为圆盘的相对厚度，D/d为圆盘的相对直径，d_e/d为等效锥形最大截面的相对直径。对于B/d一定的圆盘，随着D/d的增加，d_e/d随之增加。当D/d增加到一定阈值D_th/d时，d_e/d将不再增加，极值为d_emax/d，对应α达到极值α_max。对于不同的B/d，D_th/d、d_emax/d和α_max均不相同，其变化规律如图 7所示。

由图 7可知，D_th/d、d_emax/d随B/d分别呈近似二次多项式变化，α_max随B/d呈近似线性变化，变化规律分别表示为：

$ \frac{D_{\mathrm{th}}}{d}=0.7315\left(\frac{B}{d}\right)^2+0.3095 \frac{B}{d}+1, $

(6)

$ \frac{d_{\mathrm{emax}}}{d}=1.1764\left(\frac{B}{d}\right)^2-0.0959 \frac{B}{d}+1, $

(7)

$ \alpha_{\max }=48.609^{\circ} \frac{B}{d}+1.5382^{\circ} . $

(8)

随着B/d在0.2~0.8变化，α_max变化范围约为10°~40°，因此45°法在进行圆盘弯曲刚度等效时存在显著误差，尤其是对较薄的圆盘进行等效时。

图 6中，通常情况下，d_e < D，对应的尺寸关系如图 8a所示；然而当D/d较小时，可能存在d_e>D，对应的尺寸关系如图 8c所示，原因是当D/d接近于1时，圆盘两端的截面可以看作2个独立的轴阶梯，当采用对称锥形进行弯曲刚度等效时，假如d_e≤D，圆盘中间部分的弯曲刚度将被过低地表达，因此需要增加d_e/d才能获得更接近实际的弯曲刚度，从而导致d_e>D；此外，在特定情况下，存在d_e=D，如图 8b所示。

3 阶梯转子的转子动力学计算

下面通过算例验证本文提出的建模方法的计算精度及相对传统梁单元及45°法的优势。

3.1 计算方法

在转子小幅线性运动情况下，转子的运动方程可以描述为

$ \boldsymbol{M} \ddot{\boldsymbol{u}}+(\boldsymbol{C}+\boldsymbol{G}) \dot{\boldsymbol{u}}+\boldsymbol{K} \boldsymbol{u}=\boldsymbol{Q}. $

(9)

式中，$\boldsymbol{M}=\left[\begin{array}{cc} \boldsymbol{M}_1 & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{M}_1 \end{array}\right], \quad \boldsymbol{G}=\left[\begin{array}{cc} \boldsymbol{0} & \boldsymbol{G}_1 \\ -\boldsymbol{G}_1 & \boldsymbol{0} \end{array}\right], \quad \boldsymbol{K}=\left[\begin{array}{cc} \boldsymbol{K}_1 & \bf{0} \\ \bf{0} & \boldsymbol{K}_1 \end{array}\right], \boldsymbol{u}=\left[\begin{array}{l} \boldsymbol{u}_1 \\ \boldsymbol{u}_2 \end{array}\right], Q=\left[\begin{array}{l} \boldsymbol{Q}_1 \\ \boldsymbol{Q}_2 \end{array}\right]$。

M₁、G₁、K₁分别为单个平面上转子的质量矩阵、陀螺矩阵和刚度矩阵，u₁、u₂和Q₁、Q₂分别为2个平面上的广义坐标和广义力，在进行临界转速求解时，通常忽略阻尼矩阵C。令z=u₁+iu₂，则式(9)的齐次项表述为

$ \boldsymbol{M}_1 \ddot{\boldsymbol{z}}-\mathrm{i} \boldsymbol{G}_1 \dot{\boldsymbol{z}}+\boldsymbol{K}_1 \boldsymbol{z}=\bf{0}. $

(10)

取z=z₀e^iωt，得到频率方程：

$ \left|-\boldsymbol{M}_1 \omega^2+\boldsymbol{G}_1 \omega+\boldsymbol{K}_1\right|=\bf{0} . $

(11)

对频率方程进行求解，得到转子的临界转速。式(11)中，K₁由各个梁单元的刚度矩阵k_e及轴承等支撑的刚度矩阵组装而成。对于长度为l的等截面圆柱梁单元，EI为常数，有

$ \begin{gathered} \boldsymbol{k}_{\mathrm{e}}=\int_0^l E I \boldsymbol{N}^{\prime \prime \mathrm{T}} \boldsymbol{N}^{\prime \prime} \mathrm{d} z= \\ \frac{2 E I}{l^3}\left[\begin{array}{cccc} 6 & 3 l & -6 & 3 l \\ 3 l & 2 l^2 & -3 l & l^2 \\ -6 & -3 l & 6 & -3 l \\ 3 l & l^2 & -3 l & 2 l^2 \end{array}\right] . \end{gathered} $

(12)

式中N为形函数向量，其表达式可见文[16]的式(8.6.39)。

对于两端直径分别为d₁、d₂，长度为l的锥形梁单元，EI为随z变化的变量，有

$ \begin{gathered} \boldsymbol{k}_{\mathrm{e}}=\int_0^l(E I)(z) \boldsymbol{N}^{\prime \prime \mathrm{T}} \boldsymbol{N}^{\prime \prime} \mathrm{d} z= \\ \frac{\pi E}{1120 l^3}\left[\begin{array}{llll} k_{11} & k_{12} & k_{13} & k_{14} \\ k_{21} & k_{22} & k_{23} & k_{24} \\ k_{31} & k_{32} & k_{33} & k_{34} \\ k_{41} & k_{42} & k_{43} & k_{44} \end{array}\right] . \end{gathered} $

(13)

其中：

$ \begin{gathered} k_{11}=k_{33}=-k_{13}=-k_{31}= \\ 66 d_1^4+30 d_1^3 d_2+18 d_1^2 d_2^2+30 d_1 d_2^3+66 d_2^4, \\ k_{12}=k_{21}=-k_{23}=-k_{32}= \\ \left(47 d_1^4+22 d_1^3 d_2+9 d_1^2 d_2^2+8 d_1 d_2^3+19 d_2^4\right) l \\ k_{14}=k_{41}=-k_{34}=-k_{43}= \\ \left(19 d_1^4+8 d_1^3 d_2+9 d_1^2 d_2^2+22 d_1 d_2^3+47 d_2^4\right) l, \\ k_{22}=\left(34 d_1^4+18 d_1^3 d_2+8 d_1^2 d_2^2+4 d_1 d_2^3+6 d_2^4\right) l^2, \\ k_{24}=k_{42}=\left(13 d_1^4+4 d_1^3 d_2+d_1^2 d_2^2+4 d_1 d_2^3+13 d_2^4\right) l^2, \\ k_{44}=\left(6 d_1^4+4 d_1^3 d_2+8 d_1^2 d_2^2+18 d_1 d_2^3+34 d_2^4\right) l^2 . \end{gathered} $

采用等效刚度法对带有阶梯轴和圆盘的转子进行等效梁单元建模和转子动力学计算的流程如下：

1) 对于阶梯轴，根据阶梯尺寸，查询图 4或代入式(4)—(5)，得到对应的l_e或α。

2) 对于带有整体式圆盘的转子，根据圆盘尺寸，查询图 7或代入式(6)得到D_th/d。如果D/d≥D_th/d，直接通过查询图 7或代入式(7)—(8)得到d_emax/d或α_max；如果D/d < D_th/d，则通过查询图 6得到相应的d_e/d。

3) 根据前2步得到的锥形尺寸，通过式(13)获取k_e，代入K进行转子动力学计算。在对梁单元进行等效建模时，仅在生成K时采用等效的直径，在生成M和G时应采用原几何尺寸。

3.2 计算结果

按照上述流程，对图 9a的阶梯轴进行了等效建模，转子左右对称，阶梯轴从左至右前3处阶梯位置l_e分别是4.24、5.3、6.125 mm；按照45°法得到l_e分别为2、2、2.5 mm。同样按照上述流程对图 9c的带圆盘转子进行了等效建模，得到d_e=24.8 mm；按照45°法得到d_e=30 mm。以上转子两端均由轴承支撑，轴承宽度为10 mm，x和y方向的刚度系数为K_xx=K_yy=1×10¹⁰ N/m，忽略交叉刚度和轴承阻尼。按照式(13)求得等效梁单元刚度矩阵，代入式(11)进行转子动力学分析，得到了基于梁单元模型的临界转速。

图 9b、9d中，在ANSYS中采用三维实体单元对以上阶梯转子进行建模，通过模态分析求解得到转子的临界转速，作为评价等效梁单元模型计算精度的参考依据。主要计算单元为SOLID186，在网格划分时，对截面突变处的网格数量进行了加密，以提高计算精度。

通过不同梁单元建模方法计算得到的转子前3阶正进动临界转速与通过ANSYS计算得到的参考值的对比如表 1所示。在阶梯轴计算中，采用传统梁单元模型时，由于梁单元弯曲刚度偏大，导致计算得到的第1~3阶临界转速偏高4.90%~10.59%，相对误差较大；采用45°法进行梁单元建模，相对误差减小至3.46%~8.88%；采用等效刚度法进行梁单元建模，相对误差减小至1.35%~6.58%，其中对第1~2阶临界转速的求解精度明显提高，相对误差小于3%。

在带圆盘转子计算中，采用传统梁单元模型时，由于没有对圆盘处轴的弯曲刚度进行强化，导致计算得到的临界转速偏低，相对误差为-3.59%~-9.14%；采用45°法进行梁单元建模，由于等效尺寸过大，导致计算结果偏高2.21%~4.32%；采用等效刚度法进行梁单元建模，结果与参考值较为接近，相对误差仅为-2.23%~1.23%，其中对第1~2阶临界转速的求解精度明显提高，第1阶临界转速相对误差仅为-0.13%。

需要说明的是，采用梁单元进行转子动力学分析时的误差是由多种因素共同造成的，阶梯处的弯曲刚度只是造成计算误差的原因之一。只有当阶梯处弯曲刚度在误差中占主导作用时，采用本文提出的等效方法才能够显著降低计算误差。在部分算例中，转子同时存在阶梯轴和圆盘，在阶梯轴对弯曲刚度的削弱效应和圆盘对弯曲刚度的增强效应的共同作用下，误差部分抵消，仅对其中一项误差进行修正可能反而会增加总的计算误差。例如，在图 9a的阶梯轴的中心位置增加1个B=10 mm、D=80 mm的圆盘，ANSYS三维实体单元模型计算得到的第1阶临界转速参考值为4 333 rad/s；采用传统梁单元模型的计算结果为4 488 rad/s；如果仅针对圆盘进行等效而不对阶梯轴进行等效，计算结果为4 521 rad/s，误差增加；对圆盘和阶梯轴均进行等效的计算结果为4 367 rad/s，与参考值最为接近。此外，针对空心阶梯轴、圆盘两端直径不同等特殊情况，需要进一步研究其等效规律，而不能简单采用本文中提出的建模方法。

4 结论

通过三维实体单元对阶梯转子进行了静力学分析，在此基础上提出了基于弯曲刚度的等效梁单元建模方法，主要结论如下：

1) 与ANSYS静力学仿真结果相比，在不对阶梯轴的阶梯处进行弯曲刚度等效时，传统梁单元计算得到的弯曲刚度明显偏大；采用45°法的锥形等效梁单元可以反映弯曲刚度变化趋势，然而与参考值存在偏差。

2) 锥形梁单元符合阶梯处弯曲刚度的变化趋势，采用最小二乘法拟合得到了与参考值最为接近的等效梁单元尺寸；随着d₂/d₁在0.5~0.95增加，α呈近似二次多项式递减，范围约为40°~10°。

3) 对称的锥形梁单元符合圆盘处弯曲刚度的变化趋势，随着圆盘D/d的增加，d_e/d增加；当D/d增加到一定阈值时，d_e/d及α达到极值；D_th/d、d_emax/d及α_max与B/d相关，随着B/d在0.2~0.8增加，α_max呈近似线性递增，范围约为10°~40°。

采用本文提出的等效刚度法对转子进行了等效梁单元建模，并在此基础上进行了转子动力学分析。结果表明，传统梁单元模型的相对误差较大；采用45°法等效梁单元能够减小阶梯轴计算相对误差，计算带圆盘转子得到的临界转速偏大；采用等效刚度法进行梁单元建模得到的转子临界转速与参考值更为接近，尤其是第1~2阶临界转速的相对误差小于3%；相比于45°法，采用本文提出的建模方法能够显著提高转子动力学计算中临界转速的求解精度。

等效建模时应综合考虑阶梯轴和圆盘的作用，以取得较高的精度。