泥頁巖多場耦合流變模型

王 萍，張燕明，呂 楊，梁凌云，崔云群，屈 展，劉利鋒，李樹生

(1.西安石油大學 石油工程學院，陜西 西安 710065； 2.陜西省油氣井及儲層滲流與巖石力學重點實驗室，陜西 西安 710065；3.低滲透油氣田勘探開發(fā)國家工程實驗室，陜西 西安 710018； 4.長慶油田分公司 油氣工藝研究院，陜西 西安 710018； 5.長慶油田分公司 第五采氣廠，陜西 西安 710018)

引 言

井壁失穩(wěn)問題經(jīng)歷了從力學研究、化學研究到力學與化學耦合研究，由Chenevert、Gray和Darley[1-3]提出：井壁失穩(wěn)的一個重要原因是泥頁巖水化，而不單單是純力學問題。此后，研究人員開始對黏土礦物的水化問題如泥頁巖水化對井壁穩(wěn)定性的影響展開研究。Chenevert[1]通過總吸附水量相關法對泥頁巖進行了力-化耦合的定量分析，Hale和Mody[4-5]通過半透膜滲透壓理論，將井底化學場的影響因素等效為孔隙壓力，首次提出等效孔隙壓力法。Tan[6-9]認為巖體內流體流動的動力主要靠孔隙壓力與滲透壓力共同作用下的總水勢差導致的，依據(jù)巖石中引入水介質的連續(xù)性方程進而推導出膨脹性巖體中的流體運動方程。雖然國內外學者就泥頁巖水化問題進行了大量的研究，但是由于問題本身具有復雜性，所得出的研究成果并不一致，仍然存在著許多疑問。如有些研究者認為總吸附水量影響著所有的水化變化，事實上用“總吸附水”概括這是不準確的，因為泥頁巖中還存在著與水化反應無關的水分子，如滲透水化水、離子水化水、孔隙自由水等。總水勢模型的溶質流動參數(shù)難以確定，且沒有給出解決的方法。等效孔隙壓力法以半透膜滲透壓為出發(fā)點，引入“非理想半透膜”，在力學分析中將化學因素的影響定量化，實現(xiàn)力學與化學的耦合計算[10-12]。

在泥頁巖井壁穩(wěn)定性分析中，通常從力學、化學以及力學與化學耦合幾個方面進行研究，然而這些研究往往很少考慮巖石的流變特性和初始微缺陷造成的損傷及其耦合發(fā)生過程中的演化擴展。泥頁巖與鉆井液相互作用，在其內部將會產(chǎn)生一系列的物理化學變化，改變巖石物質內部結構與力學性能，使巖石具有明顯的流變特征。隨著含水量的增加，巖石的流變更為明顯，在一定的條件下可能會造成井眼縮徑、卡鉆、坍塌。為此，本文用經(jīng)典半透膜滲透理論方法，通過將泥頁巖與鉆井液間的化學作用過程產(chǎn)生的應力引入到有效應力中，進行力學與化學的耦合計算，然后基于多孔介質有效應力原理，從巖石應力場、化學場、滲流場耦合分析的基本方程——平衡方程、連續(xù)性方程入手，建立泥頁巖應力場、化學場與滲流場耦合作用下的流變分析模型。

1 泥頁巖多場耦合流變模型

在建立泥頁巖的流變分析模型時假設：⑴ 所研究泥頁巖屬微小變形，泊松比不變；⑵ 巖體所有初始微缺陷統(tǒng)稱為“孔隙”，將泥頁巖視為等效連續(xù)介質或具有初始損傷的連續(xù)介質；(3)鉆井液是微可壓縮流體，滲流服從達西定律；(4)結合滲透系數(shù)與流變參數(shù)之間的聯(lián)系，在數(shù)值模擬計算時，每2個時步對滲透系數(shù)修正一次[13-14]。

對于單重孔隙介質，最常采用的有效應力張量形式為

σij*=σij-αpδij，

(1)

式中：σij*是有效應力，Pa；σij是施加正應力，Pa；p為孔隙流體壓力，Pa；δij為Kronecker符號；α為比奧系數(shù)(0<α≤1)。

根據(jù)泥頁巖與鉆井液間化學勢差造成的應力,得到遠場孔隙壓力

(2)

式(2)負號可為正或負，因為近井地帶的巖石孔隙壓力可能比遠場地層中的孔隙壓力大，也可能小，要視鉆井液與泥頁巖中的水活度大小而定。

將式(2)代入式(1)得到在應力場、化學場與滲流場耦合作用條件下，巖石的有效應力[15]

(3)

(4)

在流變變形過程中，由總應變與黏彈性應變、黏塑性應變的關系可得

(5)

將式(5)代入式(4)，可得到多場耦合時的黏彈塑性本構關系：

(6)

由巖石流變中的西原模型式(如圖1)可知，在黏彈性體N//H(//表示兩個原件并聯(lián))中，應力σ是應力σH與應力σN之和，而兩者的應變相等。

圖1 巖石流變模型Fig.1 Rheological model of rock

即

(7)

整理得

(8)

若t=0時，εve=0，則在常應力σ=σs作用下求解方程

(9)

式中，E1為黏彈性體的黏彈性模量，η1為黏彈性體的黏滯系數(shù)。

則

(10)

式中，K*=K(1-DT)和G1*=G1(1-DT)為考慮初始水化損傷時，巖石蠕變的各項參數(shù)，其中，DT為巖石的初始水化損傷量，體積模量和黏性剪切模量分別為

式中，E0為彈性模量。

對于黏塑性變形部分，采用有效應力表述黏塑性蠕變。巖石黏塑性蠕變率可表示為

(11)

其中損傷參數(shù)Dc可由損傷率[11]

(12)

確定。式中，A、n為巖石材料常數(shù)，由實驗確定；Dc為巖石蠕變損傷，σs為巖石屈服應力。

對式(12)進行積分得

(13)

據(jù)文獻[16]、[17]并結合式(11)，則有

(14)

(15)

式中，F(xiàn)為巖石屈服函數(shù)，F(xiàn)0為巖石屈服函數(shù)的初始參考值。

將式(10)和式(15)代入式(6)，得到多場耦合流變本構方程為

(16)

2 流變分析的有限元格式

2.1 平衡方程

式(6)的矩陣格式為

{σ}=[D][?]{f}-[D]({εve}+{εvp})-α{M}p+α{M}β，

(17)

可得到單元內任意點的總應力矢量的近似解為

(18)

(19)

式中，[?]為算子矩陣，{σ}、{εe}分別為t時刻巖石內任意一點的應力和彈性應變，[D]為彈性矩陣，{M}=[1 1 1 0 0 0 ]T，

根據(jù)虛功原理可導出分析系統(tǒng)的平衡方程

(20)

將式(19)代入式(20)，得到

{Ru}e，簡寫為

[kuu]{δ}e+[kup]{p}e+[kuβ]={fve}e+{fvp}e+{Ru}e。

(21)

在Δtn=tn+1-tn內的位移、水壓增量分別為{Δδ}3、{Δp}e，而黏彈性應變增量和黏塑性應變增量分別為{Δεve}e和{Δεvp}e，則平衡方程(21)的增量形式為

[kuu]{Δδ}e+[kup]{Δp}e+[kuβ]={Δfve}e+{Δfvp}e+{ΔRu}e，

(22)

式中：

對所有單元建立式(22)，形成整體平衡方程的增量形式為

[Kuu]{ΔU}+[Kup]{ΔP}+[Kuβ]={ΔRu}，

(23)

式中

{ΔRu}=∑({Δfve}e+{Δfvp}e+{ΔRu}e)。

2.2 連續(xù)性方程

泥頁巖巖石是由巖石骨架、孔隙和流體所組成的多孔介質。一方面，巖石骨架與孔隙流體之間的摩擦力使孔隙流體的排出受阻，巖石變形延滯；另一方面，孔隙結合水的黏滯性使得骨架的變形有一個過程。由于巖石中水的存在和遷移使巖石變形特性既不是彈性體，也不是塑性體，而是黏彈塑性體。所以在泥頁巖巖石發(fā)生蠕變的同時必然伴隨孔隙中水的遷移。根據(jù)混合理論，分別建立巖石骨架、孔隙以及流體的控制方程，其次通過考慮方程之間的相互關聯(lián)，最后得出巖石的連續(xù)性方程。因此，連續(xù)性方程可根據(jù)流體運動方程、物性狀態(tài)方程及泥頁巖骨架和流體的質量守恒方程來建立。

當有內部或外部的流體源時，上式變?yōu)?/p>

(24)

根據(jù)水壓邊界條件p=0和流速邊界條件vn=0，將這2個邊界條件代入達西定律表達式得

(25)

若某邊界上法向流速或比流量已知，則該邊界條件為

lvx+myy+nvz=vn。

(26)

代入達西定律得

(27)

式中vn為邊界流體源矢量。將邊界條件式(26)寫成矩陣形式

(28)

(29)

由流速邊界條件式(29)，利用插值函數(shù)進行離散，由伽遼金變分原理得到

(30)

式中：

設tn到tn+1時刻單元結點的參數(shù)隨時間變化后，采用時間積分的一般格式

(31)

(32)

式中θ為時間積分因子，其值域為0～1。

對所有單元建立式(32)，形成整體連續(xù)性方程

[Kup]T{ΔU}+[Kpp]{ΔP}={ΔRp}。

(33)

其中：

[Kpp]=∑([kp]+θΔt[kpp])；

{ΔRp}=∑Δt({Rp}e+θ{ΔRp}e-[kpp]{p}ne)。

2.3 總體控制方程

建立的多場耦合流變有限元平衡方程和連續(xù)性方程均具有耦合項，需聯(lián)立求解[18]。因此，總體控制方程將由方程(23)和方程(33)組成，即

(34)

式(34)所給出的巖石流變多場耦合方程是一種應力場、化學場和滲流場耦合流變分析的統(tǒng)一形式，適用對符合西原流變模型的泥頁巖進行多場耦合分析。考慮到泥頁巖地層中的實際情況，針對不同的泥頁巖井壁圍巖耦合情況，當巖石受力狀態(tài)或變形特性發(fā)生變化時，只需改變方程便可利用式(34)中的{ΔRu}計算格式，而連續(xù)性方程(33)保持不變，便可利用式(34)進行相應的多場耦合有限元分析。

若對泥頁巖井壁進行黏彈性蠕變耦合分析，則只需取{Δεvp}=0，即有{Δfve}e=0，修正{ΔRu}，便可得到多場耦合的黏彈性分析計算公式：

{ΔRu}=∑({Δfve}e+{ΔRu}e)。

(35)

若對泥頁巖井壁進行黏彈塑性蠕變耦合分析，根據(jù)式(19)有

后得到

{ΔRu}=∑({Δfve}e+{Δfvp}e+{ΔRu}e) 。

(36)

3 泥頁巖多場耦合流變數(shù)值分析

泥頁巖多場耦合模型的總體控制方程可由離散化的平衡方程和連續(xù)性方程得到，進而得出總體控制方程的有限元格式。通過二次開發(fā)子程序引入到ABAQUS中，并進行泥頁巖井壁蠕變損傷失穩(wěn)數(shù)值模擬分析[19-21]。

取長慶油田的長7層西峰233井區(qū)某井進行計算，其地層參數(shù)如下：井深H=3 400 m； 井眼半徑rw=231 mm；地應力σH=38.5 MPa；最大水平主應力σ1=44.2 MPa；最小水平主應力σ3=35.8 MPa；地層原始內聚力C=18.7 MPa；地層原始內摩擦角φ=31.58°；地層泊松比μ=0.25；彈性模量E=31.25 GPa；吸水擴散系數(shù)Cf=0.013 4 cm2/h；有效應力系數(shù)α=0.4；吸水膨脹系數(shù)K1=0.033 3，K2=0.832。流變耦合分析采用西原模型，其中黏彈性模量E1=4 065 MPa，黏彈性系數(shù)η1=6 831 MPa·d，黏塑性系數(shù)η2=6 005 MPa·d。由實驗測得不同時間的水化損傷量，見表1。

表1 不同時間巖石的水化損傷量Tab.1 Rock hydration damage at different times

由于是井眼圍巖是軸對稱形式，取1/4進行計算分析，建立模型長200 m、寬200 m、井徑2 cm。經(jīng)分析，邊界對結果影響可以忽略，其模型和網(wǎng)格劃分如圖2。具體分析為：第一步，地應力平衡；第二步，鉆進至目的層；第三步，圍巖蠕變35 d。

圖2 圍巖蠕變模型與網(wǎng)格劃分Fig.2 Creep model of surrounding rockand its meshing

3.1 應力場變化

圖3和圖4分別為蠕變損傷分析開始到蠕變分析結束的井壁圍巖的應力分布云圖。由圖3可知，在剛鉆進至儲層時，近井地帶處于地應力平衡，此時，井內流體沒有與儲層發(fā)生接觸。從圖4可看出，泥頁巖井壁圍巖在多場耦合效應的作用下水化嚴重，出現(xiàn)垮塌現(xiàn)象。

圖3 初始蠕變地應力平衡階段Fig.3 Initial creep in-situ stress equilibrium stageS:各方向應力值

圖4 蠕變損傷應力分布云圖Fig.4 Stress nephogram of creep damageS:各方向應力值

井壁圍巖蠕變損傷導致其變形具有明顯的時效性和非線性，其蠕變損傷演化規(guī)律如圖5所示。由圖可直觀地反映初始蠕變、等速蠕變和加速蠕變各階段損傷值與時間的關系，體現(xiàn)了非線性蠕變損傷模型的特點。

圖5 蠕變損傷演化曲線Fig.5 Creep damage evolution curve

巖石的蠕變損傷與其內部微裂紋的延伸和擴展密切相關，宏觀表現(xiàn)為蠕變過程中的體積擴容。導致體積擴容的原因是由于井眼周圍應力場的變化，即周向應力和徑向應力的變化。多場耦合蠕變損傷產(chǎn)生的井周應力如圖6所示。在加載階段和初始變形階段(t≤15 d)，不考慮蠕變損傷耦合分析和考慮蠕變損傷耦合分析中的徑向應力和周向應力變化差別均不大，但在t>25 d以后，考慮蠕變損傷的井周應力變化比不考慮蠕變損傷的應力峰值大，且隨時間差別擴大。這可以解釋現(xiàn)場鉆井不是立即坍塌，而是鉆開后幾天甚至半個月才出現(xiàn)破壞，此種破壞與時間有關，對鉆井的危害也很大。

圖6 多場耦合效應井周應力變化曲線Fig.6 Stress curves around well at different time under multi-field coupling effect

近井壁處由于長時間受多場耦合效應作用，井周應力變化劇烈，考慮蠕變損傷的井周應力與不考慮蠕變損傷的井周應力具有明顯的差別，應力變化具有非線性和時效性。應力最大值不在井壁上，而是在近井壁地帶，隨著與井眼間距離的增大徑向應力和周向應力開始減小且趨于穩(wěn)定。因此在井眼附近會形成應力集中，造成井眼破壞，導致井壁失穩(wěn)。

3.2 孔隙壓力的變化

本文將泥頁巖與鉆井液間的化學作用過程產(chǎn)生的應力作用等效為孔隙壓力，即化學場的耦合效應由孔隙壓力的變化來表征。耦合效應下孔隙壓力與時間、井眼距離的變化關系如圖7、圖8所示。由圖7可知，不考慮蠕變損傷的孔隙壓力隨著時間增大明顯，后逐漸趨于平緩。而考慮蠕變損傷的孔隙壓力先隨著時間增大而后開始緩慢下降。由圖8可知，距離井眼越遠孔隙壓力越大，考慮蠕變損傷前期的孔隙壓力較之不考慮蠕變增幅較快，后期增幅減小。巖石骨架的變形，導致孔隙結構和巖石孔隙體積的改變，造成孔隙壓力發(fā)生改變。不考慮蠕變損傷時則忽略了巖石骨架變形的影響，相對于考慮蠕變損傷時孔隙壓力大。這種變形是孔隙流體的長時間作用下由多孔介質的流變效應導致的，因此隨著時間的增加孔隙壓力差別越明顯。

圖7 孔隙壓力-時間曲線Fig.7 Pore pressure-time curves

圖8 孔隙壓力-井眼距離曲線Fig.8 Pore pressure-borehole distance curves

3.3 滲透率的變化

滲透率的變化特性也可以直接或間接地反映出巖石內部滲流場的分布，如圖9所示。巖石的蠕變損傷與其內部微裂紋的延伸和擴展密切相關，宏觀表現(xiàn)為蠕變過程中的體積擴容。由圖9可知，不考慮蠕變損傷的井壁圍巖應變、滲透率隨著鉆進時間增大，這是因為泥頁巖吸水后發(fā)生膨脹，應力增大裂隙或孔隙發(fā)生擴展，導致體積應變增大，同時形成一定的滲透通道。

圖9 滲透率/應變-時間曲線Fig.9 Permeability/strain-time curves

考慮蠕變損傷，在初期蠕變變形階段，巖石內微裂紋隨機擴展增加，當微裂紋的擴展數(shù)量及長度達到一定程度后不再發(fā)生顯著變化，則使得滲透率變化趨勢變小，也說明蠕變變形進入穩(wěn)態(tài)階段，同時這也解釋了在長期地應力作用下巖石內孔隙和微裂隙受壓閉合而導致滲透率變小的現(xiàn)象。當滲透率隨著時間推移達到某一最低值后會有一段增幅期，進入了加速蠕變階段，在這一階段巖石裂紋在隨著水平應力逐漸變大的同時，其內部微裂隙繼續(xù)迸發(fā)、延伸、貫通及擴展成為宏觀裂縫。在非線性蠕變損傷變形期，隨著軸向應力增大，巖石滲透率緩慢增加，當圍巖發(fā)生失穩(wěn)破壞時出現(xiàn)較大的階躍。因此，巖石蠕變損傷(微破裂)演化引起滲透率的變化與巖體的蠕變損傷是一致的。

4 結 論

(1)井壁圍巖蠕變損傷導致其變形具有明顯的時效性和非線性特性，井周應力在整個蠕變損傷階段的變化比不考慮蠕變損傷的應力峰值大，且隨時間增加它們之間的差別擴大，此種破壞與時間有關，對鉆井的危害大。

(2)應力最大值不在井壁上，而是在近井壁地帶，隨著與井眼距離的增大，徑向應力和周向應力開始減小且趨于穩(wěn)定。在井眼附近會形成應力集中，造成井眼破壞，導致井壁失穩(wěn)。

(3)在孔隙流體的長時間作用下多孔介質產(chǎn)生流變效應，導致了巖石孔隙結構和孔隙體積的改變，造成孔隙壓力發(fā)生改變?？紫秹毫﹄S時間、距離的增長而增大后趨于平衡。

(4)在整個蠕變階段，圍巖的滲透率先減小后增大。巖石蠕變損傷發(fā)展引起滲透率的變化與巖石的裂縫擴展是一致的。