简介超声Lamb波的频散曲线是无损检测中模式识别和参数优化的关键依据。这份Matlab程序包面向从事结构健康监测、超声检测研究的工程师与研究生基于Rayleigh2Lamb方程分析无限大板中的频散特性提供相速度与群速度的数值计算思路并通过编程绘制铝板中的频散曲线便于直接对照试验结果。压缩包共9个文件以8个m脚本为核心另有1个prj工程文件m脚本覆盖模态分类、扫频计算、曲线绘图等环节prj用于统一管理工程配置整体大小仅9KB轻量且易于二次开发。已有5325人学习浏览表明该方法在Lamb波检测教学与工程验证中获得较广泛认可。借助完整代码与清晰的文件分工读者可快速复现频散曲线开展模态分析、缺陷检测参数选择或扩展为其他板结构的求解程序。1. Lamb波频散曲线方程好抄根难求用MATLAB把Lamb波的频散曲线画出来真正卡人的不是绘图函数而是“稳定地求根”。频散曲线是相速度cp或群速度cg对频厚积F的曲线它直接决定导波检测里选什么频率、按什么波速算定位。做过结构健康监测或者薄板探伤的朋友都有印象低频段S0、A0两条线干干净净频厚积一过几个MHz·mm高阶模式挤在一起漏根、跳线、单位错乱轮着来。这篇文章面向正要搭导波仿真、或者想用频散曲线设计激励参数的研究生和工程师按“方程定义—求解实现—群速度与验证—避坑—延伸应用”的顺序把这个MATLAB实现讲透。2. 频散方程与参数约定先弄懂对称和反对称再谈画图2.1 S模式与A模式两组方程一套求解逻辑Lamb波在自由板中沿板面方向传播时按照位移关于板中面的对称性分成两类对称模式记作Sn反对称模式记作An。它们各自满足Rayleigh-Lamb方程形式很像但求解时谁也不能替代谁。以相速度cp、角频率ω、波数k ω/cp定义两个辅助量p² ω²(1/cL² - 1/cp²)q² ω²(1/cT² - 1/cp²)对称模式满足tan(qd) / tan(pd) -4k²pq / (q²-k²)²反对称模式满足tan(qd) / tan(pd) -(q²-k²)² / (4k²pq)这里d是半板厚。需要提醒的是当cp小于cT甚至小于cL时p或q会变成虚数tan作用在虚宗量上会退化为tanh。很多初版代码在这里直接报NaN不是方程写错了而是没有处理复变量。下方这张表把变量含义和这套资源里推荐的单位制列清楚后面所有公式和代码都按这个约定走。变量含义推荐单位cp相速度mm/us数值等于km/scL纵波波速mm/uscT横波波速mm/usf频率MHz与mm/us搭配h板总厚度mmF频厚积F f × hMHz·mmd半板厚d h/2mmmm/us这个单位制的妙处在于数值好看钢的cL约5.9、cT约3.2画图时纵轴直接标注km/s完全不需要换算。频厚积用MHz·mm数值上等于km/s量纲也方便和文献对照。2.2 材料参数换算从E、ν、ρ得到cL和cT多数人手里只有杨氏模量、泊松比和密度写代码第一步就是把它们换算成纵波和横波波速。公式如下cL sqrt(E·(1-ν) / (ρ·(1ν)·(1-2ν)))cT sqrt(E / (2ρ·(1ν)))下面这段MATLAB代码以普通钢板参数为例输出单位为m/s实际使用时换算成mm/us% 材料参数换算E[GPa], nu[-], rho[kg/m^3] E 210; % 钢的杨氏模量 nu 0.3; % 泊松比 rho 7800; % 密度 cL sqrt(E*1e9*(1-nu)/(rho*(1nu)*(1-2*nu))); % m/s cT sqrt(E*1e9/(2*rho*(1nu))); % m/s % 转成 mm/us方便后续与 MHz、mm 搭配 cL cL * 1e-3; cT cT * 1e-3; fprintf(cL %.3f mm/us, cT %.3f mm/us\n, cL, cT);这段代码的逻辑是先按国际单位制算出波速再做一次单位缩放。之所以要在最后统一转成mm/us是因为频厚积F f × h在MHz和mm组合下数值适中残差函数中的πF再乘上速度倒数后量级接近1不容易触发数值溢出。参数方面如果你换成铝或者复合材料只需要替换E、ν、ρ三行后续求解函数完全不用改。2.3 “半板厚”与频厚积你以为的厚度可能差了一倍Rayleigh-Lamb方程里的d是半板厚不是整板厚度。早期版本如果把总厚度h直接当成d代进去算出来的频散曲线会在频厚积轴上被压缩一半S1、A1的截止位置全都对不上文献。这个错误很隐蔽因为曲线形状看起来是合理的只有和参考图对比时才发现错位。更稳妥的做法是从一开始就用频厚积F f × h作为自变量把半板厚d从公式里消掉。因为pd πF·sqrt(1/cL² - 1/cp²)正好是一个只依赖F的表达式具体推导只需要把ω 2πf、d h/2代入pd的定义约分后F自然出现。后面所有求解代码都不再显式出现d只接收F、cL、cT三个参数这样既避开了半板厚陷阱也让求解器可以在不同板厚之间复用。3. MATLAB求解实现扫频、二分、去重、跟踪的完整流程3.1 核心残差函数用sin/cos形式避开tan奇点直接按tan形式写残差函数会遇到两类麻烦。第一tan在角度接近π/2时发散符号求根法会把奇点误判为过零点第二当cp小于cT时q变成虚数tan函数在复数域的行为不如sin/cos直观。因此把方程两边同时乘以cos(pd)cos(qd)改写成sin/cos组合形式对称模式残差sin(qd)cos(pd)(q²-k²)² 4k²pq·cos(qd)sin(pd) 0反对称模式残差4k²pq·sin(qd)cos(pd) (q²-k²)²·cos(qd)sin(pd) 0进一步把ω因子从p、q、k中提取出来残差函数只保留F和cp两个变量。这样写的好处是输入参数少且在F范围0.1到20 MHz·mm内不会因为频率量级产生数值问题。function res lambResidualNorm(cp, F, cL, cT, isSym) % lambResidualNorm 归一化残差 % cp : 相速度mm/us % F : 频厚积MHz*mmF 频率(MHz)*板厚(mm) % cL,cT: 纵波/横波波速mm/us % isSym: true对称模式false反对称模式 p sqrt((1/cL)^2 - (1/cp)^2); % 可能为虚数MATLAB复数运算自动处理 q sqrt((1/cT)^2 - (1/cp)^2); pd pi * F * p; qd pi * F * q; sp sin(pd); cp2 cos(pd); sq sin(qd); cq cos(qd); q2k2 (1/cT)^2 - 2*(1/cp)^2; % (q^2-k^2)/omega^2 coeff 4 * (1/cp)^2 * p * q; % 4*k^2*p*q/omega^4 if isSym res sq .* cp2 .* (q2k2.^2) coeff .* cq .* sp; else res coeff .* sq .* cp2 (q2k2.^2) .* cq .* sp; end res real(res); % 复数运算的虚部理论上抵消浮点残余直接去掉 end这段代码看起来短但两个设计点值得说明。第一p和q保持复数形式MATLAB的sin和cos原生支持复数参数当cp cT时自动退化为双曲函数不需要手工分情况讨论。第二残差函数只依赖F而不是频率和厚度分开传入这是上一节“消去半板厚”的落地实现。参数方面cp传入的是mm/us数值如果直接传入m/s残差函数会完全失效因为1/cp的量级差了1000倍。3.2 每个频厚积下的求根器对数网格扫一遍符号变化定区间残差函数准备好之后下一步就是在给定F下找出所有满足方程的cp。常见做法是在cp搜索域内密集采样找到残差符号变化的区间再对每个区间调fzero精确收敛。搜索域需要覆盖全部物理上可能的根。A0模式在低频时相速度趋近于0所以下界不能太高高阶模式在截止频率附近相速度趋于无穷上界也不能只设到cL。实践中常用范围是0.02cT到4cL。低频加一个对数采样对A0这种低相速度根更友好。function rootsCps solveLambAtF(F, cL, cT, isSym) % solveLambAtF 在给定频厚积下求解某类模式的所有相速度根 cpMin cT * 0.02; cpMax cL * 4; N 6000; % 对数网格低频段更密照顾A0的小相速度根 cpGrid logspace(log10(cpMin), log10(cpMax), N); r arrayfun((cp) lambResidualNorm(cp, F, cL, cT, isSym), cpGrid); % 符号变化的网格段即为疑似根区间 signProd sign(r(1:end-1)) .* sign(r(2:end)); idxSeg find(signProd 0); rootsCps zeros(1, numel(idxSeg)); for i 1:numel(idxSeg) a cpGrid(idxSeg(i)); b cpGrid(idxSeg(i) 1); try rootsCps(i) fzero((cp) lambResidualNorm(cp, F, cL, cT, isSym), [a b]); catch rootsCps(i) NaN; end end % 去掉求解失败的项并按相对容差合并离得过近的根 rootsCps rootsCps(~isnan(rootsCps)); if ~isempty(rootsCps) rootsCps uniquetol(rootsCps, 1e-5, DataScale, max(rootsCps)); end end这段代码的运行逻辑是先在整个物理区间内找符号变化点再用fzero把根精确出来。之所以要加try-catch是因为网格段端点处残差可能出现恰好为0或者接近0的情况fzero的闭区间要求两端严格异号遇到等号会报错。参数N6000对常见的01到20 MHz·mm区间足够如果扫到高阶模式密集区仍发现漏根把N提到12000耗时大约增加一倍但稳定性明显更好。3.3 模式追踪不要简单排序要靠连续性“认领”根每个频率点求出的根是一组无序数。如果直接把每个F点上的根从小到大排列后连线低频段还能对上高频段会出现S1和S2互换位置的情况画出来的曲线在交叉点附近会突然跳线。解决思路是对每个模式单独跟踪从上一个频率点的根出发在当前频率的所有根里找距离最近的那个认领为本模式的延续。已经被认领的根要从候选列表里剔除防止两个模式抢同一个根。function cpMat traceModes(Fvec, cL, cT, isSym) % traceModes 沿频率轴追踪每一条模式曲线 nF numel(Fvec); cpMat NaN(nF, 20); % 预留20列实际模式数自动扩展 for iF 1:nF roots solveLambAtF(Fvec(iF), cL, cT, isSym); if iF 1 % 第一个频率点按大小顺序初始化为S0,S1,... roots sort(roots); nMode numel(roots); cpMat(iF, 1:nMode) roots; else prev cpMat(iF-1, :); prev prev(~isnan(prev)); for im 1:numel(prev) if isempty(roots), break; end [~, jMin] min(abs(roots - prev(im))); cpMat(iF, im) roots(jMin); roots(jMin) NaN; % 防止两个模式认领同一个根 end end end cpMat cpMat(:, ~all(isnan(cpMat), 1)); end这段代码的关键在于“认领后置NaN”。如果不做这一步当两个模式的根距离相近时后面的模式可能把前面已经选走的根再选一次造成两个模式共用同一个解。参数方面列数20对常见钢板厚度已经够用复合材料或者厚板上的高阶模式超过20个时可以把这个数字改大。4. 群速度计算与可视化从相速度到实际传播速度4.1 数值差分前先平滑不然群速度曲线全是毛刺相速度曲线连续求出来后群速度不能直接抄公式里的解析表达式而是要用数值差分。群速度和相速度的关系是cg dω/dk其中k ω/cp。如果直接用相邻两个频率点做差分任何细微的相速度抖动量级会被放大画出来的曲线在高频段全是毛刺。% 假设 cpMat 维度为 nF x nMFvec 为对应频厚积列向量 Fvec Fvec(:); omega 2 * pi * Fvec; % 与厚度无关的比例关系画图时足够用 kMat omega ./ cpMat; % 波数矩阵 cgMat zeros(size(kMat)); for im 1:size(kMat, 2) kSm movmean(kMat(:, im), 11); % 先做11点滑动平均 cgMat(:, im) gradient(omega) ./ gradient(kSm); end % 两端diff噪声大裁掉前后5个点再画图 Fplot Fvec(6:end-5); cgPlot cgMat(6:end-5, :);逻辑上先平滑k再差分比先差分再平滑更稳。gradient在均匀网格上计算中心差分端点本身就不准所以裁头裁尾是必要操作。参数方面11点滑动窗口适合F步长小于0.05的密集扫频如果F步长较大窗口要相应缩小否则会把群速度曲线上真实存在的凹陷也抹平。4.2 绘图与坐标习惯横轴MHz·mm纵轴km/s频散曲线的标准画法是横轴频厚积F、纵轴速度。相速度和群速度可以画在同一张图里实线表示相速度、虚线表示群速度同一条模式用同一颜色这样一眼就能看出频散强弱。figure(Color, w); hold on; colors lines(6); linestyles {-, --}; for im 1:size(cpMat, 2) plot(Fvec, cpMat(:, im), Color, colors(im, :), LineWidth, 1.4); plot(Fplot, cgPlot(:, im), --, Color, colors(im, :), LineWidth, 1.2); end xlabel(频厚积 F [MHz\cdotmm]); ylabel(速度 v [km/s]); legend(S0 cp, S0 cg, A0 cp, A0 cg);因为mm/us和km/s数值相等纵轴直接标成km/s。读图时需要特别注意A0模式它的相速度在F趋近0时也趋近0群速度曲线却不是单调的会先升后降这是弯曲波典型的频散行为。S0模式低频段接近板波速度大约5.4 km/s然后逐渐下降。4.3 极限验证低频板速和高频瑞利速代码跑通后第一件事不是看全貌而是取两个极限点和公式对比。钢板例子中S0在F趋近0时应收敛到板波速度c_plate 2·cT·sqrt(1 - (cT/cL)²)高频时所有模式相速度和群速度都收敛到瑞利波速工程近似式为cR ≈ cT·(0.87 1.12·ν) / (1 ν)用前面钢的材料参数代入c_plate约5.44 km/scR约2.98 km/s。如果代码画出的S0低频值在5.4附近、高频所有曲线尾部收在3.0附近说明求解和单位制都没有问题。这两处验证价值很高因为它们是解析值不依赖任何其他软件。5. 避坑清单漏根、跳线、单位混乱、溢出和搜索域过窄5.1 漏根符号变化法在“接触零点”处失效现象某一阶模式曲线在某个频厚积区间突然缺了一小段把网格加密后又恢复了最初以为是fzero收敛失败实际上根压根没被检测到。原因符号变化判据依赖相邻采样点乘积为负。如果网格恰好落在根上或者两个根离得太近导致网格段内残差同号符号变化法就漏检了。解决两层保险。第一在网格采样后额外检查abs(r) 1e-8的网格点直接把它们也作为候选根第二对F方向步长加密让相邻F点的相速度变化不超过搜索网格间距保证漏掉的机会变小。我一般把F范围再分成两段扫低F段用细网格高F段用粗网格漏根概率能压到很低。5.2 模式跳线根不排序直接连线曲线互相串现象单独求每个F的根再画图S1和A1在某个频段发生明显交叉两条曲线像麻花一样拧在一起但对照文献图谱发现真实曲线并不会这样。原因求解器返回的根无序按大小排序后高频段不同模式相互逼近时顺序会发生交换连线就串了。解决用上一节traceModes里的“认领”逻辑每次以上一频率点的根为基准找最近根并且同一个根不允许被认领两次。这个方案比排序稳健得多代价只是每帧需要多跑一次min查找对几千个频率点来说耗时完全可以接受。5.3 单位混乱m/s配kHz算出满屏NaN现象所有设置看起来都对但fzero提示端点函数值不异号或者画出来的曲线数量级差了1000倍。检查残差函数里1/cp发现cp如果用的是m/s而cL、cT用的是mm/us速度倒数直接差出3个数量级。原因公式里所有速度必须同一单位制。mm/us加MHz是最舒服的组合但很多人习惯先把波速换成m/s频率用Hz这时如果没同步转换方程就废了。解决在脚本最开头做一次单位断言比如检查cL是否在1到20之间因为常见材料的纵波速度在mm/us制下最多是个位数。单位制统一后不妨再验证一次S0低频板速值这一步能排查掉90%以上的单位问题。5.4 溢出现象高频厚板时残差函数直接爆掉现象板厚达到几十毫米、频率到MHz时数组函数输出Inf或NaN符号变化判断完全失效求根器一个根都返回不了。原因sin/cos复数形式中虚部很大时指数增长double超过约1e308就溢出。归一化后的残差已经比tan形式好很多但极端频厚积下依然有风险。解决把残差整体除一个尺度因子cosh(imag(pd))·cosh(imag(qd))让sin/cos项变成有界组合。关键代码就一行% 在lambResidualNorm返回前加一句防溢出缩放 scale cosh(imag(pd)) .* cosh(imag(qd)); res res ./ scale;这样处理不会改变零点位置因为除以的是一个恒正的数。高频厚板场景里这条路比直接改N有效得多。5.5 搜索域太窄cp上限定到cL高次模式全丢现象S0、A0、S1、A1都在往上S2、A2、S3怎么也出不来或者到高频末尾突然断掉。原因高阶模式在截止频率附近相速度趋于无穷搜索域上限如果只到cL它们的根根本没有进入求解区间。解决把上界放宽到4倍cL以上计算量增加不多边界效果却很显著。低频段还要注意下界A0在F很小时相速度低于0.1cT下界用0.02cT在大多数情况下够用如果扫F低于0.05建议下界降到0.005cT。6. 验证与延伸把频散曲线落到激励设计上6.1 三种低成本验证手段逐点验证极限值是第一种手段。除了前面提到的板波速度和瑞利波速还可以看A0在F趋近0时相速度和群速度同时趋近0这是弯曲波低频行为。第二种手段是和文献里的标准谱图对比重点看S1群速度曲线的“凹陷”位置这个特征类似模式指纹不同材料参数下凹陷位置变化明显如果对得上代码基本没有大问题。第三种手段是切换材料参数比如换成铝观察S0低频板速是否等于对应公式值能检验代码对材料参数的依赖是否正确。6.2 从频散曲线到激励信号设计工程上拿到频散曲线后最常用的场景是选择激励频率。做法是看群速度曲线的平坦区平坦区意味着传播时波形包络畸变小。比如要激发S0模式就选S0群速度曲线在目标频段内梯度最小的点作为中心频率。另一个常见需求是估计两个模式到达时间差以便设计时间窗分离波包。到达时间差公式为Δt L·(1/cg1 - 1/cg2)其中L是传播距离cg1和cg2是两个模式的群速度。从曲线上读出对应F下的群速度值直接代入就能算出分离窗口。这套MATLAB实现里附带的扫频和模式追踪脚本可以直接改参数跑通不需要再从零搭框架。6.3 一点个人习惯以前我调一套复合材料的参数时把搜索域上限设成cL画出的曲线少了一整支高阶模式当时盯着屏幕怀疑是求解器写错了排查了两天才发现是边界问题。从那以后我每次换材料、换板厚都会先跑一遍极限验证确认S0低频板速和高频瑞利速都落在合理范围再去看整条曲线。二维曲线好看是表象求根稳定才是这套实现真正值钱的地方。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/15 13:30:57 · 阅读 1.2k 基于MATLAB的Lamb波频散曲线求解与模式追踪
简介超声Lamb波的频散曲线是无损检测中模式识别和参数优化的关键依据。这份Matlab程序包面向从事结构健康监测、超声检测研究的工程师与研究生基于Rayleigh2Lamb方程分析无限大板中的频散特性提供相速度与群速度的数值计算思路并通过编程绘制铝板中的频散曲线便于直接对照试验结果。压缩包共9个文件以8个m脚本为核心另有1个prj工程文件m脚本覆盖模态分类、扫频计算、曲线绘图等环节prj用于统一管理工程配置整体大小仅9KB轻量且易于二次开发。已有5325人学习浏览表明该方法在Lamb波检测教学与工程验证中获得较广泛认可。借助完整代码与清晰的文件分工读者可快速复现频散曲线开展模态分析、缺陷检测参数选择或扩展为其他板结构的求解程序。1. Lamb波频散曲线方程好抄根难求用MATLAB把Lamb波的频散曲线画出来真正卡人的不是绘图函数而是“稳定地求根”。频散曲线是相速度cp或群速度cg对频厚积F的曲线它直接决定导波检测里选什么频率、按什么波速算定位。做过结构健康监测或者薄板探伤的朋友都有印象低频段S0、A0两条线干干净净频厚积一过几个MHz·mm高阶模式挤在一起漏根、跳线、单位错乱轮着来。这篇文章面向正要搭导波仿真、或者想用频散曲线设计激励参数的研究生和工程师按“方程定义—求解实现—群速度与验证—避坑—延伸应用”的顺序把这个MATLAB实现讲透。2. 频散方程与参数约定先弄懂对称和反对称再谈画图2.1 S模式与A模式两组方程一套求解逻辑Lamb波在自由板中沿板面方向传播时按照位移关于板中面的对称性分成两类对称模式记作Sn反对称模式记作An。它们各自满足Rayleigh-Lamb方程形式很像但求解时谁也不能替代谁。以相速度cp、角频率ω、波数k ω/cp定义两个辅助量p² ω²(1/cL² - 1/cp²)q² ω²(1/cT² - 1/cp²)对称模式满足tan(qd) / tan(pd) -4k²pq / (q²-k²)²反对称模式满足tan(qd) / tan(pd) -(q²-k²)² / (4k²pq)这里d是半板厚。需要提醒的是当cp小于cT甚至小于cL时p或q会变成虚数tan作用在虚宗量上会退化为tanh。很多初版代码在这里直接报NaN不是方程写错了而是没有处理复变量。下方这张表把变量含义和这套资源里推荐的单位制列清楚后面所有公式和代码都按这个约定走。变量含义推荐单位cp相速度mm/us数值等于km/scL纵波波速mm/uscT横波波速mm/usf频率MHz与mm/us搭配h板总厚度mmF频厚积F f × hMHz·mmd半板厚d h/2mmmm/us这个单位制的妙处在于数值好看钢的cL约5.9、cT约3.2画图时纵轴直接标注km/s完全不需要换算。频厚积用MHz·mm数值上等于km/s量纲也方便和文献对照。2.2 材料参数换算从E、ν、ρ得到cL和cT多数人手里只有杨氏模量、泊松比和密度写代码第一步就是把它们换算成纵波和横波波速。公式如下cL sqrt(E·(1-ν) / (ρ·(1ν)·(1-2ν)))cT sqrt(E / (2ρ·(1ν)))下面这段MATLAB代码以普通钢板参数为例输出单位为m/s实际使用时换算成mm/us% 材料参数换算E[GPa], nu[-], rho[kg/m^3] E 210; % 钢的杨氏模量 nu 0.3; % 泊松比 rho 7800; % 密度 cL sqrt(E*1e9*(1-nu)/(rho*(1nu)*(1-2*nu))); % m/s cT sqrt(E*1e9/(2*rho*(1nu))); % m/s % 转成 mm/us方便后续与 MHz、mm 搭配 cL cL * 1e-3; cT cT * 1e-3; fprintf(cL %.3f mm/us, cT %.3f mm/us\n, cL, cT);这段代码的逻辑是先按国际单位制算出波速再做一次单位缩放。之所以要在最后统一转成mm/us是因为频厚积F f × h在MHz和mm组合下数值适中残差函数中的πF再乘上速度倒数后量级接近1不容易触发数值溢出。参数方面如果你换成铝或者复合材料只需要替换E、ν、ρ三行后续求解函数完全不用改。2.3 “半板厚”与频厚积你以为的厚度可能差了一倍Rayleigh-Lamb方程里的d是半板厚不是整板厚度。早期版本如果把总厚度h直接当成d代进去算出来的频散曲线会在频厚积轴上被压缩一半S1、A1的截止位置全都对不上文献。这个错误很隐蔽因为曲线形状看起来是合理的只有和参考图对比时才发现错位。更稳妥的做法是从一开始就用频厚积F f × h作为自变量把半板厚d从公式里消掉。因为pd πF·sqrt(1/cL² - 1/cp²)正好是一个只依赖F的表达式具体推导只需要把ω 2πf、d h/2代入pd的定义约分后F自然出现。后面所有求解代码都不再显式出现d只接收F、cL、cT三个参数这样既避开了半板厚陷阱也让求解器可以在不同板厚之间复用。3. MATLAB求解实现扫频、二分、去重、跟踪的完整流程3.1 核心残差函数用sin/cos形式避开tan奇点直接按tan形式写残差函数会遇到两类麻烦。第一tan在角度接近π/2时发散符号求根法会把奇点误判为过零点第二当cp小于cT时q变成虚数tan函数在复数域的行为不如sin/cos直观。因此把方程两边同时乘以cos(pd)cos(qd)改写成sin/cos组合形式对称模式残差sin(qd)cos(pd)(q²-k²)² 4k²pq·cos(qd)sin(pd) 0反对称模式残差4k²pq·sin(qd)cos(pd) (q²-k²)²·cos(qd)sin(pd) 0进一步把ω因子从p、q、k中提取出来残差函数只保留F和cp两个变量。这样写的好处是输入参数少且在F范围0.1到20 MHz·mm内不会因为频率量级产生数值问题。function res lambResidualNorm(cp, F, cL, cT, isSym) % lambResidualNorm 归一化残差 % cp : 相速度mm/us % F : 频厚积MHz*mmF 频率(MHz)*板厚(mm) % cL,cT: 纵波/横波波速mm/us % isSym: true对称模式false反对称模式 p sqrt((1/cL)^2 - (1/cp)^2); % 可能为虚数MATLAB复数运算自动处理 q sqrt((1/cT)^2 - (1/cp)^2); pd pi * F * p; qd pi * F * q; sp sin(pd); cp2 cos(pd); sq sin(qd); cq cos(qd); q2k2 (1/cT)^2 - 2*(1/cp)^2; % (q^2-k^2)/omega^2 coeff 4 * (1/cp)^2 * p * q; % 4*k^2*p*q/omega^4 if isSym res sq .* cp2 .* (q2k2.^2) coeff .* cq .* sp; else res coeff .* sq .* cp2 (q2k2.^2) .* cq .* sp; end res real(res); % 复数运算的虚部理论上抵消浮点残余直接去掉 end这段代码看起来短但两个设计点值得说明。第一p和q保持复数形式MATLAB的sin和cos原生支持复数参数当cp cT时自动退化为双曲函数不需要手工分情况讨论。第二残差函数只依赖F而不是频率和厚度分开传入这是上一节“消去半板厚”的落地实现。参数方面cp传入的是mm/us数值如果直接传入m/s残差函数会完全失效因为1/cp的量级差了1000倍。3.2 每个频厚积下的求根器对数网格扫一遍符号变化定区间残差函数准备好之后下一步就是在给定F下找出所有满足方程的cp。常见做法是在cp搜索域内密集采样找到残差符号变化的区间再对每个区间调fzero精确收敛。搜索域需要覆盖全部物理上可能的根。A0模式在低频时相速度趋近于0所以下界不能太高高阶模式在截止频率附近相速度趋于无穷上界也不能只设到cL。实践中常用范围是0.02cT到4cL。低频加一个对数采样对A0这种低相速度根更友好。function rootsCps solveLambAtF(F, cL, cT, isSym) % solveLambAtF 在给定频厚积下求解某类模式的所有相速度根 cpMin cT * 0.02; cpMax cL * 4; N 6000; % 对数网格低频段更密照顾A0的小相速度根 cpGrid logspace(log10(cpMin), log10(cpMax), N); r arrayfun((cp) lambResidualNorm(cp, F, cL, cT, isSym), cpGrid); % 符号变化的网格段即为疑似根区间 signProd sign(r(1:end-1)) .* sign(r(2:end)); idxSeg find(signProd 0); rootsCps zeros(1, numel(idxSeg)); for i 1:numel(idxSeg) a cpGrid(idxSeg(i)); b cpGrid(idxSeg(i) 1); try rootsCps(i) fzero((cp) lambResidualNorm(cp, F, cL, cT, isSym), [a b]); catch rootsCps(i) NaN; end end % 去掉求解失败的项并按相对容差合并离得过近的根 rootsCps rootsCps(~isnan(rootsCps)); if ~isempty(rootsCps) rootsCps uniquetol(rootsCps, 1e-5, DataScale, max(rootsCps)); end end这段代码的运行逻辑是先在整个物理区间内找符号变化点再用fzero把根精确出来。之所以要加try-catch是因为网格段端点处残差可能出现恰好为0或者接近0的情况fzero的闭区间要求两端严格异号遇到等号会报错。参数N6000对常见的01到20 MHz·mm区间足够如果扫到高阶模式密集区仍发现漏根把N提到12000耗时大约增加一倍但稳定性明显更好。3.3 模式追踪不要简单排序要靠连续性“认领”根每个频率点求出的根是一组无序数。如果直接把每个F点上的根从小到大排列后连线低频段还能对上高频段会出现S1和S2互换位置的情况画出来的曲线在交叉点附近会突然跳线。解决思路是对每个模式单独跟踪从上一个频率点的根出发在当前频率的所有根里找距离最近的那个认领为本模式的延续。已经被认领的根要从候选列表里剔除防止两个模式抢同一个根。function cpMat traceModes(Fvec, cL, cT, isSym) % traceModes 沿频率轴追踪每一条模式曲线 nF numel(Fvec); cpMat NaN(nF, 20); % 预留20列实际模式数自动扩展 for iF 1:nF roots solveLambAtF(Fvec(iF), cL, cT, isSym); if iF 1 % 第一个频率点按大小顺序初始化为S0,S1,... roots sort(roots); nMode numel(roots); cpMat(iF, 1:nMode) roots; else prev cpMat(iF-1, :); prev prev(~isnan(prev)); for im 1:numel(prev) if isempty(roots), break; end [~, jMin] min(abs(roots - prev(im))); cpMat(iF, im) roots(jMin); roots(jMin) NaN; % 防止两个模式认领同一个根 end end end cpMat cpMat(:, ~all(isnan(cpMat), 1)); end这段代码的关键在于“认领后置NaN”。如果不做这一步当两个模式的根距离相近时后面的模式可能把前面已经选走的根再选一次造成两个模式共用同一个解。参数方面列数20对常见钢板厚度已经够用复合材料或者厚板上的高阶模式超过20个时可以把这个数字改大。4. 群速度计算与可视化从相速度到实际传播速度4.1 数值差分前先平滑不然群速度曲线全是毛刺相速度曲线连续求出来后群速度不能直接抄公式里的解析表达式而是要用数值差分。群速度和相速度的关系是cg dω/dk其中k ω/cp。如果直接用相邻两个频率点做差分任何细微的相速度抖动量级会被放大画出来的曲线在高频段全是毛刺。% 假设 cpMat 维度为 nF x nMFvec 为对应频厚积列向量 Fvec Fvec(:); omega 2 * pi * Fvec; % 与厚度无关的比例关系画图时足够用 kMat omega ./ cpMat; % 波数矩阵 cgMat zeros(size(kMat)); for im 1:size(kMat, 2) kSm movmean(kMat(:, im), 11); % 先做11点滑动平均 cgMat(:, im) gradient(omega) ./ gradient(kSm); end % 两端diff噪声大裁掉前后5个点再画图 Fplot Fvec(6:end-5); cgPlot cgMat(6:end-5, :);逻辑上先平滑k再差分比先差分再平滑更稳。gradient在均匀网格上计算中心差分端点本身就不准所以裁头裁尾是必要操作。参数方面11点滑动窗口适合F步长小于0.05的密集扫频如果F步长较大窗口要相应缩小否则会把群速度曲线上真实存在的凹陷也抹平。4.2 绘图与坐标习惯横轴MHz·mm纵轴km/s频散曲线的标准画法是横轴频厚积F、纵轴速度。相速度和群速度可以画在同一张图里实线表示相速度、虚线表示群速度同一条模式用同一颜色这样一眼就能看出频散强弱。figure(Color, w); hold on; colors lines(6); linestyles {-, --}; for im 1:size(cpMat, 2) plot(Fvec, cpMat(:, im), Color, colors(im, :), LineWidth, 1.4); plot(Fplot, cgPlot(:, im), --, Color, colors(im, :), LineWidth, 1.2); end xlabel(频厚积 F [MHz\cdotmm]); ylabel(速度 v [km/s]); legend(S0 cp, S0 cg, A0 cp, A0 cg);因为mm/us和km/s数值相等纵轴直接标成km/s。读图时需要特别注意A0模式它的相速度在F趋近0时也趋近0群速度曲线却不是单调的会先升后降这是弯曲波典型的频散行为。S0模式低频段接近板波速度大约5.4 km/s然后逐渐下降。4.3 极限验证低频板速和高频瑞利速代码跑通后第一件事不是看全貌而是取两个极限点和公式对比。钢板例子中S0在F趋近0时应收敛到板波速度c_plate 2·cT·sqrt(1 - (cT/cL)²)高频时所有模式相速度和群速度都收敛到瑞利波速工程近似式为cR ≈ cT·(0.87 1.12·ν) / (1 ν)用前面钢的材料参数代入c_plate约5.44 km/scR约2.98 km/s。如果代码画出的S0低频值在5.4附近、高频所有曲线尾部收在3.0附近说明求解和单位制都没有问题。这两处验证价值很高因为它们是解析值不依赖任何其他软件。5. 避坑清单漏根、跳线、单位混乱、溢出和搜索域过窄5.1 漏根符号变化法在“接触零点”处失效现象某一阶模式曲线在某个频厚积区间突然缺了一小段把网格加密后又恢复了最初以为是fzero收敛失败实际上根压根没被检测到。原因符号变化判据依赖相邻采样点乘积为负。如果网格恰好落在根上或者两个根离得太近导致网格段内残差同号符号变化法就漏检了。解决两层保险。第一在网格采样后额外检查abs(r) 1e-8的网格点直接把它们也作为候选根第二对F方向步长加密让相邻F点的相速度变化不超过搜索网格间距保证漏掉的机会变小。我一般把F范围再分成两段扫低F段用细网格高F段用粗网格漏根概率能压到很低。5.2 模式跳线根不排序直接连线曲线互相串现象单独求每个F的根再画图S1和A1在某个频段发生明显交叉两条曲线像麻花一样拧在一起但对照文献图谱发现真实曲线并不会这样。原因求解器返回的根无序按大小排序后高频段不同模式相互逼近时顺序会发生交换连线就串了。解决用上一节traceModes里的“认领”逻辑每次以上一频率点的根为基准找最近根并且同一个根不允许被认领两次。这个方案比排序稳健得多代价只是每帧需要多跑一次min查找对几千个频率点来说耗时完全可以接受。5.3 单位混乱m/s配kHz算出满屏NaN现象所有设置看起来都对但fzero提示端点函数值不异号或者画出来的曲线数量级差了1000倍。检查残差函数里1/cp发现cp如果用的是m/s而cL、cT用的是mm/us速度倒数直接差出3个数量级。原因公式里所有速度必须同一单位制。mm/us加MHz是最舒服的组合但很多人习惯先把波速换成m/s频率用Hz这时如果没同步转换方程就废了。解决在脚本最开头做一次单位断言比如检查cL是否在1到20之间因为常见材料的纵波速度在mm/us制下最多是个位数。单位制统一后不妨再验证一次S0低频板速值这一步能排查掉90%以上的单位问题。5.4 溢出现象高频厚板时残差函数直接爆掉现象板厚达到几十毫米、频率到MHz时数组函数输出Inf或NaN符号变化判断完全失效求根器一个根都返回不了。原因sin/cos复数形式中虚部很大时指数增长double超过约1e308就溢出。归一化后的残差已经比tan形式好很多但极端频厚积下依然有风险。解决把残差整体除一个尺度因子cosh(imag(pd))·cosh(imag(qd))让sin/cos项变成有界组合。关键代码就一行% 在lambResidualNorm返回前加一句防溢出缩放 scale cosh(imag(pd)) .* cosh(imag(qd)); res res ./ scale;这样处理不会改变零点位置因为除以的是一个恒正的数。高频厚板场景里这条路比直接改N有效得多。5.5 搜索域太窄cp上限定到cL高次模式全丢现象S0、A0、S1、A1都在往上S2、A2、S3怎么也出不来或者到高频末尾突然断掉。原因高阶模式在截止频率附近相速度趋于无穷搜索域上限如果只到cL它们的根根本没有进入求解区间。解决把上界放宽到4倍cL以上计算量增加不多边界效果却很显著。低频段还要注意下界A0在F很小时相速度低于0.1cT下界用0.02cT在大多数情况下够用如果扫F低于0.05建议下界降到0.005cT。6. 验证与延伸把频散曲线落到激励设计上6.1 三种低成本验证手段逐点验证极限值是第一种手段。除了前面提到的板波速度和瑞利波速还可以看A0在F趋近0时相速度和群速度同时趋近0这是弯曲波低频行为。第二种手段是和文献里的标准谱图对比重点看S1群速度曲线的“凹陷”位置这个特征类似模式指纹不同材料参数下凹陷位置变化明显如果对得上代码基本没有大问题。第三种手段是切换材料参数比如换成铝观察S0低频板速是否等于对应公式值能检验代码对材料参数的依赖是否正确。6.2 从频散曲线到激励信号设计工程上拿到频散曲线后最常用的场景是选择激励频率。做法是看群速度曲线的平坦区平坦区意味着传播时波形包络畸变小。比如要激发S0模式就选S0群速度曲线在目标频段内梯度最小的点作为中心频率。另一个常见需求是估计两个模式到达时间差以便设计时间窗分离波包。到达时间差公式为Δt L·(1/cg1 - 1/cg2)其中L是传播距离cg1和cg2是两个模式的群速度。从曲线上读出对应F下的群速度值直接代入就能算出分离窗口。这套MATLAB实现里附带的扫频和模式追踪脚本可以直接改参数跑通不需要再从零搭框架。6.3 一点个人习惯以前我调一套复合材料的参数时把搜索域上限设成cL画出的曲线少了一整支高阶模式当时盯着屏幕怀疑是求解器写错了排查了两天才发现是边界问题。从那以后我每次换材料、换板厚都会先跑一遍极限验证确认S0低频板速和高频瑞利速都落在合理范围再去看整条曲线。二维曲线好看是表象求根稳定才是这套实现真正值钱的地方。希望帮到你。本文还有配套的精品资源点击获取