Hi, thanks for providing the great library.
I would like to simulate a small 40kHz-ultrasound transducer with a flat reflector in 3D.
When I set the medium heterogeneous (air and acrylic), I cannot see any propagation even in the air. Is this because of the large change in impedance? I would appreciate it if you could tell me how to make my code work.
Thanks advance!
Best regards,
Ryuji
% =========================================================================
% sorce code
% =========================================================================
clear all;
% =========================================================================
% DEFINE LITERALS
% =========================================================================
% select which code to run
% 1: MATLAB CPU code
% 2: MATLAB GPU code
% 3: C++ code
MODEL = 2;
% scale parameter (changes the resolution of the simulation)
SC = 1;
% grid parameters
PML_SIZE = 10; % size of the perfectly matched layer [grid points]
Nx = 128*SC - 2*PML_SIZE; % number of grid points in the x direction
Ny = 64*SC - 2*PML_SIZE; % number of grid points in the y direction
Nz = 64*SC - 2*PML_SIZE; % number of grid points in the z direction
% sampling parameters
PPW = 10*SC; % points per wavelength (this defines the grid size)
PPP = 10*SC; % points per period (this defines the CFL)
T_END = 1e-3; % total compute time [s] (this must be long enough to reach steady state)
RECORD_PERIODS = 2; % number of periods to record
% medium parameters
C0 = 346; % sound speed [m/s]
RHO0 = 1.2; % density [kg/m^3]
CR = 1.43e3; % sound speed in the reflector [m/s]
RHOR = 1.18e3; % density of the reflector [kg/m^3]
% source parameters
F0 = 40e3; % source frequency [Hz]
SOURCE_RADIUS = 5e-3; % piston radius [m]
SOURCE_MAG = 0.5e6; % source pressure [Pa]
% =========================================================================
% CREATE GRID
% =========================================================================
% calculate the grid spacing based on the PPW and F0
dx = C0 / (PPW * F0); % [m]
dt = 1 / (PPP*F0); % [s]
SOURCE_RADIUS = SOURCE_RADIUS / dx;
% calculate the CFL
disp(['CFL = ' num2str(C0*dt/dx)]);
% create the computational grid
kgrid = makeGrid(Nx, dx, Ny, dx, Nz, dx);
% create the time array
kgrid.t_array = 0:dt:T_END;
% assign medium properties
% medium.sound_speed = C0;
% medium.density = RHO0;
medium.sound_speed = C0 * ones(Nx, Ny, Nz);
medium.density = RHO0 * ones(Nx, Ny, Nz);
% assign reflector properties
medium.sound_speed(Nx, :, :) = CR;
medium.density(Nx, :, :) = RHOR;
% assign the reference sound speed to the background medium
% medium.sound_speed_ref = C0;
% =========================================================================
% CREATE SOURCE
% =========================================================================
% create piston source mask
piston = makeDisc(Ny, Nz, Ny/2, Nz/2, SOURCE_RADIUS);
source.p_mask = zeros(Nx, Ny, Nz);
source.p_mask(1, :, :) = piston;
% create time varying source
source.p = SOURCE_MAG * sin(2*pi*F0*kgrid.t_array);
% filter source with up-ramp
ramp_length = round((2*pi/F0)/dt);
ramp = (-cos( (0:(ramp_length-1))*pi/ramp_length ) + 1)/2;
source.p(1:ramp_length) = ramp .* source.p(1:ramp_length);
% =========================================================================
% CREATE SENSOR MASK
% =========================================================================
% set sensor mask to record central axis, not including the source point
% sensor.mask = zeros(Nx, Ny, Nz);
% sensor.mask(2:end, Ny/2, Nz/2) = 1;
sensor.mask = [2, 1, Nz/2, Nx, Ny, Nz/2].';
% record the maximum pressure
sensor.record = {'p', 'p_rms'};
% average only the final few periods when the field is in steady state
sensor.record_start_index = kgrid.Nt - RECORD_PERIODS*PPP + 1;
% =========================================================================
% RUN k-WAVE SIMULATION
% =========================================================================
% run code
switch MODEL
case 1
% MATLAB CPU
sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, ...
'DataCast', 'single', 'PMLInside', false, ...
'DisplayMask', source.p_mask, 'PlotScale', [-1, 1]*SOURCE_MAG);
case 2
% MATLAB GPU
sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, ...
'DataCast', 'gpuArray-single', 'PMLInside', false, ...
'DisplayMask', source.p_mask, 'PlotScale', [-1, 1]*SOURCE_MAG);
case 3
% C++
sensor_data = kspaceFirstOrder3DC(kgrid, medium, source, sensor, 'PMLInside', false);
otherwise
end
% =========================================================================
% ANALYTICAL SOLUTION
% =========================================================================
% calculate the wavenumber
k = 2*pi*F0./C0;
% define radius and axis
a = SOURCE_RADIUS*dx;
x = (kgrid.x_vec(2:end, :) - min(kgrid.x_vec(:)));
% calculate the analytical solution for a piston in an infinite baffle
% for comparison (Eq 5-7.3 in Pierce)
r_a = sqrt(x.^2 + a^2);
p_ref = SOURCE_MAG*abs(-2i*exp(1i*k*(x+r_a)/2).*sin((k*r_a - k*x)/2));
% =========================================================================
% VISUALISATION
% =========================================================================
% plot the pressure along the focal axis of the piston
figure;
plot(sensor_data.p(end,:), 'bx');
img = gather(sensor_data.p_rms);
image(img, 'CDataMapping', 'scaled')
daspect([1 1 1])
colormap default
colorbar