声波方程数值模拟实验报告(共12页)

1、精选优质文档-倾情为你奉上声波方程数值模拟实验报告一基础理论知识需要的已知条件包括：11）震源函数2）地层速度（波速）3）边界条件2弹性波方程：声波方程的有限差分法数值模拟对于二维速度-深度模型，地下介质中地震波的传播规律可以近似地用声波方程描述： （4-1）是介质在点（x , z）处的纵波速度，为描述速度位或者压力的波场，为震源函数。为求式（4-1）的数值解，必须将此式离散化，即用有限差分来逼近导数，用差商代替微商。为此，先把空间模型网格化（如图4-1所示）。设x、z方向的网格间隔长度为，为时间采样步长，则有： （i为正整数） (j为正整数) (n为正整数)表示在(i,j)点，k时刻的波场值

2、。将在(i,j)点k时刻用Taylor展式展开：(4-2)将在(i,j)点k时刻用Taylor展式展开：(4-3)将上两式相加，略去高阶小量，整理得(i,j)点k时刻的二阶时间微商为：(4-4)对于空间微分，采用四阶精度差分格式，（以X方向为例）即将、分别在(i,j)点k时刻展开到四阶小量，消除四阶小量并解出二阶微分得： (4-5)同理可得： (4-6)这就实现了用网个点波场值的差商代替了偏微分方程的微商，将上三个式子代入(4-1)式中得： (4-7)式中为介质速度的空间离散值，是空间离散步长，为时间离散步长，为震源函数，关于一般使用一个理论的雷克型子波代替，即： （4-8）上式中，为时间，

3、为中心频率，一般取为20-40HZ，为控制频带宽度的参数，一般取3-5。在实际计算过程中，需把此震源函数离散，参与波场计算。确定震源位置。稳定性条件：(4-8)这里表示的是地下介质的最大波速；若地下介质网格间隔、最小速度、及时间采样间隔不符合(4-8)式时，第推求解（4-7）式，波场值会出现误差（高阶小量）累积，出现不稳定现象。频散关系式： (4-9)式中为最小速度，为Nyquist频率。一般取震源子波中的主频的2倍值参与计算，G为每个波长所占的网格点数，对于空间二阶差分、时间二阶差分取8，而对于空间为四阶差分的情况则取4方能有效减少频散。二实验步骤1、 应用声波方程作为正演模拟的波动方程，忽

4、略转换波的产生、传播；2、 将所提供震源函数离散后绘图；震源函数为雷克子波，离散绘图如下：fm=30;r=3;t=0.002; for n=1:200w(n)=exp(-(2*pi*fm/r)2*(t*n)2)*cos(2*pi*fm*t*n); endfigure(1),plot(w);3、 对于小模型，整个区域的速度值可设为常数，即只有一种介质，将震源点放在模型中间，分别记录两个时刻的波前快照（即区域内所有网格点的波场值）。第一时刻为地震波还未传播到边界上的某时刻。由稳定条件，设v为2000，可以令DT=0.001,DH=5,此时精度较高，且满足频散关系程序，图形如下：#include &

5、lt;stdio.h>#include <math.h>#include <stdlib.h>#define PI 3.#define FM 30#define R 3#define KN 200#define XN 101#define ZN 101#define DH 5#define DT 0.001void main()FILE *fp;int i,j,k,m,n; float u1XNZN,u2XNZN,u3XNZN,u4XNZN,fXNZN; /不能直接初值为0 float u5XNZN,vXNZN,wKN,uu0,uu1,uu2; for(k=0;k

6、<KN;k+) wk=exp(-(2*PI*FM/R)*(2*PI*FM/R)*(k*DT)*(k*DT)*cos(2*PI*FM*k*DT); for(i=0;i<XN;i+) /定义f函数，当且仅当i,j同时为50时，f为1，其余为0 for(j=0;j<ZN;j+) if(i=50&&j=50) fij=1; else fij=0; for(i=0;i<XN;i+) for(j=0;j<ZN;j+) u1ij=0.0; u2ij=0.0; u3ij=0.0; u4ij=0.0; vij=2000; /速度相同表示同一介质 for(k=0;k&

7、lt;KN;k+) for(i=2;i<XN-2;i+) for(j=2;j<ZN-2;j+) uu0=(vij)*(vij)*(DT/DH)*(DT/DH); uu1=-1.0/12*(u2i-2j+u2i+2j)+4.0/3*(u2i-1j+u2i+1j)-5.0/2*u2ij; uu2=-1.0/12*(u2ij-2+u2ij+2)+4.0/3*(u2ij-1+u2ij+1)-5.0/2*u2ij; u3ij=2*u2ij-u1ij+uu0*uu1+uu0*uu2+wk*fij; for(m=0;m<XN;m+) for(n=0;n<ZN;n+) u1mn=u2m

8、n; u2mn=u3mn; if(k=100) for(m=0;m<XN;m+) for(n=0;n<ZN;n+) u4mn=u3mn;/记录波前快照,中间点 if(fp=fopen("wavefront.dat","w")!=NULL) fprintf(fp,"%dn",XN); fprintf(fp,"%dn",ZN); for(i=0;i<XN;i+) for(j=0;j<ZN;j+) fprintf(fp,"%fn",u4ij); fclose(fp); 第二时刻

9、为地震波已经传播到边界上的某时刻，体会其人工边界反射；程序图形如下：#include <stdio.h>#include <math.h>#include <stdlib.h>#define PI 3.#define FM 30#define R 3#define KN 200#define XN 101#define ZN 101#define DH 5#define DT 0.001void main()FILE *fp;int i,j,k,m,n; float u1XNZN,u2XNZN,u3XNZN,u4XNZN,fXNZN; /不能直接初值为0 f

10、loat u5XNZN,vXNZN,wKN,uu0,uu1,uu2; for(k=0;k<KN;k+) wk=exp(-(2*PI*FM/R)*(2*PI*FM/R)*(k*DT)*(k*DT)*cos(2*PI*FM*k*DT); for(i=0;i<XN;i+) /定义f函数，当且仅当i,j同时为50时，f为1，其余为0 for(j=0;j<ZN;j+) if(i=50&&j=50) fij=1; else fij=0; for(i=0;i<XN;i+) for(j=0;j<ZN;j+) u1ij=0.0; u2ij=0.0; u3ij=0.0

11、; u4ij=0.0; vij=2000; /速度相同表示同一介质 for(k=0;k<KN;k+) for(i=2;i<XN-2;i+) for(j=2;j<ZN-2;j+) uu0=(vij)*(vij)*(DT/DH)*(DT/DH); uu1=-1.0/12*(u2i-2j+u2i+2j)+4.0/3*(u2i-1j+u2i+1j)-5.0/2*u2ij; uu2=-1.0/12*(u2ij-2+u2ij+2)+4.0/3*(u2ij-1+u2ij+1)-5.0/2*u2ij; u3ij=2*u2ij-u1ij+uu0*uu1+uu0*uu2+wk*fij; for(

12、m=0;m<XN;m+) for(n=0;n<ZN;n+) u1mn=u2mn; u2mn=u3mn; if(k=160) /只需改动K值，即显示边界反射 for(m=0;m<XN;m+) for(n=0;n<ZN;n+) u4mn=u3mn;/记录波前快照，边界 if(fp=fopen("wavefront.dat","w")!=NULL) fprintf(fp,"%dn",XN); fprintf(fp,"%dn",ZN); for(i=0;i<XN;i+) for(j=0;j&l

13、t;ZN;j+) fprintf(fp,"%fn",u4ij); fclose(fp); 4、 对于大模型，定义为水平层状速度模型；做两个实验，一是将震源点放在区域表层任一点，记录下某些时刻的波前快照，体会地震波在两种介质的分界面上传播规律，指出哪是反射波，哪是透射波；这时取小模型的常量，为减少频散，速度v至少为2400程序图形如下：#include <stdio.h>#include <math.h>#include <stdlib.h>#define PI 3.#define FM 30#define R 3#define KN 20

14、0#define XN 200#define ZN 100#define DH 5#define DT 0.001void main()FILE *fp;int i,j,k,m,n; float u1XNZN,u2XNZN,u3XNZN,u4XNZN,u5XNZN; /不能直接初值为0 float fXNZN,vXNZN,wKN,uu0,uu1,uu2; for(k=0;k<KN;k+) wk=exp(-(2*PI*FM/R)*(2*PI*FM/R)*(k*DT)*(k*DT)*cos(2*PI*FM*k*DT); for(i=0;i<XN;i+) for(j=0;j<ZN;

15、j+) u1ij=0.0; u2ij=0.0; u3ij=0.0; u4ij=0.0; fij=0.0; if (j<=30) vij=2400; /第一层速度为2400 else vij=3000; /第二层速度为3000 f10010=1; /定义f函数，在100，10为1 for(k=0;k<KN;k+) for(i=2;i<XN-2;i+) for(j=2;j<ZN-2;j+) uu0=(vij)*(vij)*(DT/DH)*(DT/DH); uu1=-1.0/12*(u2i-2j+u2i+2j)+4.0/3*(u2i-1j+u2i+1j)-5.0/2*u2ij

16、; uu2=-1.0/12*(u2ij-2+u2ij+2)+4.0/3*(u2ij-1+u2ij+1)-5.0/2*u2ij; u3ij=2*u2ij-u1ij+uu0*uu1+uu0*uu2+wk*fij; u5ik=u1i10; /地震记录 for(m=0;m<XN;m+) for(n=0;n<ZN;n+) u1mn=u2mn; u2mn=u3mn; if(k=100) for(m=0;m<XN;m+) for(n=0;n<ZN;n+) u4mn=u3mn;/记录波前快照 if(fp=fopen("wavefront.dat","w&q

17、uot;)!=NULL) fprintf(fp,"%dn",XN); fprintf(fp,"%dn",ZN); for(i=0;i<XN;i+) for(j=0;j<ZN;j+) fprintf(fp,"%fn",u4ij); fclose(fp); 反射波透射波二是合成一个地震记录，即记录下与震源同一深度点的各点所有时刻的波场值，并指出记录上的同向轴分别对应哪些波？这时取小模型的常量，为减少频散，速度v至少为2400程序图像如下：#include <stdio.h>#include <math.h&g

18、t;#include <stdlib.h>#define PI 3.#define FM 30#define R 3#define KN 200#define XN 200 #define ZN 200#define DH 5#define DT 0.001void main()FILE *fp;int i,j,k,m,n; float u1XNZN,u2XNZN,u3XNZN,u4XNZN,u5XNZN; /不能直接初值为0 float fXNZN,vXNZN,wKN,uu0,uu1,uu2; for(k=0;k<KN;k+) wk=exp(-(2*PI*FM/R)*(2*

19、PI*FM/R)*(k*DT)*(k*DT)*cos(2*PI*FM*k*DT); for(i=0;i<XN;i+) for(j=0;j<ZN;j+) u1ij=0.0; u2ij=0.0; u3ij=0.0; u4ij=0.0; fij=0.0; if (j<=30) vij=2400; /第一层速度为2400 else vij=3100; /第二层速度为3100 f10010=1; /定义f函数，在100，10为1 for(k=0;k<KN;k+) for(i=2;i<XN-2;i+) for(j=2;j<ZN-2;j+) uu0=(vij)*(vij)*(DT/DH)*(DT/DH); uu1=-1.0/12*(u2i-2j+u2i+2j)+4.0/3*(u2i-1j+u2i+1j)-5.0/2*u2ij; uu2=-1.0/12*(u2ij-2+u2ij+2)+4.0/3*(u2ij-1+u2ij+1)-5.0/2*u2ij; u3ij=2*u2ij-u1ij+uu0*uu1+uu0*uu2+wk*fij; u5i