于恩毅 邸 元,2) 吳 輝 曹小朋 張慶福 張傳寶

* (北京大學工學院,北京 100871)

? (北京大學地球與空間科學學院,北京 100871)

** (中國石油化工股份有限公司勝利油田分公司勘探開發(fā)研究院,山東東營 257015)

引言

二氧化碳(CO2)被認為是主要的溫室氣體之一[1].由于大量排放,CO2在大氣中的含量不斷攀升,引起全球氣候變化,嚴重威脅人類的生存環(huán)境,減少CO2的排放已成為當今世界關(guān)注的焦點問題[2-3].在可預(yù)見的未來,化石能源仍將在能源消費體系中扮演重要角色,CO2捕集、利用和封存(CO2capture,utilization and storage,CCUS)是實現(xiàn)化石能源低碳化利用、減少CO2排放最直接和最關(guān)鍵的技術(shù)途徑[4],能夠為實現(xiàn)碳中和目標提供重要支撐[5].

CO2地質(zhì)封存是指通過工程技術(shù)手段將捕集的CO2注入深部地質(zhì)儲層,實現(xiàn)CO2與大氣長期隔絕的過程[5].如圖1 所示,封存點一般選擇在深度800～1000 m 以下,滲透性好且上方存在低滲蓋層的儲層[6].儲層類型主要包括咸水層、油氣藏和煤層等.在所有封存類型中,深部咸水層分布廣泛,封存容量占比約為98%[5].油氣藏也是較為理想的CO2封存場所,因為存在詳細的地質(zhì)勘探資料和較完備的地面設(shè)施等基礎(chǔ)條件[5],并且將CO2注入油氣藏能夠利用CO2驅(qū)替達到提高原油采收率的目的[7-9].

圖1 CO2 地質(zhì)封存的途徑[10]Fig.1 Approaches of CO2 geological storage[10]

CO2在儲層巖體中的封存過程較為復(fù)雜,受到地層特性和地球化學反應(yīng)等影響[11],涉及溫度場-滲流場-力學場-化學場(T-H-M-C)的耦合作用[12],如圖2 所示.CO2注入地層后會導致較大的流體壓力積聚,導致有效應(yīng)力場發(fā)生變化,應(yīng)力場的改變會影響巖體孔隙度、滲透率和毛管壓力等,從而影響CO2的注入和流動[13-14].注入的CO2為超臨界態(tài),溫度要明顯低于周圍地層溫度,二者的溫度差導致局部區(qū)域內(nèi)的溫度場發(fā)生變化,從而改變CO2流體的密度、黏度和溶解度等,影響其滲流特性,溫度變化引起的熱應(yīng)力也會直接改變巖體的應(yīng)力狀態(tài)[13].CO2易溶于水,進而與周圍巖石礦物發(fā)生化學反應(yīng),溶解巖石或鈣化沉淀,并可能與蓋層有機質(zhì)發(fā)生反應(yīng),導致蓋層滲透率和孔隙度等性質(zhì)發(fā)生變化[11,15].在多孔介質(zhì)中的化學反應(yīng)速率又受到溫度、壓力、滲流速度和CO2擴散速率等因素的影響[16].

圖2 CO2 地質(zhì)封存多場耦合原理[13]Fig.2 Multi-field coupled principle of CO2 geological storage[13]

CO2地質(zhì)封存的主要機制包括: (1)構(gòu)造封存,即低滲的密封蓋層限制CO2的遷移[17-18];(2)毛細封存,即CO2被孔隙結(jié)構(gòu)中的毛管力束縛[19];(3)吸附封存,即CO2吸附在黏土礦物表面[20];(4)溶解封存,即CO2與地層水或油混溶[21-22];(5) 礦化封存,即CO2與巖石發(fā)生化學反應(yīng)生成碳酸鹽礦化物[23-24].當CO2注入到地層后,首先圈閉在油氣藏或深層鹽水層中,CO2的運移由多相對流過程主導,主要由注入壓力和流體密度差異引起的浮力驅(qū)動[25].隨著時間的推移,保持自由相的CO2通過巖石中的滲流通道向上遷移,剩余的CO2移動非常緩慢,直到最終被殘留在孔隙中、溶解在地層水中,或作為碳酸鹽巖礦物沉淀[10,26].

CO2地質(zhì)封存有其自身的風險.若CO2垂向運移至近地表區(qū)域,會造成淺層飲用水污染[27-30]、逃逸至大氣環(huán)境[31]等風險.所以,在CO2注入作業(yè)前,必須對儲層-蓋層系統(tǒng)進行力學分析和穩(wěn)定性評估[32-33].蓋層一般是未受干擾、低滲透、厚且橫向分布廣泛的地層,常見為頁巖、泥巖和碳酸鹽巖等[34],具有高毛細管力和突破壓力,以防止注入流體泄漏[35].在CO2封存過程中,蓋層必須承受短期的過量注入壓力和長期的浮力驅(qū)動壓力[36].斷層和裂縫帶是蓋層巖體的軟弱帶,在CO2注入過程中,蓋層中的裂縫可能發(fā)生起裂、擴展,斷層可能被活化,發(fā)生滑移[37-38].CO2注入還可能會導致儲層中的巖體發(fā)生壓裂,裂縫可向上延伸至蓋層,甚至穿透蓋層[39].天然裂縫/斷層和壓裂裂縫形成了CO2逸出的潛在通道,對蓋層完整性構(gòu)成威脅[35,40-41].

儲層-蓋層系統(tǒng)多場耦合數(shù)值模擬是分析蓋層完整性和斷層穩(wěn)定性并進行預(yù)警控制的關(guān)鍵核心技術(shù)[25,42].CO2地質(zhì)封存過程的數(shù)值模擬需要考慮大時間尺度和大空間尺度下的多物理過程、多場耦合作用、儲層非均質(zhì)性和流體相態(tài)變化等因素[43-45],需要研發(fā)計算效率高、計算能力強的CO2地質(zhì)封存數(shù)值模擬程序.本文對CO2地質(zhì)封存儲層-蓋層系統(tǒng)多場耦合數(shù)值模擬方面的研究進行了全面綜述,介紹了CO2地質(zhì)封存過程中T-H-M-C 多場耦合問題的數(shù)學模型、耦合問題數(shù)值解法和封存風險分析等方面國內(nèi)外的研究進展,對當前我國CO2地質(zhì)封存數(shù)值模擬技術(shù)面臨的主要問題進行了總結(jié).

1 數(shù)學模型

CO2地質(zhì)封存多場耦合問題的數(shù)值模擬涉及的研究對象是巖石固體骨架與孔隙流體(CO2、水、油等)組成的多相多組分系統(tǒng),其數(shù)學模型主要包括基本控制方程和相平衡計算方法.基本控制方程包括應(yīng)力平衡方程、質(zhì)量守恒方程、能量守恒方程和化學反應(yīng)等.

1.1 應(yīng)力平衡方程

根據(jù)可變形多孔介質(zhì)多相流孔隙熱彈性理論,應(yīng)力平衡方程為[12,46]

式中,ω 和 λ 為拉梅常數(shù),u為位移,P為流體壓力,Ks為巖石骨架體積模量,Km為基質(zhì)顆粒體積模量,β為骨架的體積熱膨脹系數(shù),T為溫度.本文公式中的變量均為國際單位制.

用平均應(yīng)力表述,應(yīng)力平衡方程也可寫為[47-48]

式中,v為泊松比,σm為平均應(yīng)力,α 為Biot 系數(shù).

1.2 質(zhì)量守恒方程

采用組分模型對CO2地質(zhì)封存進行數(shù)值模擬,α組分的質(zhì)量守恒方程為

對于 α 組分,其質(zhì)量累積項和傳輸項的表達式分別為

式中,下標 β 表示相(包含油相、氣(CO2) 相和水相),? 為孔隙度,ρ 為密度,S為飽和度,表示 α 組分在 β 相中的質(zhì)量分數(shù),Fadv為對流項,Fdif為擴散項,V為質(zhì)量流速,j為擴散速度.

巖石基質(zhì)中流體的流動符合達西定律[49],其表達式為

式中,k為滲透率,kr為相對滲透率,μ 為黏度,g為重力加速度,D為深度.

考慮CO2擴散,擴散項可采用 Fick 定律描述

式中,D為擴散系數(shù).

對于低速非達西滲流和高速非達西滲流,達西定律描述的滲流速度與水力梯度之間的線性關(guān)系不再成立[50].對于低滲透-致密儲層中的低速非達西滲流,只有驅(qū)動壓力梯度 ?P-ρg?D大于擬啟動壓力梯度L時,流體才能流動,其表達式為[51]

對于裂縫中的高速非達西滲流,可采用Forchheimer定律描述,其表達式為

式中,χ為高速非達西滲流系數(shù).

1.3 能量守恒方程

能量守恒方程與質(zhì)量守恒方程具有相同的形式

式中,ρR為巖石密度,CR為巖石比熱,u為比內(nèi)能.

忽略輻射傳熱,熱量的流量由熱傳導項與熱對流項組成,其表達式為

式中,κ 為熱傳導系數(shù),h為比焓.

1.4 化學反應(yīng)

CO2注入地層后發(fā)生的化學反應(yīng)主要包括CO2與地層流體之間的離子交換和溶解的CO2與地層巖石之間的礦物反應(yīng)[21].對于T-H-M-C 多場耦合問題,化學反應(yīng)過程通過質(zhì)量和能量守恒方程中的化學反應(yīng)源項R來表述,其表達式為[12,26]

式中,ξ 為化學反應(yīng)計量系數(shù),r為化學反應(yīng)速率.礦物溶解和礦化反應(yīng)的化學反應(yīng)速率計算公式為[12]

式中,A為反應(yīng)接觸面積,kR為速率常數(shù),Q為化學親和性,Keq為平衡常數(shù).化學反應(yīng)速率r若為正,表示礦化沉淀;若為負,表示礦物溶解.

根據(jù)不同的反應(yīng)機理,礦物溶解和沉淀不僅受到純H2O 的催化,而且受到H+和OH-的影響.因此,速率常數(shù)kR的計算公式包含3 項[47]

上述平衡常數(shù)和速率常數(shù)對化學反應(yīng)動力學模型的準確性至關(guān)重要,然而由于地球化學過程的復(fù)雜性,特別是相關(guān)機制的持續(xù)時間極長,準確的測定較為困難[26].

1.5 巖石應(yīng)力-應(yīng)變本構(gòu)關(guān)系

采用雙孔隙彈性理論描述巖石的應(yīng)力-應(yīng)變關(guān)系[52]

式中,G為剪切模量,為裂縫(i=1)和基質(zhì)(i=2)中的流體和熱效應(yīng)耦合系數(shù),δij為Kronecker符號.

基于等效連續(xù)介質(zhì)模型,裂縫性巖石的本構(gòu)關(guān)系為

式中,ε 為正應(yīng)變,γ 為剪應(yīng)變,σ′為有效應(yīng)力,τ 為剪應(yīng)力,GM和GF分別為基質(zhì)和裂縫的剪切模量,EM和EF分別為基質(zhì)和裂縫的楊氏模量,ECF為裂縫的壓縮模量,αM和 αF分別為巖體中基質(zhì)和裂縫的體積分數(shù).

由于裂縫性巖石中基質(zhì)占據(jù)的體積遠大于裂縫占據(jù)的體積,因此可認為 αM≈1,式(17)中的裂縫性巖石柔度矩陣C可以寫為基質(zhì)與裂縫柔度矩陣的加權(quán)平均

式中,CM為基質(zhì)柔度矩陣,CF為裂縫柔度矩陣.

裂縫的力學性質(zhì)與巖石基質(zhì)存在很大差別,其力學參數(shù)和滲流能力會隨施加在裂縫壁面上的法向載荷產(chǎn)生動態(tài)變化,需要通過裂縫柔度矩陣獲得裂縫壁面位移與巖石應(yīng)力間的關(guān)系,以對裂縫孔隙度、滲透率等參數(shù)進行更新[53].對于天然裂縫在法向加載下的變形過程,Goodman 等[54]指出裂縫在法向受壓下的變形本構(gòu)關(guān)系近似為雙曲線型模型.Bandis等[55]對大量存在天然裂縫的巖塊進行實驗,總結(jié)出裂縫受壓變形的經(jīng)驗性雙曲線型本構(gòu)公式.

1.6 相平衡計算及輔助方程

地層條件下油-水-氣(CO2)系統(tǒng)的相平衡計算是油氣藏地質(zhì)封存CO2數(shù)值模擬研究的核心問題之一.油-水-氣(CO2)的相平衡計算主要有平衡常數(shù)法、基于Rachford-Rice 函數(shù)的閃蒸計算法和Gibbs自由能最小化法3 種計算方法[56].

平衡常數(shù)法采用基于實驗的經(jīng)驗公式來計算平衡常數(shù),計算速度快、穩(wěn)定性良好,雖然避免了復(fù)雜、耗時的相平衡計算,但由于缺乏嚴密的熱力學理論基礎(chǔ),計算結(jié)果往往偏差較大.基于Rachford-Rice 函數(shù)的閃蒸計算法主要用于油-氣(CO2)兩相體系的相平衡計算,同平衡常數(shù)法相比較,有較高的精度,但計算量較平衡常數(shù)法大.適用于油-水-氣(CO2) 三相體系的閃蒸計算法還有待研究和發(fā)展.Gibbs 自由能最小化方法理論嚴密、適應(yīng)性廣、計算穩(wěn)定和收斂性好,比較適合于油-水-氣(CO2)三相的相平衡計算,但是因為需要反復(fù)迭代,計算量比平衡常數(shù)法和閃蒸計算法都要大很多.

由于滲流場、溫度場、應(yīng)力場和化學場之間存在耦合作用,在求解T-H-M-C 多場耦合問題的過程中,需要實時更新計算的物性參數(shù),主要包括毛細管力、滲透率、相對滲透率、孔隙度、流體黏度和密度等.相對滲透率曲線可以采用Brooks-Corey 模型.考慮應(yīng)力場、溫度場和化學場對孔隙度的影響[12,53]

式中,?0為0 應(yīng)力下的初始孔隙度大小,σm為平均總應(yīng)力(壓為正),Kb為巖體的體積模量,βT為熱膨脹系數(shù),T為當前時刻溫度,cs為固相s組分濃度,下標 0 表示初始狀態(tài).

基質(zhì)滲透率和毛細管力可根據(jù)孔隙度變化進行更新

式中,k0和Pc0分別為0 應(yīng)力下的初始滲透率和初始毛細管力,d為實驗測定的參數(shù).

溫度變化會造成流體黏度的變化,一般來說隨著溫度的升高流體黏度降低,可采用經(jīng)驗公式描述[53]

式中,A,B,C和D均為實驗測定參數(shù).

2 耦合問題數(shù)值解法

對于CO2地質(zhì)封存的T-H-M-C 多場耦合問題,可采用有限元和有限體積混合的方法進行離散,即對于式(1)的應(yīng)力平衡方程采用有限元法離散,對于式(3)質(zhì)量守恒方程和式(10)能量守恒方程采用有限體積法進行離散.離散后得到的方程可采用全耦合、迭代耦合、弱耦合、顯式耦合和擬耦合等數(shù)值方法進行求解,如圖3 所示.

圖3 常用的耦合問題數(shù)值解法Fig.3 Commonly used coupling solutions

全耦合方法即為全隱式的數(shù)值求解方法,該方法對離散后的應(yīng)力守恒方程、質(zhì)量守恒方程、能量守恒方程和化學反應(yīng)方程同時求解,所有參數(shù)變量需要在每個迭代步內(nèi)進行更新[13].全耦合方法的收斂性最好,計算精度最高,但由于耦合問題的方程數(shù)量龐大,全耦合方法的計算效率往往較低.

為了提高多場耦合問題數(shù)值模擬的計算速度,通常采用半隱式的方法進行求解,即優(yōu)先對一部分控制方程(一般是質(zhì)量守恒方程和能量守恒方程)進行隱式求解,之后將求解得到的未知量作為已知量代入另一部分尚未求解的控制方程中,求解其余的未知量.顯式耦合方法、弱耦合方法和迭代耦合方法均屬于半隱式求解方法.

迭代耦合方法中,在每個時間步的牛頓迭代內(nèi)首先對質(zhì)量和能量守恒方程進行隱式求解,將計算得到的平均孔隙壓力傳遞至應(yīng)力平衡方程中以求解巖石的變形情況,根據(jù)巖石變形程度對孔隙度、滲透率等關(guān)鍵物性參數(shù)進行更新,開始下一迭代循環(huán),分別對質(zhì)量守恒方程、能量守恒方程和應(yīng)力平衡方程進行計算,直至達到收斂條件進入下一時間步.與全耦合方法相比,迭代耦合方法顯著提高計算效率[57].

在顯式耦合中,滲流場、溫度場和應(yīng)力場的計算獨立進行,每個時間步中信息僅從滲流場、溫度場向應(yīng)力場單向傳遞一次,即在每一時間步中僅進行一次應(yīng)力平衡方程計算和參數(shù)更新,而應(yīng)力場的信息則不會反饋給滲流場、溫度場.顯示耦合方法是耦合程度較低的方法,計算精度有所降低,但計算效率高、收斂性好,適用于大尺度、復(fù)雜地質(zhì)條件和復(fù)雜巖體本構(gòu)模型的力學模擬[13].

弱耦合方法是在顯式耦合方法的基礎(chǔ)上進一步減少應(yīng)力平衡方程的計算頻率,在每個時間步內(nèi)對滲溫場和應(yīng)力場進行獨立計算,每隔數(shù)個時間步才進行一次應(yīng)力平衡方程的計算和參數(shù)更新,如此反復(fù)循環(huán)迭代直至結(jié)束.該方法計算結(jié)果精度有所降低,但計算速度快、收斂性好、適用范圍廣,能夠耦合不同用途的程序[13].

擬耦合方法并不進行應(yīng)力場求解,僅通過經(jīng)驗公式或數(shù)據(jù)表建立孔隙度、滲透率等與孔隙壓力的關(guān)系.該方法無法對巖體的變形過程進行模擬,因此精度有限.

3 CO2 地質(zhì)封存風險分析

分析CO2長久封存過程的安全性是多場耦合數(shù)值模擬的關(guān)鍵應(yīng)用.圖4 展示的是在CO2地質(zhì)封存中儲層-蓋層可能發(fā)生的典型地質(zhì)風險,主要包括CO2泄漏、地表變形、蓋層破壞和地震誘發(fā)等[39,58-61].蓋層破壞的常見形式包括斷層活化和蓋層壓裂(裂縫擴展),斷層活化主要與斷層剪切破壞有關(guān),蓋層壓裂(裂縫擴展)主要與蓋層拉伸破壞有關(guān)[58].

圖4 CO2 封存儲層-蓋層中的地質(zhì)風險示意圖[39]Fig.4 Geomechanical risks during CO2 geological storage[39]

3.1 CO2 泄漏

CO2通過蓋層泄漏的主要機理包括壓力驅(qū)動的流體對流[62]和濃度驅(qū)動的分子擴散[63],如圖5 所示.斷層/裂縫和巖石基質(zhì)存在巨大的滲透率差異[64],斷層核心可視為阻礙流動的屏障,而斷層兩側(cè)破碎的損傷區(qū)可視為沿斷層的流動通道[35,65].在CO2的注入過程中,流動路徑主要受裂縫網(wǎng)絡(luò)幾何形狀控制,CO2在遇到裂縫時優(yōu)先沿裂縫方向運移,在裂縫周圍擴散并集中分布在裂縫周圍,如圖6 所示.隨著CO2的積累,碳酸鹽巖礦物在裂縫交叉點沉淀,降低了裂縫網(wǎng)絡(luò)的整體連通性,低滲透帶內(nèi)的物理圈閉效應(yīng)會顯著增強[66].裂縫組合導致儲層滲透率各向異性,相對滲透率和毛管壓力效應(yīng)也會影響儲層各向異性,從而顯著影響CO2的分布[67-69].

圖5 通過斷層泄漏到地表的潛在途徑[70]Fig.5 Potential CO2 leakage path through a fault to ground surface[70]

圖6 裂縫對CO2 運移的影響(CO2 飽和度分布)[67]Fig.6 Effect of fractures on CO2 migration (CO2 saturation distribution)[67]

CO2在多孔介質(zhì)中的擴散速度相對于對流速度較慢.CO2的擴散速度對流體-巖石反應(yīng)較為敏感,流體-巖石反應(yīng)可能提高(礦物溶解孔隙度增加)或降低(礦化反應(yīng)消耗CO2)擴散速度[70].但是現(xiàn)有研究表明,擴散并非深層儲層CO2的主要泄漏機制.根據(jù)Wang 等[71]的數(shù)值模擬結(jié)果,將CO2以27 MPa的注入壓力注入儲層(儲層壓力為19 MPa),如圖7所示,穩(wěn)定注入1000 年后CO2擴散距離僅為0.93 m.此外,根據(jù)Busch 等[70]的估算,純水通過擴散透過10 m 厚度的薄蓋層(有效擴散系數(shù)為2.0×10-9m2/s),大約需要105年.

圖7 CO2 擴散運移數(shù)值模擬結(jié)果[71]Fig.7 Numerical simulation of CO2 diffusion[71]

3.2 地表變形

阿爾及利亞In Salah 二氧化碳地質(zhì)封存項目已經(jīng)觀測到由CO2注入導致的地表抬升現(xiàn)象[72-74],如圖8 所示.地表變形的分析方法主要包括解析法、半解析法和數(shù)值模擬方法[33].Fjaer 等[75]針對簡化的一維橫向薄儲層提出計算儲層垂直膨脹的解析式,Xu 等[76]針對半無限空間的問題推導了壓力源地面變形的解析解,Selvadurai[77]提出地表巖層和深部巖體之間相互作用的基本模型,并給出地表巖層隆起的半解析式.對于實際CO2地質(zhì)封存問題中復(fù)雜巖體的變形,需要采用數(shù)值模擬方法進行計算.

圖8 In Salah KB-502 井地面垂直位移的演化[72]Fig.8 Transient evolution of vertical ground displacement at the In Salah KB-502 well[72]

Rinaldi 等[72]通過數(shù)值模擬分析了KB-502 地表位移演化情況,并與監(jiān)測數(shù)據(jù)進行了比對,如圖9 所示.郝術(shù)仁等[78]模擬研究了不同流量和不同壓力注入CO2引起地表位移的變化規(guī)律.Siriwardane 等[45]模擬分析了蓋層斷層/裂縫帶對地表位移的影響,無斷層/裂縫時地表變形較小且對稱,有斷層/裂縫時地表位移較大且非對稱.裂縫帶滲透率較低時,地面垂直位移較小;蓋層裂縫帶離注入源越近,地表變形幅度越大.

圖9 In Salah KB-502 井2006 年12 月23 日地表變形模擬結(jié)果[72]Fig.9 The simulation results of ground displacement at the In Salah KB-502 well on December 23,2006[72]

3.3 剪切破壞及斷層活化

流體壓力的演化直接影響蓋層的地質(zhì)力學穩(wěn)定性,蓋層的低滲透性阻止了超壓向蓋層快速傳播,但在靠近儲層界面的蓋層部分,由于流體壓力的積聚導致有效應(yīng)力降低,應(yīng)力狀態(tài)向失效狀態(tài)演化[79-80],如圖10 所示.根據(jù)Mohr-Coulomb 強度失效準則(如下式),當作用在斷層面上的剪應(yīng)力超過斷層的剪切強度時,斷層活化[81-82]

圖10 Mohr-Coulomb 強度失效準則[59]Fig.10 Mohr-Coulomb failure criterion[59]

表1 國外部分CO2 封存儲層的力學參數(shù)[59]Table 1 Stress state at several CO2 injection sites abroad[59]

CO2注入地層后流體壓力積聚、溫度場的變化,會改變有效應(yīng)力場,化學反應(yīng)會降低巖石的強度[84].孔隙壓力的增加會導致莫爾圓半徑減小并左移(參見圖11 中的紅色莫爾圓).CO2到達注入井底部時溫度低于儲層溫度,注入井周圍的巖石發(fā)生冷卻,溫度變化誘發(fā)溫度應(yīng)力[59],應(yīng)力狀態(tài)向剪切破壞狀態(tài)轉(zhuǎn)移(參見圖11 中的藍色莫爾圓).溫度應(yīng)力主要取決于溫差、巖石基質(zhì)的體積模量和體積熱膨脹系數(shù)[85],由溫度變化 ΔT引起的孔隙壓力變化 ΔP可以根據(jù)下式進行估算[86]

圖11 應(yīng)力狀態(tài)耦合效應(yīng)示意圖[59]Fig.11 Coupled effects on stress state[59]

式中,Δ ζ 為流體應(yīng)變,βf為流體體積熱膨脹系數(shù),為巖石基質(zhì)體積熱膨脹系數(shù),Kf為流體體積模量,Kb為巖石體積模量.

斷層除非位于注入井附近,一般不會直接受到溫度差的影響.但在注入后期,儲層巖石冷卻收縮,水平應(yīng)力減小,從而引起遠場應(yīng)力的變化,導致遠離注入井的斷層可能處于不穩(wěn)定狀態(tài)(參見圖11 中的綠色莫爾圓)[87].

3.4 拉伸破壞及裂縫擴展

蓋層由于CO2注入發(fā)生拉伸破壞時,裂縫起裂和擴展方向依賴于最小主應(yīng)力方向和儲層非均質(zhì)性[88].當孔隙壓力P超過蓋層的最小主應(yīng)力 σ3時,發(fā)生拉伸破壞,裂縫開始萌生并擴展[89-90]

當儲層或蓋層巖石剪應(yīng)力超過剪切強度時,儲層或蓋層已有的裂縫會發(fā)生剪切滑移[91],如圖10 所示.如果儲層內(nèi)部形成壓裂裂縫、裂縫不向蓋層擴展,或者裂縫僅在蓋層有限的范圍內(nèi)延伸時,則不會發(fā)生泄漏,且能夠提高CO2的注入能力[90-93].

4 數(shù)值模擬方法研究需解決的主要問題

在地下能源和資源領(lǐng)域,國際上已經(jīng)研發(fā)出了諸多用于多孔介質(zhì)多相流數(shù)值模擬的程序和軟件,能夠針對CO2地質(zhì)封存進行T-H-M-C 多場耦合或部分多物理場耦合的模擬計算,表2 給出了目前國外主要的可用于CO2地質(zhì)封存的數(shù)值模擬器[25-26,32].針對CO2地質(zhì)封存,近期國際石油工程師協(xié)會發(fā)布了第11 個標準算例(the 11th SPE CSP)[94].該算例包括一個實驗室尺度的二維模型、一個油藏尺度的二維和一個油藏尺度的三維地質(zhì)模型,可用于對不同的數(shù)值模擬方法和模擬器進行對比研究.由于CO2地質(zhì)封存需要對目標封存場地儲層-蓋層體系的三維地質(zhì)結(jié)構(gòu)進行精細化建模,需要考慮大時間尺度和大空間尺度下的多場耦合作用,所以在計算 效率和計算能力方面,對數(shù)值模擬技術(shù)提出了更高的要求.除CO2地質(zhì)封存過程中復(fù)雜的相態(tài)計算和地球化學反應(yīng)涉及的本構(gòu)關(guān)系等問題以外,在數(shù)值模擬方法研究方面,還需要解決如下問題.

表2 國外部分可用于CO2 地質(zhì)封存的數(shù)值模擬器Table 2 Overview of numerical simulators for CO2 geological storage

4.1 基于角點網(wǎng)格的有限單元法

角點網(wǎng)格可以準確描述地層界面、斷層面和尖滅等復(fù)雜地質(zhì)情況以及儲層的空間展布和屬性分布特征,因此通常都采用角點網(wǎng)格進行三維地質(zhì)體精細化建模[95-96].角點網(wǎng)格系統(tǒng)中,通常存在重復(fù)節(jié)點,其節(jié)點無法傳遞力、位移和溫度等連續(xù)性變量,相鄰網(wǎng)格之間允許錯動、斷開和夾層,因此基于角點網(wǎng)格建立的地質(zhì)模型無法直接用于三維地應(yīng)力的有限元模擬.基于角點網(wǎng)格數(shù)據(jù)進行有限元網(wǎng)格的剖分,幾何處理上較為困難,且有限元網(wǎng)格數(shù)量十分巨大.因此,需要結(jié)合角點網(wǎng)格和有限元方法的優(yōu)點,建立基于角點網(wǎng)格的地應(yīng)力有限元模擬方法,這樣既可以完整準確地保留地質(zhì)模型復(fù)雜的幾何特征,也可以利用有限單元準確地計算三維應(yīng)力場[97].

4.2 多場耦合的復(fù)合介質(zhì)模型

CO2注入涉及的地層區(qū)域往往包含多尺度的裂縫系統(tǒng),從長度米級的小裂縫到公里級的大斷層[33].對于多尺度裂縫系統(tǒng),僅采用多重介質(zhì)模型,會喪失斷層或大裂縫準確的幾何信息;若僅采用離散裂縫模型,計算量極大,不適用于實際封存問題的模擬;采用嵌入式離散裂縫模型,連續(xù)介質(zhì)雖然能夠采用結(jié)構(gòu)化網(wǎng)格,但是需要細分裂縫,并對所有不同尺度的裂縫進行獨立表征,將導致網(wǎng)格數(shù)量急劇增加,造成極大的計算負擔,且不適用于各向異性的連續(xù)介質(zhì).

為了更準確表征裂縫型儲層和提高數(shù)值模擬計算效率,應(yīng)分別對不同尺度的裂縫進行建模: 小型裂縫縫長遠小于網(wǎng)格尺度,對于滲流特征的影響局限在基質(zhì)網(wǎng)格內(nèi)部,可采用等效單重介質(zhì)模型建模;中等裂縫縫長與網(wǎng)格尺度相當,集中分布常以裂縫簇的形式賦存,其方向性與連通性對集中分布區(qū)域的滲流特征造成較大影響,可采用雙重介質(zhì)模型建模;大型裂縫和斷層數(shù)量有限,但對于整體滲流場的分布具有十分重要的影響,為了準確描述每一條大型裂縫或斷層的水力特征,應(yīng)采用整體嵌入式離散裂縫模型建模[98].由此形成復(fù)合介質(zhì)模型,將離散裂縫模型、多重介質(zhì)模型和整體嵌入式離散裂縫模型結(jié)合在一起,進而在其基礎(chǔ)上建立相應(yīng)的T-H-M-C 多場耦合復(fù)合介質(zhì)數(shù)學模型.

4.3 分區(qū)域的有限元-有限體積混合數(shù)值離散方法

CO2地質(zhì)封存數(shù)值模擬需要考慮儲層的上覆蓋層、周邊地層和下部地層,建立基層-儲層-蓋層-地表全地層模型[87,99].CO2地質(zhì)封存場地尺度的地質(zhì)建模,往往達到百萬級以上的網(wǎng)格數(shù)量[100-101],模擬計算的規(guī)模十分龐大.圖12 是勝利油田G89-1 區(qū)塊的CO2地質(zhì)封存全地層模型,其網(wǎng)格數(shù)量約240 萬,待解的自由度數(shù)已達約1500 萬.

圖12 勝利油田G89-1 區(qū)塊全地層模型Fig.12 The full synthetic field model of G89-1 block in Shengli oilfield

CO2地質(zhì)封存的應(yīng)力計算區(qū)域通常遠大于滲流和溫度計算區(qū)域.CO2地質(zhì)封存數(shù)值模擬中對于滲流和溫度計算具有較高的要求,儲層內(nèi)網(wǎng)格要保證一定的細密程度.而全地層模型應(yīng)力場的計算則并不需要采用相同細密程度的網(wǎng)格劃分,由于應(yīng)力平衡方程采用有限元法計算,其主變量個數(shù)也遠大于滲流場和溫度場的計算,采用粗網(wǎng)格可以避免運算資源的浪費[54].因此,對于全地層模型的多場耦合問題,在滲流和傳熱計算區(qū)域內(nèi),采用細密網(wǎng)格進行剖分(如圖13 所示),基于有限體積法對質(zhì)量守恒方程和能量守恒方程進行數(shù)值離散;在應(yīng)力計算區(qū)域內(nèi),采用粗網(wǎng)格進行剖分,基于有限元法對應(yīng)力平衡方程進行數(shù)值離散.根據(jù)粗細網(wǎng)格間的映射關(guān)系,將離散方程聯(lián)立求解,從而形成分區(qū)域的有限元-有限體積混合數(shù)值離散方法.這一方法能夠在保證儲層滲流場和溫度場計算精度的同時,降低應(yīng)力場求解帶來的計算負擔,有效提高計算效率.

圖13 分區(qū)域網(wǎng)格劃分示意圖[54]Fig.13 Sub-regional grid division diagram[54]

4.4 高效的并行求解技術(shù)

CO2地質(zhì)封存問題涉及的時間尺度和空間尺度較大.在多場耦合作用的數(shù)值模擬方面,需要發(fā)展高效的并行求解技術(shù).采用并行計算技術(shù),可將一個復(fù)雜的大型計算問題分成若干相互獨立可以同時計算的子任務(wù)[102-103].常見的分解方法大體分為時間和空間上的并行,時間并行方法常采用流水線的方式實現(xiàn),空間上并行是目前絕大多數(shù)計算問題的并行解決方案,目前常采用區(qū)域分解算法[102].

對圖12 的全地層模型注CO2問題,采用并行技術(shù)進行多場耦合的數(shù)值模擬,計算得到注入3.2 年時的豎向地應(yīng)力如圖14 所示.這一計算雖然采用了并行技術(shù)(Intel Xeon Platinum CPU 48 核),但模擬耗時仍達40 h.

圖14 勝利油田G89-1 區(qū)塊儲層豎向有效應(yīng)力模擬結(jié)果Fig.14 The simulation results of vertical effective stress at the Shengli oilfield G89-1 block

為了大幅度提高CO2地質(zhì)封存大規(guī)模模型數(shù)據(jù)處理和數(shù)值計算的效率,可基于MPI,OpenMP,CUDA 等工具研發(fā)高效的并行求解技術(shù).其中,MPI主要適用于分布式存儲的可拓展的并行計算機和工作站機群;OpenMP 主要適用于共享內(nèi)存多處理器系統(tǒng)和多核處理器等體系結(jié)構(gòu),且與MPI 類似,適用于粗粒度任務(wù)并行計算;而CUDA 適用于包含GPU 設(shè)備的單個計算節(jié)點,或者與MPI/OpenMP 協(xié)同應(yīng)用于包含多個GPU 設(shè)備的集群,主要適用于細粒度數(shù)據(jù)并行計算[104].

5 結(jié)論

儲層-蓋層系統(tǒng)的力學分析和穩(wěn)定性評估是實現(xiàn)CO2地質(zhì)封存工業(yè)化的關(guān)鍵問題.CO2注入地層后會導致較大的流體壓力積聚、有效應(yīng)力場變化,引發(fā)CO2泄漏、地表位移、蓋層破壞和地震誘發(fā)等地質(zhì)風險.蓋層破壞的常見形式包括由剪切破壞引發(fā)的斷層活化和由拉伸破壞引發(fā)的蓋層壓裂(裂縫擴展).蓋層裂縫的起裂、擴展和斷層的活化、滑移,形成了CO2逸出的潛在通道,對蓋層完整性構(gòu)成了重大威脅.CO2泄漏的主要機理是濃度驅(qū)動的分子擴散和壓力驅(qū)動的對流,流動路徑主要受裂縫網(wǎng)絡(luò)幾何形狀控制,但也受到地層特性和地球化學反應(yīng)等影響,涉及溫度場-滲流場-力學場-化學場的多物理場耦合.

CO2地質(zhì)封存數(shù)值模擬的研究對象是巖石固體骨架與孔隙流體(CO2、水和油等)組成的多相多組分系統(tǒng),其數(shù)學模型主要包括基本控制方程和相平衡計算方法等.基本控制方程包括應(yīng)力平衡方程、質(zhì)量守恒方程、能量守恒方程和化學反應(yīng)等.對離散后的控制方程耦合求解,求解方法大體可以分為全耦合、迭代耦合、弱耦合、顯式耦合和擬耦合等算法.

為了適應(yīng)CO2地質(zhì)封存儲層-蓋層體系精細化建模、大時間尺度和大空間尺度的多物理場耦合對計算模型和計算效率的要求,還需要在裂縫性介質(zhì)數(shù)值模擬方法方面進行重點研究.目前,三維地質(zhì)體精細化建模通常都采用角點網(wǎng)格.因此,需要將角點網(wǎng)格和有限元方法的優(yōu)點相結(jié)合,建立基于角點網(wǎng)格的地應(yīng)力有限元模擬方法.將離散裂縫模型、多重介質(zhì)模型和整體嵌入式離散裂縫模型相結(jié)合,建立相應(yīng)的T-H-M-C 多場耦合復(fù)合介質(zhì)數(shù)學模型,能夠更準確地表征裂縫型儲層并提高數(shù)值模擬的計算效率.由于CO2地質(zhì)封存全地層模型的應(yīng)力計算區(qū)域遠大于滲流和溫度計算區(qū)域,采用分區(qū)域的有限元-有限體積混合數(shù)值離散方法,能夠在保證滲流場和溫度場計算精度的同時,降低應(yīng)力場求解的計算量.由于CO2地質(zhì)封存問題涉及的時間尺度和空間尺度較大,在多物理場耦合方面需要發(fā)展高效的并行求解技術(shù).