張盼, 王正楷, 鄧少貴

(中國石油大學地球科學與技術(shù)學院, 山東 青島 266580)

0 引 言

套管井電阻率測井響應(yīng)模擬的主要困難在于金屬套管和其他介質(zhì)(水泥環(huán)、地層等)的強電導率對比度所導致的計算精度問題[1-2]。傳輸線方程較好描述了套管井環(huán)境場的分布特點[3],成為過套管電阻率測井理論計算的重要方程[4]。劉福平等[5]建立了改進型傳輸線模型,把影響因素吸收到傳輸線方程的系數(shù)中,假定介質(zhì)參數(shù)分塊變化,提出了改進型傳輸線方程的遞推算法,在對各種因素的考察中獲得廣泛的應(yīng)用[6-8]。但是,當套管縱向單位長度電阻、水泥環(huán)單位厚度電阻或地層電阻率非分塊均勻時,改進型傳輸線的遞推算法變得困難。為提升改進型傳輸線的適用性和計算精度,采用自適應(yīng)hp有限元求解具有靈活性。

本文在自適應(yīng)hp有限元框架下[9-12],求解ECOS過套管電阻率的改進型傳輸線方程,自動調(diào)節(jié)網(wǎng)格剖分大小h和單元插值函數(shù)階數(shù)p,尋找適合問題解本身特征的擬優(yōu)化的網(wǎng)格。在相同自由度條件下,自適應(yīng)hp有限元能夠達到最高的計算精度。實際計算中,基于先驗信息的擬最優(yōu)網(wǎng)格的幾乎是不可行的。自適應(yīng)hp有限元從粗網(wǎng)格迭代開始,不斷尋找當前粗網(wǎng)格的細化網(wǎng)格,有限元誤差可控。

1 ECOS過套管電阻率的測井原理

ECOS的過套管電阻率測井儀器結(jié)構(gòu)如圖1所示,以測量電極N為中心上下對稱分布1對測量電極M1、M2和發(fā)射電極A1和A2,測量電極N到測量電極M1的距離為0.5 m,到發(fā)射電極A1的距離為1.8 m。測井時,儀器推靠到井壁固定,上發(fā)射電極A1發(fā)射電流強度為I1的電流,分別測量測量電極N的電壓UN(I1)、測量電極M1與N的電壓和測量電極M2與N的電壓;下發(fā)射電極A2發(fā)射電流強度為I2的電流,分別測量測量電極N的電壓UN(I2),測量電極M1與N的電壓和測量電極M2與N的電壓,根據(jù)電壓的線性疊加特性,得到過套管的電阻率表達式

(1)

式中,UN(I1)、UN(I2)分別為上供電電極和下供電電極供電時,測量電極N同套管的接觸點的電位;ΔUN(I1)、ΔUN(I2)分別為上供電電極和下供電電極供電時,電極系的2個電極M1、M2的一階差分;Δ2UN(I1)、Δ2UN(I2)分別為上供電電極和下供電電極供電時,3個測量電極同套管接觸點之間套管段上電位的二階差分;I1、I2分別為A1和A2供電電極向套管注入的電流。

2 一維自適應(yīng)hp有限元方法

2.1 過套管電阻率的傳輸線方程

設(shè)單極供電電極位于井軸線上z點,則產(chǎn)生的電位沿井軸的分布可表述為

(2)

當Rr/Rz?1時,套管邊界面上的電流近似垂直金屬套管壁流向地層[13],徑向電流可以表示為

(3)

式中,U(z)為井軸z處的電位,V;Jz(z)為經(jīng)過同一坐標的套管橫截面的電流,A;Jr(z)為由單位長度深度段井壁進入周圍地層的電流,A/m;Rr為單位厚度介質(zhì)對電流Jr(z)的電阻,Ω/m;Rz為單位長度套管對軸向電流Jz(z)的電阻,Ω/m。

以發(fā)射電極中心圓柱的中心作一個套管圓柱體,應(yīng)用恒定電流的連續(xù)性方程,即

∮sJ·ds=I

(4)

將式(4)整理得

(5)

式中,I表示在z0點施加的直流電流的強度,A。

考慮邊界條件,套管壁上的電勢滿足下列邊值問題

u|Γ=0

(6)

根據(jù)Galerkin方法,式(6)等價成下列的問題

(7)

2.2 變階插值函數(shù)

自適應(yīng)hp有限元要求求解區(qū)域單元上的形狀函數(shù)為層狀形狀函數(shù),即高階的形狀函數(shù)包含低階的形狀函數(shù),單元上的形狀函數(shù)通過拼合而構(gòu)成有限元空間的基函數(shù)[14]。通常情況下,選擇Lobatto基函數(shù)作為自適應(yīng)hp有限元的基函數(shù)。它是由2個一階函數(shù)和多個高階函數(shù)組成,k階Lobatto基函數(shù)定義為

(8)

(9)

(10)

式中,Lk(x)為k次勒讓德多項式,其定義為

(11)

2.3 投影型插值方法

投影型插值是單元頂點上的插值和單元內(nèi)部投影相結(jié)合的擬合方法,它是自適應(yīng)hp有限元自適應(yīng)策略的核心組成部分[15]。投影型插值函數(shù)定義:假設(shè)某個單元上連續(xù)函數(shù)u,定義在[a,b]上到的p階層次形狀函數(shù)為φ1,φ2,φ3,i(i=1,2,…,p-1),則投影型插值函數(shù)為

(12)

式中,u1=u(a),u2=u(b),u3,j由在單元內(nèi)部的投影系數(shù)求得下列的線性方程得到

(13)

投影型插值保證插值算子是局部的,它只與單個單元上的解的信息有關(guān);是整體連續(xù)的,每個單元的算子在整個網(wǎng)格連續(xù);是最優(yōu)的,它具有整體插值的特性,對h和p都是漸近收斂的。這些特點使得自適應(yīng)hp策略效率更高。

2.4 自適應(yīng)hp策略

自適應(yīng)hp迭代是從一個粗網(wǎng)格開始,利用參考解在粗網(wǎng)格單元多種可能細化網(wǎng)格下的投影型插值的誤差相對大小,尋找當前粗網(wǎng)格相對應(yīng)的擬最優(yōu)細化網(wǎng)格。當前粗網(wǎng)格的擬最優(yōu)細化網(wǎng)格滿足

hpopt=

(14)

實際上,真解u是未知的,采用參考解來近似。假設(shè)給定的當前粗網(wǎng)格單元尺寸為h和網(wǎng)格單元上的插值函數(shù)的階數(shù)p,將每個單元二分變成2個等尺寸單元,變成h/2,每個單元上的階數(shù)提高1階,在新的全局一致細化網(wǎng)格上求得的解叫做當前粗網(wǎng)格下的參考解,記為uh/2,p+1。利用參考解和投影型插值,式(14)的全局優(yōu)化問題可以通過粗網(wǎng)格每個單元上局部的優(yōu)化問題來求解,將每個單元上增加的單位個自由度誤差下降最大值稱為誤差下降速率,其表達式為

(15)

當前粗網(wǎng)格中每個單元的最優(yōu)化的細化單元采用枚舉的方式獲得,枚舉的細化單元的自由度個數(shù)介于當前粗網(wǎng)格單元和當前粗網(wǎng)格單元的參考單元自由度個數(shù)之間。

2.5 自適應(yīng)hp算法

給定初始的粗網(wǎng)格單元和每個單元上插值函數(shù)的階數(shù),一般取1～2階。

步驟1:計算當前粗網(wǎng)格上的有限元的解,記為uhp。

步驟2:對粗網(wǎng)格上的每個單元,二分單元h并將得到每個單元的多項式階數(shù)p提升到p+1,計算有限元的參考解uh/2,p+1,如果相對能量誤差‖uh/2,p+1-uh,p‖/‖uh/2,p+1‖小于容忍誤差,停止迭代,并輸出hp有限元的解uh/2,p+1,否則轉(zhuǎn)入步驟3。

步驟3:將參考解投影到粗網(wǎng)格中的每個單元所有可能細化單元上,利用式(15)計算誤差下降速率ri,確定每個單元誤差下降速率的h細化和p細化。

3 測井響應(yīng)計算

3.1 方法驗證

圖3 問題真解和粗網(wǎng)格有限元的解

圖4 問題真解和自適應(yīng)hp有限元的解

為了驗證自適應(yīng)hp有限元方法的有效性和程序的正確性,給出一個具有解析解的邊值問題,該解析解在中心部分劇烈變化。圖3為問題真解和粗網(wǎng)格有限元的解,圖3中的橫坐標表示求解區(qū)域,縱坐標表示解的大小,實線為實際問題的解,虛線為粗網(wǎng)格剖分下的有限元的解。圖3中低端部分表示網(wǎng)格的大小,顏色表示階數(shù)。圖3中右側(cè)部分是顏色刻度,不同的階數(shù)對應(yīng)不同的顏色。可以看出,求解該邊值問題的初始網(wǎng)格為4個單元,每個單元的插值函數(shù)都為3階,粗網(wǎng)格的有限元解與真解之間的誤差較大。圖4為問題的真解和自適應(yīng)hp有限元解,計算中相對能量誤差取為0.1%,從圖4中可以看出二者吻合相當好。細化網(wǎng)格集中在問題真解劇烈變化,說明基于參考解投影型插值的自適應(yīng)hp策略捕捉到問題真解的劇烈變化的特性,充分體現(xiàn)了自適應(yīng)hp方法的正確性和有效性(見圖4)。

3.2 金屬套管異常條件影響

為研究金屬套管電導率異常對測井響應(yīng)影響,所建研究模型的金屬套管的電導率為5×106S/m,深度分別為1 993.5～1 995.5 m、1 999.5～2 000.5 m和2 004.5～2 006.5 m的電導率為3×106S/m。圖5為在該模型下的過套管電阻率測井響應(yīng),圖5中藍線為地層的真電阻率,帶點的實線為過套管電阻率的測井響應(yīng)。從圖5中可以看出,在金屬套管電導率發(fā)生突變的地方,視電阻率出現(xiàn)類似一個周期正弦的異常,該周期的寬度約為電極系中M1和M2之間的距離(1 m),金屬套管高電導率部分對應(yīng)約為半個周期的低值部分,而金屬套管低電導率部分對應(yīng)約為半個周期的高值部分,如果金屬電導率分段異常部分小于M1和M2之間的距離,過套管電阻率異常將出現(xiàn)重疊。

圖5 金屬套管異常條件下過套管電阻率測井響應(yīng)

假設(shè)金屬套管的電導率為5×106S/m不變,在深度為1 982～1 986 m、1 990～1 994 m和1 998～2 002 m的電導率分別為3×106、2×106、4×106S/m。圖6為該金屬套管電導率條件下的過套管電阻測井響應(yīng),圖6中,除了圖5中展示的特點,類似正弦的異常值的大小與金屬套管電導率的突變值有關(guān),金屬套管的突變值越大,過套管電導率的異常值越大。

圖6 金屬套管異常條件下過套管電導率

圖7為金屬套管電導率線性變化模型,圖8為對應(yīng)的過套管電阻率測井響應(yīng),結(jié)合這2張圖可以看出,金屬套管線性變化部分對應(yīng)的過套管電阻率沒有出現(xiàn)異常,異常值出現(xiàn)在金屬套管電導率導數(shù)突變的地方,異常值隨線性變化率增大而增大。一般情況下,異常值通常很小。

圖7 金屬套管電導率線性變化模型

圖8 金屬套管電導率線性變化過套管電阻率測井響應(yīng)

對于金屬套管電導率正弦變化的情形,視電阻率也出現(xiàn)異常,異常的低值與套管的電導率的高值對應(yīng),異常的高值與套管電導率的低值對應(yīng),但電導率的異常值一般很小。

3.3 圍巖電阻率影響

圖9為分層地層的過套管電阻率測井響應(yīng)。該模型中背景電阻率為2 Ω·m,深度1 994～1 995 m、1 999～2 001 m、2 005～2 008 m之間為3層厚度不同電阻率為10 Ω·m的地層。從圖9中可以看出,ECOS過套管電阻率縱向分辨率的極限為1 m,由于計算中記錄點長度為0.1 m,在層厚為1 m的地層中最多只有一個點反映到地層的真實電阻率,其余各點都受到圍巖的影響。當記錄點長度變大時,測量到圍巖的影響就愈加明顯。記錄點長度小于1 m時,層厚大于2 m的地層在實際測井中受圍巖的影響可以不計。圖10為地層電阻率為正弦電阻率變化,視電阻率大致能夠反映原來地層的形態(tài),在高值點由于受到上下低阻圍巖的影響,視電阻率變小;在低值點由于受到上下高阻圍巖的影響,視電阻率變大。

圖9 分層地層的過套管電阻率測井響應(yīng)

圖10 地層電阻率連續(xù)變化時測井響應(yīng)

圖11 水泥環(huán)存在時過套管電阻率測井響應(yīng)

3.4 水泥環(huán)的影響

圖11為水泥環(huán)存在時過套管電阻率測井響應(yīng)。計算模型:套管的電導率1×105S/m,水泥環(huán)的厚度為5 cm,電阻率為5 Ω·m,共有10層地層(見圖11)。圖11中,水泥環(huán)電阻率大于地層電阻率,水泥環(huán)存在使視電阻率增加;水泥環(huán)電阻率等于地層電阻率,水泥環(huán)存在不影響電阻率;水泥環(huán)電阻率小于地層電阻率,水泥環(huán)的電阻率存在使視電阻率降低,且二者差異越大,影響越明顯。

4 結(jié) 論

(1) 采用一維自適應(yīng)hp有限元求解了ECOS過套管電阻率在不同條件下的測井響應(yīng),數(shù)值計算結(jié)果表明,該方法適應(yīng)性強和計算精度高,適合于介質(zhì)參數(shù)連續(xù)變化等復雜條件下問題的求解。

(2) 金屬套管電導率的連續(xù)變化對ECOS過套管電阻率視電阻率的影響與間斷變化的影響特征不同,通常條件下,連續(xù)變化的影響很小,可以忽略。

(3) 水泥環(huán)電阻率連續(xù)變化和地層電阻率連續(xù)變化對ECOS過套管電阻率的影響與相應(yīng)參數(shù)間斷變化特征類似。

參考文獻:

[1] SINGER B H, FANINI O. Through-casing Resistivity: 2-D and 3-D Distortions and Correction Techniques [C]∥ SPWLA 36th Annual Logging Symposium, June 26-29, 1995.

[2] 劉福平, 陳小安, 楊長春. 過套管電阻率測井響應(yīng)數(shù)值模擬研究進展 [J]. 地球物理學進展, 2011, 26(6): 2018-2025.

[3] KAUFMAN A A. The Electrical Field in a Borehole With a Casing [J]. Geophysics, 1990, 55(1): 29-38.

[4] KAUFMAN A A. WIGHTMAN W E. A Transmission-line Modelfor Electrical Logging Through Casing [J]. Geophysics, 1993, 58(12): 1739-1747.

[5] 劉福平, 高杰, 包德洲, 等. 實際井眼條件下過套管電阻率測井響應(yīng)的傳輸線方程的正演算法 [J]. 地球物理學報, 2007, 50(6): 1905-1913.

[6] 高杰, 劉福平. 過套管電阻率測井方法研究 [C]∥ 第四屆中俄測井國際學術(shù)交流會 [M].北京: 石油工業(yè)出版社, 2006: 339-349.

[7] 高杰, 劉福平, 保德州, 等. 非均勻套管井中的過套管電阻率測井響應(yīng) [J]. 地球物理學報, 2008, 51(4): 1255-1261.

[8] 魏寶君, 田坤, 張旭, 等. 模擬過套管電阻率測井響應(yīng)的遞推矩陣方法 [J]. 中國石油大學學報: 自然科學版, 2011, 35(6): 59-65.

[9] PARDO D, DEMKOWICZ L F, TORRES-VERDIN, et al. A Goal-oriented Hp-adaptive Finite Element Method with Electromagnetic Applications. Part I: Electrostatics, Int. J. Numer [J]. Methods Engrg. 2006, 65(8): 1269-1309.

[10] VARCIACASTILLO L E, PARDO D, GARCA-CASTILLO L E, et al. A Two-dimensional Self-adaptivehpFinite Element Method for the Characterization of Waveguide Discontinuities. Part I: Energy-norm Based Automatic Hp-adaptivity [J]. Comput. Methods Appl. Mech. Engrg., 2007, 196: 4823-4852.

[11] PARDO D, DEMKOWICZ L F, TORRES-VERES C, et al. Simulation of Resistivity Logging-while-drilling (LWD) Measurements Using a Self-adaptive Goal-orientedhp-finite Element Method [J]. SIAM J. Appl. Math, 2006, 66(6): 2085-2106.

[12] KURTZ J, DEMKOWICZ L F. A Fully Automatichp-adaptivity for Elliptic PDEs in Three Dimensions [J]. Computer Methods in Applied Mechanics and Engineering, 2007(196): 3534-3545.

[13] 劉福平, 陳小安, 趙寶成, 等. 過套管電阻率測井高電阻率層段不穩(wěn)定性分析 [J]. 測井技術(shù), 2013, 37(01): 85-89.

[14] 李輝, 劉得軍, 劉悅, 等. 基于自適應(yīng)hp-FEM的過套管電阻率測井儀器響應(yīng)數(shù)值模擬 [J]. 地球物理學進展, 2013, 28(06): 3243-3253.

[15] DEMKOWICZ L F. Computing withhp-Adaptive Finite Elements, Vol. I. One and Two Dimensional Elliptic and Maxwell Problems [M]. Taylor and Francis: Chapman & Hall/CRC Press, 2007.