简介：这份MATLAB程序包面向船舶轴系的振动与控制分析，适用于计算机、电子信息工程、数学等专业大学生的课程设计、期末大作业与毕业设计，也为需要快速开展振动仿真的初学者提供了低门槛的实践工具。压缩包共57个文件，容量仅352KB，包含34个m源码文件、10个fig图形文件、10张jpg图片、2个txt说明和1个asv备份，源码、界面图与说明文档相互对应，便于按模块查阅和二次修改。目前已有121人学习下载。程序兼容MATLAB 2014、2019a及2024a多个版本，自带可直接运行的案例数据，采用参数化编程，参数可灵活更改，注释详尽，逻辑清晰，初学者也能快速替换数据模拟不同工况。针对船舶轴系纵向振动与扭转振动两个典型方向，程序涵盖建模、仿真到结果可视化流程，可帮助使用者直观理解轴系动力学特性，是完成课程报告或毕业设计的高效辅助资源。
1. 船舶轴系的振动与控制分析MATLAB程序：试航报告上的振动超标，靠什么算清楚
一条新船试航时，转速表刚到某个值，艉轴附近就传来一阵低沉轰鸣，振动测点数据也明显超限。轮机工程师第一反应是“轴系扭振又来了”，但要把该转速下扭振固有频率、激励阶次、临界转速算明白，手边没有现成工具的话，光是列传递矩阵就要一下午。所谓“船舶轴系的振动与控制分析MATLAB程序”，就是把轴系离散建模、固有频率计算、临界转速判别、减振器设计以及主动控制仿真这一整套流程，做成能在MATLAB里直接运行的函数和脚本。它能让设计人员不用等有限元分析一周出结果，也能让研究生快速验证控制算法。适合船厂设计所、船舶与海洋工程专业学生，以及做旋转机械振动监测的人。
2. 轴系振动建模：先把连续轴段离散成矩阵，再给轴承刚度找初值
2.1 为什么选集中质量法而不是有限元：自由度少，控制好算
船舶推进轴系属于典型的细长旋转结构，从主机曲轴到螺旋桨，长度十几米甚至几十米，直径往往只有二三百毫米。如果整体用实体有限元，网格量不小，而且边界条件里轴承油膜刚度、螺旋桨附水质量都是估计值，网格再细也不一定准。工程上做扭振分析，最常用的还是集中质量模型：把每一根轴段折算成刚性圆盘，只在圆盘上保留转动惯量，轴段只保留扭转刚度；做回旋振动时，圆盘上保留质量与回转效应，轴段保留弯曲刚度。这样一条三十段左右的轴，系统自由度只有三十几，MATLAB里用普通eig也能秒出结果。
更重要的是，后续做振动控制分析时，需要把运动方程化成状态空间模型，控制器设计基于线性定常系统。集中质量模型天然就是低阶线性系统，方便直接转为A、B、C、D矩阵。而有限元模型动辄上万自由度，做主动控制前必须先做模态截断和降阶，又多一层不确定。所以我一般先用集中质量模型把轴系的固有特性算清楚，再决定要不要用有限元核对局部细节。
2.2 核心建模代码：轴段转动惯量与扭转刚度的MATLAB实现
集中质量建模第一步，是把轴系图纸拆成若干轴段，每段参数包括内径、外径、长度、材料密度、剪切模量。下面这段代码是我常用的输入方式，每一行对应一个轴段，数组shafts按从主机端到螺旋桨端排列：
% 轴段参数表：每行 [内径, 外径, 长度, 材料密度, 弹性模量, 剪切模量] shafts = [ 0.050 0.180 0.450 7850 2.06e11 0.79e11; % 主机输出端短轴 0.050 0.200 1.200 7850 2.06e11 0.79e11; % 中间轴1 0.040 0.220 2.000 7850 2.06e11 0.79e11; % 中间轴2 0.038 0.240 2.500 7850 2.06e11 0.79e11; % 艉轴 ]; rho = shafts(:,4); G = shafts(:,6); d_in = shafts(:,1); d_out = shafts(:,2); L = shafts(:,3); % 截面几何极惯性矩和圆盘转动惯量 I_p_geom = pi/32 .* (d_out.^4 - d_in.^4); % m^4 J_disk = rho .* L .* I_p_geom; % kg m^2 K_t = G .* I_p_geom ./ L; % N m/rad这段代码核心是三个量：I_p_geom是圆环截面极惯性矩，单位是m^4；J_disk是这一轴段折算到节点处的转动惯量，单位kg·m^2；K_t是轴段扭转刚度，单位N·m/rad。注意J_disk只是把整个轴段当作一个圆盘时的值，实际建模时应把相邻两个半轴段的惯量合并到共享节点上，而不是每段单独当成一个独立圆盘，否则会多算一倍惯量。
参数说明：内径填0就是实心轴；如果是空心轴，内径不能取0。G即剪切模量，钢轴通常取7.9e10 Pa，但实际会随材料牌号有差异。如果轴段有法兰、联轴器，它们的转动惯量要额外加在对应节点上，比如飞轮、齿轮箱连接法兰，通常厂家会给惯量值，直接加进去即可。
2.3 轴承刚度和螺旋桨附水质量的给出方法
轴系振动分析里，轴承不能当成刚体支撑。滑动轴承的油膜刚度与转速、载荷、间隙、油温都有关系。扭振分析中轴承一般不起作用，因为轴在轴承处仍可自由转动；但纵向振动和回旋振动分析必须在对应节点加上支撑刚度。程序包里常见的做法是预留一个轴承刚度数组，例如艉轴承取5e7~2e8 N/m，中间轴承取1e8~5e8 N/m，数值需要根据载荷和转速用雷诺方程估算，或者直接查轴承厂家的刚度曲线。
螺旋桨侧最容易被低估的，是附水质量和附水惯量。转动的螺旋桨带动周围水一起运动，等效于自身质量和转动惯量增加。对扭转振动，螺旋桨的附水转动惯量一般取桨叶自身转动惯量的1.3~1.6倍；对纵向振动，附水质量约为螺旋桨排水质量的0.2~0.5倍。这个系数是半经验值，不同资料差很多，所以我在程序里会把它单独列成一个变量prop_added_inertia_ratio，方便后面做敏感性分析。
另外，阻尼矩阵不要瞎给。扭振分析里结构阻尼通常体现为模态阻尼比，每阶取0.01~0.05。如果程序里用的是物理坐标下的质量、刚度、阻尼矩阵，则需要给定比例阻尼系数α和β，满足C = αM + βK。α和β可以通过两阶已知阻尼比反推，例如：
% 已知第1、2阶固有频率w1 w2和阻尼比zeta1 zeta2 % alpha = 2*w1*w2*(zeta1*w2 - zeta2*w1)/(w2^2 - w1^2) % beta = 2*(zeta2*w2 - zeta1*w1)/(w2^2 - w1^2)这两个系数算出来后，直接组装C = alpha*M + beta*K。这一套模型就齐了。
3. 扭转振动固有频率计算：传递矩阵扫频与临界转速判断
3.1 传递矩阵法的状态向量约定：转角与扭矩的传递关系
有了节点转动惯量和轴段刚度，求固有频率最直接的方法是组装质量矩阵和刚度矩阵，再用eig(K, M)解广义特征值问题。但船用轴系里用得更普遍的是传递矩阵法，因为它只需在两端边界条件之间来回传递，避免了大矩阵的特征值分解，程序逻辑也更贴近“扫频找共振”的物理图像。
传递矩阵法中最容易落地的是Holzer递推格式：从自由端出发，假设起始转角为1、起始扭矩为0，逐个经过圆盘和轴段，递推出末端扭矩。末端扭矩等于0时，当前频率就是固有频率。圆盘的作用是让扭矩改变-J*w^2*theta，转角不变；轴段的作用是让转角改变T/k，扭矩不变。这里的theta是扭转角，T是扭矩，w是角频率，J是圆盘转动惯量，k是轴段扭转刚度。
这个方法的巧妙之处在于边界条件自己“冒”出来。两端自由的轴系，起始端扭矩是0，末端扭矩也必须是0；只有频率恰为固有频率时，激励频率与惯性力、弹性力恰好平衡，末端才能满足零扭矩条件。其余频率下末端扭矩都不为零，残差的符号和大小就指示了离固有频率多近。
3.2 扫频找根的MATLAB代码：过零点加二分收敛
我常用的扫频代码如下，频率范围先给大，后面再细化：
% 已知节点惯量 iners（含飞轮、螺旋桨附水）和轴段刚度 ks % iners 长度比 ks 多 1 function fn = find_torsional_modes(iners, ks, fmax, N) f_scan = linspace(0.01, fmax, N); res = zeros(size(f_scan)); for i = 1:N res(i) = torsional_residual(2*pi*f_scan(i), iners, ks); end % 找符号变化的位置 idx = find(diff(sign(res)) ~= 0); fn = zeros(size(idx)); for k = 1:length(idx) fn(k) = fzero(@(f) torsional_residual(2*pi*f, iners, ks), ... [f_scan(idx(k)), f_scan(idx(k)+1)]); end end function r = torsional_residual(w, iners, ks) theta = 1.0; % 起始转角 T = 0.0; % 起始扭矩（自由端） for i = 1:length(iners) % 经过圆盘：扭矩增加 -J*w^2*theta，转角不变 T = T - iners(i) * w^2 * theta; % 经过轴段：转角变化 T/k，扭矩不变 if i < length(iners) theta = theta + T / ks(i); end end r = T; % 末端扭矩，固有频率处应过零 end逻辑说明：主函数先做粗扫频，res记录每个候选频率下的末端扭矩。固有频率附近残差必然变号，所以用diff(sign(res))找到过零区间，再交给fzero精确求根。torsional_residual函数内部严格按“圆盘—轴段—圆盘—轴段”的顺序递推，起始自由端扭矩为0，最后一个圆盘之后不再经过轴段，直接得到末端扭矩。这个递推式也便于画振型：把每段递推时的theta记下来，归一化就是模态转角形状。
参数说明：fmax要根据发动机最高转速和考虑的最高激励阶次来定。比如最高转速1800 rpm，轴频30 Hz，考虑8阶激励就到240 Hz，fmax取300 Hz足够。N粗扫点数取1000~2000，太多会增加无谓耗时。fzero要求区间两端函数值异号，这个条件已经由diff(sign(res))保证；如果不放心，可以加大N或改为在过零点附近再加密一次。
3.3 从固有频率到坎贝尔图：判断临界转速是否落在工作区
算完固有频率，下一步是判断哪些转速下会发生共振。螺旋桨产生的脉动扭矩主阶次等于桨叶数，四叶桨就是四阶；柴油机还要考虑发火阶次，比如六缸四冲程机的主简谐是3阶、6阶、9阶等。把激励频率画成随转速变化的斜线，把固有频率画成水平线，相交点对应的转速就是临界转速。MATLAB里用plot就能画：
rpm = 120:10:2100; % 转速范围 freq_h = rpm/60 * 4; % 四叶桨激励 figure; hold on; for i = 1:length(fn) yline(fn(i), '--', sprintf('固有频率 %.2f Hz', fn(i))); end plot(rpm, freq_h, 'r', 'LineWidth', 2); xlabel('曲轴转速 rpm'); ylabel('频率 Hz');运行这行脚本，就能直接看到哪几个转速落在工作转速范围里。如果交点在常用转速区，且激励频率与固有频率之间裕度不足5%，就得考虑加装减振器，或者调整轴径、改变节点惯量来移动固有频率。程序包里通常还附带一个函数，自动搜索交点并输出警告，我自己的习惯是宁可漏算高阶也要把所有低于1.2倍最高激励频率的固有频率都找出来，因为实际机桨匹配下高阶共振也不少见。
4. 振动控制仿真：从调频减振器到LQR主动控制的参数整定
4.1 开环响应仿真：施加螺旋桨脉动扭矩看轴段扭转角
固有频率算完，控制分析才有起点。控制之前先看开环响应：在螺旋桨节点上施加一个脉动扭矩激励，观察靠近主机端或飞轮端的扭转角幅值。用集中质量模型直接组装状态方程：
M_mtx = diag(iners); % 质量矩阵（转动惯量） n = length(iners); K_mtx = zeros(n, n); % 组装刚度矩阵：相邻节点用轴段刚度连接 for j = 1:length(ks) K_mtx(j, j) = K_mtx(j, j) + ks(j); K_mtx(j+1, j+1) = K_mtx(j+1, j+1) + ks(j); K_mtx(j, j+1) = K_mtx(j, j+1) - ks(j); K_mtx(j+1, j) = K_mtx(j+1, j) - ks(j); end C_mtx = 0.01*M_mtx + 0.001*K_mtx; % 比例阻尼示意 % 状态空间 A = [zeros(n), eye(n); -M_mtx\K_mtx, -M_mtx\C_mtx]; exc_vec = zeros(n,1); exc_vec(end) = 1; % 螺旋桨节点激励 B = [zeros(n,1); M_mtx\exc_vec]; C_out = zeros(1, 2*n); C_out(1) = 1; % 输出第一个节点转角 sys = ss(A, B, C_out, 0); t = 0:1e-4:2; f_exc = 30; % 轴频激励示例，实际要换激励阶次 u = sin(2*pi*f_exc*t).'; [y, t] = lsim(sys, u, t); plot(t, y);这里的状态是[转角; 角速度]，维度是2n。注意矩阵组装时K_mtx只取了简化的三对角结构，实际程序包里还要处理飞轮和螺旋桨节点的附加惯量。这里的阻尼系数是示意值，真实值要靠实验或经验取。
参数说明：exc_vec是激励位置向量，比如螺旋桨节点编号是n_prop，那么exc_vec(n_prop)=1。观察输出点若选中间轴某节点，就改C_out中对应位置为1。f_exc必须和转速匹配，比如想模拟1800 rpm时四叶桨激励，就是1800/60*4=120 Hz，不是30。上例里30 Hz是轴频，默认写了轴频，实际要按激励阶次改。
从开环响应能看到哪些转速下振幅放大明显。如果振幅超标，就要上控制。
4.2 被动控制与动力吸振器的参数设计
被动控制最常见的是在轴系上安装动力吸振器（调谐质量阻尼器），针对某一阶固有频率起作用。设计参数是吸振器固有频率f_d、质量比mu和阻尼比zeta_d。质量比通常取0.01~0.1，超过0.1体积重量就不现实。最优调谐比近似为1/(1+mu)，最优阻尼比近似为sqrt(3*mu/(8*(1+mu)^3))。程序里先把动力学方程扩展到含吸振器节点的自由度，再重新算闭环响应：
% 吸振器参数计算 mu = 0.05; % 质量比 main_inertia = sum(iners); % 主系统总惯量（简化） m_d = mu * main_inertia; % 吸振器惯量 kg m^2 f_target = 45; % 目标扭振固有频率 Hz f_d = f_target / (1 + mu); % 调谐频率 zeta_d = sqrt(3*mu/(8*(1+mu)^3)); k_d = (2*pi*f_d)^2 * m_d; % 吸振器刚度 N m/rad注意吸振器用于扭振时是“扭转吸振器”，本质上是一个在轴系节点上附加的惯性盘，通过弹性元件与轴系连接。它的惯量m_d要按kg·m^2计，刚度k_d是N·m/rad。把m_d和k_d加到对应节点上，重新组矩阵，扫频看固有频率分裂成两个。被动吸振器的缺点是只针对窄频带，转速变化大的轴系效果有限。
4.3 主动控制：状态空间建模与LQR权重调节
主动控制思路更直接：在轴系某节点（比如推力轴承处或中间轴段）施加额外主动扭矩，抵消螺旋桨激励。程序包里的控制部分通常基于状态空间模型设计LQR或H∞。LQR控制器是工程里最容易上手的一种，代价函数为J = ∫(x^T Q x + u^T R u) dt。Q矩阵决定状态收敛速度，R矩阵决定控制能耗。
MATLAB里用lqr一步就能得到反馈增益K：
% 状态空间矩阵沿用4.1的A、B Q = diag([ones(1,n)*10, ones(1,n)*1]); % 状态加权 R = 1e-4; % 控制力加权 K = lqr(A, B, Q, R); A_cl = A - B*K; sys_cl = ss(A_cl, B, C_out, 0); [y_cl, t] = lsim(sys_cl, u, t);逻辑说明：lqr函数需要满足A、B矩阵可控，且R>0。Q中对转角加权大，就会优先压转角幅值；对速度加权大，则优先抑制转速波动。R越小，控制力越强，但执行器可能饱和。实际工程中作动器输出扭矩有上限，必须先在程序里加饱和模块，再跑闭合响应，否则控制器在仿真里“完美”，实验台上一开机就抖。
另外，LQR只反馈全状态，现实中只有部分节点能放传感器，还要用lqr配合kalman设计LQG，或者用输出反馈。程序包如果只提供lqr，那只能验证理论上限，要落地必须做降维观测。我的做法是先跑全状态LQR得到性能上界，再从可测节点设计观测器，确认性能损失在可接受范围内。
5. 轴系MATLAB程序避坑指南：单位、惯量与数值溢出的5个坑
5.1 单位混用：N·m与N·mm让临界转速差一个数量级
现象：程序跑完，固有频率几十赫兹，但换算成临界转速后比试航报告的共振转速高出好几个数量级，完全对不上。
原因：输入轴段长度用了毫米，弹性模量用了MPa，导致刚度K的单位变成N·mm/rad，而转动惯量用的是kg·m^2。omega = sqrt(K/J)的根号里差了10^6量级，频率直接差出千倍。
解决：在程序入口统一转换为SI基本单位。比如用户习惯输入mm，就先除以1000。我在程序第一行写死：assert(max(L)<10, '长度请用m')，用数据范围检查单位。另外把所有扭矩、功率单位也一并规范，避免后续控制仿真里力与力矩混用。
5.2 转动惯量漏算：螺旋桨附水惯量不是可选项
现象：计算的第一阶扭振固有频率比实测高20%~40%，而且后几阶误差越来越大。
原因：只按轴段材料算了转动惯量，忘了把螺旋桨、飞轮、齿轮系的质量加上，也漏了螺旋桨推进时附水效应。轴系固有频率主要由大惯量部件决定，螺旋桨节点惯量占系统总惯量很大比例。
解决：螺旋桨自身转动惯量可以从桨几何计算，附水惯量系数取1.3~1.6，不确定时先取1.5，然后用第6章的方法做参数标定。飞轮、联轴器的惯量不要只用图纸外形估算，让设备厂家提供实测值最稳。
5.3 频率单位混用：Hz与rad/s让你找不到根
现象：扫频程序在f=100 Hz附近没有过零点，但把f改成628时又出现符号变化，判断固有频率出现在100 rad/s附近。
原因：传递矩阵点矩阵里用的是w^2，如果你直接用f去当w，相当于把扫描轴当成rad/s，最后得到的根单位是rad/s，输出时却标成Hz。
解决：代码里只保留一个角频率变量w = 2*pi*f，输入和输出都用w/2/pi转回Hz。写注释时把单位写清楚，别写w=omega=2πf又漏乘。
5.4 传递矩阵数值溢出：轴段多频率高时矩阵爆掉
现象：模型超过20个轴段，扫描到高频区时，残差曲线出现NaN，或者在某些频率突然中断。
原因：传递矩阵中1/k和J*w^2随频率增长会变得很大，连乘几十次后矩阵元素达到1e20以上，双精度浮点溢出。高频段固有频率本身也密集，数值溢出的假过零点容易误导。
解决：改用无量纲化传递矩阵：频率用w/w_ref归一化，刚度、惯量除以基准值；或者直接用eigs(K, M, n_required, 'smallestabs')解广义特征值问题，MATLAB自带的特征值求解器数值稳定性更好。程序包里如果目标是算低阶频率，我们只需要前5阶，没必要把扫描范围设到1000 Hz以上。
5.5 仿真步长设错：控制响应波形出现伪振荡
现象：主动控制闭环仿真波形毛刺多，明明控制器逻辑没问题，却看到高频抖动叠加在主响应上。
原因：lsim的输出步长或状态空间离散化步长太大，系统固有高频模态（可能是建模带来的数值模态）发生混叠，产生虚假振荡。步长太小又会让仿真时间成倍增加，控制迭代变慢。
解决：先算系统所有特征值，取最高关心频率的1/10作为仿真步长。比如系统最高固有频率200 Hz，步长取0.0005s甚至更小。用ss连续模型lsim时指定固定步长t = 0:dt:T，不要用自动步长，避免在刚性区域步长突变。另外检查A_cl的稳定性，如果出现正实部特征值，先修模型再谈步长。
6. 进阶技巧：用实测数据反调模型参数，把程序变成标定工具
程序算出的结果永远是模型的输出，模型参数不对，结果再精确也只是白算。我拿到一个新轴系程序，第一步不是直接跑设计工况，而是找一条同类型船舶的实测扭振数据，比如试航时测的轴段扭转角或轴表面应变，把固有频率实测值和计算结果对比。对不上时，优先怀疑螺旋桨附水系数和艉轴承油膜刚度，这两个参数最难直接测量。
具体做法是把模型固有频率计算函数封装成一个输入参数、输出频率向量的函数，然后用fminsearch去拟合实测频率：
fn_model = @(par) compute_frequencies(par(1), par(2)); % 附水系数, 轴承刚度 fun = @(par) sum((fn_model(par) - fn_meas).^2); par_opt = fminsearch(fun, [1.5, 1e8]);注意fn_meas至少要有两阶频率，否则参数辨识不唯一。调完后务必重新画坎贝尔图，确认临界转速位置也移动了。这一步做完，程序就从“演示用”变成了“可以信的工具”。
我的教训是第一次做这类程序时，只用了厂家给的轴段尺寸，漏掉了螺旋桨附水，结果前两阶固有频率差出12%。后来靠实测数据反调，把附水系数从1.0改成1.48、艉轴承刚度从5e8调到8e7，误差才压到1.2%。所以，程序包里最好预留参数标定接口，而不是把数据硬编码在脚本开头。希望帮到你。
本文还有配套的精品资源，点击获取
