计算声学第一课

2 i t



[
F 2


[ 2 K 0 ( f r ) ]e
i kzz
dk z ] e
d
FFT 用TRUE
1 2
i t

( ) e
d

( ) e
i t
df , df
1 2
d
Notice : The “DFP” in X2=DFP*DSQRT(DBLE(N2-1))
此外，根据为分析波场的传播机制：
第一 计算全波声场，
第二 计算分波声场。
பைடு நூலகம்
本课程的主要目的
★掌握基本的声学类问题的数值模拟程 序及有关分析手段的程序，特别是一些 基本的通用子程序。
★学会如何调试子程序，如何修改和改 进主程序，达到能够针对具体问题独立 编程的水平。
★使学者在上完该课程后的有关声学问 题的计算能力得到提高
计 算 声 学
马 俊 吉林大学物理学院声学与微波物理系
第一章 绪 论
在学习了弹性动力学、声学和声测井理论之 后，如何数值模拟声场的激发和传播？ 如何分析弹性波场的扩散、衰减以及波型的 转换？ 特别在介质分布失去轴对称时，纯解析方法 不适用时用什么方法来数值模拟声场？ 这些都是计算声学所要解决的问题，相关的 具体算法、处理手段和程序构成计算声学的主 要内容。
反变换: (我们使用)
CH=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)
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
★ CALL COMPFR(XR,XI,N2,M2,.TRUE.)
例如：
f (t )
1 2



[ F ( ) ] e
i t
d
★ CALL COMPFR(XR,XI,N2,M2,.FALSE.)
例如：
F ( )



[ f ( t ) ]e
i t
dt
柱函数子程序
I(X) K(X)
一阶导数 一阶导数
柱函数子程序组成及调用方式
★ SUBROUTINE DCBES(CZ,N,M) SUBROUTINE CPAB ★ 在主程序中 CALL CPAB ★ 之后在任意位置 CALL DCBES(CZ,N,M) 例如：要计算井壁处流体径向虚波数对应的贝塞耳 函数CZ=DCSQRT(CKZ*CKZ-CKF*CKF)*R
CALL PF2 (CP, FR(J), FI(J)) ENDDO
2
Palse of the source FFT from the spectrum
Amplitude
0
-2 0.0 0.5 1.0 1.5 2.0
Time(ms)
0.00060 0.00055 0.00050 0.00045 0.00040 0.00035 0.00030 0.00025 0.00020 0.00015 0.00010 0.00005 0.00000 -0.00005 0 1 2 3
柱函数子程序
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)
(k
z
dl , dl
1 2
dk z
Notice : The “DFL” in X1=DFL*DSQRT(DBLE(N1-1))
均匀流体介质中的全波
声源脉冲的频域响应：
FFT 用FALSE
F ( )
FFT 用TRUE
1 2



[ f ( t ) ]e
i t
dt
声源脉冲波列：简谐（单色）波的叠加
均匀流体介质中的全波
★对粘弹性介质，频率为复频率时：
kf
i
V f (1 i 2Q f )
Qf
★对理想弹性
0
介质，频率为 实频率时：
kf

Vf
Comparision between EW and DW
2
°¹ Ê ¼ Ë µ Ö ´ ³ ´ « ½ Æ ã Ä ±ï ¡ °» ²Õ ¿ Ê ¼ Ë µ Ö ´ ³ ´ ý Ö ¹ ª ½ Æ ã Ä ±ï ¡
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 )
★
Amplitude
Two times FFT Spectrum of Palse
4
5
Frequency(kHz)
PF2 (NT=2)
F20=2kHz,
NN=3
快速傅立叶变换子程序
FFT
快速傅立叶变换子程序
SUBROUTINE COMPFR(A,B,N1,M1,INV) A－－－－－－变换量实部数组 B－－－－－－变换量虚部数组 N1,M1－－－－变换量数组元素量及相关量
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 CI=DCMPLX(0.D0,1.D0)
1 2
1 cos t cos 2 f t H ( t ) H ( t t )
2 tc tc 2 tc 2 20 c
1 声源脉冲的类型 PF3
频域函数
sin( 0 . 5 t c 0 . 5 t c 0 ) sin( 0 . 5 t c 0 . 5 t c 0 ) 0 .5t 0 .5t 0 . 5 t c 0 . 5 t c 0 c c 0 sin( 0 . 5 t c 0 . 5 t c 0 ) 0 .5t 0 .5t c c 0 i 0 .5 t sin( 0 . 5 t c 0 . 5 t c 0 ) c F ( ) 0 . 25 t c e 0 . 5 t c 0 . 5 t c 0 0 .5 sin( 0 . 5 t c 0 . 5 t c 0 ) 0 . 5 t c 0 . 5 t c 0 sin( 0 . 5 t c 0 . 5 t c 0 ) 0 . 5 t c 0 . 5 t c 0
4
5
Frequency(kHz)
PF3 (NT=3)
F20=2kHz,
NN=3
2 声源脉冲的类型 PF2 变形瑞克子波
时域函数
f (t )
4 3
t e
2 2 0


0 3
t
sin( 0 t )
2 声源脉冲的类型 PF2
频域函数
3 F ( )
8 3 2 0 0
N11 2
M 11
INV－－－－逻辑变量（.TRUE. OR .FALSE.） 正变换 反变换 A, B既是输入变量又是输出变量
快速傅立叶变换子程序组成及 调用方式
★ SUBROUTINE COMPFR(XR,XI,N2,M2,.INV.)
SUBROUTINE FFT(A,B,N1,M1,KS) SUBROUTINE REORD(A,B,N1,M1,KS,REEL)
f (t )
1 2



[ F ( ) ] e
i t
d
振幅
均匀流体介质中的全波
定声压点源流体直达场★在空间-时间域的 表达式(DW)为 ik r z F e i t 1
2 2 f
PD ( r , z , t )
1 2
2


r

2
z
e
d
★在空间-频率域的表达式为（EW）
PD ( r , z , ) F e r
2 ik
f
r z
2
2
z
2
FFT 用FALSE
1 2


dk z
F 2





[ 2 K 0 ( f r ) ]e
)e
ik z z
i kzz
dk
z

(k
z
)e
ik z z

计算方法分类
根据方法本身的不同主要分四类： 第一 第二 第三 第四 纯解析算法（DW）， 摄动理论近似求解方法， 半解析方法， 纯数值求解方法（FD）。
具体根据介质类型 可以把每一类方法具体展开： 第一 第二 第三 理想弹性介质， 准（粘）弹性介质， 双相（多相）介质。
具体根据声源类型：
点源（单极和多极）， 柱源。
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 －－－－角频率间隔 P0=PAI2*F20 －－－－中心角频率
Z=1.0m
Pressure
r=0.1m
0
Vf=1600m/s
-2 0.0 0.5 1.0
Time(ms)
Fig. 01 Waves propagating in water at a center frequency of 10kHz
声源脉冲的类型
1 声源脉冲的类型 PF3 余弦包络
时域函数
f (t )
0
3
i
2

2
0
3
2

2 0 i 0 3
SUBROUTINE PF2(CP,FR,FI) IMPLICIT DOUBLE COMPLEX(C) IMPLICIT DOUBLE PRECISION(A-B,D-H,O-Z) COMMON/DATA/AA,AA2,W0 PAI=3.1415926535D0 CI=(0.D0,1.D0) CA1=3.D0*(AA-CI*CP)**2-W0**2 CA2=8.d0*AA2*W0*CA1 CA3=((AA-CI*CP)**2+W0**2)**3 CF=CA2/CA3 FR=DREAL(CF) FI=DIMAG(CF) RETURN END
声源脉冲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)
第二章 纯解析算法----全波场 的计算
§2.1 单相准弹性介质地层下单极点 源激发的井孔声波全波场的计算 本节学习在单相（准）弹性介 质流体中点源（单极源）激发与辐 射声场的数值模拟方法----离散波数 法。
一、设声源为定声压源， 则流体内点源声压直达 场的计算
均匀流体介质中的全波
定声压点源流体直达场
1
Palse of the source FFT from the spectrum
Amplitude
0
-1 0.0 0.5 1.0 1.5 2.0
Time(ms)
0.00025
Amplitude
0.00020 0.00015 0.00010 0.00005 0.00000 0 1 2 3
Two times FFT Spectrum of Palse