hi everyone，

I'm doing simulation about a HIFU transducer, but I found that when I set the grids' size smaller, the pressure of focus doesn't converge. But when I set the medium linear(delete medium.BonA=7), the results show convergency. Can anyone help me to solve this problem?

Here are my codes:

Nx = 380; % number of grid points in the x direction

Ny = 240; % number of grid points in the y direction

Nz = 240; % 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]媒介密度

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;

sensor.record = {'p_max'};

% % input arguments

input_args = {'PlotPML', true, ...

'DataCast', 'single', 'CartInterp', 'linear','PlotSim',true,'PlotLayout',true,'PMLInside',true,'PMLSize',[20,20,20]};

% 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;

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

end

%ONeil solution

% 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]');

title('Maximum Pressure');

c = colorbar;

c.Label.String = 'pressure()';