声波加pml边界条件(共11页).docx

上传人:飞****2 文档编号:13838025 上传时间:2022-05-01 格式:DOCX 页数:11 大小:119.07KB
返回 下载 相关 举报
声波加pml边界条件(共11页).docx_第1页
第1页 / 共11页
声波加pml边界条件(共11页).docx_第2页
第2页 / 共11页
点击查看更多>>
资源描述

《声波加pml边界条件(共11页).docx》由会员分享,可在线阅读,更多相关《声波加pml边界条件(共11页).docx(11页珍藏版)》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。

1、精选优质文档-倾情为你奉上 雷克子波因计算机资源有限,不便在太大的区域做PML的计算,一般只在边界上取有限宽度的区域作为PML计算区域。根据相关文献的研究,PML边界区域最少长度应为半个波长6。本文综合考虑了效果与开销等因素,选取了边界上50层作为PML的计算区域。常规计算区域与PML边界区域的如图2-4所示。衰减系数式中为PML层的厚度,为层内的点距PML与非PML的边界的距离,为纵波速度,那么在PML边界区域内,对于式(2-13)即为理论反射系数,一般取0.001较为合适,为方向的空间步长。,可看作为在常规的计算方程基础上,减去一项进行PML的阻尼修正项。因本文中只考虑各项同性介质中的地震

2、波传播规律,故可做假设。在此利用一下三个假设:因为以上三个近似精度均为时间方向上的近似,且时间精度均为二阶精度,因交错网格技术的时间精度为二阶,故以上近似不影响本式的计算精度。故可得:(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(-

3、(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:

4、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/m3 %均匀介质模型fo

5、r i=1:n for k=1:m V(i,k)=2.0e3; endendfor i=1:n for k=1:m d(2*i,2*k)=2.3e3; endend % % %层状介质模型% % for i=1:n% % for k=1:m% % if i Vmax Vmax=V(i,k); end endend%*衰减系数*% ddx、ddz 即,x方向和z方向的衰减系数R=1e-6; % 理论反射系数ddx=zeros(n,m); ddz=zeros(n,m); for i=1:n for k=1:m % 区域1 if i=1 & i=1 & k=1 & im-pml & kn-pml &

6、 i=1 & kn-pml & im-pml & k=m x=k-(m-pml);z=i-(n-pml); ddx(i,k)=-log(R)*3*Vmax*x2/(2*plx2); ddz(i,k)=-log(R)*3*Vmax*z2/(2*plz2); % 区域2 elseif ipml & kn-pml & ipml & kpml & i=n-pml & kpml & im-pml & k=m x=k-(m-pml);z=0; ddx(i,k)=-log(R)*3*Vmax*x2/(2*plx2);ddz(i,k)=0; end endend% figure(1),imagesc(ddz)

7、,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); for t=dt:dt:T p0(z0,x0)=dt*V(z0,x0)2*wave(round(t/dt); for

8、i=2:n-1 for k=2:m-1 K(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

9、)-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); end end p0=p1; pz0=pz1;px0=px1; Vz0=Vz1;Vx0=Vx1; for i=1:n-2*pml for k=1:m-2*pml p(i,k)=p1(i+pml,k+pml); end end% imagesc(p1),title(声波波场),pause(0.);endfigure(1),imagesc(p1);title(声波波场-交错网格,2阶);figure(2),imagesc(p);title(声波方程-交错网格、加边界);% figure(2),imagesc(pz1);title(z分量);% figure(3),imagesc(px1);title(x分量);%*Toc波场快照专心-专注-专业

展开阅读全文
相关资源
相关搜索

当前位置:首页 > 教育专区 > 教案示例

本站为文档C TO C交易模式,本站只提供存储空间、用户上传的文档直接被用户下载,本站只是中间服务平台,本站所有文档下载所得的收益归上传人(含作者)所有。本站仅对用户上传内容的表现方式做保护处理,对上载内容本身不做任何修改或编辑。若文档所含内容侵犯了您的版权或隐私,请立即通知淘文阁网,我们立即给予删除!客服QQ:136780468 微信:18945177775 电话:18904686070

工信部备案号:黑ICP备15003705号© 2020-2023 www.taowenge.com 淘文阁