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

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

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

1 真三軸顆粒離散元模型

1.1 建立離散元模型

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

首先,根據(jù)國際巖石力學(xué)學(xué)會對真三軸巖心尺寸要求,確定數(shù)值模型尺寸為2.5 cm×2.5 cm×10 cm,再根據(jù)巖心室內(nèi)試驗測得巖心粒度組成及力學(xué)性質(zhì),數(shù)值模型粒度組成既要確保顆粒體系能表征巖心的力學(xué)特性,又要確保模型中顆粒數(shù)量適中提高計算效率,最終確定模型顆粒半徑為0.14～3.62 mm,生成了一個包含27 000多顆粒的純顆粒模型;然后,考慮到枯竭型油氣藏儲氣庫建庫之前經(jīng)過多年的油氣開采,儲層巖石中存在大量微裂縫,通過PFC3D中縫網(wǎng)生成功能建立隨機原生裂縫模型;最后,將純顆粒模型和裂縫模型疊加,得到離散元真三軸巖心數(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é)試驗測量結(jié)果標(biāo)定模型接觸參數(shù),最終確定的模型接觸參數(shù)如表1所示。室內(nèi)試驗曲線與數(shù)值模擬曲線對比如圖2所示。由圖2可以看出,數(shù)值模擬獲得的應(yīng)力-應(yīng)變曲線與巖心室內(nèi)試驗測得的應(yīng)力-應(yīng)變曲線吻合度較好,建立的巖心真三軸離散元模型以及標(biāo)定的模型參數(shù)能夠表征儲氣庫儲層巖石力學(xué)特性,可以用來進行后續(xù)的多周期等幅值軸向循環(huán)應(yīng)力加載數(shù)值模擬研究。

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

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

圖2 室內(nèi)試驗曲線與數(shù)值模擬曲線對比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)力對儲層巖石內(nèi)部結(jié)構(gòu)的影響

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

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

多周期循環(huán)應(yīng)力加載對無黏結(jié)接觸數(shù)量影響如圖5所示。由圖5可以看出,循環(huán)應(yīng)力加載過程中模型內(nèi)顆粒間生成的無黏結(jié)接觸數(shù)量增加速率呈遞減趨勢,循環(huán)應(yīng)力加載前30個周期模型內(nèi)顆粒間無黏結(jié)接觸數(shù)量增加速度較快,循環(huán)應(yīng)力加載后30個周期模型內(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)力加載對無黏結(jié)接觸數(shù)量影響Fig.5 Effect of cyclic stress loading on number of non-bonded contacts

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

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

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

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

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

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

圖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)力下限時模型歸一化孔隙度逐漸降低,且模型歸一化孔隙度降低速率呈遞減趨勢,循環(huán)應(yīng)力加載前30個周期模型歸一化孔隙度降低速率較快,循環(huán)應(yīng)力加載后30個周期模型歸一化孔隙度降低速率減緩。其中前10個循環(huán)應(yīng)力加載周期內(nèi)歸一化孔隙度降低量約占整個循環(huán)應(yīng)力加載過程中歸一化孔隙度總降低量的30%。

圖8 循環(huán)應(yīng)力加載對模型孔隙度的影響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)力加載的進行模型內(nèi)游離顆粒運移至大孔洞和裂縫中,導(dǎo)致模型中大孔洞逐漸被填充,裂縫逐漸被壓密甚至閉合,模型在相同應(yīng)力水平下顆粒分布更加均勻,模型的孔隙度逐漸減小。隨著循環(huán)應(yīng)力加載周期增加,模型內(nèi)生成的無黏結(jié)接觸數(shù)量和微裂縫數(shù)量增加速率逐漸減緩,因此模型內(nèi)不可逆塑性變形量逐漸減小,模型孔隙度降低量也隨之減小,模型孔隙度降低速率逐漸減緩。

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

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

圖9為模型軸向和側(cè)向應(yīng)變隨著環(huán)應(yīng)力加載周期的變化。由圖9可以看出,隨著循環(huán)應(yīng)力加載周期增加,在軸向加載應(yīng)力達到循環(huán)應(yīng)力上限時,模型軸向應(yīng)變與最小主應(yīng)力方向側(cè)向應(yīng)變均逐漸增大,但兩者應(yīng)變增加速率差異明顯,最小主應(yīng)力方向側(cè)向應(yīng)變增加速率約是軸向應(yīng)變增加速率的兩倍。這表明在多周期等幅值軸向循環(huán)應(yīng)力加載的疲勞損傷過程中,模型側(cè)向膨脹的影響遠大于軸向壓縮的影響,模型側(cè)向膨脹占主導(dǎo)作用。

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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

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