2004 柱坐标系下井孔声场的有限差分模拟_丛健生

大庆石油学院学报JOU RNA L OF DA Q ING PET ROL EU M IN ST IT UT E 第28卷Vo l .28　第2期N o .2　2004年4月Apr .2004

收稿日期:2003-07-24;审稿人:蔺景龙;编辑:曾湘华 基金项目:黑龙江省自然科学基金项目(A9914)

作者简介:丛健生(1970-),男,实验师,主要从事井中声传播以及超声流量检测等方面的研究.

柱坐标系下井孔声场的有限差分模拟

丛健生1,陈德华1,宋延杰1,吕金玲2,杨大军3

(1.大庆石油学院地球科学学院,黑龙江大庆　163318;　2.临盘采油厂地质所,山东临邑　251507;　3.大庆石油

管理局测井公司,黑龙江大庆　163412)

摘　要:在弹性波动力学基础上,建立了井孔声场在柱坐标系下的速度—应力交错网格高阶差分公式,编制了相关

计算程序.算例结果表明,该模拟程序计算过程稳定,所采用的吸收边界条件吸收效果良好,消除了人为截断界面对声场

模拟的影响;当井外介质存在径向分层时,垂直地层界面对井内声场的影响随其与井壁之间距离的增大而减小,随源距

的增大而增大.

关　键　词:井孔声场;有限差分模拟;交错网格;吸收边界条件;径向分层

中图分类号:P313.3 文献标识码:A 文章编号:1000-1891(2004)02-0011-04

裸眼井中弹性波的数值模拟对石油地球物理勘探,尤其是声波测井具有重要的意义.在模拟井中弹性波场时,可采用解析方法或数值方法;但对于各向异性介质模型,传统的解析方法一般不适用.目前,解决这一问题主要有:混合法、边界积分法、边界元法、有限差分法和有限元法;其中有限差分和有限元法是模拟复杂介质中波场最直接的2种方法,它们可以给出各瞬间波场在计算区域各点处的传播规律[1,2].采用波动方程的微变公式,有限元法可处理各向异性弹性介质和黏弹性介质的复杂模型,其适合小尺度介质中地震波场的数值模拟;但该方法计算量大,而且低阶有限元法具有较弱的频散特性,常规的高阶有限元法又会出现伪波现象[3].有限差分法与有限元法相比要简单、灵活,尽管该方法也存在稳定性和频散方面的问题,但针对不同的问题可提出越来越多的差分格式以提高计算效率.采用优化的交错网格,并取高阶差分来近似微分方程,可使该方法的应用范围更广,同时也降低了空间的频散现象.笔者建立了柱坐标系下速度-应力的有限差分公式,该差分算子精度在空间域上达到8阶,时间域为2阶.应用该算法程序模拟了地层中声场演化过程及井内声波全波列,考察了径向分层地层界面对井内声场的影响.

1　基本原理

1.1　弹性力学公式

对具有轴对称性的完全弹性介质井眼模型,三维模型可简化为二维,此时柱坐标系下用应力和位移所表示的运动方程和本构方程分别为

ρ 2u t

2=1r r (r τrr )+ τrz z ,ρ 2w t 2=1r r (r τrz )+ τzz z ;(1) τrr =(λ+2μ) u r +λu r + w z ,τzz =(λ+2μ) w z +λu r + u r ,τrz =μ u z + w r ;(2)

式中:u ,w 分别为r ,z 方向上的位移分量;τrr ,τzz 分别为r ,z 方向上的正应力;τrz 为剪切应力;ρ为介质的密度;λ,μ为介质的拉梅系数,与地层纵波速度和横波速度的关系分别为v c =(λ+2μ)/ρ,v s =

μ/ρ.在式(1),(2)两端同时对时间求导数,则两式可变为

ρ v r t =1r r (r τrr )+ τrz z

,ρ v z t =1r r (r τrz )+ τzz z ;(3) τrr t =(λ+2μ) v r r +λv r r + v z z , τzz t =(λ+2μ) v z z +λv r r + v r r , τrz t =μ v r z + v z r ;(4)

式中:v r ,v z 分别为r ,z 方向上的速度分量.式(3)和式(4)即为柱坐标系下速度-应力差分方程.这些弹性动力学方程能够控制切过水平和垂直裂缝的、充满流体的井中声场的传播,模型是轴对称的,并且允许介质的弹性在径向和垂向上变化.

1.2　速度-应力差分算法的实现

传统的有限差分法在每一波长内采样点太少时,数值频散,分辨率降低;而采样点增多,又会增加运算量,降低运算速度[4]

.笔者采用高阶差分算子来近似表示空间导数,而时间导数仍采用2阶差分算子来近似,同时,在波场物理量的离散过程中采用交错网格技术来减小频散现象.

按照函数的卷积或Tay lo r 级数展开,函数f (x )对x 的导数可写为如下最佳差分算子[5]:

D +x f =

∑L /2l =1a l [f (x +l Δx )-f (x -(l -1)Δx )]≈ f (x +Δx /2) x ,(5)

D -x f =∑L /2l =1

a l [f (x +(l -1)Δx )-f (x -l Δx )]≈ f (x -Δx /2) x ;(6)式中:a l 为权重系数,是待选择的最佳差分系数,当频带宽度和相对误差已知时可以通过计算得到.L /2为差分算子的长度,当L =2,4,6,8时,式(5)和式(6)分别代表2阶、4阶、6阶和8阶差分算子.

采用交错网格技术,式(3)和式(4)可分别改写为

v n +1/2ri +1/2,j =v n -1/2ri +1/2,j +τn rr r +D -r τn rr +D -z τn rz Δt /ρ,v n +1/2zi ,j +1/2=v n -1/2z i ,j +1/2+τn rz r +D +r τn rz +D +z τn zz Δt /ρ;(7)τn +1rri ,j

=τn rri ,j +(λ+2μ)D -r v n +1/2r +λv n +1/2r r +D -z v n +1

/2z Δt ,τn +1zz i ,j

=τn zzi ,j +(λ+2μ)D -z v n +1/2z +λv n +1/2r r +D -r v n +1/2r Δt ,τn +1rzi +1/2,j +1/2=τn rzi +1/2,j +1/2+μ[D +r v n +1/2z +D +z v n +1/2r ]Δt .

(8)用f 表示式(7),(8)中各速度或应力分量,则f n i ,j =f (

i Δr ,j Δz ,n Δt );其中Δr ,Δz 分别为r ,z 方向上的图1　交错网格中各物理量的位置空间采样间隔,本文中取Δr =Δz ;Δt 为时间采样间隔.D +r ,D -r ,

D +z ,D -z 可由式(5),(6)得到.各应力以及质点速度在交错网格上的分布见图1,不同的物理量定义在不同的位置上.例如,质点的速

度分量定义在应力分量的中间,因此应力的空间导数可用质点的速

度分量迭代得到;并且每个物理量都与相邻的其它物理量联系在一

起.这种交错网格相对常规的中心差分能大大降低空间的频散现

象.交错格式的优点是采用空间和时间的高阶差分算子,在保证计

算精度的基础上,尽可能地采用较大的空间步长,以提高计算效率.1.3　吸收边界条件的选取

有限差分模拟中可选用的吸收边界条件有许多种,但大部分都受声波入射角的影响.本文中采用的阻尼吸收边界条件可以克服这一缺陷,阻尼本身可以理解成在模型边界附近加一个吸收带,波场到达吸收大　庆　石　油　学　院　学　报 第28卷　2004年