計(jì)算聲學(xué)第一課

1、計(jì) 算 聲 學(xué)1第一章 緒 論 在學(xué)習(xí)了彈性動(dòng)力學(xué)、聲學(xué)和聲測(cè)井理論之后，如何數(shù)值模擬聲場(chǎng)的激發(fā)和傳播？ 如何分析彈性波場(chǎng)的擴(kuò)散、衰減以及波型的轉(zhuǎn)換？ 特別在介質(zhì)分布失去軸對(duì)稱(chēng)時(shí)，純解析方法不適用時(shí)用什么方法來(lái)數(shù)值模擬聲場(chǎng)？ 這些都是計(jì)算聲學(xué)所要解決的問(wèn)題，相關(guān)的具體算法、處理手段和程序構(gòu)成計(jì)算聲學(xué)的主要內(nèi)容。2計(jì)算方法分類(lèi) 根據(jù)方法本身的不同主要分四類(lèi)： 第一 純解析算法（DW）， 第二 攝動(dòng)理論近似求解方法， 第三 半解析方法， 第四 純數(shù)值求解方法（FD）。3 具體根據(jù)介質(zhì)類(lèi)型 可以把每一類(lèi)方法具體展開(kāi)： 第一 理想彈性介質(zhì)， 第二 準(zhǔn)（粘）彈性介質(zhì)， 第三 雙相（多相）介質(zhì)。4 具體根

2、據(jù)聲源類(lèi)型： 點(diǎn)源（單極和多極）， 柱源。5 此外，根據(jù)為分析波場(chǎng)的傳播機(jī)制： 第一 計(jì)算全波聲場(chǎng)， 第二 計(jì)算分波聲場(chǎng)。6 本課程的主要目的掌握基本的聲學(xué)類(lèi)問(wèn)題的數(shù)值模擬程序及有關(guān)分析手段的程序，特別是一些基本的通用子程序。學(xué)會(huì)如何調(diào)試子程序，如何修改和改進(jìn)主程序，達(dá)到能夠針對(duì)具體問(wèn)題獨(dú)立編程的水平。使學(xué)者在上完該課程后的有關(guān)聲學(xué)問(wèn)題的計(jì)算能力得到提高 7第二章 純解析算法-全波場(chǎng)的計(jì)算2.1 單相準(zhǔn)彈性介質(zhì)地層下單極點(diǎn)源激發(fā)的井孔聲波全波場(chǎng)的計(jì)算本節(jié)學(xué)習(xí)在單相（準(zhǔn)）彈性介質(zhì)流體中點(diǎn)源（單極源）激發(fā)與輻射聲場(chǎng)的數(shù)值模擬方法-離散波數(shù)法。8一、設(shè)聲源為定聲壓源，則流體內(nèi)點(diǎn)源聲壓直達(dá)場(chǎng)的計(jì)算 9

3、均勻流體介質(zhì)中的全波FFT 用FALSE 定聲壓點(diǎn)源流體直達(dá)場(chǎng)在空間-頻率域的表達(dá)式為（EW） Notice : The “DFL” in X1=DFL*DSQRT(DBLE(N1-1)10均勻流體介質(zhì)中的全波聲源脈沖的頻域響應(yīng)：FFT 用FALSE 聲源脈沖波列：簡(jiǎn)諧（單色）波的疊加FFT 用TRUE 振幅 11均勻流體介質(zhì)中的全波FFT 用TRUE 定聲壓點(diǎn)源流體直達(dá)場(chǎng)在空間-時(shí)間域的表達(dá)式(DW)為 Notice : The “DFP” in X2=DFP*DSQRT(DBLE(N2-1)12均勻流體介質(zhì)中的全波對(duì)理想彈性介質(zhì)，頻率為實(shí)頻率時(shí)：對(duì)粘彈性介質(zhì)，頻率為復(fù)頻率時(shí)：13Fig.

4、01 Waves propagating in water at a center frequency of 10kHz Comparision between EW and DWZ=1.0mr=0.1mVf=1600m/s14聲源脈沖的類(lèi)型15聲源脈沖的類(lèi)型 PF3余弦包絡(luò) 時(shí)域函數(shù)161 聲源脈沖的類(lèi)型 PF3頻域函數(shù)17 SUBROUTINE PF3(CPT,FR,FI) IMPLICIT DOUBLE PRECISION(A-B,D-H,O-Z) IMPLICIT DOUBLE COMPLEX(C) COMMON/BLKDAT/TC,P0,P0T PAI=3.1415926535D0

5、CI=DCMPLX(0.D0,1.D0) A=P0T+PAI B=P0T-PAI CH=.25D0*TC*(CSIW(CPT+P0T)+CSIW(CPT-P0T) & +.5D0*(CSIW(CPT+A)+CSIW(CPT-A) & +CSIW(CPT+B)+CSIW(CPT-B) 18 正變換： CH=CH*CDEXP(CI*CPT ) 反變換: (我們使用) CH=CH*CDEXP( CI*CPT ) FR=DREAL(CH) FI=DIMAG(CH) RETURN END19其中CSIW為函數(shù)子程序 DOUBLE COMPLEX FUNCTION CSIW(CX) IMPLICIT DO

6、UBLE COMPLEX(C) IMPLICIT DOUBLE PRECISION(A-B,D-H,O-Z) X=DREAL(CX) Y=DIMAG(CX) IF(X.EQ.0.D0.AND.Y.EQ.0.D0) GOTO 111 110 CSIW=CDSIN(CX)/CX GO TO 30 111 CSIW=1. 30 RETURN END20聲源脈沖PF3的調(diào)用方式 DP=PAI2*DFP 角頻率間隔 P0=PAI2*F20 中心角頻率 W0=PAI2*F0 復(fù)頻率虛部對(duì)應(yīng)角頻率 P0T=PAI*TC*F20 TC時(shí)域脈沖波列長(zhǎng) DPT=.5D0*DP*TC WW0=0.5D0*W0*TC

7、 DO J=JFL,JFM PT=J*DPT 對(duì)（角）頻率的循環(huán) CPT=DCMPLX(PT,-WW0) CALL PF3(CPT, FR(J), FI(J) ENDDO21PF3 (NT=3) F20=2kHz, NN=322聲源脈沖的類(lèi)型 PF2變形瑞克子波時(shí)域函數(shù)232 聲源脈沖的類(lèi)型 PF2頻域函數(shù)24 SUBROUTINE PF2(CP,FR,FI)IMPLICIT DOUBLE COMPLEX(C)IMPLICIT DOUBLE PRECISION(A-B,D-H,O-Z)COMMON/DATA/AA,AA2,W0PAI=3.1415926535D0 CI=(0.D0,1.D0)C

8、A1=3.D0*(AA-CI*CP)*2-W0*2CA2=8.d0*AA2*W0*CA1CA3=(AA-CI*CP)*2+W0*2)*3CF=CA2/CA3FR=DREAL(CF)FI=DIMAG(CF)RETURNEND 25聲源脈沖PF2的調(diào)用方式 DP=PAI2*DFP 角頻率間隔 W0=PAI2*F20 中心角頻率 AA=P0/DSQRT(3) AA2=AA*AA W01=PAI2*F0 復(fù)頻率虛部對(duì)應(yīng)角頻率 DO J=JFL,JFM P=J*DP 對(duì)（角）頻率循環(huán) CP=DCMPLX(P,-W01) CALL PF2 (CP, FR(J), FI(J) ENDDO26PF2 (NT=

9、2) F20=2kHz, NN=327快速傅立葉變換子程序FFT28快速傅立葉變換子程序SUBROUTINE COMPFR(A,B,N1,M1,INV)A變換量實(shí)部數(shù)組B變換量虛部數(shù)組 N1,M1變換量數(shù)組元素量及相關(guān)量 INV邏輯變量（.TRUE. OR .FALSE.） 正變換 反變換 A, B既是輸入變量又是輸出變量 29快速傅立葉變換子程序組成及調(diào)用方式 SUBROUTINE COMPFR(XR,XI,N2,M2,.INV.) SUBROUTINE FFT(A,B,N1,M1,KS) SUBROUTINE REORD(A,B,N1,M1,KS,REEL) CALL COMPFR(XR,

10、XI,N2,M2,.TRUE.) 例如： CALL COMPFR(XR,XI,N2,M2,.FALSE.) 例如：30柱函數(shù)子程序I(X) K(X)31柱函數(shù)子程序 SUBROUTINE DCBES(CZ,N,M) CZ自變量(實(shí)部,虛部) N控制變量 N=1計(jì)算I和K；N=其它值只計(jì)算K M1計(jì)算貝塞耳函數(shù)的階數(shù) CZ是輸入變量 輸出變量通過(guò)公用塊導(dǎo)出COMMON /BLKBES/CI(0:5),CII(0:5),CK(0:5),CKK(0:5) 一階導(dǎo)數(shù) 一階導(dǎo)數(shù)32柱函數(shù)子程序組成及調(diào)用方式 SUBROUTINE DCBES(CZ,N,M) SUBROUTINE CPAB 在主程序中 C

11、ALL CPAB 之后在任意位置 CALL DCBES(CZ,N,M)例如：要計(jì)算井壁處流體徑向虛波數(shù)對(duì)應(yīng)的貝塞耳 函數(shù)CZ=DCSQRT(CKZ*CKZ-CKF*CKF)*R CKZ軸向波數(shù)，CKF流體波數(shù)33點(diǎn)源流體直達(dá)場(chǎng)的計(jì)算程序單極點(diǎn)源在流體中的輻射場(chǎng)34直達(dá)場(chǎng)計(jì)算程序SOURCE0.FOR主要流程：計(jì)算聲源及接收器有關(guān)量為計(jì)算貝塞耳函數(shù)作準(zhǔn)備 CALL CPAB計(jì)算波數(shù)CALL CKRI計(jì)算復(fù)頻率、離散點(diǎn)間隔（時(shí)頻、空波）對(duì)頻率作外循環(huán)對(duì)波數(shù)作內(nèi)循環(huán)計(jì)算波數(shù)域的聲源輻射場(chǎng)函數(shù) （注意對(duì)稱(chēng)性），每次內(nèi)循環(huán)后作FFT(波數(shù)空間)35直達(dá)場(chǎng)計(jì)算程序SOURCE0.FOR主要流程：對(duì)頻率循環(huán)

12、完成后即得到頻率空間域的聲源輻射場(chǎng)函數(shù)。并從中取出各接收器處的聲場(chǎng)。作頻率時(shí)間域FFT，并對(duì)復(fù)頻率還原 計(jì)算聲源激發(fā)的瞬態(tài)聲場(chǎng) 主程序有詳細(xì)的中文說(shuō)明36直達(dá)場(chǎng)計(jì)算程序輸入文件SOURCE0.YY $CONS R=0.1D-09, 接收器(場(chǎng)點(diǎn))的徑向坐標(biāo) Vf=1500.D0, 流體波速 Qf=1.d+100, 流體品質(zhì)因素 TC=1.0D-03, PF3聲源脈沖波列長(zhǎng) F10=0.10D+02, 要求接收器與聲源的最大距離 F20=0.20D+04, 聲源的中心頻率 DFL=0.48828125D-01, 線(xiàn)波數(shù)間隔 DFP=0.50D+02, 離散頻率點(diǎn)間隔 ZP=0.1D-06, 第一

13、接收器軸向坐標(biāo) DZ=0.1D+01, 相鄰接收器間隔 $END 37直達(dá)場(chǎng)計(jì)算程序輸入文件SOURCE0.YY $DEVP N1=8193, 波數(shù)空間數(shù)組元素量 N2=8193, 頻率時(shí)間數(shù)組元素量 M1=14, 與N1對(duì)應(yīng) M2=14, 與N2對(duì)應(yīng) NN=3, 控制頻域計(jì)算范圍的參數(shù)之一 NT=2, 聲源脈沖的類(lèi)型選擇2PF2,3-PF3 NW=5, 接收器道數(shù)5 KEY=1, 功能開(kāi)關(guān)。0聲源脈沖及頻譜計(jì)算 1-用公式直接計(jì)算空間頻率域的場(chǎng) 2-用FFT變換由頻率波數(shù)得出頻率空間域的場(chǎng) $END 38直達(dá)場(chǎng)計(jì)算程序輸入文件SOURCE0.YYDFL=1.0/(DDZ*(N-1)DFL =0.48828125D-1；0.244140625D-1,0.1220703125D-1,0.D-1N,M =1025 11；2049 12； 4097 13； 8193 14； 16385 15 39單極點(diǎn)源輻射場(chǎng)的計(jì)算結(jié)果點(diǎn)源在無(wú)限流體中激發(fā)的聲波場(chǎng)40PF2 (NT=2) F20=2kHz, NN=3Vf=1500m/s -EW - DW41PF3 (NT=3) F20=2kHz, NN=3Vf=1500m/s