Hi, Brad. Thanks for the last reply. It help a lot.
But I also confused about the attenuation is Elastic2d Model.
I am trying to simulate the wave propagation in sandstone, and I set the compression alpha coefficient as 2.7[dB/MHZ^2 cm], 5.8[dB/MHZ^2 cm] for shear. But it seems the system is not stable. I don't understand how that happened.
Thanks!!
clc
clear
close all
% =========================================================================
% DEFINE SIMULATION PROPERTIES
% =========================================================================
% define the properties used in the simulation
c0 = 3760; % compressional sound speed [m/s]
shear_c0 = 2300; % shear sound speed [m/s]
rho0 = 2650; % density [kg/m^3]
frac_c0 = 1481; % open fracture compression sound speed [m/s], water
fracshear_c0 = 1; % open fracture shear sound speed [m/s], water
frac_rho0 = 1000; % open fracture density [kg/m^3], water
points_per_wavelength = 30; % number of grid points per wavelength at f0
wavelength_separation = 10; % separation between the source and detector(farest)
pml_size = 20; % PML size
pml_alpha = 2; % PML alpha
% compute corresponding grid spacing
dx = 0.0001; % grid point spacing in the x direction [m]
dy = dx; % grid point spacing in the y direction [m]
% compute corresponding grid size
Nx = wavelength_separation*points_per_wavelength + 20; % number of grid points in t he x (row) direction
Ny = wavelength_separation*points_per_wavelength + 20; % number of grid points in the y (column) direction
material_length = (wavelength_separation*points_per_wavelength*dx)*1e3; % [mm]
grid_length = 3*dx*1e3;% [mm]fracture is three grids thick
% =========================================================================
% SIMULATION
% =========================================================================
% create the computational grid
kgrid = makeGrid(Nx, dx, Ny, dy);
% define a single sensor position an integer number of wavelengths away
sensor.mask = zeros(Nx, Ny);
y_loc = (10:60:wavelength_separation*points_per_wavelength + 10);
x_loc = (10+60:60:wavelength_separation*points_per_wavelength + 10 - 60);
sensor.mask(10, y_loc) = 1;
sensor.mask(Nx-10, y_loc) = 1;
sensor.mask(x_loc, 10) = 1;
sensor.mask(x_loc, Ny-10) = 1;
% define the crack locations
start_x = randi([10 Nx-10],1,1);
start_y = randi([10 Nx-10],1,1);
dist = randi([50,200],1,1);
angle_degree = randi([1 180],1,1);
end_x = round(start_x + cos(angle_degree*3.1415926/180).* dist );
end_y = round(start_y + sin(angle_degree*3.1415926/180).* dist );
while end_x <=10 || end_x >= Nx-10 || end_y <= 10 || end_y >= Nx-10
start_x = randi([10 Nx-10],1,1);
start_y = randi([10 Nx-10],1,1);
dist = randi([30,200],1,1);
angle_degree = randi([1 180],1,1);
end_x = round(start_x + cos(angle_degree*3.1415926/180).* dist );
end_y = round(start_y + sin(angle_degree*3.1415926/180).* dist );
end
[xx_new,yy_new] = XiaolinWu(start_x-1, start_y, end_x-1, end_y);
[x_new,y_new] = XiaolinWu(start_x, start_y, end_x, end_y);
[xxx_new,yyy_new] = XiaolinWu(start_x+1, start_y, end_x+1, end_y);
% define the absorption properties
medium.alpha_coeff_compression = 2.7; % [dB/(MHz^2 cm)]??????
medium.alpha_coeff_shear = 5.4; % [dB/(MHz^2 cm)]??????
% define the properties of the propagation medium
medium.sound_speed_compression = c0*ones(Nx, Ny); % [m/s]
medium.sound_speed_shear = shear_c0*ones(Nx, Ny); % [m/s]
medium.density = rho0*ones(Nx, Ny); % [kg/m^3]
for j = 1: length(x_new)
medium.sound_speed_compression(x_new(j), y_new(j)) = frac_c0; % [m/s]
medium.sound_speed_shear(x_new(j), y_new(j)) = fracshear_c0; % [m/s]
medium.density(x_new(j), y_new(j)) = frac_rho0;
medium.sound_speed_compression(xx_new(j), yy_new(j)) = frac_c0; % [m/s]
medium.sound_speed_shear(xx_new(j), yy_new(j)) = fracshear_c0; % [m/s]
medium.density(xx_new(j), yy_new(j)) = frac_rho0;
medium.sound_speed_compression(xxx_new(j), yyy_new(j)) = frac_c0; % [m/s]
medium.sound_speed_shear(xxx_new(j), yyy_new(j)) = fracshear_c0; % [m/s]
medium.density(xxx_new(j), yyy_new(j)) = frac_rho0;
end
figure(1)
imagesc(medium.density);
hold on
scatter (10, Nx/2, 'r.')
hold on
plot (y_loc, 10, 'k.')
hold on
plot (10, x_loc, 'k.')
hold on
plot (y_loc, 310, 'k.')
hold on
plot (310, x_loc, 'k.')
% set end time
t_end_cal = (material_length/3760)*10^6;
t_end = 15e-6; % end at 50us
% create the time array
kgrid.t_array = makeTime(kgrid, max(medium.sound_speed_compression(:)) ,[], t_end);
% define the initial pressure distribution source term
disc_magnitude = 1e5; % [Pa]...
disc_x_pos = Nx/2; % [grid points]
disc_y_pos = 10; % [grid points]
disc_radius = 1; % [grid points]
source.p0 = disc_magnitude*makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);
source.p0 = smooth(kgrid, source.p0, true);
% set the simulation options
input_args = {'PlotFreq', 20, 'PlotScale', 'auto', 'PlotLayout', false,...
'PMLInside', false, 'PMLSize', pml_size, 'PMLAlpha', pml_alpha, 'PlotPML', true};
% define the sensor to record the maximum particle velocity everywhere
sensor.record = {'p'};
% run the simulation
sensor_data = pstdElastic2D(kgrid, medium, source, sensor, input_args{:});
% =========================================================================
% VISUALISATION
% =========================================================================
%%
[t_sc, scale, prefix] = scaleSI(max(kgrid.t_array(:)));
figure(4)
for n = 1 : 6
subplot(6, 1, n), plot(kgrid.t_array*scale, sensor_data.p(n,:), 'r-');
xlabel(['Time [' prefix 's]']);
ylabel('Signal Amplitude-Pressure');
axis tight;
% ylim([-max(sensor_data(n,:)) max(sensor_data(n,:))])
end
%%
figure(5)
for n = 7 : 14
subplot(8, 1, n-6), plot(kgrid.t_array*scale, sensor_data.p(n,:), 'r-');
xlabel(['Time [' prefix 's]']);
ylabel('Signal Amplitude-Pressure');
axis tight;
% ylim([-max(sensor_data(n,:)) max(sensor_data(n,:))])
end
%%
figure(6)
for n = 15 : 20
subplot(6, 1, n-14), plot(kgrid.t_array*scale, sensor_data.p(n,:), 'r-');
xlabel(['Time [' prefix 's]']);
ylabel('Signal Amplitude-Pressure');
axis tight;
% ylim([-max(sensor_data(n,:)) max(sensor_data(n,:))])
end