應(yīng)用邊界元法分析頁巖地層井眼坍塌問題

馬天壽，陳平

(西南石油大學 油氣藏地質(zhì)及開發(fā)工程國家重點實驗室，四川 成都，610500)

應(yīng)用邊界元法分析頁巖地層井眼坍塌問題

馬天壽，陳平

(西南石油大學 油氣藏地質(zhì)及開發(fā)工程國家重點實驗室，四川 成都，610500)

基于各向異性連續(xù)介質(zhì)力學和邊界元理論，建立橫觀各向同性頁巖地層井眼坍塌定問題的基本微分方程，推導出基本方程的邊界積分方程及其離散解，并得出井周應(yīng)力和位移的邊界元離散解，結(jié)合Mohr?Coulomb準則判別井眼穩(wěn)定性，形成基于邊界元法(BEM)求解井眼坍塌問題的基本方法。建立頁巖地層井眼坍塌問題的物理模型，采用各向同性地層模型對BEM進行檢驗，并開展井周應(yīng)力分布規(guī)律研究和應(yīng)用實例分析。研究結(jié)果表明：BEM求解結(jié)果與解析解吻合較好，二者相對誤差小于2.49%；彈性模量各向異性、水平地應(yīng)力差異和鉆井液密度等對井壁應(yīng)力分布影響較大(其中水平地應(yīng)力差異的影響最大)，而泊松比各向異性的影響較?。籛201井1 530 m井段井眼穩(wěn)定性分析結(jié)果與電測結(jié)果吻合良好，BEM計算的擴徑率約為9.0%，而電測擴徑率約為10.0%。

頁巖；井眼穩(wěn)定；橫觀各向同性；邊界元(BEM)；數(shù)值方法

頁巖氣是以富有機質(zhì)頁巖為氣源巖、儲層或蓋層，在頁巖地層中不間斷供氣、連續(xù)聚集而形成的一種非常規(guī)天然氣，存在于幾乎所有的盆地中，由于埋藏深度、含氣飽和度等差異較大，分別具有不同的工業(yè)價值[1?2]。通過多年的努力，美國取得頁巖氣勘探開發(fā)技術(shù)突破，近10 a產(chǎn)量快速增長，美國已超過俄羅斯成為全球最大的天然氣生產(chǎn)國，由天然氣凈進口國轉(zhuǎn)變?yōu)閮舫隹趪?。受美國頁巖氣成功開發(fā)的啟示，我國頁巖氣開發(fā)也已成為了業(yè)界的焦點問題，并開展了頁巖氣勘探評價和先導試驗[1]。在中石油在四川南部的長寧—威遠國家級頁巖氣示范區(qū)、殼牌和中石油合作的云南昭通國家級頁巖氣示范區(qū)以及中石化在重慶東部礁石壩構(gòu)造的鉆探試采結(jié)果表明，頁巖地層鉆井過程中井眼垮塌較嚴重，導致卡鉆、埋鉆等井下復雜和事故頻頻發(fā)生[1?2]。減少井下復雜和事故發(fā)生的關(guān)鍵在于提高頁巖地層井眼的穩(wěn)定性，這使頁巖井眼穩(wěn)定問題成為頁巖氣高效開發(fā)的關(guān)鍵問題之一[3]。針對頁巖地層井眼穩(wěn)定問題的研究，國內(nèi)外已經(jīng)開展了大量研究并取得了較豐富的成果[4]，已經(jīng)形成多種井眼穩(wěn)定的計算和分析方法，如解析解方法、有限元方法(FEM)、離散元方法(DEM)、其他數(shù)值計算方法等。解析解方法為最早出現(xiàn)的井眼穩(wěn)定分析方法之一，WESTERGARD[5]把柱坐標引入井眼坐標系中，利用應(yīng)力函數(shù)計算井眼彈塑性問題，初步研究了井壁穩(wěn)定機理，并對現(xiàn)場實例和處理方法進行了討論，直接推動井眼穩(wěn)定問題進入定量化研究階段；此后逐漸形成了一系列以線彈性力學、孔隙彈性力學、彈塑性力學為基礎(chǔ)的井壁穩(wěn)定力學計算方法，由于這些方法計算簡便，且可得出沿井深的井眼坍塌、破裂壓力剖面，因而被廣泛采用。然而，解析解方法假設(shè)條件過多，計算結(jié)果往往存在較大的誤差。隨著有限元理論的快速發(fā)展，F(xiàn)EM開始被用于分析井眼穩(wěn)定問題，F(xiàn)EM主要被用于井眼穩(wěn)定彈塑性分析和多場耦合分析；考慮到井眼穩(wěn)定問題的對稱性，大量學者將其簡化為廣義平面應(yīng)變問題進行求解[6?10]，也有一些學者直接采用三維有限元模型求解分析[11?12]；FEM計算靈活性較高、可處理應(yīng)力集中和非線性問題；但FEM經(jīng)常會出現(xiàn)網(wǎng)格劃分困難、矩陣運算量過大、計算內(nèi)存大、計算較慢等問題。DEM是一種非連續(xù)介質(zhì)力學分析方法[13?15]，因此，DEM主要適用于分析節(jié)理、破碎性地層的井眼穩(wěn)定問題。有限差分方法(FDM)[16?18]計算簡單，但不適用于不規(guī)則物體，也不能處理應(yīng)力集中問題。真實破裂過程分析方法(RFPA)[19]為一種基于統(tǒng)計損傷力學的方法，在工程應(yīng)用中其基礎(chǔ)參數(shù)獲取困難。不難看出，這些方法仍存在一些難以解決的問題。WANG等[20?21]采用邊界元方法(BEM)研究了井漏中的裂縫延伸問題，指出BEM對邊界積分方程離散求解的現(xiàn)代數(shù)值分析方法，屬于近似求解方法之一，BEM適合于求解無窮域或應(yīng)力集中問題，具有降階、降維的特點[22?23]，其優(yōu)點在于依靠微分方程階數(shù)的減少以降低對近似試函數(shù)所必需的可導階數(shù)[23]。因此，BEM方法所得控制方程為積分方程，該方法降低了計算維度，主要考慮物體邊界(邊界單元處理容易)，其計算精度較高、內(nèi)存小、速度快，目前采用BEM方法研究井眼坍塌問題的文獻鮮見報道。因此，本文作者引入BEM方法來分析井眼坍塌問題，通過分析橫觀各向同性頁巖地層的BEM求解方法和求解結(jié)果，闡述BEM在井眼穩(wěn)定問題分析中的應(yīng)用，借此將BEM引入到井眼穩(wěn)定分析方法中，從而為井眼穩(wěn)定問題的分析提供更多有力的手段和方法。

1　頁巖井眼穩(wěn)定問題分析方法

1.1數(shù)學模型

頁巖的典型特點是層理、微裂隙發(fā)育，層理性地層是指有一組近于平行的層理面，層理的強度明顯比巖石本體強度低很多，CHENEVERT等[24]通過實驗證實當層理面法線與最大主應(yīng)力夾角為20o～30o時，其強度比垂直于層理面方向取心降低40%左右。層理性頁巖與橫觀各向同性材料具有極強的相似性，在層理面內(nèi)各方向的性質(zhì)基本相同，而在層理面法線方向性質(zhì)不同，因此，可以將頁巖地層視為橫觀各向同性材料，并給出相應(yīng)的數(shù)學模型。

平衡方程為

幾何方程為

橫觀各向同性頁巖的本構(gòu)方程[23]

式中：σij為總應(yīng)力；xj為坐標軸；fi為體力；εij為應(yīng)變分量；ui為位移；sijkl為柔度矩陣；i，j，k和l為循環(huán)角標。

式中：E和E′為平行和垂直橫觀各向同性面的彈性模量；v和v′為平行和垂直橫觀各向同性面的泊松比。

若假設(shè)井眼軸線方向為z軸，則井眼平面橫截面面積和形狀沿z軸不發(fā)生改變，作用的外力與z軸垂直，則井周地層的位移將主要發(fā)生在橫截面內(nèi)，因此，往往將其簡化為廣義的平面應(yīng)變問題。于是，在Oxyz直角坐標系中所有的應(yīng)變分量均與z方向無關(guān)[25]，即

將式(3)代入式(6)可得

再將式(7)代入式(3)，可得平面應(yīng)變問題本構(gòu)方程

式中：σxx，σyy，σzz，τxy，τyz和τxz為應(yīng)力分量；εxx，εyy，εzz，γxy，γyz和γxz為應(yīng)變分量。

在Γ1上，靜力邊界條件為

在Γ2上，位移邊界條件為

式中：pi為邊界應(yīng)力；為邊界上給定的應(yīng)力分量，上方記號“”表示該值已知，下同；ui為邊界位移；為邊界上給定的位移分量；Γ1和Γ2為邊界。

將式(2)代入式(8)中，再將得到的結(jié)果代入式(1)中，即可得到由位移表示的應(yīng)力平衡方程，即本文求解的基本方程。

1.2邊界元求解方法

為了求解橫觀各向同性頁巖地層井眼穩(wěn)定問題的基本方程，首先要將基本方程轉(zhuǎn)化為邊界積分方程。但橫觀各向同性材料基本解的數(shù)學描述較復雜，為此，采用一種適用于任一各向異性材料的方法，即參照各向同性材料基本解并以迭代法逼近正確解[22]。若忽略體力，則基本微分方程的加權(quán)余量表達式為

若考慮采用參照彈性模型的各向同性材料基本解，此參照彈性模型是以各向異性材料常數(shù)的平均值建立的。因此，原有彈性常數(shù)張量dijkl[22]表示為

式中：dijkl為實際彈性常數(shù)；為參照態(tài)彈性常數(shù)；為實際與參照態(tài)彈性常數(shù)的差值。因此，本構(gòu)方程可表示為

將式(11)進行分部積分得

將式(13)代入式(14)，并將式(15)帶入，進行分部積分后可得

由于基本解同時滿足域內(nèi)和邊界，將式(16)針對邊界上的力作用點m，可改寫為

式(17)中不僅出現(xiàn)了邊界積分項，還出現(xiàn)了域積分項，可將域積分項變換成體力項和邊界積分項[22]再計算。對式(17)右端第2項采用分部積分可得

式(18)中右端第一項可等效處理為體力項，即

將式(19)代入式(17)可得

將平面問題位移、表面力、等效面力和等效體力矢量可表示為

為了將邊界積分方程離散化，采用常單元分割求解域及邊界，則u和p在單元內(nèi)是常量，等于單元中點(節(jié)點)的值。則式(20)可寫成離散形式，即

式中：ui為基本解多用點i(邊界或域內(nèi))的位移；ci為2×2階的常數(shù)矩陣，取決于該節(jié)點i所在的位置及幾何形狀；為與i節(jié)點和j節(jié)點有關(guān)的邊界積分項，稱為影響系數(shù)[22]，經(jīng)過積分計算后成為2×2階的矩陣和Gij。

而等效體力項的數(shù)值積分需按下式計算

式中：wq為權(quán)函數(shù)；Ωk為域元面積。數(shù)值積分計算后得到Bik的2個分量。函數(shù)在積分節(jié)點q上計算，q的取值為1～t。

據(jù)此，將式(24)代入式(23)可得

若令

則式(25)可寫成

對邊界上所有節(jié)點應(yīng)用式(27)，可得矩陣方程

式中：H和G均為2N×2N階的矩陣，N為邊界上的節(jié)點總數(shù)；B為體力項列向量。

將已知邊界條件代入式(28)，并將未知量及其系數(shù)移至方程的左邊，其余項移至方程的右邊可得

式中：A為與H和G相關(guān)的矩陣；x為未知數(shù)組，包括位移和面力；F為已知邊界條件和對應(yīng)H和G相乘所得的列向量。

求解式(29)可以獲得全部邊值。當獲得全部邊界值(位移和表面力)后，便可計算域內(nèi)任一點的位移和應(yīng)力。為了得到求解域內(nèi)的位移分布，可采用Somigliana公式求解[22,25]，將Somigliana公式表示的邊界積分方程寫成離散形式：

為了得到求解域內(nèi)點的應(yīng)力分布，可通過對該點位移的求導和本構(gòu)方程加以確定，求得域內(nèi)點應(yīng)力分布的邊界積分方程離散形式[22,25]為

其中：Dkl和Skl為系數(shù)矩陣；Dkij和Skij為與邊界條件有關(guān)的系數(shù)；δij為克羅內(nèi)克函數(shù)；r為矢徑長度。

式(31)即為采用BEM方法求解井眼周圍應(yīng)力分布的方法，通過BEM計算得出井周應(yīng)力分布后，結(jié)合強度判別準則即可進行井眼穩(wěn)定性分析，本文采用Mohr?Coulomb準則進行井眼穩(wěn)定性分析。

1.3強度判別準則

在進行井眼穩(wěn)定性判別前，需要將廣義平面二向應(yīng)力狀態(tài)結(jié)果轉(zhuǎn)換為空間主應(yīng)力狀態(tài)，即計算出空間應(yīng)力狀態(tài)下的主應(yīng)力。采用BEM求解得到的應(yīng)力張量σ為

因此，井眼局部圓柱坐標系下空間主應(yīng)力為：

式中：1σ和3σ分別井周地層的最大和最小主應(yīng)力；γ為井壁最大主應(yīng)力與井眼軸線之間的夾角。

獲得層理性頁巖井周最大、最小主應(yīng)力后，代入強度準則(Mohr?Coulomb準則)，即可進行穩(wěn)定性判別，Mohr?Coulomb準則可表示為

當采用Mohr?Coulomb準則判別井眼不穩(wěn)定時，假設(shè)頁巖地層強度系數(shù)(strength factor，F(xiàn)S)為

當FS＞1時，井周頁巖地層穩(wěn)定；當FS＜1時，井周頁巖地層將發(fā)生失穩(wěn)。因此，F(xiàn)S越小，則表示井眼失穩(wěn)的風險越大。

若將最大主應(yīng)力σ1和最小主應(yīng)力σ3代入式(38)，可得Mohr?Coulomb準則的另一種表達形式為

式中：τ為剪切面上剪應(yīng)力；σn為剪切面上正應(yīng)力；c0為巖石內(nèi)聚力；?為巖石內(nèi)摩擦角；pp為地層孔隙壓力；α為Biot系數(shù)。

2　邊界元求解模型的建立與驗證

2.1求解模型的建立

對于頁巖地層中所鉆的井眼來說，沿軸向方向的應(yīng)變可近似認為是一常量，即可將復雜的三維空間應(yīng)力分析問題轉(zhuǎn)化為廣義平面應(yīng)變問題來考慮。因此，當采用BEM求解時，所建立的井周地層計算網(wǎng)絡(luò)模型如圖1所示，分別為求解域計算網(wǎng)絡(luò)模型與井壁周圍離散后的計算網(wǎng)絡(luò)模型局部視圖，本文后續(xù)分析均采用這種計算網(wǎng)絡(luò)模型，圖1(b)中的a→b→c→d→a為應(yīng)力分布取值路徑。在后續(xù)計算和分析中，所采用的基礎(chǔ)參數(shù)如表1所示。

圖1　井周地層計算網(wǎng)絡(luò)模型Fig. 1 Mesh model of formation around borehole

表1　井周地層計算基礎(chǔ)參數(shù)Table 1 Basic parameters of formation around borehole

2.2模型驗證

為了檢驗本文模型的準確性，取各向同性均質(zhì)地層直井井眼模型來檢驗BEM的準確性，分別采用常規(guī)井壁穩(wěn)定力學分析方法(解析解)[4]和BEM計算井周應(yīng)力分布情況。此時，取各向異性系數(shù)為1.0(即nE=E′/E=1.0，nυ=υ′/υ=1.0)，其他基礎(chǔ)參數(shù)如表1所示。首先，采用BEM計算出各向同性均質(zhì)地層井周應(yīng)力分布情況，提取出一倍井眼范圍內(nèi)的井周應(yīng)力分布，結(jié)果如圖2所示。為了對比解析解方法和BEM計算結(jié)果，取井壁處a→b→c→d→a路徑的應(yīng)力(徑向應(yīng)力rσ、切向應(yīng)力θσ、軸向應(yīng)力zσ)計算結(jié)果進行對比，如圖3所示。

圖2　邊界元方法解出的井周應(yīng)力分布圖Fig. 2 Stress distribution of wellbore by BEM

圖3　井壁(r=R)應(yīng)力分布圖Fig. 3 Stress distribution of borehole wall (r=R)

由圖3可以看出，井周應(yīng)力分布的求解結(jié)果在徑向、切向和軸向(z軸)3個方向的計算結(jié)果均吻合較好，BEM計算結(jié)果與解析解計算結(jié)果之間的相對誤差小于2.49%。另外，井周應(yīng)力分布在最小水平地應(yīng)力方向上的應(yīng)力差異最大(井壁切向應(yīng)力遠高于徑向應(yīng)力)，在最大水平地應(yīng)力方向上的差異較小(井壁切向、軸向和徑向應(yīng)力十分接近)，使各向同性地層中井眼坍塌通常發(fā)生在最小水平地應(yīng)力方向上。對于二者計算誤差，主要原因是由于BEM是一種近似求解的數(shù)值方法，加之劃分網(wǎng)格單元的不規(guī)則性，使其計算結(jié)果與解析解計算結(jié)果有些差異，但總體上誤差較小，能夠滿足工程所需精度，驗證了BEM的正確性和準確性。

3　計算結(jié)果討論

3.1彈性模量對井壁應(yīng)力的影響

定義橫觀各向同性頁巖地層彈性模量各向異性系數(shù)nE=E′/E，則nE表征了頁巖地層巖石的各向異性程度，其值越小表示各向異性程度越高，nE=1表示各向同性，通過計算得到圖4所示的結(jié)果。

由圖4可見：頁巖地層井壁切向應(yīng)力分布規(guī)律近似“啞鈴”狀，nE對井壁切向應(yīng)力的影響比較顯著；各向異性程度越高(nE越小)，對井壁切向應(yīng)力分布的影響越大，井壁切向應(yīng)力分布由規(guī)則“啞鈴”狀逐漸變成不規(guī)則形狀；在不同各向異性程度的地層中，井壁切向應(yīng)力的極大值和極小值不同，各向異性程度越高則井壁切向應(yīng)力極值越高(如nE=0.25時切向應(yīng)力極值為97.68 MPa，而各向同性地層為94.32 MPa)，且切向應(yīng)力極小值也越小；井壁切向應(yīng)力分布規(guī)律還與井周角和地應(yīng)力方位有關(guān)，井壁切向應(yīng)力極大值點在最小水平地應(yīng)力方向附近，隨著各向異性程度增加，應(yīng)力極大值點開始向最大水平地應(yīng)力方向偏轉(zhuǎn)。

圖4　nE對井壁切向應(yīng)力分布的影響Fig. 4 Influence of nEon tangential stress distribution at borehole wall

3.2泊松比對井壁應(yīng)力的影響

定義橫觀各向同性頁巖地層泊松比各向異性系數(shù)nυ=υ′/υ，則nυ為頁巖地層巖石的各向異性程度，其值越小表示各向異性程度越高，nυ=1為各向同性，通過計算得圖5所示的結(jié)果。由圖5可以看出：頁巖地層井壁切向應(yīng)力分布規(guī)律近似“啞鈴”狀，nυ對井壁切向應(yīng)力幾乎沒有影響；井壁切向應(yīng)力極大值約為94.32 MPa，其位于最小水平地應(yīng)力方向；井壁切向應(yīng)力極小值約為22.43 MPa，其位于最大水平地應(yīng)力方向。

圖5　nυ對井壁切向應(yīng)力分布的影響Fig. 5 Influence of nυon tangential stress distribution at borehole wall

3.3地應(yīng)力對井壁應(yīng)力的影響

定義頁巖地層水平地應(yīng)力比mσ=σh/σH，為此，可固定最大水平地應(yīng)力σH，取水平地應(yīng)力比mσ=0.25，0.50，0.75，1.00共4種情況分別計算，以研究地應(yīng)力對井壁應(yīng)力分布的影響?？紤]到泊松比各向異性對井壁應(yīng)力分布的影響較小，彈性模量各向異性對井壁應(yīng)力分布的影響較大，此處分別計算了在nE為0.25和1.00這2種情況下，不同地應(yīng)力比下的井壁應(yīng)力分布情況，其結(jié)果如圖6所示。

圖6　水平地應(yīng)力比mσ對井壁切向應(yīng)力分布的影響Fig. 6 Influence of ratio of horizontal in-situ stress (mσ) on tangential stress distribution at borehole wall

由圖6可以看出：在均勻地應(yīng)力下(mσ=1.00)，井壁切向應(yīng)力分布形狀為圓形(即均勻分布)，在非均勻地應(yīng)力下，井壁切向應(yīng)力分布形狀由圓形逐漸向“啞鈴”狀變化(即井壁切向應(yīng)力差異增加)，若水平地應(yīng)力比mσ進一步增大，井壁切向應(yīng)力分布的“啞鈴”狀更加明顯，這是導致井眼失穩(wěn)的主要原因之一；地應(yīng)力和彈性模量各向異性系數(shù)nE對井壁切向應(yīng)力均有較大影響，但地應(yīng)力的影響更為顯著；水平地應(yīng)力差異越大(mσ越小)，井壁切向應(yīng)力越大，均質(zhì)地層(nE=1.00)中mσ=0.25與均勻地應(yīng)力下的井壁切向應(yīng)力差異為36.13 MPa，橫觀各向同性地層(nE=0.25)中mσ=0.25與均勻地應(yīng)力下的井壁切向應(yīng)力差異為42.88 MPa，可見由于地應(yīng)力差異造成的井壁切向應(yīng)力差異達36.13 MPa，而由于地層彈性模量各向異性造成的井壁切向應(yīng)力差異僅為6.75 MPa，地應(yīng)力差異的影響比地層各向的影響更為顯著。

3.4井筒內(nèi)鉆井液密度對井壁應(yīng)力的影響

鉆開地層形成井眼后，鉆井液液柱壓力取代了所鉆巖層提供的支撐，破壞了地層原有的應(yīng)力平衡，引起井眼周圍巖石的應(yīng)力重新分布。如果重新分布的應(yīng)力超過巖石所能承受的最大載荷，將會導致井眼失穩(wěn)。為此，分析了井筒內(nèi)鉆井液密度(即井筒內(nèi)壓力)對井壁應(yīng)力的影響，分析結(jié)果如圖7所示。

由圖7可見：無論是均質(zhì)各向同性地層，還是橫觀各向同性地層，井筒內(nèi)鉆井液密度對井壁切向應(yīng)力的分布均有一定的影響，鉆井液密度增加，將導致井壁切向應(yīng)力降低；同時，隨著井內(nèi)鉆井液密度的增加，徑向應(yīng)力也顯著增加(井壁處徑向應(yīng)力等于井內(nèi)鉆井液液柱壓力)；在切向應(yīng)力降低、徑向應(yīng)力增加的綜合作用下，井周分布的應(yīng)力差異降低，這對維持井周地層穩(wěn)定(防止剪切破壞)具有重要作用，因此，增加鉆井液密度有利于維持井眼穩(wěn)定性。

4　應(yīng)用實例

以四川南部長寧—威遠頁巖氣示范區(qū)的一口評價井(W201井)為例進行分析，該井1 530 m井段為龍馬溪組頁巖地層，該地層中層理、頁理發(fā)育，產(chǎn)狀約為60°∠5°，井斜角為3.3°、井斜方位為313.1°，說明井眼與層理相交關(guān)系近似垂直關(guān)系。因此，采用橫觀各向同性BEM方法進行計算和分析較合理，基礎(chǔ)參數(shù)如表1所示，其中，彈性模量各向異性系數(shù)nE=0.82，泊松比各向異性系數(shù)nυ=0.85，計算得到的結(jié)果分別如圖8和圖9所示。

圖7　鉆井液密度對井壁切向應(yīng)力分布的影響Fig. 7 Influence of drilling fluid density on tangential stress distribution at borehole wall

由圖8和圖9可以看出：

1) 若以各向同性介質(zhì)分析W201井龍馬溪組頁巖地層井眼穩(wěn)定問題，采用密度為1.33 g/cm3的鉆井液，能夠維持井眼穩(wěn)定。井周地層強度系數(shù)大于1.03，強度系數(shù)最低的位置在最小水平地應(yīng)力方向上，即沿圖1所示路徑a→b之間井周角48°位置，分析結(jié)果與前述解析解分析結(jié)果一致，進一步驗證了方法的準確性和正確性。

2) 若以橫觀各向同性介質(zhì)分析W201井龍馬溪組頁巖地層井眼穩(wěn)定問題，采用密度為1.33 g/cm3的鉆井液，不能夠維持井眼穩(wěn)定，如圖8(b)和圖8(c)所示的強度系數(shù)分布情況。井周地層強度系數(shù)最低為0.83，根據(jù)強度系數(shù)定義，強度系數(shù)越低，井眼越不穩(wěn)定，強度系數(shù)最低的位置仍然在最小水平地應(yīng)力方向上，即路徑a→b之間井周角48°的位置。為此，通過提取計算結(jié)果，提取出強度系數(shù)小于1.0的區(qū)域，并繪制成井眼不穩(wěn)定區(qū)域圖(見圖9)。由圖9可以看出：井眼不穩(wěn)定區(qū)域不止發(fā)生在井壁上，而是發(fā)生在進入地層一定深度處，不安全區(qū)域進入地層的深度約為1.09倍井眼半徑處，即井眼擴徑失穩(wěn)后將形成9.0%的擴徑率，而不穩(wěn)定區(qū)域的井周夾角約為115°。

3) W201井1 530 m井段在實鉆過程中，采用密度為1.28 g/cm3的聚磺體系鉆井液，完鉆電測井徑表明該井段井眼擴徑率為10%左右。采用各向同性介質(zhì)計算結(jié)果表明，1.28 g/cm3的聚磺體系鉆井液能夠維持井眼穩(wěn)定；而采用橫觀各向同性介質(zhì)計算結(jié)果表明，1.28 g/cm3的聚磺體系鉆井液鉆進將發(fā)生約10%的井眼擴徑率。說明本文采用的BEM能夠滿足井眼穩(wěn)定問題分析的需求，進一步檢驗了BEM方法的精度。

圖8　井周地層強度系數(shù)分布圖Fig. 8 Strength factor distribution around borehole

圖9　橫觀各向同性頁巖井眼不穩(wěn)定區(qū)域Fig. 9 Instability area around borehole in transverse isotropic shale formation

5　結(jié)論

1) 基于各向異性連續(xù)介質(zhì)力學和邊界元理論，建立了橫觀各向同性頁巖地層的井眼穩(wěn)定問題基本微分方程，給出了基本方程的邊界積分方程及其離散形式解以及井周應(yīng)力和位移的邊界元離散形式解，結(jié)合Mohr?Coulomb準則判別井眼穩(wěn)定性，形成了基于BEM求解頁巖井眼穩(wěn)定問題的基本方法。

2) 建立了頁巖地層井眼穩(wěn)定問題的BEM求解模型，并以各向同性均質(zhì)地層為例檢驗?zāi)Ｐ?，BEM求解結(jié)果與解析解吻合較好，二者的相對誤差小于2.49%，驗證了BEM的正確性與準確性。

3) 在橫觀各向同性頁巖地層中，彈性模量各向異性系數(shù)nE、水平地應(yīng)力比m和井筒內(nèi)鉆井液密度等對井壁應(yīng)力分布產(chǎn)生顯著的影響，而泊松比各向異性系數(shù)nυ幾乎沒有影響。

4) W201井1530井段實例分析表明，采用橫觀各向同性頁巖地層計算結(jié)果顯示井眼失穩(wěn)擴徑率約9.0%，以各向同性均質(zhì)地層計算結(jié)果顯示該井段井眼穩(wěn)定性較好，而實鉆電測資料顯示該段井眼擴徑率為10.0%左右，說明采用BEM方法并考慮頁巖地層的各向異性能夠更加準確的判斷井眼穩(wěn)定情況，也進一步驗證了本文方法的精度。

[1] 郭彤樓, 張漢榮. 四川盆地焦石壩頁巖氣田形成與富集高產(chǎn)模式[J]. 石油勘探與開發(fā), 2014, 41(1): 28?36. GUO Tonglou, ZHANG Hanrong. Formation and enrichment mode of Jiaoshiba shale gas field, Sichuan Basin[J]. Petroleum Exploration and Development, 2014, 41(1): 28?36.

[2] 馬天壽, 陳平. 基于CT掃描技術(shù)研究頁巖水化細觀損傷特性[J]. 石油勘探與開發(fā), 2014, 41(2): 227?233. MA Tianshou, CHEN Ping. Studies of meso-damage characteristics of shale hydration based on CT scanning technology[J]. Petroleum Exploration and Development, 2014, 41(2): 227?233.

[3] 王紅巖, 劉玉章, 董大忠, 等. 中國南方海相頁巖氣高效開發(fā)的科學問題[J]. 石油勘探與開發(fā), 2013, 40(5): 574?579. WANG Hongyan, LIU Yuzhang, DONG Dazhong, et al. Scientific issues on effective development of marine shale gas in southern China[J]. Petroleum Exploration and Development, 2013, 40(5): 574?579.

[4] 陳勉, 金衍, 張廣清. 石油工程巖石力學[M]. 北京：科學出版社, 2008: 102?166. CHEN Mian, JIN Yan, ZHANG Guangqing. Rock mechanics on petroleum engineering[M]. Beijing: Science Press, 2008: 102?166.

[5] WESTERGARD H M. Plastic state of stress around a deep well[J]. Journal of Boston Society Civil Engineers, 1940, 27: 387?391.

[6] ROSHAN H, RAHMAN S S. Effects of ion advection and thermal convection on pore pressure changes in high permeable chemically active shale formations[J]. Petroleum Science and Technology, 2013, 31(7): 727?737.

[7] 盛金昌, 劉繼山, 許孝臣, 等. 頁巖鉆孔過程中的流固熱化學耦合模型[J]. 工程力學, 2009, 26(12): 240?245. SHENG Jinchang, LIU Jishan, XU Xiaochen, et al. A coupled porochemothermoelastic model for a borehole in shales[J]. Engineering Mechanics, 2009, 26(12): 240?245.

[8] 王倩, 周英操, 王剛, 等. 泥頁巖井壁穩(wěn)定流固化耦合模型[J].石油勘探與開發(fā), 2012, 39(4): 475?480. WANG Qian, ZHOU Yingcao, WANG Gang, et al. A fluid-solid-chemistry coupling model for shale wellbore stability[J]. Petroleum Exploration and Development, 2012, 39(4): 475?480.

[9] 王倩, 周英操, 唐玉林, 等.泥頁巖井壁穩(wěn)定影響因素分析[J].巖石力學與工程學報, 2012, 31(1): 171?179. WANG Qian, ZHOU Yingcao, TANG Yulin, et al. Analysis of effect factor in shale wellbore stability[J]. Chinese Journal of Rock Mechanics and Engineering, 2012, 31(1): 171?179.

[10] 賈善坡, 鄒臣頌, 王越之, 等. 基于熱?流?固耦合模型的石油鉆井施工過程數(shù)值分析[J]. 巖土力學, 2012, 33(S2): 321?328. JIA Shanpo, ZOU Chensong, WANG Yuezhi, et al. Numerical analysis of construction process of petroleum drilling based on thermal?hydro?mechanical coupling[J]. Rock and Soil Mechanics, 2012, 33(S2): 321?328.

[11] 李軍, 陳勉, 金衍, 等. 大位移井井壁穩(wěn)定三維彈塑性有限元分析[J]. 巖石力學與工程學報, 2004, 23(14): 2385?2389. LI Jun, CHEN Mian, JIN Yan, et al. Three-dimensional elastoplastic FEM analysis on borehole stability[J]. Chinese Journal of Rock Mechanics and Engineering, 2004, 23(14): 2385?2389.

[12] ROSHAN H, RAHMAN S S. A fully coupled chemo-poroelastic analysis of pore pressure and stress distribution around a wellbore in water active rocks[J]. Rock Mechanics and Rock Engineering, 2011, 44(2): 199?210.

[13] 李嗣貴, 鄧金根, 李明志. 節(jié)理破碎地層井壁穩(wěn)定的離散元分析[J]. 巖石力學與工程學報, 2002, 21(S1): 2139?2143. LI Sigui, DENG Jingen, LI Mingzhi. Discrete element analysis on wellbore stability in jointed broken formation[J]. Chinese Journal of Rock Mechanics and Engineering, 2002, 21(S1): 2139?2143.

[14] 屈平, 申瑞臣, 付利, 等. 三維離散元在煤層水平井井壁穩(wěn)定中的應(yīng)用[J]. 石油學報, 2011, 32(1): 153?157. QU Ping, SHEN Ruichen, FU Li, et al. Application of the 3D discrete element method in the wellbore stability of coal-bed horizontal wells[J]. Acta Petrolei Sinica, 2011, 32(1): 153?157.

[15] 尹虎, 李黔, 郭亮, 等. 離散單元法評價煤層井壁穩(wěn)定性[J].天然氣工業(yè), 2012, 32(11): 59?63. YIN Hu, LI Qian, GUO Liang, et al. Evaluation of wellbore stability in coal seams by use of the universal discrete element code (UDEC) method[J]. Natural Gas Industry, 2012, 32(11): 59?63.

[16] 鄧金根, 郭東旭, 周建良, 等. 泥頁巖井壁應(yīng)力的力學?化學耦合計算模式及數(shù)值求解方法[J]. 巖石力學與工程學報, 2003, 22(S1): 2250?2253. DENG Jingen, GUO Dongxu, ZHOU Jianliang, et al. Mechanics?chemistry coupling calculation model of borehole stress in shale formation and its numerical solving method[J]. Chinese Journal of Rock Mechanics and Engineering, 2003, 22(S1): 2250?2253.

[17] 程遠方, 張鋒, 王京印, 等. 泥頁巖井壁坍塌周期分析[J]. 中國石油大學學報(自然科學版), 2007, 31(1): 63?66. CHENG Yuanfang, ZHANG Feng, WANG Jingyin, et al. Analysis of borehole collapse cycling time for shale[J]. Journal of China University of Petroleum (Science & Technology Edition), 2007, 31(1): 63?66.

[18] 鄧金根, 劉書杰, 石得勤, 等. 軟泥巖井眼彈塑性變形的拉格朗日元法計算[J]. 地質(zhì)力學學報, 1999, 5(1): 33?37. DENG Jingen, LIU Shujie, SHI Deqin, et al. Calculation of elastoplastic deformation of wellbore in soft mudstone using lagrangian method[J]. Journal of Geomechanics, 1999, 5(1): 33?37.

[19] 鄒靈戰(zhàn), 鄧金根, 徐顯廣, 等. 山前高陡構(gòu)造節(jié)理圍巖的井壁失穩(wěn)機制研究[J]. 巖石力學與工程學報, 2008, 27(S1): 2733?2733. ZOU Lingzhan, DENG Jingen, XU Xianguang, et al. Study on wellbore collapse mechanism in jointed rock masses of high-dip-structures before mountains[J]. Chinese Journal of Rock Mechanics and Engineering, 2008, 27(S1): 2733?2733.

[20] WANG H, SOLIMAN M Y, SHAN Z, et al. Understanding the effects of leakoff tests on wellbore strength[J]. SPE Drilling & Completion, 2011, 26(4): 531?539.

[21] BASSEY A, DOSUNMU A, IORKAM A, et al. A new (3D MUDSYSTTM Model) approach to wellbore strengthening while drilling in depleted sands: a critical application of LCM and stress caging model[C]// Nigeria Annual International Conference and Exhibition. Lagos, Nigeria: SPE, 2012: 1?18.

[22] 申光憲. 邊界元法[M]. 北京: 機械工業(yè)出版社, 1998: 138?224. SHEN Guangxian. Boundary element method[M]. Beijing: China Machine Press, 1998: 138?224.

[23] 張志增, 李仲奎. 橫觀各向同性巖體中圓形巷道反分析的惟一性[J]. 巖土力學, 2011, 32(7): 2066?2072. ZHANG Zhizeng, LI Zhongkui. Uniqueness of displacement back analysis of a circular tunnel in transversely isotropic rock mass[J]. Rock and Soil Mechanics, 2011, 32(7): 2066?2072.

[24] CHENEVERT M E, GATLIN C. Mechanical anisotropies of laminated sedimentary rocks[J]. SPE Journal, 1965, 5(1): 67?77.

[25] 馬天壽, 陳平. 邊界元法在深部地層井眼穩(wěn)定中的應(yīng)用[J].科技導報, 2015, 33(7): 49?54. MA Tianshou, CHEN Ping. Application of boundary element method to borehole stability problem in deep formation[J]. Science & Technology Review, 2015, 33(7): 49?54.

(編輯 劉錦偉)

Boundary element method and its application to borehole collapse problems in shale formations

MA Tianshou, CHEN Ping
(State Key Laboratory of Oil & Gas Reservoir Geology and Exploitation, Southwest Petroleum University, Chengdu 610500, China)

Based on the anisotropic continuum mechanics and the boundary element theory, the basic differential equations of the borehole collapse problem in transverse isotropic shale formation were established. The boundary integral equation of the basic equations and its discrete solution were deduced, and the discrete solutions of boundary element for stress and displacement were also deduced. The Mohr?Coulomb criterion was used for determining the stability of wellbore, and a basic method for solving borehole collapse by using the boundary element method (BEM) was formed the above processes. The physical model of borehole stability was established for shale formations, the isotropic model was used to verify the correctness of BEM, and the stress distribution around borehole and application example were analyzed by BEM. The results show that the BEM’s solutions are in good agreement with analytical solutions and their relative error is less than 2.49%. The anisotropy of elastic modulus, the differences of horizontal in-situ stress and the density of drilling fluid have significant influence on the stress distribution (the influence of the in-situ stress is the largest), and the influence of the anisotropy of poisson’s ratio is low. The analysis results of 1 530 m section in W201 well are in good agreement with the electric log data, the hole enlargement rate is about 9.0% by BEM, and the hole enlargement rate is about 10.0% by the electric logging data.

shale; borehole stability; transverse isotropy; boundary element method (BEM); numerical method

TE21

A

1672?7207(2016)03?0839?11

10.11817/j.issn.1672-7207.2016.03.017

2015?03?12；

2015?05?20

國家重點基礎(chǔ)研究發(fā)展計劃(973計劃)項目(2013CB228003)；四川省國際科技合作與交流項目(2016HH0001) (Project (2013CB228003) supported by the National Basic Research Program (973 Program) of China; Project(2016HH0001) supported by the Scientific Research Foundation of International Cooperation and Exchanges of Sichuan Province)

馬天壽，博士，講師，從事油氣井工程方面的教學與科研；E-mail: matianshou@126.com