当前位置：首页 > news >正文

超声波数值模拟与声场传播的MATLAB实现

news 2025/8/18 2:41:05

理论基础

声波传播方程

超声波在介质中的传播可以用波动方程描述：

\(∇²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');