基于廣義Kelvin模型的非定常鹽巖蠕變模型

韓偉民，閆怡飛，閆相禎

(1.中國石油大學(華東)儲運與建筑工程學院，山東青島，266580；2.中國石油大學(華東)機電工程學院，山東青島，266580)

四川盆地作為油氣資源的重要賦存地，分布著大量厚度不一的鹽巖地層。而位于川東北三疊系嘉陵江組鹽巖地層的多口油氣田生產(chǎn)井出現(xiàn)了套管變形損壞情況[1]，因此，研究該地層段鹽巖的流變特性就顯得尤為重要。當應力水平低于巖石長期強度時，其蠕變曲線僅存在衰減蠕變階段和穩(wěn)態(tài)蠕變階段；當應力水平高于巖石長期強度時，其蠕變曲線呈現(xiàn)完整的3個階段，即衰減蠕變階段、穩(wěn)態(tài)蠕變階段和加速蠕變階段。室內試驗結果表明鹽巖的蠕變曲線也存在上述3個階段[2]，但對于工程問題而言，處于高應力狀態(tài)下的深層鹽巖體蠕變過程中并不會進入加速蠕變階段而發(fā)生損傷破壞，因此，研究衰減蠕變階段和穩(wěn)態(tài)蠕變階段鹽巖的非線性蠕變特征對工程實際更有意義。其中，衰減蠕變階段的存在時間相對短暫，在恒定載荷作用下穩(wěn)態(tài)蠕變階段則更為持久。國內外學者針對鹽巖的蠕變特征開展了大量研究[3-6]，但由于地質環(huán)境、礦物成分及含水率等影響因素的存在，不同地區(qū)的鹽巖彈性參數(shù)及流變特性存在明顯差異，因此，單一的蠕變本構模型形式并不具有普適性。目前針對鹽巖的非線性流變問題，主要有以下幾種常見研究方法[7]：1)以組合的黏彈性或黏塑性元件模型為基礎，在其后串聯(lián)加入非線性元件，蠕變參數(shù)可根據(jù)室內試驗結果來確定。劉開云等[8]將Maxwell模型與非線性黏塑性體串聯(lián)組成新的非線性黏塑性模型以描述泥巖的蠕變過程。一般的工程問題可采用該方法進行近似處理，但對于復雜的巖石蠕變本構，組合的黏彈性或黏塑性元件模型并不能反映與蠕變時間和應力水平相關的非線性蠕變特征，此時，該方法的近似處理結果并不適用。2)將元件模型中的線性元件視為時間的函數(shù)，并以非線性元件進行替代，改進后的非線性元件組合模型的蠕變參數(shù)可根據(jù)室內試驗結果進行確定。王軍保等[9-10]用非線性黏壺代替常規(guī)Burgers模型的線性黏壺，通過考慮損傷因子和蠕變對黏滯系數(shù)的影響等，建立了改進的非線性元件組合模型以描述鹽巖的非線性蠕變特性。ZHOU等[11-12]用分數(shù)階Abel黏壺元件替代西原模型的線性黏壺，建立了改進的西原模型以反映鹽巖流變的3個階段。該方法雖較為理想，但模型中蠕變參數(shù)辨識的計算處理相對繁瑣和困難，且改進后的非線性元件模型并不能反映溫度因素對巖石蠕變過程的影響。3)基于室內試驗結果，考慮損傷和斷裂力學，建立經(jīng)驗本構關系式。高小平等[13]應用不同應力和溫度條件下的鹽巖變形機制，得到了鹽巖的穩(wěn)態(tài)蠕變應變率與偏應力、能量與溫度之間的函數(shù)關系式。楊春和等[14-15]通過不同應力水平作用下鹽巖蠕變變形與時間的關系，建立了鹽巖瞬態(tài)與穩(wěn)態(tài)蠕變的耦合本構方程。王安明等[16-17]采用Norton Power鹽巖蠕變本構模型分別對鹽穴儲氣庫及鹽巖井段井眼縮徑進行了計算分析。該方法可充分利用試驗數(shù)據(jù)進行回歸分析得到較為可靠的本構關系式，但由于巖石取芯費用較高且室內蠕變時間較長，短期內較難獲得大量可靠的室內試驗數(shù)據(jù)。鑒于此，為了避免上述幾種非線性蠕變模型構建方法的局限性，根據(jù)川東北油氣田鹽巖三軸蠕變試驗結果，并參考其他學者在研究鹽巖蠕變特性中采用的方法[7,10-12,18]，本文作者將改進后的非線性元件組合模型與非線性元件進行串聯(lián)，構建新的適用于川東北油氣田嘉陵江組鹽巖的非線性蠕變本構模型，以便能夠更好地描述蠕變時間、溫度和應力水平對鹽巖衰減蠕變階段和穩(wěn)態(tài)蠕變階段的非線性蠕變特征的影響。首先，在廣義Kelvin模型的基礎上用非定常黏壺替代線性黏壺構建非定常廣義Kelvin模型，將其與能夠描述鹽巖穩(wěn)態(tài)蠕變特征的Heard模型進行串聯(lián)組合，構建新的NGKH非線性鹽巖蠕變模型，以描述鹽巖的瞬態(tài)蠕變和穩(wěn)態(tài)蠕變過程；然后，基于套損現(xiàn)象集中發(fā)生的嘉陵江組鹽巖地層的巖心室內蠕變試驗結果，采用Levenberg-Marquardt算法分別對不同應力水平下的NGKH模型、非定常廣義Kelvin模型和廣義Kelvin模型參數(shù)進行辨識，通過比較驗證新模型的適用性；最后，基于FLAC3D提供的接口程序對NGKH模型進行二次開發(fā)，并通過數(shù)值試驗對其正確性進行驗證。

1 鹽巖蠕變本構模型

1.1 非定常廣義Kelvin模型的構建

廣義Kelvin模型是由胡克體與Kelvin模型串聯(lián)而成，如圖1所示。

圖1 廣義Kelvin模型Fig.1 Generalized Kelvin model

廣義Kelvin模型的蠕變方程為[18]：

式中：E0為串聯(lián)胡克體的彈性模量；E1為Kelvin體的黏彈性模量；η1為Kelvin體的初始黏滯系數(shù)。

由于廣義Kelvin模型只能描述巖石的線性蠕變特征，對于蠕變特性復雜的鹽巖來說并不適用。將Kelvin體中牛頓黏壺的黏滯系數(shù)視為時間的函數(shù)，建立非定常黏壺[19]：

式中：η(t)為非定常黏壺的黏滯系數(shù)；a為待定系數(shù)。

對式(2)求導可得

由式(3)可知：當a>0時，˙(t)≤0，黏滯系數(shù)隨時間衰減；當a=0時，η(t)=η，非定常黏壺退化為線性黏壺；當a<0時，˙(t)≥0，黏滯系數(shù)隨時間增大。因此，可通過變化a來實現(xiàn)非定常Kelvin體的非線性蠕變特征。用上述非定常黏壺替換圖1中廣義Kelvin模型的線性黏壺，構建而成的非定常廣義Kelvin蠕變模型如圖2所示。

圖2 非定常廣義Kelvin模型Fig.2 Non-stationary generalized Kelvin model

非定常Kelvin體的狀態(tài)方程為

式中：˙為蠕變速率。

分離變量求定積分，并由ε=J(t)σ可得其蠕變柔量JK為

由式(1)和式(5)可得非定常廣義Kelvin模型的蠕變柔量為

則非定常廣義Kelvin模型的一維本構方程為

在三維應力狀態(tài)下，巖石的應力張量可分解為偏應力張量Sij和球應力張量σm。通常考慮應力偏量僅引起蠕變變形，球應力引起彈性變形，同時假定巖石為各向同性體，由式(1)可得廣義Kelvin模型的三維蠕變方程為[19]

式中：εij為應變張量；K和G0分別為串聯(lián)胡克體的體積模量和剪切模量；G1為Kelvin體的剪切模量。

在等圍壓三軸壓縮試驗中σ2=σ3，代入式(8)可得常規(guī)三軸壓縮試驗中廣義Kelvin模型的軸向應變方程：

同理，由式(7)和式(8)可得非定常廣義Kelvin模型的三維蠕變方程和軸向應變方程分別為[19]：

1.2 Heard模型的構建

研究鹽巖的穩(wěn)態(tài)蠕變本構方程即確定其穩(wěn)態(tài)蠕變速率與溫度、偏應力之間的變化關系，鹽巖的蠕變機制和蠕變本構方程隨其所處地質環(huán)境的不同而變化[20]。對于川東北油氣田的深部鹽巖地層來說，其蠕變機制屬于位錯滑移機制[6]，符合Heard本構模型，本構方程如式(12)所示。曾德智等[21-22]利用Heard蠕變模型對鹽巖地層套管蠕變外載、鹽巖層井壁穩(wěn)定性及井眼縮徑量進行了計算分析；LI等[6]通過對鹽巖和鹽膏巖進行的室內三軸蠕變試驗，分析了不同礦物成分和圍壓對鹽巖蠕變速率的影響，并采用Heard模型對試驗結果進行了擬合。

式中：為穩(wěn)態(tài)蠕變速率；Q為激活能；R為摩爾氣體常量，R=8.314 J/(mol·K)；σ為偏應力；T為熱力學溫度，K；A和B為流變常數(shù)。

按照GB/T 50266—2013“工程巖體試驗方法標準”，將取自川東北某氣田嘉陵江組地層埋深3 838～3 858 m處的巖石巖芯，加工成高為100mm、直徑為50mm的標準圓柱體巖石試樣共2組。該地層段鹽巖層厚度為174m，地層巖性以灰白色硬石膏巖、白色鹽巖、灰白色膏鹽巖為主，夾灰色灰質白云巖、深灰色、灰色白云巖、膏質灰?guī)r、泥灰?guī)r。本試驗采用MTS材料試驗機進行，該裝置由軸向加壓、側向加壓、溫控設備和微機系統(tǒng)組成。對第一組巖石試樣，在室溫條件下將圍壓加載到15MPa并保持不變，軸向采用單級加載方式分別施加30，35，40和45MPa的均布載荷并保持不變，加載時間持續(xù)40 h左右。由前人研究成果可知：溫度因素對鹽巖的蠕變特性影響非常顯著，因此，對于第二組巖石試樣，首先將三軸室的溫度分別加到50℃和100℃并保持3 h左右，再將圍壓加載至15MPa并保持不變，軸向仍然采用單級加載方式施加30MPa均布載荷并保持不變，加載時間持續(xù)30 h左右。圖3(a)和圖3(b)所示分別為圍壓15MPa時，不同偏應力及不同溫度下的鹽巖應變-時間曲線?；贖eard蠕變模型，采用Origin分析軟件對圖3中不同條件下的穩(wěn)態(tài)蠕變率進行非線性擬合，得到的參數(shù)A，B和Q見表1。

圖3 鹽巖室內試驗蠕變曲線Fig.3 Creep curves of salt rock

表1 擬合得到的Heard模型蠕變參數(shù)Table1 Creep parameters of fitted Heard model

1.3 NGKH模型的構建

王軍保等[18]基于芒硝的室內三軸蠕變壓縮試驗，將分數(shù)階廣義Kelvin模型與未考慮溫度因素的Heard模型相結合，構建了能夠反映芒硝非線性蠕變特性的蠕變本構模型。而在不同地層深度和地層溫度下，鹽巖的蠕變過程受溫度影響較為顯著，因此，本文采用的Heard模型中仍然考慮溫度影響因素。

上述建立的非定常廣義Kelvin模型能夠描述鹽巖的瞬時變形及衰減蠕變狀態(tài)，而Heard蠕變模型能夠反映鹽巖的偏應力、溫度與穩(wěn)態(tài)蠕變速率的關系，因此，參考巖石非線性蠕變模型的經(jīng)典構建思路，將非定常廣義Kelvin模型與Heard蠕變模型進行串聯(lián)，構成新的鹽巖非線性蠕變模型，簡稱NGKH蠕變模型(見圖4)，下式即為NGKH蠕變模型的軸向應變方程：

式(13)中第1項和第2項反映初始蠕變階段的瞬時應變，第3項反映衰減蠕變階段的黏彈性應變，第4項反映穩(wěn)態(tài)蠕變階段的黏性流動應變。

圖4 NGKH蠕變模型Fig.4 NGKH creep model

1.4 蠕變參數(shù)敏感性分析

在應力狀態(tài)、溫度相同的條件下，對式(13)中的待定系數(shù)a和A進行敏感性分析，圖5(a)和圖5(b)所示分別為系數(shù)a和A對NGKH模型蠕變曲線的影響。對式(3)進行分析可知a的變化能改變Kelvin體中非線性黏壺的蠕變性質。由圖5(a)可以看出a的變化對蠕變曲線的形狀尤其是衰減蠕變階段影響較為顯著：當a>0時，隨著a增大，衰減蠕變階段的蠕變速率衰減速度加快，進入穩(wěn)態(tài)蠕變階段的時間提前，此后各蠕變曲線幾乎重合；當a<0時，隨著a減小，衰減蠕變階段的蠕變速率衰減速度同樣加快，進入穩(wěn)態(tài)蠕變階段的時間同樣提前，但此后各蠕變曲線近似平行，軸向應變呈近似線性減小。從圖5(b)可知：A的變化主要影響穩(wěn)態(tài)蠕變階段的曲線形態(tài)，對衰減蠕變階段則影響很小。隨著A增大，穩(wěn)態(tài)蠕變率明顯增大，軸向應變也明顯增大。蠕變參數(shù)B對蠕變曲線的影響規(guī)律與A的影響規(guī)律相似。

圖5 蠕變參數(shù)a和A對NGKH模型蠕變曲線的影響Fig.5 Influence of parameters a and A on creep curves of NGKH model

將圍壓固定為15MPa，在其他蠕變參數(shù)不變的情況下，分別觀察應力σ1的變化對NGKH模型及廣義Kelvin模型蠕變曲線的影響，如圖6所示。

圖6 軸向載荷對NGKH模型和廣義Kelvin模型蠕變曲線的影響Fig.6 Effect of axial load on creep curves of NGKH model and generalized Kelvin model

從圖6可以看出：隨著軸向載荷增大，NGKH模型及廣義Kelvin模型的軸向應變均隨之增大，但廣義Kelvin模型的蠕變曲線呈線性變化，而NGKH模型的穩(wěn)態(tài)蠕變速率和蠕變曲線均呈非線性變化。因此，在不同應力狀態(tài)下，鹽巖的蠕變曲線形態(tài)可通過變化NGKH模型中a，A和B來調整，使NGKH模型更加適用于巖石的流變分析。

2 蠕變模型參數(shù)辨識

針對本文提出的NGKH蠕變模型，式(13)中需要給出的參數(shù)分別是K，G0，G1，a，η1，A，B和Q，其中A，B和Q已經(jīng)通過室內蠕變試驗曲線擬合分析確定。令A1=(σ1+2σ3)/(9K)+(σ1-σ3)/(3G0)，再利用K，G0與彈性模量E0和泊松比λ0之間的關系求得E0，進而反求出K和G0，此時，需要擬合確定的參數(shù)變?yōu)锳1，G1，a和η1。根據(jù)圖3中不同應力狀態(tài)下鹽巖室內三軸蠕變試驗結果，運用Origin分析軟件的Levenberg-Marquardt算法[23]，通過自定義函數(shù)模型，對A1，G1，a和η1進行辨識。表2、表3和表4所示分別為NGKH模型、非定常廣義Kelvin模型和廣義Kelvin模型的蠕變參數(shù)辨識結果。將擬合得到的各模型蠕變曲線與室內試驗結果進行對比，如圖7所示。其中圖7(a)～7(d)所示分別為不同偏應力下的蠕變曲線對比，圖7(a)、圖7(e)和圖7(f)所示分別為不同溫度下的蠕變曲線對比。從表2～4可見：在相同應力狀態(tài)下，不同蠕變模型的K、G0和G1較為接近，而a和η1卻存在明顯差異。從圖7(a)和圖7(b)可知：a和η1組成的非定常黏壺對衰減蠕變階段起著控制作用，而Heard體的加入直接影響了穩(wěn)態(tài)蠕變階段的蠕變曲線走向，同時反映出NGKH模型的靈活性和適用性。

表2 NGKH模型蠕變參數(shù)辨識Table2 Creep parameters identification for NGKH model

表3 非定常廣義Kelvin模型參數(shù)辨識Table3 Creep parameters identification for non-stationary generalized Kelvin model

表4 廣義Kelvin模型參數(shù)辨識Table4 Creep parameters identification for generalized Kelvin model

從圖7可知：與非定常廣義Kelvin模型及常規(guī)廣義Kelvin模型相比，由于引入了非定常黏壺及Heard蠕變體，本文提出的NGKH模型的擬合效果相對較好，理論計算結果與室內試驗結果吻合度明顯更高，能夠較好地對不同應力狀態(tài)下鹽巖的瞬時應變、瞬態(tài)蠕變階段和穩(wěn)態(tài)蠕變階段進行描述，更能反映鹽巖的非線性蠕變特征，由此驗證了該模型的合理性和可靠性。

3 蠕變模型二次開發(fā)

FLAC3D軟件內置的蠕變模型均有特定的使用范圍，在一定程度上影響了FLAC3D的廣泛應用[24]。因此，部分學者利用其提供的接口程序對新的蠕變模型進行了二次開發(fā)，例如熊良宵等[25-27]就基于FLAC3D分別對改進Burgers模型和河海模型進行了二次開發(fā)。為了能夠將本文提出的NGKH蠕變模型程序化，首先將其本構方程改寫成三維差分形式。

3.1 NGKH模型的有限差分形式

3.1.1 模型各部分的偏應力與偏應變關系式

由圖4和式(13)可知：NGKH模型各部分之間為串聯(lián)關系，因此應力相等，應變相加，則有[24,28]

對于串聯(lián)彈簧胡克體，偏應力與偏應變速率的關系為

對于非定常Kelvin體，偏應力Sij和偏應變有如下關系[28]：

式中：為非定常Kelvin體的黏滯系數(shù)，且

對于串聯(lián)Heard體，有

引入應力強度q，且q=在FLAC3D中可通過相應指針讀取應力張量的各分量，然后求出q。由經(jīng)典彈塑性力學可知因此，對于Heard體，偏應力Sij與偏應變速率之間的關系為

3.1.2 模型各部分的增量形式

為了利用FLAC3D軟件進行二次開發(fā)，將NGKH本構方程的上述內容改寫成有限差分形式。將式(15)寫成增量形式：

式(16)的增量形式可寫為

將式(17)寫成增量形式：

圖7 NGKH模型、非定常廣義Kelvin模型和廣義Kelvin模型擬合曲線對比Fig.7 Comparison of fitting curves of NGKH model,non-stationary generalized Kelvin model and generalized Kelvin model

將式(22)整理后，得到非定常Kelvin體第i步偏應變更新公式為

將式(23)進一步整理后可得：

式中：γ為蠕變乘子，隨應力和蠕變應變的變化而變化。

由此式(19)的增量形式可寫為

將式(21)、式(24)和式(27)代入式(20)，有

整理可得NGKH模型第i步應力更新公式為

3.2 二次開發(fā)程序的驗證

基于C++語言，利用Visual Studio2010平臺將上述NGKH蠕變模型本構方程的差分形式進行編譯寫入動態(tài)鏈接庫.dll文件中，利用FLAC3D的UDM接口程序進行二次開發(fā)。圖8所示為NGKH模型的二次開發(fā)流程。

圖8 自定義NGKH本構模型二次開發(fā)流程Fig.8 Secondary development process of self-defined NGKH constitutive model

為了驗證開發(fā)NGKH模型程序的正確性，建立與室內試驗規(guī)格相同的直徑為50mm、高度為100mm的鹽巖圓柱體數(shù)值模型，共劃分1 440個單元和1 595個節(jié)點，如圖9所示。模型底面施加法向約束，頂面分別施加30，35，40和45MPa的均布載荷，側面施加15MPa圍壓，分別進行三軸蠕變壓縮試驗。

3.2.1 黏彈性特征驗證

圖9 鹽巖的有限差分數(shù)值試樣Fig.9 Finite difference numerical sample for salt rock

由式(11)可知：當a=0時，非定常廣義Kelvin模型退化為廣義Kelvin模型，因此，在NGKH模型中將a和A同時設置為0，此時NGKH模型即退化為廣義Kelvin模型。在FLAC3D程序內置的常規(guī)Burgers模型中，對Maxwell體串聯(lián)牛頓黏壺的黏性系數(shù)不進行賦值，即退化為廣義Kelvin模型。退化后的NGKH模型和退化后的Burgers模型均屬于三元件模型，按照表4所示的廣義Kelvin參數(shù)分別對這2個退化模型進行賦值，計算不同軸向載荷下的鹽巖試樣蠕變曲線。圖10(a)～10(d)所示分別為偏應力25MPa時由退化NGKH模型計算得到的蠕變10，20，30和40 h后的鹽巖數(shù)值試樣位移分布云圖，圖11所示為不同偏應力下退化NGKH模型與退化Burgers模型的蠕變曲線對比圖。

由圖10可以看出：鹽巖試樣的頂部軸向位移最大，往底部方向逐漸減小；隨著蠕變時間的增加，鹽巖試樣的最大軸向位移隨之增大，但在20 h后軸向應變增速開始放緩。從圖11可知：在非定常參數(shù)a及Heard體不起作用的情況下，退化NGKH模型與退化Burgers模型的瞬時應變一致，曲線吻合良好，表明開發(fā)的NGKH模型程序進行的黏彈性計算結果是可靠的。

3.2.2 非線性特征驗證

按照表2中蠕變參數(shù)對NGKH蠕變模型賦值并進行三軸蠕變數(shù)值計算，圖12所示為開發(fā)的NGKH模型程序模擬結果與室內試驗結果的對比。其中圖12(a)所示為不同偏應力下數(shù)值模擬結果與室內試驗結果對比，圖12(b)所示為不同溫度下數(shù)值模擬結果與室內試驗結果對比。從圖12可以看出開發(fā)NGKH模型的程序模擬結果與室內試驗結果吻合良好，證明該模型的非線性蠕變計算結果是可靠的。

圖10 不同蠕變時刻退化Burgers模型的軸向位移分布Fig.10 Axial deformation distributions of degenerated Burgers model at different creep time

圖11 退化NGKH模型和退化Burgers模型的黏彈性蠕變曲線對比Fig.11 Comparison of viscoelastic creep curves between degraded NGKH model and degraded Burgers model

圖12 NGKH模型的非線性計算結果與試驗結果的對比Fig.12 Comparison between nonlinear computation results from NGKH model and creep test results

上述分析驗證了本文提出的將包含非定常黏壺的非定常廣義Kelvin模型與穩(wěn)態(tài)蠕變Heard體進行串聯(lián)，最終構建NGKH蠕變模型以描述鹽巖的非線性蠕變特征思路的正確性以及該蠕變模型在FLAC3D中二次開發(fā)的正確性。開發(fā)得到的FLAC3D蠕變本構模型為鹽巖地層井壁圍巖穩(wěn)定性分析及套損問題研究提供了參考。

4 結論

1)將牛頓黏壺的黏滯系數(shù)視為時間的函數(shù)，采用非定常黏壺替換常規(guī)廣義Kelvin模型中的線性牛頓黏壺，建立了非定常廣義Kelvin模型，將其與非線性Heard體進行串聯(lián)，組成能夠描述鹽巖瞬時蠕變階段和穩(wěn)態(tài)蠕變階段的四元件非線性NGKH蠕變模型。

2)基于鹽巖的室內三軸蠕變試驗結果，對NGKH模型的蠕變參數(shù)進行辨識。擬合結果與室內試驗曲線吻合良好，表明與廣義Kelvin模型及非定常廣義Kelvin模型相比，本文提出的NGKH蠕變模型更加適合描述川東北油氣田嘉陵江組鹽巖的非線性蠕變特征。

3)基于FLAC3D二次開發(fā)NGKH模型程序的黏彈性計算結果與非線性計算結果都是可靠的，從而證實了開發(fā)NGKH模型的正確性與合理性。