简介本资源是一套专用于计算与绘制空心圆管中超声导波频散曲线的MATLAB开源程序集面向无损检测、结构健康监测领域的工程师、高校科研人员及研究生解决导波模态分析、相速度/群速度求解及频散特性可视化等核心问题。压缩包为RAR格式共24个文件全部为.m脚本涵盖核心求解器pcdisp.m、矩阵行列式计算pcmatdet.m、频散曲线绘制pcplotmatdet2D.m、模态命名pcmodename.m、激励信号建模pcwaveform.m及数值积分intsimpson.m等关键模块总大小仅51KB轻量易部署。已有1421人学习下载代码注释详尽便于理解导波频散方程求解逻辑、参数敏感性分析及二次开发。读者可直接运行复现典型频散曲线修改管径、壁厚、材料参数后快速适配不同工况亦可基于源码拓展模态识别、群速度提取或与实验数据比对功能。1. 管中导波频散曲线为什么先看它再定检测方案前阵子有同行问我3英寸管做导波检测激励频率定多少合适用 T(0,1) 还是 L(0,2)缺陷回波按什么波速算时差。我说你先去看频散曲线。他愣了一下说曲线倒是会用软件画但 PCdisp 算出来的那堆线他不知道怎么对应到实际检测参数上。这个问题其实非常典型。管中导波频散曲线是所有超声导波检测方案设计的起点。你选的模态、定的激励频率、算的传播时差本质上都是在跟这张图打交道。而 PCdisp 这类计算程序恰好就是把波动方程的解变成工程上能直接用的曲线的最顺手工具——它算的是空心圆柱波导里导波模态的相速度和群速度随频率的变化关系也就是频散曲线。这篇文章我会从计算原理、源程序编译运行、读图取参、数值坑和代码整理几个角度把整个链路讲透。1.1 波速不是常数检测就必须先查表很多人刚开始接触导波时习惯性认为波在固体里传播的速度就是材料声速测到回波时间除以 2 再乘声速就是距离。这个思路在常规超声直探头场景下基本成立但在管道导波场景下会出问题因为导波的速度不是常数。频散dispersion的本质是波在有限几何截面比如管壁中传播时不同频率成分感受到的边界约束不同导致相速度随频率变化。打个比方白光经过三棱镜会散成彩带是因为不同波长的光在玻璃里折射率不同。导波在管道里也一样不同频率的折射率不同传播速度就不同。这里必须区分两个速度相速度 cp ω/k某个频率成分的等相位面移动速度群速度 cg dω/dk整个波包脉冲的能量包络移动速度。实际检测中我们发射的是一个有限持续时间的脉冲信号能量和信息是以群速度传播的。所以缺陷定位、时差换算全都用群速度而模态识别、波长计算用相速度。频散曲线必须把这两条都算出来这就是 PCdisp 这类程序的核心输出。1.2 三类模态族L、T、F 怎么记管道导波模态按圆周方向和谐振结构分为三个族模态族名称特征典型用途L(0,m)纵向模态轴对陈质点位移主要在轴向和径向长距离检测L(0,2) 常用T(0,m)扭转模态轴对称质点位移主要在周向长距离检测首选T(0,1) 最常用F(n,m)弯曲模态非轴对称n 为周向阶数m 为模阶短距离、缺陷精细表征括号里第一个数字是周向阶数 n0 表示轴对称第二个数字 m 是同一族里按截止频率从低到高的模态序号。L 是 longitudinalT 是 torsionalF 是 flexural。其中 T(0,1) 是管道导波检测里公认的劳模模态。它在很宽的频率范围内几乎不频散群速度接近剪切波速钢大约 3220 m/s没有截止频率对管内介质的敏感性也相对可控。很多长距离管道筛查系统直接以 T(0,1) 为工作模态就冲它那条近乎水平的群速度曲线。1.3 频散曲线的三个用途把频散曲线拿到手主要干三件事选模态。看要检测什么类型的缺陷、能接受多大频散、管道有没有包覆层决定激励哪个模态族选频率。在选定模态的群速度曲线上找平坦段避开陡峭段这样信号包络不会很快被拉散算时间。用选定频率下的精确群速度值做时差-距离换算而不是拿一个笼统的材料声速凑合。这三件事每件都离不开准确的频散数据。用手算或用近似公式在简单平板情形下还能对付一旦进入管道这种圆柱几何模态耦合复杂必须靠数值求解这就是 PCdisp 存在的意义。2. PCdisp 计算原理从波动方程到可用的曲线2.1 圆柱波导的数学描述PCdisp 计算的物理模型是弹性空心圆柱管道材料假设为各向同性线弹性。控制方程是纳维方程在柱坐标下把位移做亥姆霍兹分解引入标量势函数和矢量势函数然后设解沿轴向以 exp(i(kz-ωt)) 形式传播。将势函数代入后径向部分会得到贝塞尔方程解由第一类和第二类贝塞尔函数组合而成。管道的内外表面都是自由边界这就要求应力分量在 r 内半径 和 r 外半径 处同时为零。把位移、应变、应力全部用势函数表示代入四个边界条件内外表面各三个应力分量共六个但轴对称时部分解耦就得到一个关于未知系数的齐次方程组。这个方程组有非零解的条件是系数矩阵的行列式等于零。行列式为零就构成了频散方程——一个包含频率 ω、波数 k、材料参数和几何参数的超越方程。PCdisp 做的事情本质上就是在给定频率下求解这个超越方程得到对应的波数再换算成相速度。对于 T(0,m) 扭转模态问题会明显简化因为扭转波只涉及周向位移边界条件退化为独立的单方程求解稳定。而 L(0,m) 和 F(n,m) 需要处理 2x2 或更高阶的行列式数值上要小心。2.2 特征方程与零点的求法超越方程没有解析通解只能数值求根。PCdisp 这类程序典型做法是对每个周向阶数 n在感兴趣的频率范围内对每个频率点扫描波数 k或先扫描相速度区间计算行列式值通过符号变化判断是否存在零点用二分法或牛顿迭代法精确定位零点频率步进推进时根据上一频率点的模态顺序做分支跟踪避免把不同模态的曲线接错。这里最考验程序功底的是模态跟踪。频率足够高时模态数量多不同模态的曲线可能交叉或接近行列式零点之间的间距很小如果频率步长太大很容易漏模态或跳支。后面我会专门讲这些坑。2.3 群速度怎么算群速度 cg dω/dk有两条路线解析路线对频散方程关于 k 求偏导解出 dω/dk 的显式表达式数值路线在频散曲线上取相邻两个频率点对应的波数用差分近似 dω/dk。数值路线实现简单但曲线斜率变化剧烈的区域比如截止频率附近误差会放大。PCdisp 的常见实现是混合策略用数值差分得到初值再结合局部插值平滑。我们自己用源程序时如果发现群速度曲线出现锯齿状抖动优先怀疑差分步长是否合适。2.4 输入参数与输出内容使用 PCdisp 前需要明确几类输入参数几何参数管道外半径、内半径或壁厚材料参数杨氏模量 E、泊松比 ν、密度 ρ有的版本直接要求拉梅常数或纵波、横波声速计算参数起始频率、终止频率、频率步长或点数、最大周向阶数 n、需要的模态族输出控制是否导出文本、是否同时输出相速度和群速度。输出通常是每个模态一组数据包含频率、相速度、群速度有的版本还有衰减或波结构成分形式为文本文件方便后续用 Python、MATLAB 或者 Origin 绘图。PCdisp 的定位就是算然后交给用户自己处理所以它的核心竞争力在求解器的稳定性和完备性而不是界面有多漂亮。3. 源程序编译与运行从拿到代码到跑出第一张图3.1 源程序结构长什么样PCdisp 的源程序在流传过程中有多个版本文件组织不完全一样但大体逃不出下面这几类模块我以常见版本结构为例说明主程序文件负责读输入、循环频率、调用求解器、写输出特征方程模块根据模态族构造行列式或方程组贝塞尔函数计算模块第一类和第二类贝塞尔函数及其导数求根模块二分法、牛顿法或组合算法输入输出模块解析输入文件、格式化输出结果。拿到源码后第一件事不是急着编译而是先看 README 或注释里的版本说明搞清楚它依赖哪些库、预期用什么编译器。老版本程序很多是为 Fortran 77/90 写的新的可能用 C/C 封装了界面。不同版本差异很大这里我只能说一个通用的处理路径。3.2 编译环境如果拿到的是 Fortran 源程序Linux 或 Windows 下用 gfortran 编译通常最省事。一个典型的最小化编译命令长这样gfortran -O2 -o pcdisp main.f90 dispersion.f90 bessel.f90 rootfinder.f90 -lmWindows 下如果没配过命令行环境可以用 Code::Blocks 加 MinGW或者直接装 WSL 在 Ubuntu 里编。编译报错时九成问题出在隐性类型和数组越界上老 Fortran 代码对编译器检查选项很敏感。我自己习惯编译时加-ffixed-line-length-none -fdefault-real-8gfortran 选项前者避免固定格式行宽报错后者把默认实数提升为双精度。导波频散计算对精度要求高尤其是接近截止频率时单精度很容易让求根过程漂移。3.3 输入文件准备以 3 英寸 Schedule 40 钢管为例外径 88.9 mm壁厚 5.49 mm我们算 10~100 kHz 范围内的前几个模态。输入文件大致是这种形式不同版本字段顺序可能不同以实际源码注释为准# PCdisp input # Material: steel 210.0e9 # Youngs modulus (Pa) 0.3 # Poisson ratio 7800.0 # density (kg/m3) # Geometry: pipe (m) 44.45e-3 # outer radius 38.96e-3 # inner radius # Frequency range (Hz) and steps 10000.0 # fmin 100000.0 # fmax 500 # number of frequency steps # circumferential order 3 # max n单位是这里最容易翻车的地方。源码里几何如果写的是米你填毫米就会得到完全离谱的曲线如果材料参数用的是 GPa 而不是 Pa频率轴也会全偏。拿到任何一套源程序先确认内部单位制再准备输入文件。为什么选 10~100 kHz这是管道导波检测的典型激励频段。频率太高时波长太短对管道表面状况和包覆层过于敏感衰减严重频率太低时波长比管径大太多模态选择性差。具体定多少正是要用频散曲线来论证的。3.4 运行与结果验证编译通过、输入文件准备好后命令行运行./pcdisp input.txt跑完会生成若干输出文件命名通常是 L01.txt、T01.txt、F11.txt 这样的形式每行对应一个频率点列分别是频率、相速度、群速度。拿到结果先别急着画图做两个 sanity checkT(0,1) 模态在整个频段内的群速度应该稳定在剪切波速附近对钢来说约 3220 m/s。如果偏出去很多说明材料参数或单位有问题L(0,1) 模态在低频极限应接近一维杆的纵波速度 sqrt(E/ρ)对钢大约 5190 m/s。这是从弹性理论可以直接推出来的低频极限是很好的验证锚点。这两个检查通过了基本可以认定程序跑通了。3.5 一个够用的绘图脚本PCdisp 的文本输出用 Python 处理很方便。下面是读取 L(0,1) 和 T(0,1) 群速度并绘图的脚本框架import numpy as np import matplotlib.pyplot as plt def load_mode(filename): data np.loadtxt(filename) return data[:, 0], data[:, 1], data[:, 2] # f, cp, cg f1, cp1, cg1 load_mode(L01.txt) f2, cp2, cg2 load_mode(T01.txt) plt.figure(figsize(8, 5)) plt.plot(f1 / 1e3, cg1 / 1e3, labelL(0,1) group velocity) plt.plot(f2 / 1e3, cg2 / 1e3, labelT(0,1) group velocity) plt.xlabel(Frequency (kHz)) plt.ylabel(Group velocity (m/s)) plt.grid(True, alpha0.3) plt.legend() plt.show()把 L、T、F 各模态都画在同一张图上就是完整的频散曲线图。图中那些从低频往高频延伸的线每条都代表一个模态线越平缓说明该频段频散越弱、越适合检测。4. 把曲线变成检测参数读图与时间换算4.1 选模态选频率的工程逻辑拿到频散曲线图后选模态的经典逻辑是做长距离筛查优先 T(0,1)。看它的群速度曲线在 20~100 kHz 是不是一条平线是的话说明在这个频段脉冲不会因为频散被拉散回波识别容易想用纵波模态看 L(0,2)它在较高频段群速度接近纵波速度传播效率高但低频端曲线坡度大选频要避开截止频率附近的陡峭段做缺陷定位和尺寸评估才去考虑 F 族模态用非轴对称波的波结构的周向分布来判断缺陷方位。选频率的核心指标是平坦度。群速度曲线越平信号的波形保真度越好。一个实用的经验法则是在工作频段内群速度变化量控制在 5% 以内回波时差换算的系统误差可接受。如果曲线在候选频率附近斜率过大稍微偏离中心频率波包的不同频率分量就会以显著不同的速度传播把回波拉成一大片。4.2 缺陷回波时间怎么算假设选定了 T(0,1)从频散曲线读出激励频率下的 cg 3220 m/s。被测管道长 50 m缺陷在距离激励端 30 m 处那么回波往返时间t 2z / cg 2 × 30 / 3220 ≈ 18.6 ms这个 18.6 ms 是判定回波位置的基准。如果现场仪器读到的回波在 18.6 ms 附近对应的反射界面就在 30 m 左右。同理如果回波时间偏差 1 ms折算到距离偏差约 1.6 m。这个灵敏度直接决定了检测定位精度前提是 cg 读数准确而不是拿一个拍脑袋的 5000 m/s 去算。用 L(0,2) 时更要小心。它的群速度在低频段变化大你必须用激励频率处曲线上的精确数值做换算不能取一个平均速度。这也是为什么现场工程师手里都会带一张频散数据的表格。4.3 把频散数据落到检测方案里把 PCdisp 算出的文本数据做成一个查询表是工程上很实用的做法。例如导出 CSV 后用几行代码就能按频率查群速度import numpy as np data np.loadtxt(T01.txt, delimiter,) freqs, cgs data[:, 0], data[:, 2] def group_velocity_at(f_target): idx np.argmin(np.abs(freqs - f_target)) return freqs[idx], cgs[idx] print(group_velocity_at(32000)) # 32 kHz 处的群速度这个小函数在实测数据处理时非常有用。从采集软件里读出回波时间调用它拿到准确群速度距离直接算出来整个过程不依赖网络不依赖在线工具原始数据和程序都在你手里随时可以复核。5. 数值计算中的坑模态漏算、跳支与伪根5.1 模态漏算与跳支PCdisp 这类求根程序最常见的毛病是模态漏算。频率升高后一个周向阶数下会出现多个模态它们的曲线在某些频段会非常接近甚至交叉。如果频率步长取得太大行列式符号变化被跨过去零点就丢了或者把 A 模态的尾段错误接续到 B 模态的头上曲线整体跳支。排查办法很直接逐步加密频率步长观察曲线是否发生突变。正常模态曲线应该是光滑连续的任何尖角、突然折断、速度跳变都值得怀疑。另外可以把不同 n 阶的曲线画在一起看相邻阶数的模态曲线在交叉区域往往有规律的排列如果某条线明显插到了别人的位置就要查跟踪逻辑。5.2 截止频率附近的伪根每个高阶模态都有截止频率。在截止频率附近相速度趋于无穷大行列式的零点和其他物理量的数值特性会变得很差。数值程序容易在这里找到伪根——行列式确实很小但不是真实的物理模态。伪根的特征是曲线只在很窄的频率范围内出现随后消失或者算出的群速度为负、超过材料纵波速度等不合理值。出现这种情况先看程序输出里有没有对行列式归一化处理再看截止频率附近是否用了足够细的频率步长。实在不行把该频段的结果和文献上的公开频散曲线对比人工剔除异常支。5.3 材料参数敏感性与单位错误频散曲线的形状对剪切波速 cT 特别敏感。钢的标准参数 E 210 GPa、ν 0.3、ρ 7800 kg/m³ 时cT ≈ 3220 m/s。如果实际管道材料的 cT 偏差 5%所有模态曲线的频率轴位置都会有明显偏移。所以算曲线前最好用实测声速横波和纵波反推材料参数而不是直接套用书本值。单位错误是另一类高频事故。几何参数用 mm 还是 m、频率用 Hz 还是 kHz、弹性模量用 Pa 还是 GPa任何一处不一致都会让结果完全失真。我见过有人把外半径 44.45 直接按 44.45 m 填进去算出来的频散曲线在几赫兹处出现一堆模态——这显然不对。5.4 对不上时的排查清单如果 PCdisp 的曲线和实测数据或文献曲线对不上按顺序检查检查项说明单位制几何、频率、模量单位是否与源码内部一致材料参数是否用了与该管道实际状态匹配的声速几何尺寸内外半径是否填写正确壁厚是否算对频率范围是否超出了该程序设计的有效范围模态命名同一符号在不同文献里可能指不同模态数值参数频率步长、迭代精度是否足够这套清单我从科研到现场项目用过很多次绝大多数曲线对不上问题都能在最后两项之前解决。6. 源程序整理与鉴材一次热搜引发的代码管理提醒6.1 为什么这类补材料情况屡见不鲜最近在行业热搜里看到一句话提交的源程序鉴别材料中出现大量与已登记或申请的软件源程序相同的内容。这其实是很多用开源或学术共享代码做项目的人都会撞上的问题——PCdisp 源程序本身是公开流传的很多人拿过来改改用到提交成果材料时直接把原版代码或大段未修改的模块交了上去。我不谈具体流程和规定只从做工程的角度说一句任何一套代码如果你打算作为自己的成果材料交出去那它必须是你真正常握、能逐行解释、并且确实做了二次开发的版本。直接复制原版源程序当交付内容与其说是流程问题不如说是代码管理意识的问题——一来没有体现你的工作二来一旦被要求补正重新整理的成本远高于一开始就规范。6.2 代码怎么整理才算自己的结合我这些年折腾各种开源计算程序的经验把 PCdisp 这类代码变成自己的成果至少要做到以下几点逐行读懂。不要留任何这段没看明白但能跑的代码。读不懂就加注释读懂了用自己的话重写注释重构模块。把原版的单文件大函数拆成有明确职责的小模块变量名按你的习惯重命名消除全局变量泛滥的问题写清物理注释。比如某段是在算 T 模态的行列式注释里写明对应的边界条件和应力分量公式这既是整理也是加深理解保留验证用例。把 3 英寸钢管的输入、输出和验证结果存成一个 test case能随时复现证明程序行为符合预期明确区分来源。哪些模块是沿用原版思路重写的哪些是全新加的最好在 README 里写清楚。这不是推卸而是负责任地说明你的增量贡献。我自己习惯在新代码开头加一段版本变更记录注明基于哪个版本、改了什么、解决了哪些数值问题。这不仅是给别人看的半年后你自己回来看也能快速找回上下文。6.3 代码管理的日常习惯从这次热搜还折射出一个更普遍的问题很多人不把为项目写的脚本和配置当代码管理。其实 PCdisp 的输入文件、后处理脚本、测试数据这些和源程序一样重要。我建议至少做这几件事用 Git 管理整个目录包括输入文件、输出样例、绘图脚本每次修改输入参数或源码都提交一次commit message 写清楚改了什么README 里记录编译命令、运行方式、已知问题和对应解决办法把验证用的文献对比图保存起来作为程序可信度的证据。这些习惯花不了多少时间但能在需要交材料、复现结果、或者换电脑重装环境时替你省下大量补材料的时间。尤其是像 PCdisp 这种从学术圈流出来的程序版本多、依赖杂、注释风格各异不好好整理放三个月你自己都未必能重新跑通。我的个人习惯是拿到任何一套新源码先建一个干净目录把原始版本原封不动存一份再另建工作目录开始改。原始版本是基准线出了问题可以随时 diff 回来看改动。另外每次跑完一批曲线我会顺手把输入参数、输出文件和当时用的材料参数一起打包存档文件名带上日期。这些做法看起来很笨但在你需要向别人说明这条曲线是用什么参数、什么版本代码、什么步骤算出来的时这套记录就是最快的答案。本文还有配套的精品资源点击获取 从频散曲线到检测方案：PCdisp源程序解析与工程应用
简介本资源是一套专用于计算与绘制空心圆管中超声导波频散曲线的MATLAB开源程序集面向无损检测、结构健康监测领域的工程师、高校科研人员及研究生解决导波模态分析、相速度/群速度求解及频散特性可视化等核心问题。压缩包为RAR格式共24个文件全部为.m脚本涵盖核心求解器pcdisp.m、矩阵行列式计算pcmatdet.m、频散曲线绘制pcplotmatdet2D.m、模态命名pcmodename.m、激励信号建模pcwaveform.m及数值积分intsimpson.m等关键模块总大小仅51KB轻量易部署。已有1421人学习下载代码注释详尽便于理解导波频散方程求解逻辑、参数敏感性分析及二次开发。读者可直接运行复现典型频散曲线修改管径、壁厚、材料参数后快速适配不同工况亦可基于源码拓展模态识别、群速度提取或与实验数据比对功能。1. 管中导波频散曲线为什么先看它再定检测方案前阵子有同行问我3英寸管做导波检测激励频率定多少合适用 T(0,1) 还是 L(0,2)缺陷回波按什么波速算时差。我说你先去看频散曲线。他愣了一下说曲线倒是会用软件画但 PCdisp 算出来的那堆线他不知道怎么对应到实际检测参数上。这个问题其实非常典型。管中导波频散曲线是所有超声导波检测方案设计的起点。你选的模态、定的激励频率、算的传播时差本质上都是在跟这张图打交道。而 PCdisp 这类计算程序恰好就是把波动方程的解变成工程上能直接用的曲线的最顺手工具——它算的是空心圆柱波导里导波模态的相速度和群速度随频率的变化关系也就是频散曲线。这篇文章我会从计算原理、源程序编译运行、读图取参、数值坑和代码整理几个角度把整个链路讲透。1.1 波速不是常数检测就必须先查表很多人刚开始接触导波时习惯性认为波在固体里传播的速度就是材料声速测到回波时间除以 2 再乘声速就是距离。这个思路在常规超声直探头场景下基本成立但在管道导波场景下会出问题因为导波的速度不是常数。频散dispersion的本质是波在有限几何截面比如管壁中传播时不同频率成分感受到的边界约束不同导致相速度随频率变化。打个比方白光经过三棱镜会散成彩带是因为不同波长的光在玻璃里折射率不同。导波在管道里也一样不同频率的折射率不同传播速度就不同。这里必须区分两个速度相速度 cp ω/k某个频率成分的等相位面移动速度群速度 cg dω/dk整个波包脉冲的能量包络移动速度。实际检测中我们发射的是一个有限持续时间的脉冲信号能量和信息是以群速度传播的。所以缺陷定位、时差换算全都用群速度而模态识别、波长计算用相速度。频散曲线必须把这两条都算出来这就是 PCdisp 这类程序的核心输出。1.2 三类模态族L、T、F 怎么记管道导波模态按圆周方向和谐振结构分为三个族模态族名称特征典型用途L(0,m)纵向模态轴对陈质点位移主要在轴向和径向长距离检测L(0,2) 常用T(0,m)扭转模态轴对称质点位移主要在周向长距离检测首选T(0,1) 最常用F(n,m)弯曲模态非轴对称n 为周向阶数m 为模阶短距离、缺陷精细表征括号里第一个数字是周向阶数 n0 表示轴对称第二个数字 m 是同一族里按截止频率从低到高的模态序号。L 是 longitudinalT 是 torsionalF 是 flexural。其中 T(0,1) 是管道导波检测里公认的劳模模态。它在很宽的频率范围内几乎不频散群速度接近剪切波速钢大约 3220 m/s没有截止频率对管内介质的敏感性也相对可控。很多长距离管道筛查系统直接以 T(0,1) 为工作模态就冲它那条近乎水平的群速度曲线。1.3 频散曲线的三个用途把频散曲线拿到手主要干三件事选模态。看要检测什么类型的缺陷、能接受多大频散、管道有没有包覆层决定激励哪个模态族选频率。在选定模态的群速度曲线上找平坦段避开陡峭段这样信号包络不会很快被拉散算时间。用选定频率下的精确群速度值做时差-距离换算而不是拿一个笼统的材料声速凑合。这三件事每件都离不开准确的频散数据。用手算或用近似公式在简单平板情形下还能对付一旦进入管道这种圆柱几何模态耦合复杂必须靠数值求解这就是 PCdisp 存在的意义。2. PCdisp 计算原理从波动方程到可用的曲线2.1 圆柱波导的数学描述PCdisp 计算的物理模型是弹性空心圆柱管道材料假设为各向同性线弹性。控制方程是纳维方程在柱坐标下把位移做亥姆霍兹分解引入标量势函数和矢量势函数然后设解沿轴向以 exp(i(kz-ωt)) 形式传播。将势函数代入后径向部分会得到贝塞尔方程解由第一类和第二类贝塞尔函数组合而成。管道的内外表面都是自由边界这就要求应力分量在 r 内半径 和 r 外半径 处同时为零。把位移、应变、应力全部用势函数表示代入四个边界条件内外表面各三个应力分量共六个但轴对称时部分解耦就得到一个关于未知系数的齐次方程组。这个方程组有非零解的条件是系数矩阵的行列式等于零。行列式为零就构成了频散方程——一个包含频率 ω、波数 k、材料参数和几何参数的超越方程。PCdisp 做的事情本质上就是在给定频率下求解这个超越方程得到对应的波数再换算成相速度。对于 T(0,m) 扭转模态问题会明显简化因为扭转波只涉及周向位移边界条件退化为独立的单方程求解稳定。而 L(0,m) 和 F(n,m) 需要处理 2x2 或更高阶的行列式数值上要小心。2.2 特征方程与零点的求法超越方程没有解析通解只能数值求根。PCdisp 这类程序典型做法是对每个周向阶数 n在感兴趣的频率范围内对每个频率点扫描波数 k或先扫描相速度区间计算行列式值通过符号变化判断是否存在零点用二分法或牛顿迭代法精确定位零点频率步进推进时根据上一频率点的模态顺序做分支跟踪避免把不同模态的曲线接错。这里最考验程序功底的是模态跟踪。频率足够高时模态数量多不同模态的曲线可能交叉或接近行列式零点之间的间距很小如果频率步长太大很容易漏模态或跳支。后面我会专门讲这些坑。2.3 群速度怎么算群速度 cg dω/dk有两条路线解析路线对频散方程关于 k 求偏导解出 dω/dk 的显式表达式数值路线在频散曲线上取相邻两个频率点对应的波数用差分近似 dω/dk。数值路线实现简单但曲线斜率变化剧烈的区域比如截止频率附近误差会放大。PCdisp 的常见实现是混合策略用数值差分得到初值再结合局部插值平滑。我们自己用源程序时如果发现群速度曲线出现锯齿状抖动优先怀疑差分步长是否合适。2.4 输入参数与输出内容使用 PCdisp 前需要明确几类输入参数几何参数管道外半径、内半径或壁厚材料参数杨氏模量 E、泊松比 ν、密度 ρ有的版本直接要求拉梅常数或纵波、横波声速计算参数起始频率、终止频率、频率步长或点数、最大周向阶数 n、需要的模态族输出控制是否导出文本、是否同时输出相速度和群速度。输出通常是每个模态一组数据包含频率、相速度、群速度有的版本还有衰减或波结构成分形式为文本文件方便后续用 Python、MATLAB 或者 Origin 绘图。PCdisp 的定位就是算然后交给用户自己处理所以它的核心竞争力在求解器的稳定性和完备性而不是界面有多漂亮。3. 源程序编译与运行从拿到代码到跑出第一张图3.1 源程序结构长什么样PCdisp 的源程序在流传过程中有多个版本文件组织不完全一样但大体逃不出下面这几类模块我以常见版本结构为例说明主程序文件负责读输入、循环频率、调用求解器、写输出特征方程模块根据模态族构造行列式或方程组贝塞尔函数计算模块第一类和第二类贝塞尔函数及其导数求根模块二分法、牛顿法或组合算法输入输出模块解析输入文件、格式化输出结果。拿到源码后第一件事不是急着编译而是先看 README 或注释里的版本说明搞清楚它依赖哪些库、预期用什么编译器。老版本程序很多是为 Fortran 77/90 写的新的可能用 C/C 封装了界面。不同版本差异很大这里我只能说一个通用的处理路径。3.2 编译环境如果拿到的是 Fortran 源程序Linux 或 Windows 下用 gfortran 编译通常最省事。一个典型的最小化编译命令长这样gfortran -O2 -o pcdisp main.f90 dispersion.f90 bessel.f90 rootfinder.f90 -lmWindows 下如果没配过命令行环境可以用 Code::Blocks 加 MinGW或者直接装 WSL 在 Ubuntu 里编。编译报错时九成问题出在隐性类型和数组越界上老 Fortran 代码对编译器检查选项很敏感。我自己习惯编译时加-ffixed-line-length-none -fdefault-real-8gfortran 选项前者避免固定格式行宽报错后者把默认实数提升为双精度。导波频散计算对精度要求高尤其是接近截止频率时单精度很容易让求根过程漂移。3.3 输入文件准备以 3 英寸 Schedule 40 钢管为例外径 88.9 mm壁厚 5.49 mm我们算 10~100 kHz 范围内的前几个模态。输入文件大致是这种形式不同版本字段顺序可能不同以实际源码注释为准# PCdisp input # Material: steel 210.0e9 # Youngs modulus (Pa) 0.3 # Poisson ratio 7800.0 # density (kg/m3) # Geometry: pipe (m) 44.45e-3 # outer radius 38.96e-3 # inner radius # Frequency range (Hz) and steps 10000.0 # fmin 100000.0 # fmax 500 # number of frequency steps # circumferential order 3 # max n单位是这里最容易翻车的地方。源码里几何如果写的是米你填毫米就会得到完全离谱的曲线如果材料参数用的是 GPa 而不是 Pa频率轴也会全偏。拿到任何一套源程序先确认内部单位制再准备输入文件。为什么选 10~100 kHz这是管道导波检测的典型激励频段。频率太高时波长太短对管道表面状况和包覆层过于敏感衰减严重频率太低时波长比管径大太多模态选择性差。具体定多少正是要用频散曲线来论证的。3.4 运行与结果验证编译通过、输入文件准备好后命令行运行./pcdisp input.txt跑完会生成若干输出文件命名通常是 L01.txt、T01.txt、F11.txt 这样的形式每行对应一个频率点列分别是频率、相速度、群速度。拿到结果先别急着画图做两个 sanity checkT(0,1) 模态在整个频段内的群速度应该稳定在剪切波速附近对钢来说约 3220 m/s。如果偏出去很多说明材料参数或单位有问题L(0,1) 模态在低频极限应接近一维杆的纵波速度 sqrt(E/ρ)对钢大约 5190 m/s。这是从弹性理论可以直接推出来的低频极限是很好的验证锚点。这两个检查通过了基本可以认定程序跑通了。3.5 一个够用的绘图脚本PCdisp 的文本输出用 Python 处理很方便。下面是读取 L(0,1) 和 T(0,1) 群速度并绘图的脚本框架import numpy as np import matplotlib.pyplot as plt def load_mode(filename): data np.loadtxt(filename) return data[:, 0], data[:, 1], data[:, 2] # f, cp, cg f1, cp1, cg1 load_mode(L01.txt) f2, cp2, cg2 load_mode(T01.txt) plt.figure(figsize(8, 5)) plt.plot(f1 / 1e3, cg1 / 1e3, labelL(0,1) group velocity) plt.plot(f2 / 1e3, cg2 / 1e3, labelT(0,1) group velocity) plt.xlabel(Frequency (kHz)) plt.ylabel(Group velocity (m/s)) plt.grid(True, alpha0.3) plt.legend() plt.show()把 L、T、F 各模态都画在同一张图上就是完整的频散曲线图。图中那些从低频往高频延伸的线每条都代表一个模态线越平缓说明该频段频散越弱、越适合检测。4. 把曲线变成检测参数读图与时间换算4.1 选模态选频率的工程逻辑拿到频散曲线图后选模态的经典逻辑是做长距离筛查优先 T(0,1)。看它的群速度曲线在 20~100 kHz 是不是一条平线是的话说明在这个频段脉冲不会因为频散被拉散回波识别容易想用纵波模态看 L(0,2)它在较高频段群速度接近纵波速度传播效率高但低频端曲线坡度大选频要避开截止频率附近的陡峭段做缺陷定位和尺寸评估才去考虑 F 族模态用非轴对称波的波结构的周向分布来判断缺陷方位。选频率的核心指标是平坦度。群速度曲线越平信号的波形保真度越好。一个实用的经验法则是在工作频段内群速度变化量控制在 5% 以内回波时差换算的系统误差可接受。如果曲线在候选频率附近斜率过大稍微偏离中心频率波包的不同频率分量就会以显著不同的速度传播把回波拉成一大片。4.2 缺陷回波时间怎么算假设选定了 T(0,1)从频散曲线读出激励频率下的 cg 3220 m/s。被测管道长 50 m缺陷在距离激励端 30 m 处那么回波往返时间t 2z / cg 2 × 30 / 3220 ≈ 18.6 ms这个 18.6 ms 是判定回波位置的基准。如果现场仪器读到的回波在 18.6 ms 附近对应的反射界面就在 30 m 左右。同理如果回波时间偏差 1 ms折算到距离偏差约 1.6 m。这个灵敏度直接决定了检测定位精度前提是 cg 读数准确而不是拿一个拍脑袋的 5000 m/s 去算。用 L(0,2) 时更要小心。它的群速度在低频段变化大你必须用激励频率处曲线上的精确数值做换算不能取一个平均速度。这也是为什么现场工程师手里都会带一张频散数据的表格。4.3 把频散数据落到检测方案里把 PCdisp 算出的文本数据做成一个查询表是工程上很实用的做法。例如导出 CSV 后用几行代码就能按频率查群速度import numpy as np data np.loadtxt(T01.txt, delimiter,) freqs, cgs data[:, 0], data[:, 2] def group_velocity_at(f_target): idx np.argmin(np.abs(freqs - f_target)) return freqs[idx], cgs[idx] print(group_velocity_at(32000)) # 32 kHz 处的群速度这个小函数在实测数据处理时非常有用。从采集软件里读出回波时间调用它拿到准确群速度距离直接算出来整个过程不依赖网络不依赖在线工具原始数据和程序都在你手里随时可以复核。5. 数值计算中的坑模态漏算、跳支与伪根5.1 模态漏算与跳支PCdisp 这类求根程序最常见的毛病是模态漏算。频率升高后一个周向阶数下会出现多个模态它们的曲线在某些频段会非常接近甚至交叉。如果频率步长取得太大行列式符号变化被跨过去零点就丢了或者把 A 模态的尾段错误接续到 B 模态的头上曲线整体跳支。排查办法很直接逐步加密频率步长观察曲线是否发生突变。正常模态曲线应该是光滑连续的任何尖角、突然折断、速度跳变都值得怀疑。另外可以把不同 n 阶的曲线画在一起看相邻阶数的模态曲线在交叉区域往往有规律的排列如果某条线明显插到了别人的位置就要查跟踪逻辑。5.2 截止频率附近的伪根每个高阶模态都有截止频率。在截止频率附近相速度趋于无穷大行列式的零点和其他物理量的数值特性会变得很差。数值程序容易在这里找到伪根——行列式确实很小但不是真实的物理模态。伪根的特征是曲线只在很窄的频率范围内出现随后消失或者算出的群速度为负、超过材料纵波速度等不合理值。出现这种情况先看程序输出里有没有对行列式归一化处理再看截止频率附近是否用了足够细的频率步长。实在不行把该频段的结果和文献上的公开频散曲线对比人工剔除异常支。5.3 材料参数敏感性与单位错误频散曲线的形状对剪切波速 cT 特别敏感。钢的标准参数 E 210 GPa、ν 0.3、ρ 7800 kg/m³ 时cT ≈ 3220 m/s。如果实际管道材料的 cT 偏差 5%所有模态曲线的频率轴位置都会有明显偏移。所以算曲线前最好用实测声速横波和纵波反推材料参数而不是直接套用书本值。单位错误是另一类高频事故。几何参数用 mm 还是 m、频率用 Hz 还是 kHz、弹性模量用 Pa 还是 GPa任何一处不一致都会让结果完全失真。我见过有人把外半径 44.45 直接按 44.45 m 填进去算出来的频散曲线在几赫兹处出现一堆模态——这显然不对。5.4 对不上时的排查清单如果 PCdisp 的曲线和实测数据或文献曲线对不上按顺序检查检查项说明单位制几何、频率、模量单位是否与源码内部一致材料参数是否用了与该管道实际状态匹配的声速几何尺寸内外半径是否填写正确壁厚是否算对频率范围是否超出了该程序设计的有效范围模态命名同一符号在不同文献里可能指不同模态数值参数频率步长、迭代精度是否足够这套清单我从科研到现场项目用过很多次绝大多数曲线对不上问题都能在最后两项之前解决。6. 源程序整理与鉴材一次热搜引发的代码管理提醒6.1 为什么这类补材料情况屡见不鲜最近在行业热搜里看到一句话提交的源程序鉴别材料中出现大量与已登记或申请的软件源程序相同的内容。这其实是很多用开源或学术共享代码做项目的人都会撞上的问题——PCdisp 源程序本身是公开流传的很多人拿过来改改用到提交成果材料时直接把原版代码或大段未修改的模块交了上去。我不谈具体流程和规定只从做工程的角度说一句任何一套代码如果你打算作为自己的成果材料交出去那它必须是你真正常握、能逐行解释、并且确实做了二次开发的版本。直接复制原版源程序当交付内容与其说是流程问题不如说是代码管理意识的问题——一来没有体现你的工作二来一旦被要求补正重新整理的成本远高于一开始就规范。6.2 代码怎么整理才算自己的结合我这些年折腾各种开源计算程序的经验把 PCdisp 这类代码变成自己的成果至少要做到以下几点逐行读懂。不要留任何这段没看明白但能跑的代码。读不懂就加注释读懂了用自己的话重写注释重构模块。把原版的单文件大函数拆成有明确职责的小模块变量名按你的习惯重命名消除全局变量泛滥的问题写清物理注释。比如某段是在算 T 模态的行列式注释里写明对应的边界条件和应力分量公式这既是整理也是加深理解保留验证用例。把 3 英寸钢管的输入、输出和验证结果存成一个 test case能随时复现证明程序行为符合预期明确区分来源。哪些模块是沿用原版思路重写的哪些是全新加的最好在 README 里写清楚。这不是推卸而是负责任地说明你的增量贡献。我自己习惯在新代码开头加一段版本变更记录注明基于哪个版本、改了什么、解决了哪些数值问题。这不仅是给别人看的半年后你自己回来看也能快速找回上下文。6.3 代码管理的日常习惯从这次热搜还折射出一个更普遍的问题很多人不把为项目写的脚本和配置当代码管理。其实 PCdisp 的输入文件、后处理脚本、测试数据这些和源程序一样重要。我建议至少做这几件事用 Git 管理整个目录包括输入文件、输出样例、绘图脚本每次修改输入参数或源码都提交一次commit message 写清楚改了什么README 里记录编译命令、运行方式、已知问题和对应解决办法把验证用的文献对比图保存起来作为程序可信度的证据。这些习惯花不了多少时间但能在需要交材料、复现结果、或者换电脑重装环境时替你省下大量补材料的时间。尤其是像 PCdisp 这种从学术圈流出来的程序版本多、依赖杂、注释风格各异不好好整理放三个月你自己都未必能重新跑通。我的个人习惯是拿到任何一套新源码先建一个干净目录把原始版本原封不动存一份再另建工作目录开始改。原始版本是基准线出了问题可以随时 diff 回来看改动。另外每次跑完一批曲线我会顺手把输入参数、输出文件和当时用的材料参数一起打包存档文件名带上日期。这些做法看起来很笨但在你需要向别人说明这条曲线是用什么参数、什么版本代码、什么步骤算出来的时这套记录就是最快的答案。本文还有配套的精品资源点击获取