声波加pml边界条件.docxVIP

雷克子波 因计算机资源有限，不便在太大的区域做PML的计算，一般只在边界上取有限宽度的区域作为PML计算区域。根据相关文献的研究，PML边界区域最少长度应为半个波长[6]。本文综合考虑了效果与开销等因素，选取了边界上50层作为PML的计算区域。常规计算区域与PML边界区域的如图2-4所示。 衰减系数 式中为PML层的厚度，为层内的点距PML与非PML的边界的距离，为纵波速度，那么在PML边界区域内，对于式（2-13）即为理论反射系数，一般取0.001较为合适，为方向的空间步长。 ，可看作为在常规的计算方程基础上，减去一项进行PML的阻尼修正项。因本文中只考虑各项同性介质中的地震波传播规律，故可做假设。在此利用一下三个假设： 因为以上三个近似精度均为时间方向上的近似，且时间精度均为二阶精度，因交错网格技术的时间精度为二阶，故以上近似不影响本式的计算精度。 故可得： （2-18a） 同理： （2-18b） （2-18c） 此为在PML边界区域内的弹性声波应力-速度方程组。 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%交错网格---非均匀介质二维声波方程(一阶压力--速度)、2阶时间差分、2阶空间差分精度 %%加上边界 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% close all;clear,clc tic %%***********************震源为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; % 纵向网格大小，单位m dx=5; % 横向网格大小，单位m dt=0.0001; % 时间步长，单位s T=0.5; % 波动传播时间，单位s wave(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; % 纵波最大速度 %%Setting Velocity Density zt=-750-plz:dz/2:750+plz; xt=-1000-plx:dx/2:1000+plx; % 速度与密度采样区间 nt=length(zt); mt=length(xt); % 速度与密度采样点数目 V=zeros(n,m); % 介质速度,m/s d=zeros(nt,mt); % 介质密度,kg/m^3 %%均匀介质模型 for i=1:n for k=1:m V(i,k)=2.0e3; end end for i=1:n for k=1:m d(2*i,2*k)=2.3e3; end end % % %%层状介质模型 % % for i=1:n % % for k=1:m % % if i round(n/3) % % V(i,k)=2.3e3; % % else % % V(i,k)=3.0e3; % % end % % end % % end for i=1:n-1 for k=1:m-1 d(2*i+1,2*k)=(d(2*i,2*k)+d(2*(i+1),2*k))/2;