Hello, ask for a question:
1. I want to simulate the sound wave from the first medium obliquely incident into the second medium, reflect and refract at the interface of the two layers of media, then the sound wave re-enters the third medium from the second medium, and produces reflection at the interface between the two media. And refraction. I tried to program, but it didn't work. I don't know if I can simulate three layers of media using this software, or there are other ways to solve my problem.
I am looking forward to the expert's answer, thank you very much.
2.Regarding the PML problem, I have not understood thoroughly whether one of the four boundaries of the square region can be set to be non-PML, so that when the sound wave is incident on the bottom surface of the medium, it is reflected back instead of being absorbed.
I look forward to your answer, thank you very much.
k-Wave
A MATLAB toolbox for the time-domain
simulation of acoustic wave fields
Ultrasonic three-layer medium propagation simulation
(4 posts) (2 voices)-
Posted 6 years ago #
-
Hi zhangwei198202,
(1) Yes this is possible. Take a look at the examples included in k-Wave, in particular, the "Snell's Law And Critical Angle Reflection Example".
(2) It is not currently possible to turn off just one PML. In any case, k-Wave uses a Fourier collocation spectral method, so turning off the PML results in wave wrapping rather than reflections.
Brad.
Posted 6 years ago # -
Hello, Bradley Treeby, thank you for your advice. I carefully refer to the example "Shear Waves And Critical Angle Reflection Example" in the software, because this example is very similar to the model I want to do. But my ability is limited, I can't fully understand the program. I added a layer of air medium to the program, placed it on the bottom as the third layer medium, and set all the parameters of the air to zero. Although the program can run, I know that it is definitely not correct. I sent the program to you. Please help me check it.
What I don't understand is that the difference in acoustic impedance between the second medium and the third layer is obvious, why the reflected wave is not visible in the second medium.
I hope that you can take the time to give your answer, thank you very much!
Here is the program:clearvars;
% =========================================================================
% SIMULATION PARAMETERS
% =========================================================================% change scale to 2 to reproduce higher resolution figures in help file
scale = 1;% create the computational grid
PML_size = 10; % [grid points]
Nx = 128*scale - 2*PML_size; % [grid points]
Ny = 192*scale - 2*PML_size; % [grid points]
dx = 0.5e-3/scale; % [m]
dy = 0.5e-3/scale; % [m]
kgrid = kWaveGrid(Nx, dx, Ny, dy);% define the medium properties for the top layer
cp1 = 1540; % compressional wave speed [m/s]
cs1 = 0; % shear wave speed [m/s]
rho1 = 1000; % density [kg/m^3]
alpha0_p1 = 0.1; % compressional absorption [dB/(MHz^2 cm)]
alpha0_s1 = 0.1; % shear absorption [dB/(MHz^2 cm)]% define the medium properties for the bottom layer
cp2 = 3000; % compressional wave speed [m/s]
cs2 = 1400; % shear wave speed [m/s]
rho2 = 1850; % density [kg/m^3]
alpha0_p2 = 1; % compressional absorption [dB/(MHz^2 cm)]
alpha0_s2 = 1; % shear absorption [dB/(MHz^2 cm)]% define the third layer medium properties for the bottom layer
cp3 = 0; % compressional wave speed [m/s]
cs3 = 0; % shear wave speed [m/s]
rho3 = 0; % density [kg/m^3]
alpha0_p3 = 0; % compressional absorption [dB/(MHz^2 cm)]
alpha0_s3 = 0; % shear absorption [dB/(MHz^2 cm)]% create the time array
cfl = 0.1;
t_end = 60e-6;
kgrid.makeTime(cp1, cfl, t_end);% define position of heterogeneous slab
slab = zeros(Nx, Ny);
slab(Nx/3:2/3*Nx, :) = 1;% define position of heterogeneous slab
air = zeros(Nx, Ny);
air(2/3*Nx:end, :) = 2;% define the source geometry in SI units (where 0, 0 is the grid center)
arc_pos = [-15e-3, -25e-3]; % [m]
focus_pos = [5e-3, 5e-3]; % [m]
radius = 25e-3; % [m]
diameter = 20e-3; % [m]% define the driving signal
source_freq = 500e3; % [Hz]
source_strength = 1e6; % [Pa]
source_cycles = 3; % number of tone burst cycles% define the sensor to record the maximum particle velocity everywhere
sensor.record = {'u_max_all'};% set the input arguments
input_args = {'PMLSize', PML_size, 'PMLAlpha', 2, 'PlotPML', false, ...
'PMLInside', false, 'PlotScale', [-1, 1]*source_strength, ...
'DisplayMask', 'off', 'DataCast', 'single'};% =========================================================================
% FLUID SIMULATION
% =========================================================================% assign the medium properties
medium.sound_speed = cp1*ones(Nx, Ny);
medium.sound_speed(slab == 1) = cp2;
medium.density = rho1*ones(Nx, Ny);
medium.density(slab == 1) = rho2;
medium.alpha_coeff = alpha0_p1*ones(Nx, Ny);
medium.alpha_coeff(slab == 1) = alpha0_p2;
medium.alpha_power = 2;% convert the source parameters to grid points
arc_pos = round(arc_pos / dx) + [Nx/2, Ny/2];
focus_pos = round(focus_pos / dx) + [Nx/2, Ny/2];
radius = round(radius / dx);
diameter = round(diameter / dx);% force the diameter to be odd
if ~rem(diameter, 2)
diameter = diameter + 1;
end% generate the source geometry
source_mask = makeArc([Nx, Ny], arc_pos, radius, diameter, focus_pos);% assign the source
source.p_mask = source_mask;
source.p = source_strength*toneBurst(1/kgrid.dt, source_freq, source_cycles);% run the fluid simulation
sensor_data_fluid = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});% =========================================================================
% ELASTIC SIMULATION
% =========================================================================% define the medium properties
clear medium
medium.sound_speed_compression = cp1*ones(Nx, Ny);
medium.sound_speed_compression(slab == 1) = cp2;
medium.sound_speed_compression(air == 1) = cp3;medium.sound_speed_shear = cs1*ones(Nx, Ny);
medium.sound_speed_shear(slab == 1) = cs2;
medium.sound_speed_shear(air == 1) = cs3;medium.density = rho1*ones(Nx, Ny);
medium.density(slab == 1) = rho2;
medium.density(air == 1) = rho3;medium.alpha_coeff_compression = alpha0_p1*ones(Nx, Ny);
medium.alpha_coeff_compression(slab == 1) = alpha0_p2;
medium.alpha_coeff_compression(air == 1) = alpha0_p3;medium.alpha_coeff_shear = alpha0_s1*ones(Nx, Ny);
medium.alpha_coeff_shear(slab == 1) = alpha0_s2;
medium.alpha_coeff_shear(air == 1) = alpha0_s3;% assign the source
clear source
source.s_mask = source_mask;
source.sxx = -source_strength*toneBurst(1/kgrid.dt, source_freq, source_cycles);
source.syy = source.sxx;% run the elastic simulation
sensor_data_elastic = pstdElastic2D(kgrid, medium, source, sensor, input_args{:});% =========================================================================
% VISUALISATION
% =========================================================================% define plot vector
x_vec = kgrid.x_vec * 1e3;
y_vec = kgrid.y_vec * 1e3;% calculate square of velocity magnitude
u_e = sensor_data_elastic.ux_max_all.^2 + sensor_data_elastic.uy_max_all.^2;
u_f = sensor_data_fluid.ux_max_all.^2 + sensor_data_fluid.uy_max_all.^2;% plot layout
figure;
imagesc(y_vec, x_vec, double(source_mask | slab));xlabel('y [mm]');
ylabel('x [mm]');
axis image;
colormap(flipud(gray));% plot beam patterns for the fluid simulation
figure;
subplot(2, 1, 1);
imagesc(y_vec, x_vec, 20*log10(u_f./max(u_f(:))));
xlabel('y [mm]');
ylabel('x [mm]');
axis image;
h = colorbar;
caxis([-50, 0]);
title('Fluid Model');
ylabel(h, '[dB]');% plot beam patterns for the elastic simulation
subplot(2, 1, 2);
imagesc(y_vec, x_vec, 20*log10(u_e./max(u_e(:))));
xlabel('y [mm]');
ylabel('x [mm]');
axis image;
h = colorbar;
caxis([-50, 0]);
title('Elastic Model');
colormap(jet(256));
ylabel(h, '[dB]');scaleFig(1, 1.5);
Posted 6 years ago # -
Hi zhangwei198202,
You set:
air(2/3*Nx:end, :) = 2;But later use:
medium.sound_speed_compression(air == 1) = cp3;(etc)You should use either both 1 or both 2.
Brad
Posted 6 years ago #
Reply
You must log in to post.