張盼, 王正楷, 鄧少貴

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

0 引 言

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

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

1 ECOS過(guò)套管電阻率的測(cè)井原理

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

(1)

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

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

2.1 過(guò)套管電阻率的傳輸線方程

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

(2)

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

(3)

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

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

∮sJ·ds=I

(4)

將式(4)整理得

(5)

式中,I表示在z0點(diǎn)施加的直流電流的強(qiáng)度,A。

考慮邊界條件,套管壁上的電勢(shì)滿足下列邊值問(wèn)題

u|Γ=0

(6)

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

(7)

2.2 變階插值函數(shù)

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

(8)

(9)

(10)

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

(11)

2.3 投影型插值方法

投影型插值是單元頂點(diǎn)上的插值和單元內(nèi)部投影相結(jié)合的擬合方法,它是自適應(yīng)hp有限元自適應(yīng)策略的核心組成部分[15]。投影型插值函數(shù)定義:假設(shè)某個(gè)單元上連續(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)

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

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

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

hpopt=

(14)

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

(15)

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

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

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

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

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

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

3 測(cè)井響應(yīng)計(jì)算

3.1 方法驗(yàn)證

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

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

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

3.2 金屬套管異常條件影響

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

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

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

圖6 金屬套管異常條件下過(guò)套管電導(dǎo)率

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

圖7 金屬套管電導(dǎo)率線性變化模型

圖8 金屬套管電導(dǎo)率線性變化過(guò)套管電阻率測(cè)井響應(yīng)

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

3.3 圍巖電阻率影響

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

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

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

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

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

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

4 結(jié) 論

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

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

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

參考文獻(xiàn):

[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] 劉福平, 陳小安, 楊長(zhǎng)春. 過(guò)套管電阻率測(cè)井響應(yīng)數(shù)值模擬研究進(jìn)展 [J]. 地球物理學(xué)進(jìn)展, 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] 劉福平, 高杰, 包德洲, 等. 實(shí)際井眼條件下過(guò)套管電阻率測(cè)井響應(yīng)的傳輸線方程的正演算法 [J]. 地球物理學(xué)報(bào), 2007, 50(6): 1905-1913.

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

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

[8] 魏寶君, 田坤, 張旭, 等. 模擬過(guò)套管電阻率測(cè)井響應(yīng)的遞推矩陣方法 [J]. 中國(guó)石油大學(xué)學(xué)報(bào): 自然科學(xué)版, 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] 劉福平, 陳小安, 趙寶成, 等. 過(guò)套管電阻率測(cè)井高電阻率層段不穩(wěn)定性分析 [J]. 測(cè)井技術(shù), 2013, 37(01): 85-89.

[14] 李輝, 劉得軍, 劉悅, 等. 基于自適應(yīng)hp-FEM的過(guò)套管電阻率測(cè)井儀器響應(yīng)數(shù)值模擬 [J]. 地球物理學(xué)進(jìn)展, 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.