聲波方程數(shù)值模擬實驗報告

1、1.2)3)>2u2fcu=Vp 2t:x22w 2 w 廠二 V （T :t；x聲波方程的有限差分法數(shù)值模擬對于二維速度-深度模型，地下介質(zhì)中地震波的傳播規(guī)律可以近似地用聲波方程描述:2.彈性波方程:二b s(t)一 z22_w)；z2)聲波方程數(shù)值模擬實驗報告.基礎(chǔ)理論知識需要的已知條件包括：1）震源函數(shù)地層速度（波速） 邊界條件.2.2.2(4-1)葦=v2(諾二2)S(t)一 t:x:zv（x, Z）是介質(zhì)在點（x , z ）處的縱波速度，u為描述速度位或者壓力的波場，s（t）為震源函數(shù)。為求式（4-1）的數(shù)值解，必須將此式離散化，即用有限差分來逼近導(dǎo)數(shù)，用差商代替微商。為此，先

2、把空間模型網(wǎng)格化（如圖4-1所示）。；2，i 1,L=2J/i +1, j 2i +2, j /id, j 一i 1, j.J1i +1, ji +2, j d一，j/1,ji丿，j寸i +1, j十2, jj1Fi 4, j勺i 1,j也J, j +1i +1, j +1i +2, j勺i 4, j七i 1,j七.J, j七i +1,i +2i +2, j七設(shè)x、z方向的網(wǎng)格間隔長度為h，厶t為時間采樣步長，則有:X二ih（i為正整數(shù)）z = j.）h（j為正整數(shù)） n ：t（n為正整數(shù)）u：j表示在（i,j）點，k時刻的波場值。k將ui,j在（i,j）點k時刻用Taylor展式展開：(4-

3、2)k 1 kjUui,j 7,j ct將Uikj在(i,j)點k時刻用Taylor展式展開:k 1 k；'UUi,jUi,2*72 ；:t2(4-3)將上兩式相加，略去高階小量,整理得(i,j)點k時刻的二階時間微商為:2k 1kk J；：u ui,j -2ui,j ui,j.:t242(4-4)對于空間微分，采用四階精度差分格式,(以X方向為例)即將Ui*；j、*Tj分別在(i,j)點k時刻展開到四階小量，消除四階小量并解出二階微分得:-2 u 11 k k 4 k k 5 kT7ui 二 j +uid2,j+：ui4j +Ui*,j =ui,j1232:x lx2(4-5)同理可

4、得:2；=u 1.1kk4 kk 5 k一 2 人 2一 石Ui,j，+Ui,j2+：Ui,j+Ui,j；Ui,j :z-z1232(4-6)這就實現(xiàn)了用網(wǎng)個點波場值的差商代替了偏微分方程的微商，將上三個式子代入(4-1)式中得:k 1Ui,jkk 4= 2ui,j Ui,jV :-11 k k 4 k kpu=j Ui/ yijj：h25 kH-Ui,jh2 12k k 4 k k 5 kUi2j ui 2,j3ui 4,j ui1,j ?Ui,jS(t)*、(i -i°)*(j - j°)(4-7)式中v(i, j)為介質(zhì)速度的空間離散值，：h是空間離散步長，=t為時間

5、離散步長，s(k)為震源函數(shù)，關(guān)于 s(k) 一般使用一個理論的雷克型子波代替，即:s(t)二 e(-2 f / )2t2cos2 二 ft(4-8)上式中，t為時間，f為中心頻率，一般取為20-40HZ，為控制頻帶寬度的參數(shù),般取3-5。在實際計算過程中，需把此震源函數(shù)離散， 參與波場計算。（i 一 i0）"* （i - jo）確定震源位置。(4-8)穩(wěn)定性條件：Vmax * At / Ah £8這里Vmax表示的是地下介質(zhì)的最大波速；若地下介質(zhì)網(wǎng)格間隔、最小速度、及時間采 樣間隔不符合（4-8）式時，第推求解（4-7）式，波場值會出現(xiàn)誤差（高階小量）累積，出現(xiàn) 不穩(wěn)定現(xiàn)

6、象。頻散關(guān)系式：-Vmin /（GfN）（4-9）式中Vmin為最小速度，fN為Nyquist頻率。一般取震源子波中的 主頻f的2倍值參與 計算，G為每個波長所占的網(wǎng)格點數(shù)，對于空間二階差分、時間二階差分G取8,而對于空間為四階差分的情況則 G取4方能有效減少頻散。忽略轉(zhuǎn)換波的產(chǎn)生、傳播二.實驗步驟1、應(yīng)用聲波方程作為正演模擬的波動方程,2、將所提供震源函數(shù)離散后繪圖； 震源函數(shù)為雷克子波，離散繪圖如下：fm=30;r=3;t=0.002;for n=1:200w( n)=exp(-(2*pi*fm/r)A2*(t* n) A2)*cos(2*pi*fm*t* n); endfigure(1)

7、,plot(w);3、對于小模型，整個區(qū)域的速度值可設(shè)為常數(shù)，即只有一種介質(zhì)，將震源點放在模型中間,分別記錄兩個時刻的波前快照(即區(qū)域內(nèi)所有網(wǎng)格點的波場值)。第一時刻為地震波還未傳播到邊界上的某時刻。由穩(wěn)定條件，設(shè) v為2000，可以令DT=0.001,DH=5,此時精度較高，且滿足頻散關(guān)系 程序，圖形如下：#i nclude <stdio.h>#in elude <math.h>#i nclude <stdlib.h>#defi ne PI 3.141593#defi ne FM 30#defi ne R 3#defi ne KN 200#defi ne

8、XN 101#defi ne ZN 101#defi ne DH 5#defi ne DT 0.001void mai n()不能直接初值為0FILE *fp; int i,j,k,m, n;float u1XNZN,u2XNZN,u3XNZN,u4XNZN,fXNZN; /float 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函數(shù)，當(dāng)且僅當(dāng)i,j同時為50時，f為1,其

9、余為0for(j=0;j<ZN;j+)if(i=5 0&&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;/速度相同表示同一介質(zhì)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)

10、-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+)u1m n=u2m n;u2m n=u3m n;if(k=100)for(m=0;m<XN;m+)for(n=0; n< ZN; n+)u4mn=u3mn; 記錄波前快照，中間點if(fp=fope n( "wavefro nt.dat","w&q

11、uot;)!=NULL)fpri ntf(fp,"%dn",XN);fprin tf(fp,"%dn",ZN);for(i=0;i<XN;i+)for(j=0;j<ZN;j+)fprin tf(fp,"%fn",u4ij);fclose(fp);Wavefronrt of therelicaF vefield第二時刻為地震波已經(jīng)傳播到邊界上的某時刻，體會其人工邊界反射; 程序圖形如下：#i nclude <stdio.h>#in elude <math.h>#i nclude <stdlib.

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

13、T)*(k*DT)*cos(2*PI*FM*k*DT);for(i=0;i<XN;i+)定義f函數(shù)，當(dāng)且僅當(dāng)i,j同時為50時，f為1,其余為0for(j=0;j<ZN;j+)if(i=5 0&&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;/速度相同表示同一介質(zhì)for(k=0;k<KN;k+)for(i=2;i<XN-2;i+)for(j=2;j<ZN-2;j+)uu0=(vij)*

14、(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+)u1m n=u2m n;u2m n=u3m n;if(k=160)/只需改動K值，即顯示邊界反射for(m=0;m<XN;m+)for(n=

15、0; n< ZN; n+)u4mn=u3mn;記錄波前快照，邊界if(fp=fope n( "wavefro nt.dat","w")!=NULL)fpri ntf(fp,"%dn",XN);fprin tf(fp,"%dn",ZN);for(i=0;i<XN;i+)for(j=0;j<ZN;j+)fprin tf(fp,"%fn",u4ij);fclose(fp);4、對于大模型，定義為水平層狀速度模型；做兩個實驗，一是將震源點放在區(qū)域表層任一 點，記錄下某些時刻的波前快照，

16、體會地震波在兩種介質(zhì)的分界面上傳播規(guī)律，指出哪 是反射波，哪是透射波；這時取小模型的常量，為減少頻散，速度v至少為2400程序圖形如下：#i nclude <stdio.h>#in elude <math.h>#i nclude <stdlib.h>#defi ne PI 3.141593#defi ne FM 30#defi ne R 3#defi ne KN 200#defi ne XN 200#defi ne ZN 100#defi ne DH 5#defi ne DT 0.001void mai n()FILE *fp;int i,j,k,m, n;

17、float u1XNZN,u2XNZN,u3XNZN,u4XNZN,u5XNZN; /不能直接初值為 0float fXNZN,vXNZN,wKN,uuO,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;j+) u1ij=O.O;u2ij=0.0;u3ij=0.0;u4ij=O.O;fij=0.0;if (j<=30)vij=2400;第一層速度為 2400elsevij=3000;/

18、第二層速度為 3000f10010=1;定義 f 函數(shù)，在 100, 10 為 1for(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;u

19、5ik=u1i10;/ 地震記錄for(m=0;m<XN;m+)for(n=0;n <ZN; n+)u1m n=u2m n;u2m n=u3m n;if(k=100)for(m=0;m<XN;m+)for(n=0; n< ZN; n+)u4mn=u3mn; 記錄波前快照if(fp=fope n("wavefro nt.dat","w")!=NULL)fpri ntf(fp,"%dn",XN);fprin tf(fp,"%dn",ZN);for(i=0;i<XN;i+) for(j=0;

20、j<ZN;j+)fprin tf(fp,"%fn",u4ij); fclose(fp);Wsvefrorit of heretical 'A'avefield二是合成一個地震記錄，即記錄下與震源同一深度點的各點所有時刻的波場值，并指出記錄上的同向軸分別對應(yīng)哪些波？這時取小模型的常量，為減少頻散，速度v至少為2400程序圖像如下：#i nclude <stdio.h>#in clude <math.h>#i nclude <stdlib.h>#defi ne PI 3.141593#defi ne FM 30#defi

21、 ne R 3#defi ne KN 200#defi ne XN 200#defi ne ZN 200#defi ne DH 5#define DT 0.001 void mai n()FILE *fp;int i,j,k,m, n;float u1XNZN,u2XNZN,u3XNZN,u4XNZN,u5XNZN; /不能直接初值為0float 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<X

22、N;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;第一層速度為 2400elsevij=3100;/ 第二層速度為 3100f10010=1;定義 f 函數(shù)，在 100, 10 為 1for(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*(u