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)