Page 21 - 201806

第 37 卷 第 6 期 王颖等： 孔隙介质与液体界面上波动方程数值模拟的改进方案——MPML 边界条件的应用 851

 k+ 1 1 [ ]
k
τ 2 = τ k− 2 + ∆t (λ + 2µ) D L x v + λ i,j D L z v k ,
 xx,i,j xx,i,j i,j x,i,j x z,i,j z



 [ ]
 1 1
 k+ k− k k
τ 2 = τ 2 + ∆t λ i,j D L x v + (λ + 2µ) D L z v ,

 zz,i,j zz,i,j x,i,j x i,j z,i,j z



 1 1 [ ]
k+ k− k k
τ 2 1 = τ 2 1 + ∆tµ ∗ 1 1 D L x 1 v + D L z 1 v x , (2)
z
1
1
1
1
xz,i+ ,j+ xz,i+ ,j+ i+ ,j+ x,i+ ,j+ z,i+ ,j+
 2 2 2 2 2 2 2 2 2 2

 ( )

 k+ 1 k+ 1
 k+1 k L x L z
v = v + b 1 ∆t D 1 τ xx 2 + D τ xz 2 ,

1
1
 x,i+ ,j x,i+ ,j x,i+ ,j x,i+ ,j 1
 2 2 2 2 z,i+ , j
 2

 ( )
 k+ 1 k+ 1
 k+1 k L x 2 L z 2
v 1 = v 1 + b 1 ∆t D 1 τ xz + D 1 τ zz ,
z,i,j+ z,i,j+ z,i,j+ x,i,j+ z,i,j+
2 2 2 2 2
其中，D x 、D z 分别表示对 x 和z 的一阶微分算子，i、 而在MPML中，三个正交方向上的衰减因子同为坐
j 为分别为 x 和 z 方向的空间网格位置，k 为时间层 标x的函数，即
位置，∆t 为时间步长。而且对 µ、b x 、b z 等参数的设 (x) (x) (x)
d x = d x (x) , d y = d y (x) , d z = d z (x) ,
置遵循了Zeng等 [10] 改进的真空法，如下所示：
d (x) (x) = p (y/x) (x) (x) ,
d
y
x
µ ∗ 1 = (x) (z/x) (x)
1
i+ ,j+ d (x) = p d (x) , (7)
2 2 z x
1 1 1 1 1 其中，p (y/x) 、p (z/x) 分别为y、z 方向衰减因子的比例
 [ ( )] −1
 + + + ,


 4 µ i,j µ i+1,j µ i,j+1 µ i+1,j+1 系数，通过调整该系数可达到改善计算结果稳定性


( )
if µ i,j µ i+1,j µ i,j+1 µ i+1,j+1 ̸= 0 , (3) 的目的 [9] 。当比例系数均为零时，MPML 方程 (7)



 退化为传统的 PML形式方程 (6)。在高泊松比介质


0, otherwise,
中，MPML能保证数值模拟的稳定性。
b x,i+ ,j =
1
2 2.2 数值模拟和算例分析

2

 , if (ρ i+1,j ̸= 0 or ρ i,j ̸= 0) , 本文中设计了两个模型：具有弹性海底的海洋
ρ i+1,j + ρ i,j (4)
 环境模型以及充填液体的井孔模型。对方程 (2) 迭
 0, if (ρ i+1,j = 0, ρ i,j = 0) ,
代计算，即可得到不同时刻的波场值。采用 MPML
b z,i,j+ 1 =
2 作为吸收边界条件，即在模型的正演模拟中，将模型

2 ( )
 的左右边界及下边界设为MPML层，以有效地压制
 , if ρ i,j+1 ̸= 0 or ρ i,j ̸= 0 ,
ρ i,j + ρ i,j+1 (5) 边界处的人工反射，同时与采用了 PML 层的模型

 0, if (ρ i,j+1 = 0, ρ i,j = 0) ,
正演模拟结果对比，验证MPML在长时间数值模拟
通过式(3)∼式(5)保证了模型上边界处自由边界条 方面的稳定性优势。
件的实现。而且通过该方法，在海水与海底交界处， 图 1 为具有弹性海底的海洋环境模型，海水
对应的液 -固界面处剪应力为零的边界条件也能自 深度 50 m，数值模拟时震源设置在水深 45 m 处
动实现。 (图 1 中星形所示)，在震源右侧距离分别为 500 m、
在波场模拟中由于计算区域是有限的，为了压 1500 m、2500 m、3500 m、5000 m、10000 m 的六个
制截断边界造成的人工边界反射，本论文中采用
0
MPML 作为吸收边界条件。该边界条件作为传统 v p=1500 m/s v s=0
PML边界条件的扩展，在多个正交方向同时引入衰 ງए/m 50
减因子，通过调整不同方向衰减因子的比例系数达 100 v p=3200 m/s v s=1800 m/s
150
到改善 MPML 稳定性的目的。以 x 方向为例，传统 0 50 100 150 200 250 300 350 400 450 500
ᡰሏ/m
的PML方程中，只有一个衰减因子是空间变量x的
函数，即 图 1 含液 -固界面狭长模型 (以源距 500 m 为例)
Fig. 1 Elongated model with liquid-solid interface
d x = d x (x) , d y = 0, d z = 0, (6)
(take source-receiver interval of 500 m for example)