韓偉民，閆怡飛，閆相禎

(1.中國(guó)石油大學(xué)(華東)儲(chǔ)運(yùn)與建筑工程學(xué)院，山東青島，266580；2.中國(guó)石油大學(xué)(華東)機(jī)電工程學(xué)院，山東青島，266580)

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

1 鹽巖蠕變本構(gòu)模型

1.1 非定常廣義Kelvin模型的構(gòu)建

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

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

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

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

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

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

對(duì)式(2)求導(dǎo)可得

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

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

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

式中：˙為蠕變速率。

分離變量求定積分，并由ε=J(t)σ可得其蠕變?nèi)崃縅K為

由式(1)和式(5)可得非定常廣義Kelvin模型的蠕變?nèi)崃繛?/p>

則非定常廣義Kelvin模型的一維本構(gòu)方程為

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

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

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

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

1.2 Heard模型的構(gòu)建

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

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

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

圖3 鹽巖室內(nèi)試驗(yàn)蠕變曲線Fig.3 Creep curves of salt rock

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

1.3 NGKH模型的構(gòu)建

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

FLAC3D軟件內(nèi)置的蠕變模型均有特定的使用范圍，在一定程度上影響了FLAC3D的廣泛應(yīng)用[24]。因此，部分學(xué)者利用其提供的接口程序?qū)π碌娜渥兡Ｐ瓦M(jìn)行了二次開發(fā)，例如熊良宵等[25-27]就基于FLAC3D分別對(duì)改進(jìn)Burgers模型和河海模型進(jìn)行了二次開發(fā)。為了能夠?qū)⒈疚奶岢龅腘GKH蠕變模型程序化，首先將其本構(gòu)方程改寫成三維差分形式。

3.1 NGKH模型的有限差分形式

3.1.1 模型各部分的偏應(yīng)力與偏應(yīng)變關(guān)系式

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

對(duì)于串聯(lián)彈簧胡克體，偏應(yīng)力與偏應(yīng)變速率的關(guān)系為

對(duì)于非定常Kelvin體，偏應(yīng)力Sij和偏應(yīng)變有如下關(guān)系[28]：

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

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

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

3.1.2 模型各部分的增量形式

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

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

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

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

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

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

式中：γ為蠕變乘子，隨應(yīng)力和蠕變應(yīng)變的變化而變化。

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

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

整理可得NGKH模型第i步應(yīng)力更新公式為

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

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

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

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

3.2.1 黏彈性特征驗(yàn)證

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

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

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

3.2.2 非線性特征驗(yàn)證

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

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

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

圖12 NGKH模型的非線性計(jì)算結(jié)果與試驗(yàn)結(jié)果的對(duì)比Fig.12 Comparison between nonlinear computation results from NGKH model and creep test results

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

4 結(jié)論

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

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

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