基于柱坐標(biāo)系復(fù)頻移完全匹配層的隨鉆套管井聲場有限差分模擬

楊書博, 趙琪琪, 王 磊, 倪衛(wèi)寧, 李 新

(1.頁巖油氣富集機(jī)理與有效開發(fā)國家重點實驗室,北京 100101; 2.中國石化石油工程技術(shù)研究院,北京 100101;3.中國石油大學(xué)(北京)克拉瑪依校區(qū),新疆克拉瑪依 834000)

近年來,隨著國內(nèi)外油氣勘探開發(fā)程度的不斷提高,常規(guī)油氣藏的發(fā)現(xiàn)逐漸減少,深水深層及非常規(guī)油氣藏逐步成為勘探開發(fā)的主要對象,大斜度井和水平井的業(yè)務(wù)量也大幅度增加。在這些復(fù)雜井況下進(jìn)行電纜測井存在成本高、難度大、大量占用井架時間等問題,而隨鉆測井可在下放或上提儀器時進(jìn)行滑眼測量,既可以降低深水和深層環(huán)境中的測量時間和成本,又可以彌補(bǔ)大位移井和水平井中電纜測井爬行器無法深入測量的不足[1-5]。國內(nèi)外均已經(jīng)針對隨鉆固井質(zhì)量評價開展了一系列的研究,并在現(xiàn)場測試中取得了很好的應(yīng)用效果[5-10]。然而,由于隨鉆測井和電纜測井在測量環(huán)境中存在巨大差異,不能直接將電纜測井的結(jié)論直接應(yīng)用到隨鉆測井中。目前基于隨鉆聲波測井的固井質(zhì)量評價研究尚存在模式波傳播特征不清楚的情況,有必要開展相關(guān)數(shù)值模擬方面的工作。時域有限差分算法(finite difference time domain,FDTD)是現(xiàn)階段地球物理領(lǐng)域應(yīng)用最為廣泛的數(shù)值模擬方法[11-13]。為了消除來自人為截斷邊界的反射,通常需要在主研究區(qū)域外增加完全匹配層(perfectly matched layer,PML)。傳統(tǒng)的PML吸收邊界條件僅采用一個衰減參數(shù),在掠入射、高波阻抗差異以及各向異性的情況下,反射現(xiàn)象十分嚴(yán)重,甚至?xí)霈F(xiàn)不穩(wěn)定的現(xiàn)象[14-16]。針對這種情況,Roden和Gendney提出了更通用的復(fù)頻移完全匹配層(CFS-PML)吸收邊界條件[17]。CFS-PML吸收邊界條件采用3個參數(shù),通過將常規(guī)的復(fù)伸展坐標(biāo)系的實軸放大,并把虛軸的零點位置平移到負(fù)半軸,極大地提高了適用性[18-19]。然而,上述研究工作大多數(shù)是基于直角坐標(biāo)系開展,能否引入到柱坐標(biāo)系中尚需進(jìn)一步研究。隨鉆條件下固井質(zhì)量評價測井的數(shù)值模擬對于吸收邊界條件的要求很高,其原因主要包括以下兩個方面:第一,模型中層狀介質(zhì)的徑向長度與軸向長度之比很小。例如,流體環(huán)隙的厚度通常小于1 cm,而接收器源距通常大于1 m,即模型中流體環(huán)隙橫縱比通常小于1/100。第二,模型包含多個波阻抗差異很大的流固界面。例如,流體和鉆鋌或套管之間的波阻抗差異超過30倍。傳統(tǒng)的提高吸收邊界效果的方法是增加PML厚度,然而這會造成計算量急劇增大。筆者基于柱坐標(biāo)系彈性波一階速度-應(yīng)力方程組,提出一種適用于井孔聲場有限差分模擬的CFS-NPML吸收邊界條件。利用改進(jìn)的有限差分算法對比不同吸收邊界條件在隨鉆套管井聲場計算中的吸收效果,分析不同環(huán)隙位置和膠結(jié)質(zhì)量條件下充液井孔內(nèi)各模式波的傳播特征。

1 理論推導(dǎo)

基于彈性波一階速度-應(yīng)力方程組,推導(dǎo)CFS-NPML吸收邊界條件在二維柱坐標(biāo)系有限差分算法中的實現(xiàn)過程[16]。由于三維柱坐標(biāo)系和二維柱坐標(biāo)系中吸收邊界條件的實現(xiàn)過程完全相同,在此不再贅述。

1.1 柱坐標(biāo)系下CFS-NPML吸收邊界條件的推導(dǎo)

二維柱坐標(biāo)系下頻率域一階速度-應(yīng)力方程組滿足:

(1)

式中,Vr、Vθ和Vz為頻率域速度分量;Trr、Tθθ、Tzz、Trθ、Trz和Tθz為頻率域應(yīng)力分量;λ和μ為介質(zhì)的彈性系數(shù);m為聲源階數(shù),取0、1和2分別代表單極子、偶極子和四極子聲源。

(2)

(3)

sp(ω)為復(fù)頻移伸展函數(shù),對于CFS-NPML吸收邊界條件[19],其表達(dá)式為

(4)

式中,κp為比例因子,主要影響瞬逝波的吸收效果;αp為頻移因子,主要影響低頻成分的吸收效果;dp為衰減因子,主要影響波場在PML內(nèi)呈指數(shù)衰減的效果。故

(5)

將坐標(biāo)系變換回時間域可得

ζp(t).

(6)

式中,δ(t)為狄拉克函數(shù);H(t)為階躍函數(shù)。因此,將式(3)轉(zhuǎn)換到時間域可得

(7)

式中,*為卷積符號。在離散情況下,nΔt時刻的卷積可寫為

(8)

式中,n為時間索引;Δt為時間步長。利用梯形積分公式可得

(9)

其中

(10)

其中

(11)

其中

1.2 參數(shù)設(shè)置

理論上,對于具有相同彈性介質(zhì)參數(shù)而衰減參數(shù)不同的PML介質(zhì),波均將無反射地衰減傳播。但是,有限差分的數(shù)值離散會在不同PML層界面處引起數(shù)值反射。為了減小這種不期望的反射,必須使PML衰減參數(shù)平滑地變化[17,19]。因此本文中使用下式中的多項式函數(shù)來設(shè)置各吸收邊界參數(shù):

(12)

其中

式中,dp,max、κp,max和αp,max分別為dp、κp和αp的最大值;lp為到PML內(nèi)邊界的距離;lp,max為PML的厚度;nd、nκ和nα為多項式的階數(shù),通常取nd=2、nκ=2和nα=1;vf為流體的縱波速度;R為Dirichlet條件下平面縱波垂直入射到PML外層邊界時的理論反射系數(shù),通常取R=0.000 001;W0為該數(shù)值格式中每個波長所需要的最少采樣點數(shù),通常取W0=12;Δmax為空間網(wǎng)格的最大值;f0為聲源主頻。

上述CFS-NPML吸收邊界條件的各參數(shù)值為參考值,可以滿足絕大多數(shù)聲波測井?dāng)?shù)值模擬的需求。但是,在一些極特殊的情況下,仍然可以通過調(diào)整上述CFS-NPML吸收邊界條件的各參數(shù)值來獲得更好的吸收效果。

2 數(shù)值模擬分析

2.1 計算模型

計算模型如圖1所示。模型關(guān)于井軸對稱,尺寸為0.5 m×4 m。z軸與井軸重合,r軸與井軸垂直。模型從內(nèi)到外共有6層介質(zhì),即鉆鋌內(nèi)流體、鉆鋌、鉆鋌外流體、套管、水泥和地層,所有介質(zhì)層的界面都是以井軸為中心的圓柱面。聲源和接收器均位于鉆鋌外壁處的流體中,聲源中心點在z方向的坐標(biāo)為0.5 m,各接收器到聲源中心點的距離(源距)分別為2.0、2.2、2.4、2.6、2.8、3.0、3.2和3.4 m,分別記為R1、R2、R3、R4、R5、R6、R7和R8。計算模型中各層介質(zhì)的參數(shù)如表1所示。

圖1 計算模型示意圖Fig.1 Schematic diagram of calculation model

表1 計算模型參數(shù)

聲源所施加的激勵信號的峰值頻率為10 kHz;采用雷克子波作為聲源函數(shù);采用爆炸源作為單極子聲源的加載方式;采用完全匹配層邊界條件對在主研究區(qū)域邊界處的波場進(jìn)行吸收,采用并行算法對計算進(jìn)行加速[20]。r方向的空間步長為2.5 mm,z方向的空間步長均為5 mm;時間步長為0.2 μs。

2.2 吸收邊界條件對比

分別利用基于NPML和CFS-NPML吸收邊界條件的有限差分算法對隨鉆套管井模型進(jìn)行模擬。設(shè)置NPML參數(shù)滿足:lr,max=lz,max=20,dr,max=dz,max=621 700;CFS-NPML參數(shù)滿足:lr,max=lz,max=20,dr,max=dz,max=621 700,κr,max=κz,max=5,αr,max=αz,max=πf0。圖2為不同吸收邊界條件下各時刻的波場快照。由圖2可知,在1 ms之前,NPML和CFS-NPML吸收邊界條件的吸收效果都較好,吸收邊界處均未出現(xiàn)明顯的反射現(xiàn)象,計算過程十分穩(wěn)定;在2 ms時刻,NPML吸收邊界處出現(xiàn)不穩(wěn)定現(xiàn)象,而CFS-NPML吸收邊界處未出現(xiàn)不穩(wěn)定現(xiàn)象和明顯的反射現(xiàn)象;繼續(xù)計算至12 ms時刻,NPML吸收邊界處的計算噪聲持續(xù)傳播并污染整個計算區(qū)域,而CFS-NPML吸收邊界條件依舊具有良好的吸收效果。圖3為不同吸收邊界條件下2 m源距處的接收波形。由圖3可知,在計算穩(wěn)定的時刻(約2.5 ms之前),NPML吸收邊界條件和CFS-NPML吸收邊界條件的吸收效果基本一致,接收器均未接收到明顯的反射波;當(dāng)計算出現(xiàn)不穩(wěn)定之后,NPML吸收邊界條件所產(chǎn)生的計算噪聲的幅度持續(xù)增長,嚴(yán)重影響后續(xù)波形的計算。

圖2 不同吸收邊界條件下各時刻的波場快照Fig.2 Snapshots with different absorption boundary conditions at each time

圖3 不同吸收邊界條件下2 m源距處的接收波形Fig.3 Received waveforms with different absorption boundary conditions with the offset of 2 m

2.3 聲場響應(yīng)分析

為了分析環(huán)隙位置對接收波形的影響,分別計算完全膠結(jié)、Ⅰ界面膠結(jié)不好、Ⅱ界面膠結(jié)不好和完全不膠結(jié)4種情況下的隨鉆套管井聲場。計算模型如圖4所示,與完全膠結(jié)模型相比,Ⅰ界面膠結(jié)不好模型是在套管和水泥之間存在一個厚度1 cm的流體層,水泥的厚度也相應(yīng)地減小1 cm;Ⅱ界面膠結(jié)不好模型是在水泥和地層之間存在一個厚度1 cm的流體層,水泥的厚度也相應(yīng)地減小1 cm;完全不膠結(jié)模型是將水泥完全替代成流體。圖5～8分別為4種模型中的聲場響應(yīng)。將陣列波形數(shù)據(jù)通過時間-慢度相關(guān)法進(jìn)行處理,可以得到各模式波的時間-慢度相關(guān)(STC)圖[5]。觀察可知,4種模型的接收波形中均主要包含套管波(CA)、鉆鋌波(CO)、滑行橫波(S)和斯通利波(ST),滑行縱波被掩蓋。套管波和鉆鋌波的幅度小于滑行橫波和斯通利波的幅度。套管波到時最早,斯通利波到時最晚。由圖5可知,完全膠結(jié)模型中4種模式波的波速分別約為4 970、4 500、2 280和1 350 m/s。由圖6可知, Ⅰ 界面膠結(jié)不好模型中4種模式波的波速分別約為5 530、4 540、2 280和1 350 m/s。由圖7可知, Ⅱ 界面膠結(jié)不好模型中4種模式波的波速分別約為5 120、4 750、2 280和1 350 m/s。由圖8可知,完全不膠結(jié)模型中4種模式波的波速分別約為5 540、4 450、2 280和1 350 m/s。觀察可知,就套管波波速而言,完全不膠結(jié)模型約等于 Ⅰ 界面膠結(jié)不好模型,大于 Ⅱ 界面膠結(jié)不好模型,大于完全膠結(jié)模型,這與電纜測井的情況基本一致[5,21-25]。就鉆鋌波波速而言,4種模型有一定差異,但差異不大,即鉆鋌波受套管井環(huán)境影響相對較小。就滑行橫波和斯通利波而言,4種模型完全一樣。由于套管波和鉆鋌波波速相近,在時域上難以分離(特別是小源距情況下),因此共同構(gòu)成了首波。

圖4 不同環(huán)隙位置模型Fig.4 Calculation models with annuluses at different positions

圖5 完全膠結(jié)模型中的聲場響應(yīng)Fig.5 Responses in well-cemented cased-hole model

圖6 Ⅰ界面膠結(jié)不好模型中的聲場響應(yīng)Fig.6 Responses in cased-hole model with a weak casing-cement interface

圖7 Ⅱ界面膠結(jié)不好模型中的聲場響應(yīng)Fig.7 Responses in cased-hole model with a weak cement-formation interface

圖8 完全不膠結(jié)模型中的聲場響應(yīng)Fig.8 Responses in bad-cemented cased-hole model

對不同環(huán)隙位置模型中2 m源距處的首波波形進(jìn)行對比,結(jié)果如圖9所示。觀察可知,4種情況下首波幅度差異較大,從大到小依次為完全不膠結(jié)模型、Ⅰ 界面膠結(jié)不好模型、Ⅱ 界面膠結(jié)不好模型和完全膠結(jié)模型,這與電纜測井的情況也基本一致[5,21-25]。

圖9 不同環(huán)隙位置模型2 m源距處的接收波形對比Fig.9 Comparison of received waveforms in models with different annulus positions for the offset of 2 m

為了分析膠結(jié)質(zhì)量對接收波形的影響,保持水泥充滿套管與地層間的環(huán)空,改變水泥參數(shù)以模擬膠結(jié)質(zhì)量從0到1的變化,分別計算了不同水泥參數(shù)情況下的隨鉆套管井聲場。如表2所示,當(dāng)水泥參數(shù)分別滿足水泥100、水泥80、水泥60、水泥40、水泥20和水泥0時,代表完全膠結(jié)、80%膠結(jié)、60%膠結(jié)、40%膠結(jié)、20%膠結(jié)和完全不膠結(jié)的情況。各種膠結(jié)質(zhì)量情況下,水泥彈性參數(shù)的設(shè)置規(guī)則為

(13)

表2 水泥參數(shù)

式中,vp,i、vs,i和ρi分別為第i種膠結(jié)質(zhì)量情況下水泥的縱波波速、橫波波速和密度;Ii為第i種膠結(jié)質(zhì)量情況下的膠結(jié)指數(shù);vp,100%、vs,100%和ρ100%分別為完全膠結(jié)情況下水泥的縱波波速、橫波波速和密度;vp,0%、vs,0%和ρ0%分別為完全不膠結(jié)情況下水泥的縱波波速、橫波波速和密度。圖10為不同膠結(jié)質(zhì)量模型中2 m源距處的聲場響應(yīng)。通過統(tǒng)計不同情況下的首波峰峰值,并除以完全不膠結(jié)情況下的首波峰峰值,可以得到首波相對幅度隨膠結(jié)指數(shù)的變化關(guān)系。分析可知,當(dāng)膠結(jié)指數(shù)小于0.5時,隨著膠結(jié)指數(shù)的增大,首波幅度快速減小;當(dāng)膠結(jié)指數(shù)大于0.5時,隨著膠結(jié)指數(shù)的增大,首波幅度基本不變。造成這種現(xiàn)象的原因是,當(dāng)膠結(jié)指數(shù)較小時,套管波的幅度遠(yuǎn)遠(yuǎn)大于鉆鋌波的幅度,在首波中占主導(dǎo)地位,首波變化規(guī)律與電纜測井基本相同[5,21-25];當(dāng)膠結(jié)指數(shù)較大時,套管波的幅度與鉆鋌波的幅度相近甚至更小,首波變化規(guī)律與電纜測井不再相同[5,21-25],這也說明鉆鋌波的幅度與膠結(jié)指數(shù)關(guān)系不大。

圖10 不同膠結(jié)質(zhì)量模型2 m源距處的聲場響應(yīng)對比Fig.10 Comparison of received waveforms in models with different bond indexes for the offset of 2 m

3 結(jié) 論

(1)與NPML吸收邊界條件相比,CFS-NPML吸收邊界條件在隨鉆套管井聲場有限差分模擬中具有更好的吸收效果和更高的計算穩(wěn)定性,且CFS-NPML吸收邊界條件對數(shù)值模擬的計算量、存儲量和編程復(fù)雜度沒有太大影響。

(2)當(dāng)環(huán)隙位置不同時,套管波在隨鉆條件下的變化規(guī)律與電纜條件下的變化規(guī)律基本一致,而鉆鋌波受套管井環(huán)境的影響相對較小。套管波波速和鉆鋌波波速相近,短源距條件下時域波形會發(fā)生重疊,共同構(gòu)成首波。

(3)當(dāng)膠結(jié)指數(shù)小于0.5時,隨著膠結(jié)指數(shù)的增大,首波幅度快速減小;當(dāng)膠結(jié)指數(shù)大于0.5時,隨著膠結(jié)指數(shù)的增大,首波幅度基本不變。該特點使得常規(guī)電纜聲波固井質(zhì)量評價方法不適用于隨鉆測井環(huán)境,需要開展新方法的研究。