Hi,
I am trying to run simulations to model the ultrasound propagation of a focused transducer in 3D.
My issue is that my source creates non-physical edge waves at the edges of the source. (My source is modeled as a bowl.) The edge waves also appear in a linear simulation and in 2D simulations. However, in 2D simulations, the effect is much weaker. I don't believe that the problem is a discretization issue because I have not perceived any improvements when reducing the grid spacing.
Is there a solution to this?
I have also attached the code for more detailed information.
Thank you!
clearvars; % Clears all variables
D = 64e-3; % Aperture diameter [m]
C = 98e-3; % Radius of Curvature [m]
c = 1482; % Speed of Sound [m/s]
f = 1.1e6; % Frequency [Hz]
mag = 3e5; % Magnitude [Pa]
rho = 998.2; % Density [kg/m^3]
cfl = 0.3;
% Domain
x_dim = 80e-3; % [m]
y_dim = 160e-3; % [m]
z_dim = 80e-3;
% Size of the domain (grid points)
Nx = 864;
Ny = 2048;
Nz = 864;
% =========================================================================
% SIMULATION
% =========================================================================
% create the computational grid
dx = x_dim/Nx; % grid point spacing in the x direction [m]
dy = y_dim/Ny; % grid point spacing in the y direction [m]
dz = z_dim/Nz; % grid point spacing in the z direction [m]
% Create the grid
kgrid = kWaveGrid(Nx, dx, Ny, dy,Nz,dz);
% define the properties of the propagation medium
medium.sound_speed = c ; % [m/s]
medium.density = rho ; % [kg/m^3]
medium.alpha_coeff =2.17e-3; % [dB/(MHz^y cm)] for water
medium.alpha_power = 2; %water
medium.BonA = 5.0; % solve nonlinear equations
distanceTraveledAtTEnd = 160e-3;
tEnd = distanceTraveledAtTEnd/ min(medium.sound_speed(:));
kgrid.makeTime(medium.sound_speed,cfl,tEnd);
grid_size = [Nx,Ny,Nz];
bowl_pos = [Nx/2,50,Nz/2]; % position of the center of the bowl%50
radius = round(C/dx); % in grid spaces
diameter = 2*floor(D/(2*dx))+1; % in grid spaces
focus_pos = [Nx/2,400,Nz/2]; % position of the focus
bowl_1 = makeBowl(grid_size, bowl_pos, radius, diameter, focus_pos, 'Plot', false);
source.p_mask = bowl_1;
% define a time varying sinusoidal pressure burst
source_freq = f; % [Hz]
source_mag = mag;
number_cycles = 4;
sinePressure = source_mag * sin(2 * pi * source_freq * kgrid.t_array);
% ensure smooth transition from mag=0 to mag=source_mag and the reverse
% model one pressure burst via Tukey window
burst_duration = number_cycles/source_freq;
t_mask = kgrid.t_array(kgrid.t_array <= burst_duration);
num_samples_in_burst = length(t_mask);
window = [transpose(tukeywin(num_samples_in_burst, 0.5)), ...
zeros(1, length(kgrid.t_array) - num_samples_in_burst)];
source.p = sinePressure .* window; % Apply the Tukey window to the pressure signal
%filter the source to remove any high frequencies not supported by the grid
source.p = filterTimeSeries(kgrid, medium, source.p);
% create sensor - line along transducer axis
% sensor.mask = [Nx/2, 70, Nz/2, Nx/2,1999,Nz/2].';
sensor.mask = zeros(Nx, Ny,Nz);
sensor.mask(1:2:Nx, 1:2:799,1:2:Nz) = 1;
sensor.record = {'p_final'};%, 'p_max', 'p_min'};
% run the simulation
input_args = {'PMLSize', 30, 'PMLInside',true, 'SaveToDisk', 'input'};
sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});