地球物理学报  2012, Vol. 55 Issue (8): 2787-2797   PDF    
自适应hp-FEM在随钻电阻率测井仪器响应数值模拟中的应用
李辉 , 刘得军 , 刘彦昌 , 郭智勇 , 邢晓楠     
中国石油大学(北京)地球物理与信息工程学院,北京 102249
摘要: 随钻电阻率测井仪器响应数值模拟在电法测井仪器研制和资料处理与解释中具有重要的意义.本文利用自适应hp有限元算法(hp-FEM)模拟和分析了地层各向异性以及相对倾斜对仪器接收信号幅值比和相位差的影响,并对仪器响应特性进行了评价;同时,利用HERMES有限元软件包与VC++开发平台相结合,可以对随钻电阻率测井仪器响应数值模拟过程中,仪器产生的感应电场在地层中和地层分界面处的分布进行显示,从而为钻井地质导向提供有效的信息.计算结果表明,自适应hp-FEM算法收敛速度快且计算精度高,非常适合于对随钻电阻率测井仪器响应进行数值模拟计算.
关键词: 随钻电阻率测井      数值模拟      自适应hp有限元      HERMES      幅值比      相位差     
Application of self-adaptive hp-FEM in numerical simulation of resistivity logging-while-drilling
LI Hui, LIU De-Jun, LIU Yan-Chang, GUO Zhi-Yong, XING Xiao-Nan     
College of Geophysics and Information Engineering, China University of Petroleum, Beijing 102249, China
Abstract: Numerical simulation of resistivity logging-while-drilling (LWD) is important to the electrical well logging instrument design and logging data applications. In this paper, self-adaptive hp-FEM method is used to simulate and analyze the anisotropy of the stratum, and the influence of different dipping angles which can change the amplitude ratio and phase difference of the receive signal from the instrument, the response characteristic is also evaluated. Meanwhile, the finite element analysis software HERMES and VC++ development platform are combined as a tool to reflect the electric field distribution in the stratum and at the boundary, which can provide the effective information in drilling geosteering. Numerical illustration results show that the self-adaptive hp-FEM method has the ability of fast convergence and high computing accuracy. Thus, self-adaptive hp-FEM is an effective method for numerical simulation of electrical logging..
Key words: Resistivity logging-while-drilling      Numerical simulation      Self-adaptive hp-FEM      HERMES      Amplitude ratio      Phase difference     
1 引 言

在石油钻井领域,随钻电阻率测井技术已经成为实时井场数据采集、解释和现场决策以及指导地质导向钻井的关键技术,由于其在地质导向、实时地层对比评价中具有重要应用价值而得到广泛应用[1].如今,随钻电阻率测井技术已成为开发复杂油藏的关键技术,同时对随钻电阻率测井仪器测量精度提出了更高要求.依据随钻电阻率测井仪器响应模拟得到的结果建立系统、准确的随钻电阻率测井资料解释理论,对开发高性价比随钻电阻率测井仪器具有重要意义[2-3].

目前,随钻电阻率测井已成为开发垂直井和大斜度井中重要的测井方法,开展随钻电阻率测井仪器响应数值模拟是研制高精度随钻电阻率测井仪器的重要手段.近年来,国内外有关学者对随钻电阻率仪器响应数值模拟进行了大量的研究并提出了许多卓有成效的数值模拟算法.例如,Zhou 等[4]研究了随钻电阻率测井仪器在水平多层介质中的解析解及数值计算结果;Lovell等[5]提出将有限元算法应用于随钻电阻率测井响应研究;Jing等[6]利用柱坐标传输线矩阵法对随钻电阻率测井仪器响应进行了模拟并取得了较好的效果;Zhang等[7]在应用共轭梯度-快速傅里叶变换算法的基础上提出了使用稳定型双共轭梯度-快速傅里叶变换算法来对三维非均匀介质的电场情况进行数值模拟;Wang 等[8]采用了时域有限差分法模拟随钻电阻率测井仪器响应并获得了一定的应用效果;魏宝君等[9]在柱坐标系下推导出径向成层介质中磁流源并矢格林函数的矢量本征函数展开式,从而得出研究随钻电阻率测井仪器的电磁场分布、接收线圈响应、刻度方法所必需的解析计算式;刘福平等[10]使用传输线矩阵法对过套管电阻率测井仪器响应进行了研究;邢光龙等[11]对电磁传播电阻率测井仪的响应函数进行了分析,系统地考察了其接收信号的相位差和幅度比随不同电导率和仪器参数的变化规律;Nie等[12]提出采用体面线积分方程法来提高井眼条件下随钻电阻率测井响应的计算效率和建模准确度.以上随钻电阻率测井仪器响应数值模拟算法都取得了一定的应用效果.

有限元(FEM)算法是一种研究工程结构问题的有效计算方法,适用于具有复杂边界形n或边界条件、含有复杂媒质的问题,在分析井眼周围复杂地层结构时具有明显的优势.自适应有限元算法(hp-FEM)是从经典有限元算法中派生出来的一种现代有限元算法,由Babuˇska[13]等学者于20世纪80年代首次提出;1998年Demkowicz等[14]提出使用hp-FEM 算法对电阻率测井仪器响应进行模拟并得到了二维电阻率测井响应模拟曲线;2007 年Pardo等[15]结合多重网格技术对hp-FEM 进行了改进,提出了自适应hp-FEM 算法,模拟了简单边界条件下地层电场的分布情况并得到电阻率测井响应仿真曲线;2011年陈晓晖等[16]利用高精度自适应hp-FEM算法模拟和分析了水平层n各向同性均匀地层中随钻电阻率测井电场分布情况并取得了一定的应用效果.由于自适应hp-FEM 算法能够有效地提高有限元分析的效率及可靠性,从而为许多工程技术问题开辟了道路,对于随钻电阻率测井仪器响应数值模拟,自适应hp-FEM 算法具有一定的优势.当井眼周围地层倾斜时,由于井周地质条件的复杂性以及地层几何结构的多变性,就需要数值模拟模型能够进行自适应的调整.由于自适应hp-FEM 算法包含多种技术,主要有:自适应网格调整、非线性问题中载荷增量的自适应选取以及瞬态问题中时间步长的自适应调整,等.因此,自适应hp-FEM 算法适合于对随钻电阻率测井仪器响应进行数值模拟计算.同时与其他随钻电阻率测井仪器响应数值模拟算法相比,自适应hp-FEM 算法还具有收敛速度快、计算精度高、鲁棒性好等特点.因此,本文基于HERMES有限元软件包[17]和VC+ + 开发平台,将自适应hp-FEM 算法应用于随钻电阻率测井仪器响应[18-20]数值模拟.利用自适应hp-FEM 算法计算和分析了随钻电阻率测井仪器响应特性,同时考察了地层倾角、各向异性以及接收线圈倾斜对测井仪器响应的影响.

2 原理和方法 2.1 时谐麦克斯韦方程

假设随钻电阻率测井仪器工作在高频环境下,为单发双收结构且发射源频率为2 MHz,如图 1所示.其中,T为发射线圈,R1R2 为接收线圈,TR1 的距离为L1R1R2 之间的距离为L2.发射线圈和接收线圈的半径均为r,匝数为ns为接收线圈面积,发射线圈通以交变电流IT=Ieiωt.由于接收线圈感生电动势是由交变电流产生,即麦克斯韦方程组中的场量是单频的谐振函数,因此求解的电磁场为时谐场.由于发射线圈内部存在交变电流,在地层中必然产生随时间变化的磁场,所以接收线圈附近的感应磁场为时谐电磁场.根据电磁场原理,麦克斯韦方程可表示为

图 1 随钻电阻率测井仪器模型 Fig. 1 Model of LWD instrumen

(1)

其中,HE分别代表磁场和电场;Mimp是外加磁场;εμσ 分别代表介电常数、磁导率和介质电导率;ρ 代表场强分布;ω (ω≠0)代表发射源角频率;Jimp代表外加电流.将式(1)中法拉第定律代入安培定律,可得时谐麦克斯韦方程在求解域Ω 的波动方程,如(2)式所示,

(2)

在矢量场中,任意旋度场一定是无散场,存在

(3)

(4)

将(1)式中高斯电场定律和(3)式代入(4)式可得,

(5)

则(5)式为电场连续方程.

2.2 边界条件

(1) 理想导体边界条件理想导体是一种理想化的高导电性媒质,存在σ → ∞.依据安培定律可知,在理想导体某一封闭区域内,存在E→ 0.同理,在时变环境下,存在磁场H→0.理论上,任何边界处电场强度的切向分量及磁感应强度的法向分量都是连续的,因此在理想导体上不存在电场切向分量及磁场法向分量.理想导体导电壁表面边界条件定义为

(6)

(2) 发射线圈边界条件

利用外加源电流Jimp 建立发射线圈模型.依据等效原理[15],可将外加源电流Jimp 等效为面电流Jsimp ,以此来建立发射线圈模型.发射线圈边界条件定义为

(7)

(3) 非理想导体边界条件

与理想导体不同,非理想导体存在损耗且σ≠0.因此,非理想导体导电壁表面边界条件定义为

(8)

2.3 变分方程

在高频情况下采用矢量有限元法在H(curl)空间对电场进行求解,可以保证各单元相邻表面之间电场切向的连续性,以及电场旋度零空间的正确表示.同时,在H(curl)空间对电场进行求解,能够使求解过程中不会出现伪解.H(curl)空间定义为

(9)

假设随钻电阻率测井问题的有限元求解区域为Ω,求解区域边界∂Ω 为理想边界,则

(10)

其中,F 为任意测试函数;当,其中,E1E2 为接收线圈处电场强度,x1x2 为计算域,由于求解区域边界是理想边界,则求解电场E所用的矢量求解空间H可定义为

(11)

将封闭区域边界条件作为约束方程代入波动方程(2)式,同时再依据法拉第定律,并引入任意测试函数F,在Ω 域积分可得

(12)

再依据安培定律,

(13)

其中,,为在Г2 边界条件下任意测试函数F 沿切线方向的矢量,n为单位法线矢量.将(13)式代入(12)式可以得到H空间下基于时谐麦克斯韦方程的有限元标准变分方程:

(14)

其中,Г1 为理想导体边界条件;Г2 为发射线圈边界条件.当外加磁场Mimp为零时,(14)式可写成一个双线性函数和一个线性函数的组合,即:

(15)

(16)

(17)

利用伽辽金法可将求解空间H离散为若干个有限维子空间,并对其分别进行求解组成各自的刚度矩阵,最后将各刚度矩阵组集为总刚度矩阵,其矩阵形式为

(18)

其中,A为总刚度矩阵,B为施加条件数,x为需要求解的未知变量.由于(18)式为大型稀疏矩阵,而自适应hp-FEM 算法在求解区域内的网格剖分是自适应的且求解区域较大,因此其施加条件数很大.此时,为提高计算速度,本文基于HERMES有限元软件包,利用UMFPACK[19-20]求解器对(18)式进行求解.

用初始网格τhpΩ 划分为多个互不重叠的单元k1k2,…,kn.其中,各单元尺寸为h1h2,…,hn>0,多项式阶次为p1p2,…,pn≥1.此时在求解空间H下有限元求解子空间定义为

(19)

其中,Ehpki中依据初始网格得到的电场E的近似解.由H(curl)空间特性可知,此时空间Hhp中的矢量场分布是不连续的,但该场量的切向分量在整个求解区域各单元表面上保持连续变化,满足高频情况下求解随钻电阻率测井过程中电场分布的要求.

2.4 自适应hp-FEM 算法求解过程

自适应hp-FEM 算法求解过程如下:

(1) 设初始网格τhp由多个互不重叠的单元组成,求解过程中最大容许误差MPE 大于0,每步hp自适应能增加的最大自由度数为D

(2) 选择一个初始网格K作为待细化对象,计算此时电场近似解EhpHhp

(3) 计算参考解ErefHrefHref=Hh/2,p+1

(4) 将k细化为多个子单元ki(i=1,2,…,n),根据约束情况和最小规则对每条边的多项式阶次进行调整;

(5) 对网格中每个子单元计算其能量范数,根据能量范数和Eref构造子单元误差公式ERRi.然后根据子单元误差求解全局误差ERR.其中,

(6) 当ERR≤MPE 时,停止计算转入后处理.否则将所有子单元按ERRi递减顺序加入到序列l中;

(7) 如果该步增加的自由度数量小于D,则在l中选取误差最大的单元k(l中排列第一的单元)进行hp细化;

(8) 继续第(2)步,直至满足ERR≤MPE.

2.5 建模过程

自适应hp-FEM 算法可根据模型的实际情况不同选择多种方法进行细化:(1)直接在初始网格的各单元中插入一些节点,在空间上将每个单元分裂成多个子单元,这种方式称为h细化.h细化适用于场域分布变化剧烈的区域,无需节点间连接信息,易于实现;(2)不改变初始网格的各单元尺寸,只提高其多项式阶次p,这种方式称为p细化.p细化适用于场域分布变化平滑、网格密集的区域,它的自由度都集中在节点上,无需考虑元素间的协调性;(3)将hp细化相结合,根据地层的实际情况和对求解区域中场量分布求解得到的误差指示选择分别或同时使用hp细化策略对网格进行剖分密度的调整.本文利用HERMES 有限元软件包,并依据自适应hp-FEM 算法建立随钻电阻率测井仪器响应数值模拟模型,建模过程如图 2所示.

图 2 建模过程 Fig. 2 Modeling process
3 数值算例

随钻电阻率测井仪器工作原理是利用电磁波在地层介质中传播时,不同介质对电磁波的吸收差异,将不同地层的电导率大小转变为接收线圈之间感应电动势的幅值比及相位差,再根据输出响应曲线来解释各地层的电阻率大小.随钻电阻率测井仪器的基本结构如图 1所示,由两接收线圈和一个发射线圈所组成.发射线圈T与接收线圈R1 的距离为L1=0.35m (源距),接收线圈R1 与接收线圈R2 之间的距离为L2=0.16m (间距),发射频率f=2 MHz,发射线圈和接收线圈带有磁缓冲器,磁缓冲器电阻率为1×104 Ωm,钻杆电阻率为1×106 Ωm,钻杆周围充满泥浆,泥浆电阻率为0.2Ωm.通过考察接收线圈R1R2 所接收到电磁波信号的幅值比和相位差来反映地层电导率的大小.其中,幅值比(VEATT)和相位差(ΔФ)定义为

(20)

3.1 模型Ⅰ

模型Ⅰ为各向同性均匀地层模型,随钻电阻率测井仪器结构和地层模型如图 1 所示.设三层地层(上、下围岩和目的层)均为各向同性,各层水平分量和垂直分量电导率相等,上围岩(设为低阻层)电导率为σ0=2S/m,目的层(设为高阻层)电导率σ1=0.1S/m,下围岩(设为低阻层)电导率σ2=2S/m.该各向同性均匀地层模型其水平方向厚度为Lh=2.773m,垂直方向深度为Lv=7.6 m,目的层在垂直方向层厚为2m.当接收线圈无倾斜时,考察仪器穿过地层倾角α 分别为0°,30°,45°和60°的地层时,两接收线圈所接收到电磁波信号的相位差和幅值比的变化情况.图 3为仪器行进至目的层,地层垂直方向深度Lv=3.3 m,地层倾角为30°时依据自适应hp-FEM 算法建立模型其电场分布情况;图 4 为依据自适应hp-FEM 算法对目的层部分求解区域进行网格剖分、细化和定阶.限于篇幅,不再给出基于自适应hp-FEM 算法,当α 为0°、45°和60°仪器行进至目的层时的电场分布与网格细化及阶次图.图 5为仪器穿过地层模型时,两接收线圈所接收到电磁波信号的幅值比和相位差曲线.

图 3 电场分布(α = 30°) Fig. 3 Distribution of electric field (α = 30°)
图 4 网格细化及阶次图(α=30°) (a)网格细化;(b)阶次分布. Fig. 4 Illustration ofmesh (α = 30°) Mesh subdivision; (b) Order distribution.
图 5 模型I (a)幅值比和(b)相位差θ=0°) Fig. 5 Amplitude ratio (a) and phase difference (b) formodel I (θ=0°)

通过图 345可以看出,当地层倾角较低,仪器穿过地层分界面时电场变化不大且连续性较好,地层分界面处网格剖分稀疏,幅值比和相位差曲线平滑;当地层倾角逐渐变大,仪器穿过地层分界面时电场变化较强且出现局部电场不连续,此时幅值比和相位差曲线同样会出现波动,仪器响应表现为角峰,也就是极化角.但无论仪器由上围岩行进至目的层或由目的层行进至下围岩,所测幅值比和相位差曲线表现出的极化角类型相同,此时只能判断地层边界的存在却无法准确判断仪器穿越地层为高阻层还是低阻层.同时,由于极化角宽度较窄,只有当仪器靠近地层分界面甚至穿过分界面时才会出现极化角,因此无法及时提供地质导向信息.

当接收线圈按一定角度倾斜时(θ=45°),考察地层倾角α 分别为0°,30°,45°和60°仪器穿过地层模型时,两接收线圈所接收到电磁波信号的幅值比和相位差变化情况.图 6为接收线圈倾斜45°角仪器以不同角度穿过地层模型时,两接收线圈所接收到电磁波信号的幅值比和相位差曲线.

图 6 模型Ⅰ(a)幅值比和(b)相位差(θ=45°) Fig. 6 Amplitude ratio(a)and phase difference(b)for modelI(θ=45°)

通过图 5图 6可以看出,接收线圈有倾角和无倾角时,两接收线圈之间的幅值比和相位差会出现明显差异.通过图 6可以看出,接收线圈倾斜45°时仪器响应表现为,当仪器从低阻层进入高阻层时,幅值比和相位差曲线无极化角出现而仪器从高阻层进入到低阻层时幅值比和相位差曲线出现既宽又深的极化角,此时如果仪器以相对倾角60°从高阻层穿越至低阻层,由图 6 中幅值比曲线可知所产生的极化角沿仪器行进方向最大宽度约为0.5m,

如果仪器沿垂直方向行进0.5 m,则仪器的实际穿越厚度约为0.5/cos60°=1m[21].利用地层分界面处极化角出现的不对称性,可以判断仪器在行进过程中钻入地层为低阻层或高阻层.

为了验证自适应hp-FEM 算法在该模型下的精确性和快速性,对于模型Ⅰ,本文再利用自适应h-FEM 算法和自适应p-FEM 算法计算接收线圈无倾角且地层倾角α 分别为0°,30°,45°和60°时,两接收线圈所接收到电磁波信号的幅值比和相位差随地层电导率变化的关系.

图 7a为仪器行进至目的层,地层垂直方向深度Lv=3.3 m,地层倾角为30°时依据自适应h-FEM算法建立模型对目的层进行网格剖分,图 7b为依据自适应p-FEM 算法建立模型对目的层进行网格定阶.以下不再给出基于自适应h-FEM 和自适应p-FEM 算法,α为0°、45°和60°仪器行进至目的层时的电场分布和网格细化以及网格阶次图.

图 7 >(a)h-FEM 建模和(b) p-FEM 建模(α=30°) Fig. 7 h-FEM model (a) and p-FEM model (b) (α=30°)

表 1为自适应hp-FEM、自适应h-FEM 和自适应p-FEM 三种算法在接收线圈无倾斜且地层倾角在30°时的仪器响应.其中,z代表地层深度,D代表自由度,T代表计算时间,N代表迭代次数,E代表全局误差.通过表 1 可以看出,自适应hp-FEM 算法所需要的自由度范围从19052 个至23814 个,平均使用自由度为21277个;计算时间范围从1.61s至2.54s,平均计算时间约为2.13s;迭代次数最大为13 次;全局误差最大值为0.899%,最小值为0.534%,全局误差平均值约为0.734%.自适应h-FEM 算法平均使用自由度为55782个,平均计算时间约为3.54s,全局误差平均值约为0.839%.自适应p-FEM 算法平均使用自由度为19539 个,平均计算时间约为2.77s,全局误差平均值约为0.82%.在地层分界面处(2.0m 和4.0 m)自适应hp-FEM算法所需要的自由度、计算时间、迭代次数及全局误差都偏大,其他两种算法表现相同,主要原因在于感应电流在穿过地层交界面时由于界面两侧的地层介质电导率不同,根据欧姆定律,界面两侧的电场必然不同,界面处电场垂向分量不连续造成界面电荷积累,从而仪器响应表现为在地层交界面处出现极化角,同时在地层交界面处为保证计算精度自适应网格剖分更细多项式阶次更高从而导致计算时间、自由度和迭代次数的增加.与自适应h-FEM 和自适应p-FEM 算法相比,自适应hp-FEM 算法具有较快的计算速度和较高的计算精度.因此,自适应hp-FEM 算法优于其他两种算法.以下本文将利用自适应hp-FEM 算法对目的层为各向异性和上下围岩及目的层均为各向异性的两种模型进行探讨.

表 1 仪器响应 Table 1 Instrument response
3.2 模型Ⅱ

模型Ⅱ为上、下围岩各向同性目的层为各向异性地层模型.其中,上、下围岩为各向同性且电导率相同,σ0=σ2 =2S/m,视为低阻层.目的层为各向异性,水平电导率为σh1=0.2S/m,垂直电导率为σv1=0.04S/m,视为高阻层,该地层模型水平方向厚度为Lh=2.773m,垂直方向深度为Lv=7.6m,目的层在垂直方向层厚为2m,仪器参数与模型Ⅰ相同.当仪器接收线圈无倾斜时,考察仪器穿过地层倾角α 分别为0°,30°,45°和60°的地层模型时,两接收线圈所接收到电磁波信号的幅值比和相位差,如图 8a图 8b所示.利用三种算法(自适应hp-FEM、自适应h-FEM 和自适应p-FEM)对模型Ⅱ进行数值模拟,限于篇幅对目的层电场分布、网格细化及网格定阶图不再给出.

图 8 模型Ⅱ (a)幅值比和(b)相位差(θ=0°) Fig. 8 Amplitude ratio (a) and phase difference (b) formodel Ⅱ(θ=0°)

通过图 8可以看出,随着地层倾角的增大幅值比和相位差呈下降趋势,仪器响应曲线表现为沿仪器行进方向地层电导率作用较为强烈.当地层倾角大于45°时相位差曲线在地层分界面处出现极化角,地层倾角大于60°时极化角较为明显.同时,电场在地层分界面处沿地层分界面分布较密,而在离地层分界面较远的地方电场变化平缓且出现减弱趋势,其主要原因在于,电磁波在地层中传播时会产生衰减,地层电阻率越低,衰减速度越快.当接收线圈倾斜45°,利用自适应hp-FEM 算法考察仪器穿过地层倾角α 分别为0°,30°,45°和60°的地层时,两接收线圈所接收到电磁波信号的幅值比和相位差,如图 9a9b所示.

通过图 8图 9可以看出,当接收线圈倾斜45°时,两接收线圈之间的幅值比和相位差出现差异.图 9中幅值比和相位差都出现极化角而图 8a中极化角则不太明显.当感应电流穿过地层交界面时,因界面两侧介质电导率不同,导致界面两侧感应电场出现异常,由于目的层存在各向异性且受垂直电导率作用明显,因此在地层分界面处电荷会不断积累从而形成极化层;当仪器接收线圈经过地层分界面时会受到极化层的作用,从而导致接收线圈接收到的电磁波信号出现失真,此时接收到电磁波信号的幅值比和相位差曲线就会出现极化角,并随地层倾斜角的不断增大而发生改变.同时,由于电磁波在地层中传播时会产生衰减,在低阻层中衰减较快在高阻层中衰减较慢,从而导致界面两侧感应电场的强度发生明显改变;而仪器接收线圈倾斜一定角度后仪器源距将产生一定程度的变化,源距越大受地层倾斜程度的影响就越大,必然在一定程度上增强了极化角的幅度.因此,当仪器从低阻层进入高阻层时极化角较浅,而当仪器从高阻层进入低阻层时极化角则十分明显.

图 9 模型Ⅱ(a)幅值比和(b)相位差(θ=45°) Fig. 9 Amplitude ratio (a) and phase difference (b) formodel Ⅱ (θ=45°)

为了验证自适应hp-FEM 算法在该模型下的精确性和快速性,利用自适应h-FEM 算法和自适应p-FEM 算法计算接收线圈无倾角且地层倾角α为45°时的仪器响应.通过计算可得,自适应hp-FEM 算法所需要的自由度范围从18810 个至21605个,平均使用自由度为22123个;计算时间范围从1.46s至2.29s,平均计算时间约为1.99s;迭代次数最大为10次;全局误差最大值为0.98%,最小值为0.681%,全局误差平均值约为0.734%.自适应h-FEM 算法平均使用自由度为92105个,平均计算时间约为2.44s,全局误差平均值约为0.84%.自适应p-FEM 算法平均使用自由度为31200 个,平均计算时间约为2.01s,全局误差平均值约为0.93%.三种算法在地层分界面处所需要的自由度、计算时间、迭代次数及全局误差都会显着增大,主要原因不仅因为交界面两侧的地层介质电导率不同,同时与目的层存在各向异性,其水平和垂直方向电导率不同有关.由于目的层存在各向异性,随着垂向深度的增加,三种算法表现出的仪器响应其各项参数都随之而变化,但都表现为沿垂直方向电导率作用较为强烈.

3.3 模型Ⅲ

模型Ⅲ为上、下围岩及目的层都为各向异性地层模型.上围岩水平电导率σh0=2S/m,垂直电导率σv0=0.5S/m,视为低阻层;目的层水平电导率σh1=0.1S/m,垂直电导率σv1= 0.01S/m,视为高阻油层;下围岩水平电导率σh2=2S/m,垂直电导率σv2=0.5S/m,视为低阻层;该地层模型水平方向厚度为Lh=2.773m,垂直方向深度为Lv=7.6 m,目的层在垂直方向层厚为2m,仪器参数与模型Ⅰ相同.接收线圈无倾斜时,考察仪器穿过地层倾角α 分别为0°,30°,45°和60°的地层模型时两接收线圈所接收到电磁波信号的幅值比和相位差,如图 10a10b所示.

图 10 模型Ⅲ(a)幅值比和(b)相位差(θ=0°) Fig. 10 Amplitude ratio(a)and phase difference(b)formodel Ⅲ (θ=0°)

当仪器行进至目的层(高阻油层),地层垂直方向深度Lv=3.3 m,地层倾角为45°,利用自适应hp-FEM 算法对模型Ⅲ 进行数值模拟时,目的层部分求解区域的电场分布、网格细化及网格定阶图,如图 11所示.通过图 11可以看出,由于上、下围岩和目的层均为各向异性介质,其水平电导率和垂直电导率各不相同,发射线圈产生的电场分别沿水平和垂直方向产生衰减,在地层分界面处衰减较快并伴有电场不连续.此时,自适应hp-FEM 算法在地层分界面处网格细化程度明显加强,为保证计算精度网格剖分越细计算精度则越高,但是网格过度细化会影响计算时间、增加未知量个数以及增大计算机内存的使用量,因此网格的细化只能达到一定的程度,并非越细越密越好.在不改变网格尺寸的条件下,提高计算区域多项式阶次同样可以达到提高求解精度的要求.在地层分界面处,网格细化程度增强同时求解区域阶次升高的主要原因是由于电场在地层分界面处分布不连续而引起的.当接收线圈倾斜45°时,利用自适应hp-FEM 算法考察仪器穿过地层倾角α 分别为0°,30°,45°和60°的地层模型时,两接收线圈所接收到电磁波信号的幅值比和相位差,如图 12a12b所示.

图 11 电场分布与网格细化及阶次分布图(α = 45°) (a)电场分布;(b)网格细化;(c)阶次分布. Fig. 11 Distribution of electric field, mesh and order (α = 45°) (a) Distribution of electric field; (b) Mesh subdivision: (c) Order distribution.
图 12 模型Ⅲ (a)幅值比和(b)相位差(θ=45°) Fig. 12 Amplitude ratio (a) and phase difference (b) formodel Ⅲ (θ=45°)

通过图 12可以看出,当接收线圈倾斜45°时两接收线圈之间的幅值比和相位差与图 10 相比变化较大,在地层分界面处极化角更加明显.当仪器以不同的倾角穿过地层时,地层的水平分量电导率和垂直分量电导率对仪器测量结果将产生不同程度的影响,仪器测量曲线表现为沿仪器行进方向电导率作用比较明显.设目的层(高阻层)为含油层,通过图 12可以看出,当仪器从低阻层(上围岩)钻入高阻层(含油层)时,幅值比和相位差曲线产生的极化角较浅,无法及时判断仪器在行进过程中钻入地层是否为高阻含油层;当仪器接收线圈倾斜后,仪器测井曲线所产生的极化角较深,可以准确及时地判断出仪器在行进过程中所处的地层位置以及地层性质.根据接收线圈倾斜后随钻电阻率测井仪器测井曲线(幅值比和相位差)中的极化角,可预先指示出高阻层和低阻层边界,从而提供准确的地质导向信息并及时调整测井仪器行进方向,使仪器始终保持在高阻含油层内.此外,极化角的出现不仅可以对随钻电阻率测井过程中仪器导向提供参考,还可以对地质构造解释和地层界面判断提供依据.

为了验证自适应hp-FEM 算法在模型Ⅲ 下的精确性和快速性,利用自适应h-FEM 算法和自适应p-FEM 算法计算了接收线圈无倾角且地层倾角α 为45°时的仪器响应.通过计算可得,自适应hp-FEM 算法所需要的自由度范围从18786 个至21601个,平均使用自由度为20281个;计算时间范围从1.46s至2.25s,平均计算时间约为1.79s;全局误差最大值为0.871%,最小值为0.678%,全局误差平均值约为0.776%.自适应h-FEM 算法平均使用自由度为36475个,平均计算时间约为1.99s,全局误差平均值约为0.87%.自适应p-FEM 算法平均使用自由度为29075 个,平均计算时间约为1.88s,全局误差平均值约为0.91%.因此,自适应hp-FEM 算法在计算时间和计算精度等方面优于其他两种算法,将该方法应用于模型Ⅲ有效可行.

利用自适应hp-FEM 算法,在三种不同类型地层模型下,对随钻电阻率测井仪器响应进行了数值模拟,对不同地层模型下仪器的响应特性进行了比较和评价,以此来验证该方法在随钻电阻率测井仪器响应数值模拟中的有效性.数值模拟结果显示,自适应hp-FEM 算法在三种地层模型中,其计算精度可达0.001%,全局误差可保持在1%以内,并具有较快的计算速度.因此,自适应hp-FEM 算法在随钻电阻率测井仪器响应数值模拟中行之有效.需要指出,以上三种模型中仪器接收线圈倾斜角相同且发射线圈无倾斜,实际工程应用中接收线圈倾角可选择其他角度,但应避免线圈倾斜角与地层倾斜角相似;对于发射线圈倾斜而接收线圈无倾斜的情况可以得到上述类似结论.

4 结 论

本文利用自适应hp-FEM 算法对不同地质环境下随钻电阻率测井仪器响应进行了数值模拟,自适应hp-FEM 算法计算精度高、收敛速度快,能够有效地对随钻电阻率测井过程中仪器响应特性进行评价.

(1) 依据自适应hp-FEM 算法建立模型,对地层倾斜时不同地质条件下随钻电阻率测井仪器的响应情况进行了探讨,同时对自适应hp-FEM 算法应用于不同模型时的精确性和快速性进行了分析,为高精度随钻电阻率测井仪器的研制提供参考.

(2) 基于自适应hp-FEM 算法,利用HERMES有限元软件包和VC+ + 开发平台,可以对随钻电阻率测井仪器响应数值模拟过程中,仪器所处的地层位置,以及由仪器产生的感应电场在地层中和地层分界面处的分布情况进行显示,将其与随钻电阻率测井仪器测量结果相结合可及时准确地提供钻井地质导向信息.

(3) 自适应hp-FEM 算法能够根据地层模型的实际情况自动调整计算策略,具有不必求解整个计算域内的精确解就能使有用信号快速准确地逼近其真值.在随钻电阻率测井仪器响应数值模拟过程中,该算法具有较快的计算速度和较高的计算精度,将其应用于随钻电阻率测井仪器探测特性评价以及仪器参数标定将具有十分重要的意义.此外,本文依据自适应hp-FEM 算法对随钻电阻率测井仪器响应进行数值模拟时未考虑其他因素的影响,实际中有很多因素会影响随钻电阻率测井仪器的测量结果,如:井斜角的影响、地层电阻率对比度的影响、泥浆滤液侵入的影响以及仪器结构和仪器精度的影响,等.因此,下一步将利用自适应hp-FEM 算法重点研究井斜情况下,上述各影响因素对随钻电阻率测井仪器响应及测量结果的影响.

参考文献
[1] 时鹏程. 随钻测井技术在我国石油勘探开发中的应用. 测井技术 , 2002, 26(6): 441–445. Shi P C. M/LWD technology plays an important role in China oilfield development. Well Logging Technology (in Chinese) , 2002, 26(6): 441-445.
[2] 孙向阳, 聂在平, 赵延文, 等. 用矢量有限元方法模拟随钻测井仪在倾斜各向异性地层中的电磁响应. 地球物理学报 , 2008, 51(5): 1600–1606. Sun X Y, Nie Z P, Zhao Y W, et al. The electromagnetic modeling of logging-while-drilling tool in tilted anisotropic formations using vector finite element method. Chinese J. Geophys (in Chinese) , 2008, 51(5): 1600-1606.
[3] Chen Q, Pardo D, Li H B, et al. New post-processing method for interpretation of through casing resistivity (TCR) measurements. J. Appl. Geophys , 2011, 74(1): 19-25. DOI:10.1016/j.jappgeo.2011.02.011
[4] Zhou Q, Hilliker D J. MWD resistivity tool response in a layered medium. Geophysics , 1991, 56(11): 1738-1748. DOI:10.1190/1.1442986
[5] Lovell J R. Finite element methods in resistivity logging [Ph. D. thesis]. The Netherlands: Delft University of Technology, 1993.
[6] Li J, Liu C. A three-dimensional transmission line matrix method (TLM) for simulations of logging tools. IEEE Transactions on Geoscience and Remote Sensing , 2000, 38(4): 1522-1529. DOI:10.1109/36.851952
[7] Zhang Z Q, Liu Q H. Applications of the BCGS-FFT method to 3-D induction well logging problems. IEEE Transactions on Geoscience and Remote Sensing , 2003, 41(5): 998-1004. DOI:10.1109/TGRS.2003.811547
[8] Wang T, Signorelli J. Finite-difference modeling of electromagnetic tool response for logging while drilling. Geophysics , 2004, 69(1): 152-160. DOI:10.1190/1.1649383
[9] 魏宝君, 田坤, 张旭, 等. 用并矢Green 函数的矢量本征函数展开式评价偏心对随钻电磁波电阻率测井响应的影响. 中国石油大学学报 (自然科学版) , 2010, 34(5): 57–62. Wei B J, Tian K, Zhang X, et al. Evaluating influence of eccentricity on response of electromagnetic wave resistivity logging-while-drilling by vector eigenfunction expansion formulae for dyadic Green's functions. Univ. Petro. (Nat. Sci|) (in Chinese) , 2010, 34(5): 57-62.
[10] 刘福平, 高杰, 孙宝佃, 等. 实际井眼条件下过套管电阻率测井响应的传输线方程正演算法. 地球物理学报 , 2007, 50(6): 1905–1913. Liu F P, Gao J, Sun B D, et al. Forward calculation of the resistivity logging response through casing by transmission line equation for multi-layer formations. Chinese J. Geophy (in Chinese) , 2007, 50(6): 1905-1913.
[11] 邢光龙, 王宏建, 杨善德. 二维轴对称介质中电磁波测井的响应函数. 地球物理学报 , 2008, 51(3): 924–932. Xing G L, Wang H J, Yang S D. The response functions of electromagnetic wave logs in the 2-D axis-symmetric formation. Chinese J. Geophys (in Chinese) , 2008, 51(3): 924-932.
[12] Nie X C, Yuan N, Liu C R. Simulation of LWD tool response using a fast integral equation method. IEEE Transactions on Geoscience and Remote Sensing , 2010, 48(1): 72-81. DOI:10.1109/TGRS.2009.2027112
[13] Babuka I, Suri M. The h-p version of the finite element method with quasiuniform meshes. Mathematical Modelling and Numerical Analysis , 1987, 21(2): 199-238.
[14] Demkowicz L, Vardapetyan L. Modeling of electro magnetic absorption/scattering problems using hp-adaptive finite elements. Computer Methods in Applied Mechanics and Engineering , 1998, 152(1-2): 103-124. DOI:10.1016/S0045-7825(97)00184-9
[15] Pardo D, Demkowicz L, Torres-Verdín C, et al. A self- adaptive goal-oriented hp finite element method with electromagnetic applications. Computer Methods in Applied Mechanics and Engineering , 2007, 196: 3585-3597. DOI:10.1016/j.cma.2006.10.016
[16] 陈晓晖, 刘得军, 马中华. 基于高精度自适应hp-FEM的随钻电阻率测井电场数值模拟. 计算物理 , 2011, 28(1): 50–56. Chen X H, Liu D J, Ma Z H. Numerical simulation of electric field in resistivity LWD using high accuracy self- adaptive hp-FEM. Comput. Phys (in Chinese) , 2011, 28(1): 50-56.
[17] Tomá V, Pavel S, Martin Z, et al. Modular hp-FEM system HERMES and its application to Maxwell's equations. Mathematics and Computers in Simulation , 2007, 76(1-3): 223-228. DOI:10.1016/j.matcom.2007.02.001
[18] Demkowicz L, Rachowicz W, Devloo P. A fully automatic hp-adaptivity. J. Sci. Comput , 2002, 17(1-4): 117-142.
[19] olín P, erveny J, Doleel I. Arbitrary-level hanging nodes and automatic adaptivity in the hp-FEM. Mathematics and Computers in Simulation , 2008, 77(1): 117-132. DOI:10.1016/j.matcom.2007.02.011
[20] olín P, Dubcova L, Kruis J. Adaptive hp-FEM with dynamical meshes for transient heat and moisture transfer problems. J. Comput. Appl. Math , 2010, 233(12): 3103-3112. DOI:10.1016/j.cam.2009.07.025
[21] 杨锦舟, 魏宝君, 林楠. 倾斜线圈随钻电磁波电阻率测量仪器基本原理及其在地质导向中的应用. 中国石油大学学报 (自然科学版) , 2009, 33(1): 44–49. Yang J Z, Wei B J, Lin N. Basic theory of electromagnetic wave resistivity measurement while drilling tool with tilted antennas and its application for geo-steering. Chin. Univ. Petrol. (Nat. Sci|) (in Chinese) , 2009, 33(1): 44-49.