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

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

中圖分類號： O241. 82 文獻(xiàn)標(biāo)志碼： A DOI： 10. 19907/j. 0490-6756. 2024. 041002

1 引言

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

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

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

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

2 模型

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

準(zhǔn)確模擬多孔介質(zhì)中單相流耦合地應(yīng)力問題需要研究流體流動與地應(yīng)力的耦合問題，以及巖石結(jié)構(gòu)因地應(yīng)力而產(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 所示，對應(yīng)的巖石流體參數(shù)如表1 所示. 假設(shè)該區(qū)域有5 口生產(chǎn)井存在，坐標(biāo)分別為［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 所示，速度及剪切應(yīng)力如圖3 所示，最大主應(yīng)力及最小主應(yīng)力如圖4 所示，x 方向位移及y 方向位移如圖5 所示.

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

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

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

6 結(jié)論

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

參考文獻(xiàn)：

［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.

（責(zé)任編輯： 周興旺）

基金項目： 國家自然科學(xué)基金（11971337）