[2018年最新整理]声波加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,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; % 纵波最大速度%%Setting Velocity 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%%均匀介质模型for i=1:nfor k=1:m V(i,k)=2.0e3;endendfor i=1:nfor k=1:m d(2*i,2*k)=2.3e3;endend% % %%层状介质模型% % 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% % endfor i=1:n-1for k=1:m-1 d(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;endendfor i=1:nfor k=1:mif V(i,k) Vmax Vmax=V(i,k);endendend%%**********************衰减系数************