声波方程数值模拟实验报告文档格式.docx
《声波方程数值模拟实验报告文档格式.docx》由会员分享,可在线阅读,更多相关《声波方程数值模拟实验报告文档格式.docx(34页珍藏版)》请在冰点文库上搜索。
2)地层速度(波速)
3)边界条件
2.弹性波方程:
声波方程的有限差分法数值模拟
对于二维速度-深度模型,地下介质中地震波的传播规律可以近似地用声波方程描述:
(4-1)
是介质在点(x,z)处的纵波速度,
为描述速度位或者压力的波场,
为震源函数。
为求式(4-1)的数值解,必须将此式离散化,即用有限差分来逼近导数,用差商代替微商。
为此,先把空间模型网格化(如图4-1所示)。
设x、z方向的网格间隔长度为
,
为时间采样步长,则有:
(i为正整数)
(j为正整数)
(n为正整数)
表示在(i,j)点,k时刻的波场值。
将
在(i,j)点k时刻用Taylor展式展开:
(4-2)
(4-3)
将上两式相加,略去高阶小量,整理得(i,j)点k时刻的二阶时间微商为:
(4-4)
对于空间微分,采用四阶精度差分格式,(以X方向为例)即将
、
分别在(i,j)点k时刻展开到四阶小量,消除四阶小量并解出二阶微分得:
(4-5)
同理可得:
(4-6)
这就实现了用网格点波场值的差商代替了偏微分方程的微商,将上三个式子代入(4-1)式中得:
}
(4-7)
式中
为介质速度的空间离散值,
是空间离散步长,
为时间离散步长,
为震源函数,关于
一般使用一个理论的雷克型子波代替,即:
(1)
(4-8)
(2)
(4-9)
上式中,
为时间,
为中心频率,一般取为20-40HZ,
为控制频带宽度的参数,一般取3-5。
在实际计算过程中,需把此震源函数离散,参与波场计算。
确定震源位置。
稳定性条件:
(4-10)
这里
表示的是地下介质的最大波速;
若地下介质网格间隔、最大速度及时间采样间隔不符合(4-8)式时,递推求解(4-7)式,波场值会出现误差(高阶小量)累积,出现不稳定现象。
频散关系式:
(4-11)
为最小速度,
为Nyquist(奈奎斯特)频率。
一般取震源子波中的主频
的2倍值参与计算,G为每个波长所占的网格点数,对于空间二阶差分、时间二阶差分
取8,而对于空间为四阶差分的情况则
取4方能有效减少频散。
同理:
对于空间微分,采用二阶精度差分格式为
(4-12)
考虑Reynolds边界条件(详细推到见文献“基于Marmousi模型的声波方程有限差分正演算法”),差分格式如下:
(4-13)
注:
其中参数的设置:
借用以前实验的数据,当然可以根据限制条件4-10、4-11计算得到;
至于震源放于[101*5,101*5]的矩形中心,根据速度与位移可以计算传播到边界时的时间,此处忽略。
二.实验步骤
1、应用声波方程作为正演模拟的波动方程,忽略转换波的产生、传播;
2、将所提供震源函数离散后绘图;
震源函数为雷克子波,离散绘图如下:
fm=30;
r=3;
t=0.002;
forn=1:
50
w(n)=exp(-(2*pi*fm/r)^2*(t*n)^2)*cos(2*pi*fm*t*n);
end
subplot(211)
plot(w)
t1=0.002;
w1(n)=(1-2*(pi*fm*(t1*n-1/fm))^2)*exp(-(pi*fm*(t1*n-1/fm))^2);
subplot(212)
plot(w1)
3、对于小模型,整个区域的速度值可设为常数,即只有一种介质,将震源点放在模型中间,分别记录两个时刻的波前快照(即区域内所有网格点的波场值)。
第一时刻为地震波还未传播到边界上的某时刻。
由稳定条件,设v为2000,可以令DT=0.001,DH=5,此时精度较高,且满足频散关系
程序,图形如下:
#include<
stdio.h>
math.h>
stdlib.h>
#definePI3.141593
#defineFM30
#defineR3
#defineKN200
#defineXN101
#defineZN101
#defineDH5
#defineDT0.001
voidfunction(constintflag1,constintflag2)
{
FILE*fp;
inti,j,k,m,n;
floatu1[XN][ZN],u2[XN][ZN],u3[XN][ZN],u4[XN][ZN],f[XN][ZN];
//不能直接初值为0
floatv[XN][ZN],w[KN],uu0,uu1,uu2;
for(i=0;
i<
XN;
i++)//定义f函数,当且仅当i,j同时为50时,f为1,其余为0
for(j=0;
j<
ZN;
j++)
{if(i==50&
&
j==50)
f[i][j]=1;
else
f[i][j]=0;
}
i++)
for(j=0;
{
u1[i][j]=0.0;
u2[i][j]=0.0;
u3[i][j]=0.0;
u4[i][j]=0.0;
v[i][j]=2000;
//速度相同表示同一介质
}
if(0==flag1)
{
for(k=0;
k<
KN;
k++)
w[k]=exp(-(2*PI*FM/R)*(2*PI*FM/R)*(k*DT)*(k*DT))*cos(2*PI*FM*k*DT);
elseif(1==flag1)
w[k]=(1-2*pow((PI*FM*(k*DT-1./FM)),2))*exp(-2*pow((PI*FM*(k*DT-1./FM)),2));
for(k=0;
k++)
if(0==flag2)//时间二阶差分格式、空间四阶差分格式
{
for(i=2;
XN-2;
for(j=2;
ZN-2;
{
uu0=(v[i][j])*(v[i][j])*(DT/DH)*(DT/DH);
uu1=-1.0/12*(u2[i-2][j]+u2[i+2][j])+4.0/3*(u2[i-1][j]+u2[i+1][j])-5.0/2*u2[i][j];
uu2=-1.0/12*(u2[i][j-2]+u2[i][j+2])+4.0/3*(u2[i][j-1]+u2[i][j+1])-5.0/2*u2[i][j];
u3[i][j]=2*u2[i][j]-u1[i][j]+uu0*uu1+uu0*uu2+w[k]*f[i][j];
}
}
elseif(1==flag2)//时间二阶差分格式、空间二阶差分格式
uu1=u2[i+1][j]-2*u2[i][j]+u2[i-1][j];
uu2=u2[i][j+1]-2*u2[i][j]+u2[i][j-1];
for(m=0;
m<
m++)
for(n=0;
n<
n++)
u1[m][n]=u2[m][n];
u2[m][n]=u3[m][n];
if(k==100)
for(m=0;
u4[m][n]=u3[m][n];
//记录波前快照,中间点
if(0==flag1&
0==flag2)
fp=fopen("
wavefront14.dat"
"
w"
);
elseif(1==flag1&
wavefront24.dat"
elseif(0==flag1&
1==flag2)
wavefront12.dat"
wavefront22.dat"
if(fp!
=NULL)
for(i=0;
for(j=0;
fprintf(fp,"
%d%d%f\n"
i,j,u4[i][j]);
fclose(fp);
voidmain()
function(0,0);
//时间二阶差分格式、空间四阶差分格式、震源1
function(0,1);
//时间二阶差分格式、空间二阶差分格式、震源1
function(1,0);
//时间二阶差分格式、空间四阶差分格式、震源2
function(1,1);
//时间二阶差分格式、空间二阶差分格式、震源2
(1)时间二阶差分格式、空间四阶差分格式、震源1
(2)时间二阶差分格式、空间二阶差分格式、震源1
(3)时间二阶差分格式、空间四阶差分格式、震源2
(4)时间二阶差分格式、空间二阶差分格式、震源2
第二时刻为地震波已经传播到边界上的某时刻,体会其人工边界反射(此处用Reynolds边界条件);
程序图形如下:
floatv[XN][ZN],w[KN],uu0,uu1,uu2,uu3;
m++)//波向前传播记录连续三个时刻的值
for(n=0;
{
u1[m][n]=u2[m][n];
//Reynolds边界条件
uu3=v[i][j]*(DT/DH);
u3[0][j]=u2[0][j]+u2[1][j]-u1[1][j]+uu3*(u2[1][j]-u2[0][j]-u1[2][j]+u1[1][j]);
//左边界
u3[XN][j]=u2[XN][j]+u2[XN-1][j]-u1[XN-1][j]+uu3*(u2[XN-1][j]-u2[XN][j]-u1[XN-2][j]+u1[XN-1][j]);
//右边界
u3[i][0]=u2[i][0]+u2[i][1]-u1[i][1]+uu3*(u2[i][1]-u2[i][0]-u1[i][2]+u1[i][1]);
//下边界
u3[i][ZN]=u2[i][ZN]+u2[i][ZN-1]-u1[i][ZN-1]+uu3*(u2[i][ZN-1]-u2[i][ZN]-u1[i][ZN-2]+u1[i][ZN-1]);
//上边界
if(k==180)//只需改动K值,即显示边界反射
for(m=0;
for(n=0;
u4[m][n]=u3[m][n];
//记录波前快照,边界
4、对于大模型,定义为水平层状速度模型;
做两个实验,一是将震源点放在区域表层任一点,记录下某些时刻的波前快照,体会地震波在两种介质的分界面上传播规律,指出哪是反射波,哪是透射波;
这时取小模型的常量,为减少频散,速度v至少为2400
#defineXN200
#defineZN100
floatu1[XN][ZN],u2[XN][ZN],u3[XN][ZN],u4[XN][ZN];
floatf[XN][ZN],v[XN][ZN],w[KN],uu0,uu1,uu2;
f[i][j]=0.0;
if(j<
=30)
v[i][j]=2400;
//第一层速度为2400
else
v[i][j]=3000;
//第二层速度为3000
f[100][10]=1;
//定义f函数,在100,10为1
//记录波前快照
wa