交错网格PML弹性波.docx

透射p波透射s波反射s波反射p波直达P波上图为本程序运行波场快照代码如下：//空间4阶,时间2阶,pml吸收边界条件，P-SV波交错网格数值模拟 version 3.0//程序并不完善，但应该算比较有用的/* 网格交错方式如下txx,tzz------vx-------x||||||vz---------txz|||z*/#include stdio.h#include stdlib.h#include string.h#include math.h#include time.h#define PI 3.1415926535float **space2d(int nr, intnc);void free_space2d(float **a, intnx);voidwfile(char filename[], float **data, int nr, intnc);voidcreate_model(float **vp, float **vs, float **rho, float **lamda, float **lamda2u,float **mu,int nr, intnc);float **extmodel(float **init_model,intnr,intnc,intnp);int main(){//给定参数int NX=1000;//x方向网格点数int NZ=1000;//z方向网格点数int NP=20;//pml层网格点数intsx=500+NP;//震源坐标点号intsz=500+NP;intNX_ext=NX+2*NP;intNZ_ext=NZ+2*NP; int NT=1000; //时间层数double H=1.0;//空间步长double RC=0.0001;double DP=NP*H; double DT=0.0002;//时间步长double F0=30.0;//震源主频double T0=1.2/F0; double Vpmax=2500.0; //模型最大纵波速度,用于稳定性计算double Vpmin=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;double tt,x,z,xoleft,xoright,d0;double v0,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