声波加pml边界条件

雷克子波因计算机资源有限，不便在太大的区域做PML的计算，一般只在边界上取有限宽度的区域作为PML计算区域。根据相关文献的研究，PML边界区域最少长度应为半个波长[6]。本文综合考虑了效果与开销等因素，选取了边界上50层作为PML的计算区域。常规计算区域与PML边界区域的如图2-4所示。xxxxxPvaPKtxzzzzzPvaPKtzxzxxvvPaPKtxz衰减系数2331log2xzPvRL式中L为PML层的厚度，为层内的点距PML与非PML的边界的距离，Pv为纵波速度，R那么在PML边界区域内，对于式（2-13）1xxxvPvtx即为理论反射系数，一般取0.001较为合适，为方向的空间步长。1xxxvPvtx，可看作为在常规的计算方程基础上，减去一项进行PML的阻尼修正项。因本文中只考虑各项同性介质中的地震波传播规律，故可做xz假设。在此利用一下三个假设：2121,2121,21,21kjixkjixkjixvvv21,2121,21,2121kjizkjizkjizvvv1,1,21,21kjikjikjiPPP因为以上三个近似精度均为时间方向上的近似，且时间精度均为二阶精度，因交错网格技术的时间精度为二阶，故以上近似不影响本式的计算精度。故可得：NmkmjikmjimkjikjixkjixPPaxtvttv11,,21,2121,2121,15.015.011（2-18a）同理：NmkjmikjmimkjikjizkjizPPaztvttv1,1,21,21,2121,2115.015.011（2-18b）NmkmjizkmjixmNmkjmizkjmizmkjikjikjivvaxtvvaztKPttP12121,2121,121,2121,2121,,1,5.015.011（2-18c）此为在PML边界区域内的弹性声波应力-速度方程组。%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%交错网格---非均匀介质二维声波方程(一阶压力--速度)、2阶时间差分、2阶空间差分精度%%加上边界%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%closeall;clear,clctic%%***********************震源为Ricker子波*********dtt=0.0001;tt=-0.06:dtt:0.06;fm=30;A=0.01;wave=A*(1-2*(pi*fm*tt).^2).*exp(-(pi*fm*tt).^2);plot(wave),title('震源子波--Ricker子波');%%***********************************************%%模型参数设置dz=5;%纵向网格大小，单位mdx=5;%横向网格大小，单位mdt=0.0001;%时间步长，单位sT=0.5;%波动传播时间，单位swave(round(T/dt))=0;%将子波后面部分补零%%%研究区域%z=-750:dz:750;x=-1000:dz:1000;pml=50;%吸收层的网格数plx=pml*dx;%上下吸收层的厚度plz=pml*dz;%左右吸收层的厚度z=-750-plz:dz:750+plz;x=-1000-plx:dx:1000+plx;%采样区间n=length(z);m=length(x);%采样点数z0=round(n/2);x0=round(m/2);%震源位置Vmax=0;%纵波最大速度%%SettingVelocity&Densityzt=-750-plz:dz/2:750+plz;xt=-1000-plx:dx/2:1000+plx;%速度与密度采样区间nt=length(zt);mt=length(xt);%速度与密度采样点数目V=zeros(n,m);%介质速度,m/sd=zeros(nt,mt);%介质密度,kg/m^3%%均匀介质模型fori=1:nfork=1:mV(i,k)=2.0e3;endendfori=1:nfork=1:md(2*i,2*k)=2.3e3;endend%%%%层状介质模型%%fori=1:n%%fork=1:m%%ifiround(n/3)%%V(i,k)=2.3e3;%%else%%V(i,k)=3.0e3;%%end%%end%%endfori=1:n-1fork=1:m-1d(2*i+1,2*k)=(d(2*i,2*k)+d(2*(i+1),2*k))/2;d(2*i,2*k+1)=(d(2*i,2*k)+d(2*i,2*(k+1)))/2;endendfori=1:nfork=1:mifV(i,k)VmaxVmax=V(i,k);endendend%%**********************衰减系数************************%%ddx、ddz即，x方向和z方向的衰减系数R=1e-6;%理论反射系数ddx=zeros(n,m);ddz=zeros(n,m);fori=1:nfork=1:m%%区域1ifi=1&i=pml&k=1&k=pmlx=pml-k;z=pml-i;ddx(i,k)=-log(R)*3*Vmax*x^2/(2*plx^2);ddz(i,k)=-log(R)*3*Vmax*z^2/(2*plz^2);elseifi=1&i=pml&km-pml&k=mx=k-(m-pml);z=pml-i;ddx(i,k)=-log(R)*3*Vmax*x^2/(2*plx^2);ddz(i,k)=-log(R)*3*Vmax*z^2/(2*plz^2);elseifin-pml&i=n&k=1&k=pmlx=pml-k;z=i-(n-pml);ddx(i,k)=-log(R)*3*Vmax*x^2/(2*plx^2);ddz(i,k)=-log(R)*3*Vmax*z^2/(2*plz^2);elseifin-pml&i=n&km-pml&k=mx=k-(m-pml);z=i-(n-pml);ddx(i,k)=-log(R)*3*Vmax*x^2/(2*plx^2);ddz(i,k)=-log(R)*3*Vmax*z^2/(2*plz^2);%%区域2elseifi=pml&kpml&km-pml+1x=0;z=pml-i;ddx(i,k)=0;ddz(i,k)=-log(R)*3*Vmax*z^2/(2*plz^2);elseifin-pml&i=n&kpml&k=m-pmlx=0;z=i-(n-pml);ddx(i,k)=0;ddz(i,k)=-log(R)*3*Vmax*z^2/(2*plz^2);%%区域3elseifipml&i=n-pml&k=pmlx=pml-k;z=0;ddx(i,k)=-log(R)*3*Vmax*x^2/(2*plx^2);ddz(i,k)=0;elseifipml&i=n-pml&km-pml&k=mx=k-(m-pml);z=0;ddx(i,k)=-log(R)*3*Vmax*x^2/(2*plx^2);ddz(i,k)=0;endendend%figure(1),imagesc(ddz),title('z方向衰减系数');%figure(2),imagesc(ddx),title('x方向衰减系数');%%**************************************************%%**********************波场模拟********************p0=zeros(n,m);p1=zeros(n,m);px0=zeros(n,m);px1=zeros(n,m);pz0=zeros(n,m);pz1=zeros(n,m);K=zeros(n,m);Vx1=zeros(nt,mt);Vx0=zeros(nt,mt);Vz1=zeros(nt,mt);Vz0=zeros(nt,mt);fort=dt:dt:Tp0(z0,x0)=dt*V(z0,x0)^2*wave(round(t/dt));fori=2:n-1fork=2:m-1K(i,k)=d(2*i,2*k)*V(i,k)^2;Vz1(2*i+1,2*k)=((1-0.5*dt*ddz(i,k))*Vz0(2*i+1,2*k)-dt*(p0(i+1,k)-p0(i,k))/(d(2*i+1,2*k)*dz))/(1+0.5*dt*ddz(i,k));Vx1(2*i,2*k+1)=((1-0.5*dt*ddx(i,k))*Vx0(2*i,2*k+1)-dt*(p0(i,k+1)-p0(i,k))/(d(2*i,2*k+1)*dx))/(1+0.5*dt*ddx(i,k));pz1(i,k)=((1-0.5*dt*ddz(i,k))*pz0(i,k)-K(i,k)*dt*(Vz1(2*i+1,2*k)-Vz1(2*i-1,2*k))/dz)/(1+0.5*dt*ddz(i,k));px1(i,k)=((1-0.5*dt*ddx(i,k))*px0(i,k)-K(i,k)*dt*(Vx1(2*i,2*k+1)-Vx1(2*i,2*k-1))/dx)/(1+0.5*dt*ddx(i,k));p1(i,k)=px1(i,k)+pz1(i,k);endendp0=p1;pz0=pz1;px0=px1;Vz0=Vz1;Vx0=Vx1;fori=1:n-2*pmlfork=1:m-2*pmlp(i,k)=p1(i+pml,k+pml);endend%imagesc(p1),title('声波波场'),pause(0.0000001);endfigure(1),imagesc(p1);title('声波波场--交错网格，2阶');figure(2),imagesc(p);title('声波方程--交错网格、加边界');%figure(2),imagesc(pz1);title('z分量');%figure(3),imagesc(px1);title('x分量');%%*************************************************Toc波场快照