计算声学第一课

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