三維感應(yīng)測(cè)井儀器刻度響應(yīng)有限元法數(shù)值模擬研究

俞燕明，肖加奇，柯式鎮(zhèn)，李強(qiáng)

(1.中國(guó)石油長(zhǎng)城鉆探工程有限公司測(cè)井技術(shù)研究院，北京 100124;2.中國(guó)石油大學(xué)，北京 100124)

0 引言

三維感應(yīng)測(cè)井目前商業(yè)化的儀器主要是斯倫貝謝公司的 Rt-Scanner[1]和阿特拉斯公司的 3DEX[1]。三維感應(yīng)測(cè)井的數(shù)值模擬研究已經(jīng)開(kāi)展多年，主要方法有均勻場(chǎng)解析法[1-2]、傳輸線(xiàn)理論法[3](TLM)、有限差分方法[4-6](FDM)、有限元法[7-9](FEM)、數(shù)值模式匹配法[10](NMM)、混合勢(shì)理論法[11]和有限體積法[12](FVM)等。國(guó)內(nèi)外的數(shù)值模擬研究主要針對(duì)地層模型開(kāi)展算法研究，對(duì)于刻度響應(yīng)及其影響的計(jì)算未見(jiàn)到報(bào)道，相關(guān)文獻(xiàn)主要是針對(duì)陣列感應(yīng)測(cè)井儀器采用解析解的方法進(jìn)行刻度環(huán)參數(shù)及刻度響應(yīng)的計(jì)算[13-14]。本文針對(duì)三維感應(yīng)測(cè)井儀器刻度響應(yīng)進(jìn)行數(shù)值模擬研究，建立了三維有限元數(shù)值模擬計(jì)算程序，并在單一傾斜刻度環(huán)下的刻度響應(yīng)計(jì)算與物理實(shí)驗(yàn)結(jié)果進(jìn)行了驗(yàn)證，計(jì)算考察了刻度環(huán)半徑、刻度環(huán)電阻、刻度環(huán)傾斜角度、刻度環(huán)傾斜方位、刻度環(huán)水平放置位置對(duì)儀器響應(yīng)的影響關(guān)系，為三維感應(yīng)測(cè)井儀器刻度裝置設(shè)計(jì)提供理論指導(dǎo)。

1 三維感應(yīng)測(cè)井原理

與常規(guī)感應(yīng)測(cè)井方法不同，三維感應(yīng)測(cè)井采用3個(gè)相互正交的發(fā)射線(xiàn)圈(Tx、Ty和Tz)和3個(gè)相互正交的接收線(xiàn)圈(Rx、Ry和Rz)進(jìn)行發(fā)射和接收(見(jiàn)圖1)，以便在各個(gè)深度點(diǎn)上同時(shí)測(cè)量9個(gè)磁場(chǎng)分量(Hxx、Hxy、Hxz、Hyx、Hyy、Hyz、Hzx、Hzy和 Hzz)，再根據(jù)這9個(gè)分量計(jì)算由9個(gè)電導(dǎo)率分量組成的電導(dǎo)率張量、各向異性指標(biāo)參數(shù)及裂縫參數(shù)等。

圖1 三維感應(yīng)測(cè)井原理示意圖

對(duì)于電性各向異性地層，電導(dǎo)率張量的形式為

對(duì)于砂泥巖薄互層，可以等效為軸向各向異性地層，電導(dǎo)率張量簡(jiǎn)化為

式中，σxx=σyy。

2 有限元計(jì)算方法

對(duì)于感應(yīng)測(cè)井，其電磁場(chǎng)問(wèn)題可以由Maxwell方程組描述

式中，E為電場(chǎng)強(qiáng)度;B為磁場(chǎng)強(qiáng)度;D為電通密度;H為磁通密度;J為導(dǎo)電介質(zhì)體內(nèi)的電流密度;Js為電流源提供的電流密度;ρ為電荷密度。狀態(tài)方程

式中，ε為電容率;σ為電導(dǎo)率;μ為磁導(dǎo)率。

令

設(shè)整個(gè)求解空間被離散化成e0個(gè)單元N0個(gè)節(jié)點(diǎn)，對(duì)于任一節(jié)點(diǎn)i，在直角坐標(biāo)系中，其矢勢(shì)Ai的3個(gè)分量分別記為Axi、Ayi和Azi。每個(gè)單元有n0個(gè)節(jié)點(diǎn)，且相應(yīng)于節(jié)點(diǎn)i的形狀函數(shù)記為，則每個(gè)單元中任意點(diǎn)的磁矢勢(shì)可表示為

根據(jù)泛函的極值求解原理可得

式中，對(duì)單元內(nèi)的泛函偏導(dǎo)數(shù)的計(jì)算為

這樣，對(duì)于某個(gè)單元可以計(jì)算得到單元?jiǎng)偠染仃嚍?/p>

式中，

對(duì)應(yīng)于單元上i節(jié)點(diǎn)的矩陣方程右邊的常數(shù)項(xiàng)按式(20)計(jì)算

將整個(gè)求解空間的所有單元的所有節(jié)點(diǎn)裝配完后形成了有限元求解問(wèn)題的總剛度矩陣[K]及其矩陣方程

對(duì)于超大規(guī)模矩陣方程(21)的求解可以有多種算法。如直接解法中有LU分解法、LDLT分解法、GAUSS法和帶狀矩陣算法等。迭代法分有共軛梯度法和復(fù)雙共軛梯度法等。經(jīng)研究發(fā)現(xiàn)當(dāng)線(xiàn)性方程組系數(shù)矩陣的條件數(shù)較差時(shí)LU分解方法最為穩(wěn)定。當(dāng)系數(shù)矩陣條件數(shù)中等或中等以上時(shí)，采用預(yù)條件的共軛梯度法計(jì)算速度較快，可以用于三維感應(yīng)測(cè)井響應(yīng)的數(shù)值計(jì)算中。

3 計(jì)算實(shí)例與分析

根據(jù)上述算法在ANSYS平臺(tái)上實(shí)現(xiàn)了三維感應(yīng)測(cè)井儀器響應(yīng)的數(shù)值模擬計(jì)算，其過(guò)程是，①選擇電磁功能模塊。②進(jìn)入預(yù)處理階段。根據(jù)儀器參數(shù)建立有限元模型(見(jiàn)圖2)，其中圓柱塊為三維感應(yīng)去直耦線(xiàn)圈單元和接收線(xiàn)圈單元。最右側(cè)的6個(gè)線(xiàn)圈組成了三維感應(yīng)發(fā)射線(xiàn)圈單元。中間大的傾斜環(huán)為刻度環(huán)。其他的小線(xiàn)圈為普通感應(yīng)接收線(xiàn)圈。整個(gè)線(xiàn)圈系中有2套由發(fā)射線(xiàn)圈T、去直耦線(xiàn)圈B和接收線(xiàn)圈R三線(xiàn)圈組成三維感應(yīng)測(cè)量線(xiàn)圈系即淺探測(cè)和深探測(cè)線(xiàn)圈系，其中淺探測(cè)線(xiàn)圈系的發(fā)射到去直耦線(xiàn)圈的距離TB=0.5453 m，發(fā)射到接收線(xiàn)圈的距離TR=0.7620 m，深探測(cè)的TB=1.2700 m，TR=1.6002 m，發(fā)射電流和頻率均為0.45 A和20 kHz，線(xiàn)圈匝數(shù)見(jiàn)表1。除了儀器模型外，還要給出地層模型，通常是圓柱體。本文因?yàn)橹饕怯?jì)算刻度環(huán)響應(yīng)，因而介質(zhì)模型中是刻度環(huán)，其半徑可變，小環(huán)半徑0.25 m，大環(huán)半徑0.55 m，本文采用小環(huán)。刻度環(huán)之外的其他空間是半徑10 m，高度20 m的圓柱體空氣介質(zhì)。③對(duì)所建的計(jì)算模型進(jìn)行網(wǎng)格剖分和屬性賦值。網(wǎng)格采用變尺度自動(dòng)劃分，由內(nèi)到外依次變大?？諝鈱傩允请娮杪薀o(wú)窮大，線(xiàn)圈屬性是電阻率無(wú)窮小，不考慮渦流?？潭拳h(huán)當(dāng)導(dǎo)電介質(zhì)，其電阻率由刻度電阻換算成電阻率賦給。④邊界條件約束和激勵(lì)加載。邊界條件約束主要是無(wú)窮遠(yuǎn)邊界施加矢量磁位A=0。激勵(lì)加載則是將電流密度加到發(fā)射線(xiàn)圈中。1次計(jì)算只加載1個(gè)方向的電流，分3次計(jì)算才能完成所有9個(gè)分量的計(jì)算。⑤方程組求解。給定分析類(lèi)型和加載的信號(hào)頻率，選擇求解器，然后求解。⑥后處理。根據(jù)計(jì)算結(jié)果輸出的參數(shù)，計(jì)算儀器測(cè)量信號(hào)，完成整個(gè)數(shù)值模擬計(jì)算過(guò)程。本文運(yùn)用該方法開(kāi)展了三維感應(yīng)儀器刻度響應(yīng)的計(jì)算、驗(yàn)證及影響因素考察。

圖2 三維感應(yīng)儀器刻度原理示意圖

表1 線(xiàn)圈系中的各線(xiàn)圈匝數(shù)

3.1 刻度環(huán)水平放置位置的影響計(jì)算與驗(yàn)證

不改變刻度環(huán)的其他參數(shù)值，只改變其水平放置位置Zcal，計(jì)算出的儀器響應(yīng)隨刻度環(huán)水平放置位置的變化關(guān)系如圖3和圖4所示。由結(jié)果可以看出，對(duì)于不同探測(cè)深度的線(xiàn)圈系，其刻度響應(yīng)取得極大值的位置不同。對(duì)于2個(gè)ZZ分量，進(jìn)行了物理試驗(yàn)，實(shí)測(cè)的結(jié)果如圖3中的實(shí)測(cè)Vr1zz和實(shí)測(cè)Vr2zz所示。Vr1zz、Vr2zz為淺、深三維感應(yīng)儀器響應(yīng)水平分量;比較圖3的2組曲線(xiàn)可以看出，雖然它們?cè)诹恐瞪洗嬖诓町?主要由于所用的頻率和電流密度值存在差異引起)，但曲線(xiàn)的形狀非常一致。

圖3 三維感應(yīng)儀器刻度響應(yīng)水平分量隨刻度環(huán)水平位置的變化及與實(shí)測(cè)結(jié)果對(duì)比

圖4 三維感應(yīng)儀器刻度響應(yīng)非水平分量隨刻度環(huán)水平位置的變化

3.2 刻度環(huán)半徑的影響

刻度環(huán)傾斜角度、水平放置位置和刻度環(huán)電阻不變，只改變刻度環(huán)半徑RADcal并計(jì)算儀器響應(yīng)，得到淺、深三維感應(yīng)測(cè)井儀器響應(yīng)的ZZ分量隨刻度環(huán)半徑的關(guān)系如圖5所示。由結(jié)果表明，一個(gè)刻度環(huán)不管其半徑如何取，均無(wú)法使得2個(gè)及以上不同探測(cè)深度的線(xiàn)圈系同時(shí)獲得最大響應(yīng)，只能取折衷的半徑值。

圖5 三維感應(yīng)儀器刻度響應(yīng)隨刻度環(huán)半徑的變化

3.3 刻度環(huán)電阻的影響

刻度環(huán)傾斜角度、水平放置、半徑等參數(shù)不變，只改變刻度環(huán)電阻Rcal值，計(jì)算出的儀器響應(yīng)隨刻度環(huán)電阻的變化關(guān)系如圖6和圖7所示。圖6、圖7中各分量均取絕對(duì)值進(jìn)行作圖，可以看出，存在使儀器各分量響應(yīng)共同達(dá)到最大值的刻度電阻Rcal值。

圖6 三維感應(yīng)儀器刻度響應(yīng)水平分量隨刻度電阻的變化

圖7 三維感應(yīng)儀器刻度響應(yīng)垂直和交叉分量隨刻度電阻的變化

3.4 刻度環(huán)傾斜角度的影響

其他刻度環(huán)參數(shù)不變，只改變刻度環(huán)的傾斜角度ALFA，計(jì)算出的儀器響應(yīng)隨ALFA的變化關(guān)系如圖8所示。由圖8可以看出，存在一個(gè)使XZ和ZX分量取得極大(小)值的ALFA。

圖8 三維感應(yīng)儀器刻度響應(yīng)隨刻度環(huán)傾斜角度的變化

3.5 刻度環(huán)傾斜方位的影響

刻度環(huán)的其他參數(shù)不變，只改變刻度環(huán)的傾斜方位XETA，計(jì)算出的儀器響應(yīng)隨XETA的變化關(guān)系如圖9所示。由圖9可以看出，XX、YY、YZ和ZY分量隨方位角呈現(xiàn)余弦變化，XY和YX分量隨方位角呈現(xiàn)正弦變化，ZX和XZ則呈現(xiàn)為半周期的余弦變化。所有分量均存在使其取得極大值的XETA值。

圖9 三維感應(yīng)儀器刻度響應(yīng)隨刻度環(huán)傾斜方位的變化

4 結(jié)論

(1)對(duì)于淺、深探測(cè)的三維感應(yīng)線(xiàn)圈系，其刻度響應(yīng)達(dá)到最大的刻度環(huán)的半徑是不同的，因而無(wú)法確定使得這2個(gè)不同探測(cè)深度的線(xiàn)圈系同時(shí)達(dá)到最大響應(yīng)的刻度環(huán)半徑值，只能折中考慮。為此，可以考慮采用雙環(huán)設(shè)計(jì)。

(2)對(duì)于不同分量取得刻度響應(yīng)的極大值其傾斜角和方位角是不同的，因而要取得最大刻度響應(yīng)需要針對(duì)不同的分量改變不同的刻度環(huán)傾斜角和傾斜方位來(lái)配合。

(3)對(duì)于給定半徑和位置的刻度環(huán)，存在一個(gè)使其刻度響應(yīng)達(dá)到極值的刻度電阻值。

(4)對(duì)于深、淺線(xiàn)圈系，儀器刻度響應(yīng)取得極值的位置不同，刻度時(shí)應(yīng)予以考慮。

[1]張國(guó)艷，肖加奇，郝永杰.三維感應(yīng)測(cè)井?dāng)?shù)值計(jì)算與理論分析[J].測(cè)井技術(shù)，2012，36(1):15-19.

[2]白彥，仵杰.三分量感應(yīng)線(xiàn)圈系在均勻地層中的響應(yīng)特性分析[J].電子測(cè)試，2010，215(12):22-25.

[3]Wang H N，Tao H G，Yao J J，et al.Fast Multi-parameter Reconstruction of Multi-component Induction Well Logging Datum in Deviated Well in a Horizontally Stratified Anisotropic Formation[C]∥ IEEE Trans on Geosci and Remote Sensing，2008，46(5):1525-1534.

[4]汪功禮，張庚驥，崔鋒修，等.三維感應(yīng)測(cè)井響應(yīng)計(jì)算的交錯(cuò)網(wǎng)格有限差分法[J].地球物理學(xué)報(bào)，2003，46(4):561-567.

[5]王昌學(xué)，楊韋華，儲(chǔ)昭坦，等.多分量感應(yīng)測(cè)井響應(yīng)的交錯(cuò)網(wǎng)格有限差分法模擬[J].石油大學(xué)學(xué)報(bào):自然科學(xué)版，2005，29(3):35-40.

[6]沈金松.用有限差分法計(jì)算各向異性介質(zhì)中多分量感應(yīng)測(cè)井的響應(yīng)[J].地球物理學(xué)進(jìn)展，2004，19(1):101-107.

[7]李飛虎，張中慶，王卓遠(yuǎn).用矢量棱邊元素法模擬三維感應(yīng)測(cè)井響應(yīng)[J].復(fù)旦學(xué)報(bào):自然科學(xué)版，2009，48(5):560-566.

[8]黃臨平，戴世坤.復(fù)雜條件下3D電磁場(chǎng)有限元計(jì)算方法[J].地球科學(xué)，2002，27(6):775-779.

[9]張繼鋒，湯井田，王燁，等.被動(dòng)源電磁測(cè)深三維有限元數(shù)值模擬[J].成都理工大學(xué)學(xué)報(bào):自然科學(xué)版，2010，137(6):654-659.

[10]汪宏年，陶宏根，姚敬金，等.用模式匹配算法研究層狀各向異性?xún)A斜地層中多分量感應(yīng)測(cè)井響應(yīng)[J].地球物理學(xué)報(bào)，2008，51(5):1591-1599.

[11]Wang H N，So P，Yang S W，et al.Nummerical Modeling of Multi-component Induction Well Logging Tools in the Cylindrically Stratified Anisotropic Media and its Response[J].IEEE Trans on Geosci and Remote Sensing，2003，46(4):1134-1147.

[12]張燁，汪宏年，陶宏根，等.基于耦合標(biāo)勢(shì)與矢勢(shì)的有限體積法模擬非均勻各向異性地層中多分量感應(yīng)測(cè)井三維響應(yīng)[J].地球物理學(xué)報(bào)，2012，55(6):2141-2152.

[13]仵杰，劉春雅，張?zhí)鹛穑?陣列感應(yīng)測(cè)井儀刻度系數(shù)特性研究[J].測(cè)井技術(shù)，2006，30(5):400-403.

[14]張群華，范宜仁，譚寶海，等.陣列感應(yīng)測(cè)井儀刻度研究[J].測(cè)井技術(shù)，2010，34(6):576-584.