多孔介質(zhì)中單相流地應力耦合問題的二重網(wǎng)格混合有限元法

摘要： 本文提出了多孔介質(zhì)中單相流耦合地應力問題的二重網(wǎng)格混合有限元法. 本文首先給出了多孔介質(zhì)中的單相流的控制方程及地應力問題的控制方程，推導了多孔介質(zhì)中的微可壓縮流體的單相流耦合地應力問題的非線性控制方程. 然后，為數(shù)值求解該方程本文提出了二重網(wǎng)格混合有限元法，將問題轉(zhuǎn)化為在粗網(wǎng)格上求解小規(guī)模非線性問題、在細網(wǎng)格上求解大規(guī)模線性問題，以提高計算效率. 最后，數(shù)值算例驗證了方法的有效性.

關(guān)鍵詞： 多孔介質(zhì)； 單相流耦合地應力問題； 二重網(wǎng)格法； 混合有限元

中圖分類號： O241. 82 文獻標志碼： A DOI： 10. 19907/j. 0490-6756. 2024. 041002

1 引言

油氣藏數(shù)值模擬是一個復雜交叉學科問題，需要地質(zhì)工程、滲流學、計算數(shù)學、熱力學和計算機科學等多學科共同研究. 油氣藏數(shù)值模擬需要通過計算機軟件來模擬地下油氣在不同條件下的變化，以便為制定開采方案提供參考，達到油氣高產(chǎn)的目標或預測油氣產(chǎn)量.

目前，多數(shù)油氣藏數(shù)值模擬都只關(guān)注多孔介質(zhì)中的油氣流動. 但是，在實際的頁巖油氣儲層中，流場和地應力場之間存在重要的相互作用，因而在構(gòu)建模型及數(shù)值模擬時考慮地應力的影響很重要. 對于這個問題，不同領(lǐng)域?qū)W者均有研究，這些研究共同推動人們對地應力變化規(guī)律及流場與地應力場的耦合相互作用的認識［1-5］. Gudala 等［1］使用單相流與地質(zhì)力學耦合模型研究了縫隙多孔介質(zhì)中單相流的流動，并對模型進行了數(shù)值求解.Li 等［4］把巖體看作多孔介質(zhì)，計入溫度和孔隙率的影響，建立了高溫高壓及高孔隙壓環(huán)境下的巖體多孔介質(zhì)流變模型. Chen 等［5］基于計算流體力學和巖土力學，考慮流場和應力場的相互作用建立了流-固耦合模型，等.

在數(shù)學模型基礎(chǔ)上，數(shù)值模擬還需要合適的計算方法與算法設(shè)計. 其中，用混合有限元法數(shù)值求解耦合地應力問題已有一些研究［6-9］. Zhang 等［7］使用多尺度混合有限元法研究了可壓縮流體在多孔介質(zhì)中的滲流. Li 等［8］建立了單相滲流問題的有限元方程，并針對復雜邊界情況用有限元法進行了模擬計算，給出了混合邊界油藏壓力分布的規(guī)律. 進一步，鑒于單相流耦合地應力模型的非線性可能導致計算效率低，并考慮到地應力模擬時間較長、數(shù)值模擬較慢，研究者提出了二重網(wǎng)格混合有限元法［10-13］. Zhang 等［10］提出了一種求解半線性反應擴散方程的二重網(wǎng)格混合有限元法. Rui 等［11］提出了多孔介質(zhì)Darcy-Forchheimer 流動的二重網(wǎng)格塊中心差分法，等.

本文研究多孔介質(zhì)中單相流的耦合地應力問題，提出了一個二重網(wǎng)格混合有限元格式，構(gòu)造了格式對應的算法，并用數(shù)值算例驗證了方法的有效性.

2 模型

在推導問題對應的控制模型前，我們假設(shè)：（i）在耦合流動過程中溫度恒定；（ii）考慮基巖和流體的壓縮性時，基巖與裂縫中的流體流動均滿足達西定律；（iii）巖石變形為微小線彈性變形；（iv）不考慮慣性力，并忽略毛細管力與重力作用.

準確模擬多孔介質(zhì)中單相流耦合地應力問題需要研究流體流動與地應力的耦合問題，以及巖石結(jié)構(gòu)因地應力而產(chǎn)生的線彈性變形問題. 式（1）為多孔介質(zhì)中單相流的控制方程.

得到線性問題（即方程（13））的二重網(wǎng)格混合有限元法的解（ vn + 1h ，pn + 1h ，un + 1h ）∈ Vh × Wh × Uh.

5 數(shù)值算例

考慮10 m × 10 m 的二維區(qū)域，區(qū)域的右邊界和上邊界牽引力為5 N. 對該區(qū)域采用64 × 64 的一致三角形單元離散，如圖1 所示，對應的巖石流體參數(shù)如表1 所示. 假設(shè)該區(qū)域有5 口生產(chǎn)井存在，坐標分別為［5，5］，［2. 5，2. 5］，［2. 5，7. 5］，［7. 5，2. 5］及［7. 5，7. 5］，生產(chǎn)速率為 2 × 10-7 s-1，模擬時間為100 d，時間步長為1440 min，即1 d.

下面我們先使用牛頓迭代法模擬運行100 d，模擬運行100 d 后區(qū)域中的壓力及孔隙率如圖2 所示，速度及剪切應力如圖3 所示，最大主應力及最小主應力如圖4 所示，x 方向位移及y 方向位移如圖5 所示.

下面我們使用二重網(wǎng)格混合有限元方法再次進行計算. 首先，在粗網(wǎng)格上進行8 × 8 的一致三角形單元離散，再劃分同圖1 的64 × 64 細網(wǎng)格.我們先在粗網(wǎng)格上求解非線性問題，再在細網(wǎng)格上求解線性問題，得到數(shù)值解. 同樣，我們模擬100 d，時間步長取為1 d.

模擬運行100 d 后，區(qū)域中的壓力及孔隙率如圖6 所示，速度及剪切應力如圖7 所示，最大主應力及最小主應力如圖8 所示，x 方向位移及y 方向位移如圖9 所示.

對比2 種方法的結(jié)果，可以看到二者幾乎相同. 另外，圖10 展示了2 種方法的計算時間，可見2種方法所需計算時間與模擬天數(shù)線性相關(guān)，但二重網(wǎng)格混合有限元法的效率明顯高于牛頓迭代法.

6 結(jié)論

本文提出用二重網(wǎng)格混合有限元法來求解多孔介質(zhì)單相流耦合地應力問題，此類問題涉及的模型為單相流耦合地應力模型，模擬的時間尺度大，計算量也大，計算耗時非常長. 我們分別用牛頓迭代法和二重網(wǎng)格混合有限元法對同一問題進行求解，并對比2 種方法的結(jié)果和計算時間，發(fā)現(xiàn)二重網(wǎng)格混合有限元方法的效率明顯高于牛頓迭代法.

參考文獻：

［1］ Gudala M， Govindarajan S K. Numerical modellingof coupled single-phase fluid flow and geomechanicsin a fractured porous media ［J］. J Petrol Sci Eng，2020， 191： 107215.

［2］ Gao C， Gray K E. A coupled geomechanics and res?ervoir simulator with a staggered grid finite differencemethod［ J］. J Petrol Sci Eng， 2022， 209： 109818.

［3］ Shojaei M J. Pore-scale dynamics of foam flow in porousmedia［ D］. Manchester： University of Manchester，2020.

［4］ Li J， Wang Y， Wang H. A study on the rheologicalmodel of porous media in deep rock mass ［J］. RockSoil Mech， 2008， 29： 2355.

［5］ Chen W， Pei W， Li S， et al. Numerical model andengineering application of fluid structure coupling insaturated porous media ［J］. Chin J Rock Mech Eng，2013（z2）： 3346.

［6］ Zidane A， Firoozabadi A. An implicit numericalmodel for multicomponent compressible two-phaseflow in porous media ［J］. Adv Water Res， 2015，85： 64.

［7］ Zhang N， Yao J. Research on multiscale mixed finiteelement numerical method for compressible fluid flow［J］. Chin J Comput Mech， 2017， 34： 226.

［8］ Li Y， Qu X， Mao Y. Finite element simulation ofseepage pressure distribution in complex boundary oilreservoirs ［J］. J Shandong Inst Petrol Chem Technol，2007， 21： 1.

［9］ Yu J， Ren Y， Cao W， et al. Extended mixed elementmethod for compressible porous media flow［J］. J Shandong Univ（ Nat Sci Ed）， 2017， 52： 25.

［10］ Zhang J， Han H， Yu Y， et al. A new two-gridmixed finite element analysis of semi-linear reaction–diffusion equation ［J］. Comput Math Appl， 2021，92： 172.

［11］ Rui H， Liu W. A two-grid block-centered finite differencemethod for Darcy-Forchheimer flow in porousmedia ［J］. SIAM J Numer Anal， 2015， 53：1941.

［12］ Hu H. Two-grid method for compressible miscibledisplacement problem by mixed finite element methodsand expanded mixed finite element method ofcharacteristics ［J］. Numer Algorithms， 2022，89： 611.

［13］ Gao M. A double grid finite element method for solvingnonlinear porous elastic problems ［D］. Jinan：Shandong University， 2022.

［14］ Rutqvist J， Tsang C F. A study of caprock hydromechanicalchanges associated with CO2-injection into abrine formation［ J］. Environ Geol， 2002， 42： 296.

（責任編輯： 周興旺）

基金項目： 國家自然科學基金（11971337）