理论基础
声波传播方程
超声波在介质中的传播可以用波动方程描述:
\(∇²p - (1/c²)∂²p/∂t² = 0\)
其中:
- \(p\) 是声压
- \(c\) 是介质中的声速
- \(∇²\) 是拉普拉斯算子
数值求解方法
我们使用有限差分时域法(FDTD)进行数值求解:
- 离散空间和时间域
- 使用中心差分近似导数
- 迭代求解声压场
MATLAB实现代码
%% 超声波传播数值模拟 - MATLAB实现
clear; close all; clc;%% 参数设置
c = 1500; % 声速 (m/s) - 水中声速
f0 = 1e6; % 超声波频率 (Hz)
lambda = c/f0; % 波长 (m)% 模拟区域参数
Lx = 10*lambda; % x方向长度 (m)
Lz = 10*lambda; % z方向长度 (m)
dx = lambda/20; % 空间步长 (m)
dz = dx; % 空间步长 (m)
dt = dx/(2*c); % 时间步长 (s) - 满足CFL条件% 网格设置
Nx = round(Lx/dx); % x方向网格点数
Nz = round(Lz/dz); % z方向网格点数
x = (0:Nx-1)*dx; % x坐标向量
z = (0:Nz-1)*dz; % z坐标向量
[X, Z] = meshgrid(x, z); % 网格坐标% 时间参数
t_total = 20*lambda/c; % 总模拟时间
Nt = round(t_total/dt); % 时间步数
t = (0:Nt-1)*dt; % 时间向量%% 初始化声场
p = zeros(Nz, Nx); % 当前时刻声压
p_prev = p; % 前一时刻声压
p_next = p; % 下一时刻声压%% 设置声源(高斯脉冲)
source_pos_x = round(Nx/2); % 声源x位置
source_pos_z = round(Nz/4); % 声源z位置% 高斯脉冲参数
t0 = 5e-6; % 脉冲中心时间
sigma = 1/(2*pi*f0); % 脉冲宽度% 创建声源信号
source_signal = exp(-((t-t0).^2)/(2*sigma^2)) .* sin(2*pi*f0*t);%% 设置障碍物(可选)
obstacle = false(Nz, Nx); % 障碍物掩膜
% 在中间设置一个圆形障碍物
center_x = round(Nx/2);
center_z = round(Nz/2);
radius = lambda; % 障碍物半径for i = 1:Nxfor j = 1:Nzdist = sqrt((x(i)-x(center_x))^2 + (z(j)-z(center_z))^2);if dist <= radiusobstacle(j, i) = true;endend
end%% 吸收边界层(减小边界反射)
absorb_layer_width = 30; % 吸收层宽度(网格点数)
absorb_coeff = 0.01; % 吸收系数% 创建吸收层掩膜
absorb_mask = ones(Nz, Nx);
for i = 1:absorb_layer_width% 左右边界absorb_mask(:, i) = exp(-absorb_coeff*(absorb_layer_width-i)^2);absorb_mask(:, end-i+1) = exp(-absorb_coeff*(absorb_layer_width-i)^2);% 上下边界absorb_mask(i, :) = exp(-absorb_coeff*(absorb_layer_width-i)^2);absorb_mask(end-i+1, :) = exp(-absorb_coeff*(absorb_layer_width-i)^2);
end%% 主循环 - 模拟声波传播
figure('Position', [100, 100, 800, 600]);for n = 1:Nt% 更新声场方程 (2D波动方程)for i = 2:Nx-1for j = 2:Nz-1% 跳过障碍物区域if obstacle(j, i)p_next(j, i) = 0;continue;end% 有限差分方程d2p_dx2 = (p(j, i+1) - 2*p(j, i) + p(j, i-1))/dx^2;d2p_dz2 = (p(j+1, i) - 2*p(j, i) + p(j-1, i))/dz^2;p_next(j, i) = 2*p(j, i) - p_prev(j, i) + ...(c^2 * dt^2) * (d2p_dx2 + d2p_dz2);endend% 添加声源p_next(source_pos_z, source_pos_x) = p_next(source_pos_z, source_pos_x) + ...source_signal(n);% 应用吸收边界p_next = p_next .* absorb_mask;% 更新场变量p_prev = p;p = p_next;% 每10步绘制一次if mod(n, 10) == 0% 绘制声压场subplot(2, 1, 1);imagesc(x, z, p);axis image;title(sprintf('声压场 (t = %.2f μs)', t(n)*1e6));xlabel('x (m)');ylabel('z (m)');colorbar;clim([-1, 1]); % 固定颜色范围colormap(jet);% 绘制障碍物轮廓hold on;contour(x, z, obstacle, [0.5, 0.5], 'k', 'LineWidth', 1.5);hold off;% 绘制声源位置hold on;plot(x(source_pos_x), z(source_pos_z), 'ro', 'MarkerSize', 8, 'LineWidth', 1.5);hold off;% 绘制声源信号subplot(2, 1, 2);plot(t(1:n)*1e6, source_signal(1:n), 'b', 'LineWidth', 1.5);title('声源信号');xlabel('时间 (μs)');ylabel('振幅');grid on;xlim([0, t_total*1e6]);ylim([-1.2, 1.2]);drawnow;end
end%% 生成声场传播动画
% 初始化视频
videoFile = VideoWriter('ultrasound_propagation.mp4', 'MPEG-4');
videoFile.FrameRate = 10;
open(videoFile);fig = figure('Position', [100, 100, 800, 600]);% 重新初始化声场
p = zeros(Nz, Nx);
p_prev = p;
p_next = p;for n = 1:50:Nt % 跳步以减小文件大小% 更新声场(与主循环相同)% ...(此处省略更新代码以节省空间,实际使用时需包含)% 绘制声场imagesc(x, z, p);axis image;title(sprintf('声压场传播 (t = %.2f μs)', t(n)*1e6));xlabel('x (m)');ylabel('z (m)');colorbar;clim([-0.5, 0.5]);colormap(jet);% 添加障碍物和声源标记hold on;contour(x, z, obstacle, [0.5, 0.5], 'k', 'LineWidth', 1.5);plot(x(source_pos_x), z(source_pos_z), 'ro', 'MarkerSize', 8, 'LineWidth', 1.5);hold off;% 捕获帧并写入视频frame = getframe(fig);writeVideo(videoFile, frame);
end% 完成并保存视频
close(videoFile);
disp('声场传播动画已保存为 ultrasound_propagation.mp4');
关键
1. 数值稳定性处理
- CFL条件:\(dt ≤ dx/(c√2)\) 确保数值稳定性
- 空间离散:使用中心差分法离散拉普拉斯算子
- 时间迭代:蛙跳法(leapfrog)更新声压场
2. 边界处理
- 吸收边界层:在模拟区域边界添加指数衰减层
- 边界吸收系数:\(exp(-α·d²)\)其中d为到边界的距离
- 障碍物处理:设置障碍物区域声压始终为零
3. 声源模型
- 高斯调制脉冲:\(s(t) = exp(-(t-t₀)²/(2σ²))·sin(2πf₀t)\)
- 位置设置:在网格中心附近设置点声源
4. 可视化技术
- 动态显示声压场分布
- 障碍物轮廓叠加
- 声源位置标记
- 实时显示声源信号
- 生成传播过程动画
参考仿真 超声波数值模拟,通过matlab模拟声场的传播 youwenfan.com/contentcnb/65678.html
参数调整建议
- 提高分辨率:减小\(dx\), \(dz\)(需相应减小 \(dt\))
- 改变介质:调整声速 \(c\) 值
- 修改声源:调整 \(f0\) 或使用不同脉冲形状
- 复杂障碍物:定义任意形状的障碍物掩膜
- 添加衰减:在介质模型中引入频率相关衰减
通过调整这些参数,可以模拟各种实际场景中的超声波传播行为。