I have the problem with 3D simulation of Lamb wave in plate wich is embended in artificial fluid (which should simulate air).
My simulations were based on existing in help code for simulations of shearwave and Snells law.
I have performed simulations in 2D and the results are very good and in agreement with theoretical (analitical) model. Then I tried to perform such anaysis in 3D and in fact simulations is working but the results are bad (I receive NAN in the values of UX, uy, uz, p). Below please find well working 2D code and the 3D code I have developed based on 2D code. I would appreciate for any sugestions
Kind Regards Michal
1. 2D code - GOOD
clear all;
clc
close all
% =========================================================================
% SIMULATION PARAMETERS
% =========================================================================
scale = 1;
% create the computational grid
PML_size = 20; % size of the PML in grid points
Nx = 59*scale - 2*PML_size; % number of grid points in the x direction
Ny = 300*scale - 2*PML_size; % number of grid points in the y direction
% Nz = 300*scale - 2*PML_size; % number of grid points in the y direction
dx = 0.5e-3/scale; % grid point spacing in the x direction [m]
dy = 0.5e-3/scale; % grid point spacing in the y direction [m]
% dz=dx;
kgrid = makeGrid(Nx, dx, Ny, dy);
% kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz);
% define the medium properties
cp1 = 600; % compressional wave speed [m/s]
cs1 = 0; % shear wave speed [m/s]
rho1 = 200; % density [kg/m^3]
alpha0_p1 = 5.0; % compressional wave absorption [dB/(MHz^2 cm)]
alpha0_s1 = 5.0; % shear wave absorption [dB/(MHz^2 cm)]
cp2 = 2656; % compressional wave speed [m/s]
cs2 = 1078; % shear wave speed [m/s]
rho2 = 1140; % density [kg/m^3]
alpha0_p2 = 0.5; % compressional wave absorption [dB/(MHz^2 cm)]
alpha0_s2 = 0.5; % shear wave absorption [dB/(MHz^2 cm)]
% create the time array
cfl = 0.1;
t_end = 200e-6;
kgrid.t_array= makeTime(kgrid, cp1, cfl, t_end);
fs=1/kgrid.dt;
% define position of heterogeneous slab
% slab = ones(Nx, Ny, Nz);
% slab (1,:,:)=0 ;
% slab (19,:,:)=0 ;
slab = ones(Nx, Ny);
slab (1,:)=0 ;
slab (19,:)=0 ;
% define the source properties
source_freq = 0.05e6; % [Hz]
source_strength = 2e6; % [Pa]
source_cycles = 2; % number of tone burst cycles
%3D
% source_mask=zeros(Nx,Ny,Nz);
% source_mask(2,30,6)=1;
% source_mask(18,30,6)=1;
%2D
source_mask=zeros(Nx,Ny);
source_mask(2,6)=1;
source_mask(18,6)=1;
% assign the source
source.s_mask = source_mask;
source.sxx(1,:) = source_strength*toneBurst(1/kgrid.dt, source_freq, source_cycles);
source.sxx(2,:) = -source_strength*toneBurst(1/kgrid.dt, source_freq, source_cycles);
% define the sensor mask
sensor.mask = zeros(Nx, Ny);
for k=171:1:221
sensor.mask(2,k:k+10) = 1; %2D
end
sensor.record = {'p','u'};
% set the input arguments
input_args = {'PMLSize', PML_size, 'PMLAlpha', 2, 'PlotPML', false, ...
'PMLInside', false, 'PlotScale', [-1, 1]*source_strength/2, ...
'DisplayMask', 'off', 'DataCast', 'single'};
% =========================================================================
% ELASTIC SIMULATION
% =========================================================================
% define the medium properties
medium.sound_speed_compression = cp1*ones(Nx, Ny);
medium.sound_speed_compression(slab == 1) = cp2;
medium.sound_speed_shear = cs1*ones(Nx, Ny);
medium.sound_speed_shear(slab == 1) = cs2;
medium.density = rho1*ones(Nx, Ny);
medium.density(slab == 1) = rho2;
medium.alpha_coeff_compression = alpha0_p1*ones(Nx, Ny);
medium.alpha_coeff_compression(slab == 1) = alpha0_p2;
medium.alpha_coeff_shear = alpha0_s1*ones(Nx, Ny);
medium.alpha_coeff_shear(slab == 1) = alpha0_s2;
% run the elastic simulation
sensor_data_elastic = pstdElastic2D(kgrid, medium, source, sensor, input_args{:});
p=sensor_data_elastic.p;
subplot(2,2,1),imagesc(p),title('Preassure')
ux=sensor_data_elastic.ux;
subplot(2,2,2),imagesc(ux),title('Ux')
uy=sensor_data_elastic.uy;
subplot(2,2,3),imagesc(uy),title('Uy')
save('SIM_Results2D.mat','p','ux','uy','fs','source_freq')
2. 3D code - BAD
clear all;
clc
close all
% =========================================================================
% SIMULATION PARAMETERS
% =========================================================================
scale = 1;
% create the computational grid
PML_size = 20; % size of the PML in grid points
Nx = 59*scale - 2*PML_size; % number of grid points in the x direction
Ny = 100*scale - 2*PML_size; % number of grid points in the y direction
Nz = 300*scale - 2*PML_size; % number of grid points in the y direction
dx = 0.5e-3/scale; % grid point spacing in the x direction [m]
dy = 0.5e-3/scale; % grid point spacing in the y direction [m]
dz=dx;
kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz);
% define the medium properties
cp1 = 600; % compressional wave speed [m/s]
cs1 = 0; % shear wave speed [m/s]
rho1 = 200; % density [kg/m^3]
alpha0_p1 = 5.0; % compressional wave absorption [dB/(MHz^2 cm)]
alpha0_s1 = 5.0; % shear wave absorption [dB/(MHz^2 cm)]
cp2 = 2656; % compressional wave speed [m/s]
cs2 = 1078; % shear wave speed [m/s]
rho2 = 1140; % density [kg/m^3]
alpha0_p2 = 0.5; % compressional wave absorption [dB/(MHz^2 cm)]
alpha0_s2 = 0.5; % shear wave absorption [dB/(MHz^2 cm)]
% create the time array
cfl = 0.1;
t_end = 100e-6;
kgrid.t_array= makeTime(kgrid, cp1, cfl, t_end);
fs=1/kgrid.dt;
% define position of heterogeneous slab
slab = ones(Nx, Ny, Nz);
slab (1,:,:)=0 ;
slab (19,:,:)=0 ;
% define the source properties
source_freq = 0.05e6; % [Hz]
source_strength = 2e6; % [Pa]
source_cycles = 2; % number of tone burst cycles
source_mask=zeros(Nx,Ny,Nz);
source_mask(2,30,6)=1;
source_mask(18,30,6)=1;
% assign the source
source.s_mask = source_mask;
source.sxx(1,:) = source_strength*toneBurst(1/kgrid.dt, source_freq, source_cycles);
source.sxx(2,:) = -source_strength*toneBurst(1/kgrid.dt, source_freq, source_cycles);
% define the sensor mask
sensor.mask = zeros(Nx, Ny, Nz);
for k=171:1:221
sensor.mask(2,30,k:k+10) = 1;
end
sensor.record = {'p','u'};
% set the input arguments
input_args = {'PMLSize', PML_size, 'PMLAlpha', 2, 'PlotPML', false, ...
'PMLInside', false, 'PlotScale', [-1, 1]*source_strength/2, ...
'DisplayMask', 'off', 'DataCast', 'single'};
% =========================================================================
% ELASTIC SIMULATION
% =========================================================================
% define the medium properties
clear medium
medium.sound_speed_compression = cp1*ones(Nx, Ny, Nz);
medium.sound_speed_compression(slab == 1) = cp2;
medium.sound_speed_shear = cs1*ones(Nx, Ny, Nz);
medium.sound_speed_shear(slab == 1) = cs2;
medium.density = rho1*ones(Nx, Ny, Nz);
medium.density(slab == 1) = rho2;
medium.alpha_coeff_compression = alpha0_p1*ones(Nx, Ny, Nz);
medium.alpha_coeff_compression(slab == 1) = alpha0_p2;
medium.alpha_coeff_shear = alpha0_s1*ones(Nx, Ny, Nz);
medium.alpha_coeff_shear(slab == 1) = alpha0_s2;
% run the elastic simulation
sensor_data_elastic = pstdElastic3D(kgrid, medium, source, sensor, input_args{:});
p=sensor_data_elastic.p;
subplot(2,2,1),imagesc(p),title('Preassure')
ux=sensor_data_elastic.ux;
subplot(2,2,2),imagesc(ux),title('Ux')
uy=sensor_data_elastic.uy;
subplot(2,2,3),imagesc(uy),title('Uy')
save('SIM_Results3D.mat','p','ux','uy','fs','source_freq')