OUC地震波场论数值模拟实验报告.docVIP

声波方程数值模拟试验报 姓名: XXX 专业年级:XXXXXXXXXXXX 学号:XXXXXXXXXXX 第一部分：实验结果展示 图1：零相位的雷克子波 图2：单层介质中的波前快照 图3：经过人工界面反射后的波前快照 图4：两层速度模型震源为任意一点某时刻的波前快照 图5：与震源同一深度点的各点所有时刻的波场值 第二部分：实验原理 对于二维速度-深度模型，地下介质中地震波的传播规律可以近似地用声波方程描述： （4-1） 是介质在点（x , z）处的纵波速度， 为描述速度位或者压力的波场， 为震源函数。 图4-1 差分网格划分示意图 网格间隔长度 Δh 时间采样步长 Δt X=i Δh Z=j Δh T=k Δt 表示(i,j)点k时刻的波场值 时间二阶、空间二阶差分格式推导如下： 将上两式相加，略去高阶小量，整理得(i,j)点k时刻的二阶时间微商为： 同理可得(i,j)点k时刻的二阶空间微商分别为： 这就实现了用网个点波场值的差商代替了偏微分方程的微商，将上三个式子代入(4-1)式中得： 式中 同理可得空间四阶精度的差分格式为 FM=20; 同理可得时间二阶、空间四阶精度的声波方程差分格式为： 一般使用一个理论的雷克型子波代替，即： 第三部分：程序源代码及解释 1、雷克子波源程序如下： Xn =100; Zn =100; R =3.0; dt =0.001 Dh= 5.0; V=2500.0; n=1:400; f=exp(-(2*pi*FM*n*dt/R).*(2*pi*FM*n*dt/R)).*cos(2*pi*FM*n*dt); Plot(n,f); 2、小模型代码 对于小模型，输入时间是只要输入不同的值，就可以观察不同时刻的波前值，当n较小时为没到达边界时的值，n较大时为人工反射的波前值，第一部分中的两幅图为其中的两个代表。 C程序代码： #include stdio.h #include math.h #define PI 3.1415926 #define FM 20 #define Xn 100 #define Zn 100 #define R 3.0 #define dt 0.001 #define dh 5.0 #define v 2500.0 #define i_0 50 #define j_0 50 int delta(int x) //定义δ函数 { if(x==i_0) return 1; else return 0; } double s(int n) //定义雷克子波 { double f; f=exp(-pow(2*PI*FM*n*dt/R,2))*cos(2*PI*FM*n*dt); return f; } void main() { int i, j, k, l,m,n,x,y; double p0[Xn][Zn],p1[Xn][Zn],p2[Xn][Zn]; FILE *fp; for(i=0;iXn;i++) //对数组p0,p1,p2初始化为零 for(j=0;jZn;j++) { p0[i][j]=0; p1[i][j]=0; p2[i][j]=0; } printf(请输入时间n(n为dt的整数倍)n=\n); //由用户输入想要得到那一时刻的波前值 scanf(%d,n); for(k=0;kn;k++) { for(i=2;iXn-2;i++) for(j=2;jZn-2;j++) p2[i][j]=2*p1[i][j]-p0[i][j]+pow(v*dt/dh,2)*(-1/12*(p1[i-2][j]+p1[i+2][j])+4/3*(p1[i-1][j]+p1[i+1][j])-5/2*p1[i][j])+pow(v*dt/dh,2)*(-1/12*(p1[i][j-2]+p1[i][j+2])+4/3*(p1[i][j-1]+p1[i][j+1])-5/2*p1[i][j])+s(k)*delta(i)*delta(j); for(x=0;xZn;x++) for(y=0;yZn;y++)// printf(%f\n,p2[i][j]); { p0[x][y]=p1[x][y]; p1[x][y]=p2[x][y]; } } if((fp=fopen(wavefront.dat,w))!=NULL) { fprintf