简介面向航空发动机转子动力学分析与设计验证场景这份资源提供基于传递矩阵法的临界转速计算工具核心是一个轻量MATLAB脚本。脚本通过输入轴段长度、截面参数、材料属性与支承刚度等模型信息即可建立传递矩阵并求解转子固有频率进而得到各阶临界转速帮助工程师在概念设计阶段快速评估共振风险。传递矩阵法将复杂连续轴系离散为若干单元编码简洁、计算开销小脚本完整展示了从模型离散、矩阵构建到频率求解的代码实现对学习转子动力学和MATLAB数值方法均有参考价值。压缩包共1个文件文件类型为m脚本整个资源包仅1KB轻巧易用。目前已有237人学习适合航空发动机转子、透平机械等旋转机械的临界转速快速估算也可作为课程设计或进一步开发轴承-转子系统分析程序的基础模板无论是初学者理解传递矩阵流程还是工程师快速获得基准结果都能提供直观支持。1. 传递矩阵法算转子临界转速先搞懂这一个矩阵乘法航空发动机转子的临界转速计算听起来是有限元的事但工程一线最常见的资料包就是那类rotor_criticalspeed_calculation.rar解开后一个 Fortran 文件或 MATLAB 脚本核心算法是传递矩阵法。原因是设计初期要反复改盘的位置、轴承刚度、轴径有限元建模一次要 20 分钟传递矩阵法把转子切成几十段每段一个 4×4 矩阵乘法连乘一次只要几毫秒扫一遍频率就能找到各阶临界转速。本文不依赖某个商业软件直接讲清状态向量、场矩阵、点矩阵怎么拼给一套 Python 实现参数设置和漏根排错也会展开。适合刚接转子动力学任务的强度工程师以及想把手算例题程序化的学生。2. 传递矩阵法建模状态向量与场/点矩阵2.1 转子截面的 4 个状态分量怎么排传递矩阵法把一个转子从小端到大端切成若干段每一段的连接截面用一个状态向量描述。对横向弯曲状态向量取[ Z [y,\ \theta,\ M,\ T]^T ]其中 y 是挠度θ 是转角M 是弯矩T 是剪力。整个转子的振动问题变成从一个端面出发逐段把状态向量乘到另一端最后满足两端边界条件的频率就是临界转速。这个思路和弹塑性梁的初参数法同源只是把初参数法的双曲函数解换成数值矩阵。正负号约定是代码对不上的第一个原因。弯矩和剪力的符号必须和材料力学里梁的内力符号约定一致弯矩逆时针为正剪力使左段向上右段向下为正。如果从某份源码里抄的矩阵把 T 的符号反了算出来的临界转速频率规律是对的但振型会上下镜像而且和支承耦合时误差很隐蔽。拿到任何传递矩阵代码第一步永远是把符号约定写进注释再拿一个单盘悬臂转子的解析解做标定。切段规则按转动节点来。两端的轴段、盘、轴承都落在离散节点上盘和轴承只占据一个节点的位置不占用一段轴长。对航空发动机转子压气机盘和涡轮盘比同段轴至少重一个量级这样集中质量的处理足够准确叶片离心刚化效应本来也不在临界转速计算的基本公式里后续修正时再单独加项。2.2 场传递矩阵无质量轴段的 4×4 传递关系节点 i 到节点 i1 之间有一段等截面圆轴长度 l抗弯刚度 EI。先忽略这段轴本身的分布质量它只提供弹性变形。从截面 i 到截面 i1状态向量满足[ Z_{i1} F_i , Z_i ][ F_i \begin{bmatrix} 1 l \dfrac{l^2}{2EI} \dfrac{l^3}{6EI}\ 0 1 \dfrac{l}{EI} \dfrac{l^2}{2EI}\ 0 0 1 l\ 0 0 0 1 \end{bmatrix} ]这个矩阵的每一列都对应一个初始状态分量对末端各状态量的影响。例如第一列表示当截面 i 处 y1、其他状态量为 0 时经过这段悬臂式变形末端产生挠度、转角、弯矩和剪力分别是矩阵第一列的四个元素。第二列以上依次类推。矩阵是低三角的说明弯矩只由剪力和外载荷产生挠度则由上一截面的所有状态累积而成符合简单的梁变形几何关系。注意这个 F 矩阵里没有频率 ω也没有质量项。这是把轴段质量“砍掉”后的简化场矩阵。相比直接采用带分布质量的精确解里面会出现 cos/sinh 组合分段集中质量的做法在段数足够时误差只推向高频段而实际关心的低阶临界转速正好在低频段。因此工程代码里大量使用这种无质量场矩阵加集中质量点矩阵的组合计算速度快物理意义清楚。如果你要算第三阶以上建议把每段再加密一倍而不是换更复杂的场矩阵。2.3 点传递矩阵把盘和轴承压进一个 4×4 矩阵轴段经过节点时节点上的集中质量 m 和轴承刚度 k 会对状态向量产生一个突变。设节点左侧状态为 (Z_l)右侧为 (Z_r)挠度、转角、弯矩在节点两侧连续只有剪力因惯性力和支承反力发生跳跃。假设转子以频率 ω 做同步正进动惯性力幅值为 (m\omega^2 y)方向与挠度一致轴承反力为 (-k y)。于是[ T_r T_l m\omega^2 y - k y ]写成矩阵[ P \begin{bmatrix} 1 0 0 0\ 0 1 0 0\ 0 0 1 0\ m\omega^2 - k 0 0 1 \end{bmatrix} ]点在矩阵右下角是 1表示剪力从节点左侧平移到右侧左下角那个 (m\omega^2-k) 是响应项。没有盘的节点直接让 m0没有支承的节点让 k0。把某一小段轴和其后置的节点合在一起就是一级基本传递[ T_i P_i F_i ]从最左端开始连乘到最右端[ T_{total} P_n F_n \cdots P_1 F_1 ]如果左端状态为 (Z_0)右端就是 (Z_n T_{total} Z_0)。左端边界条件会消掉两个未知数右端边界条件又给出两个方程最后只留下一个关于 ω 的残余量。残余量为零时ω 就是临界转速。这一整套过程里没有求解大型线性方程组矩阵永远只有 4×4 连乘这就是它快的原因。3. 实现用 Python 写一个能跑的转子临界转速计算程序3.1 把转子的几何和材料拆成 segments/masses 两个表开始写代码前先把模型数据整理成两张表。第一张表是轴段列表每段含长度 l 和抗弯刚度 EI第二张表是节点属性每个节点含集中质量 m 和支承刚度 k。以下代码以一个单盘悬臂转子和一个双支承双盘转子为例。import numpy as np def field_matrix(l, EI): 无质量轴段的场传递矩阵 return np.array([ [1.0, l, l * l / (2.0 * EI), l ** 3 / (6.0 * EI)], [0.0, 1.0, l / EI, l * l / (2.0 * EI)], [0.0, 0.0, 1.0, l], [0.0, 0.0, 0.0, 1.0] ]) def point_matrix(m, w2, k0.0): 节点上的点传递矩阵w2 为 omega^2 P np.eye(4) P[3, 0] m * w2 - k return Pfield_matrix对应 2.2 节的 F 矩阵point_matrix对应 P 矩阵。注意point_matrix里参数名是w2调用时要传 (\omega^2)不要传成 (\omega)否则量纲全错。轴段的数据结构用列表里的元组表示每个元组是(l, EI, m, k)顺序从左端到右端。m 和 k 是该段右侧节点的质量与轴承刚度。这种数据排法最直观也方便从 Excel 或旧 Fortran 数据文件转换。3.2 扫频找残余行列式过零点的完整代码以一个左端自由、右端自由的转子为例。左端没有弯矩和剪力所以 (M_00, T_00)未知的只有 (y_0) 和 (\theta_0)。右端同样自由要求 (M_n0, T_n0)。由 (Z_n T_{total} Z_0)右端弯矩、剪力分别是总矩阵第三、四行的前两列与 (y_0, \theta_0) 的组合def total_matrix(rotor, w): 连乘所有段的场矩阵和点矩阵 T np.eye(4) for l, EI, m, k in rotor: T point_matrix(m, w * w, k) field_matrix(l, EI) T # 先过场矩阵再过点矩阵顺序和节点布置一致 return T def residual(rotor, w): 自由-自由边界条件下的残余行列式 T total_matrix(rotor, w) # 右端弯矩、剪力仅依赖于左端 y0, theta0 # 非零解条件为 2x2 行列式等于 0 return T[2, 0] * T[3, 1] - T[2, 1] * T[3, 0]残余量随 ω 变化每过一次零点就对应一个临界转速。扫频函数从 0 开始小步长推进符号变化后原地二分细化把根锁定到 10^{-6} rad/s 精度def find_critical_speeds(rotor, w_max1000.0, dw1.0): roots [] w 0.0 r0 residual(rotor, w) while w w_max: r1 residual(rotor, w dw) if r0 0.0: roots.append(w) elif r0 * r1 0: # 跨越零点的粗判断 a, b w, w dw for _ in range(60): mid 0.5 * (a b) rm residual(rotor, mid) if rm * r0 0: b mid else: a mid r0 rm roots.append(0.5 * (a b)) w dw r0 residual(rotor, w) return rootsr0 * r1 0是常见的异号判断但如果残余量在某个频率附近有上下尖峰而不真正过零这种判断会丢根或误报。更稳的做法是先画一条残余量曲线看形状再决定步长。代码里的dw1.0意味着每 1 rad/s 采样一次对通常几百 rad/s 的低阶临界转速够密对频率间隔小于 1 rad/s 的密集模态则不够后面第 4 章会讲步长怎么选。3.3 临界转速落点参数表与计算精度给一个双支承双盘转子的算例模拟简化后的发动机高压转子总长 1.2 m两段轴各 0.6 m轴径 50 mm材料弹性模量按 GH4169 取 210 GPa抗弯刚度约为 6.45×10^4 N·m²。两个盘的质量分别是 12 kg 和 18 kg两个轴承刚度都取 2×10^7 N/m。参数数值说明轴段长度0.6 / 0.6 m两段抗弯刚度 EI6.45e4 N·m²按实心圆轴计算盘质量12 / 18 kg集中质量轴承刚度2e7 N/m滚动轴承量级扫频范围0~2000 rad/s覆盖前三阶弯曲用上述代码跑出来的结果大致是一阶临界转速 650~750 rad/s二阶 1400~1600 rad/s 量级具体值随段数和刚度变化。段数从每根轴 1 段加密到 10 段时一阶值的变化会小于 1%这就是收敛的标志。计算前先确认单位长度用米质量用千克力用牛ω 用 rad/s得到的转速要转成 rpm 时乘 (60 / (2\pi))。4. 参数怎么设根才不漏步长、端部条件与陀螺效应4.1 扫频步长和漏根的三个对策临界转速漏根是最常见的问题。dw太大时两个相邻根之间残余量两次变号有可能被一步跨过残差符号没有变化根就丢了。三个对策第一扫频前先想清楚要几阶。发动机转子工作转速往往在 3000~15000 rpm换算是 314~1570 rad/s至少要保证一阶弯曲、二阶弯曲在扫频范围内所以上限取到 2000 rad/s 比较合适步长先取 1 rad/s。第二符号判断加密度检查。如果某一段里残差绝对值出现很小值比如小于最大残差的 1e-6但符号不变要把它当作可疑根重新细化。因为真根两侧残差必然异号如果符号连续接近零多半是模态密集或数值溢出需要加密该区间。第三扫频结果和细化结果互相验证。二分细化时保留 60 次迭代理论上能让区间缩小到 (1/2^{60})实际不用那么多但可以顺便输出细化区间的大小如果最终区间仍然大于 1e-3说明该处不是干净的单根要怀疑是数值震荡。4.2 端部条件怎么选悬臂、简支还是弹性支承传递矩阵法里端部条件决定了残余量的具体表达式。上面的代码用的是自由-自由边界适用于风扇转子在地面试车时两端不约束的状态。但发动机实际安装在飞机上时转子通过轴承座连接到静子机匣此时两端边界应当取弹性支承。具体做法是给两端节点配置轴承刚度而不是把边界条件设成自由。对所有节点统一给定轴承刚度后残余行列式表达式不变右端自由状态自然由大刚度或小刚度区分出来。支承刚度的数量级对临界转速影响极大。一个常用的工程参考滚动轴承刚度在 1e7~1e8 N/m挤压油膜阻尼器等效刚度大约在 1e6~1e7 N/m。算之前如果不能确定刚度就给一组低刚度、一组高刚度分别计算得到临界转速的“散布带”再和测振数据对照。支承刚度取低了临界转速整体下降取高了转子趋向刚性支承临界转速逼近轴本身固有频率。这是设计阶段的常规灵敏度分析和传递矩阵法的计算效率正好匹配。4.3 陀螺力矩让临界转速怎么漂轴上的盘不仅有质量还有极转动惯量。盘偏转时会产生陀螺力矩等效为一个和转角相关的弯矩作用于弯矩平衡方程。这个效应在低速下可以忽略但航空发动机转速高、盘径大必须考虑。最直接的做法是给点矩阵加一项def point_matrix_gyro(m, Jp, w, Omega, k0.0, sign1.0): 考虑盘极转动惯量 Jp 的点矩阵 Omega 为自转转速w 为涡动频率 同步正进动时 Omega wsign 取 1.0 P np.eye(4) P[3, 0] m * w * w - k P[2, 1] sign * Jp * Omega * w # 陀螺力矩项 return P把P[2,1]从 0 改成 (J_p \Omega \omega)表示弯矩受转角影响。同步正进动时转子自转方向和涡动方向一致陀螺力矩使转子表现为“变硬”临界转速比不考虑陀螺项时升高反进动时则相反临界转速降低。这就是为什么坎贝尔图里每一条“临界转速线”实际上是两条斜线交叉出来的。实际发动机转子是多盘结构每个盘的极转动惯量必须从三维模型或称重摆测得到。圆盘极转动惯量近似为 (J_p m R^2/2)如果手里只有盘质量和半径可以先用这个估算。加入陀螺项后残余量不再只是 ω 的函数而是 (\omega) 和 (\Omega) 的二元函数计算时先固定自转转速 (\Omega) 为工作转速再扫频找 (\omega)两者相等时对应同步正进动临界转速。4.4 数值陷阱病态矩阵与单位一致性传递矩阵法最隐蔽的问题是数值尺度。矩阵里同时存在 (l^3/EI) 和 (m \omega^2)当轴很长或 ω 很大时矩阵元素数量级差可达 10^{10}浮点截断会把小量吃掉表现出来就是一阶根能算对、二阶根误差大、三阶根乱跳。对策有三个。一是把单位放在毫米、吨、千牛、千分之秒上让 EI、质量、频率的数量级尽量贴近 1再换算回工程单位。二是加密分段但不要无限加密40~100 段足够超过 200 段后矩阵连乘的点数增加逐步累积误差反而比网格加密的收益大。三是如果高频段始终不稳定改用 Riccati 传递矩阵或阻抗耦合法它们把矩阵划成子块损失更小。单位一致性要单独检查。常见错误包括弹性模量用 GPa10^9 Pa乘惯性矩时直接相乘得到错误数量级质量用克而密度用千克每立方米转速用 rpm 却当成 rad/s 传入扫频函数。建议在程序开头把所有输入转换成统一的 SI 基准单位再进入矩阵计算转换过程写成函数而不是散落在各段代码里出问题能一次定位。5. 用解析解和坎贝尔图验证传递矩阵法算出的临界转速单盘悬臂转子是验证传递矩阵代码最好的标定算例。一根长度 L0.8 m、抗弯刚度 EI 的无质量轴末端一个质量 m20 kg不考虑轴质量时一阶临界转速有精确解[ \omega_c \sqrt{\frac{3EI}{m L^3}} ]假设 EI 取 6.45×10^4 N·m²那么理论值为 ( \sqrt{3 \times 64500 / (20 \times 0.8^3)} ) ≈ 137.6 rad/s换算成转速约 1314 rpm。把上面的rotor列表写成[(0.8, 64500, 20, 0)]算出来的根应当和理论值一致到十进制精度。如果发现差了很多先检查左端边界条件是否把自由端设成了悬臂固定端再检查残余行列式是否取了正确的行和列。第二个验证手段是加密分段。对任意工程转子把每轴段从 1 段改成 4、8、16、32 段各算一次一阶临界转速。正常情况是段数加密后结果单调收敛相邻两次加密差值不断缩小如果 8 段到 16 段的变化大于 1%说明模型还没收敛要么继续加密要么检查有没有轴段跨过直径突变处。直径突变的位置务必设置为节点否则场矩阵用错 EI收敛速度会明显变慢。验证完数值还要验证物理趋势把盘质量增加 10%临界转速应该大致下降 (\sqrt{1/1.1} \approx 0.953) 倍把轴承刚度减半弹性支承转子的临界转速下降下降幅度越大说明该阶模态受支承影响越强。这个灵敏度测试不需要额外解析解却能快速暴露矩阵拼装错误。如果关心转速变化下的临界转速变化可以画坎贝尔图。固定自转转速 (\Omega) 为一组值比如 0 到 2000 rad/s 分 21 档对每个 (\Omega) 用第 4.3 节的二元扫频找同步正进动频率把得到的 (\omega) 与 (\Omega) 画在横竖坐标相同的图上和 45 度线相交的点就是实际临界转速。整理出来的座标直接用于发动机适航取证中的转速裕度评估。最后一个实用技巧算完临界转速后把左右支承刚度同时降到原来的一半再算一遍比较前两阶变化。如果移动量超过 5%说明转子由支承主导调整支承刚度比改轴径更有效如果几乎不动说明转子由自身弯曲刚度主导应该优先改轴径、壁厚或材料。这一条比任何理论判据都更快地指导结构修改方向。本文还有配套的精品资源点击获取 
长河君 · 2026/9/16 1:46:05 · 阅读 1.2k 传递矩阵法计算转子临界转速：从矩阵乘法到Python实现
简介面向航空发动机转子动力学分析与设计验证场景这份资源提供基于传递矩阵法的临界转速计算工具核心是一个轻量MATLAB脚本。脚本通过输入轴段长度、截面参数、材料属性与支承刚度等模型信息即可建立传递矩阵并求解转子固有频率进而得到各阶临界转速帮助工程师在概念设计阶段快速评估共振风险。传递矩阵法将复杂连续轴系离散为若干单元编码简洁、计算开销小脚本完整展示了从模型离散、矩阵构建到频率求解的代码实现对学习转子动力学和MATLAB数值方法均有参考价值。压缩包共1个文件文件类型为m脚本整个资源包仅1KB轻巧易用。目前已有237人学习适合航空发动机转子、透平机械等旋转机械的临界转速快速估算也可作为课程设计或进一步开发轴承-转子系统分析程序的基础模板无论是初学者理解传递矩阵流程还是工程师快速获得基准结果都能提供直观支持。1. 传递矩阵法算转子临界转速先搞懂这一个矩阵乘法航空发动机转子的临界转速计算听起来是有限元的事但工程一线最常见的资料包就是那类rotor_criticalspeed_calculation.rar解开后一个 Fortran 文件或 MATLAB 脚本核心算法是传递矩阵法。原因是设计初期要反复改盘的位置、轴承刚度、轴径有限元建模一次要 20 分钟传递矩阵法把转子切成几十段每段一个 4×4 矩阵乘法连乘一次只要几毫秒扫一遍频率就能找到各阶临界转速。本文不依赖某个商业软件直接讲清状态向量、场矩阵、点矩阵怎么拼给一套 Python 实现参数设置和漏根排错也会展开。适合刚接转子动力学任务的强度工程师以及想把手算例题程序化的学生。2. 传递矩阵法建模状态向量与场/点矩阵2.1 转子截面的 4 个状态分量怎么排传递矩阵法把一个转子从小端到大端切成若干段每一段的连接截面用一个状态向量描述。对横向弯曲状态向量取[ Z [y,\ \theta,\ M,\ T]^T ]其中 y 是挠度θ 是转角M 是弯矩T 是剪力。整个转子的振动问题变成从一个端面出发逐段把状态向量乘到另一端最后满足两端边界条件的频率就是临界转速。这个思路和弹塑性梁的初参数法同源只是把初参数法的双曲函数解换成数值矩阵。正负号约定是代码对不上的第一个原因。弯矩和剪力的符号必须和材料力学里梁的内力符号约定一致弯矩逆时针为正剪力使左段向上右段向下为正。如果从某份源码里抄的矩阵把 T 的符号反了算出来的临界转速频率规律是对的但振型会上下镜像而且和支承耦合时误差很隐蔽。拿到任何传递矩阵代码第一步永远是把符号约定写进注释再拿一个单盘悬臂转子的解析解做标定。切段规则按转动节点来。两端的轴段、盘、轴承都落在离散节点上盘和轴承只占据一个节点的位置不占用一段轴长。对航空发动机转子压气机盘和涡轮盘比同段轴至少重一个量级这样集中质量的处理足够准确叶片离心刚化效应本来也不在临界转速计算的基本公式里后续修正时再单独加项。2.2 场传递矩阵无质量轴段的 4×4 传递关系节点 i 到节点 i1 之间有一段等截面圆轴长度 l抗弯刚度 EI。先忽略这段轴本身的分布质量它只提供弹性变形。从截面 i 到截面 i1状态向量满足[ Z_{i1} F_i , Z_i ][ F_i \begin{bmatrix} 1 l \dfrac{l^2}{2EI} \dfrac{l^3}{6EI}\ 0 1 \dfrac{l}{EI} \dfrac{l^2}{2EI}\ 0 0 1 l\ 0 0 0 1 \end{bmatrix} ]这个矩阵的每一列都对应一个初始状态分量对末端各状态量的影响。例如第一列表示当截面 i 处 y1、其他状态量为 0 时经过这段悬臂式变形末端产生挠度、转角、弯矩和剪力分别是矩阵第一列的四个元素。第二列以上依次类推。矩阵是低三角的说明弯矩只由剪力和外载荷产生挠度则由上一截面的所有状态累积而成符合简单的梁变形几何关系。注意这个 F 矩阵里没有频率 ω也没有质量项。这是把轴段质量“砍掉”后的简化场矩阵。相比直接采用带分布质量的精确解里面会出现 cos/sinh 组合分段集中质量的做法在段数足够时误差只推向高频段而实际关心的低阶临界转速正好在低频段。因此工程代码里大量使用这种无质量场矩阵加集中质量点矩阵的组合计算速度快物理意义清楚。如果你要算第三阶以上建议把每段再加密一倍而不是换更复杂的场矩阵。2.3 点传递矩阵把盘和轴承压进一个 4×4 矩阵轴段经过节点时节点上的集中质量 m 和轴承刚度 k 会对状态向量产生一个突变。设节点左侧状态为 (Z_l)右侧为 (Z_r)挠度、转角、弯矩在节点两侧连续只有剪力因惯性力和支承反力发生跳跃。假设转子以频率 ω 做同步正进动惯性力幅值为 (m\omega^2 y)方向与挠度一致轴承反力为 (-k y)。于是[ T_r T_l m\omega^2 y - k y ]写成矩阵[ P \begin{bmatrix} 1 0 0 0\ 0 1 0 0\ 0 0 1 0\ m\omega^2 - k 0 0 1 \end{bmatrix} ]点在矩阵右下角是 1表示剪力从节点左侧平移到右侧左下角那个 (m\omega^2-k) 是响应项。没有盘的节点直接让 m0没有支承的节点让 k0。把某一小段轴和其后置的节点合在一起就是一级基本传递[ T_i P_i F_i ]从最左端开始连乘到最右端[ T_{total} P_n F_n \cdots P_1 F_1 ]如果左端状态为 (Z_0)右端就是 (Z_n T_{total} Z_0)。左端边界条件会消掉两个未知数右端边界条件又给出两个方程最后只留下一个关于 ω 的残余量。残余量为零时ω 就是临界转速。这一整套过程里没有求解大型线性方程组矩阵永远只有 4×4 连乘这就是它快的原因。3. 实现用 Python 写一个能跑的转子临界转速计算程序3.1 把转子的几何和材料拆成 segments/masses 两个表开始写代码前先把模型数据整理成两张表。第一张表是轴段列表每段含长度 l 和抗弯刚度 EI第二张表是节点属性每个节点含集中质量 m 和支承刚度 k。以下代码以一个单盘悬臂转子和一个双支承双盘转子为例。import numpy as np def field_matrix(l, EI): 无质量轴段的场传递矩阵 return np.array([ [1.0, l, l * l / (2.0 * EI), l ** 3 / (6.0 * EI)], [0.0, 1.0, l / EI, l * l / (2.0 * EI)], [0.0, 0.0, 1.0, l], [0.0, 0.0, 0.0, 1.0] ]) def point_matrix(m, w2, k0.0): 节点上的点传递矩阵w2 为 omega^2 P np.eye(4) P[3, 0] m * w2 - k return Pfield_matrix对应 2.2 节的 F 矩阵point_matrix对应 P 矩阵。注意point_matrix里参数名是w2调用时要传 (\omega^2)不要传成 (\omega)否则量纲全错。轴段的数据结构用列表里的元组表示每个元组是(l, EI, m, k)顺序从左端到右端。m 和 k 是该段右侧节点的质量与轴承刚度。这种数据排法最直观也方便从 Excel 或旧 Fortran 数据文件转换。3.2 扫频找残余行列式过零点的完整代码以一个左端自由、右端自由的转子为例。左端没有弯矩和剪力所以 (M_00, T_00)未知的只有 (y_0) 和 (\theta_0)。右端同样自由要求 (M_n0, T_n0)。由 (Z_n T_{total} Z_0)右端弯矩、剪力分别是总矩阵第三、四行的前两列与 (y_0, \theta_0) 的组合def total_matrix(rotor, w): 连乘所有段的场矩阵和点矩阵 T np.eye(4) for l, EI, m, k in rotor: T point_matrix(m, w * w, k) field_matrix(l, EI) T # 先过场矩阵再过点矩阵顺序和节点布置一致 return T def residual(rotor, w): 自由-自由边界条件下的残余行列式 T total_matrix(rotor, w) # 右端弯矩、剪力仅依赖于左端 y0, theta0 # 非零解条件为 2x2 行列式等于 0 return T[2, 0] * T[3, 1] - T[2, 1] * T[3, 0]残余量随 ω 变化每过一次零点就对应一个临界转速。扫频函数从 0 开始小步长推进符号变化后原地二分细化把根锁定到 10^{-6} rad/s 精度def find_critical_speeds(rotor, w_max1000.0, dw1.0): roots [] w 0.0 r0 residual(rotor, w) while w w_max: r1 residual(rotor, w dw) if r0 0.0: roots.append(w) elif r0 * r1 0: # 跨越零点的粗判断 a, b w, w dw for _ in range(60): mid 0.5 * (a b) rm residual(rotor, mid) if rm * r0 0: b mid else: a mid r0 rm roots.append(0.5 * (a b)) w dw r0 residual(rotor, w) return rootsr0 * r1 0是常见的异号判断但如果残余量在某个频率附近有上下尖峰而不真正过零这种判断会丢根或误报。更稳的做法是先画一条残余量曲线看形状再决定步长。代码里的dw1.0意味着每 1 rad/s 采样一次对通常几百 rad/s 的低阶临界转速够密对频率间隔小于 1 rad/s 的密集模态则不够后面第 4 章会讲步长怎么选。3.3 临界转速落点参数表与计算精度给一个双支承双盘转子的算例模拟简化后的发动机高压转子总长 1.2 m两段轴各 0.6 m轴径 50 mm材料弹性模量按 GH4169 取 210 GPa抗弯刚度约为 6.45×10^4 N·m²。两个盘的质量分别是 12 kg 和 18 kg两个轴承刚度都取 2×10^7 N/m。参数数值说明轴段长度0.6 / 0.6 m两段抗弯刚度 EI6.45e4 N·m²按实心圆轴计算盘质量12 / 18 kg集中质量轴承刚度2e7 N/m滚动轴承量级扫频范围0~2000 rad/s覆盖前三阶弯曲用上述代码跑出来的结果大致是一阶临界转速 650~750 rad/s二阶 1400~1600 rad/s 量级具体值随段数和刚度变化。段数从每根轴 1 段加密到 10 段时一阶值的变化会小于 1%这就是收敛的标志。计算前先确认单位长度用米质量用千克力用牛ω 用 rad/s得到的转速要转成 rpm 时乘 (60 / (2\pi))。4. 参数怎么设根才不漏步长、端部条件与陀螺效应4.1 扫频步长和漏根的三个对策临界转速漏根是最常见的问题。dw太大时两个相邻根之间残余量两次变号有可能被一步跨过残差符号没有变化根就丢了。三个对策第一扫频前先想清楚要几阶。发动机转子工作转速往往在 3000~15000 rpm换算是 314~1570 rad/s至少要保证一阶弯曲、二阶弯曲在扫频范围内所以上限取到 2000 rad/s 比较合适步长先取 1 rad/s。第二符号判断加密度检查。如果某一段里残差绝对值出现很小值比如小于最大残差的 1e-6但符号不变要把它当作可疑根重新细化。因为真根两侧残差必然异号如果符号连续接近零多半是模态密集或数值溢出需要加密该区间。第三扫频结果和细化结果互相验证。二分细化时保留 60 次迭代理论上能让区间缩小到 (1/2^{60})实际不用那么多但可以顺便输出细化区间的大小如果最终区间仍然大于 1e-3说明该处不是干净的单根要怀疑是数值震荡。4.2 端部条件怎么选悬臂、简支还是弹性支承传递矩阵法里端部条件决定了残余量的具体表达式。上面的代码用的是自由-自由边界适用于风扇转子在地面试车时两端不约束的状态。但发动机实际安装在飞机上时转子通过轴承座连接到静子机匣此时两端边界应当取弹性支承。具体做法是给两端节点配置轴承刚度而不是把边界条件设成自由。对所有节点统一给定轴承刚度后残余行列式表达式不变右端自由状态自然由大刚度或小刚度区分出来。支承刚度的数量级对临界转速影响极大。一个常用的工程参考滚动轴承刚度在 1e7~1e8 N/m挤压油膜阻尼器等效刚度大约在 1e6~1e7 N/m。算之前如果不能确定刚度就给一组低刚度、一组高刚度分别计算得到临界转速的“散布带”再和测振数据对照。支承刚度取低了临界转速整体下降取高了转子趋向刚性支承临界转速逼近轴本身固有频率。这是设计阶段的常规灵敏度分析和传递矩阵法的计算效率正好匹配。4.3 陀螺力矩让临界转速怎么漂轴上的盘不仅有质量还有极转动惯量。盘偏转时会产生陀螺力矩等效为一个和转角相关的弯矩作用于弯矩平衡方程。这个效应在低速下可以忽略但航空发动机转速高、盘径大必须考虑。最直接的做法是给点矩阵加一项def point_matrix_gyro(m, Jp, w, Omega, k0.0, sign1.0): 考虑盘极转动惯量 Jp 的点矩阵 Omega 为自转转速w 为涡动频率 同步正进动时 Omega wsign 取 1.0 P np.eye(4) P[3, 0] m * w * w - k P[2, 1] sign * Jp * Omega * w # 陀螺力矩项 return P把P[2,1]从 0 改成 (J_p \Omega \omega)表示弯矩受转角影响。同步正进动时转子自转方向和涡动方向一致陀螺力矩使转子表现为“变硬”临界转速比不考虑陀螺项时升高反进动时则相反临界转速降低。这就是为什么坎贝尔图里每一条“临界转速线”实际上是两条斜线交叉出来的。实际发动机转子是多盘结构每个盘的极转动惯量必须从三维模型或称重摆测得到。圆盘极转动惯量近似为 (J_p m R^2/2)如果手里只有盘质量和半径可以先用这个估算。加入陀螺项后残余量不再只是 ω 的函数而是 (\omega) 和 (\Omega) 的二元函数计算时先固定自转转速 (\Omega) 为工作转速再扫频找 (\omega)两者相等时对应同步正进动临界转速。4.4 数值陷阱病态矩阵与单位一致性传递矩阵法最隐蔽的问题是数值尺度。矩阵里同时存在 (l^3/EI) 和 (m \omega^2)当轴很长或 ω 很大时矩阵元素数量级差可达 10^{10}浮点截断会把小量吃掉表现出来就是一阶根能算对、二阶根误差大、三阶根乱跳。对策有三个。一是把单位放在毫米、吨、千牛、千分之秒上让 EI、质量、频率的数量级尽量贴近 1再换算回工程单位。二是加密分段但不要无限加密40~100 段足够超过 200 段后矩阵连乘的点数增加逐步累积误差反而比网格加密的收益大。三是如果高频段始终不稳定改用 Riccati 传递矩阵或阻抗耦合法它们把矩阵划成子块损失更小。单位一致性要单独检查。常见错误包括弹性模量用 GPa10^9 Pa乘惯性矩时直接相乘得到错误数量级质量用克而密度用千克每立方米转速用 rpm 却当成 rad/s 传入扫频函数。建议在程序开头把所有输入转换成统一的 SI 基准单位再进入矩阵计算转换过程写成函数而不是散落在各段代码里出问题能一次定位。5. 用解析解和坎贝尔图验证传递矩阵法算出的临界转速单盘悬臂转子是验证传递矩阵代码最好的标定算例。一根长度 L0.8 m、抗弯刚度 EI 的无质量轴末端一个质量 m20 kg不考虑轴质量时一阶临界转速有精确解[ \omega_c \sqrt{\frac{3EI}{m L^3}} ]假设 EI 取 6.45×10^4 N·m²那么理论值为 ( \sqrt{3 \times 64500 / (20 \times 0.8^3)} ) ≈ 137.6 rad/s换算成转速约 1314 rpm。把上面的rotor列表写成[(0.8, 64500, 20, 0)]算出来的根应当和理论值一致到十进制精度。如果发现差了很多先检查左端边界条件是否把自由端设成了悬臂固定端再检查残余行列式是否取了正确的行和列。第二个验证手段是加密分段。对任意工程转子把每轴段从 1 段改成 4、8、16、32 段各算一次一阶临界转速。正常情况是段数加密后结果单调收敛相邻两次加密差值不断缩小如果 8 段到 16 段的变化大于 1%说明模型还没收敛要么继续加密要么检查有没有轴段跨过直径突变处。直径突变的位置务必设置为节点否则场矩阵用错 EI收敛速度会明显变慢。验证完数值还要验证物理趋势把盘质量增加 10%临界转速应该大致下降 (\sqrt{1/1.1} \approx 0.953) 倍把轴承刚度减半弹性支承转子的临界转速下降下降幅度越大说明该阶模态受支承影响越强。这个灵敏度测试不需要额外解析解却能快速暴露矩阵拼装错误。如果关心转速变化下的临界转速变化可以画坎贝尔图。固定自转转速 (\Omega) 为一组值比如 0 到 2000 rad/s 分 21 档对每个 (\Omega) 用第 4.3 节的二元扫频找同步正进动频率把得到的 (\omega) 与 (\Omega) 画在横竖坐标相同的图上和 45 度线相交的点就是实际临界转速。整理出来的座标直接用于发动机适航取证中的转速裕度评估。最后一个实用技巧算完临界转速后把左右支承刚度同时降到原来的一半再算一遍比较前两阶变化。如果移动量超过 5%说明转子由支承主导调整支承刚度比改轴径更有效如果几乎不动说明转子由自身弯曲刚度主导应该优先改轴径、壁厚或材料。这一条比任何理论判据都更快地指导结构修改方向。本文还有配套的精品资源点击获取