Hi All,
I am trying to infer some quantitative pressure values that reaches a transducer from a multilayered solid. To start simple, I am now testing the reflected pressure from a layer interface. However, there is a discrepancy between the theoretical reflection coefficient (0.86) and the ratio of the incident and reflected waves (0.69).
For that, designed a simple circular layout. I emit a tone-burst signal from a circular transducer and record pressures along the vertical radius. So if you check the sensor position 166 >> then the amplitude ratio of incident and reflected pressures gives 0.69. Is this difference expected? or Am I missing something?
Here is the code I am using
% =========================================================================
% MATERIAL PROPERTIES
% =========================================================================
mat1.speed_of_sound = 1482; %[m/s]
mat1.density = 994; %[kg/m3]
mat2.speed_of_sound = 8432; %[m/s]
mat2.density = 2329; %[kg/m3]
% =========================================================================
% SIMULATION GRID PARAMETERS
% =========================================================================
dx = 16e-6;
dy = dx;
Nx = 1024; % [grid points]
Ny = 1024;
kgrid = kWaveGrid(Nx, dx, Ny, dy);
PML_size_x = Nx/16;
PML_size_y = Ny/16;
PML_size = [PML_size_x PML_size_y];
PML_alpha = 3;
% =========================================================================
% Layout Design
% =========================================================================
r_source = 6e-3; % [m]
r_layer1 = 10e-3; % [m]
source_region_mask = ((kgrid.x).^2 + (kgrid.y).^2) < r_source^2;
layer1_mask = ((kgrid.x).^2 + (kgrid.y).^2) < r_layer1^2;
dilated_source_region = imdilate(source_region_mask, strel('disk', 1));
source_ring = ~source_region_mask & dilated_source_region;
% =========================================================================
% MEDIUM PROPERTY ASSIGNMENT
% =========================================================================
medium.sound_speed = ones(Nx,Ny) * mat1.speed_of_sound;
medium.density = ones(Nx,Ny) * mat1.density;
medium.sound_speed(~layer1_mask) = mat2.speed_of_sound;
medium.density(~layer1_mask) = mat2.density;
medium.alpha_coeff = zeros(Nx,Ny);
medium.alpha_power = 1.005;
medium.sound_speed_ref = mat1.speed_of_sound;
% =========================================================================
% SENSOR DEFINITION
% =========================================================================
% define the sensor to record the pressure
sensor.mask = ~ones(Nx, Ny);
sensor.mask(1:Nx/2, Ny/2) = 1;
sensor.record = {'p'};
% =========================================================================
% SOURCE DEFINITION
% =========================================================================
% assign the source
clearvars source;
source.p_mask = source_ring;
% define the driving signal
source_freq = 5e6; % [Hz]
source_strength = 1e6; % [Pa]
source_cycles = 3; % number of tone burst cycles
% create the time steps and the time array
CFL = 0.1; % stability criteria for convergence
t_end = 4e-6; % (tof)
kgrid.makeTime(max(medium.sound_speed(:)), CFL, t_end);
sigma = 0.6/source_freq;
mu = 3 * sigma;
source_sig = -source_strength*exp(-((kgrid.t_array-mu).^2./(2*(sigma^2)))).*sin(2*pi*source_freq*kgrid.t_array);
source.p = source_sig;
% =========================================================================
% ELASTIC SIMULATION
% =========================================================================
% set the input arguments
input_args = { 'PMLSize', PML_size, ...
'PMLAlpha', PML_alpha, ...
'PMLInside', true, ...
'DataCast', 'gpuArray-single', ...
'DataRecast', true, ...
'PlotSim', true, ...
'PlotFreq', 100, ...
'PlotScale', source_strength *[-1 1], ...
'RecordMovie', true, ...
};
% run the simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});
refl_coeff = (mat1.speed_of_sound*mat1.density - mat2.speed_of_sound * mat2.density) / ...
(mat1.speed_of_sound*mat1.density + mat2.speed_of_sound * mat2.density);
sensor_pressure_env = abs(hilbert(plot(sensor_data.p(166,:))));