應用邊界元法分析頁巖地層井眼坍塌問題

馬天壽，陳平

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

應用邊界元法分析頁巖地層井眼坍塌問題

馬天壽，陳平

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

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

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

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

1　頁巖井眼穩定問題分析方法

1.1數學模型

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

平衡方程為

幾何方程為

橫觀各向同性頁巖的本構方程[23]

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

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

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

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

再將式(7)代入式(3)，可得平面應變問題本構方程

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

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

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

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

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

1.2邊界元求解方法

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

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

式中：dijkl為實際彈性常數；為參照態彈性常數；為實際與參照態彈性常數的差值。因此，本構方程可表示為

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

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

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

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

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

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

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

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

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

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

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

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

若令

則式(25)可寫成

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

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

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

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

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

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

其中：Dkl和Skl為系數矩陣；Dkij和Skij為與邊界條件有關的系數；δij為克羅內克函數；r為矢徑長度。

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

1.3強度判別準則

在進行井眼穩定性判別前，需要將廣義平面二向應力狀態結果轉換為空間主應力狀態，即計算出空間應力狀態下的主應力。采用BEM求解得到的應力張量σ為

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

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

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

當采用Mohr?Coulomb準則判別井眼不穩定時，假設頁巖地層強度系數(strength factor，FS)為

當FS＞1時，井周頁巖地層穩定；當FS＜1時，井周頁巖地層將發生失穩。因此，FS越小，則表示井眼失穩的風險越大。

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

式中：τ為剪切面上剪應力；σn為剪切面上正應力；c0為巖石內聚力；?為巖石內摩擦角；pp為地層孔隙壓力；α為Biot系數。

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

2.1求解模型的建立

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

圖1　井周地層計算網絡模型Fig. 1 Mesh model of formation around borehole

表1　井周地層計算基礎參數Table 1 Basic parameters of formation around borehole

2.2模型驗證

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

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

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

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

3　計算結果討論

3.1彈性模量對井壁應力的影響

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

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

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

3.2泊松比對井壁應力的影響

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

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

3.3地應力對井壁應力的影響

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

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

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

3.4井筒內鉆井液密度對井壁應力的影響

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

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

4　應用實例

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

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

由圖8和圖9可以看出：

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

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

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

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

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

5　結論

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

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

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

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

[1] 郭彤樓, 張漢榮. 四川盆地焦石壩頁巖氣田形成與富集高產模式[J]. 石油勘探與開發, 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掃描技術研究頁巖水化細觀損傷特性[J]. 石油勘探與開發, 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] 王紅巖, 劉玉章, 董大忠, 等. 中國南方海相頁巖氣高效開發的科學問題[J]. 石油勘探與開發, 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] 王倩, 周英操, 王剛, 等. 泥頁巖井壁穩定流固化耦合模型[J].石油勘探與開發, 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] 王倩, 周英操, 唐玉林, 等.泥頁巖井壁穩定影響因素分析[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] 賈善坡, 鄒臣頌, 王越之, 等. 基于熱?流?固耦合模型的石油鉆井施工過程數值分析[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] 李軍, 陳勉, 金衍, 等. 大位移井井壁穩定三維彈塑性有限元分析[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] 李嗣貴, 鄧金根, 李明志. 節理破碎地層井壁穩定的離散元分析[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] 屈平, 申瑞臣, 付利, 等. 三維離散元在煤層水平井井壁穩定中的應用[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] 尹虎, 李黔, 郭亮, 等. 離散單元法評價煤層井壁穩定性[J].天然氣工業, 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] 鄧金根, 郭東旭, 周建良, 等. 泥頁巖井壁應力的力學?化學耦合計算模式及數值求解方法[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]. 地質力學學報, 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] 鄒靈戰, 鄧金根, 徐顯廣, 等. 山前高陡構造節理圍巖的井壁失穩機制研究[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]. 北京: 機械工業出版社, 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] 馬天壽, 陳平. 邊界元法在深部地層井眼穩定中的應用[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

國家重點基礎研究發展計劃(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