交错网格PML弹性波.docx
- 文档编号:18247258
- 上传时间:2023-08-14
- 格式:DOCX
- 页数:19
- 大小:55.76KB
交错网格PML弹性波.docx
《交错网格PML弹性波.docx》由会员分享,可在线阅读,更多相关《交错网格PML弹性波.docx(19页珍藏版)》请在冰点文库上搜索。
交错网格PML弹性波
上图为本程序运行波场快照
代码如下:
//空间4阶,时间2阶,pml吸收边界条件,P-SV波交错网格数值模拟version3.0
//程序并不完善,但应该算比较有用的
/*网格交错方式如下
txx,tzz------vx------->x
||
||
||
vz---------txz
|
|
|z
*/
#include
#include
#include
#include
#include
#definePI3.1415926535
float**space2d(intnr,intnc);
voidfree_space2d(float**a,intnx);
voidwfile(charfilename[],float**data,intnr,intnc);
voidcreate_model(float**vp,float**vs,float**rho,float**lamda,float**lamda2u,float**mu,intnr,intnc);
float**extmodel(float**init_model,intnr,intnc,intnp);
intmain()
{
//给定参数
intNX=1000;//x方向网格点数
intNZ=1000;//z方向网格点数
intNP=20;//pml层网格点数
intsx=500+NP;//震源坐标点号
intsz=500+NP;
intNX_ext=NX+2*NP;
intNZ_ext=NZ+2*NP;
intNT=1000;//时间层数
doubleH=1.0;//空间步长
doubleRC=0.0001;
doubleDP=NP*H;
doubleDT=0.0002;//时间步长
doubleF0=30.0;//震源主频
doubleT0=1.2/F0;
doubleVpmax=2500.0;//模型最大纵波速度,用于稳定性计算
doubleVpmin=2000.0;//模型最小纵波速度,用于控制数值频散
doubleVsmax=1300.0;
doubleVsmin=1100.0;
float**vs;//初始模型横波速度
float**vp;//初始模型纵波速度
float**rho;//初始模型密度
float**mu;//剪切模量
float**lamda;
float**lamda2u;
float**sis_x;//地震记录x分量
float**sis_z;//地震记录z分量
float**vx_x;
float**vx_z;
float**vx;
float**vz_x;
float**vz_z;
float**vz;
float**txx_x;
float**txx_z;
float**tzz_x;
float**tzz_z;
float**txz_x;
float**txz_z;
float**dxi,**dxi2;
float**dzj,**dzj2;
doublett,x,z,xoleft,xoright,d0;
doublev0,y1,y2,y3,y4,y5,y6,y7,y8,y9,y10,rho_tempx,rho_tempz,muxz,best_dt;
intix,iz,it;
vs=space2d(NZ,NX);
vp=space2d(NZ,NX);
rho=space2d(NZ,NX);
mu=space2d(NZ,NX);
lamda=space2d(NZ,NX);
lamda2u=space2d(NZ,NX);
create_model(vp,vs,rho,lamda,lamda2u,mu,NZ,NX);
charvs_name[]="vs_ext.dat";
float**vs_ext;
float**vp_ext;
float**rho_ext;
float**mu_ext;
float**lamda_ext;
float**lamda2u_ext;
vs_ext=extmodel(vs,NZ,NX,NP);
vp_ext=extmodel(vp,NZ,NX,NP);
rho_ext=extmodel(rho,NZ,NX,NP);
mu_ext=extmodel(mu,NZ,NX,NP);
lamda_ext=extmodel(lamda,NZ,NX,NP);
lamda2u_ext=extmodel(lamda2u,NZ,NX,NP);
wfile(vs_name,vs_ext,NZ_ext,NX_ext);
vx_x=space2d(NZ_ext,NX_ext);
vx_z=space2d(NZ_ext,NX_ext);
vx=space2d(NZ_ext,NX_ext);
vz_x=space2d(NZ_ext,NX_ext);
vz_z=space2d(NZ_ext,NX_ext);
vz=space2d(NZ_ext,NX_ext);
txx_x=space2d(NZ_ext,NX_ext);
txx_z=space2d(NZ_ext,NX_ext);
tzz_x=space2d(NZ_ext,NX_ext);
tzz_z=space2d(NZ_ext,NX_ext);
txz_x=space2d(NZ_ext,NX_ext);
txz_z=space2d(NZ_ext,NX_ext);
dxi=space2d(NZ_ext,NX_ext);
dxi2=space2d(NZ_ext,NX_ext);
dzj=space2d(NZ_ext,NX_ext);
dzj2=space2d(NZ_ext,NX_ext);
xoleft=DP;
xoright=(NX_ext-1)*H-DP;
for(iz=0;iz { for(ix=0;ix { x=ix*H; if(x { v0=vp_ext[iz][0]; d0=3.0*v0*log(1.0/RC)/(2.0*DP); dxi[iz][ix]=d0*pow(((xoleft-x)/DP),2); dxi2[iz][ix]=d0*pow(((xoleft-x-0.5*H)/DP),2); } elseif(x>=0.9999*xoright) { v0=vp_ext[iz][NX_ext-1]; d0=3.0*v0*log(1.0/RC)/(2.0*DP); dxi[iz][ix]=d0*pow(((x-xoright)/DP),2); dxi2[iz][ix]=d0*pow(((x+0.5*H-xoright)/DP),2); } else { dxi[iz][ix]=0.;dxi2[iz][ix]=0.; } } } xoright=(NZ_ext-1)*H-DP; for(iz=0;iz { z=iz*H; for(ix=0;ix { if(z { v0=vp_ext[0][ix]; d0=3.0*v0*log(1.0/RC)/(2.0*DP); //d0=3.0*v0*(8.0/15.0-3.0/100.0*NP+NP*NP/1500.0)/H; dzj[iz][ix]=d0*pow(((xoleft-z)/DP),2); dzj2[iz][ix]=d0*pow(((xoleft-z-0.5*H)/DP),2); } elseif(z>=0.9999*xoright) { v0=vp_ext[NZ_ext-1][ix]; d0=3.0*v0*log(1.0/RC)/(2.0*DP); dzj[iz][ix]=d0*pow(((z-xoright)/DP),2); dzj2[iz][ix]=d0*pow(((z+0.5*H-xoright)/DP),2); } else { dzj[iz][ix]=0.;dzj2[iz][ix]=0.; } } } //开始时间递推 doublea1=9.0/8.0; doublea2=-1.0/24.0; //检查稳定性条件 best_dt=6.0*H/(7.0*sqrt(2.0)*Vpmax); if(DT>=best_dt) printf("时间步长过大,应该小于%f\n",best_dt); //控制网格频散 if(Vsmin/(F0*H)<15) printf("空间步长太大,可能引起明显的网格频散\n"); clock_tstart,finish; start=clock(); sis_x=space2d(NX,NT);sis_z=space2d(NX,NT); for(it=0;it { tt=it*DT; if(it%100==0)printf("it=%d\n",it); //计算质点速度v for(iz=2;iz for(ix=2;ix { rho_tempx=0.5*(rho_ext[iz][ix]+rho_ext[iz][ix+1]); rho_tempz=0.5*(rho_ext[iz][ix]+rho_ext[iz+1][ix]); y1=(DT/(rho_tempx*H))*(a1*(txx_x[iz][ix+1]-txx_x[iz][ix]+txx_z[iz][ix+1]-txx_z[iz][ix])+ a2*(txx_x[iz][ix+2]-txx_x[iz][ix-1]+txx_z[iz][ix+2]-txx_z[iz][ix-1])); y2=(DT/(rho_tempx*H))*(a1*(txz_x[iz][ix]-txz_x[iz-1][ix]+txz_z[iz][ix]-txz_z[iz-1][ix])+ a2*(txz_x[iz+1][ix]-txz_x[iz-2][ix]+txz_z[iz+1][ix]-txz_z[iz-2][ix])); y3=(DT/(rho_tempz*H))*(a1*(txz_x[iz][ix]-txz_x[iz][ix-1]+txz_z[iz][ix]-txz_z[iz][ix-1])+ a2*(txz_x[iz][ix+1]-txz_x[iz][ix-2]+txz_z[iz][ix+1]-txz_z[iz][ix-2])); y4=(DT/(rho_tempz*H))*(a1*(tzz_x[iz+1][ix]-tzz_x[iz][ix]+tzz_z[iz+1][ix]-tzz_z[iz][ix])+ a2*(tzz_x[iz+2][ix]-tzz_x[iz-1][ix]+tzz_z[iz+2][ix]-tzz_z[iz-1][ix])); vx_x[iz][ix]=(1.0/(1.0+0.5*DT*dxi2[iz][ix]))*((1-0.5*DT*dxi2[iz][ix])*vx_x[iz][ix]+y1); vx_z[iz][ix]=(1.0/(1.0+0.5*DT*dzj[iz][ix]))*((1-0.5*DT*dzj[iz][ix])*vx_z[iz][ix]+y2); vz_x[iz][ix]=(1.0/(1.0+0.5*DT*dxi[iz][ix]))*((1-0.5*DT*dxi[iz][ix])*vz_x[iz][ix]+y3); vz_z[iz][ix]=(1.0/(1.0+0.5*DT*dzj2[iz][ix]))*((1-0.5*DT*dzj2[iz][ix])*vz_z[iz][ix]+y4); } //计算应力txy,tzy for(iz=2;iz for(ix=2;ix { muxz=1/(0.25*(1/mu_ext[iz][ix]+1/mu_ext[iz][ix+1]+1/mu_ext[iz+1][ix]+1/mu_ext[iz+1][ix+1])); //muxz=0.25*(mu_ext[iz][ix]+mu_ext[iz][ix+1]+mu_ext[iz+1][ix]+mu_ext[iz+1][ix+1]); y5=(DT*lamda2u_ext[iz][ix]/H)*(a1*(vx_x[iz][ix]-vx_x[iz][ix-1]+vx_z[iz][ix]-vx_z[iz][ix-1])+ a2*(vx_x[iz][ix+1]-vx_x[iz][ix-2]+vx_z[iz][ix+1]-vx_z[iz][ix-2])); y6=(lamda_ext[iz][ix]*DT/H)*(a1*(vz_x[iz][ix]-vz_x[iz-1][ix]+vz_z[iz][ix]-vz_z[iz-1][ix])+ a2*(vz_x[iz+1][ix]-vz_x[iz-2][ix]+vz_z[iz+1][ix]-vz_z[iz-2][ix])); y7=(lamda_ext[iz][ix]*DT/H)*(a1*(vx_x[iz][ix]-vx_x[iz][ix-1]+vx_z[iz][ix]-vx_z[iz][ix-1])+ a2*(vx_x[iz][ix+1]-vx_x[iz][ix-2]+vx_z[iz][ix+1]-vx_z[iz][ix-2])); y8=(lamda2u_ext[iz][ix]*DT/H)*(a1*(vz_x[iz][ix]-vz_x[iz-1][ix]+vz_z[iz][ix]-vz_z[iz-1][ix])+ a2*(vz_x[iz+1][ix]-vz_x[iz-2][ix]+vz_z[iz+1][ix]-vz_z[iz-2][ix])); y9=(muxz*DT/H)*(a1*(vz_x[iz][ix+1]-vz_x[iz][ix]+vz_z[iz][ix+1]-vz_z[iz][ix])+ a2*(vz_x[iz][ix+2]-vz_x[iz][ix-1]+vz_z[iz][ix+2]-vz_z[iz][ix-1])); y10=(muxz*DT/H)*(a1*(vx_x[iz+1][ix]-vx_x[iz][ix]+vx_z[iz+1][ix]-vx_z[iz][ix])+ a2*(vx_x[iz+2][ix]-vx_x[iz-1][ix]+vx_z[iz+2][ix]-vx_z[iz-1][ix])); txx_x[iz][ix]=(1.0/(1.0+0.5*DT*dxi[iz][ix]))*((1-0.5*DT*dxi[iz][ix])*txx_x[iz][ix]+y5); txx_z[iz][ix]=(1.0/(1.0+0.5*DT*dzj[iz][ix]))*((1-0.5*DT*dzj[iz][ix])*txx_z[iz][ix]+y6); tzz_x[iz][ix]=(1.0/(1.0+0.5*DT*dxi[iz][ix]))*((1-0.5*DT*dxi[iz][ix])*tzz_x[iz][ix]+y7); tzz_z[iz][ix]=(1.0/(1.0+0.5*DT*dzj[iz][ix]))*((1-0.5*DT*dzj[iz][ix])*tzz_z[iz][ix]+y8); txz_x[iz][ix]=(1.0/(1.0+0.5*DT*dxi2[iz][ix]))*((1-0.5*DT*dxi2[iz][ix])*txz_x[iz][ix]+y9); txz_z[iz][ix]=(1.0/(1.0+0.5*DT*dzj2[iz][ix]))*((1-0.5*DT*dzj2[iz][ix])*txz_z[iz][ix]+y10); } //加震源 //vz_x[sz][sx]=vz_x[sz][sx]-2*PI*PI*F0*F0*(tt-T0)*exp(-PI*PI*F0*F0*(tt-T0)*(tt-T0));//高斯函数的一阶导函数,对质点速度加集中力源 tzz_x[sz][sx]=tzz_x[sz][sx]-2*PI*PI*F0*F0*(tt-T0)*exp(-PI*PI*F0*F0*(tt-T0)*(tt-T0));//对正应力加爆炸源 txx_x[sz][sx]=txx_x[sz][sx]-2*PI*PI*F0*F0*(tt-T0)*exp(-PI*PI*F0*F0*(tt-T0)*(tt-T0)); for(iz=0;iz { vx_x[iz][0]=0.;vx_x[iz][NX_ext-1]=0.;vx_z[iz][0]=0.;vx_z[iz][NX_ext-1]=0.; vz_x[iz][0]=0.;vz_x[iz][NX_ext-1]=0.;vz_z[iz][0]=0.;vz_z[iz][NX_ext-1]=0.; } for(ix=0;ix { vx_x[NZ_ext-1][ix]=0.;vx_z[NZ_ext-1][ix]=0.; vz_x[NZ_ext-1][ix]=0.;vz_z[NZ_ext-1][ix]=0.; vx_x[0][ix]=0;vx_z[0][ix]=0; vz_x[0][ix]=0;vz_z[0][ix]=0; } //记录地震记录 for(ix=NP;ix { sis_z[ix-NP][it]=vz_x[NP][ix]+vz_z[NP][ix]; sis_x[ix-NP][it]=vx_x[NP][ix]+vx_z[NP][ix]; } } finish=clock(); printf("%fseconds\n",(double)(finish-start)/CLOCKS_PER_SEC); //输出波场快照 charvzname[]="vz.dat"; charvxname[]="vx.dat"; for(iz=0;iz for(ix=0;ix { vx[iz][ix]=vx_x[iz][ix]+vx_z[iz][ix]; vz[iz][ix]=vz_x[iz][ix]+vz_z[iz][ix]; } wfile(vxname,vx,NZ_ext,NX_ext); wfile(vzname,vz,NZ_ext,NX_ext); //输出地震记录 charsisxname[]="sisx.dat"; charsiszname[]="sisz.dat"; wfile(sisxname,sis_x,NX,NT); wfile(siszname,sis_z,NX,NT); free_space2d(tzz_x,NZ_ext); free_space2d(tzz_z,NZ_ext); free_space2d(txx_x,NZ_ext); free_space2d(txx_z,NZ_ext); free_space2d(txz_x,NZ_ext); free_space2d(txz_z,NZ_ext); free_space2d(vx_x,NZ_ext); free_space2d(vx_z,NZ_ext); free_space2d(vx,NZ_ext); free_space2d(vz_x,NZ_ext); free_space2d(vz_z,NZ_ext); free_space2d(vz,NZ_ext); free_space2d(vs_ext,NZ_ext); free_space2d(vp_ext,NZ_ext); free_space2d(rho_ext,NZ_ext); free_space2d(vs,NZ); free_space2d(vp,NZ); free_space2d(mu,NZ); free_space2d(mu_ext,NZ_ext); free_space2d(lamda,NZ); free_space2d(lamda2u,NZ); free_space2d(lamda_ext,NZ_ext); free_space2d(lamda2u_ext,NZ_ext); free_space2d(rho,NZ); free_space2d(sis_x,NX); free_space2d(sis_z,NX); free_space2d(dxi,NZ_ext); free_space2d(dxi2,NZ_ext); free_space2d(dzj,NZ_ext); free_space2d(dzj2,NZ_ext); return0; } //申请二维动态数组 float**space2d(intnr,intnc) { float**a; inti; a=(float**)calloc(nr,sizeof(float*)); for(i=0;i
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 交错 网格 PML 弹性