简介这是面向铁路声学与结构动力学方向的MATLAB源码包聚焦轨道频散曲线的计算与物理意义分析。源码采用传递矩阵法将钢轨分别视为欧拉梁与Timoshenko梁通过构建周期单元间的传递矩阵计算不同频率下波的传播速度并绘制频散曲线可帮助研究者理解声子晶体的带隙特征、横波传播规律及频散特性适用于周期性轨道结构振动与噪声控制相关课题。压缩包共2个文件均为.m脚本大小约2KB其中oula_guangyi.m处理欧拉梁模型timo_guangyi.m处理Timoshenko梁模型代码包含参数设置、矩阵构建、迭代求解与绘图模块便于直接运行或二次开发。该资源已有535人学习适合具备MATLAB基础和结构波动理论背景的工程师、研究生快速上手借助源码验证理论推导并开展参数化分析。1. 轨道频散曲线为什么钢轨里的波速不是一个定值同样一根 60 kg/m 的钢轨在低频段弯曲波相速度可能不到 1000 m/s而在 5 kHz 以上却可能超过 2000 m/s如果换用纵向波速度又会稳定在 5000 m/s 上下。这里的关键不是“钢轨材料变了”而是波在钢轨这个受限波导中传播时频率直接决定了波速和波形形态。把“波在不同频率下跑多快、能传多远、以什么模态传”画成曲线就是轨道频散曲线。对做钢轨探伤、轮轨噪声分析、传感器布点和结构健康监测的工程师来说读不懂这张图就很难解释“为什么同一个激励信号在钢轨里传播 50 m 后波形被拉得像一串尾巴”。下面从频散曲线的物理意义讲起再给出用 MATLAB 计算轨道频散曲线的完整脚本最后说清楚它在检测频率选择和群速度补偿里的实际用法。2. 频散曲线的物理意义把相速度、群速度和截止频率一次讲透2.1 相速度与群速度一个在“找相位”一个在“搬能量”先给两个基本定义。相速度描述某个等相位面沿钢轨推进的速度表达式为 c_p ω/k群速度描述整个波包的能量传播速度表达式为 c_g dω/dk。在非频散介质中两者相同在钢轨这种受限波导中两者明显偏离。工程上最容易踩的误区是把检测系统标定的“波速”直接当成缺陷定位里的距离换算率。实际传感器记录到的波包前沿由群速度决定而波包内部的相位则按相速度移动。在 5 kHz 附近弯曲波的相速度与群速度可能相差 300600 m/s如果标距按相速度计算波包实际到达时间就会出现明显偏差。更麻烦的是不同频率成分的波包速度各不相同叠加之后波形在时域上被逐渐拉长。这就是“频散”最直观的工程后果。2.2 频散从哪里来波动方程里的一个指数代入频散的来源不是材料非线性而是波导截面的几何约束。把钢轨垂向弯曲的位移设为 w(x,t) W·exp[i(kx−ωt)]代入考虑剪切变形的 Timoshenko 梁方程组后可以得到频散方程k⁴ − (ρ/E ρ/(κG))k²ω² − ρA/(EI)·ω² ρ²/(EκG)·ω⁴ 0其中 k 是波数ω 是角频率E 是弹性模量G 是剪切模量ρ 是密度A 是钢轨截面积I 是截面惯性矩κ 是剪切修正系数。这个方程把 k 和 ω 绑定在一起给定频率不能直接写出一个常数速度而要先求出对应的波数 k再通过 ω/k 得到相速度。方程中的每一项都有明确物理角色第一项来自弯曲刚度主导低频段行为第二项来自转动惯量与剪切变形的联合修正负责在高频段把速度增长趋势“拉平”第三项是惯性项决定波的低频起始状态第四项是剪切惯性耦合项是第二条传播分支出现的根源。把这些项逐频点求解并整理成曲线就是后面用 MATLAB 实现的核心内容。2.3 钢轨里的三类主导模态与模态截止钢轨不是简单一维梁而是一个开放的厚壁波导。实际传播的导波按质点运动形态分成三大类纵向波、弯曲波、扭转波。三类模态的速度范围和频散强弱差异很大在检测中必须区分对待。模态质点运动方向典型速度频散程度常见激发源纵向 L沿钢轨轴线约 5100 m/s弱车轴冲击、轨头缺陷弯曲 B垂直于轨面约 10003000 m/s强轮轨垂向力扭转 T绕纵向轴旋转约 2800 m/s中等车轮偏磨、扣件不平顺纵向波在低频段近似不频散传播远、波速稳定所以长距离探伤通常会优先考虑纵向激励弯曲波频散最严重但轮轨垂向力最容易激发出它因此现场信号里往往既包含有用的弯曲波也夹杂大量干扰成分。每种模态还存在一个“截止频率”概念某个高阶模态能被激发的最低频率。低于截止频率时该模态的波数出现虚部波沿钢轨呈指数衰减不再形成有效的传播波。这也是为什么在同一传感器布置下当激励频带升高时时域波形中的振型数量会突然增加。2.4 物理意义落地频散对检测“看到什么”的决定性影响把频散曲线翻译成工程语言曲线上的每一个点都代表“这个频率下的波包在钢轨中以什么速度前进”多条曲线在同一个频率重叠时代表同一频率同时存在多个波速的传播模态。检测系统在时域里收到的波形是多个频率、多个波速分量叠加后的结果。因此频散曲线的物理意义不能只停留在“波速随频率变化”这个层面。它直接决定距离分辨率、最大可检测距离和信号的可解释性。频散越强波形在传播过程中被拉伸得越明显两个距离很近的缺陷反射波在时域上就越难分开频散曲线越平坦的频率区间才是传感器布点和激励频带设计的最优区间。这一点在后面的群速度补偿中还会再次用到。3. 用 MATLAB 计算轨道频散曲线从 Timoshenko 梁模型起步3.1 为什么先选 Timoshenko 梁而不是 Euler–Bernoulli 梁Euler–Bernoulli 梁只考虑弯曲变形得到的频散关系是 c_p (EI/ρA)^(1/4)·√ω频率越高速度无限上升。这个趋势在钢轨低频段大体可用但超过 12 kHz 就明显失真钢轨轨头与轨底之间的剪切变形比普通梁更突出转动惯量也不能忽略。Timoshenko 梁在弯曲方程中同时引入剪切修正系数 κ 和转动惯量 ρI频散方程里因此多出 k²ω² 与 ω⁴ 两项使相速度在高频段收敛到接近瑞利波速的有限值。换来的代价是方程从关于 k 的四次方程变成关于 k² 的二次方程这正好适合用 MATLAB 的roots函数逐频点求根。对只需要轨道弯曲波频散曲线的场景Timoshenko 模型是精度和实现成本之间的一个平衡点。3.2 钢轨截面参数与剪切修正系数取值以 60 kg/m 钢轨为基准垂向弯曲计算常用的参数如下。特别强调截面惯性矩 I 应取绕水平中性轴的垂向弯曲惯性矩而不是扭转惯性矩用错数值低频段相速度会整体偏差约 10% 以上。参数含义取值单位E弹性模量2.1×10¹¹PaG剪切模量8.1×10¹⁰Paρ密度7850kg/m³A截面面积7.686×10⁻³m²I垂向弯曲惯性矩2.16×10⁻⁵m⁴κ剪切修正系数0.400.50——剪切修正系数 κ 是这条简化路径里唯一需要经验的参数。对轨头—轨底形状明显不对称的钢轨截面0.45 是常用的折中值如果只关心低频段取 0.50 对结果影响小于 3%。但频率到了 6 kHz 以上κ 从 0.40 变到 0.50第二支传播分支的起振频率可能移动数百赫兹处理高频问题时需要把 κ 当作一个可调参数来校对。3.3 频散方程的数值化把四阶方程压成 k² 的二次式将前面的频散方程改写为关于 x k² 的二次方程x² − (ρ/E ρ/(κG))ω²·x − ρA/(EI)·ω² ρ²/(EκG)·ω⁴ 0对每个给定频率求解两个根 x₁ 和 x₂。实正根对应传播波数负实根对应衰减波复数根说明该频率下没有有效传播波。在 MATLAB 中roots直接返回两个根筛选条件是“虚部接近 0 且实部大于 0”最后开方得到波数 k。需要留意浮点误差理论上为 0 的虚部在 MATLAB 里可能残留 1e-10 量级的噪声所以判断时用abs(imag(r)) 1e-8比直接比较imag(r) 0更稳妥。3.4 可运行的 MATLAB 脚本与相速度/群速度输出下面这段脚本可直接复制到 MATLAB R2016 以上版本运行。它扫描 20 Hz10 kHz输出两条弯曲传播分支的相速度与群速度。% rail_dispersion_timoshenko.m % 用 Timoshenko 梁模型计算 60 kg/m 钢轨弯曲波频散曲线 clear; clc; close all; % ---- 参数 ---- E 2.1e11; % 弹性模量 [Pa] G 8.1e10; % 剪切模量 [Pa] rho 7850; % 密度 [kg/m^3] A 7.686e-3; % 截面积 [m^2] I 2.16e-5; % 垂向弯曲惯性矩 [m^4] kappa 0.45; % 剪切修正系数 % ---- 频率扫描 ---- f linspace(20, 10000, 800); w 2 * pi * f; % ---- 频散方程组合系数 ---- p1 rho/E rho/(kappa*G); p2 rho*A / (E*I); p3 rho^2 / (E*kappa*G); % 两条传播分支的波数 k1 zeros(size(f)); k2 zeros(size(f)); for i 1 : numel(f) a -p1 * w(i)^2; % 一次项系数 b -p2 * w(i)^2 p3 * w(i)^4; % 常数项 r roots([1, a, b]); % 解 x k^2 % 只取正的实根对应传播波 r r(abs(imag(r)) 1e-8 real(r) 0); r sort(real(r)); if ~isempty(r) k1(i) sqrt(r(1)); end if numel(r) 1 k2(i) sqrt(r(2)); end end代码执行逻辑是先计算三个组合系数 p1、p2、p3再对 800 个频率点逐个构造二次方程并求根。abs(imag(r)) 1e-8用于剔除浮点虚部噪声real(r) 0保证只保留传播波对应的正波数。第二支在没有正根时保持为 0绘图时不会显示。接下来把波数换算成相速度和群速度% ---- 相速度与群速度 ---- idx1 k1 0; kp1 k1(idx1); wp1 w(idx1); fp1 f(idx1); cp1 wp1 ./ kp1; % 第1支相速度 % 群速度中心差分: cg dw/dk cg1 zeros(size(kp1)); cg1(2:end-1) (wp1(3:end) - wp1(1:end-2)) ./ ... (kp1(3:end) - kp1(1:end-2)); % 第2支 idx2 k2 0; kp2 k2(idx2); wp2 w(idx2); fp2 f(idx2); cp2 wp2 ./ kp2; cg2 zeros(size(kp2)); if numel(kp2) 2 cg2(2:end-1) (wp2(3:end) - wp2(1:end-2)) ./ ... (kp2(3:end) - kp2(1:end-2)); end % ---- 画图 ---- figure(Color, w, Position, [100 100 900 380]); subplot(1,2,1); hold on; grid on; box on; plot(fp1/1000, cp1/1000, LineWidth, 1.6); plot(fp2/1000, cp2/1000, --, LineWidth, 1.4); xlabel(频率 [kHz]); ylabel(相速度 [km/s]); legend(第1支, 第2支, Location, best); subplot(1,2,2); hold on; grid on; box on; plot(fp1(2:end-1)/1000, cg1(2:end-1)/1000, LineWidth, 1.6); plot(fp2(2:end-1)/1000, cg2(2:end-1)/1000, --, LineWidth, 1.4); xlabel(频率 [kHz]); ylabel(群速度 [km/s]); legend(第1支, 第2支, Location, best);这里没有使用 MATLAB 自带的gradient因为频率均匀分布但波数并不均匀gradient隐含的等间距假设会引入额外误差。中心差分用三个相邻点计算斜率两端点保持 0 只用于绘图。如果后续要做频散补偿建议把[fp1, cg1]作为速度表插值到 DFT 频率格点。3.5 结果怎么读两条弯曲分支各自代表什么脚本跑完后会看到两组关键结果。第一支弯曲波在 1 kHz 以下近似按 √f 上升这与 Euler–Bernoulli 梁的预期一致5 kHz 以后增长速度变缓相速度逐渐接近约 3000 m/s 量级的钢轨表面波速度。第二支在约 6.5 kHz 附近才出现对应剪切惯性耦合项的阈值频率该值可由 f_t 1/(2π)·√(κGA/(ρI)) 估计代入表内参数约为 6470 Hz。第二支代表钢轨截面剪切变形参与形成的高阶弯曲传播波。第二支出现后同一个频率下存在两个不同的群速度。也就是说一个宽频激励经过钢轨传播后时域波形中会出现两组速度不同的波包。这不一定是故障信号而是钢轨波导的固有传播特征通常叫做模式分叉。识别模式分叉是后续选择激励频带、解释波形时绕不开的一步。4. 从简化梁到真实钢轨截面用 SAFE 方法补上高频段的缺口4.1 梁模型在 8 kHz 以上还剩下哪些坑Timoshenko 梁能解释前几 kHz 的弯曲行为但超过 8 kHz 后钢轨横截面不再保持刚性平截面轨头、轨底会发生局部压扁、剪切畸变等截面变形。这些振型在梁模型里根本不存在却在高频轮轨噪声和扣件区域振动中占据主导。另一个常被忽略的因素是钢轨底部与轨枕扣件的约束它会改变局部边界条件使截面模态提前出现。实测中这类截面模态在 515 kHz 区间相当密集每张频散图上可能叠着几十条曲线。继续用单根梁模型算到这里误差已经超过工程可接受范围。要往更高频走常见做法是改用半解析有限元方法即 SAFE。4.2 SAFE 方法的基本思路横向有限元加轴向解析波SAFE 的核心思想是做方向分治钢轨横截面使用二维有限元离散长度方向保留解析简谐波形式 exp[i(kx−ωt)]。这样就把三维波动问题压缩到截面上的二维自由度计算规模远小于全三维有限元。离散后得到控制方程[K₀ kK₁ k²K₂ − ω²M] U 0其中 K₀ 对应面内应变刚度K₁、K₂ 分别来自轴向应变与转动耦合M 是质量矩阵U 是截面上所有节点自由度。给定波数 k 时这是关于 ω² 的广义特征值问题直接用eig(Km, M)求解给定频率求波数则需要解二次特征值问题相对更复杂。提示用 SAFE 前先想清楚目标频率。只需 10 kHz 以下且只关心弯曲波Timoshenko 梁足够10 kHz 以上或需要纵向、扭转、截面模态的完整频谱再上 SAFE。SAFE 的调试成本主要在网格划分和模态排序不适合只做一次粗略验证的场景。4.3 MATLAB 中拼接频散方程的一般套路下面给出已经验证过求解结构的骨架代码省略单元组装细节。最省事的建模方式是先用 MATLAB PDE 工具箱对钢轨截面做三角剖分再根据每个单元的节点坐标组装出 K₀、K₁、K₂、M 四组矩阵。% safe_dispersion_skeleton.m % 假设截面网格与节点自由度已经生成, 矩阵为 K0, K1, K2, M k_scan linspace(0.1, 30, 300); % 波数扫描 [1/m] n_modes_wanted 6; % 只关注前6阶模态 f_safe zeros(numel(k_scan), n_modes_wanted); for i 1 : numel(k_scan) k k_scan(i); Km K0 k*K1 k^2*K2; % 当前波数下的刚度矩阵 [~, D] eig(Km, M); % 广义特征值问题 omega2 real(diag(D)); omega2 omega2(omega2 0); % 丢弃伪零/负频率 omega sort(sqrt(omega2)); f_safe(i, :) omega(1:n_modes_wanted) / (2*pi); end plot(k_scan, f_safe / 1000, .); xlabel(波数 [1/m]); ylabel(频率 [kHz]);SAFE 扫描推荐扫 k 而不是扫 f原因有两点一是 k 给定时特征值问题保持线性对称结构调用eig稳定二是频散曲线通常呈现为频率随波数单调增长的形式从低频到高频逐条读取更方便。当矩阵规模较大时可以改用eigs(Km, M, n_modes_wanted, smallestabs)来减少计算量。4.4 与 Timoshenko 结果对比时的注意点拿到 SAFE 结果之后第一件事不是看高频而是验证低频段的弯曲支是否与 Timoshenko 曲线重合。如果 SAFE 的弯曲支在 2 kHz 以下与梁模型偏差超过 5%多半是网格密度不够或者边界条件处理不当。钢轨底部在没有扣件建模时应使用自由边界近似与真实扣件约束的差异在低频段影响较小。对比项Timoshenko 梁SAFE 截面模型自由度梁挠度 截面转角截面节点二维位移弯曲支覆盖频率约 08 kHz0几十 kHz是否能出纵向/扭转模态不能可以主要误差来源截面刚性假设、剪切修正网格密度、边界约束MATLAB 计算代价秒级分钟级取决于网格第二个注意点是模态排序。SAFE 在同一波数下会解出纵向、弯曲、扭转、截面模态的混合频谱eig不会按物理意义自动排序。常见处理方式是在后处理中按位移场方向余弦分类纵向分量占优的归为 L 支面内横向分量占优的归为 B 或 T 支。这个步骤比较繁琐最简单的方法是直接看振型图手动确认模态转换点附近的曲线归属。5. 把频散曲线用起来检测频率选取、群速度补偿与激励设计5.1 检测距离和频率的权衡别只看相速度很多人拿到频散曲线后习惯只看相速度然后选一个“波速最快”的频率工作。但检测距离和分辨率真正争夺的其实是群速度曲线。当群速度曲线在某段频率内保持平坦时波包形状在传播过程中基本稳定检测系统可以使用简单的阈值或互相关定位。对 60 kg/m 钢轨弯曲波这个平坦区间通常落在 13 kHz群速度变化幅度约在 3% 以内一旦超过 5 kHz群速度下降速度变快同样的传播距离会让波包多展宽 10%20%。在安排传感器间距时可以把传播距离控制在激励主频对应群速度的 30 倍波长以内。例如 2 kHz、群速度约 3900 m/s 时波长约 1.95 m30 倍波长约 58 m这个范围内反射信号的首波仍可分辩超过这个距离再谈高分辨率就需要引入补偿算法。5.2 用群速度做频散补偿把“拉长”的波形还原如果必须在频散较强的频段工作可以用群速度曲线对采集信号做一次频域相位补偿。原理是每个频率分量在钢轨中传播距离 L 后额外滞后 τ(f)L/c_g(f)补偿就是把这个滞后在频域中移回去。% compensate_dispersion.m % x: 采集时域信号; L: 探头到目标距离; f: FFT频率向量 X fft(x); cgc interp1(f_disp, cg_disp, f, linear, extrap); Xc X .* exp(1i * (2*pi * f .* L ./ cgc)); xc real(ifft(Xc));f_disp和cg_disp来自前面频散计算得到的群速度表interp1把离散群速度曲线插值到 DFT 频率格点。补偿后波形中原本被拉长的波包会重新聚拢但注意不要把这个补偿加到信噪比不足的频段否则只会放大噪声。实际工程处理时可以先对信号做带通滤波只保留频散曲线中线性度较好的频率范围再做相位补偿。5.3 激励波形设计的两个实用参数激励设计要同时控制中心频率和带宽最简单常用的是汉宁窗调制的正弦脉冲s(t) sin(2πf_ct)·0.5(1−cos(2πf_ct/N))。其中周期数 N 越小时域脉冲越窄、带宽越大N 越大频带越窄但时域越长。窄带脉冲传播距离更远距离分辨率差宽带脉冲分辨率高但更容易在频散区间被拉散需要配合上面的补偿代码使用。结合频散曲线选参数的常见做法是探测距离 30 m 内中心频率取 35 kHz、N 取 23得到约 1 kHz 的带宽探测距离 50 m 以上降到 1.52 kHz、N 取 35让主能量落在群速度曲线平坦区。这样既保证反射波有足够的峰值强度又不会让频散把波形破坏到无法识别。具体折中值取决于钢轨扣件间隔、传感器灵敏度和现场噪声底。本文还有配套的精品资源点击获取