一阶弹性波方程交错网格差分法(时间2阶空间4阶,带边界条件).pdf

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 一阶弹性波方程交错网格差分法(时间2阶空间4阶，带边界条件) % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all; t=0; dt=0.001; snaptime=0.4; X=200; Z=200; dx=10; dz=10; Vp=3000; Vs=1800; D=2000;%密度 %拉梅系数 M=D*Vs^2; L=D*Vp^2-2*D*Vs^2; %差分权系数 c1=9/8; c2=-1/24; U=zeros(X,Z); W=zeros(X,Z); P=zeros(X,Z); Q=zeros(X,Z); S=zeros(X,Z); Up=zeros(X,Z); Wp=zeros(X,Z); Pp=zeros(X,Z); Qp=zeros(X,Z); Sp=zeros(X,Z); Ux=zeros(X,Z);%位移 Uz=zeros(X,Z); for t=0:dt:snaptime %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 应力 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i=3:X-2 for j=3:Z-2 P(i,j)=Pp(i,j)+(dt/dx)*D*Vp*Vp*(1.125*(Up(i+1,j)-Up(i,j))-0(Up(i+2,j)-Up( i-1,j)))+(dt/dz)*D*(Vp*Vp-2.0*Vs*Vs)*(1.125*(Wp(i,j+1)-Wp(i,j))-0(Wp(i,j+ 2)-Wp(i,j-1))); Q(i,j)=Qp(i,j)+(dt/dx)*D*(Vp*Vp-2.0*Vs*Vs)*(1.125*(Up(i+1,j)-Up(i,j))-0(U p(i+2,j)-Up(i-1,j)))+(dt/dz)*D*Vp*Vp*(1.125*(Wp(i,j+1)-Wp(i,j))-0(Wp(i,j+ 2)-Wp(i,j-1))); S(i,j)=Sp(i,j)+(dt/dx)*D*Vs*Vs*(1.125*(Wp(i,j)-Wp(i-1,j))-0(Wp(i+1,j)-Wp( i-2,j)))+(dt/dz)*D*Vs*Vs*(1.125*(Up(i,j)-Up(i,j-1))-0(Up(i,j+1)-Up(i,j-2)) ); end end %四周降阶处理 % i=2; % for j=3:Z-2 % P(i,j)=Pp(i,j)+(dt/dx)*D*Vp*Vp*(1.125*(Up(i+1,j)-Up(i,j))-0(Up(i+2,j)-Up( i-1,j)))+(dt/dz)*D*(Vp*Vp-2.0*Vs*Vs)*(1.125*(Wp(i,j+1)-Wp(i,j))-0(Wp(i,j+ 2)-Wp(i,j-1))); % Q(i,j)=Qp(i,j)+(dt/dx)*D*(Vp*Vp-2.0*Vs*Vs)*(1.125*(Up(i+1,j)-Up(i,j))-0(U p(i+2,j)-Up(i-1,j)))+(dt/dz)*D*Vp*Vp*(1.125*(Wp(i,j+1)-Wp(i,j))-0(Wp(i,j+ 2)-Wp(i,j-1))); % S(i,j)=Sp(i,j)+(dt/dx)*D*Vs*Vs*(Wp(i,j)-Wp(i-1,j))+(dt/dz)*D*Vs*Vs*(1.125*(Up(i,j)- Up(i,j-1))-0(Up(i,j+1)-Up(i,j-2))); % end % % i=X-1; % for j=3:Z-2 % P(i,j)=Pp(i,j)+(dt/dx)*D*Vp*Vp*(Up(i+1,j)-Up(i,j))+(dt/dz)*D*(Vp*Vp-2.0