Hello,
I'm trying to model scattering from a plane shear wave (signal=wavelet) from a hole. My material is steel, and I set the hole material velocity and density to be less than steel. I get several problems:
1. I can make a compression plane wave emanating from the left edge of the grid, by setting source.ux=0 and source.uy=signal. But if I reverse them, to make a shear wave, it gives nonsense results.
2. To simulate a hole (air inside), I set the density very low. However the "hole" seems to generate its own energy. I tried to kill that by setting attenuation high, but the program ground to a halt.
Thanks for any thoughts!
Code:
clear source sensor medium kgrid
% define variables
freq=5e6;
epw=12;
cp=5850;
cs=3230;
rho=8030;
SDHdia=1.5e-3;
%==define simulation area (remember, horizontal is Y direction)==================
xlen=20e-3;
ylen=16e-3;
% ==================== setup grid =========================================
lambda=cp/freq;
dx=lambda/epw;
dy=lambda/epw;
Nx= ceil(xlen/dx);
Ny= ceil(ylen/dy);
kgrid = makeGrid(Nx, dx, Ny, dy);
%====set up material regions=========================================
medium.sound_speed_compression =ones(kgrid.Nx, kgrid.Ny)*cp;
medium.sound_speed_shear =ones(kgrid.Nx, kgrid.Ny)*cs;
medium.alpha_coeff_compression = ones(kgrid.Nx, kgrid.Ny)*.03/(freq/1e6); %to give dB/(MHz^2 cm)
medium.alpha_coeff_shear = ones(kgrid.Nx, kgrid.Ny)*.03/(freq/1e6)*4;
medium.density=ones(kgrid.Nx, kgrid.Ny)*rho;
SDHloc(1)=0;
SDHloc(2)=2e-3;
SDHdia=SDHdia/dx;
[SDHcenterGridx,SDHcenterGridy]=cartPt2Grid(kgrid,[SDHloc(1) SDHloc(2)]);
disc = makeDisc(Nx, Ny, SDHcenterGridx, SDHcenterGridy, SDHdia/2);
medium.sound_speed_compression(disc == 1) = cp*.9;
medium.sound_speed_shear(disc == 1) = cs*.9;
medium.density(disc==1)=60;
medium.alpha_coeff_compression(disc==1) = 1/(freq/1e6); %to give dB/(MHz^2 cm)
medium.alpha_coeff_shear(disc==1) = 4/(freq/1e6)*4;
t_end=.3*2*ylen/cs;
%if ReverseFire==1; t_end=t_end/2;end
dt=.3*dx/(cp); %set time steps (C should be Cmax)
kgrid.t_array = 0:dt:t_end;
Nt=length(kgrid.t_array);
PML=20*dx;
%======================= setup source ===================================================
signal = wavelet(freq,4,0,2/freq,kgrid.t_array);
source.u_mask = zeros(kgrid.Nx,kgrid.Ny); % or source.p_mask
source.u_mask(:,1)=1;
source.ux=signal/(rho*cp);
source.uy=0;
%====setup sensor =====================================================
sensor.mask=zeros(kgrid.Nx, kgrid.Ny);
%sensor.mask(1:5:kgrid.Nx,1:5:kgrid.Ny)=1;
sensor.mask(1:2:kgrid.Nx,1:2:kgrid.Ny)=1;
sensor.record = {'u_max_all','u'}; %{'p','p_max'};
%===== display_mask =========================================
display_mask=disc;
%==================== RUN =================================================
figpos=get(0,'DefaultFigurePosition');
set(0,'DefaultFigurePosition',[340,275,980,735]);
input_args = {'DisplayMask', display_mask, 'PMLInside', true, 'PlotPML', false,'PlotScale','auto','DataCast', 'single'};%,'RecordMovie',true};
sensor_data = pstdElastic2D(kgrid, medium, source, sensor, input_args{:});
%sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});
set(0,'DefaultFigurePosition',figpos);
kmodel=sensor_data;
clear sensor_data;
kmodel.freq=freq;
kmodel.Nx=Nx;kmodel.Ny=Ny;
kmodel.x_vec=kgrid.x_vec;kmodel.y_vec=kgrid.y_vec;
kmodel.Nt=length(kgrid.t_array);
kmodel.epw=epw;
save kmodel.mat kmodel -v7.3
uall=sqrt(kmodel.ux_max_all.^2+kmodel.uy_max_all.^2);
u=sign(kmodel.ux.*kmodel.uy).*sign(kmodel.ux).*sqrt(kmodel.ux.^2+kmodel.uy.^2);
%uh=abs(hilbert(u'));uh=uh';
%nox=ceil(kmodel.Nx/5);noy=ceil(kmodel.Ny/5);
nox=ceil(kmodel.Nx/2);noy=ceil(kmodel.Ny/2);
dat=zeros(nox,noy);%may need to reverse this is Nx, Ny are in wrong order
for t=1:size(u,2);
for j=1:nox
for k=1:noy;
%dat(j,k,t)=uh((k-1)*nox+j,t);
dat(j,k,t)=u((k-1)*nox+j,t);
%daty(j,k,t)=kmodel.uy((k-1)*nox+j,t);
end
end
end
% Record the movie
%maxdat=max(max(max(dat)));
%ax=[0 maxdat/500];
nframes=50;
step=floor((kmodel.Nt-1)/nframes);
figure;imagesc(kmodel.y_vec,kmodel.x_vec,(dat(:,:,200)));axis equal;%caxis(ax)
set(gca,'nextplot','replacechildren');
ax=caxis;
for j= 1:nframes;
t=1+(j-1)*step;
imagesc(kmodel.y_vec,kmodel.x_vec,(dat(:,:,t)));caxis(ax)
F(j) = getframe;
end
% Play the movie
movie(F,1,1)