基于時(shí)域有限差分亞網(wǎng)格與共形網(wǎng)格技術(shù)的隨鉆電磁波測(cè)井響應(yīng)數(shù)值模擬

范宜仁, 李煒,3, 李虎, 王磊, 李栗

(1.中國(guó)石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院, 山東 青島 266580; 2.中國(guó)石油大學(xué)CNPC測(cè)井重點(diǎn)實(shí)驗(yàn)室, 山東 青島 266580; 3.中海油田服務(wù)股份有限公司物探事業(yè)部, 天津 300451)

0 引 言

隨鉆電磁波測(cè)井具有實(shí)時(shí)性、探測(cè)范圍大等特征,被廣泛應(yīng)用于隨鉆地質(zhì)導(dǎo)向和水平井儲(chǔ)層評(píng)價(jià)等方向[1]。Tsili Wang等使用時(shí)域有限差分(FDTD)模擬了傾斜各向同性地層及儀器偏心條件下的測(cè)井響應(yīng)[2],Hwa Ok Lee等使用柱坐標(biāo)系下的FDTD和PML吸收邊界模擬了各向異性地層的隨鉆電磁響應(yīng)[3]。董建黨采用FDTD模擬了三維各向異性地層條件下隨鉆電磁波測(cè)井響應(yīng)[4],其地層模型較為簡(jiǎn)單。范宜仁等研究了儀器偏心條件下三維傾斜地層的測(cè)井響應(yīng)[5]。使用FDTD方法進(jìn)行傳統(tǒng)電測(cè)井?dāng)?shù)值模擬時(shí),網(wǎng)格剖分往往采用階梯近似方法,具有易于實(shí)現(xiàn)、計(jì)算速度快等特點(diǎn),但在處理三維復(fù)雜地層邊界、儀器偏心、線圈、井眼等情況時(shí),為保證計(jì)算精度,需不斷加密網(wǎng)格,計(jì)算速度與存儲(chǔ)量急劇增加。

柴玫、耿軒等研究了二維時(shí)域有限差分亞網(wǎng)格技術(shù),并將其應(yīng)用于帶涂層的局部機(jī)翼、有刻槽的金屬長(zhǎng)方柱等模擬[6-7],大大節(jié)約了計(jì)算內(nèi)存和時(shí)間。張玉、曾碧能等采用共形網(wǎng)格技術(shù)或改進(jìn)的共形網(wǎng)格技術(shù)對(duì)PBG仿真或涂層三角面元的散射等進(jìn)行了模擬,極大地提高了模擬精度[8-9]。但中國(guó)國(guó)內(nèi)對(duì)于將亞網(wǎng)格和共形網(wǎng)格技術(shù)應(yīng)用于隨鉆電磁波測(cè)井響應(yīng)的數(shù)值模擬研究很少。本文把時(shí)域有限差分亞網(wǎng)格和共形網(wǎng)格技術(shù)應(yīng)用于復(fù)雜地層隨鉆電磁波測(cè)井響應(yīng)的數(shù)值模擬,并與單純使用時(shí)域有限差分算法的模擬結(jié)果進(jìn)行對(duì)比,提高了其計(jì)算速度和模擬精度。

1 三維圓柱坐標(biāo)系下FDTD差分格式

模擬的電性各向同性地層為有耗介質(zhì),其本構(gòu)方程為

D=ε·E,J=σ·E

(1)

地層磁性為

B=μ·H,Jm=σm·H

(2)

此時(shí),在無(wú)源區(qū)Maxwell方程微分形式為

(3)

(4)

式中,ε為地層介電常數(shù);μ為地層磁導(dǎo)率;σ為地層電導(dǎo)率;σm為地層導(dǎo)磁率。

將圓柱坐標(biāo)系Yee元胞中的E、H的某一分量設(shè)為f(r,φ,z,t),表示為

f(r,φ,z,t)=f(iΔr,jΔφ,kΔz,nΔt)=fn(i,j,k)

(5)

式中,Δr、Δφ、Δz分別為空間步長(zhǎng);Δt為時(shí)間步長(zhǎng)。采用二階中心差分近似進(jìn)行離散

(6)

按照Yee元胞可以得出FDTD迭代格式(以Er分量為例)

(7)

2 亞網(wǎng)格技術(shù)

FDTD作為模擬隨鉆電磁波傳播的方法其模擬精度主要取決于剖分網(wǎng)格的大小。網(wǎng)格尺寸δ的取值范圍必須滿(mǎn)足色散條件,即δ≤λ/10,其中λ為FDTD模擬區(qū)域內(nèi)介質(zhì)尺寸最小的波長(zhǎng)。δ的大小還需考慮目標(biāo)區(qū)域的某些細(xì)微結(jié)構(gòu),比如隨鉆電磁波測(cè)井儀器的線圈以及井眼,這時(shí)需要對(duì)這些細(xì)微區(qū)域采用較小的網(wǎng)格,否則會(huì)影響模擬結(jié)果的準(zhǔn)確性。對(duì)于這種情況,可以在線圈、井眼等細(xì)微結(jié)構(gòu)處采用亞網(wǎng)格技術(shù),其他FDTD模擬區(qū)域采用粗網(wǎng)格技術(shù)。這樣,可以實(shí)現(xiàn)在保證計(jì)算精度的前提下大量節(jié)省計(jì)算機(jī)的運(yùn)行內(nèi)存和時(shí)間。

2.1 亞網(wǎng)格算法

以軸對(duì)稱(chēng)圓柱坐標(biāo)系中r—φ平面為例,研究三維軸對(duì)稱(chēng)坐標(biāo)系時(shí)域有限差分(FDTD)亞網(wǎng)格技術(shù)。采用粗網(wǎng)格和細(xì)網(wǎng)格同時(shí)計(jì)算的方式,將粗網(wǎng)格和細(xì)網(wǎng)格在時(shí)間步上交替計(jì)算。粗細(xì)網(wǎng)格分界面在切向電場(chǎng)面上,如圖1(a)、(b)所示,粗網(wǎng)格和細(xì)網(wǎng)格的尺寸為奇數(shù)比,這里設(shè)為3∶1,則粗細(xì)網(wǎng)格比nf=3。因?yàn)榇志W(wǎng)格和細(xì)網(wǎng)格的比例為奇數(shù),所以對(duì)于此r—φ平面布局,粗網(wǎng)格的磁場(chǎng)和細(xì)網(wǎng)格的磁場(chǎng)在時(shí)間和空間上均重合,不需要采用插值。

圖1 柱坐標(biāo)系下的亞網(wǎng)格示意圖

設(shè)粗網(wǎng)格的時(shí)間步長(zhǎng)為Δtf,尺寸大小為Δr和Δφ,細(xì)網(wǎng)格的時(shí)間步Δtf=Δtc/nf。不論是粗網(wǎng)格還是亞網(wǎng)格,在網(wǎng)格內(nèi)部,傳統(tǒng)FDTD的計(jì)算都適用。因此,亞網(wǎng)格算法的重點(diǎn)是粗細(xì)網(wǎng)格邊界上E和H的計(jì)算,采用空間上與時(shí)間上的內(nèi)插插值法進(jìn)行計(jì)算。對(duì)于粗網(wǎng)格直接采取FDTD計(jì)算,局部細(xì)網(wǎng)格采取蛙跳方式進(jìn)行計(jì)算。

(8)

式中

(9)

(10)

式中,E0、E1、E2由上步得到。粗細(xì)網(wǎng)格邊界處其他與粗網(wǎng)格不重合點(diǎn)的細(xì)網(wǎng)格電場(chǎng)值亦可通過(guò)式(10)得到。為了提高計(jì)算時(shí)電場(chǎng)的穩(wěn)定性,對(duì)于粗細(xì)網(wǎng)格邊界內(nèi)一層的電場(chǎng)進(jìn)行修正,采用加權(quán)平均的方式。針對(duì)圖1所示,采用公式

(11)

(12)

圖2 亞網(wǎng)格效果圖

測(cè)井儀器、井眼、發(fā)射線圈、接收線圈等處需要精細(xì)刻畫(huà),以保證電磁波傳播模擬更加精準(zhǔn),所以需要在這些地方采用細(xì)網(wǎng)格,在地層等其他地方采用粗網(wǎng)格,其中粗網(wǎng)格與細(xì)網(wǎng)格尺寸比例為3∶1。具體橫截面和徑向的網(wǎng)格剖分見(jiàn)圖2(a)、(b)。圖2(a)中,在r—φ平面上,對(duì)線圈、井眼處采用了細(xì)網(wǎng)格,地層其他地方采用了粗網(wǎng)格;圖2(b)中,在r—z平面上,對(duì)發(fā)射線圈和接收線圈處,采用了細(xì)網(wǎng)格,地層其他地方采用了粗網(wǎng)格。

2.2 計(jì)算速度/精度對(duì)比與實(shí)例

通過(guò)在均勻地層條件下模擬單發(fā)雙收隨鉆電磁波儀器的測(cè)井響應(yīng),對(duì)比分析了傳統(tǒng)FDTD算法和亞網(wǎng)格算法的計(jì)算時(shí)間和精度。取發(fā)射頻率為2 MHz,線圈距為[24 30] in*非法定計(jì)量單位,1 ft=12 in=0.304 8 m,下同,模擬地層電阻率分別取1、10 Ω·m,2種算法的對(duì)比結(jié)果見(jiàn)表1、表2。

表1 R=1 Ω·m模擬結(jié)果精度與時(shí)間

表2 R=10 Ω·m模擬結(jié)果精度與時(shí)間

由傳統(tǒng)FDTD和亞網(wǎng)格技術(shù)的計(jì)算時(shí)間和模擬精度對(duì)比可以發(fā)現(xiàn),亞網(wǎng)格技術(shù)計(jì)算時(shí)間僅為傳統(tǒng)FDTD的19%,亞網(wǎng)格技術(shù)模擬相對(duì)誤差均控制在0.8%之內(nèi)。圖3(a)和(b)為傳統(tǒng)FDTD和亞網(wǎng)格技術(shù)計(jì)算的發(fā)射頻率2 MHz、線圈距[24 30] in條件下的相位差、幅度比刻度圖版,可以看出兩者吻合很好,驗(yàn)證了亞網(wǎng)格技術(shù)的正確性。

圖3 利用亞網(wǎng)格技術(shù)計(jì)算刻度圖版與FDTD對(duì)比

3 共形網(wǎng)格技術(shù)

圖4 傾斜地層切割Yee元胞示意圖

圖5 儀器偏心橫截面示意圖

使用FDTD模擬三維復(fù)雜地層隨鉆電磁波測(cè)井響應(yīng)時(shí),對(duì)于復(fù)雜地層比如傾斜地層(見(jiàn)圖4)或儀器偏心(見(jiàn)圖5)等會(huì)產(chǎn)生階梯誤差。剖分的網(wǎng)格尺寸越小,模擬過(guò)程中產(chǎn)生的階梯誤差越小。但是減少網(wǎng)格尺寸勢(shì)必造成網(wǎng)格總數(shù)的增加和時(shí)間步長(zhǎng)的減少,導(dǎo)致計(jì)算機(jī)計(jì)算時(shí)間和存儲(chǔ)空間的急劇增加。對(duì)于復(fù)雜地層,這種減少網(wǎng)格尺寸的做法并不可取。這一難題可以用共形網(wǎng)格技術(shù)予以解決。

3.1 共形網(wǎng)格技術(shù)原理

圖6 r—z平面共形網(wǎng)格中地層參數(shù)及其等效

設(shè)地層1的電磁參數(shù)為ε1、σ1、μ1、σm1;地層2的電磁參數(shù)為ε2、σ2、μ2、σm2。電場(chǎng)分量位于元胞棱邊的中點(diǎn),根據(jù)Maxwell方程積分形式,電磁參數(shù)等效值可以由加權(quán)平均得到,因此,圖6中A、C處4個(gè)電場(chǎng)分量中的等效ε和等效σ為

(13)

(14)

同理,F處的磁場(chǎng)分量可由不同地層所占的元胞面積加權(quán)平均得到,F處磁場(chǎng)分量中的等效μ和等效σm為

(15)

共形網(wǎng)格中處的FDTD格式與傳統(tǒng)FDTD相同,只需將等效電磁參數(shù)替代普通地層電磁參數(shù),代入FDTD差分格式即可。共形網(wǎng)格技術(shù)的這種簡(jiǎn)單近似,應(yīng)用到介質(zhì)突變分界面等情況,可以很好地減少階梯近似產(chǎn)生的誤差,提高模擬的精度,且不會(huì)增加計(jì)算時(shí)間。圖7(a)、(b)為3層45 °傾斜地層共形網(wǎng)格和階梯近似對(duì)比圖,可以看出傾斜地層的層界面處的電阻率為上下地層電阻率的加權(quán)平均。

圖7 45 °的3層傾斜地層電阻率

3.2 計(jì)算速度/精度對(duì)比與實(shí)例

模擬單發(fā)雙收隨鉆電磁波儀器偏心時(shí)的測(cè)井響應(yīng),對(duì)比分析了傳統(tǒng)FDTD階梯近似和共形網(wǎng)格算法的計(jì)算時(shí)間和精度。取發(fā)射頻率為2 MHz,線圈距為[24 30] in,模擬地層電阻率為0.1 Ω·m,并與階梯近似模擬結(jié)果、Hue一維半數(shù)值解進(jìn)行對(duì)比(見(jiàn)表3,以一維半數(shù)值解為參照)。圖8為共形網(wǎng)格、階梯近似、一維半數(shù)值解3種方法模擬的儀器偏心條件下隨鉆電磁波測(cè)井儀器響應(yīng)特征對(duì)比。

圖8 利用亞網(wǎng)格技術(shù)計(jì)算刻度圖版與FDTD對(duì)比

由表3和圖8可知,儀器偏心條件下隨鉆電磁波測(cè)井響應(yīng)數(shù)值模擬,共形網(wǎng)格技術(shù)在保證計(jì)算時(shí)間的前提下,減少了相對(duì)誤差,精確度較高。

表3 儀器偏心條件下模擬結(jié)果精度與時(shí)間

4 結(jié) 論

(1) 研究了三維柱坐標(biāo)系下亞網(wǎng)格技術(shù)與共形網(wǎng)格技術(shù),將其應(yīng)用于復(fù)雜地層條件隨鉆電磁波測(cè)井響應(yīng)的數(shù)值模擬。亞網(wǎng)格技術(shù)解決了對(duì)井眼、線圈等精細(xì)結(jié)構(gòu)網(wǎng)格剖分過(guò)小、導(dǎo)致CPU計(jì)算時(shí)間過(guò)多的難題;精細(xì)結(jié)構(gòu)處用細(xì)網(wǎng)格剖分,其他結(jié)構(gòu)采用粗網(wǎng)格,其計(jì)算時(shí)間遠(yuǎn)遠(yuǎn)小于傳統(tǒng)FDTD,且相對(duì)誤差控制在0.8%以?xún)?nèi)。

(2) 共形網(wǎng)格技術(shù)解決了復(fù)雜地層邊界或儀器偏心等情況下采用階梯近似產(chǎn)生較大誤差的難題,在保證計(jì)算時(shí)間的條件下大大減少了階梯誤差,提高了精度。

(3) 這2種技術(shù)可用于模擬三維復(fù)雜地層隨鉆電磁波測(cè)井響應(yīng),為大斜度井/水平井評(píng)價(jià)和隨鉆快速反演提供理論指導(dǎo)和技術(shù)支持。

參考文獻(xiàn):

[1] 楊震. 非均勻復(fù)雜地層隨鉆電磁波測(cè)井響應(yīng)研究 [D]. 青島: 中國(guó)石油大學(xué)(華東), 2009.

[2] Wang T, Signorelli J. Finite-difference Modeling of Electromagnetic Tool Response for Logging While Drilling

[J]. Geophysics, 2004, 69(1): 152-160.

[3] Hue Y K, Teixeiral F L, San L E, et al. Modeling of EM Logging Tools in Arbitrary 3-D Borehole Geometries Using PML FDTD [J]. IEEE Geosci Remote Sens Lett, 2005, 2(1): 78-81.

[4] 董建黨. 各向異性地層電磁測(cè)井響應(yīng)的時(shí)域有限差分方法研究 [D]. 青島: 中國(guó)石油大學(xué)(華東), 2008.

[5] 范宜仁, 胡云云, 李虎, 等. 隨鉆電磁波測(cè)井儀器偏心條件下響應(yīng)模擬與分析 [J]. 中國(guó)石油大學(xué)學(xué)報(bào): 自然科學(xué)版, 2014, 38(2): 59-66.

[6] 柴玫, 閻玉波, 葛德彪. 一種可穿越介質(zhì)的二維時(shí)域有限差分亞網(wǎng)格技術(shù) [J]. 西安電子科技大學(xué)學(xué)報(bào): 自然科學(xué)版, 2000, 27(5): 581-584.

[7] 耿軒, 徐金平. R-FDTD與亞網(wǎng)格相結(jié)合技術(shù)及其應(yīng)用 [J]. 東南大學(xué)學(xué)報(bào): 自然科學(xué)版, 2004, 34(2): 161-165.

[8] 張玉, 宋健, 梁昌洪. 并行共形FDTD算法及其在PBG結(jié)構(gòu)仿真中的應(yīng)用 [J]. 電子學(xué)報(bào), 2003, 12(A): 2142-2144.

[9] 曾碧能. 基于FDTD算法的非均勻網(wǎng)格及共形網(wǎng)格技術(shù)的實(shí)現(xiàn) [D]. 成都: 電子科技大學(xué), 2007.