Hi,
I am trying to simulate the pressure distribution in a viscoelastic medium using pstdElastic3D. When I am using a velocity source with a rectangular source mask and getting exceptionally high-pressure values (P_max). I am not sure where I am going wrong. Please help. Code is pasted below.
clear; close all; clc;
%Define the medium and grids (Sampling)
L_x=(20*10^-3);% Unit:m
L_y=(20*10^-3);% Unit:m
L_z=(20*10^-3);% Unit:m
C_L=1540; % compressional sound speed; unit: m/s
Frequency_US=1*10^6; %unit: Hz (Maximum frequency)
Wavelength_Min=C_L/Frequency_US; % Unit: m
dx=Wavelength_Min/5;
dy=dx;
dz=dx;
Nx=ceil(L_x/dx); % Number of grids in the x direction; ceil rounds the number to next integer.
Ny=ceil(L_y/dy); % Number of grids in the y direction
Nz=ceil(L_z/dz); % Number of grids in the z direction
kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz); % kWavegrid class makes and store grid points
%Assigning the temporal grids
CFL=0.5;
dt=CFL*dx/C_L;% CFL=c*dt/dx
simulation_time=40*10^-6; % unit: seconds
Nt=ceil(simulation_time/dt);% round-off to the next integer
kgrid.setTime(Nt, dt);% it sets the time grid and makes time array
%medium properties for pstdElastic3D
medium.sound_speed_compression=C_L; %m/s
medium.sound_speed_shear=1.7; %Unit: m/s
medium.density=980; %in kg/m^3
medium.alpha_coeff_compression=0.1; %dB/MHz^2-cm
medium.alpha_coeff_shear=0.3; %dB/MHz^2-cm
%Defining the source input
% define properties of the input signal
source_strength = 1e-6; % [m/s]
tone_burst_freq = Frequency_US; % [Hz]
tone_burst_cycles = 5;
input_signal = source_strength*toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles);
Transducer_height=46;% in grid points
transducer_width = 32;
source.u_mask = zeros(Nx, Ny, Nz);
source.u_mask(11, ceil(Ny/2 - transducer_width/2) + 1:ceil(Ny/2 + transducer_width/2), ceil(Nz/2 - Transducer_height/2) + 1:ceil(Nz/2 + Transducer_height/2)) = 1
source.u_mode='dirichlet';
focus_position=[0, 0, 0];
focus_signal = focus(kgrid, input_signal, source.u_mask, focus_position, C_L);
% plot(focuse_signal((93*0)+101,:))
source.ux=focus_signal; %focused_signal;
%% define the sensor mask
sensor.mask = zeros(Nx, Ny, Nz);
%taking the sensor mask at the X-Y plane at the center of the medium (z=0
%plane)
sensor.mask(:,:,round(Nz/2)) = 1;
sensor.record = {'p', 'p_max', 'u', 'I', 'I_avg'};
sensor.record_start_index=1; %(default=1)
%set input values for PML boundary layer
alpha_coeff_PML=2.5; % dB/MHz-cm (default=2)
PML_size=[10,10,10]; %(default=10)
% set the optional input as input arguements
input_args = {'PlotPML', false, 'PMLAlpha', alpha_coeff_PML, 'PMLSize', PML_size};
% RUN the elastic simulation
%use pstdElastic3D function to simulate propagation of elastic wave in this
%homogenous viscoelastic medium
sensor_data=pstdElastic3D(kgrid, medium, source, sensor, input_args{:});