计算声学第一课

《计算声学第一课》由会员分享，可在线阅读，更多相关《计算声学第一课（43页珍藏版）》请在装配图网上搜索。

1、计 算 声 学马 俊吉林大学物理学院声学与微波物理系第一章 绪 论 在学习了弹性动力学、声学和声测井理论之后，如何数值模拟声场的激发和传播？如何分析弹性波场的扩散、衰减以及波型的转换？特别在介质分布失去轴对称时，纯解析方法不适用时用什么方法来数值模拟声场？这些都是计算声学所要解决的问题，相关的具体算法、处理手段和程序构成计算声学的主要内容。计算方法分类 根据方法本身的不同主要分四类：第一 纯解析算法（DW），第二 摄动理论近似求解方法，第三 半解析方法，第四 纯数值求解方法（FD）。具体根据介质类型 可以把每一类方法具体展开：第一 理想弹性介质，第二 准（粘）弹性介质，第三 双相（多相）介质。

2、具体根据声源类型：点源（单极和多极），柱源。此外，根据为分析波场的传播机制：第一 计算全波声场，第二 计算分波声场。本课程的主要目的掌握基本的声学类问题的数值模拟程序及有关分析手段的程序，特别是一些基本的通用子程序。学会如何调试子程序，如何修改和改进主程序，达到能够针对具体问题独立编程的水平。使学者在上完该课程后的有关声学问题的计算能力得到提高 第二章 纯解析算法-全波场的计算 2.1 单相准弹性介质地层下单极点源激发的井孔声波全波场的计算 本节学习在单相（准）弹性介质流体中点源（单极源）激发与辐射声场的数值模拟方法-离散波数法。一、设声源为定声压源，则流体内点源声压直达场的计算 均匀流体介质

3、中的全波zzkiFzrikDdkerfKzreFzrPzf 0222)(2 ),(22FFT 用FALSE zzikzzzikzdkdldlekdkekzz21 ,)()(21定声压点源流体直达场在空间-频率域的表达式为（EW）Notice:The“DFL”in X1=DFL*DSQRT(DBLE(N1-1)均匀流体介质中的全波dtetfFit 21)()(声源脉冲的频域响应：FFT 用FALSE deFtfit 21)()(声源脉冲波列：简谐（单色）波的叠加FFT 用TRUE 振幅 均匀流体介质中的全波dedkerfKdezreFtzrPtizzkiFtizrikDzf)(2 ),(0221

4、222122FFT 用TRUE ddfdfedetiti21 ,)()(21定声压点源流体直达场在空间-时间域的表达式(DW)为 Notice:The“DFP”in X2=DFP*DSQRT(DBLE(N2-1)均匀流体介质中的全波ffVk对理想弹性介质，频率为实频率时：对粘弹性介质，频率为复频率时：)21(fffQiVik0fQFig.01 Waves propagating in water at a center frequency of 10kHz Comparision between EW and DW0.00.51.0-202PressureTime(ms)按公式计算的直达场 按

5、积分展开式计算的直达场Z=1.0mr=0.1mVf=1600m/s声源脉冲的类型1声源脉冲的类型 PF3余弦包络 时域函数)()(2coscos1)(2202221ctttttHtHtfttfccc1 声源脉冲的类型 PF3频域函数cticccccccccccccccccccccccccetttttttttttttttttttttttttF5.00000000000005.05.0)5.05.0sin(5.05.0)5.05.0sin(5.05.0)5.05.0sin(5.05.0)5.05.0sin(5.05.05.0)5.05.0sin(5.05.0)5.05.0sin(25.0)(SUB

6、ROUTINE 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 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)正变换：CH=CH*CDEXP(CI*CPT)反变换:(我们使用)CH

7、=CH*CDEXP(CI*CPT)FR=DREAL(CH)FI=DIMAG(CH)RETURN END其中CSIW为函数子程序 DOUBLE COMPLEX FUNCTION CSIW(CX)IMPLICIT DOUBLE 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 END声源脉冲PF3的调用方式 DP=PAI2*DFP 角频

8、率间隔 P0=PAI2*F20 中心角频率 W0=PAI2*F0 复频率虚部对应角频率 P0T=PAI*TC*F20 TC时域脉冲波列长 DPT=.5D0*DP*TC WW0=0.5D0*W0*TC DO J=JFL,JFM PT=J*DPT 对（角）频率的循环 CPT=DCMPLX(PT,-WW0)CALL PF3(CPT,FR(J),FI(J)ENDDO 0123450.000000.000050.000100.000150.000200.00025AmplitudeFrequency(kHz)Two times FFT Spectrum of Palse0.00.51.01.52.0-1

9、01 Palse of the source FFT from the spectrumAmplitudeTime(ms)PF3(NT=3)F20=2kHz,NN=32声源脉冲的类型 PF2变形瑞克子波时域函数)sin()(02203430tettft2 声源脉冲的类型 PF2频域函数3202020230203833)(0iiF 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.

10、D0,1.D0)CA1=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 声源脉冲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)ENDDO PF2(NT=2

11、)F20=2kHz,NN=3012345-0.000050.000000.000050.000100.000150.000200.000250.000300.000350.000400.000450.000500.000550.00060AmplitudeFrequency(kHz)Two times FFT Spectrum of Palse0.00.51.01.52.0-202 Palse of the source FFT from the spectrumAmplitudeTime(ms)快速傅立叶变换子程序FFT快速傅立叶变换子程序SUBROUTINE COMPFR(A,B,N1,M

12、1,INV)A变换量实部数组B变换量虚部数组 N1,M1变换量数组元素量及相关量 INV逻辑变量（.TRUE.OR.FALSE.）正变换 反变换 A,B既是输入变量又是输出变量 11211MN快速傅立叶变换子程序组成及调用方式 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,XI,N2,M2,.TRUE.)例如：CALL COMPFR(XR,XI,N2,M2,.FALSE.)例如：deFtfit 21)()(dtet

13、fFit )()(柱函数子程序I(X)K(X)柱函数子程序 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)一阶导数 一阶导数柱函数子程序组成及调用方式 SUBROUTINE DCBES(CZ,N,M)SUBROUTINE CPAB 在主程序中 CALL CPAB 之后在任意位置 CALL DCBES(CZ,N,M)例如：要计算井壁处流体径向虚波数对应的贝塞耳 函数

14、CZ=DCSQRT(CKZ*CKZ-CKF*CKF)*R CKZ轴向波数，CKF流体波数点源流体直达场的计算程序单极点源在流体中的辐射场直达场计算程序SOURCE0.FOR 主要流程：计算声源及接收器有关量为计算贝塞耳函数作准备 CALL CPAB计算波数CALL CKRI计算复频率、离散点间隔（时频、空波）对频率作外循环对波数作内循环计算波数域的声源辐射场函数 （注意对称性），每次内循环后作FFT(波数空间)直达场计算程序SOURCE0.FOR 主要流程：对频率循环完成后即得到频率空间域的声源辐射场函数。并从中取出各接收器处的声场。作频率时间域FFT，并对复频率还原 计算声源激发的瞬态声场

15、主程序有详细的中文说明直达场计算程序输入文件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,第一接收器轴向坐标 DZ=0.1D+01,相邻接收器间隔$END 直达场计算程序输入文件SOURCE0.YY$DEVP N1=8193,波数空间数组元素量

16、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 直达场计算程序输入文件SOURCE0.YY DFL=1.0/(DDZ*(N-1)DFL=0.48828125D-1；0.244140625D-1,0.1220703125D-1,0.06102515625D-1N,M=1025 11；2049 12；4097 13；819

17、3 14；16385 15 单极点源辐射场的计算结果点源在无限流体中激发的声波场01234-2.00E+0070.00E+0002.00E+007Y Axis TitleTime(ms)z=10-7,r=10-10m01234-202z=1m01234-101z=2m01234-0.50.00.5z=3m01234-0.250.000.25z=4mThe radiating acoustic wave field of monipolePF2(NT=2)F20=2kHz,NN=3Vf=1500m/s -EW -DWPF3(NT=3)F20=2kHz,NN=3Vf=1500m/s -EW -DW

18、01234-2.00E+0070.00E+0002.00E+007Time(ms)z=10-7,r=10-10m01234-202z=1m01234-101z=2m01234-0.50.00.5z=3m01234-0.250.000.25z=4mThe radiating acoustic wave field of monipole单极点源辐射场的计算作业 1、点源脉冲PF3的考察 F20=10,15kHz,NN=1,2,3,4,5，TC=0.5ms,1.0ms,2.0ms 2、点源脉冲PF2的考察 F20=5,20kHz,NN=1,2,3,4,5 3、直达场的计算 任选一种声源脉冲：计算F20=20kHz，(1)ZP=1m,dz=0.5m (2)ZP=10m,dz=1m 4、思考题：如何根据要求适当选取DFL,DFP？