Page 60 - 《应用声学》2020年第1期

56 2020 年 1 月

其中， 并行计算集群，结合上述并行计算及网格划分方法，

包含仪器提升的偶极三维远探测波场模拟的大计
 C 11 = λ + 2µ,


 算量问题得到了有效解决。
C 12 = λ, (3)



C 44 = µ. 2 数值计算

式 (1) ∼ (2) 中，v i (i = x, y, z)、τ ij (i, j = x, y, z)
2.1 溶洞模型
和 g ij (i, j = x, y, z) 分别为速度分量、应力分量和
首先建立包含井孔的三维溶洞模型。图 1 为三
体积源分量。式 (3) 中，λ 和 µ 是介质的拉梅系数，
维数值模拟模型的xz 方向切面示意图，模型大小为
与纵横波速度 V p 和 V s 的关系为 λ = ρ(V − 2V )，
2
2
s
p
2
µ = ρV 。对于液体介质，比如井孔流体和储层流 x = 16 m、y = 8 m、z = 16 m，井孔半径为0.1 m，采
s
体，被认为是一种具有零剪切速度的特殊各向同性 用非均匀网格，在井孔附近网格大小为 0.01 m，地
层中网格大小为0.02 m，时间步长为16.0 × 10 −3 s。
介质。速度 -应力方程采用交错网格在空间 -时间域
如图 1 所示，溶洞为一球体，球内充满水，地层为碳
上进行离散。各速度分量和应力分量交错分布在空
酸盐岩地层，地层参数如表 1所示，球心距离井轴为
间网格上，相邻分量之间相差半个网格；在时间上各
L，球体直径为d。
分量也是交错分布的，相邻的速度分量和应力分量
之间相差半个时间步长。以 v x 和 t xx 为例，其离散
方程为 [9]
ڡࡏ
∆t ( n
n+1/2
n−1/2
v = v + δ x t
x i+1/2,j,k x i+1/2,j,k xx i,j,k
ρ i+1/2,j,k
̌ߘ d
)
+ δ x t n + δ x t n , (4)
xy i+1/2,j+1/2,k xz i+1/2,j,k+1/2 L
[
δ
τ n+1 = τ n + ∆t (λ + 2µ) i,j,k x v n+1/2
xx i,j,k xx i,j,k x i+1/2,j,k
]
+ λ i,j,k δ y v n+1/2 + λ i,j,k δ z v n+1/2 , (5)
y i,j+1/2,k z i,j,k+1/2
图 1 溶洞模型示意图
其中，δ i (i = x, y, z) 是差分算子，∆t 是时间步长。 Fig. 1 Schematic diagram of karst cave model
为了减弱数值频散，空间和时间步长应该满足：
表 1 地层参数
max(∆x, ∆y, ∆z) < C min /10f max , (6)
Table 1 Formation Parameters
√
∆t < 1.0/C max ∆x −2 + ∆y −2 + ∆z −2 , (7)
纵波速度/ 横波速度/ 密度/
其中，C min 和 C max 分别是计算模型的最小速度和 介质 (m·s −1 ) (m·s −1 ) (kg·m )
3
最大速度，f max 是声源的最高频率，∆x、∆y 和 ∆z 井孔流体 1500 — 1000
为 x、y 和 z 方向的空间步长。本文中采用非分裂完 裂缝 1500 — 1000
全匹配层作为吸收边界条件 [10] ，无需分解波场的各
溶洞 1500 — 1000
个变量，而是直接进行复数扩展代换，虽然在模拟中
地层 5200 3200 2700
引入了卷积运算，但是各个物理量的完整性得到了
保证。 图 2 是距离井轴 4 m、直径为 1.2 m 的溶洞 SH
由于三维计算模型较大，整个有限差分计算采 反射波场，其中图 2(a)、图 2(b) 和图 2(c) 为共偏移
用非均匀交错网格并行算法来提高计算效率。在井 距道集 (COG)，D 为源距；图 2(d)、图 2(e) 和图 2(f)
孔附近对网格进行精细划分，再在一般地层介质中 为共源道集 (CSG)，P 为声源提升高度。从反射波
则使用大网格以节约计算资源 [11] 。采用消息传递 到时分析可以得出，红色方框内能量较弱的为来自
和共享存储的方式，进程间采用消息传递并行编程 溶洞近井壁界面的 SH 反射波，后续能量较强的为
MPI，各进程内部采用共享内存的 OpenMP 并行， 来自溶洞远井壁界面的反射波，初至 SH 反射波能
这样的设计可以使得共享了进程内存地址空间的 量较弱，由于溶洞尺寸较小，跟后续振幅较强波场
多个线程共同执行相应进程的指令。基于天河二号 较难区分(主要是指球内聚焦波，见下文)，这也会给