MATLAB 编程计算lamb波频散曲线。 有限元算lamb波频散曲线 代码可以得到lamb波的频散曲线和群速度曲线。 完整MATLAB程序。 可运行。 有限元可以得到频散曲线 相速度曲线 #频散曲线 #MATLAB程序 不同要求可议价

最近在搞超声导波检测，发现Lamb波的频散曲线计算是个绕不过去的坎。这玩意儿在板状结构损伤识别里太重要了，直接决定探头频率选多少合适。今天咱们用MATLAB撸个有限元计算的脚本，手把手教你怎么出频散曲线。

先上硬货——全套代码扔这儿了，后面慢慢拆解：

% 材料参数设置
E = 70e9;         % 弹性模量 (Pa)
nu = 0.33;        % 泊松比
rho = 2700;       % 密度 (kg/m^3)
thickness = 0.003;% 板厚 (m)

% 频率扫描设置
freq_range = linspace(1e3, 500e3, 50);  % 1kHz到500kHz取50个点
c_phase = zeros(length(freq_range), 4); % 存储前四个模态的相速度
c_group = zeros(size(c_phase));         % 群速度矩阵

% 有限元模型初始化
model = createpde('structural','modal-solid');
geometryFromEdges(model, @squareg);  % 创建平板几何
structuralProperties(model,'YoungsModulus',E, ...
                        'PoissonsRatio',nu,...
                        'MassDensity',rho);
generateMesh(model,'Hmax',thickness/5); % 网格尺寸不大于板厚1/5

这段代码设定了铝板的材料参数，频率扫描范围从1kHz到500kHz取50个点。有限元部分用到了MATLAB自带的PDE工具箱，这里注意网格尺寸要足够小，不然高频段的结果会飘。一般来说网格尺寸取波长的1/6~1/8，这里直接按板厚的五分之一来划分，够用了。

接下来是核心计算部分：

for i = 1:length(freq_range)
    omega = 2*pi*freq_range(i);
    % 求解特征频率问题
    RF = solve(model,'FrequencyRange',[0.8 1.2]*omega/(2*pi));
    
    % 提取有效模态
    valid_modes = find(RF.NaturalFrequencies >= 0.9*omega/(2*pi) & ...
                       RF.NaturalFrequencies <= 1.1*omega/(2*pi));
    
    % 计算相速度和群速度
    for j = 1:min(length(valid_modes),4)
        k = sqrt(rho*(omega^2)*thickness^2/(E*(1-nu^2))); % 波数估计
        c_phase(i,j) = omega/k; 
        % 群速度数值微分计算
        if i > 1
            c_group(i,j) = (c_phase(i,j)-c_phase(i-1,j))/(freq_range(i)-freq_range(i-1));
        end
    end
end

这里采用扫频法逐个频率点求解。有个坑要注意——solve函数返回的模态可能包含非传播模态，所以用valid_modes筛选出频率附近的真实Lamb波模态。群速度计算用了最朴素的向前差分，实际工程中建议用中心差分提高精度。

最后画图部分来个双Y轴显示：

figure('Position',[200 200 800 400])
yyaxis left
plot(freq_range/1e3,c_phase(:,1:2),'LineWidth',1.5)
ylabel('Phase Velocity (m/s)')

yyaxis right
plot(freq_range/1e3,c_group(:,1:2),'--','LineWidth',1.5)
ylabel('Group Velocity (m/s)')

xlabel('Frequency (kHz)')
title('Lamb Wave Dispersion Curves')
legend('S0 Phase','A0 Phase','S0 Group','A0 Group')
grid on
set(gca,'FontSize',12)

跑出来的图形应该能看到典型的S0和A0模态曲线。注意看高频段S0模态的群速度会趋向于材料纵波速度，这个可以当作验证计算正确性的一个check点。

几个常见翻车点提醒：

1. 网格太粗导致高频段数据跳变——试着把Hmax参数调小
2. 出现异常模态——检查valid_modes的筛选条件
3. 群速度曲线毛刺太多——改用更精细的频率步长或者样条插值

这个代码框架已经可以处理各向同性材料的对称/反对称模态计算。如果需要算更高阶模态或者各向异性材料，改改材料参数输入和模态筛选条件就行。代码里故意留了几个可调参数（比如频率点数、模态数量限制），方便各位根据自己需求魔改。