朱帥，李明飛，姜麗萍，竇益華，王智勇

1.陜西鐵路工程職業(yè)技術(shù)學(xué)院（陜西 渭南 714000）

2.西安石油大學(xué) 機(jī)械工程學(xué)院（陜西 西安 710065）

0 引言

巖石的流變包括蠕變、松弛和彈性后效，巖石的流變影響油氣井等巖石工程結(jié)構(gòu)的穩(wěn)定性和安全性，因此，巖石流變特性的研究引起了人們的重視[1-2]。目前巖石流變研究的一般做法是：根據(jù)巖石流變的機(jī)理，建立描述巖石流變的模型和方程，實(shí)驗(yàn)測量模型和方程中的關(guān)鍵參數(shù)，再將包含這些參數(shù)的模型和方程用于巖石流變分析[3-4]。

目前，不同工程領(lǐng)域的學(xué)者結(jié)合自身學(xué)科特點(diǎn)，有針對性地開展巖石流變特征研究。肖欣宏等人針對水環(huán)境下的軟巖工程安全問題，研究了低應(yīng)力作用下的蠕變特性和該特性相同水壓條件下的瞬時(shí)應(yīng)變，得出水壓高低和應(yīng)力大小對蠕變特性的影響[5]。劉穎等人從工程施工及運(yùn)行實(shí)際中發(fā)現(xiàn)不同含水狀態(tài)下的巖石蠕變特性問題，對砂巖工程試樣進(jìn)行了蠕變力學(xué)試驗(yàn)，對比干燥狀態(tài)下的砂巖變化得出其受含水量的影響情況[6]。貴州大學(xué)的王唯等人為了探索巖石全過程蠕變特性，改進(jìn)傳統(tǒng)的Bingham模型，對比天然頁巖蠕變試驗(yàn)結(jié)果，驗(yàn)證其模型的適用性和可行性[7]。黃海軍對隧道圍巖在不同凍融循環(huán)次數(shù)下的蠕變特性進(jìn)行試驗(yàn)研究[8]。安徽理工大學(xué)的王游等為了克服傳統(tǒng)元件模型組合難以對巖石非線性蠕變特性具體描述的問題，改進(jìn)廣義Kelvin 模型進(jìn)行參數(shù)確定和驗(yàn)證，為類似問題研究提供思路和借鑒[9]。很多學(xué)者通過改變流變元件組合，改進(jìn)了巖石蠕變模型以驗(yàn)證具體工程領(lǐng)域的巖石蠕變問題[10-11]。高子璐等人給改進(jìn)模型的參數(shù)反演提供思路和方法[12]。國內(nèi)外相關(guān)專家學(xué)者也采用大邊界地層條件仿真方法對巖石流變現(xiàn)象進(jìn)行推演和性能研究[13-14]。復(fù)雜的工程背景下巖石的流變特性是復(fù)雜多變的，而西原體模型可較好地模擬巖石流變過程中的蠕變階段，這對研究流變地層巖石變化很有意義。

借鑒先進(jìn)的研究經(jīng)驗(yàn)，以流變巖層與井筒的相互作用為工程背景，針對流變巖層與井筒相互作用時(shí)的力學(xué)行為變化問題，基于改進(jìn)西原體模型的巖石蠕變，繼續(xù)探索巖石流變特性的后續(xù)松弛過程。

1 巖石流變本構(gòu)關(guān)系及其改進(jìn)

1.1 傳統(tǒng)西原體模型

巖石流變學(xué)理論主要研究巖石的應(yīng)力-應(yīng)變狀態(tài)及其隨時(shí)間增長的變化規(guī)律，建立恰當(dāng)?shù)谋緲?gòu)模型顯得尤為重要。如圖1 所示，傳統(tǒng)巖石力學(xué)西原體模型包括胡克體，黏彈性體（凱爾文體），理想黏塑性體，能夠反映巖石的彈性、黏彈性、黏塑性[1]。

圖1 傳統(tǒng)的西原體力學(xué)模型

文章公式中涉及到的應(yīng)力單位均為MPa，應(yīng)變無量綱，時(shí)間以d計(jì)。圖1中，σ為模型應(yīng)力；ε為應(yīng)變；σs為巖石的屈服應(yīng)力；k為彈簧的彈性系數(shù)；η為黏性元件的牛頓黏性系數(shù)；t為時(shí)間。其本構(gòu)方程和蠕變方程如下[1]：

1）本構(gòu)方程

當(dāng)σ＜σs時(shí)，

當(dāng)σ≥σs時(shí)，

2）蠕變方程

當(dāng)σ＜σs時(shí)，

當(dāng)σ≥σs時(shí)，

西原體模型可反映出當(dāng)應(yīng)力水平較低時(shí)，開始模型變形較快，一段時(shí)間后逐漸趨于穩(wěn)定，發(fā)生穩(wěn)定蠕變；當(dāng)應(yīng)力水平達(dá)到或超過材料某一臨界應(yīng)力值后，漸漸轉(zhuǎn)化為不穩(wěn)定蠕變。因此，該模型在巖石流變特性研究應(yīng)用中使用廣泛，特別在軟巖流變特性的描述中適用。傳統(tǒng)西原體模型的元件為理想的線性元件，只反映出蠕變的前兩個(gè)階段，不易完整描述巖石的其他流變特性。因此，為更全面地描述巖石的流變過程，在傳統(tǒng)西原體模型的基礎(chǔ)上，引入能夠表征巖石材料變形停止、應(yīng)力隨時(shí)間增長而下降的參數(shù)[3]。為便于說明問題，以鹽膏巖地層中油氣井井筒應(yīng)力變化為例進(jìn)行研究。

1.2 改進(jìn)后的流變模型本構(gòu)關(guān)系

馬克斯威爾體模型是由彈簧和阻尼器串聯(lián)組合成的一種黏彈性體，該模型能對流變特征中的蠕變現(xiàn)象和松弛現(xiàn)象較好表達(dá)，符合鉆井初期井壁巖石的流變規(guī)律。為此，結(jié)合前人對西原模型的改進(jìn)，將馬克斯威爾體與西原體模型結(jié)合起來，在模型的建立中串聯(lián)阻尼器，建立改進(jìn)的巖石流變西原體模型，如圖2 所示。該模型阻尼器觸發(fā)條件是當(dāng)模型應(yīng)變達(dá)到某一應(yīng)變值εA時(shí)，巖石進(jìn)入加速蠕變階段對應(yīng)的起始應(yīng)變。

圖2 改進(jìn)的西原體模型

圖2 中，σ為模型施加的總應(yīng)力,MPa，σs為巖石的屈服應(yīng)力，MPa；k為彈簧的彈性系數(shù)；η為黏性元件的牛頓黏性系數(shù)；ε為應(yīng)變。依據(jù)文獻(xiàn)[3]，得到改進(jìn)后的模型在加速蠕變階段的本構(gòu)方程為[3]：

其中τ=t- |tε=εA,t｜ε=εA為巖石進(jìn)入加速蠕變階段的時(shí)刻。

當(dāng)t=0 時(shí)，施加應(yīng)力σ0，并保持σ=σ0不變，得到該過程下的蠕變方程為：

為了驗(yàn)證改進(jìn)西原體模型分析巖石流變的準(zhǔn)確性，假設(shè)應(yīng)力大于摩擦片屈服應(yīng)力（σ＞σs），取彈簧元件的彈性系數(shù)k1、k2為10 MPa，牛頓黏性體的黏性系數(shù)η1、η2為10 MPa·d ，η3為50 MPa·d ，σ0為10 MPa，模型屈服應(yīng)力分別取90、80、70、60和50 MPa，得到不同屈服壓力下改進(jìn)西原體蠕變曲線，如圖3 所示。從蠕變曲線變化趨勢上看，巖石蠕變發(fā)生至第20 d時(shí)進(jìn)入加速階段，直至第26 d蠕變加速結(jié)束。整體看，改進(jìn)后的西原體模型完成了從蠕變減速到等速，再到加速的全過程，證明改進(jìn)模型對蠕變過程描述的適用性，可進(jìn)行后續(xù)流變特性的探索[15]。

圖3 改進(jìn)后的西原體模型蠕變曲線示意圖

1.3 改進(jìn)后流變模型的松弛方程及其特征

當(dāng)蠕變發(fā)生至一定程度時(shí)，巖石的變形將會保持在某一時(shí)刻的大小并且變化很小，此時(shí)認(rèn)為應(yīng)變ε為定值，此時(shí)，將模型的本構(gòu)方程進(jìn)行一階微分后得到：

等式兩邊同時(shí)積分可得：

當(dāng)t=0時(shí)，模型受到瞬時(shí)應(yīng)力σ0，可得常數(shù)計(jì)算式：

得到的模型松弛方程表達(dá)為：

為了與巖石蠕變現(xiàn)象銜接，從蠕變過程結(jié)束后對松弛方程進(jìn)行準(zhǔn)確性驗(yàn)證，保持模型相同賦值條件，對方程進(jìn)行模擬。模型分別賦予50、60、70、80和90 MPa 等屈服應(yīng)力值[15]。模擬出的松弛結(jié)果如圖4 所示，結(jié)合鉆井初期井壁巖石后續(xù)流變特性發(fā)展，對蠕變發(fā)生后一段時(shí)間內(nèi)的松弛方程進(jìn)行模擬，從蠕變?nèi)^程結(jié)束起認(rèn)為變形停止，從曲線中應(yīng)力隨時(shí)間的變化趨勢上看，模型受到的應(yīng)力從27 d 開始逐漸減小，符合松弛現(xiàn)象發(fā)生的應(yīng)力變化情況，說明該過程為巖石松弛階段。

圖4 改進(jìn)后的西原體模型松弛模擬結(jié)果

結(jié)合改進(jìn)模型后的松弛方程和模擬松弛曲線結(jié)果分析看，方程中的巖石應(yīng)力值在發(fā)生蠕變之后開始減小，松弛曲線中的應(yīng)力隨時(shí)間的增長呈下降趨勢。結(jié)合鉆井初期井壁巖石的流變規(guī)律可以推斷，當(dāng)巖石的應(yīng)變值達(dá)到松弛發(fā)生條件時(shí)刻對應(yīng)的應(yīng)變值時(shí)，此后巖石流變特征就進(jìn)入了松弛階段。在此之前，由圖3 可知模型處于蠕變階段的各個(gè)過程，在鉆井初期的流變地層中，由于該段地層巖石在鉆井后會發(fā)生應(yīng)力變形，達(dá)到一定程度會貼合井壁，可以看成經(jīng)歷了從蠕變開始到蠕變結(jié)束的過程，即受巖層自身地應(yīng)力下的應(yīng)變逐漸減小直至停止的過程?？梢耘袛嗔髯兊貙訋r石在蠕變之后產(chǎn)生松弛現(xiàn)象。為進(jìn)一步驗(yàn)證推斷，對松弛方程式（11）利用數(shù)學(xué)極限思想進(jìn)行求解計(jì)算，過程如下。

隨著時(shí)間的推移，當(dāng)t→∞時(shí)，松弛方程為“型”極限函數(shù)，隨即對方程進(jìn)行求極限得：

對于鉆井初期的井壁巖石流變的工程問題，鉆井導(dǎo)致的地層破壞，打破了流變地層原始地應(yīng)力平衡，為平衡受力，應(yīng)力則集中在井筒周圍，該狀態(tài)下的巖石就會發(fā)生蠕變現(xiàn)象。隨著蠕變?nèi)^程結(jié)束，地層與井筒緊密貼合，巖石停止蠕變變形，接著表現(xiàn)出應(yīng)力隨時(shí)間的增長而減小的松弛現(xiàn)象，巖石流變特征會因?qū)嶋H工程條件的改變發(fā)生變化。

2 地應(yīng)力計(jì)算及數(shù)值模擬

2.1 地應(yīng)力計(jì)算驗(yàn)證松弛方程

利用松弛方程進(jìn)行應(yīng)力計(jì)算時(shí)，存在黏性系數(shù)和彈性系數(shù)等5 個(gè)未知量，必須確定松弛方程中涉及到的這些系數(shù)。通常以實(shí)驗(yàn)的方式對巖石流變本構(gòu)關(guān)系中的方程進(jìn)行參數(shù)識別，現(xiàn)有文獻(xiàn)對該方面的實(shí)驗(yàn)研究較多，為方便計(jì)算，選用最貼近鉆井初期巖石流變特性的蠕變數(shù)據(jù)對松弛方程進(jìn)行參數(shù)識別，取長山巖的蠕變數(shù)據(jù)作為計(jì)算條件，改進(jìn)后得到的蠕變方程存在5 個(gè)參數(shù)[15]，相應(yīng)取5 組數(shù)據(jù)，其中100、200、300、350 和400 h 時(shí)的應(yīng)變值分別為0.051、0.054、0.058、0.063 和0.067 mm，初始應(yīng)力為14.72 MPa，屈服應(yīng)力為11.34 MPa，代入公式（11）可得到彈性系數(shù)和黏滯系數(shù)值分別為：

得到含具體參數(shù)的松弛方程：

為了論證松弛方程結(jié)論的準(zhǔn)確性，以某油田井區(qū)鹽膏巖地層中巖石力學(xué)參數(shù)為例，結(jié)合地應(yīng)力理論公式和有限差分?jǐn)?shù)值模擬兩種方法進(jìn)行該地層地應(yīng)力計(jì)算，結(jié)果與推導(dǎo)的流變松弛方程結(jié)果對比并分析誤差。根據(jù)某油田實(shí)際地質(zhì)數(shù)據(jù)和測井資料可知，該井區(qū)東西埋深在3 522～4 468 m，南北埋深在3 768～4 428 m 的鹽膏巖蠕變地層，符合巖石流變特性的地層最大厚度約220 m，最小厚度約142 m，大部分厚度平均在142～168 m，符合條件的巖層包括泥巖夾層和泥膏巖夾層。該段地層巖石平均密度為2.186 g/cm3，屈服強(qiáng)度測算為65.4 MPa。其具體巖石力學(xué)參數(shù)見表1。

表1 某油田地層巖石力學(xué)參數(shù)

根據(jù)上覆巖層壓力計(jì)算公式計(jì)算地層初始地應(yīng)力，可得所取地層垂直地應(yīng)力為112.05 MPa。結(jié)合我國地理?xiàng)l件現(xiàn)狀，根據(jù)平均水平應(yīng)力計(jì)算公式得116.4 MPa，水平主應(yīng)力系數(shù)比值取0.7[1]，得到的水平最大主應(yīng)力和水平最小主應(yīng)力分別為136.94、95.86 MPa。

結(jié)合該地層巖石力學(xué)參數(shù)及理論計(jì)算的地應(yīng)力結(jié)果，對含參數(shù)的松弛方程公式（13）進(jìn)行賦值計(jì)算地應(yīng)力。取屈服強(qiáng)度65.4 MPa，初始應(yīng)力116.4 MPa，得到應(yīng)力隨時(shí)間變化情況如圖5 所示，在26 h左右時(shí)應(yīng)力值與理論計(jì)算結(jié)果大小相對應(yīng)，水平應(yīng)力平均值隨即呈現(xiàn)出下降趨勢，直至30.9 h時(shí)，平均應(yīng)力值逐漸穩(wěn)定在113.08 MPa，與地應(yīng)力理論計(jì)算結(jié)果對比，平均誤差為3.83%。

圖5 松弛方程應(yīng)力隨時(shí)間的變化趨勢

2.2 數(shù)值模擬驗(yàn)證松弛方程

為了多角度驗(yàn)證流變松弛方程的準(zhǔn)確性和可行性，再進(jìn)行數(shù)值模擬計(jì)算。FLAC3D 有限差分軟件在大邊界模型下的數(shù)值計(jì)算較其他數(shù)值模擬軟件精度更高，計(jì)算更為準(zhǔn)確，故選用該軟件對上述井區(qū)流變地層簡化進(jìn)行建模計(jì)算。設(shè)定模型在水平方向上長度均為100 m，依據(jù)表1 中巖石力學(xué)參數(shù)，厚度方向自上而下按照20 m 泥巖層，30 m 鹽膏巖層、30 m 泥膏巖層進(jìn)行劃分，模型總厚度為80 m。對整體模型統(tǒng)一劃分網(wǎng)格，單元格個(gè)數(shù)為1 000，共1 000個(gè)正方體網(wǎng)格單元。

選用軟件自帶流變模型中的cwipp鹽巖變形模型代替鹽膏巖層，彈性模量取20 GPa，泊松比取0.31，體積模量取17.54 GPa，剪切模量取7.63 GPa，地層中的重力加速度9.8 m/s2，與表1保持一致。設(shè)置模型上表面壓力為112.05 MPa，水平最大主應(yīng)力為136.94 MPa，水平最小主應(yīng)力為95.86 MPa。采用深埋工程初始地應(yīng)力場的S-B 法計(jì)算生成，得到X方向上的水平應(yīng)力云圖，如圖6所示[15]。

圖6的計(jì)算結(jié)果顯示，模型的x方向水平應(yīng)力最大值為144.38 MPa；模型的y 方向水平應(yīng)力最大值為99.62 MPa。利用FLAC3D 有限差分?jǐn)?shù)值模擬結(jié)果對流變本構(gòu)關(guān)系的計(jì)算結(jié)果進(jìn)行校驗(yàn)，水平方向應(yīng)力平均誤差為7.25%。整體誤差較小，驗(yàn)證了巖石流變本構(gòu)松弛的準(zhǔn)確性，可以用來計(jì)算流變地層的應(yīng)力值。

圖6 初始地應(yīng)力場x方向應(yīng)力云圖

3 結(jié)論

1）改進(jìn)的西原體模型適用于流變地層中的巖石變化規(guī)律。

2）在巖石發(fā)生蠕變之后，推導(dǎo)出的松弛方程體現(xiàn)出應(yīng)力值隨著時(shí)間的增長而下降，說明巖石進(jìn)入松弛階段。

3）將流變地層物理特征代入松弛方程進(jìn)行計(jì)算，與地應(yīng)力計(jì)算結(jié)果相比誤差為3.83%；與FLAC3D有限差分軟件數(shù)值模擬計(jì)算結(jié)果相比，平均誤差為7.25%，推導(dǎo)的松弛方程可計(jì)算流變地層的應(yīng)力值。