做复合材料导波无损检测那阵子，我几乎每周都要和COMSOL里的频散曲线打交道。当时项目里用的是碳纤维增强复合材料层合板，换能器贴上去之后，激励信号一宽频，出来的波包又散又乱。要把信号读懂，就得先把这块板自身的频散曲线算清楚。COMSOL官方案例库里有用特征频率法求解各向同性板频散曲线的案例，但换到复合材料之后，层合板铺层方向一旋转，很多现成设置就失效了，网上能查到的成体系教程也少。这篇文章就把我在COMSOL里从建模、边界条件、网格到后处理的完整流程写下来，包括我踩过的几个坑和对应排查思路，希望能给正在做类似事情的同行省点时间。
1. 频散曲线到底在算什么：从一根波导到一块复材层合板
1.1 频散的本质：波速不是常数，而是频率的函数
频散曲线的“频散”两个字，说白了就是波的传播速度随频率变化的现象。想象你在水池边按下一块木板，水波传出去，长波和短波的传播速度不一样，波形从边缘开始慢慢拉开——这就是最简单的频散。对于固体里的导波来说，情况类似。当Lamb波或者板波在结构里传播时，波数k和角频率ω之间存在一个非线性关系，于是相速度cp = ω/k和群速度cg = dω/dk都会随着频率变化。
这个关系为什么重要？因为无损检测里的导波信号几乎都是宽频激励。你发射一个包含多个频率成分的脉冲，每个频率按自己的速度跑，传到接收位置时，波形被拉长、变形、甚至分裂成多个波包。如果不知道频散曲线，你就很难判断哪个波包对应哪个模态、缺陷反射信号该在什么时候到达。频散曲线画出来，通常就是一张以频率为横轴、相速度或群速度为纵轴的图，上面一条条曲线对应结构可以传播的不同模态，A0、S0、A1、S1这些名字就是它们的代号。
1.2 复合材料凭什么比金属麻烦得多
金属板为什么好算？因为它是各向同性的，任何一个方向材料性质都一样，Lamb波的理论解可以直接用Rayleigh-Lamb方程算，一个脚本就能把频散曲线扫出来。复合材料层合板就不一样了。
第一个麻烦是各向异性。碳纤维增强复合材料单向层里，沿纤维方向的弹性模量可以到130多GPa，垂直纤维方向可能不到10GPa，差了一个量级。波沿着纤维方向跑得飞快，垂直方向慢吞吞，速度自然跟传播方向有关。
第二个麻烦是铺层旋转。工程上的层合板不可能所有层都一个方向，通常按0°/90°/±45°交叉铺设。每一层有自己材料坐标系里的刚度，但求解时必须变换到全局坐标。铺层一多，层与层之间的边界还会造成波模式转换，面内位移和面外位移耦合在一起。等效出来的本构矩阵不再是简单正交各向异性，而是各个分量都可能有耦合的满阵。
第三个麻烦是模态数量。厚度方向分层越多，可传播的模态就越多。金属板低频段可能就A0和S0两条曲线，复材层合板一上来就是一大束曲线挤在一起，模态间距小、交叉点多，后处理配对相当头疼。
1.3 为什么COMSOL特征频率法能解这个问题
讲白了，频散曲线问题最后都归结为一个特征值问题。导波在传播方向上有一个波数k，假设位移场在传播方向呈简谐变化u(x,y,t)=U(y)e^{i(ωt−kx)}，代入弹性动力学方程之后，固定一个k值，就能把它变成一个x方向被“冻结”的二维截面上的特征值问题，解出来一系列对应不同模态的角频率ω。再把k从0扫到几千甚至上万rad/m，每扫一次解一批ω，最后把所有(k,ω)点画出来，频散曲线就出来了。
COMSOL做这件事的核心优势，就是它内置了Floquet周期边界，可以直接把“传播方向简谐变化”这个数学条件变成两个边界上的位移相位关系，不需要用户自己写有限元程序去处理复指数边界条件。固体力学接口里设好周期条件，特征值研究里把k设为扫描参数，其他交给求解器。
2. 两条实现路线对比：COMSOL里的特征频率扫描与频域周期边界法
2.1 路线一：二维截面模型加特征频率扫描
我印象中最常用的路子，就是把三维波导问题压成一个二维截面。比如一个层合板，波沿x方向传播，厚度方向是y，那么我只需要在COMSOL里建一个x方向取一小段长度、y方向等于板厚的矩形截面。这个截面在传播方向的两个端面上施加速度和位移的周期条件，选择“Floquet周期性”，把波数k填进去。然后在“特征值”研究中固定这个k，求解特定频率范围内的所有特征模态。
这条路的好处非常明显：二维模型自由度少，网格可以画得很细，特征值求解很快。一个十层复材板的截面模型，映射网格画出来可能也就一两万个自由度，单个k点求解几秒到几十秒，扫几百个k值完全可行。而且后处理可以直接看厚度方向的振型分布，判断模态是对称还是反对称，这对模态标签非常有帮助。
要注意的是二维模型必须在平面应变和平面应力之间选一个。无限宽板沿面外方向是无限延伸、面外应变为零，对应平面应变假设。大部分无限板导波文献都是用这个假设。如果试件是条状、面外方向很窄，那可能得用平面应力近似，但那种情况下三维截面模型更靠谱。
2.2 路线二：三维单胞加频域周期边界
第二个路线适合更复杂的情况。比如声子晶体、周期性超材料，或者结构在面内两个方向都有周期性，这时候二维截面模型就不够用了。建一个三维单胞，六个面全部施加周期边界，两个方向的波数分量kx、ky分别扫描，特征频率求解得到完整的二维色散面。
三维路线能处理的问题类型确实多，但代价也大。单胞网格一加密，自由度轻松上百万，一次特征值求解就得好几分钟。如果还要做二维波矢扫描，比如kx扫50步、ky扫50步，那就是2500次特征值求解，计算量直接爆炸。所以三维周期性模型通常只在必须考虑面内二维波传播的时候才用，而且一般会用扫掠网格或者降阶方法来控制规模。
还有一条中间路线，就是二维广义平面应变。COMSOL里有“广义平面应变”接口，允许面外应变存在，这在处理不对称层合板的面外耦合时比纯平面应变更精确，同时计算量又远低于三维建模。如果你遇到非对称铺层导致的面内面外耦合比较明显，可以用这个接口做二维截面分析。
2.3 我的选型原则：先二维，需要时再上三维
我的个人习惯是：所有平板类波导问题一律先上二维截面模型。层合板、涂层板、夹芯板、带金属薄层的结构，二维截面都能给出可靠的频散关系。除非有以下几种情况，才会考虑换三维：
- 面内存在二维周期性，比如周期性排布的打孔板、声子晶体。
- 结构面外方向几何急剧变化，无法用单一截面代表。
- 需要考虑面外有限宽度带来的边界效应。
如果只是为了复材层合板算个频散曲线，二维截面模型一定是第一选择。算完有疑问再往三维推，不要一上来就三维建模，时间成本完全不一样。
3. 建模仿真全流程：以碳纤维增强复合层合板为例的设置清单
3.1 几何简化和坐标系设定
拿一个常见的对称层合板举例，假设铺层顺序是[0/90/90/0]，四层，每层0.5mm，总厚度2mm。在COMSOL里新建一个二维模型，x方向取一个周期长度，我通常取1mm到5mm都无所谓，只要网格能覆盖波传播方向的简谐变化就行。y方向就是厚度方向，从0到2mm，按四层分成四个矩形域。
这里有个容易忽略的点：材料坐标系必须单独定义。COMSOL默认材料坐标系和全局坐标系重合，但复合材料的纤维方向通常不是全局x方向。在“材料”节点下添加“材料坐标系”，选择旋转，旋转轴取z轴，旋转角度就是铺层角。0°层不转，90°层绕z轴转90°。这样设置之后，定义正交各向异性材料时，第一主轴才真正对应纤维方向。
3.2 材料参数：从工程常数到弹性矩阵的换算
典型碳纤维/环氧单向层参数可以这样填：
| 参数 | 数值 | 说明 |
|---|---|---|
| E1 | 135 GPa | 纤维方向弹性模量 |
| E2 | 9.8 GPa | 垂直纤维方向弹性模量 |
| E3 | 9.8 GPa | 厚度方向弹性模量 |
| G12 | 4.7 GPa | 面内剪切模量 |
| G13 | 4.7 GPa | 沿厚度剪切模量 |
| G23 | 3.4 GPa | 横向剪切模量 |
| ν12 | 0.30 | 主泊松比 |
| ν13 | 0.30 | 沿厚度泊松比 |
| ν23 | 0.35 | 横向泊松比 |
| ρ | 1600 kg/m³ | 密度 |
在COMSOL材料属性里，可以选择直接输入工程常数，程序会自动换算成刚度矩阵，这个够用。如果你需要对结果做精细分析，还可以自己把柔度矩阵S算出来求逆得到C矩阵，再填进去。两种方式结果一致，但要注意工程常数之间必须满足热力学稳定性条件，否则C矩阵不正定，求解器会报错。最常见的问题就是E1、E2填反——纤维方向模量永远最大，别拿90°方向的模量当E1填进去。
3.3 边界条件与特征值研究设置
几何建好、材料属性设好后，关键是边界条件。在“定义”节点添加周期条件，源边界和目标边界分别选中截面x方向的两个端面，类型选“Floquet周期性”。这时COMSOL会要求输入波矢分量(kx, ky, kz)，对于沿x方向传播的板波，就是(kx = k, ky = 0, kz = 0)。这里k就是后面要扫描的参数。
研究类型选择“特征值”。在特征值设置里，要指定搜索区域。我习惯设定搜索频率范围0到1MHz（或者根据换能器中心频率定），并让求解器返回前N个特征频率。这里的N不能太小，低频段模态较少，随着频率升高模态增多，为了覆盖高频段需要把N设到100甚至200。更稳妥的做法是设定频率搜索上界，让求解器自动找出该范围内所有特征值，再用一个参数化扫描来遍历k的步长。
k的取值范围怎么定？最简单的标准是，用你关心的最高频率fmax和最小相速度cmin估算：kmax = 2πfmax/cmin。比如fmax = 1MHz，cmin = 1000m/s，kmax大约6280 rad/m。为了保证曲线平滑，k从0扫到kmax，步长取kmax/150到kmax/300。这个步长在高频段可能略粗，但可以先粗扫一遍看趋势，再在模态密集区细化。
3.4 网格划分：多少个单元才能不骗自己
网格是频散曲线计算里最容易导致结果偏差的一环。核心原则是：厚度方向必须能分辨高阶模态的空间变化。我的做法是每层至少3个四边形单元，推荐4个。四层板就是厚度方向16个单元。x方向因为是Floquet周期边界，平面内波长由k决定，理论上网格只要满足每个波长有10个单元以上即可。
接着按最大k值对应的最短波长校准一下：如果kmax=6280 rad/m，最短波长大约1mm，那么x方向的单元尺寸不应超过0.1mm。但实际二维模型x方向只取1mm，用映射网格划分的话，x方向20个单元就能满足10单元每波长的要求。厚度方向16个单元、x方向20个单元，总共320个单元，自由度大概几千到一万，扫300个k点完全无压力。
网格画完必须先做一次收敛性检查：选一个较大的k值，比如k=4000 rad/m，用当前网格解一次特征频率；再把网格细化一倍，重新解一次。比较前10阶特征频率，如果相对变化都在0.5%以内，说明网格够了。如果高频模态明显漂移，那就是网格太粗，先把厚度方向单元数提上去。
4. 后处理：把一堆特征频率整理成可读的频散曲线
4.1 从COMSOL里批量导出特征频率数据
扫描完几百个k点之后，COMSOL会生成一个很大的表格。最直接的导出方式是在“派生值”里选“全局计算”，表达式填“freq”，然后把参数勾选到所有k值，计算出结果后右键导出为CSV。这时候表格大概有几百行，每行有k值和对应的多个特征频率。
如果你用COMSOL with MATLAB接口，也可以写个简单循环逐k提取，然后一次性合并。这里我更推荐前者，省去脚本时间。导出后的文件格式一般是每步k对应一行，多个特征频率横向排列。接下来就用Python或者MATLAB做规整。
4.2 相速度和群速度的换算
导出数据后，相速度的计算很直接：
cp = ω / k = 2πf / k
注意频率单位是Hz，k单位是rad/m，算出来的cp就是m/s。
群速度稍微绕一点，是角频率对波数的导数：
cg = dω / dk = 2π * df / dk
实际数据处理时，用邻近k点的差分来近似df/dk。我习惯写成如下Python脚本片段：
import numpy as np import pandas as pd data = pd.read_csv('freq_data.csv') freq_cols = [c for c in data.columns if c.startswith('freq')] rows = [] for idx, row in data.iterrows(): k = row['k'] for col in freq_cols: f = row[col] if f > 0: rows.append({'k': k, 'f': f}) df = pd.DataFrame(rows) df['cp'] = 2 * np.pi * df['f'] / df['k'] df_sorted = df.sort_values(['k', 'f']).reset_index(drop=True) df['cg'] = df.groupby('f_col')['f'].transform(lambda x: 2 * np.pi * np.gradient(x, df['k']))实际操作中有一点很容易踩坑：同一个k值下，COMSOL返回的特征频率顺序并不保证和相邻k值下的模态一一对应。模态A在k=100处排第5，在k=101处可能排第7，因为另一个模态插进来了。因此直接差分群速度很可能得到错误值。要先做模态配对，按照“频率随k连续变化”的原则把各条分支理出来，再算群速度。
4.3 曲线绘制与模态标记
画图时我会画两张图：一张是f-k散点图，把所有模态点画在同一张图上；另一张是f-cp散点图。f做横轴，cp做纵轴，一眼就能看出哪条是A0、哪条是S0。对于对称层合板，还能通过查看特征模态的位移云图区分对称和反对称模态，给曲线打上标签。这一步别省，后面阶段还要靠模态标签和实验信号对上号。
如果分支之间有交叉或接近耦合，后处理连线时会乱。我的办法是手动按振型相似度补充识别；自动匹配的算法在模态交换区域很不可靠，宁可保留散点也不强行连线。频散曲线最终呈现给大家的往往是散点图，或者只在趋势很清晰的区域画线。
4.4 验证：拿已知解对照，别自嗨
COMSOL算出来的频散曲线到底靠不靠谱，必须验证。如果你算的是各向同性铝板或者钢板，直接用Rayleigh-Lamb方程自己解一遍A0和S0，跟COMSOL结果对比，前几阶模态误差应该在1%以内。如果误差偏大，大概率是边界条件或网格设置出了问题。
复合材料没有现成解析解，但有经典层合板理论（CLT）可以验证低频极限。比如[A0]模态在频率趋近于零时的波速，接近于由层合板面内等效刚度和密度决定的一个极限值。用CLT算出等效刚度A11和密度ρ，再用cp0 = sqrt(A11 / (ρ * h))做近似对照。如果COMSOL算出的低频相速度方向和参考值差很多，先回头查材料坐标系。
5. 踩坑记录：伪模态、材料方向错配、扫掠步长与网格依赖的完整排查链路
5.1 伪模态：白算了一大堆之后才发现全是假的
特征频率扫描最容易出问题的是伪模态。最常见的一种是在k非常小、接近零的地方，特征频率几乎为零。这些模态对应的其实不是导波，而是刚体平动——整个截面在x方向或y方向没有约束地移动，频率为零。COMSOL特征值求解器不会自动帮你把这些0频率解剔除，它们会混在结果表里，后处理时让低频段看起来有一堆莫名其妙的“曲线”。
检查方法很简单：挑出低频小k的几个特征模态，看位移云图。如果整个截面又平移又转动的，明显不符合导波形态，那就过滤掉。处理办法是把特征值搜索区域的下限从0设为一个很小的正值，比如100Hz，或者在后处理脚本里直接把f<1000Hz的点丢出去。
还有一种伪模态和周期边界有关。当k=0时，Floquet边界退化为普通的连续周期边界，会有一些无关的振动模态混进来，特别是厚度方向均匀拉伸或者压缩的模式。这些模态在频散曲线里通常表现为一条水平直线，看着特别突兀。碰到这种曲线，先看振型确认物理意义，再决定是不是删掉。
5.2 材料坐标系没对准：曲线整体偏了还找不到原因
这是复材仿真最容易犯的低级错误，但排查起来极其折磨。有一次我算一个[0/90/0]对称层合板的频散曲线，算出来低频相速度比文献值低了30%，整体曲线都往低频方向压。折腾了一个晚上，最后发现是90°层材料坐标系旋转方向填反了，x方向的等效刚度算小了很多。
如果材料坐标系的旋转轴和旋转方向设置错误，结果就是整条频散曲线的频率轴偏移、模态顺序错乱。排查路径如下：
- 先建一个只有0°单向层的模型，算低频A0速度，和经典层合板理论的准静态值对比。
- 再建一个只有90°单向层的模型，同样算一次，速度应该明显低于0°单向层。
- 两层结果都对，再回到层合板模型，多半就是层间坐标系定义或者铺层顺序的问题。
用这个办法逐步缩小范围，五分钟就能定位到具体的层。别直接在多层板模型里瞎试，效率太低。
5.3 扫掠步长太大：分支连错线，白画一堆图
频散曲线后处理连线时最坑的是模态配对。COMSOL特征频率求解是按频率从小到大返回的，但这个顺序在每个k点之间并不稳定。高频段、模态密集交叉的区域，频率相近的两个模态很容易在相邻k点之间互换排位，后处理自动连线就容易把A1的一部分接到S1上，曲线跨越交错，图根本没法看。
解决方法是控制k步长。我踩过这个坑之后，一般采取两步走：先用粗步长扫一遍，画出散点图，确定模态密集交叉区的位置；然后在那些区域加密k步长，比如从全局200步加密到局部400步。另外不要依赖自动连线，我现在的习惯是直接用散点图出结果，人工做分支识别，尤其在模态密集区，手动判断比任何自动聚类算法都可靠。
5.4 网格密度和频率上限的相互拉扯
频率越高，模态的空间波长越短，对网格的分辨率要求越高。网格够不够，有一个简单的经验公式：在最高频率fmax下，目标模态的相速度cp_min，取最短波长λmin = cp_min / fmax，网格尺寸必须小于λmin/10。比如fmax=1MHz，cp_min=800m/s，λmin=0.8mm，网格需要小于0.08mm。
有时候你会发现网格加密之后，高频模态的数量突然变多了，之前没算出来的模态全都冒出来。这不是坏事，而是之前网格太粗把真实模态过滤掉了。所以高频段的频散曲线如果出现“莫名其妙少了一段”，先别急着怀疑求解器，回去检查厚度方向的网格分辨率。
另外，特征频率求解的模态数量也要跟着频率上限走。设固定数量的模态数，到高频时可能覆盖不到全部需要的模态；设搜索频率上界，求解器会自动找到该范围内的所有特征值，更稳妥。如果上界内的模态数量太多，求解速度会变慢，那就分频段处理，比如先扫0-500kHz，再扫500kHz-1MHz，最后拼接起来。
总结一下我的工作流
如果从头开始算一条复合材料层合板的频散曲线，我的固定流程是：先用二维截面模型，建好层几何和材料坐标系，输入工程常数，设置Floquet周期条件，特征值研究里把k作为扫描参数，扫300个点左右，映射网格保证厚度方向每层至少3个单元。导出全部特征频率之后，用脚本把频率换算成相速度，手动匹配模态分支，最后和CLT低频极限值做交叉验证。这个流程在COMSOL 6.x各版本里都稳定跑通，耗时一般一两个晚上，比一开始就上三维模型快得多。
最后分享一个小技巧：在COMSOL里做参数化扫描时，最好把扫描参数k本身也作为一个“全局表达式”在结果表里输出。这样导出的CSV每一行都带着对应的k值，后处理直接用pandas按k分组处理，省去自己手工对应参数的麻烦。这个小习惯我第一次注意到是因为来回对参数太累了，养成之后后处理效率提升非常明显。
