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

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

0 引言

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

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

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

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

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

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

式中，σxx=σyy。

2 有限元計算方法

對于感應(yīng)測井，其電磁場問題可以由Maxwell方程組描述

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

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

令

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

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

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

這樣，對于某個單元可以計算得到單元剛度矩陣為

式中，

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

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

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

3 計算實例與分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

刻度環(huán)的其他參數(shù)不變，只改變刻度環(huán)的傾斜方位XETA，計算出的儀器響應(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)對于淺、深探測的三維感應(yīng)線圈系，其刻度響應(yīng)達到最大的刻度環(huán)的半徑是不同的，因而無法確定使得這2個不同探測深度的線圈系同時達到最大響應(yīng)的刻度環(huán)半徑值，只能折中考慮。為此，可以考慮采用雙環(huán)設(shè)計。

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

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

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

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

[2]白彥，仵杰.三分量感應(yīng)線圈系在均勻地層中的響應(yīng)特性分析[J].電子測試，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)測井響應(yīng)計算的交錯網(wǎng)格有限差分法[J].地球物理學(xué)報，2003，46(4):561-567.

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

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

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

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

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

[10]汪宏年，陶宏根，姚敬金，等.用模式匹配算法研究層狀各向異性傾斜地層中多分量感應(yīng)測井響應(yīng)[J].地球物理學(xué)報，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]張燁，汪宏年，陶宏根，等.基于耦合標勢與矢勢的有限體積法模擬非均勻各向異性地層中多分量感應(yīng)測井三維響應(yīng)[J].地球物理學(xué)報，2012，55(6):2141-2152.

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

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