Hi!
I am trying to model acoustic scattering in air from rigid surface. I've tried couple of methods before with small luck. I tried to set almost 0 sound speed for scattering surface which was fine and fast but there was significant phase shift. Second approach was to use density, sound speed and alpha coeff of real material. Works nice, but this way time step have to be quite small to avoid "explosion". And while I need to make thousands of evaluations this scheme is too slow.
Finally I discovered idea to set dirichlet conditions for velocity to 0. This scheme works fast and seems to give realistic enough scattering with correct phase.
But after some tests I discovered 2 strange things:
1) After some simulation time passes dirichlet conditions go away and waves just propagate from that moment as if u_mask was zero matrix.
2) Some strange vibrations occur on the left side of scattering surface. Strange thing that this effect is not symmetrical as there is nothing like that on the right side of the sample.
Here is my code:
% create the computational grid
Nx = 160; % number of grid points in the x (row) direction
Ny = 160; % number of grid points in the y (column) direction
dx = 0.025; % grid point spacing in the x direction [m]
dy = 0.025; % grid point spacing in the y direction [m]
kgrid = makeGrid(Nx, dx, Ny, dy);
% define the properties of the propagation medium
medium.sound_speed = 345*ones(Nx,Ny); % [m/s] Air
medium.alpha_coeff = 1.1795*ones(Nx,Ny); % [dB/(MHz^y cm)]
medium.alpha_power = 1.4484; %Air (?)
medium.density = 1.1864; %kg/m3 Air
PML_x=20;
PML_y=20;
% create the time array
dt = 1.5e-5;
total_time = 0.05;
kgrid.t_array = 0:dt:total_time;
% define source point(s)
source.p_mask = zeros(Nx, Ny);
source.p_mask(Nx/2,Ny/2) = 1;
% define Ricker source
fcent = 1500;
t0 = 0.002;
source_amp = 0.18e-2;
source.p = source_amp*-sqrt(2)*pi*fcent*((sqrt(2)*pi*fcent*(kgrid.t_array-t0)).^2 - 1) ...
.*exp(-0.5*(sqrt(2)*pi*fcent*(kgrid.t_array-t0)).^2);
% Filter
source_notfiltered = source.p;
medium_filter.sound_speed=345;
source.p = filterTimeSeries(kgrid, medium_filter, source.p);
source.u_mask = zeros(Nx, Ny);
source.u_mode='dirichlet';
source.u_mask(PML_x+1:PML_x+8,PML_y+1:Ny-PML_y-1)=1;
source.ux(1:nnz(source.u_mask))=0;
source.uy(1:nnz(source.u_mask))=0;
sensor.mask = [0;0];
sensor.record = {'p_rms'};
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, 'DataCast', 'gpuArray-single', ...
'PlotScale', [-0.1, 0.1], 'PlotSim', true, 'DataRecast', true, 'LogScale', true, 'PMLSize',[PML_x,PML_y]);