CO2 地質封存風險分析的多場耦合數值模擬技術綜述1)

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

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

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

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

引言

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

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

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

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

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

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

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

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

1 數學模型

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

1.1 應力平衡方程

根據可變形多孔介質多相流孔隙熱彈性理論,應力平衡方程為[12,46]

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

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

式中,v為泊松比,σm為平均應力,α 為Biot 系數.

1.2 質量守恒方程

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

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

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

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

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

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

式中,D為擴散系數.

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

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

式中,χ為高速非達西滲流系數.

1.3 能量守恒方程

能量守恒方程與質量守恒方程具有相同的形式

式中,ρR為巖石密度,CR為巖石比熱,u為比內能.

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

式中,κ 為熱傳導系數,h為比焓.

1.4 化學反應

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

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

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

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

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

1.5 巖石應力-應變本構關系

采用雙孔隙彈性理論描述巖石的應力-應變關系[52]

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

基于等效連續介質模型,裂縫性巖石的本構關系為

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

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

式中,CM為基質柔度矩陣,CF為裂縫柔度矩陣.

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

1.6 相平衡計算及輔助方程

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

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

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

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

基質滲透率和毛細管力可根據孔隙度變化進行更新

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

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

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

2 耦合問題數值解法

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

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

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

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

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

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

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

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

3 CO2 地質封存風險分析

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

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

3.1 CO2 泄漏

CO2通過蓋層泄漏的主要機理包括壓力驅動的流體對流[62]和濃度驅動的分子擴散[63],如圖5 所示.斷層/裂縫和巖石基質存在巨大的滲透率差異[64],斷層核心可視為阻礙流動的屏障,而斷層兩側破碎的損傷區可視為沿斷層的流動通道[35,65].在CO2的注入過程中,流動路徑主要受裂縫網絡幾何形狀控制,CO2在遇到裂縫時優先沿裂縫方向運移,在裂縫周圍擴散并集中分布在裂縫周圍,如圖6 所示.隨著CO2的積累,碳酸鹽巖礦物在裂縫交叉點沉淀,降低了裂縫網絡的整體連通性,低滲透帶內的物理圈閉效應會顯著增強[66].裂縫組合導致儲層滲透率各向異性,相對滲透率和毛管壓力效應也會影響儲層各向異性,從而顯著影響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在多孔介質中的擴散速度相對于對流速度較慢.CO2的擴散速度對流體-巖石反應較為敏感,流體-巖石反應可能提高(礦物溶解孔隙度增加)或降低(礦化反應消耗CO2)擴散速度[70].但是現有研究表明,擴散并非深層儲層CO2的主要泄漏機制.根據Wang 等[71]的數值模擬結果,將CO2以27 MPa的注入壓力注入儲層(儲層壓力為19 MPa),如圖7所示,穩定注入1000 年后CO2擴散距離僅為0.93 m.此外,根據Busch 等[70]的估算,純水通過擴散透過10 m 厚度的薄蓋層(有效擴散系數為2.0×10-9m2/s),大約需要105年.

圖7 CO2 擴散運移數值模擬結果[71]Fig.7 Numerical simulation of CO2 diffusion[71]

3.2 地表變形

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

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

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

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

3.3 剪切破壞及斷層活化

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

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

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

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

圖11 應力狀態耦合效應示意圖[59]Fig.11 Coupled effects on stress state[59]

式中,Δ ζ 為流體應變,βf為流體體積熱膨脹系數,為巖石基質體積熱膨脹系數,Kf為流體體積模量,Kb為巖石體積模量.

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

3.4 拉伸破壞及裂縫擴展

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

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

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

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

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

4.1 基于角點網格的有限單元法

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

4.2 多場耦合的復合介質模型

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

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

4.3 分區域的有限元-有限體積混合數值離散方法

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

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

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

圖13 分區域網格劃分示意圖[54]Fig.13 Sub-regional grid division diagram[54]

4.4 高效的并行求解技術

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

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

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

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

5 結論

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

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

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