• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    有限元感應(yīng)測井模擬的背景場選擇方法研究

    2015-02-18 08:00:25王健陳浩王秀明張雷
    地球物理學(xué)報(bào) 2015年6期
    關(guān)鍵詞:計(jì)算誤差總場電導(dǎo)率

    王健, 陳浩, 王秀明, 張雷

    1 中國科學(xué)院聲學(xué)研究所 聲場與聲信息國家重點(diǎn)實(shí)驗(yàn)室, 北京 100190 2 中國科學(xué)院大學(xué), 北京 100049

    ?

    有限元感應(yīng)測井模擬的背景場選擇方法研究

    王健1,2, 陳浩1, 王秀明1, 張雷1

    1 中國科學(xué)院聲學(xué)研究所 聲場與聲信息國家重點(diǎn)實(shí)驗(yàn)室, 北京 100190 2 中國科學(xué)院大學(xué), 北京 100049

    在感應(yīng)測井的有限元模擬中,為了消除源的奇異性,一般將總場分解為背景場和散射場.本文定量研究了不同背景場選擇方法對計(jì)算精度的影響.首先,分析了在均勻和徑向分層介質(zhì)中不同背景電導(dǎo)率對長短源距線圈系有限元模擬結(jié)果的影響;其次,利用三層介質(zhì)模型對比了選擇源點(diǎn)附近和線圈系中點(diǎn)附近地層電導(dǎo)率作為背景電導(dǎo)率的結(jié)果.研究結(jié)果表明,如果在選擇背景場時(shí),只考慮源附近的散射場梯度而不同時(shí)考慮源和接收器附近的散射場梯度,計(jì)算誤差明顯增大.在此基礎(chǔ)上提出一種利用Gianzero幾何因子計(jì)算的視電導(dǎo)率作為背景電導(dǎo)率的新的背景場方法.該方法綜合考慮了圍巖、井眼、線圈距等因素,特別是在介質(zhì)分界面處,可有效減小計(jì)算誤差,并取得了滿意的精度.本研究為復(fù)雜環(huán)境下的感應(yīng)測井模擬的背景場選擇提供了指導(dǎo)和依據(jù).

    數(shù)值模擬; 感應(yīng)測井; 有限元; 背景場; 幾何因子

    1 引言

    隨著石油、天然氣及礦藏勘探和開發(fā)需求的增加及測井技術(shù)的發(fā)展,需要應(yīng)用更加精確和快速的數(shù)值模擬技術(shù)對大斜度井或傾斜地層、井洞、裂縫、溶洞和礦體等復(fù)雜條件下的感應(yīng)測井響應(yīng)進(jìn)行數(shù)值模擬,以便對測量到的視電阻率進(jìn)行有效的校正和可靠的反演.

    在三維數(shù)值模擬中,目前主要的方法包括有限差分法(finite difference method:FDM)(Newman et al.,1995;Weiss and Newman, 2002; 汪功禮等, 2003; 沈金松, 2004)、有限體積法(Finite volume method:FVM)(Weiss,2013)和有限元方法(Finite element method :FEM)(Anderson,1983; Biroand Preis,1989; Badea et al.,2001a;Onegova,2010;Mukherjee and Everett,2011; Zhang et al.,2011;王昌學(xué)等, 2013).感應(yīng)測井的激勵(lì)源附近的場強(qiáng)變化非常劇烈,給計(jì)算建模帶來極大的困難(Badea et al., 2001b).為了解決此問題,通常將總場分解成背景場和散射場分別進(jìn)行計(jì)算(Newman and Alumbaugh, 1995).背景場一般可解析求解,它反映了總場中變化比較劇烈的部分,散射場則需要通過數(shù)值模擬進(jìn)行計(jì)算.這樣分解的主要優(yōu)點(diǎn)是:源可以由背景場簡單方便地引入而不需要模擬源的具體形狀;在較少網(wǎng)格條件下比直接模擬源計(jì)算的結(jié)果精度更高.Chang和Anderson(1984)在軸對稱模型中選用泥漿電導(dǎo)率作為背景電導(dǎo)率.Badea等(2001b)比較了有井模型時(shí)分別選擇均勻介質(zhì)和徑向分層介質(zhì)作為背景介質(zhì)對計(jì)算精度的影響.Newman和Alumbaugh(2002)選擇離發(fā)射線圈最近的介質(zhì)為背景介質(zhì).王昌學(xué)等(2006)認(rèn)為背景電導(dǎo)率要盡量接近整個(gè)區(qū)域的平均電導(dǎo)率.總之,已有的模擬計(jì)算中,一般選擇源所在介質(zhì)的電導(dǎo)率作為背景電導(dǎo)率.但是隨著考慮的地層愈加復(fù)雜,這種選擇方法已經(jīng)無法滿足感應(yīng)測井高精度模擬的要求.目前背景場的選擇對數(shù)值模擬的精度和速度影響的系統(tǒng)研究據(jù)我們所知尚未見報(bào)道.尤其是當(dāng)源與接收器之間是非均勻介質(zhì)時(shí),如何綜合考慮場的分布、選擇一個(gè)最佳的背景電導(dǎo)率是一個(gè)亟待解決的問題.

    本文從感應(yīng)測井有限元模擬和背景場的基本原理出發(fā),詳細(xì)研究了不同介質(zhì)模型下,不同的背景場對有限元模擬精度和速度的影響,分析了影響背景場選擇的因素.提出了利用幾何因子得到的視電導(dǎo)率作為有限元模擬的背景電導(dǎo)率及其所具備的優(yōu)點(diǎn).

    2 基礎(chǔ)理論

    由于感應(yīng)測井的頻率比較低,位移電流可以忽略不計(jì)(Moran and Kunz ,1962).電場和磁場所滿足的微分Maxwell方程為

    (1)

    上式中,ω是角頻率,μ0是自由空間的磁導(dǎo)率,σ是地層電導(dǎo)率.電流項(xiàng)J由電流源Js和感應(yīng)電流項(xiàng)σE組成.當(dāng)把介質(zhì)的電導(dǎo)率分解為背景電導(dǎo)率與差值電導(dǎo)率σ=σb+Δσ,而對應(yīng)的場分解為背景場和散射場E=Eb+Es,并將它們代入公式(1)中,得到散射場方程(NewmanandAlumbaugh,1995):

    (2)

    其中Δσ是地層電導(dǎo)率和背景場電導(dǎo)率的差值.Eb是背景場,Es是散射場,E為總場.對于半徑為a放置于z=zs處的載流線圈,它在均勻介質(zhì)中產(chǎn)生的電場的解析式為

    ×J1(λa)J1(λρ)λdλ,

    (3)

    (4)

    n為計(jì)算區(qū)域的外法向方向.

    在本文中,Es可通過有限元方法計(jì)算得到.本文的全部算例用已經(jīng)開發(fā)的矢量有限元程序計(jì)算完成.由于電磁場的有限元模擬具體實(shí)施過程并不是這篇文章的主要目的,因此不詳細(xì)闡述,相關(guān)研究見文獻(xiàn)(Jin, 2002).

    為了說明背景場對數(shù)值模擬結(jié)果的影響,我們考慮一個(gè)放置于二層介質(zhì)模型中的垂直磁偶極子(VerticalMagneticDipole:VMD)產(chǎn)生的電場強(qiáng)度分布.如圖1所示,將磁偶極子源放置在z軸上,且在介質(zhì)分界面(z=0)下方1.5 m處,磁偶極子源的頻率為20 kHz,其中上層介質(zhì)電導(dǎo)率σup=1 S/m,下層介質(zhì)電導(dǎo)率σlow=10 S/m.相應(yīng)的計(jì)算公式由文獻(xiàn)(HardmanandShen, 1986)給出.圖2是所計(jì)算的VMD在二層介質(zhì)模型中激發(fā)的電場的總場和散射場的等高線.從圖中我們可以看到源附近的總場等高線密集,在源附近電場梯度比較大,用有限元模擬電場在網(wǎng)格不夠精細(xì)或者插值函數(shù)階數(shù)不足(尤其是線性插值)時(shí),無法得到足夠精確的結(jié)果.而通過選擇下層介質(zhì)作為背景場得到的散射場在源附近等高線變得稀疏,電場梯度比較小,因此可以采用較為稀疏的網(wǎng)格或低階的插值函數(shù).

    圖1 二層介質(zhì)中的垂直磁偶極子Fig.1 Vertical magnetic dipole in two-layer medium

    圖2 VMD在二層介質(zhì)中激發(fā)的總場和散射場電場強(qiáng)度實(shí)部等高線Fig.2 Contours of real signal for the total and scattered fields excited by VMD in two-layer medium

    雙線圈系是感應(yīng)測井儀器的基本結(jié)構(gòu),因此下面以雙線圈系為例研究背景電導(dǎo)率的不同選擇對模擬精度的影響.模擬中,電流源的頻率f=20 kHz,線圈系的半徑a=0.05 m,長源距的線圈距L=1 m,短源距的線圈距L=0.3 m.有限元模擬中,采用四面體網(wǎng)格,計(jì)算區(qū)域大約為5個(gè)趨膚深度.為了方便對比,相同探測條件下背景電導(dǎo)率改變時(shí)的網(wǎng)格劃分完全一致.

    3 模擬結(jié)果和分析

    3.1 均勻介質(zhì)模型的背景場

    取地層電導(dǎo)率為0.1 S/m,而背景電導(dǎo)率為地層電導(dǎo)率的0.001~1000倍,計(jì)算結(jié)果如圖3所示.從圖3可以看出,背景電導(dǎo)率越接近地層電導(dǎo)率,有限元模擬計(jì)算的精度越高.當(dāng)背景電導(dǎo)率比地層電導(dǎo)率小時(shí),誤差變化平緩且保持在相對較低的水平;當(dāng)背景電導(dǎo)率大于地層電導(dǎo)率時(shí),誤差迅速地增加.因此當(dāng)模型電導(dǎo)率變化不大時(shí),應(yīng)盡可能選擇與模型接近的值作為背景電導(dǎo)率,而模型電導(dǎo)率變化比較劇烈時(shí),取較小的電導(dǎo)率為背景值要相比取大的電導(dǎo)率作為背景電導(dǎo)率的效果好.這樣可以避免引入過高梯度的背景場而可能導(dǎo)致的散射場梯度的增加.從圖3還可以看出,短源距的計(jì)算誤差要比長源距的小,這可能是因?yàn)殚L源距的地方網(wǎng)格比較稀疏.

    3.2 有井存在的徑向分層

    圖4a和b分別給出了短源距和長源距的雙線圈系在淡水泥漿中背景場對有限元計(jì)算結(jié)果的影響.其中井徑為0.2 m,淡水泥漿電導(dǎo)率σmf=0.01 S/m,無限厚均勻地層的電導(dǎo)率σt=0.005~15 S/m.圖中圓圈和三角分別是選用泥漿和地層作為背景場介質(zhì)計(jì)算的相對誤差.由圖4可以看到,當(dāng)?shù)貙与妼?dǎo)率等于泥漿電導(dǎo)率即介質(zhì)是均勻時(shí)誤差等于0,此時(shí)有限元計(jì)算的結(jié)果等同于解析解.短源距時(shí),用泥漿作背景場的計(jì)算誤差小于地層作背景場的計(jì)算誤差.長源距時(shí),用泥漿作背景場的計(jì)算誤差總體要大于用地層作背景場的計(jì)算誤差.這是因?yàn)槎淘淳嗟木€圈系受井眼影響大,源在電導(dǎo)率為泥漿電導(dǎo)率的均勻介質(zhì)中產(chǎn)生的背景場更接近于總場;長源距線圈系受地層的影響更大,地層背景場更接近于總場.從圖中還可以看出,長源距地層背景場的計(jì)算誤差要小于短源距井眼背景場的計(jì)算誤差.這是因?yàn)榧词乖诙淘淳嗲闆r下,由于井眼半徑比較小,泥漿電導(dǎo)率比較低,來自于地層的影響仍然占很大的比重,井眼背景場與總場之間仍然存在著相當(dāng)大的差異.而在長源距情況下,幾乎所有的影響都來自于地層,地層背景場十分接近總場,因而用有限元計(jì)算得到更高精度結(jié)果.

    圖3 無限均勻介質(zhì)中不同背景電導(dǎo)率時(shí)有限元計(jì)算結(jié)果Fig.3 FEM solutions with different background conductivities in infinite uniform medium

    圖4 淡水泥漿時(shí)不同背景場有限元計(jì)算結(jié)果Fig.4 FEM solutions with different background conductivities in fresh mud

    圖5中a和b分別給出了短源距和長源距的雙線圈系在咸水泥漿下背景場對有限元計(jì)算結(jié)果的影響.咸水泥漿電導(dǎo)率σmf=5 S/m,地層電導(dǎo)率σt=0.01~15 S/m.與淡水泥漿的情況類似,短源距時(shí)用泥漿作背景場的誤差更小,而長源距用地層作背景場的誤差更小.同時(shí),長源距以泥漿為背景場的計(jì)算誤差在地層電導(dǎo)率較低時(shí)變得很大,這是因?yàn)殚L源距受到地層的影響比較大,地層電導(dǎo)率比較低,而背景電導(dǎo)率比較高.如圖3在均勻介質(zhì)中,當(dāng)背景電導(dǎo)率大于目標(biāo)電導(dǎo)率時(shí),計(jì)算誤差會(huì)迅速上升.圖4和圖5的結(jié)果表明,在感應(yīng)測井模擬中,需要根據(jù)源距的長短以及井眼大小選擇背景電導(dǎo)率.

    3.3 水平三層模型的背景場方法研究

    為了考察背景場對用有限元模擬感應(yīng)測井儀器穿越介質(zhì)分界面時(shí)的計(jì)算精度的影響,我們分別計(jì)算了長源距雙線圈系在模型1和模型2中的響應(yīng)(圖6).分別選擇圍巖和中間層作為背景場介質(zhì).圖7和圖8分別給出了對應(yīng)于模型1和模型2的視電導(dǎo)率曲線和相對誤差曲線.圖7中,從左圖可以看出二種背景場計(jì)算的視電導(dǎo)率都有較高的精度.但從右圖的相對誤差曲線可以看到,選擇圍巖作為背景電導(dǎo)率時(shí),在圍巖處結(jié)果的相對誤差較小,一般為10-5;而在中間層處計(jì)算的誤差非常大,最大相對誤差超過了10%.這是由于當(dāng)線圈系處于圍巖且距離目標(biāo)層比較遠(yuǎn)時(shí),來自于目標(biāo)層的影響幾乎可以忽略掉,此時(shí)圍巖背景場接近總場.當(dāng)線圈系處于界面或者目標(biāo)層內(nèi)時(shí),目標(biāo)層對其的影響增加,圍巖背景場與總場的差異增大,導(dǎo)致結(jié)果誤差的增大.選擇中間層電導(dǎo)率作為背景電導(dǎo)率時(shí),結(jié)果在圍巖處誤差較大,而在中間層處誤差較小.從圖中還可以看出,圍巖背景場誤差曲線變化劇烈,而目標(biāo)層背景場誤差曲線變化平緩.這主要是因?yàn)閲鷰r背景場的背景電導(dǎo)率較高,目標(biāo)層背景場的背景電導(dǎo)率較低的原因,而較低的背景電導(dǎo)率要好于較高的背景電導(dǎo)率,參見圖3.圖8展示的規(guī)律和圖7基本一致,唯一的區(qū)別是由于中間層電導(dǎo)率較高,圍巖電導(dǎo)率較低,因此中間層背景場誤差變化的更劇烈,而圍巖背景場變化平緩.因此,對感應(yīng)測井需要模擬隨深度變化的地層而言,取固定的背景電導(dǎo)率效果并不好.如果一定要取固定的背景電導(dǎo)率,應(yīng)該取模型中最低的電導(dǎo)率作為背景值.有限元計(jì)算的結(jié)果總是在低阻地層精度較高,在高阻地層的精度較低.這是因?yàn)樵诘妥璧貙訒r(shí)來自于高阻地層的渦流影響很弱,而在高阻地層時(shí),來自于低阻地層的渦流相對強(qiáng)烈些.因此當(dāng)線圈系位于低阻地層時(shí),計(jì)算所選擇的背景場更接近總場,計(jì)算精度也就更高.

    圖6 水平三層模型示意圖Fig.6 Horizontal three-layer model

    圖5 咸水泥漿時(shí)不同背景場有限元計(jì)算結(jié)果Fig.5 FEM solutions with different background conductivities in saltwater mud

    圖7 高阻三層模型視電導(dǎo)率和相對誤差曲線Fig.7 The apparent conductivity and relative error curves of three-layer model with less conductive medium in the middle

    圖8 低阻三層模型視電導(dǎo)率和相對誤差曲線Fig.8 The apparent conductivity and relative error curves of three-layer model with more conductive medium in the middle

    圖9 高阻三層模型線圈系中點(diǎn)背景場和發(fā)射線圈背景場計(jì)算誤差對比Fig.9 Comparison between solutions obtained with the midpoint of transmitter-receiver spacing and transmitter background fields in three-layer model with less conductive medium in the middle

    圖10 低阻三層模型線圈系中點(diǎn)背景場和發(fā)射線圈背景場計(jì)算誤差對比Fig.10 Comparison between solutions obtained with midpoint of transmitter-receiver spacing and transmitter background fields in three-layer model with more conductive medium in the middle

    圖9和圖10分別給出了選用線圈系中點(diǎn)所在介質(zhì)作為背景電導(dǎo)率和發(fā)射線圈所在介質(zhì)作為背景電導(dǎo)率的相對誤差曲線.采用的模型為模型1和模型2.同圖7和圖8對比,我們可以看到這種隨線圈系移動(dòng)改變背景電導(dǎo)率的方法要遠(yuǎn)優(yōu)于選擇固定的某種介質(zhì)作為背景電導(dǎo)率.同時(shí),當(dāng)線圈系穿越地層時(shí),選擇線圈系中點(diǎn)的介質(zhì)作為背景電導(dǎo)率要好于發(fā)射線圈所在介質(zhì)作為背景電導(dǎo)率.這是因?yàn)楫?dāng)線圈中點(diǎn)在界面上時(shí),線圈系受到上下層介質(zhì)的影響是均等的.當(dāng)線圈系繼續(xù)進(jìn)入上一層時(shí),上層介質(zhì)的影響增大,這時(shí)選用上層介質(zhì)作為背景電導(dǎo)率獲得的結(jié)果精度更高.然而當(dāng)線圈系跨越界面時(shí),無論選取那種介質(zhì)作為背景電導(dǎo)率得到的背景場都不可能接近總場,導(dǎo)致數(shù)值模擬結(jié)果在這一區(qū)間誤差不穩(wěn)定,當(dāng)介質(zhì)對比度增大時(shí)會(huì)出現(xiàn)極大的誤差.

    源附近的總場是由源在均勻介質(zhì)中形成的場和介質(zhì)交界面處的反射場和折射場組成.當(dāng)源距離邊界較遠(yuǎn)時(shí),以源在附近均勻介質(zhì)中產(chǎn)生的場作為背景場會(huì)和總場相差不大,且接收器距離源越近,在接收器附近的散射場的梯度也會(huì)越小.然而在感應(yīng)測井問題中,源和接收器要經(jīng)常穿越界面,特別是當(dāng)?shù)貙雍穸缺容^小、線圈距比較大時(shí),源附近的總場與均勻介質(zhì)的場會(huì)相差很大.當(dāng)接收器距離源較遠(yuǎn)時(shí),接收器的場與源的場也有很大的差異,采用源附近的介質(zhì)作為背景場很可能使接收器的散射場有較大的梯度而使誤差變大.因此在感應(yīng)測井?dāng)?shù)值模擬選用背景場時(shí),要同時(shí)考慮源和接收器的散射場分布.正如我們在上面幾個(gè)實(shí)例中看到的,選擇離源最近的介質(zhì)作為背景電導(dǎo)率的方法很多情況下并不是最好的選擇.

    3.4 Gianzero幾何因子背景場

    在上面的多個(gè)算例中,我們選擇的背景電導(dǎo)率都是模型中原有的介質(zhì),這種背景場介質(zhì)的模型是階躍式的.當(dāng)模擬的模型含有多層介質(zhì)時(shí)(尤其是薄層),這種階躍式的背景場會(huì)產(chǎn)生非常大的誤差.由于影響總場的因素很多,包括層厚、井徑、線圈距、以及圍巖、目標(biāo)層和泥漿的電導(dǎo)率差異.這意味著選擇模型中原有介質(zhì)作為背景電導(dǎo)率在整條測井曲線上總會(huì)出現(xiàn)誤差較大的計(jì)算點(diǎn).為了通過引入背景場能夠同時(shí)使得源和接收器附近的散射場梯度減小,即需要同時(shí)考慮源和接收器之間的場的分布并對二者進(jìn)行折中,由此我們提出了用幾何因子計(jì)算得到的視電導(dǎo)率值作為有限元計(jì)算的背景電導(dǎo)率.本文主要選用Gianzero幾何因子,這是因?yàn)榕cGianzero幾何因子相比,Doll幾何因子沒有考慮電磁波的傳播效應(yīng),只能在低電導(dǎo)率、低頻率條件下,反應(yīng)感應(yīng)測井的響應(yīng)特征,限制了方法的適用范圍;Moran幾何因子計(jì)算視電導(dǎo)率時(shí)需要對電導(dǎo)率進(jìn)行積分運(yùn)算,步驟較為復(fù)雜;而Born幾何因子不適用于在均勻地層中描述測井響應(yīng)(仵杰等,2001).Gianzero幾何因子如下式所示(Gianzero and Anderson, 1981):

    (5)

    (6)

    圖11 高阻三層模型Gianzero背景場和線圈中點(diǎn)背景場的計(jì)算誤差對比Fig.11 Comparsion between solutions obtained with Gianzero and midpoint of transmitter-receiver spacing background fields in three-layer model with less conductive medium in the middle

    圖11和圖12分別給出了用Gianzero幾何因子背景場與線圈中點(diǎn)背景場在模型1和模型2中的誤差比較.可以看到,Gianzero幾何因子背景場的計(jì)算誤差在低阻地層比線圈中點(diǎn)背景場稍大,在高阻地層和分界面處要明顯好于線圈中點(diǎn)背景場,最大誤差小于1%,而線圈中點(diǎn)背景場最大誤差為2.4%.雖然Gianzero幾何因子背景場在低阻地層誤差較大,但是其計(jì)算誤差變化平緩,不同位置下都保持在較低的水平,總體可以控制在1%以內(nèi).而目標(biāo)層背景場則在分界面和高阻地層有著非常大的誤差.

    圖13給出在二層介質(zhì)模型中用Gianzero幾何因子背景場和發(fā)射線圈背景場的有限元計(jì)算結(jié)果與解析解的對比.上層介質(zhì)電導(dǎo)率為10S/m,下層介質(zhì)的電導(dǎo)率為0.001S/m.可以看到,在介質(zhì)電導(dǎo)率對比度比較大時(shí),采用發(fā)射線圈背景場的有限元解在界面處誤差非常大,而用Gianzero幾何因子背景場的有限元解誤差很小.因此與傳統(tǒng)的背景場方法相比,這種方法在大電導(dǎo)率對比度下優(yōu)勢更為明顯.

    利用圖14的軸對稱模型,圖15比較Gianzero幾何因子背景場與Badea提出(Badeaetal.,2001b)的余弦背景場的誤差.余弦背景場的背景電導(dǎo)率模型是由井眼和外部地層組成的同心圓柱地層,其計(jì)算得到的背景場是徑向分層地層的解析解(詳細(xì)的計(jì)算公式參見附錄,由于這種背景場需要計(jì)算余弦變換,所以又稱為余弦背景場).本文取模式匹配法(譚茂金等,2007;張雷等,2012)的結(jié)果為標(biāo)準(zhǔn)值,模擬線圈系長度為1m的長源距.從圖15中可以看到,雖然Gianzero幾何因子背景場在圍巖處誤差相對較大,但在目標(biāo)層和分界面處優(yōu)于余弦背景場且總體誤差都較小(一般小于1%),而余弦背景場在目標(biāo)層的誤差大于3%,因此在線圈系探測范圍內(nèi)存在縱向分層時(shí),Gianzero幾何因子計(jì)算的背景場明顯好于余弦背景場.正如前面在分析徑向分層時(shí)所提到的,源距比較大時(shí)井眼的影響遠(yuǎn)不如外部地層的,當(dāng)線圈系在目標(biāo)層時(shí),受到圍巖的顯著影響,此時(shí)考慮了這種影響的Gianzero幾何因子背景場獲得了更高的精度.而余弦背景場只考慮了井的影響卻沒有考慮圍巖的影響.在圍巖處,遠(yuǎn)離目標(biāo)層時(shí)由于余弦背景場相對于Gianzero幾何因子背景場更接近總場,因此精度更高.雖然在線圈系探測范圍內(nèi)不存在縱向分層時(shí),余弦背景場計(jì)算得到的是精確解,但是采用基于Gianzero幾何因子背景場的有限元模擬結(jié)果對多種復(fù)雜地層環(huán)境(如砂泥巖薄互儲(chǔ)層)更加具有普適性.

    圖12 低阻三層模型Gianzero背景場和線圈中點(diǎn)背景場的計(jì)算誤差對比Fig.12 Comparison between solutions obtained with Gianzero and midpoint of transmitter-receiver spacing background fields in three-layer model with more conductive medium in the middle

    圖13 Gianzero背景場與發(fā)射線圈背景場計(jì)算的有限元結(jié)果對比Fig.13 Comparison solutions obtained with Gianzeroand transmitter background fields

    圖14 軸對稱井眼模型圖Fig.14 Axisymmetric model with borehole

    Gianzero幾何因子背景場相對于余弦背景場的優(yōu)勢不僅僅在于計(jì)算的精度,其計(jì)算速度也得到較大的提升.計(jì)算余弦背景場時(shí)無論是采用801點(diǎn)的快速漢克爾函數(shù)變換(Anderson,1983)還是序列外推積分方法(QuadratureIntegrationwithSequenceExtrapolation:QWE)(Key,2012)都需要計(jì)算大量的貝塞爾函數(shù),而影響計(jì)算效率.表1給出了有限元程序分別采用余弦和Gianzero背景場在單個(gè)儀器記錄點(diǎn)所花費(fèi)的時(shí)間.程序采用Gianzero背景場的計(jì)算效率要比用余弦背景場高4倍.使用的計(jì)算機(jī)配置為4核3.1GHz的CPU,內(nèi)存12G.例如,按照表1中單個(gè)記錄點(diǎn)所需要耗費(fèi)的時(shí)間,模擬一段厚度為200m的非均勻模型,深度采樣間隔為0.05m,則總的計(jì)算點(diǎn)數(shù)為4000.采用余弦背景場比采用Gianzero背景場約多花費(fèi)10.5 h.Gianzero方法由于其解可以被認(rèn)為是偶極子源在均勻介質(zhì)中產(chǎn)生的場,有著簡單的解析表達(dá)式而不需要計(jì)算漢克爾函數(shù)變換,因此其所花費(fèi)的時(shí)間相對應(yīng)用快速漢克爾函數(shù)變換所用的時(shí)間完全可以忽略不計(jì).一般情況下,這種假設(shè)在電磁勘探問題中可以滿足精度的要求.

    圖15 Gianzero背景場和余弦背景場計(jì)算的有限元結(jié)果對比Fig.15 Comparison solutions obtained with Gianzero and cosine background fields

    四面體數(shù)自由度總時(shí)間網(wǎng)格剖分求解和后處理剛度矩陣集成(含背景場計(jì)算)余弦背景場30170236056212.712s0.822s2.25s9.642sGianzero背景場3017023605623.262s0.822s2.37s0.203s

    4 結(jié)論

    背景電導(dǎo)率的大小對電磁場計(jì)算精度有較大的影響,背景電導(dǎo)率選擇不當(dāng)會(huì)使計(jì)算結(jié)果嚴(yán)重偏離真實(shí)值,如選擇背景電導(dǎo)率遠(yuǎn)大于介質(zhì)電導(dǎo)率時(shí),計(jì)算誤差會(huì)急劇增加.對于非均勻復(fù)雜的計(jì)算模型,如砂泥巖薄互層,選擇固定背景電導(dǎo)率難以保證源-接收器移動(dòng)到不同地層,尤其當(dāng)源-接收器之間存在多個(gè)薄層時(shí)的精度.

    理想的背景場選擇要同時(shí)考慮源和接收器附近的散射場分布,這樣可以避免單純地考慮源附近散射場梯度減小而可能導(dǎo)致的接收器附近散射場的梯度增大.文中徑向分層和水平分層的例子驗(yàn)證了這種考慮的重要性,由此提出了基于Gianzero幾何因子的背景場電導(dǎo)率計(jì)算方法.與以往的階躍式的背景場不同,這種幾何因子背景場更加平滑.它同時(shí)考慮了源和接收器的間距、井眼和圍巖等影響因素.盡管在一些高導(dǎo)地層,Gianzero計(jì)算的的背景場的計(jì)算誤差較其他方法大,但維持在很低的水平.當(dāng)源和接收器處于分界面附近時(shí),這種背景場更加接近于總場,避免了介質(zhì)分界面處出現(xiàn)較大的誤差,因此它的計(jì)算誤差曲線更加平緩.數(shù)值模擬的例子證明了這種方法的可靠性和優(yōu)越性.本文所提出的背景場方法可以有效地提高感應(yīng)測井儀器正演模擬的精度和速度,尤其適用于砂泥巖薄互層等復(fù)雜的多層介質(zhì)模型,進(jìn)而為將來的快速反演提供堅(jiān)實(shí)的基礎(chǔ).針對于這種方法在水平井,大斜度井和傾斜地層中的應(yīng)用效果的精度評估有待于進(jìn)一步驗(yàn)證.

    附錄A 有限尺寸載流線圈源在井外無限厚地層中電場強(qiáng)度公式推導(dǎo)

    井眼半徑b,泥漿電導(dǎo)率σ0,外部地層電導(dǎo)率σ1,(ρ,z)是考察點(diǎn)坐標(biāo),電流強(qiáng)度等于I,載流線圈半徑為a.

    (A1)

    a≤ρ≤b,

    Anderson W L. 1983. Fourier Cosine and Sine Transforms Using Lagged Convolutions in Double-precision (Subprograms DLAGF0/DLAGF1).//US Geological SurveyOpenFile Report, 83-320.Badea E A, EverettM E, NewmanG A, et al. 2001a. Finite-element analysis of controlled-source electromagnetic induction using Coulomb-gauged potentials.Geophysics, 66(3): 786-799.

    Badea E A, Everett M E,ShenL C,et al. 2001b. Effect of background fields on three-dimensional finite element analysis of induction logging.RadioScience, 36(4): 721-729.

    Biro O,Preis K. 1989. On the use of the magnetic vector potential in the finite-element analysis of three-dimensional eddy currents.IEEETransactionsMagnetics,25(4): 3145-3159.

    Chang S K, AndersonB. 1984. Simulation of induction logging by the finite-element method.Geophysics,49(11): 1943-1958.

    Druskin V L, KnizhnermanL A, Lee P. 1999. New spectral Lanczos decomposition method for induction modeling in arbitrary 3-D geometry.Geophysics, 64(3): 701-706.

    Gianzero S, AndersonB I. 1981. A new look at skin effect. ∥SPWLA 22nd Annual Logging Symposium, Society of Petrophysicists and Well-Log Analysts.Guptasarma D, Singh B. 1997. New digital linear filters for HankelJ0andJ1transforms.GeophysicalProspecting, 45(5): 745-762.

    Hardman R H, ShenL C. 1986. Theory of induction sonde in dipping beds.Geophysics, 51(3): 800-809.

    Jin J M. 2002. The Finite Element Method in Electromagnetic. New York: Wiley.

    Key, K. 2012. Is the fast Hankel transform faster than quadrature?.Geophysics, 77(3): F21-F30.

    Onegova E V. 2010. Effect of multicoil electromagnetic tool eccentricity on measured signals.RussianGeologyandGeophysics, 51(4): 423-427.Moran J, Kunz K S. 1962. Basic theory of induction logging and application to study of two-coil sondes.Geophysics, 27(6): 829-858.Mukherjee S, Everett M E. 2011. 3D controlled-source electromagnetic edge-based finite element modeling of conductive and permeable heterogeneities.Geophysics, 76(4): F215-F226.Newman G A,Alumbaugh D L. 1995. Frequency-domain modelling of airborne electromagnetic responses using staggered finite differences.GeophysicalProspecting, 43(8): 1021-1042.

    Newman G A, Alumbaugh D L. 2002. Three-dimensional induction logging problems, Part 2: A finite-difference solution.Geophysics, 67(2): 484-491.Shen J S. 2004. Modeling of the multi-component induction log in anisotropic medium by using finite difference method.ProgressinGeophysics(in Chinese), 19(1): 101-107.

    Tan M J, Zhang G J, Yun H Y, et al. 2007. 3-D numerical mode-matching (NMM) method for resistivity logging responses in nonsymmetricconditions.ChineseJournalofGeophysics(in Chinese), 50(3): 939-945.

    Wang C X, Zhou C C, Chu Z T, et al. 2006. Modeling of

    electromagnetic response in frequency domain to electrical anisotropic formations.ChineseJournalofGeophysics(in Chinese), 49(6): 1873-1883.

    Wang C X, Chu Z T, Xiao C W, et al. 2013. The analysis of effect of the borehole environment on response of array induction logging.ChineseJournalofGeophysics(in Chinese), 56(4): 1392-1403.

    Wang G L, Zhang G J, Cui F X, et al. 2003. Application of staggered Grid Finite Difference Method to the Computation of 3-D Induction Logging Response.ChineseJournalofGeophysics(in Chinese), 46(4): 561-568.Weiss C J, NewmanG A. 2002. Electromagnetic induction in a fully 3-D anisotropic earth.Geophysics, 67(4): 1104-1114.Weiss C J. 2013. Project APhiD: A Lorenz-gauged A-Φ decomposition for parallelized computation of ultra-broadband electromagnetic induction in a fully heterogeneous.Earth.Computers&Geosciences, 58: 40-52.

    Zhang L, Chen H, Wang X M. 2012. Numerical modeling of responses to a tilted-coil antenna in a transversely isotropic formation.ChineseJournalofGeophysics(in Chinese), 55(10):3493-3505.

    Zhang Z Q, Zhang X, Mu L X. 2011. Simulation of electromagnetic logging-while-drilling tools using vector finite element methods. ∥ 2011 IEEE International Symposium on Antennas and Propagation (APSURSI).Spokane, WA: IEEE, 2011: 2499-2502.

    附中文參考文獻(xiàn)

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

    譚茂金, 張庚驥, 運(yùn)華云等. 2007. 非軸對稱條件下用三維模式匹配法計(jì)算電阻率測井響應(yīng). 地球物理學(xué)報(bào), 50(3): 939-945.

    王昌學(xué), 周燦燦, 儲(chǔ)昭坦等. 2006. 電性各向異性地層頻率域電磁響應(yīng)模擬. 地球物理學(xué)報(bào), 49(6): 1873-1883.

    王昌學(xué), 儲(chǔ)昭坦, 肖承文等. 2013. 井環(huán)境對陣列感應(yīng)測井響應(yīng)的影響分析. 地球物理學(xué)報(bào), 56(4): 1392-1403.

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

    仵杰,龐巨豐,徐景碩. 2001. 感應(yīng)測井幾何因子理論及其應(yīng)用研究.測井技術(shù), 25(6):417-422

    張雷, 陳浩, 王秀明. 2012. 橫向各向同性地層中傾斜線圈系響應(yīng)特征的數(shù)值模擬. 地球物理學(xué)報(bào), 55(10): 3493-3505.

    (本文編輯 胡素芳)

    Research on selection method of background field for finite element simulation of induction logging

    WANG Jian1,2, CHEN Hao1, WANG Xiu-Ming1, ZHANG Lei1

    1StateKeyLaboratoryofAcoustics,InstituteofAcoustics,ChineseAcademyofSciences,Beijing100190,China2UniversityofChineseAcademyofSciences,Beijing100049,China

    With the development of petroleum industry, accurate and fast numerical simulation methods are needed to achieve reliable inversion of the measured response of induction logging in complex media. Thereinto,the choice of the background field plays an important role in the forward modeling of electromagnetic.Thus, the objectives are to quantitatively compare the effects of different background fields on the calculation accuracy.The finite element method is employed to carry out simulation of induction logging. In order to eliminate the singularities, total field is usually decomposed into two parts, background field and scattered field. Analytical method is generally available for the former one, which reflects the drastic change of the field. While for the latter one, due to its character of gentle change, finite element method with sparse grid would be enough. We performed comparisons between results obtained from different background conductivities in several models to find out how background field affects the numerical accuracy.Then, we propose a new method in which use the apparent conductivity calculated by Gianzero geometry factor as the background conductivity.The main results are as following: (1) in the homogeneous media, error becomes smaller as the background conductivity matches the media′s. And taking smaller background conductivity usually gains better results; (2) in the radial layered medium, when the mud is fresh and using short spacing, mud and formation are both used as background field in the calculation. And results turn out the error of former is smaller than the latter′s. When the mud is salty, similar results are obtained except the scaleof error.(3) in the three layered model, we choose formation conductivity near the source and near the midpoint of transmitter-receiver as background conductivity respectively. The results reveal that error greatly increases without taking scattering field gradient near the source and receiver into consideration simultaneously. (4) we compare the results from Gianzero geometry factor background and cosine background, and the test shows that Gianzero background field has higher precision with long spacing or thin formation.Moreover, Gianzero method is more efficiency due to its concise expression and no need to calculate hankel function transform.In the simulation of induction logging, background field exerts a significant influence upon the calculation accuracy, and the results can deviate markedly due to the improper conductivity value. For the complex inhomogeneous earth model, such as sand-shale interbedded layers, constant background conductivity cannot ensure the accuracy when the transmitter and receiver are not in the same layer, especially there are thin layers.Besides, ideal background field should take the distribution of field near the transmitter and receiver into account simultaneously. Only in this way can the increase of gradient of scattered field near the receiver be avoided, which is caused by only considering the decrease of gradient of field near the transmitter.The main differences between Gianzero geometry factor background and other step type background are that the Gianzero geometry factor background is more smooth and gentle, and can take many factors into account, such as surrounding rock, borehole, transmitter-receiver interval. When the tool intersects bed boundary ,the results of Gianzero background field are more accurate.

    Numerical simulation; Induction logging; Finite element method; Background field; Geometrical factor

    10.6038/cjg20150630.

    國家重大科研裝備研制項(xiàng)目"深部資源探測核心裝備研發(fā)"(ZDYZ2012-1-07)資助.

    王健,男,1987年生,博士生,主要從事電磁波場數(shù)值模擬研究. E-mail: wangjianshinian@163.com

    10.6038/cjg20150630

    P631

    2014-06-04,2014-12-23收修定稿

    王健, 陳浩, 王秀明等. 2015. 有限元感應(yīng)測井模擬的背景場選擇方法研究.地球物理學(xué)報(bào),58(6):2177-2187,

    Wang J, Chen J, Wang X M, et al. 2015. Research on selection method of background field for finite element simulation of induction logging.ChineseJ.Geophys. (in Chinese),58(6):2177-2187,doi:10.6038/cjg20150630.

    猜你喜歡
    計(jì)算誤差總場電導(dǎo)率
    炭黑填充天然橡膠超彈性本構(gòu)方程的適用性分析
    綜合施策打好棉花田管“組合拳”
    水尺計(jì)重中密度測量與計(jì)算誤差分析及相關(guān)問題的思考
    水尺計(jì)重中密度測量與計(jì)算誤差分析及相關(guān)問題的思考
    基于比較測量法的冷卻循環(huán)水系統(tǒng)電導(dǎo)率檢測儀研究
    低溫脅迫葡萄新梢電導(dǎo)率和LT50值的研究
    前向雷達(dá)目標(biāo)回波成分與特性分析
    石總場早播棉花出苗顯行
    強(qiáng)度折減法中折減參數(shù)對邊坡穩(wěn)定性計(jì)算誤差影響研究
    高電導(dǎo)率改性聚苯胺的合成新工藝
    国产精品嫩草影院av在线观看| 少妇高潮的动态图| 亚洲成色77777| 一级,二级,三级黄色视频| 亚洲美女视频黄频| 一二三四中文在线观看免费高清| 欧美日韩亚洲高清精品| 在线天堂中文资源库| 91午夜精品亚洲一区二区三区| 久久久久视频综合| 大片免费播放器 马上看| 婷婷成人精品国产| 捣出白浆h1v1| 一区在线观看完整版| 久久ye,这里只有精品| 老司机影院成人| 国产日韩欧美亚洲二区| 你懂的网址亚洲精品在线观看| 国产国语露脸激情在线看| 天堂8中文在线网| 免费黄色在线免费观看| 中国美白少妇内射xxxbb| 日本黄大片高清| 国产精品秋霞免费鲁丝片| 啦啦啦中文免费视频观看日本| 日韩 亚洲 欧美在线| 国产老妇伦熟女老妇高清| 久久99一区二区三区| 日本免费在线观看一区| 欧美人与善性xxx| 精品人妻在线不人妻| 久久精品国产鲁丝片午夜精品| 欧美丝袜亚洲另类| 熟妇人妻不卡中文字幕| 国产高清三级在线| 国产午夜精品一二区理论片| 日韩,欧美,国产一区二区三区| 亚洲成人av在线免费| a级毛片在线看网站| 美女主播在线视频| 色婷婷av一区二区三区视频| 18禁裸乳无遮挡动漫免费视频| 精品一区在线观看国产| 91在线精品国自产拍蜜月| 又粗又硬又长又爽又黄的视频| 男人舔女人的私密视频| 熟女人妻精品中文字幕| 精品人妻在线不人妻| 伊人亚洲综合成人网| 又粗又硬又长又爽又黄的视频| 一本大道久久a久久精品| 亚洲少妇的诱惑av| 欧美xxⅹ黑人| 97在线视频观看| 国产精品久久久久久精品电影小说| 国产精品麻豆人妻色哟哟久久| 99国产综合亚洲精品| 日韩制服骚丝袜av| 2021少妇久久久久久久久久久| 午夜福利视频在线观看免费| 午夜福利网站1000一区二区三区| 亚洲内射少妇av| 男人爽女人下面视频在线观看| 久久国产精品男人的天堂亚洲 | 大片免费播放器 马上看| 久久精品国产a三级三级三级| 亚洲中文av在线| 中文字幕人妻熟女乱码| 亚洲精品av麻豆狂野| 18禁动态无遮挡网站| 国产一区二区激情短视频 | 天美传媒精品一区二区| 亚洲精品国产av蜜桃| 成年人免费黄色播放视频| 免费观看a级毛片全部| 看免费av毛片| 热re99久久精品国产66热6| 国产又色又爽无遮挡免| 国产精品成人在线| 99久久综合免费| 国产精品蜜桃在线观看| 你懂的网址亚洲精品在线观看| 精品久久蜜臀av无| 久久人人爽人人片av| 久久亚洲国产成人精品v| 伊人亚洲综合成人网| 国产亚洲最大av| 欧美成人午夜精品| 亚洲经典国产精华液单| 五月玫瑰六月丁香| 国产 精品1| 赤兔流量卡办理| 蜜臀久久99精品久久宅男| 人人澡人人妻人| 黄色 视频免费看| www.av在线官网国产| 26uuu在线亚洲综合色| 毛片一级片免费看久久久久| 国产又爽黄色视频| 1024视频免费在线观看| 日日爽夜夜爽网站| 在线观看一区二区三区激情| 国产精品成人在线| 午夜老司机福利剧场| 日本91视频免费播放| 永久网站在线| 国产日韩一区二区三区精品不卡| 亚洲精品456在线播放app| 久久精品人人爽人人爽视色| 久久人人爽人人爽人人片va| 热re99久久精品国产66热6| 黄色一级大片看看| 女人精品久久久久毛片| 狂野欧美激情性bbbbbb| 久久久久国产网址| 五月玫瑰六月丁香| 欧美亚洲日本最大视频资源| 久久久久久久久久人人人人人人| 黑人猛操日本美女一级片| 人妻系列 视频| 久久久久久久久久成人| 精品一区二区免费观看| 少妇人妻精品综合一区二区| 国产精品.久久久| 成人国产av品久久久| 国产高清国产精品国产三级| 亚洲国产成人一精品久久久| 亚洲欧美一区二区三区黑人 | 欧美精品亚洲一区二区| 日韩在线高清观看一区二区三区| 免费在线观看黄色视频的| 国产在视频线精品| 热99久久久久精品小说推荐| 国产熟女午夜一区二区三区| 亚洲,一卡二卡三卡| 乱人伦中国视频| 免费女性裸体啪啪无遮挡网站| 国产国拍精品亚洲av在线观看| 免费黄色在线免费观看| 老司机影院毛片| 国产一区亚洲一区在线观看| 另类精品久久| 一级毛片黄色毛片免费观看视频| 久久影院123| 中文字幕亚洲精品专区| 男女下面插进去视频免费观看 | 90打野战视频偷拍视频| 夫妻午夜视频| 日韩制服丝袜自拍偷拍| 日韩av不卡免费在线播放| 日韩成人av中文字幕在线观看| 成年av动漫网址| 你懂的网址亚洲精品在线观看| 欧美国产精品一级二级三级| 久久人人97超碰香蕉20202| 中文字幕人妻丝袜制服| 一边亲一边摸免费视频| 纯流量卡能插随身wifi吗| 熟女电影av网| xxxhd国产人妻xxx| 男人添女人高潮全过程视频| 日韩中字成人| 亚洲欧美中文字幕日韩二区| 五月开心婷婷网| 爱豆传媒免费全集在线观看| 成人二区视频| 寂寞人妻少妇视频99o| www.色视频.com| 精品国产露脸久久av麻豆| 精品国产一区二区三区四区第35| 中国美白少妇内射xxxbb| 人成视频在线观看免费观看| 观看av在线不卡| 少妇人妻久久综合中文| 最近中文字幕高清免费大全6| 欧美少妇被猛烈插入视频| 国产成人av激情在线播放| 9191精品国产免费久久| 中文字幕亚洲精品专区| 久久久久精品人妻al黑| 国产精品成人在线| 久久久久人妻精品一区果冻| 下体分泌物呈黄色| 精品国产一区二区三区四区第35| 免费观看在线日韩| 色婷婷av一区二区三区视频| 99久久综合免费| 久久人人97超碰香蕉20202| 在线天堂中文资源库| 亚洲色图 男人天堂 中文字幕 | 一区二区三区乱码不卡18| 免费少妇av软件| 2022亚洲国产成人精品| 大片免费播放器 马上看| 秋霞伦理黄片| 日本vs欧美在线观看视频| 欧美精品亚洲一区二区| 国产精品久久久久久久电影| 高清不卡的av网站| 热99久久久久精品小说推荐| 精品午夜福利在线看| 激情视频va一区二区三区| 波野结衣二区三区在线| 9色porny在线观看| 人妻人人澡人人爽人人| 国产爽快片一区二区三区| 成人毛片a级毛片在线播放| 在线观看美女被高潮喷水网站| 高清在线视频一区二区三区| 亚洲精品,欧美精品| 国语对白做爰xxxⅹ性视频网站| av视频免费观看在线观看| 五月伊人婷婷丁香| 日韩一区二区三区影片| 永久网站在线| 王馨瑶露胸无遮挡在线观看| 少妇被粗大的猛进出69影院 | 国产成人aa在线观看| 看十八女毛片水多多多| 青春草亚洲视频在线观看| 久久精品国产亚洲av天美| 母亲3免费完整高清在线观看 | 久久久久网色| 国产av码专区亚洲av| 制服人妻中文乱码| 成年动漫av网址| 人人妻人人澡人人看| 97人妻天天添夜夜摸| 亚洲精品美女久久久久99蜜臀 | 久久久国产欧美日韩av| 亚洲欧洲国产日韩| 在线免费观看不下载黄p国产| 老司机影院毛片| 亚洲国产精品一区三区| 天天影视国产精品| 日本vs欧美在线观看视频| 少妇熟女欧美另类| 亚洲欧美一区二区三区黑人 | 狠狠精品人妻久久久久久综合| 大话2 男鬼变身卡| 一区二区三区乱码不卡18| 午夜免费男女啪啪视频观看| 日韩三级伦理在线观看| 高清毛片免费看| 久久99一区二区三区| 草草在线视频免费看| 91在线精品国自产拍蜜月| 中文字幕免费在线视频6| 国产免费现黄频在线看| 搡女人真爽免费视频火全软件| 国产免费视频播放在线视频| 99热6这里只有精品| 女人被躁到高潮嗷嗷叫费观| 尾随美女入室| 欧美精品一区二区免费开放| 午夜免费鲁丝| av在线老鸭窝| freevideosex欧美| 最后的刺客免费高清国语| 精品国产一区二区三区四区第35| 天堂中文最新版在线下载| 亚洲国产欧美在线一区| 中文字幕精品免费在线观看视频 | 69精品国产乱码久久久| 一区二区三区四区激情视频| 日韩电影二区| 日韩不卡一区二区三区视频在线| 狂野欧美激情性bbbbbb| 久久久久国产网址| 少妇的逼好多水| 中文字幕免费在线视频6| 美女大奶头黄色视频| 黄色毛片三级朝国网站| 人妻 亚洲 视频| 精品一区二区免费观看| 国产成人aa在线观看| 欧美日本中文国产一区发布| 精品国产一区二区久久| 咕卡用的链子| 欧美97在线视频| 大香蕉久久成人网| 一本久久精品| 国产毛片在线视频| 国产成人免费观看mmmm| 国产乱来视频区| 国产熟女欧美一区二区| 丰满乱子伦码专区| 久久精品国产a三级三级三级| 久久久亚洲精品成人影院| 天美传媒精品一区二区| 免费人成在线观看视频色| 看非洲黑人一级黄片| 久久精品夜色国产| 亚洲成色77777| 国产精品免费大片| a级毛色黄片| 十八禁网站网址无遮挡| 日韩制服骚丝袜av| 韩国高清视频一区二区三区| 免费高清在线观看日韩| 成人漫画全彩无遮挡| 久久久a久久爽久久v久久| 制服诱惑二区| 国产成人aa在线观看| videosex国产| 有码 亚洲区| 久久久久久久大尺度免费视频| 国产有黄有色有爽视频| 亚洲精品国产色婷婷电影| 女人被躁到高潮嗷嗷叫费观| 国产免费一级a男人的天堂| 午夜免费鲁丝| 欧美精品高潮呻吟av久久| 国产精品国产三级专区第一集| 18禁在线无遮挡免费观看视频| 久久精品熟女亚洲av麻豆精品| 91国产中文字幕| 99久久综合免费| 国产伦理片在线播放av一区| 美女xxoo啪啪120秒动态图| av视频免费观看在线观看| 国产男人的电影天堂91| 高清在线视频一区二区三区| 中文欧美无线码| 亚洲欧美日韩另类电影网站| 国产精品国产三级专区第一集| 日韩 亚洲 欧美在线| 满18在线观看网站| 男人爽女人下面视频在线观看| 亚洲第一区二区三区不卡| 大片免费播放器 马上看| 精品酒店卫生间| 啦啦啦中文免费视频观看日本| 色5月婷婷丁香| 在线观看国产h片| 午夜免费男女啪啪视频观看| 久久午夜综合久久蜜桃| 18禁国产床啪视频网站| 免费大片18禁| 欧美人与善性xxx| 9色porny在线观看| 不卡视频在线观看欧美| 黄色毛片三级朝国网站| 久久女婷五月综合色啪小说| 国产av国产精品国产| 捣出白浆h1v1| 国产免费视频播放在线视频| 少妇高潮的动态图| 久久午夜福利片| av有码第一页| 久久精品国产综合久久久 | 亚洲四区av| 免费观看av网站的网址| 亚洲精品日本国产第一区| 中文字幕人妻丝袜制服| 色婷婷久久久亚洲欧美| 国产永久视频网站| 下体分泌物呈黄色| 亚洲精品日韩在线中文字幕| 在线免费观看不下载黄p国产| 久久久精品免费免费高清| 国产亚洲午夜精品一区二区久久| 国产成人午夜福利电影在线观看| 久久久久人妻精品一区果冻| 99re6热这里在线精品视频| 欧美日韩视频高清一区二区三区二| 免费在线观看完整版高清| 欧美精品av麻豆av| 久久av网站| 国产日韩一区二区三区精品不卡| 亚洲精品美女久久久久99蜜臀 | 亚洲av日韩在线播放| 老司机影院成人| 午夜福利网站1000一区二区三区| 久久精品夜色国产| 街头女战士在线观看网站| 青春草国产在线视频| 蜜臀久久99精品久久宅男| 三级国产精品片| 日产精品乱码卡一卡2卡三| 少妇 在线观看| 亚洲av国产av综合av卡| 建设人人有责人人尽责人人享有的| 2022亚洲国产成人精品| 亚洲精品456在线播放app| 中文字幕人妻熟女乱码| 国产乱人偷精品视频| 欧美精品一区二区免费开放| 飞空精品影院首页| 久久精品久久久久久噜噜老黄| 自线自在国产av| 视频区图区小说| 内地一区二区视频在线| 国产精品国产三级国产av玫瑰| 在线观看人妻少妇| 亚洲内射少妇av| 赤兔流量卡办理| 在线观看国产h片| 人人妻人人爽人人添夜夜欢视频| 日日摸夜夜添夜夜爱| 亚洲av电影在线进入| 18禁在线无遮挡免费观看视频| 亚洲av欧美aⅴ国产| 久久午夜福利片| 欧美日韩av久久| 中文字幕免费在线视频6| 一区二区三区精品91| 下体分泌物呈黄色| 1024视频免费在线观看| 曰老女人黄片| 永久免费av网站大全| 国产在视频线精品| 免费观看av网站的网址| 热99国产精品久久久久久7| 制服丝袜香蕉在线| 在线 av 中文字幕| 欧美精品人与动牲交sv欧美| 少妇人妻精品综合一区二区| 街头女战士在线观看网站| 99热6这里只有精品| 免费高清在线观看视频在线观看| 国产精品一二三区在线看| 成人手机av| 欧美丝袜亚洲另类| 久久精品国产自在天天线| 精品99又大又爽又粗少妇毛片| 亚洲欧美成人综合另类久久久| 18禁观看日本| 极品人妻少妇av视频| 高清毛片免费看| 最近最新中文字幕免费大全7| av在线老鸭窝| 一级a做视频免费观看| 80岁老熟妇乱子伦牲交| 久久综合国产亚洲精品| 少妇高潮的动态图| 亚洲综合精品二区| 天天操日日干夜夜撸| 欧美xxxx性猛交bbbb| 免费高清在线观看视频在线观看| 国产欧美日韩一区二区三区在线| 亚洲av成人精品一二三区| 97在线人人人人妻| 国产成人精品婷婷| 蜜桃在线观看..| 免费大片18禁| 成人免费观看视频高清| 一本大道久久a久久精品| 午夜福利在线观看免费完整高清在| www日本在线高清视频| 男的添女的下面高潮视频| 亚洲五月色婷婷综合| 丝袜人妻中文字幕| 另类精品久久| 免费在线观看黄色视频的| 乱人伦中国视频| 亚洲av成人精品一二三区| 久久久久国产精品人妻一区二区| 精品99又大又爽又粗少妇毛片| 18禁动态无遮挡网站| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 午夜av观看不卡| 精品国产一区二区三区久久久樱花| 日韩大片免费观看网站| 日韩 亚洲 欧美在线| 人妻一区二区av| 久久这里只有精品19| 亚洲精品第二区| 国产欧美亚洲国产| 久久久国产一区二区| 三上悠亚av全集在线观看| 亚洲国产精品一区二区三区在线| 综合色丁香网| 亚洲第一av免费看| 最新的欧美精品一区二区| 国产黄色视频一区二区在线观看| 91aial.com中文字幕在线观看| 久久久久久久国产电影| 777米奇影视久久| 国产深夜福利视频在线观看| www.色视频.com| 久久婷婷青草| 在线观看美女被高潮喷水网站| 亚洲中文av在线| 夜夜骑夜夜射夜夜干| 22中文网久久字幕| 亚洲av综合色区一区| 草草在线视频免费看| 国产成人精品一,二区| 国产片内射在线| 免费看光身美女| 丰满少妇做爰视频| 国产永久视频网站| 亚洲精品国产av蜜桃| 亚洲一级一片aⅴ在线观看| 国产精品国产av在线观看| 日本爱情动作片www.在线观看| 90打野战视频偷拍视频| 久久久久网色| 欧美精品一区二区免费开放| 伊人久久国产一区二区| 高清av免费在线| 亚洲在久久综合| av国产久精品久网站免费入址| 纵有疾风起免费观看全集完整版| 精品国产国语对白av| 18+在线观看网站| 亚洲国产av影院在线观看| www日本在线高清视频| www.av在线官网国产| 亚洲欧美成人综合另类久久久| 成年人午夜在线观看视频| 国产一区二区激情短视频 | 久久狼人影院| 欧美亚洲 丝袜 人妻 在线| 99热全是精品| 欧美日韩成人在线一区二区| 精品一区在线观看国产| 欧美97在线视频| 国产亚洲欧美精品永久| 午夜福利视频在线观看免费| 国产黄频视频在线观看| 2021少妇久久久久久久久久久| 久久精品国产亚洲av涩爱| 美女xxoo啪啪120秒动态图| 熟女av电影| 青春草亚洲视频在线观看| 国产成人午夜福利电影在线观看| 久久ye,这里只有精品| 秋霞伦理黄片| 亚洲欧洲日产国产| 两个人免费观看高清视频| 国产片内射在线| 亚洲,欧美精品.| 久久久精品94久久精品| 午夜免费观看性视频| 日韩三级伦理在线观看| 精品亚洲成国产av| 国产一区二区在线观看日韩| 丝袜脚勾引网站| 免费在线观看完整版高清| 欧美日韩精品成人综合77777| 午夜福利视频精品| 欧美老熟妇乱子伦牲交| 国产精品人妻久久久影院| av线在线观看网站| 久久久久精品久久久久真实原创| 久久久久久人妻| 国产女主播在线喷水免费视频网站| av国产精品久久久久影院| 两个人免费观看高清视频| 成人黄色视频免费在线看| 五月天丁香电影| 狠狠婷婷综合久久久久久88av| 国产av国产精品国产| 熟女av电影| 中国国产av一级| 亚洲第一av免费看| 日韩大片免费观看网站| 成年av动漫网址| 99香蕉大伊视频| 亚洲精品456在线播放app| 免费高清在线观看视频在线观看| 9热在线视频观看99| 欧美日韩一区二区视频在线观看视频在线| 一区二区三区四区激情视频| 日韩在线高清观看一区二区三区| 免费av中文字幕在线| 秋霞伦理黄片| 一级片'在线观看视频| 内地一区二区视频在线| 一区二区av电影网| 久久久久久人人人人人| 男人操女人黄网站| 一二三四在线观看免费中文在 | 久久青草综合色| 国产av码专区亚洲av| 欧美成人午夜精品| 一级毛片 在线播放| 亚洲性久久影院| 国产精品久久久久成人av| 香蕉国产在线看| 国产又爽黄色视频| 在线观看免费高清a一片| 人成视频在线观看免费观看| 亚洲精品国产av蜜桃| 亚洲精品国产色婷婷电影| 你懂的网址亚洲精品在线观看| 欧美丝袜亚洲另类| 一边亲一边摸免费视频| 成人免费观看视频高清| 三级国产精品片| 免费不卡的大黄色大毛片视频在线观看| 伦精品一区二区三区| 热99国产精品久久久久久7| 又大又黄又爽视频免费| 少妇被粗大的猛进出69影院 | 大香蕉久久网| 99热全是精品| 亚洲av福利一区| 黄色一级大片看看| 国产黄色视频一区二区在线观看| av免费在线看不卡| 少妇被粗大的猛进出69影院 | 国产免费现黄频在线看| 免费黄色在线免费观看| 美女国产视频在线观看| 午夜福利,免费看| 在现免费观看毛片| 一区二区av电影网| 高清视频免费观看一区二区| 亚洲av免费高清在线观看| 少妇被粗大猛烈的视频| kizo精华| 夜夜骑夜夜射夜夜干|