船舶轴系振动MATLAB程序：传递矩阵法贯穿扭振计算与控制

简介面向船舶轴系振动与控制分析的MATLAB程序包适合船舶、机械、自动化及计算机等专业学生用于课程设计、期末作业与毕业设计也可供相关工程技术人员进行仿真验证。程序支持MATLAB 2014、2019a与2024a等版本采用参数化编程参数便于修改注释详细并附赠可直接运行的案例数据替换数据即可快速开展不同工况下的振动分析。包体共57个文件以34个m源码文件为核心覆盖纵向振动与扭转振动等分析模块另有10个fig界面文件、10张jpg示意图、2个txt说明文档及1个asv备份文件压缩包整体约352KB结构清晰便于检索。目前已有121人学习下载。借助该程序用户可直观学习船舶轴系动力学建模、振动特性计算与减振控制思路在Simulink等工具环境下快速上手是兼顾理论理解与工程实践的高性价比学习资料。1. 船舶轴系振动与控制分析MATLAB程序机舱工程师的扭振计算工具箱船舶轴系的振动计算在实船项目里从来不是“可做可不做”的选做题。主机厂、设计所、船级社三方审图时扭振计算报告是轮机部分必查项而纵向振动和回旋振动在大型低转速螺旋桨和桨-艉部耦合激励越来越强的今天也常被要求补算。这套MATLAB程序就是把轴系离散化建模、传递矩阵法求固有频率、强迫响应计算和主动控制仿真打包在一起能直接拿轴系参数换出来一组可审图的数值结果。程序不依赖额外工具箱代码量也不大适合轮机工程专业做课程设计和毕业设计也适合船厂及设计所的工程师拿来校核外协报告。2. 轴系离散化与传递矩阵建模先弄懂三类振动和程序的数据结构2.1 三类振动先分清扭转、纵向、回旋的求解路径不同轴系振动有个老规矩先分方向再谈方法。扭振发生在回转方向激励源主要是气缸内燃气压力和往复惯性力引起的周期性转矩路径从曲轴、弹性联轴器一路传到螺旋桨纵向振动沿轴线方向激励来自螺旋桨脉动推力轴系作为弹性杆件被反复拉伸压缩回旋振动则是轴系在自身重力、螺旋桨偏心力和轴承油膜力作用下产生的横向涡动常和艉轴承磨损、轴颈偏心联系起来。这三类振动在程序里对应三种不同的数学模型。扭振用集中质量-无质量弹性轴段模型每个节点有一个转动惯量每两个节点之间用一个扭转刚度连接纵向振动按轴向杆件离散每个节点有质量和轴向刚度回旋振动则要引入轴承支撑刚度把轴系当成转子-轴承系统处理。选传递矩阵法而不是直接开有限元是因为轴系是典型的一维链式结构传递矩阵法迭代快、物理意义清楚、溯源容易有限元能做但建模工作量翻倍程序也重得多。对课程设计和工程复核来说传递矩阵法足够交代清楚问题。2.2 程序的数据结构轴段、节点、激励表的组织方式打开这套程序第一个要读的脚本是类似input_axis.m的参数输入文件。程序里把轴系划分成从主机飞轮端或曲轴自由端到螺旋桨的若干轴段每个节点对应一个集中惯量或质量每段对应一个刚度。这个划分不是等分就完事而是要按轴系实际结构“停在关键位置”弹性联轴器两端、齿轮啮合点、推力轴承位置、螺旋桨锥部、每段中间法兰这些位置都必须单独设节点。% 轴系节点与轴段参数定义SI单位制 % J: 节点转动惯量 kg*m^2K: 相邻节点间扭转刚度 N*m/rad nodes struct( ... name, {Flywheel, Coupling_Shaft, Coupling_Eng, Shaft_1, Shaft_2, Propeller}, ... J, [12.8, 0.15, 1.6, 2.1, 2.4, 38.6], ... % 螺旋桨带附水惯量 K, [0, 4.2e5, 1.08e5, 6.7e5, 7.1e5, 0]); % 首尾两端的K用0占位这里J的单位是kg*m^2K的单位是N*m/rad。程序中所有输入必须保持同一套单位制我一般统一用 SI 制避免中间手工换算。注释里特意标了“螺旋桨带附水惯量”这个值通常是螺旋桨干惯量的 1.15~1.3 倍随船舶状态营运/压载还会变程序不写死方便按工况修改。激励输入是另一张表。扭振强迫响应需要的是各转速下的谐次激励力矩幅值。一般把每缸气体压力曲线做傅里叶展开取主要谐次通常取 0.5、1.0、1.5、2.0、3.0 谐次低速机还要看到 6 谐次幅值乘以作用半径得到力矩相位按发火顺序错开。这些在程序里也是一个矩阵行是转速列是谐次。注意激励相位差不是可选项同一谐次下各缸力矩如果相位设置错误合成激励会互相抵消或虚增后文避坑章会单独讲。2.3 传递矩阵主循环从自由端起算的MATLAB实现传递矩阵法的套路是把轴系从一端比如自由端开始逐步把状态向量角位移和力矩往下传递另一端的边界条件成立时对应的频率就是固有频率。程序里的核心循环本质上是把每个轴段的“场传递矩阵”和每个节点的“点传递矩阵”顺序乘起来。% 扭振传递矩阵法计算给定频率 omega 下的总传递矩阵 function T_total torsional_transfer_matrix(nodes, omega) n numel(nodes.J); T eye(2); % 状态向量 [theta; Torque] for i 2:n K nodes.K(i); J nodes.J(i); % 场传递矩阵轴段柔度引起的转角变化 Tf [1, 1/K; 0, 1]; % 点传递矩阵集中惯量在简谐运动下产生的惯性扭矩 -J*omega^2 Tp [1, 0; -J*omega^2, 1]; T T * Tf * Tp; end T_total T; end这个循环对应的是自由端-自由端边界条件。Tf里的1/K是轴段柔度表示“这段轴在单位扭矩下的转角变化”Tp里的-J*omega^2是集中质量在简谐运动下产生的惯性扭矩方向与加速度相反。循环从i2开始是因为首节点只作为边界点提供惯量不参与传递。调用时给不同的omega扫描找到让边界条件Torque0或theta0的频率就是固有频率。这个写法同样适用于纵向振动不过状态向量变成轴向位移和轴向力Tf中的柔度变成L/(E*A)Tp中的惯性项变成-m*omega^2其中m是节点集中质量E*A是轴段轴向刚度。回旋振动则要把两个平面水平和垂直耦合起来状态向量扩成四阶代码结构不变矩阵阶数变高。边界条件的选择要跟着轴系实际约束走。曲轴自由端通常按自由端处理螺旋桨端也近似自由端如果轴系里有齿轮箱作为固定点或者推力轴承在纵振模型里表现为大刚度弹簧对应节点的边界条件就要从“自由”改成“弹性约束”。程序默认是自由-自由用其他边界条件时要在主循环里补一个约束矩阵。3. 固有频率与强迫响应计算结果与船级社许用值的对比逻辑3.1 用特征值法求扭振固有频率eig与Holzer校验传递矩阵法扫描频率直观但频率点要自己找找没找到某个根、漏没漏根心里没底。程序里通常还会配一个特征值法做交叉验证把系统的质量矩阵和刚度矩阵组出来直接调用eig。这个矩阵不大——节点数一般 20~50 个手写组装不要多少时间。% 组装质量矩阵与刚度矩阵用特征值法求固有频率 J nodes.J(2:end-1); % 去掉首尾占位 Ktotal nodes.K(2:end-1); M diag(J); K zeros(numel(J)); for i 1:numel(J)-1 k Ktotal(i1); K(i,i) K(i,i) k; K(i1,i1) K(i1,i1) k; K(i,i1) K(i,i1) - k; K(i1,i) K(i1,i) - k; end [V, D] eig(K, M); % 广义特征值问题 f_natural sqrt(abs(diag(D))) / (2*pi); % 单位 Hzeig(K, M)求解的是广义特征值问题特征值的平方根对应圆频率除以2*pi才是工程上常用的 Hz。程序里这一步算出来的应该是前几阶低于激振频率范围的扭振模态实际审图时主要关心 1 节点和 2 节点模态更高阶模态要么远离工作转速范围要么被阻尼压得不重要了。Holzer 表是传统手算方法程序里通常保留这个函数作为“物理可解释”的对照。每给定一个假设频率从自由端逐段累加惯量力矩和轴段转角最后看末端封闭条件是否近似成立。和传递矩阵法相比它更像逐级迭代的试算法当eig和 Holzer 结果差超过 2% 时我会先怀疑边界条件设错了而不是急着改程序。eig内部用的是广义特征值分解对矩阵规模不大但刚性比较大的系统直接用不带平移的默认算法就可能丢根我会先做一次对角缩放再求解。3.2 强迫响应计算激励力矩谱与放大系数固有频率算完紧接着算强迫响应。程序里对每一个工作转速把各谐次的激励力矩施加在对应节点上做频率扫描求解稳态响应幅值。这里面有个工程上容易忽略的点激励源在曲轴端响应点可能在螺旋桨端两者之间的相位差靠传递矩阵的复数形式体现所以程序里的刚度、惯量都按复数动态刚度处理而不是静力关系。% 在给定转速 n_rpm 下求螺旋桨端扭振幅值 % excite_node: 激励作用节点编号谐次 h激振力矩幅值 T_amp function theta_p forced_response(nodes, n_rpm, h, T_amp) omega 2*pi*(n_rpm/60)*h; % 激励圆频率 T_total torsional_transfer_matrix(nodes, omega); r T_total(1,1) / T_total(1,2); % 自由端边界消元 theta_p T_amp / (r * nodes.K(end-1)); end这段代码最核心的是T_total(1,1) / T_total(1,2)这一行。它的物理意义是在自由端边界下系统在激励频率处的动柔度单位力矩产生的转角。当激励频率等于固有频率时动柔度趋于无穷响应发散加了阻尼之后分母出现虚部响应被限制在有限值。程序中阻尼通常不在主循环里出现而是在最后一步以损失系数或对数衰减率折算进柔度。很多人第一次算强迫响应会问为什么共振点幅值那么大甚至超过许用应力因为程序初始阻尼取的是材料内阻尼而轴系里弹性联轴器、减振器、液力耦合器的阻尼远大于材料阻尼。这些部件的阻尼要按厂家给出的损失系数单独加进对应节点的复刚度程序输入文件里有专门一个字段loss_factor就是干这个的。频率扫描的范围要覆盖主机最低稳定转速到 105% 额定转速步长取每 5 rpm 一个点共振峰附近自动加密到 1 rpm否则峰值容易被步长漏过去。3.3 结果整理如何对齐船级社许用应力算完响应幅值还要换算成剪应力才谈得上校核。扭振剪应力等于节点两侧扭矩差除以该轴段的抗扭截面模量公式是tau T_t / W_t其中W_t pi*D^3*(1 - (d/D)^4)/16。程序在后处理脚本里自动完成这个过程并和船级社规范的许用值表对比。校核位置许用应力来源常见取值区间MPa主机曲轴自由端主机厂/船级社规范连续运转 ±30~±50中间轴船级社规范标称值按轴径换算螺旋桨轴锥部船级社规范加应力集中系数略低于中间轴弹性联轴器两侧联轴器厂家许用按转速和扭振附加应力这块最容易翻车的是应力集中系数。轴径变化、键槽、法兰根部这些位置的局部应力远高于名义应力程序里给的名义剪应力不能直接对标许用值要么乘 2~3 的应力集中系数要么直接用有限元子模型算局部应力。程序包的说明文档里明确提醒了这一点但不少用户还是直接把后处理表格里的应力拿去对标造成误判。转速扫描之后程序会输出“各转速下最大剪应力包络线”这组数据要逐点查不能只看峰值因为某谐次在低于额定转速时也可能越过许用线。4. 纵向振动与回旋振动两个容易被忽略的维度4.1 纵向振动模型螺旋桨脉动推力激励下的位移响应扭转振动是整个审图流程的主角但这几年纵向振动越来越被重视。大型集装箱船和 VLCC 采用高效率低转速螺旋桨桨叶进入艉部不均匀伴流场时脉动推力直接作用在轴系上激起纵向振动。纵向振动如果和轴系固有频率合拍轻则推力轴承磨损重则艉部结构局部疲劳。程序里的纵向振动模型和扭振几乎同构只是物理量从转角换成位移。% 纵向振动简化多自由度模型M*x C*x K*x F(t) % M: 节点质量含附水质量C: 轴向阻尼K: 轴段轴向刚度 m_vec [950, 120, 180, 210, 420]; % kg螺旋桨节点已含附水质量 k_vec [8.5e8, 9.2e8, 9.0e8, 7.8e8]; % N/m M diag(m_vec); K zeros(5); for i 1:4 K(i,i) K(i,i) k_vec(i); K(i1,i1) K(i1,i1) k_vec(i); K(i,i1) K(i,i1) - k_vec(i); K(i1,i) K(i1,i) - k_vec(i); end F_prop 2.4e5; % 桨叶脉动推力幅值 N按伴流分数估算纵向振动的激励频率是螺旋桨叶频轴频乘以叶片数比如五叶桨的轴频是 1.2 Hz叶频就是 6 Hz高转速时可能更高。程序里输入推力幅值按叶频做简谐激励响应结果输出成推力轴承处的动态力。校核时看动态力比例经验上动态推力幅值不超过平均推力的 5%~10% 是常见目标。纵向振动里附水质量同样存在。螺旋桨在水中前后振荡时周围水体跟着运动等效附水质量可达到螺旋桨自身质量的 30%~50%。程序输入文件里m_vec末尾的420就是干质量加附水质量的合计值。如果这里只填干质量纵向固有频率会被高估 10% 以上叶频共振区的位置就判断错了。实测数据反而不那么依赖附水估算因为实船测的是推力轴承根部加速度频响曲线峰直接从数据里读。4.2 回旋振动临界转速轴承刚度与转子模型的简化回旋振动whirling是轴系作为转子系统的横向振动在艉轴承和后轴承之间的悬臂段尤其明显。轴承磨损、间隙增大、艉轴承后部支撑刚度下降都会让回旋振动的临界转速掉进工作转速范围。程序里的回旋模型做了简化把螺旋桨当成悬臂端的集中质量把艉轴承当成一个横向弹簧轴段用欧拉-伯努利梁单元离散。% 回旋振动临界转速单盘-双轴承简化模型 % EI: 轴段抗弯刚度 N*m^2L: 艉轴承到螺旋桨重心距离 m % k_bearing: 艉轴承支撑刚度 N/mm_prop: 螺旋桨质量 kg EI 2.3e7; L 2.8; k_bearing 3.2e8; m_prop 3500; % 悬臂段等效刚度 k_cantilever 3*EI / L^3; k_total 1 / (1/k_bearing 1/k_cantilever); omega_c sqrt(k_total / m_prop); % 临界角速度 rad/s n_c omega_c * 60 / (2*pi) % 对应转速 rpm这里的重点是k_total的串联合并艉轴承油膜刚度和轴段悬臂弯曲刚度是串联关系总刚度一定小于两者中较小者。实际艉轴承油膜刚度本身还随转速变化程序里提供的是额定转速下的常值要更精细就得耦合约 8 次方量级的油膜 Reynolds 方程那就超出这套程序的范围了。作为快速估算常值刚度加 20% 上下浮动扫描一遍看临界转速落在哪个区间比强行求解油膜来得稳。回旋振动的关键判断是临界转速要避开工作转速的 ±15%或者临界转速对应的振动能量被轴承阻尼有效衰减。如果程序算出来的n_c落在工作转速附近优先调整的是艉轴承支撑刚度而不是轴径因为改轴径对悬臂段弯曲刚度的影响是三次方关系牵一发动全身。另外注意回旋振动还有正进动和反进动的区别程序给的是简化的同步正进动解船用轴系在这个模型下已经够用研究级精度需要把陀螺力矩项加回去。5. 避坑从单位制到控制模型降阶的六个常见坑5.1 单位制混用导致数量级崩盘现象程序跑完固有频率算出来是几千赫兹明显不合理再查某个中间量发现是扭矩单位从 N·m 换成 kN·m 后忘记在计算里换算。原因轴系参数来自不同的原始资料。主机厂的飞轮惯量表给的是kg·m^2联轴器厂家给的扭转刚度是kN·m/rad而轴段几何尺寸用mm一混就乱。解决程序里所有输入文件强制统一到 SI 制进入计算主循环之前加一个assert检查刚度值落在 1e4~1e8 之间、惯量值落在 1e-2~1e3 之间就通过否则直接报错提示“请检查单位制”。我拿到新数据的第一件事就是把所有单位写在同一张表上逐项换算后再填入程序。5.2 分段数太少导致高阶模态缺失现象特征值法只算出 2~3 阶固有频率Holzer 表也找不到更高阶的根或者程序报错“矩阵奇异”。原因轴系离散化节点数太少尤其长中间轴整段只设一个节点等价于把几十米长的弹性轴压成了一个刚性体高阶模态被架空了。解决把长轴段按“每 1~2 米一个节点”的原则重新分段至少保证每根轴段中间有一个节点。分段加密后同一台机器上算出的低阶频率通常变化不大但高阶频率才会出现。对传递矩阵法来说分段数从 10 加到 40计算时间几乎不变所以起步就往密了分。5.3 阻尼设零导致强迫响应发散现象强迫响应在某谐次共振点处幅值冲上天花板画出来的频率响应曲线呈尖峰状无法和实测数据对照。原因程序默认阻尼为零而实际轴系里弹性联轴器、减振器、缸内阻尼共同作用把共振峰值压得很低如果不额外加阻尼算出来的是“理论共振曲线”工程上不可用。解决把联轴器厂家的损失系数loss_factor填进对应节点程序里把该节点的刚度和1 1j*loss_factor相乘形成复刚度。没有厂家数据时先用经验值 0.05~0.1 试算看共振幅值是否降到可以接受的范围。比完后记得把复刚度去掉再算一次固有频率因为阻尼对固有频率的偏移虽然小但会影响 2% 以内的对比精度。5.4 螺旋桨附水惯量没计入导致固有频率偏高现象扭振固有频率比实船测试值高出 8%~15%而且转速越高偏差越大。原因螺旋桨在水中旋转时周围水体的附加惯量没有计入。程序节点表里螺旋桨节点用的是干惯量数值比实际参加振动的惯量偏低固有频率自然偏高。解决给螺旋桨节点乘上附水系数商船螺旋桨通常取 1.15~1.3具体按螺旋桨盘面比和艉部流场取。更稳妥的做法是把压载和满载两种状态的附水系数都算一遍看固有频率的变化范围是否仍然满足避开工作转速区间的要求。5.5 激励相位设置错误导致合成激励失真现象强迫响应算出来的某谐次幅值明显偏小甚至接近零和同型船报告差距很大。原因各缸激励力矩的相位没有按发火顺序错开。同一谐次下各缸贡献的力矩向量本应互相部分抵消相位错一位抵消关系就完全变了。解决程序里建一个发火顺序表按曲柄夹角计算每缸相对相位再叠加成合成激励。最常见的正确做法是第 k 缸相对第 1 缸的相位差等于发火间隔角乘以谐次数。低速机发火间隔不均匀时要逐缸单独算相位不能直接用等间隔公式。5.6 Simulink控制模型与振动计算结果对不上现象把轴系模型搭进 Simulink 做主动控制仿真控制器输出导致系统不稳定或者边仿真边发散调试了很久找不到原因。原因振动计算用的是频域传递矩阵控制仿真用的是时域状态空间两者模型阶数不一样。直接拿 40 阶的完整轴系模型做控制器设计控制器还没调好系统就振荡了而如果随意截断高阶模态模型又丢失了真实的弯扭耦合动态。解决先做模型降阶保留前 2~3 阶主导模态再设计控制器控制器阶数也压到 2~4 阶实现起来才有工程意义。Simulink 里验证通过后还要回到频域的传递矩阵法交叉校验一下降阶模型在某些关键频率点的响应是否一致这一步能排除很多“仿真很漂亮、实船不对”的情况。6. 从振动计算到主动控制Simulink闭环仿真与验证技巧6.1 状态空间模型与LQR控制器的搭建程序的最后一部分是把轴系模型转换成状态空间形式做主动减振控制仿真。适用范围包括纵向振动的主动推力减振以及扭振半主动控制的可行性验证。目标是把选定节点处的振动幅值压下来同时控制输入不超出执行机构的物理限制。% 由质量/刚度/阻尼矩阵生成状态空间模型设计LQR控制器 % 状态量 x [位移; 速度] A [zeros(n,n), eye(n,n); -M\K, -M\C]; B [zeros(n,nu); M\Bu]; % Bu: 控制力作用位置矩阵 C eye(2*n); D zeros(2*n, nu); Q diag([1e6*ones(1,n), 1e2*ones(1,n)]); % 位置权重远大于速度 R 1e-4 * eye(nu); [K_lqr, S, e] lqr(A, B, Q, R);lqr的输出K_lqr是状态反馈增益矩阵e是闭环极点。Q 矩阵里位置项的权重比速度项高四个量级目的是优先抑制位移峰值R 取小值意味着允许控制机构输出较大力矩。实际试算时我会先把R调到 1e-3看控制量峰值是否超过执行机构限制如果超过就把 R 调大牺牲一点减振效果换工程可行性。6.2 验证技巧用频域结果反向校验控制模型主动控制仿真跑通之后还有个容易被跳过的环节把 Simulink 里的闭环模型和前面频域计算的结果做交叉比对。具体做法是给 Simulink 模型加一个扫频正弦激励扫过固有频率附近区域记录响应幅值画成幅频曲线和传递矩阵法算出的强迫响应曲线叠在一起看。% 对降阶模型做扫频验证对比频域传递矩阵结果 f_vec 0.5:0.1:20; % Hz resp_sim zeros(size(f_vec)); for i 1:numel(f_vec) w 2*pi*f_vec(i); resp_sim(i) freqresp(sys_red, w); % sys_red为降阶状态空间模型 end semilogy(f_vec, abs(resp_sim), b-); hold on; semilogy(f_vec, abs(resp_tm), r--); % resp_tm来自传递矩阵法两条曲线在固有频率附近幅值偏差控制在 5% 以内就说明降阶模型保留了足够的动态精度控制器设计可以继续往下走。偏差大了优先检查降阶时是否保留了正确的模态参与因子而不是怀疑控制算法本身。从那以后我每次拿到一套轴系参数都强制走一遍“解析计算→传递矩阵交叉验证→Simulink闭环仿真→实船报告对比”四步流程确认前三步量级一致才敢往下做控制参数整定。这份程序包虽然界面朴素但每一步都看得见、改得了希望帮到你。本文还有配套的精品资源点击获取