Hi，
I'm doing bowl focused transducer simulation, with signal of 1Mhz and 543.36e3 pa. But when I use smaller dx, dy, dz ,the pressure of focus(max pressure) is getting higher. And when I use 0.3mm of grid size, the max pressure is about 30 Mpa, and 0.2mm of grid size, the max pressure is about 40 Mpa! While the same parameters I put in focusedBowlONeil() to calculate, the max pressure at focus is just about 20 Mpa, almost the same 20 Mpa results when I use 0.5mm of grid size(dx,dy,dz). Am I use so large amplitude?
My English isn't very well , hope you understand my confusion and help me to solve this problem.
Thanks!!!
Here are my codes:
clear;
clf;
close all;
addpath('D:\k-wave-toolbox-version-1.4\k-Wave');
Nx = 130*3; % number of grid points in the x direction
Ny = 80*3; % number of grid points in the y direction
Nz = 80*3; % number of grid points in the z direction
dx =0.3e-3; % grid point spacing in the x direction [m]
dy =0.3e-3; % grid point spacing in the y direction [m]
dz =0.3e-3; % grid point spacing in the z direction [m]
x_size=dx*Nx;
y_size=dy*Ny;
z_size=dz*Nz;
sound_speed=1500; % [m/s]媒介声速
density = 1000 ; % [kg/m^3]媒介密度
% radius = 63*2; %换能器曲率半径
% diameter = 63*2 - 1; %换能器直径（奇数）
kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz);
% define the properties of the propagation medium
medium.BonA=7; %水里的超声非线性参数
medium.sound_speed = sound_speed * ones(Nx, Ny, Nz); % [m/s]
medium.density = density * ones(Nx, Ny, Nz); % [kg/m^3]
kgrid.makeTime(medium.sound_speed,0.3,80e-6);
% define a centered circular sensor
% define a series of Cartesian points to collect the data
sensor.mask = zeros(Nx, Ny, Nz);
sensor.mask(:, :, Nz/2) = 1;
% set the record mode such that only the max and peak values are stored
sensor.record = {'p_max'};
% % input arguments
input_args = {'PlotPML', true, ...
'DataCast', 'single', 'CartInterp', 'linear','PlotSim',true,'PlotLayout',true,'PMLInside',true,'PMLSize',[20,20,20]};
% grid_size = [Nx, Ny, Nz];%网格空间
% bowl_pos = [1, Ny/2,Nz/2];%bowl的中心点
% %radius = Nx/2;%曲率半径
% %diameter =Nx-1;%圆的直径
% focus_pos = [Nx/2, Ny/2, Nz/2 ];%与bowl_pos连线为对称轴（确定方向） 和法向量有一点点像 决定换能器姿态
% source.p_mask=makeBowl(grid_size, bowl_pos, radius, diameter, focus_pos, 'Plot', true);
% define a time varying sinusoidal source
source_mag = 543.36e3; % [Pa]
% create empty array
karray = kWaveArray('SinglePrecision',true);
%position = [1, Ny/2, Nz/2];
position = [-x_size/2+20e-3,0,0];
% define arc properties
radius = 63.2e-3; % [m]
diameter = 64e-3; % [m]
focus_pos = [x_size/2, 0, 0];
karray.addBowlElement(position, radius, diameter, focus_pos);
source.p_mask = karray.getArrayBinaryMask(kgrid);
voxelPlot(double(source.p_mask));
source_freq = 1e6; % [Hz]声源频率
sig1 = source_mag * sin(2 * pi * source_freq * kgrid.t_array);
% smooth the source
sig1 = filterTimeSeries(kgrid, medium, sig1,'PlotSignals',true,'PlotSpectrums',true);
source_signal = zeros(1,length(sig1));
source_signal(1,1:length(sig1)) = sig1;
source.p = karray.getDistributedSourceSignal(kgrid, source_signal);
sensor_data = kspaceFirstOrder3DG(kgrid, medium, source, sensor, input_args{:});
sensor_data.p_max = reshape(sensor_data.p_max, [Nx, Ny]);
sensor_data.p_max = reshape(sensor_data.p_max, [Nx, Ny]);
figure;
imagesc(kgrid.x_vec, kgrid.y_vec,sensor_data.p_max);
xlabel('x [m]');
ylabel('y [m]');
%zlabel('z [m]');
title('Maximum Pressure');
c = colorbar;
c.Label.String = 'pressure()';
max = sensor_data.p_max(1,1);
maxi = [1,1];
[max, maxi];
for a=floor(0.3*Nx):Nx
for b=1:Ny
if sensor_data.p_max(a,b)>max
max= sensor_data.p_max(a,b);
maxi=[a,b];
end
end
end
focus_x = 0;
focusx_max = 0;
focusy_max = 0;
focus_y = 0;
focus_l = 0;
focus_r = 0;
focus_u = 0;
focus_d = 0;
[focus_x,focus_y];
for a=floor(0.3*Nx):(Nx - 1)
for b=1:(Ny - 1)
if sensor_data.p_max(a,b)<0.71*max&&sensor_data.p_max(a,b+1)>0.71*max
%sensor_data.p_max(a,b)=0;
focus_l = b+1;
end
if sensor_data.p_max(a,b)>0.71*max&&sensor_data.p_max(a,b+1)<0.71*max
%sensor_data.p_max(a,b)=0;
focus_r = b+1;
end
if sensor_data.p_max(a,b)<0.71*max&&sensor_data.p_max(a+1,b)>0.71*max
focus_u = a+1;
end
if sensor_data.p_max(a,b)>0.71*max&&sensor_data.p_max(a+1,b)<0.71*max
%sensor_data.p_max(a,b)=0;
focus_d = a+1;
end
end
focus_x = focus_r - focus_l;
focus_y = focus_d - focus_u;
if focus_x > focusx_max
focusx_max = focus_x;
end
if focus_y > focusy_max
focusy_max = focus_y;
end
% if focus_d > 0
% break;
% end
end
% define transducer parameters
radius = 63.2e-3; % [m]
diameter = 64e-3; % [m]
frequency = 1e6; % [Hz]
sound_speed = 1500; % [m/s]
density = 1000; % [kg/m^3]
Z = density*sound_speed;
p = 543.36e3;
velocity = p/Z; % [m/s]
% define position vectors
axial_position = 0:5e-4:150e-3; % [m]
lateral_position = -15e-3:5e-4:15e-3; % [m]
% evaluate pressure
[p_axial, p_lateral] = focusedBowlONeil(radius, diameter, velocity, ...
frequency, sound_speed, density, axial_position, lateral_position);
% plot
figure;
subplot(4, 1, 1);
plot(axial_position .* 1e3, p_axial .* 1e-6, 'k-');
xlabel('Axial Position [mm]');
ylabel('Pressure [MPa]');
subplot(4, 1, 2);
plot((0:dx:(Nx-1)*dx) .* 1e3, sensor_data.p_max(:,Ny/2).* 1e-6, 'k-');
xlabel('Axial Position [mm]');
ylabel('Pressure [MPa]');
subplot(4, 1, 3);
plot(lateral_position .* 1e3, p_lateral .* 1e-6, 'k-');
xlabel('Lateral Position [mm]');
ylabel('Pressure [MPa]');
subplot(4, 1, 4);
plot((0:dy:(Ny-1)*dy).* sensor_data.p_max(Ny/2,:) .* 1e-6, 'k-');
xlabel('Lateral Position [mm]');
ylabel('Pressure [MPa]');
figure;
imagesc(kgrid.x_vec, kgrid.y_vec,sensor_data.p_max);
xlabel('x [m]');
ylabel('y [m]');
%zlabel('z [m]');
title('Maximum Pressure');
c = colorbar;
c.Label.String = 'pressure()';