1、計 算 聲 學1第一章 緒 論 在學習了彈性動力學、聲學和聲測井理論之后，如何數值模擬聲場的激發和傳播？ 如何分析彈性波場的擴散、衰減以及波型的轉換？ 特別在介質分布失去軸對稱時，純解析方法不適用時用什么方法來數值模擬聲場？ 這些都是計算聲學所要解決的問題，相關的具體算法、處理手段和程序構成計算聲學的主要內容。2計算方法分類 根據方法本身的不同主要分四類： 第一 純解析算法（DW）， 第二 攝動理論近似求解方法， 第三 半解析方法， 第四 純數值求解方法（FD）。3 具體根據介質類型 可以把每一類方法具體展開： 第一 理想彈性介質， 第二 準（粘）彈性介質， 第三 雙相（多相）介質。4 具體根

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

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

4、01 Waves propagating in water at a center frequency of 10kHz Comparision between EW and DWZ=1.0mr=0.1mVf=1600m/s14聲源脈沖的類型15聲源脈沖的類型 PF3余弦包絡 時域函數161 聲源脈沖的類型 PF3頻域函數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為函數子程序 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的調用方式 DP=PAI2*DFP 角頻率間隔 P0=PAI2*F20 中心角頻率 W0=PAI2*F0 復頻率虛部對應角頻率 P0T=PAI*TC*F20 TC時域脈沖波列長 DPT=.5D0*DP*TC WW0=0.5D0*W0*TC

7、 DO J=JFL,JFM PT=J*DPT 對（角）頻率的循環 CPT=DCMPLX(PT,-WW0) CALL PF3(CPT, FR(J), FI(J) ENDDO21PF3 (NT=3) F20=2kHz, NN=322聲源脈沖的類型 PF2變形瑞克子波時域函數232 聲源脈沖的類型 PF2頻域函數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的調用方式 DP=PAI2*DFP 角頻率間隔 W0=PAI2*F20 中心角頻率 AA=P0/DSQRT(3) AA2=AA*AA W01=PAI2*F0 復頻率虛部對應角頻率 DO J=JFL,JFM P=J*DP 對（角）頻率循環 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變換量實部數組B變換量虛部數組 N1,M1變換量數組元素量及相關量 INV邏輯變量（.TRUE. OR .FALSE.） 正變換 反變換 A, B既是輸入變量又是輸出變量 29快速傅立葉變換子程序組成及調用方式 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柱函數子程序I(X) K(X)31柱函數子程序 SUBROUTINE DCBES(CZ,N,M) CZ自變量(實部,虛部) N控制變量 N=1計算I和K；N=其它值只計算K M1計算貝塞耳函數的階數 CZ是輸入變量 輸出變量通過公用塊導出COMMON /BLKBES/CI(0:5),CII(0:5),CK(0:5),CKK(0:5) 一階導數 一階導數32柱函數子程序組成及調用方式 SUBROUTINE DCBES(CZ,N,M) SUBROUTINE CPAB 在主程序中 C

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

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

13、接收器軸向坐標 DZ=0.1D+01, 相鄰接收器間隔 $END 37直達場計算程序輸入文件SOURCE0.YY $DEVP N1=8193, 波數空間數組元素量 N2=8193, 頻率時間數組元素量 M1=14, 與N1對應 M2=14, 與N2對應 NN=3, 控制頻域計算范圍的參數之一 NT=2, 聲源脈沖的類型選擇2PF2,3-PF3 NW=5, 接收器道數5 KEY=1, 功能開關。0聲源脈沖及頻譜計算 1-用公式直接計算空間頻率域的場 2-用FFT變換由頻率波數得出頻率空間域的場 $END 38直達場計算程序輸入文件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單極點源輻射場的計算結果點源在無限流體中激發的聲波場40PF2 (NT=2) F20=2kHz, NN=3Vf=1500m/s -EW - DW41PF3 (NT=3) F20=2kHz, NN=3Vf=1500m/s