clear all
close all
clc
Nx1=400;
Nz1=400; % 差分网格400*400
fm=20; % 子波主频20Hz
space_x=0.25;
space_z=0.25; % 网格步长0.25*0.25m
sample_t=0.1/1000; % 时间步长0.1ms
T=4000*0.1/1000; % 模拟时间400ms
K=T/sample_t; % 时间网格点数
layer=40*space_x; % 吸收层厚度40
Nx=Nx1+2*layer/space_x;
Nz=Nz1+layer/space_x; % 加吸收层后的差分网格
coefficientR=0.0001; % 理论反射系数
K1=50;
V1=zeros(2000,K1);
V2=zeros(2000,Nz,Nx);
V0=zeros(Nz,Nx);
ux=zeros(2000,Nz,Nx);
uz=zeros(2000,Nz,Nx);
%% 建立均匀半空间地质模型 %%
VP=zeros(Nx,Nz);
VS=zeros(Nx,Nz);
Density=zeros(Nx,Nz);
lame1=zeros(Nx,Nz);
lame2=zeros(Nx,Nz);
%% 计算区域 %%
for i=1:Nx1
for j=1:Nz1
VP(i+layer/space_x,j)=600;
VS(i+layer/space_x,j)=300;
Density(i+layer/space_x,j)=1.8*1000;
end
end
%% 吸收层 %%
for i=layer/space_x+1:Nx-layer/space_x
for j=Nz-layer/space_x+1:Nz
VP(i,j)=VP(i,Nz-layer/space_x);
VS(i,j)=VS(i,Nz-layer/space_x);
Density(i,j)=Density(i,Nz-layer/space_x);
end
end
for i=1:layer/space_x
for j=1:Nz
VP(i,j)=VP(1+layer/space_x,j);
VS(i,j)=VS(1+layer/space_x,j);
Density(i,j)=Density(1+layer/space_x,j);
end
end
for i=Nx-layer/space_x+1:Nx
for j=1:Nz
VP(i,j)=VP(Nx-layer/space_x,j);
VS(i,j)=VS(Nx-layer/space_x,j);
Density(i,j)=Density(Nx-layer/space_x,j);
end
end
%% 吸收层 %%
for i=1:Nx
for j=1:Nz
lame1(i,j)=Density(i,j)*(VP(i,j)^2-2*VS(i,j)^2);
lame2(i,j)=Density(i,j)*VS(i,j)^2;
end
end
%%%%%空间4阶精度差分系数%%%%%%%
C1=1.125;
C2=-0.0416667;
%%%%%初始条件%%%%%
for i=1:Nx
for j=1:Nz
U_x(i,j)=0;
V_x(i,j)=0;
P_x(i,j)=0;
Q_x(i,j)=0;
R_x(i,j)=0;
U_y(i,j)=0;
V_y(i,j)=0;
P_y(i,j)=0;
Q_y(i,j)=0;
R_y(i,j)=0;
U(i,j)=0;
V0(j,i)=0;
V(i,j)=0;
P(i,j)=0;
Q(i,j)=0;
R(i,j)=0;
vy(i,j)=0;
vyx(i,j)=0;
vyz(i,j)=0;
txy(i,j)=0;
tzy(i,j)=0;
dx(i)=0;
dy(j)=0;
end
end
% vy=zeros(Nx,Nz);
% vyx=zeros(Nx,Nz);
% vyz=zeros(Nx,Nz);
% txy=zeros(Nx,Nz);
% tzy=zeros(Nx,Nz);
%%%%%AEA自由表面边界条件%%%%%%%%%%%%%%%
for j=1:4
for i=1:Nx
Density(i,j)=0.5*Density(i,j);
lame1(i,j)=0;
lame2(i,j)=lame2(i,j);
end
end
%%%%%%%%%%波动方程差分计算%%%%%%%%%%%%%%%%
for k=1:K
mibinbin=k
%%%%AEA设置%%%%
for i=1:Nx
Q(i,4)=0;
end
%%%%计算更新各点速度%%%%%%
for i=4:(Nx-4)
for j=4:(Nz-4)
if i<=4+layer/space_x
dx(i)=(2*VP(i,j)/layer)*log(1/coefficientR)*((4+layer/space_x-i)*space_x/layer)^4;
end
if i>=Nx-4-layer/space_x
dx(i)=(2*VP(i,j)/layer)*log(1/coefficientR)*((i-(Nx-4-layer/space_x))*space_x/layer)^4;
end
if j>=Nz-4-layer/space_x
dy(j)=(2*VP(i,j)/layer)*log(1/coefficientR)*((j-(Nz-4-layer/space_x))*space_x/layer)^4;
end
vyx(i,j)=(vyx(i,j)*(1-sample_t*dx(i)/2)+sample_t/Density(i,j)*((C1*(txy(i,j)-txy(i-1,j))+C2*(txy(i+1,j)-txy(i-2,j)))/space_x))/(1+sample_t*dx(i)/2);
vyz(i,j)=(vyz(i,j)*(1-sample_t*dy(j)/2)+sample_t/Density(i,j)*((C1*(tzy(i,j)-tzy(i,j-1))+C2*(tzy(i,j+1)-tzy(i,j-2)))/space_z))/(1+sample_t*dy(j)/2);
vy(i,j)=vyx(i,j)+vyz(i,j);
end
end
%%%%计算更新各点应力%%%%%%%
for i=4:(Nx-4)
for j=4:(Nz-4)
if i<=4+layer/space_x
dx(i)=(2*VP(i,j)/layer)*log(1/coefficientR)*((4+layer/space_x-i)*space_x/layer)^4;
end
if i>=Nx-4-layer/space_x
dx(i)=(2*VP(i,j)/layer)*log(1/coefficientR)*((i-(Nx-4-layer/space_x))*space_x/layer)^4;
end
if j>=Nz-4-layer/space_x
dy(j)=(2*VP(i,j)/layer)*log(1/coefficientR)*((j-(Nz-4-layer/space_x))*space_x/layer)^4;
end
%% 震源设置 %%
vy(200+layer/space_x,4)=100*(1-2*(pi*fm*(k-500)*sample_t)^2)*exp(-(pi*fm*(k-500)*sample_t)^2);
txy(i,j)=(lame2(i,j)*sample_t*((C1*(vy(i+1,j)-vy(i,j))+C2*(vy(i+2,j)-vy(i-1,j)))/space_x)+txy(i,j)*(1-sample_t*dx(i)/2))/(1+sample_t*dx(i)/2);
tzy(i,j)=(lame2(i,j)*sample_t*((C1*(vy(i,j+1)-vy(i,j))+C2*(vy(i,j+2)-vy(i,j-1)))/space_z)+tzy(i,j)*(1-sample_t*dy(j)/2))/(1+sample_t*dy(j)/2);
end
end
% for i=1:Nx
% for j=1:Nz
% V0(j,i)=V(i,j);
% end
% end
% V2(k,:,:)=V0(:,:);
% for i=1:K1
% V1(k,i)=V(1+layer/space_x+1*(i-1),4); % 最小偏移距为7m，道间距为0.5m
% end
% figure(2)
% imagesc(V0)
% axis square
end
% figure(1)
% wigb(V1)
% axis([0,51,0,200])
% xlabel('道数')
% ylabel('t/ms')
评论0