隋義勇, 林堂茂, 劉 翔, 董長(zhǎng)銀, 程 威, 張廣明, 廖正毅

(1.中國(guó)石油大學(xué)(華東)石油工程學(xué)院,山東青島 266580; 2.中國(guó)石油勘探開發(fā)研究院,北京 100083)

中國(guó)越來越多的枯竭型油氣藏儲(chǔ)氣庫(kù)投入建設(shè)和運(yùn)行,而對(duì)儲(chǔ)氣庫(kù)建設(shè)仍處于初級(jí)技術(shù)階段[1-4]。針對(duì)循環(huán)應(yīng)力對(duì)儲(chǔ)氣庫(kù)巖石微觀結(jié)構(gòu)和力學(xué)性質(zhì)的影響,馬新華等[5]通過室內(nèi)巖心三軸加卸載試驗(yàn)研究巖石變形破壞特征,分析了累積塑性變形的變化規(guī)律。鄭得文等[6]通過室內(nèi)試驗(yàn)研究?jī)?chǔ)氣庫(kù)生產(chǎn)過程中的剪切破壞、拉伸破壞及斷層剪切滑移失穩(wěn)風(fēng)險(xiǎn)等問題。孫巖等[7]根據(jù)流體物性參數(shù)和相應(yīng)生產(chǎn)動(dòng)態(tài)數(shù)據(jù),優(yōu)化了多周期注采水侵量計(jì)算模型。Wang等[8]利用有限元法研究?jī)?chǔ)氣庫(kù)循環(huán)注采過程中的疲勞損傷、彈塑性變形及孔隙度變化。Guo等[9]研究循環(huán)應(yīng)力作用下黏土彈性模量變化規(guī)律及永久累積變形特征。Wang等[10]研究循環(huán)注采過程中產(chǎn)生的微地震與應(yīng)力歷史對(duì)巖石內(nèi)部損傷及破壞的影響。Sun等[11]建立了一種儲(chǔ)氣庫(kù)綜合地質(zhì)力學(xué)模型,評(píng)估儲(chǔ)氣庫(kù)長(zhǎng)期運(yùn)行過程中蓋層完整性破壞和故障泄漏風(fēng)險(xiǎn)。陳金平等[12]使用有限元法研究?jī)?chǔ)氣庫(kù)地應(yīng)力場(chǎng)和內(nèi)壓作用下,不同長(zhǎng)短軸比橢球形儲(chǔ)氣庫(kù)腔體極限壓力。Dai等[13]使用離散元法研究花崗巖在加載卸載過程中應(yīng)變和損傷演變情況。前人研究主要集中在宏觀層次,缺少對(duì)微觀結(jié)構(gòu)變化和機(jī)制的研究,巖石微觀結(jié)構(gòu)變化是影響巖石宏觀力學(xué)性質(zhì)變化的根源[14-16]。筆者采用顆粒離散元數(shù)值模擬方法研究循環(huán)應(yīng)力對(duì)儲(chǔ)氣庫(kù)巖石微觀結(jié)構(gòu)及力學(xué)性質(zhì)的影響。

1 真三軸顆粒離散元模型

1.1 建立離散元模型

離散元法將巖體看作是各離散單元的黏結(jié)集合體,假設(shè)兩離散單元之間存在若干彈簧連接,通過兩者重疊量以及彈簧剛度確定接觸力,離散單元的運(yùn)動(dòng)規(guī)律符合宏觀牛頓定律,離散單元間碰撞符合胡克定律以及彈塑性力學(xué)規(guī)律[17-21],離散元法將非連續(xù)介質(zhì)與連續(xù)介質(zhì)力學(xué)問題有機(jī)結(jié)合起來。

首先,根據(jù)國(guó)際巖石力學(xué)學(xué)會(huì)對(duì)真三軸巖心尺寸要求,確定數(shù)值模型尺寸為2.5 cm×2.5 cm×10 cm,再根據(jù)巖心室內(nèi)試驗(yàn)測(cè)得巖心粒度組成及力學(xué)性質(zhì),數(shù)值模型粒度組成既要確保顆粒體系能表征巖心的力學(xué)特性,又要確保模型中顆粒數(shù)量適中提高計(jì)算效率,最終確定模型顆粒半徑為0.14～3.62 mm,生成了一個(gè)包含27 000多顆粒的純顆粒模型;然后,考慮到枯竭型油氣藏儲(chǔ)氣庫(kù)建庫(kù)之前經(jīng)過多年的油氣開采,儲(chǔ)層巖石中存在大量微裂縫,通過PFC3D中縫網(wǎng)生成功能建立隨機(jī)原生裂縫模型;最后,將純顆粒模型和裂縫模型疊加,得到離散元真三軸巖心數(shù)值模型,如圖1所示。

圖1 巖心離散元模型建立過程Fig.1 Modeling process of core discrete element model

1.2 模型接觸參數(shù)及應(yīng)力加載參數(shù)

根據(jù)巖體的非連續(xù)特性及顆粒間力學(xué)接觸特性,選用PFC3D中自帶的線性平行黏結(jié)接觸模型模擬巖體顆粒間接觸摩擦特性與接觸黏結(jié)特性。選用PFC3D中自帶的光滑節(jié)理接觸模型模擬巖體原生裂縫內(nèi)顆粒間的接觸特性。根據(jù)巖心的室內(nèi)力學(xué)試驗(yàn)測(cè)量結(jié)果標(biāo)定模型接觸參數(shù),最終確定的模型接觸參數(shù)如表1所示。室內(nèi)試驗(yàn)曲線與數(shù)值模擬曲線對(duì)比如圖2所示。由圖2可以看出,數(shù)值模擬獲得的應(yīng)力-應(yīng)變曲線與巖心室內(nèi)試驗(yàn)測(cè)得的應(yīng)力-應(yīng)變曲線吻合度較好,建立的巖心真三軸離散元模型以及標(biāo)定的模型參數(shù)能夠表征儲(chǔ)氣庫(kù)儲(chǔ)層巖石力學(xué)特性,可以用來進(jìn)行后續(xù)的多周期等幅值軸向循環(huán)應(yīng)力加載數(shù)值模擬研究。

表1 巖心顆粒接觸模型參數(shù)

根據(jù)儲(chǔ)氣庫(kù)儲(chǔ)層巖石的原位地應(yīng)力狀態(tài)確定數(shù)值模型的初始應(yīng)力狀態(tài),已知新疆某儲(chǔ)氣庫(kù)儲(chǔ)層3個(gè)方向的原位地應(yīng)力梯度,結(jié)合儲(chǔ)層埋深,確定儲(chǔ)氣庫(kù)儲(chǔ)層3個(gè)方向的原位地應(yīng)力狀態(tài)。利用伺服原理控制側(cè)向無摩擦剛性墻體對(duì)數(shù)值模型施加水平方向的圍壓,其中水平最小主應(yīng)力方向圍壓為30 MPa,水平最大主應(yīng)力方向圍壓為48 MPa。數(shù)值模型上、下墻體為加載板,通過控制上、下墻體緩慢勻速往復(fù)運(yùn)動(dòng)對(duì)數(shù)值模型施加軸向循環(huán)加載應(yīng)力,結(jié)合儲(chǔ)氣庫(kù)生產(chǎn)上、下限庫(kù)內(nèi)壓力和儲(chǔ)層上覆地層壓力,確定數(shù)值模型軸向循環(huán)加載應(yīng)力上限為65.23 MPa,軸向循環(huán)加載應(yīng)力下限為49.23 MPa,儲(chǔ)氣庫(kù)通常要求安全平穩(wěn)生產(chǎn)50～60 a,因此決定對(duì)數(shù)值模型進(jìn)行等應(yīng)力幅值軸向循環(huán)加載60個(gè)周期,研究循環(huán)應(yīng)力加載對(duì)儲(chǔ)氣庫(kù)巖石微觀結(jié)構(gòu)和力學(xué)性質(zhì)的影響。循環(huán)應(yīng)力加載示意圖如圖3所示。

圖2 室內(nèi)試驗(yàn)曲線與數(shù)值模擬曲線對(duì)比Fig.2 Comparison of laboratory experimental curves and numerical simulation curves

圖3 循環(huán)應(yīng)力加載示意圖Fig.3 Cyclic stress loading diagram

2 循環(huán)應(yīng)力對(duì)儲(chǔ)層巖石內(nèi)部結(jié)構(gòu)的影響

2.1 循環(huán)應(yīng)力對(duì)模型內(nèi)部黏結(jié)狀態(tài)的影響

不同加載周期模型內(nèi)生成的無黏結(jié)接觸分布如圖4所示,紅色部分顯示的是循環(huán)應(yīng)力加載1、30、60個(gè)周期時(shí)模型內(nèi)顆粒間黏結(jié)接觸斷裂后生成的無黏結(jié)接觸分布情況。由圖4可以看出,隨著循環(huán)應(yīng)力加載周期增加,模型內(nèi)顆粒間生成的無黏結(jié)接觸數(shù)量逐漸增加,對(duì)比原始裂縫分布與無黏結(jié)接觸分布位置,發(fā)現(xiàn)循環(huán)應(yīng)力加載過程中模型內(nèi)顆粒間生成的無黏結(jié)接觸主要集中在原始裂縫分布密集的區(qū)域。

多周期循環(huán)應(yīng)力加載對(duì)無黏結(jié)接觸數(shù)量影響如圖5所示。由圖5可以看出,循環(huán)應(yīng)力加載過程中模型內(nèi)顆粒間生成的無黏結(jié)接觸數(shù)量增加速率呈遞減趨勢(shì),循環(huán)應(yīng)力加載前30個(gè)周期模型內(nèi)顆粒間無黏結(jié)接觸數(shù)量增加速度較快,循環(huán)應(yīng)力加載后30個(gè)周期模型內(nèi)顆粒間無黏結(jié)接觸數(shù)量增加速度減緩。

圖4 不同加載周期模型內(nèi)生成的無黏結(jié)接觸分布Fig.4 Non-bonded contacts distribution in model with different loading cycles

圖5 多周期循環(huán)應(yīng)力加載對(duì)無黏結(jié)接觸數(shù)量影響Fig.5 Effect of cyclic stress loading on number of non-bonded contacts

循環(huán)應(yīng)力加載對(duì)模型內(nèi)顆粒間黏結(jié)狀態(tài)的影響與模型內(nèi)顆粒間接觸強(qiáng)度分布和裂縫分布有關(guān)。初始時(shí)模型內(nèi)顆粒間黏結(jié)接觸強(qiáng)度服從Weibull分布,模型內(nèi)各區(qū)域顆粒間黏結(jié)接觸強(qiáng)度分布不均勻,循環(huán)應(yīng)力加載時(shí)應(yīng)力上限雖未達(dá)到模型破壞強(qiáng)度,但已經(jīng)超過模型內(nèi)顆粒間弱黏結(jié)接觸的破壞強(qiáng)度,導(dǎo)致模型內(nèi)顆粒間黏結(jié)接觸發(fā)生斷裂。模型內(nèi)顆粒在加載外力、不平衡力和顆粒間接觸力的共同作用下發(fā)生分離運(yùn)移并產(chǎn)生新接觸,力和位移共同作用導(dǎo)致模型內(nèi)顆粒間強(qiáng)黏結(jié)接觸逐漸弱化并發(fā)生斷裂,因此隨著循環(huán)應(yīng)力加載周期增加,模型內(nèi)顆粒間產(chǎn)生的無黏結(jié)接觸數(shù)量逐漸增加。另外,模型內(nèi)原始裂縫分布密集區(qū)域存在軟弱界面,在加載外力作用下應(yīng)力會(huì)在界面處集中導(dǎo)致顆粒沿裂縫面發(fā)生滑移,因此模型內(nèi)顆粒間產(chǎn)生的無黏結(jié)接觸集中在原始裂縫分布密集區(qū)域。循環(huán)應(yīng)力加載前30個(gè)周期模型內(nèi)顆粒的位移量較大,顆粒間黏結(jié)接觸斷裂產(chǎn)生的無黏結(jié)接觸數(shù)量較多。循環(huán)應(yīng)力加載后30個(gè)周期模型內(nèi)顆粒逐漸被壓密而分布較均勻,模型內(nèi)部應(yīng)力分布逐漸趨向均勻,模型內(nèi)顆粒位移量減小,顆粒間黏結(jié)接觸斷裂產(chǎn)生的無黏結(jié)接觸數(shù)量減少,因此隨循環(huán)應(yīng)力加載周期增加模型內(nèi)顆粒間產(chǎn)生的無黏結(jié)接觸數(shù)量增加速率遞減。

2.2 循環(huán)應(yīng)力對(duì)模型內(nèi)部微裂縫發(fā)育的影響

圖6為不同應(yīng)力加載周期時(shí)模型內(nèi)生成的微裂縫分布與原始裂縫分布對(duì)比,其中紅色部分代表循環(huán)應(yīng)力加載后模型內(nèi)生成的微裂縫分布,綠色部分代表原始裂縫分布。由圖6可以看出,隨著循環(huán)應(yīng)力加載周期增加,模型內(nèi)生成的微裂縫數(shù)量和微裂縫分布范圍逐漸增加,且新生成的微裂縫分布受原始裂縫分布的控制明顯,新生成的微裂縫主要分布在原始裂縫分布密集區(qū)域,在原始裂縫的邊緣和尖端萌生后逐漸延展并相互貫穿連通。

圖6 不同加載周期模型內(nèi)微裂縫分布與原始裂縫分布對(duì)比Fig.6 Comparison of micro-cracks distribution with original crack distribution with different loading cycles

循環(huán)應(yīng)力加載對(duì)模型內(nèi)微裂縫發(fā)育的影響如圖7所示。由圖7可以看出,隨著循環(huán)應(yīng)力加載周期增加,模型內(nèi)生成的微裂縫數(shù)量增加速率呈遞減趨勢(shì),循環(huán)應(yīng)力加載前30個(gè)周期模型內(nèi)生成的微裂縫數(shù)量增加速度較快,循環(huán)應(yīng)力加載后30個(gè)周期模型內(nèi)生成的微裂縫數(shù)量增加速度減緩。

圖7 循環(huán)應(yīng)力加載對(duì)模型內(nèi)微裂縫發(fā)育的影響Fig.7 Effect of cyclic stress loading on development of micro-cracks in the model

循環(huán)應(yīng)力加載過程中,模型內(nèi)原始裂縫在張應(yīng)力的作用下克服模型內(nèi)顆粒間摩擦力往復(fù)開合,導(dǎo)致原始裂縫分布區(qū)域及其周邊區(qū)域的顆粒承受往復(fù)的張應(yīng)力以及剪應(yīng)力作用,因此原始裂縫分布區(qū)域及其周邊區(qū)域更容易產(chǎn)生微裂縫。同時(shí)在加載外力、不平衡力和顆粒間接觸力的共同作用下,原始裂縫周圍顆粒易沿原始裂縫面產(chǎn)生較大位移量的滑移。這有利于微裂縫的延展和貫穿連通,因此模型內(nèi)生成的微裂縫主要集中在原始裂縫分布區(qū)域及其周邊區(qū)域。循環(huán)應(yīng)力加載前30個(gè)周期,模型內(nèi)新生成的微裂縫主要受原始裂縫控制,因此生成的微裂縫數(shù)量增加速率較快。循環(huán)應(yīng)力加載后30個(gè)周期,模型內(nèi)顆粒體系重新分布導(dǎo)致模型內(nèi)應(yīng)力趨向均衡,模型內(nèi)新生成的微裂縫受原始裂縫控制減弱,模型內(nèi)生成的微裂縫數(shù)量增加速率減緩,但模型內(nèi)生成的微裂縫分布更加廣泛。

2.3 循環(huán)應(yīng)力對(duì)模型孔隙度的影響

圖8為模型歸一化孔隙度φn/φ1(φn和φ1分別為第n和第1循環(huán)應(yīng)力加載周期模型孔隙度，%)隨循環(huán)應(yīng)力加載周期增加的變化。由圖8可以看出,隨著循環(huán)應(yīng)力加載周期增加,軸向加載應(yīng)力處于循環(huán)應(yīng)力下限時(shí)模型歸一化孔隙度逐漸降低,且模型歸一化孔隙度降低速率呈遞減趨勢(shì),循環(huán)應(yīng)力加載前30個(gè)周期模型歸一化孔隙度降低速率較快,循環(huán)應(yīng)力加載后30個(gè)周期模型歸一化孔隙度降低速率減緩。其中前10個(gè)循環(huán)應(yīng)力加載周期內(nèi)歸一化孔隙度降低量約占整個(gè)循環(huán)應(yīng)力加載過程中歸一化孔隙度總降低量的30%。

圖8 循環(huán)應(yīng)力加載對(duì)模型孔隙度的影響Fig.8 Effect of cyclic stress loading on model porosity

循環(huán)應(yīng)力加載過程中模型的變形主要有可恢復(fù)的彈性變形和不可逆的塑性變形,模型孔隙度降低主要是不可逆塑性變形引起的。隨著循環(huán)應(yīng)力加載周期增加,模型內(nèi)生成的無黏結(jié)接觸數(shù)量和微裂縫數(shù)量逐漸增加,模型內(nèi)產(chǎn)生大量游離顆粒。隨著循環(huán)應(yīng)力加載的進(jìn)行模型內(nèi)游離顆粒運(yùn)移至大孔洞和裂縫中,導(dǎo)致模型中大孔洞逐漸被填充,裂縫逐漸被壓密甚至閉合,模型在相同應(yīng)力水平下顆粒分布更加均勻,模型的孔隙度逐漸減小。隨著循環(huán)應(yīng)力加載周期增加,模型內(nèi)生成的無黏結(jié)接觸數(shù)量和微裂縫數(shù)量增加速率逐漸減緩,因此模型內(nèi)不可逆塑性變形量逐漸減小,模型孔隙度降低量也隨之減小,模型孔隙度降低速率逐漸減緩。

3 循環(huán)應(yīng)力對(duì)儲(chǔ)層巖石力學(xué)性質(zhì)的影響

3.1 循環(huán)應(yīng)力對(duì)儲(chǔ)層巖石應(yīng)變的影響

圖9為模型軸向和側(cè)向應(yīng)變隨著環(huán)應(yīng)力加載周期的變化。由圖9可以看出,隨著循環(huán)應(yīng)力加載周期增加,在軸向加載應(yīng)力達(dá)到循環(huán)應(yīng)力上限時(shí),模型軸向應(yīng)變與最小主應(yīng)力方向側(cè)向應(yīng)變均逐漸增大,但兩者應(yīng)變?cè)黾铀俾什町惷黠@,最小主應(yīng)力方向側(cè)向應(yīng)變?cè)黾铀俾始s是軸向應(yīng)變?cè)黾铀俾实膬杀?。這表明在多周期等幅值軸向循環(huán)應(yīng)力加載的疲勞損傷過程中,模型側(cè)向膨脹的影響遠(yuǎn)大于軸向壓縮的影響,模型側(cè)向膨脹占主導(dǎo)作用。

圖9 循環(huán)應(yīng)力加載對(duì)軸向應(yīng)變和側(cè)向應(yīng)變影響Fig.9 Effects of cycle stress loading on axial strains and lateral strains

3.2 循環(huán)應(yīng)力對(duì)儲(chǔ)層巖石彈性模量的影響

圖10為模型彈性模量隨循環(huán)應(yīng)力加載周期的變化。

圖10 循環(huán)應(yīng)力加載對(duì)彈性模量影響Fig.10 Effect of cyclic stress loading on elastic modulus

由圖10可以看出,隨著循環(huán)應(yīng)力加載周期增加,模型彈性模量逐漸減小,模型彈性模量降低速率呈遞減趨勢(shì),循環(huán)應(yīng)力加載前30個(gè)周期模型彈性模量降低速率較快,后30個(gè)周期模型彈性模量降低速率減緩,循環(huán)應(yīng)力加載60個(gè)周期模型彈性模量較原始彈性模量降低0.33 GPa,降幅為22%。

3.3 循環(huán)應(yīng)力對(duì)儲(chǔ)層巖石泊松比的影響

圖11為模型泊松比隨循環(huán)應(yīng)力加載周期的變化。由圖11可以看出,隨著循環(huán)應(yīng)力加載周期增加,模型泊松比逐漸增大,模型泊松比增加速率呈遞減趨勢(shì)。循環(huán)應(yīng)力加載前30個(gè)周期模型泊松比增加速率較快,后30個(gè)周期模型泊松比增加速率減緩,循環(huán)應(yīng)力加載60個(gè)周期模型泊松比較原始泊松比增加0.157,增幅為78%。

圖11 循環(huán)應(yīng)力加載對(duì)泊松比影響Fig.11 Effect of cycle stress loading on Poissons ratio

3.4 循環(huán)應(yīng)力對(duì)儲(chǔ)層巖石損傷量的影響

圖12為模型損傷量隨循環(huán)應(yīng)力加載周期的變化。由圖12可以看出,隨著循環(huán)應(yīng)力加載周期增加,模型損傷量逐漸增大,則模型抵抗變形和破壞的能力降低,模型損傷量增加速率呈遞減趨勢(shì)。循環(huán)應(yīng)力加載前30個(gè)周期模型損傷量增加速率較快,后30個(gè)周期模型損傷量增加速率減緩,循環(huán)應(yīng)力加載60個(gè)周期模型損傷量增加22%。

圖12 循環(huán)應(yīng)力加載對(duì)模型損傷量影響Fig.12 Effect of cycle stress loading on model damage amount

循環(huán)應(yīng)力加載過程中模型內(nèi)新生成的大量無黏結(jié)接觸無法抵抗張應(yīng)力,只能抵抗壓應(yīng)力,同時(shí)模型內(nèi)新生成的微裂縫數(shù)量和模型損傷量逐漸增加,模型整體強(qiáng)度降低,這導(dǎo)致模型抵抗變形的能力降低。因此隨著循環(huán)應(yīng)力加載周期增加,相同應(yīng)力水平下模型側(cè)向應(yīng)變及軸向應(yīng)變量均逐漸增加,模型彈性模量逐漸降低。根據(jù)巖石力學(xué)理論可知,巖石抵抗張應(yīng)力的能力遠(yuǎn)低于抵抗壓應(yīng)力的能力,而且模型內(nèi)新生成的無法抵抗張應(yīng)力的無黏結(jié)接觸數(shù)量逐漸增多。因此隨著循環(huán)應(yīng)力加載周期增加,相同應(yīng)力水平下模型側(cè)向應(yīng)變?cè)黾铀俾蚀笥谳S向應(yīng)變?cè)黾铀俾?循環(huán)應(yīng)力加載過程中模型側(cè)向變形占主導(dǎo)地位,模型泊松比逐漸增加。在循環(huán)應(yīng)力加載過程中,模型內(nèi)新生成的無黏結(jié)接觸數(shù)量和微裂縫數(shù)量增加速率呈遞減趨勢(shì),因此模型損傷量增加速率同樣呈遞減趨勢(shì)。

4 結(jié) 論

(1)隨著循環(huán)應(yīng)力加載周期增加,新生成的無黏結(jié)接觸數(shù)量和微裂縫數(shù)量均逐漸增加,且新生成的無黏結(jié)接觸分布和微裂縫分布均受原始裂縫分布控制明顯,而孔隙度逐漸降低,且無黏結(jié)接觸數(shù)量、微裂縫數(shù)量和孔隙度變化速率均呈遞減趨勢(shì)。

(2)隨著循環(huán)應(yīng)力加載周期增加,相同應(yīng)力水平下側(cè)向應(yīng)變和軸向應(yīng)變均增大,但側(cè)向應(yīng)變占主導(dǎo)地位,彈性模量逐漸減小,而泊松比和損傷量逐漸增加,且彈性模量、泊松比和損傷量變化速率均呈遞減趨勢(shì)。