做超声导波检测的朋友应该都有过这种经历换能器贴上去、A扫信号出来了可一看波形就懵了——同一时刻冒出好几个包络怎么数都数不对。尤其是碰上碳纤维复合材料波形更是乱成一团。这时候你翻文献十篇里有八篇都在说“依据频散曲线可以识别模态”可你连频散曲线怎么来的都不知道更别说自己算了。今天这篇就给大家讲清楚一件事怎么用COMSOL把一个复合材料板的频散曲线算出来并且画成你能直接拿去对模态的图。这个算例是我做复合材料Lamb波检测时总结下来的整套流程不算复杂但细节坑不少。适合正在做导波无损检测、压电换能器激励仿真或者复合材料结构健康监测的工程师和学生。你看完以后应该能把“频散曲线”从论文里的名词变成自己手头可以用的工具。1. 为什么复合材料频散曲线这么难算1.1 从“波速随频率变”这个现象说起先说什么叫频散。你往水面上扔一块石头能看到一圈圈波纹往外扩但如果你仔细看会发现大波纹走得快、小波纹走得慢这就是最简单的频散现象——不同频率的成分传播速度不一样波形在传播过程中就会“散开”。在固体板里这种频散就更明显了。板里传播的导波有两个最基本的家族对称模态S0、S1、S2……和反对称模态A0、A1、A2……。同一个频率下S0和A0同时存在各自的相速度和群速度完全不同同一个模态在不同频率下速度又会变。于是“频率—波数—速度”三条信息相互纠缠必须画成曲线才能讲清楚。对金属各向同性板兰姆波频散方程有现成解析解拿MATLAB写几十行就能画出来。但到了复合材料这里麻烦就来了。复合材料板是各向异性的最典型的是碳纤维增强树脂基复合材料CFRP。沿着纤维方向和垂直纤维方向弹性模量能差出十倍以上纵波、横波在不同方向上的速度都不一样板内传播的导波不再能简单分成纯对称、纯反对称而是变成准对称、准反对称模态频散关系还依赖传播方向。你再想用解析公式去推那推导量能让人写到崩溃。所以工程上基本都改用数值方法算COMSOL就是其中很好用的一类工具。1.2 算频散曲线的三条技术路线说到数值算频散曲线行业内目前常见的有三条路。第一条是半解析有限元法也就是SAFESemi-Analytical Finite Element。这个方法把波的传播方向用解析波数表示只在截面方向做有限元离散计算效率极高尤其适合任意截面形状的波导。但SAFE要自己写程序或者用专门的工具箱不是每个人都有精力去搞。第二条是用通用有限元软件做三维动态显式仿真比如ABAQUS里的频散分析。做法是建一段有限长的板在端部加扫频激励提取沿程多个点的响应做二维傅里叶变换2D-FFT得到频散图。这条路思路直白但要注意的是它本质是模拟“真实传播过程”计算量大而且为了分辨不同模态需要足够长的传播距离和足够密的测点网格跑一次要很久。第三条就是在COMSOL里用“特征频率研究 周期性边界条件”来算。这也是我要详细讲的算例路线。它的核心思想是导波在无限长板中沿一个方向传播时位移场沿传播方向呈简谐分布可以把这个问题转化为“垂直于传播方向的截面上的特征值问题”。这个思路和SAFE有相似之处但COMSOL帮你把有限元部分全部包装好了建模直观、后处理方便计算结果也已经在不少论文里和解析解、试验数据对过精度是可信的。我选择这条路线还有一个很实际的原因后面如果要继续做压电换能器激励导波的完整仿真可以在同一个COMSOL工程里把压电器件和复合板一起建进去多物理场耦合衔接很顺不需要换软件换模型。1.3 算例要解决的工程问题咱们这个算例要解决的工程问题很具体一块厚度为2毫米的单向碳纤维复合材料板我想知道在0到3兆赫兹频率范围内沿纤维方向传播的导波都有哪些模态每个模态的相速度和群速度是多少。有了这条频散曲线我后续做实验时就能预先知道某个频率下激励出来的波包应该出现在哪个时间位置模态识别就有依据了。单向板是最简单也是最基础的复合材料铺层形式理解了它你再去做多向铺层或夹芯板也只是材料参数和层结构的替换思路完全一致。2. 几何与材料建模一开始就要做对的两件事2.1 几何模型不必追求三维二维截面就够了很多人一上手就想建三维板模型觉得越接近实物越准确。但算频散曲线这件事三维模型不仅没必要而且会给自己找麻烦。我们要算的是沿一个固定方向传播的导波模态这个问题的本质是二维平面应变问题取一个截面建模就完全够了。我建的这个算例几何非常简单一个矩形长条宽度取2毫米即板的厚度方向长度取1毫米即波的传播方向的一个代表性胞元。厚度方向必须严格等于真实板厚这是所有结果正确的前提。长度方向取1毫米是出于两方面的考虑一方面满足周期性边界条件的要求另一方面1毫米对于我要扫到的最高波数对应最短波长来说网格已经足够细分不会出现波长内节点数不足的问题。直接在COMSOL里用矩形画就行不需要从SolidWorks导入STEP文件。可能有朋友习惯用CAD建模再导入我提醒一下导入STEP后COMSOL常常会提示一些几何警告比如面缺失、缝隙之类虽然不是每次都致命但会干扰网格划分。能原生的就原生这个算例的几何太简单了没必要走导入这条路。2.2 周期性边界条件到底在干什么这一步是整个算例的灵魂我建议你一定要理解它而不是仅仅照猫画虎点几下。我们的核心假设是板在传播方向上是无限长的。但计算机没法算无限大的东西所以必须截取一段有限长度的胞元同时在这个胞元的左右两端施加周期性边界条件。周期性边界条件的数学本质是Bloch定理它规定一个边界上的位移和另一个边界上的位移之间相差一个固定的相位因子。用公式表示就是u_right u_left * exp(-i·k·L)其中k是我们设定的波数L是胞元长度。所谓“特征频率扫描”就是把这个k当作已知参数然后求解在这个相位约束下板截面能够存在的特征振动频率。这样做的物理意义很直观你给定了一个空间波长由k决定然后问“在这个波长约束下板在厚度方向能产生哪些驻波形状对应的频率是多少”。求解出来的一对频率和波数就构成了频散曲线上的一点。把k从一个范围扫过去整条曲线就出来了。2.3 材料参数别想当然正交各向异性必须输全复合材料板在COMSOL里要用“正交各向异性”材料来定义不能简单填一个各项同性的弹性模量和泊松比否则算出来的频散曲线和实际测试值会差得非常远。单向碳纤维板的材料参数我实测下来用下面这组工程常数作为参考值是比较靠谱的参数数值说明密度ρ1600 kg/m³复合材料密度E1135 GPa沿纤维方向弹性模量E2 E310 GPa垂直纤维方向弹性模量G12 G135 GPa面内剪切模量G233.5 GPa面外剪切模量ν12 ν130.3主泊松比ν230.4面外泊松比材料坐标系要特别注意COMSOL里默认的x、y、z坐标轴对应材料的主方向。在这个算例里我把纤维方向设为x方向厚度方向设为y方向这样波沿x方向传播正好对应“沿纤维方向传播”这一个工况。如果你以后要算不同铺层角的板最简单的方法不是改材料而是建立一个旋转坐标系或者修改铺层方向角。COMSOL里都有现成的功能但那是另外一个话题了这里先不展开。3. 特征频率扫描实操从网格到求解器3.1 网格划分的收敛性验证网格划分这件事我吃了不少亏必须单独拿出来说。频散曲线计算中高频模态对应短波长如果网格太粗高频段的频率值会偏高曲线会出现“往上翘”的假象。如果网格太细计算量成倍增长参数扫描几十步下来一整天就搭进去了。所以网格密度必须有一个合理的折中。我的经验是厚度方向至少要划分6层单元传播方向也就是1毫米的长度方向至少要保证最短波长内有8到10个节点。在这个算例里我用二维自由四边形网格整体最大单元尺寸设为0.05毫米。按这个密度模型总共有大约几千个自由度单次特征频率求解非常快参数扫描60步也只需要几分钟。但这里有一个必须做的操作正式扫描之前先做一次网格收敛性验证。具体做法是固定一个波数比如k 5 rad/mm分别用粗网格0.1mm和细网格0.025mm计算前10阶特征频率对比两者差异。如果高频段的频率偏差超过1%说明网格还不够细需要加密如果偏差很小就用较粗的网格跑省时间。这个验证花不了十分钟但能避免你后面跑完一整轮扫描才发现曲线不连续、白干活。3.2 研究步骤和参数扫描设置物理场选“固体力学”研究类型选“特征频率”。这里的关键在于波数k要定义成一个全局参数然后在研究中用“参数化扫描”来扫描它。参数k的扫描范围怎么定这取决于你想看的频率范围和模态速度范围。对于我这个算例板厚2毫米感兴趣的频率范围是0到3MHz。在这个频率范围内最低的模态速度大概是A0模态在低频段的速度大约1000m/s左右最高的模态速度接近S0模态在低频段的速度大约5400m/s左右。由关系式k 2πf/v可以估算k的扫描范围大致是0.05到20 rad/mm。低频时只需要很小的k高频时需要较大的k为了低频段和高频段都有足够的分辨率我建议用指数分布或者人为将扫描点加密在低频段。具体扫描序列我是这样设的先扫描k从0.05到1步长0.05再扫1到20步长0.5。前后加起来大约60个扫描步每步求解特征频率。还有个细节特征频率研究需要设置“搜索特征值范围”。我在COMSOL里把搜索区间设为0到5MHz确保所有感兴趣的模态都能被找到又不会因为范围太宽而引入大量无关的高阶模态减慢求解速度。3.3 特征模态的筛选方法特征频率求解器算出来的特征值里混着一些我们不要的东西。最典型的是刚体模态频率接近0对应板的整体平移或转动这种模态没有物理意义直接忽略。怎么区分有效模态和无效模态最可靠的方法就是看位移场。在COMSOL的后处理里把变形云图打开逐个看一眼低频模态的形状。有效模态在厚度方向有明显的变形特征准对称模态S0、S1等上下表面位移同相中面没有横向位移准反对称模态A0、A1等上下表面位移反相中面横向位移最大。我一般习惯在同一个k值下把前10阶模态的振型导出来排成一排看能很快把模态家族理清楚。这一步看着笨但确实管用。另外COMSOL特征频率求解器会因为对称性同时给出正负方向传播的模态表现为两个相近的频率值。这是正常的不是错误。你在后处理提取数据的时候可以直接取一个或者求平均对曲线结果没有影响。4. 后处理与频散曲线绘制4.1 从特征频率换算相速度和群速度求解完参数扫描后CHFEComsol Multiphysics特征频率研究会给出每个k值对应的一系列特征频率f。有了f和k相速度的计算是直接的cp ω / k 2πf / k群速度的定义是cg dω / dk。在离散数据点上用差分近似cg Δω / Δk具体操作上我推荐先把数据导出来再算群速度比在COMSOL里直接定义表达式更灵活。COMSOL的“派生值”功能可以把每个扫描步的特征频率整理成表格我把这个表格导出为CSV文件然后用Python脚本做后续计算。这里有个细节导出的数据是按k值分组排序的但每个k值下的特征频率不是自动按模态顺序排好的。如果你直接做差分计算群速度很可能会把S0的频率和A0的频率做差结果完全错误。必须先把频率按“同一物理模态”分组再做差分。4.2 用Python脚本整理数据和绘制曲线我习惯把COMSOL的导出数据处理和绘图都扔给Python脚本的思路大致是这样。首先读取CSV文件文件里应该包含k、特征频率序号、频率值这三列。然后把相同k值下的频率值按从小到大排列提取每个模态的频率数组。接下来按模态频率升序对同一k处的所有频率排序人工根据频率连续性和振型特征给每个模态打上标签A0、S0、A1、S1……。这里关键是模态追踪。如果一个k值下前6个频率是模态1到6那么在下个k值处频率顺序可能发生交换例如S1和A1的f-k曲线在某个k处靠得很近甚至交叉简单按顺序分组就会出错。我用来追踪模态的方法是用上一个k值的特征向量作为初值求解下一个k值然后通过模态置信准则MACModal Assurance Criterion判断新旧模态的对应关系。但这套操作在COMSOL里要写脚本相对麻烦。有更省事的方式算完曲线后先用不排序的散点数据画一张“f-k散点图”你会看到曲线天然地连成一条条分支只是数据点之间偶尔有跳变。用肉眼确认哪些点属于同一条分支然后手动在脚本里按区间截取再分别拟合成光滑曲线。这个方法看起来土但实际工程中非常稳速度也快。4.3 算例图示的内容和读法算例的核心图示通常包含三张图第一张是f-k曲线横轴是波数k纵轴是频率f。这张图最能直观展示模态结构每条从低频延伸到高频的连续分支就是一个模态。各向异性板里A0分支从原点出发斜率随频率变大而变小S0分支在低频段近似直线频率升高后逐渐弯曲。第二张是相速度-频率曲线cp-f横轴频率纵轴相速度。这张图是工程中最常用的直接用来查“某个频率下某个模态的速度是多少”。A0在低频段相速度很低高频逐渐逼近表面波速度S0在低频段速度最高高频逐渐下降。第三张是群速度-频率曲线cg-f这张图比相速度曲线更关键因为实际检测中你看到的波包以群速度传播走时计算必须用群速度而不是相速度。A0的群速度在某一低频处会出现一个极小值这是复合材料板中很典型的现象对应着能速慢、能量堆积严重的频率点实际检测时要特别避开这个区域否则信号衰减严重。画这三张图时我都会把不同模态用不同颜色区分并把模态标签直接标注在曲线旁边。这样即使读者没做过仿真也能一眼看出“这条是S0那条是A0”。5. 常见问题与排查技巧实录5.1 模态排序出错的判定与修正这是所有做频散曲线的人都会遇到的问题我也不例外。第一次跑完参数扫描我导出的频率数据按k分组后直接按频率升序编号就去做群速度差分结果群速度图上一片跳变完全没法看。后来我意识到不同模态的频率曲线在接近时会发生交叉或回避简单的升序编号根本不能满足要求。我采用的解决方法是把数据散点先画出来观察分支的连续性再结合特征模态形状辅助判断。你在自己的算例里如果也遇到群速度曲线乱跳第一反应不要以为是求解错了大概率是模态排序问题。5.2 周期性边界条件设置不生效COMSOL的周期性条件在“固体力学”物理场接口中有专门的设置。创建周期性边界时源边界和目标边界的位移自由度要创建“配对”关系。我遇到过一次设置完成后结果和解析解完全对不上的情况排查发现是目标边界和源边界的法向方向不一致导致相位因子的符号反了相当于波传播方向反相。检查方法是先跑一个简单的k值比如k 0.1对比S0模态在低频段的相速度是否和理论值接近。如果差很多基本就是周期性边界条件的问题。5.3 高频段曲线不平滑高频段频散曲线出现锯齿状抖动几乎都是网格不够细导致的。尤其在纤维方向弹性模量很高的复合材料里S1、A1等高阶模态在高频段的位移场变化剧烈网格节点数不足时特征频率解会偏大而且不稳定。解决方法是回到3.1的收敛性验证环节把高频段对应的最大k值单独固定做一次网格加密对比。如果加密后曲线明显变平滑说明原网格确实不够重新划分网格后再跑一次全扫描。高频段还有一个容易出现的问题是扫出的频率点明显偏离曲线趋势这种点往往是求解器把某个局部高频模态当成全局模态输出。判断标准仍然是看位移场如果在厚度方向只有局部几个单元在振动那基本就是伪模态直接剔除就好。5.4 数据导出时绘图为空在COMSOL里做数据导出时如果碰上“绘图为空”的情况通常是因为当前数据集没有选对。特征频率研究的解存储在一个默认数据集中但你查看参数化扫描中某一步的结果时COMSOL可能会切换到一个不包含解的数据集。解决办法是在数据集列表里手动选中对应扫描步的数据集名称然后再做导出或绘图。这个问题不复杂但很容易让人误以为是求解失败了。6. 几个很值得做的扩展方向6.1 扫描传播角度得到方向频散曲线复合材料各向异性的特点决定了不同传播方向的频散曲线差异很大。我这个算例只算了一个方向如果你实际工程中需要在多个方向布置传感器建议把传播角度作为第二个扫描参数。做法是在材料定义里加入一个旋转角度参数然后用二维参数化扫描同时扫描k和角度一次求解就能得到整个平面的频散数据。6.2 与压电换能器仿真联动COMSOL里做压电换能器激励导波是它比较拿手的场景。你可以把压电片和复合板建在一起压电材料的电极接地用频域扫频求解器得到导波激励响应然后在板表面布置探针点提取时间信号。这就把频散曲线从“理论计算”变成了“激励响应”的完整对接。做这类仿真时需要给板端施加吸收边界或者足够长的阻尼层来模拟无限长板边界否则边界反射波会造成严重的干涉伪影。6.3 用解析解或文献数据做精度验证数值算例最好都要有一个验证过程。对单层单向板你可以在低频段用解析解对比比如A0模态在低频极限时的相速度可以通过经典的板弯曲波速公式大致估算也可以在文献里找到相同材料参数下的频散曲线数据做点对点对比。这个验证不用做得很复杂低频前两三个模态能对得上基本就可以确认模型没有系统性错误。7. 一点个人经验频散曲线这个东西很多新手一开始以为是“查表”就能拿到但真正上手做才发现坑不少。用COMSOL算频散曲线的整个过程建模只占两成时间后面的模态识别、数据整理、曲线修正占了八成时间。不要指望一次跑完就得到完美的图反复迭代几次是很正常的。我在实际做的过程中还有一个体会就是尽量保持模型简单。能算二维的不要上三维能用单向板的不要一开始就上多向铺层先把物理本质理解透再逐步增加复杂度。COMSOL求解频散曲线这个算例真正做到位你对导波在复合材料的传播规律理解会有质的提升后面看实验信号也能更有底气。