基于滲流-溫度雙場耦合的油藏型儲氣庫數(shù)值模擬

何海燕,劉先山,耿少陽,孫軍昌,孫彥春,賈倩

（1.中國石油冀東油田公司，河北 唐山 063200；2.成都理工大學(xué)能源學(xué)院，四川 成都 610059；3.東北石油大學(xué)環(huán)渤海能源研究院，秦皇島 河北 066004）

國家能源局再次強(qiáng)調(diào)要持續(xù)大力推進(jìn)天然氣“產(chǎn)供儲銷”體系建設(shè)，儲氣庫作為其中重要一環(huán)備受關(guān)注。“十四五”規(guī)劃預(yù)計2025 年建成350×108m3調(diào)峰及儲備能力，但目前中國儲氣庫建設(shè)優(yōu)質(zhì)資源缺乏，因此，油藏被納入建庫選址范疇。

中國大陸經(jīng)歷多期次構(gòu)造運動，其東部形成一系列復(fù)雜斷塊，中西部形成褶皺構(gòu)造，導(dǎo)致建庫油藏斷塊多、非均質(zhì)性強(qiáng)、邊底水及人工注水侵入地層后流體關(guān)系復(fù)雜[1-2]。因此，目前無論是已建庫的京58 儲氣庫，還是正在建庫的冀東油田南堡1-29儲氣庫，都是復(fù)雜斷塊油藏型儲氣庫。復(fù)雜斷塊油藏改建儲氣庫后，多周期高速注采過程均為油氣水三相流動，存在注采周期短、氣體流速高、壓力波及范圍小等特征。注氣周期注入的冷氣會擾動儲層溫度場，油氣水的黏度、氣油比等高壓物性參數(shù)受溫度影響十分嚴(yán)重。在油藏衰竭開發(fā)階段，有較多學(xué)者針對注氣、水、壓裂液等對儲層產(chǎn)生的溫度場擾動開展了研究。王增林等[3]使用Fluent 軟件模擬均質(zhì)儲層條件下使用不同管柱的油藏溫度場變化。鄭少婧等[4]通過實驗探究了儲氣庫交替注采工況下儲層滲透率溫度敏感性，基于實驗結(jié)果建立了考慮滲透率溫度敏感性的氣井產(chǎn)能方程。郭肖等[5]應(yīng)用熱應(yīng)力理論，推導(dǎo)了滲透率隨溫度變化的理論模型，并采用變圍壓、變內(nèi)壓應(yīng)力敏感實驗進(jìn)行了驗證。前人通過實驗和理論證實了儲層中溫度場被擾動后，會對巖石滲透率、彈性模量、泊松比及流體黏度、體積系數(shù)等高壓物性產(chǎn)生較大影響。此外，不同于氣藏衰竭開發(fā)階段中流體的低速流動只需克服黏滯阻力產(chǎn)生的壓降，儲氣庫運行階段采氣速度是氣藏衰竭開發(fā)階段的20~30 倍，井筒附近天然氣高速流動產(chǎn)生的慣性力附加壓降已不容忽視。EL-ZEHAIRY 等[6]基于XCT 數(shù)據(jù)研究了多孔介質(zhì)微觀非均質(zhì)性對慣性流的影響，研發(fā)了孔隙網(wǎng)絡(luò)模型（PNM）來模擬多孔介質(zhì)中非達(dá)西流動，與均勻多孔介質(zhì)相比，非均質(zhì)多孔介質(zhì)中由于連接較差，孔喉中存在更多停滯區(qū)域，減少了流體流動的有效面積，樣品內(nèi)從停滯區(qū)到吼道連接處速度分布由低到高，更容易觀察到流體慣性效應(yīng)。WANG 等[7]學(xué)者指出毛細(xì)管力和邊界層效應(yīng)是油藏中非達(dá)西流動產(chǎn)生的主要原因，提出了一種油藏低速注水過程中非達(dá)西流體動力學(xué)表征方法，通過該方法，給定毛細(xì)管力和邊界層厚度即可計算出單個毛細(xì)管內(nèi)流體的流速。NIE 等[8]學(xué)者也研究了流體在未固結(jié)介質(zhì)（中等強(qiáng)度支撐劑填充介質(zhì)）與固結(jié)介質(zhì)（露頭、巖心等）中流動時的非達(dá)西紊流因子測定方法及表征公式。

傳統(tǒng)油藏工程方法和常規(guī)數(shù)值模擬將儲層視為恒溫[9-11]，巖心實驗雖能模擬流體高速非達(dá)西效應(yīng)及溫度變化對儲層巖石和流體高壓物性的影響，但受限于尺度，僅能代表儲層中一個點。然而，在復(fù)雜斷塊油藏中，同一斷塊內(nèi)儲層存在非均質(zhì)性，不同斷塊中儲層具有不同溫壓及流體系統(tǒng)，在改建儲氣庫后更是存在短周期、高氣體流速以及井控范圍小等特征。因此，傳統(tǒng)油藏工程方法、巖心實驗和常規(guī)數(shù)值模擬難以同時刻畫復(fù)雜斷塊油藏型儲氣庫交替注采工況、溫度場擾動及高速非達(dá)西效應(yīng)附加壓降，導(dǎo)致儲氣庫多周期運行過程調(diào)峰、單井注采能力等生產(chǎn)動態(tài)指標(biāo)預(yù)測精度低，最終造成新鉆井?dāng)?shù)量、投資預(yù)測等誤差大。因此，結(jié)合高速非達(dá)西實驗及流體黏溫實驗，重點研究復(fù)雜斷塊油藏型儲氣庫周期注采過程儲層溫度交替變化及高速非達(dá)西附加壓力損失對儲層流體滲流及生產(chǎn)動態(tài)的影響。研究成果可以為復(fù)雜斷塊油藏建庫方案設(shè)計提供理論指導(dǎo)，為儲氣庫安全高效運行奠定基礎(chǔ)。

1 工區(qū)概況

儲氣庫由帶氣頂弱邊底水的飽和油藏改建而成，儲層平均凈毛比為0.35，平均孔隙度為0.08，平均滲透率為19×10-3μm2。儲層自上而下有3 個不同的流體系統(tǒng)，其中氣層中束縛水飽和度為0.28，另外2套油層中束縛水飽和度分別為0.36 和0.4，油層之下為含水層。工區(qū)被斷層分為7個斷塊，各斷塊氣油水界面不統(tǒng)一，流體分布（圖1）較為復(fù)雜。

圖1 不同斷塊中油氣水3相初始分布Fig.1 Initial distribution of oil,gas and water in different fault blocks

2 儲氣庫滲流-溫度耦合數(shù)學(xué)模型

為描述復(fù)雜斷塊油藏型儲氣庫儲層中低溫天然氣注入后流體高壓物性變化，以及溫度場擾動對滲流場的影響，基于多孔介質(zhì)中滲流-溫度耦合作用機(jī)理，建立了滲流-溫度耦合數(shù)學(xué)模型。模型基本假設(shè)條件如下：①枯竭油藏型儲氣庫中構(gòu)造由7個斷塊組成，采用笛卡爾網(wǎng)格剖分；②儲層中流體由油氣水三相組成，密度、黏度是溫度和壓力的函數(shù)；③天然氣在儲層中高速流動時，除了受孔喉壁面及流體間黏滯阻力影響外，還受高速慣性力產(chǎn)生的附加壓降影響。

描述了儲氣庫多孔介質(zhì)中油氣水三相流動的控制方程，主要包含質(zhì)量守恒方程和能量守恒方程。每一組方程都有3個相的子集，加在一起就組成了完整控制方程。質(zhì)量守恒方程和能量守恒方程適用于所有相。因此，控制方程可以簡化為單相形式。

多孔介質(zhì)中流體質(zhì)量守恒方程：

式中：φ為儲層孔隙度；ρm為流體密度，單位kg/m3，m表示油氣水三相；t為時間，單位s；v→m為流體速度，單位m/s；qm為流體產(chǎn)量，單位m3/s。

考慮高速非達(dá)西效應(yīng)的流動方程采用經(jīng)典的Forchheimer方程：

式中：p為儲層壓力，單位MPa；μm為流體黏度，單位mPa·s；k為滲透率，單位10-3μm2；βm為地層流體體積系數(shù)，單位m3/m3。

地層流體密度是壓力與溫度的函數(shù)，表達(dá)式為：

式中：pi為參考壓力，單位MPa；Ti為參考溫度，單位K；ρi為參考密度，單位kg/m3；pm為儲層中油氣水的壓力，單位MPa；Tm為儲層中油氣水的溫度，單位K；ρm為儲層壓力pm和儲層溫度Tm對應(yīng)的油氣水的密度，單位kg/m3；ρi為參考壓力pi和參考溫度Ti對應(yīng)的油氣水的密度，單位kg/m3；CL為流體彈性壓縮系數(shù)；αL為流體熱膨脹系數(shù)。

式中：μm為儲層壓力pm和儲層溫度Tm對應(yīng)的油、氣、水的黏度，單位mPa·s；μi為參考壓力pi和參考溫度Ti對應(yīng)的油、氣、水的黏度，單位mPa·s；γm與ηm均為儲層中流體的黏度方程系數(shù)。

流體飽和度方程：

式中：nm為流量數(shù)量；Sm為儲層流體飽和度。

將毛細(xì)管壓力定義為單位面積的力，可以表示為：

式中：pc為毛細(xì)管壓力，單位MPa；r為毛細(xì)管半徑，單位m；σ為表面張力，單位N/m；θ為接觸角，（°）。

毛細(xì)管壓力與重力作用過程，向上與向下的力是平衡的，可以得到：

式中：g為重力加速度，單位m/s2；h為潤濕相在毛細(xì)管中上升的高度，單位m。

油氣間毛管壓力可以表示為：

式中：pcgo為油氣間的毛細(xì)管壓力，單位MPa；pg為氣相壓力，單位MPa；po為油相壓力，單位MPa。

飽和度和毛細(xì)管壓力之間關(guān)系：

式中：pcow為油水間的毛細(xì)管壓力，單位MPa；pw為油相壓力，單位MPa；Sw為含水飽和度。

通常使用LEVERETT[12]提出的J函數(shù)：

為了模擬多相流動，需要確定儲層內(nèi)初始飽和度場分布。如果我們知道油水界面的位置，就可以通過結(jié)合式（7）和式（10）來確定地層中飽和度隨深度的分布：

初始條件：

式中：z為縱向距離，單位m；z0為參考垂向距離，單位m；p0為參考垂向距離z0對應(yīng)的儲層壓力，單位MPa。

為了表征流體與流體、流體與巖石間的熱量交換關(guān)系，引入能量守恒方程：

式中：hm為流體的焓，單位J；λm為流度系數(shù)，單位10-3μm2/（mPa·s），其中，λm=k·kr,m/μm，kr,m為儲層流體相對滲透率，單位10-3μm2；（ρU）eff為單位質(zhì)量的有效內(nèi)能，單位J；Λeff為有效導(dǎo)熱系數(shù)，單位W/（m·K）。

式中：Um為流體的內(nèi)能，單位J；Ur為巖石的內(nèi)能，單位J；ρr為巖石密度，單位kg/m3。

式中：Λm為流體導(dǎo)熱系數(shù)，單位W/（m·K）；Λr為巖石導(dǎo)熱系數(shù)，單位W/（m·K）。

基于MATLAB 軟件中油藏數(shù)值模擬工具箱（Matlab Reservoir Simulation Toolbox，簡稱MRST），非線性方程離散采用有限體積法（Finite Volume Method，簡稱FVM），在空間上采用兩點通量近似（Two-Point-Flux-Approximation，簡稱TPFA）有限體積格式，在時間上采用后向（隱式）歐拉格式對方程進(jìn)行耦合離散求解[13]，井模型采用Peaceman 模型[14]進(jìn)行離散求解。

3 流體高速非達(dá)西及高壓物性實驗

3.1 高速非達(dá)西實驗

結(jié)合Forchheimer 方程（式16）可擬合紊流因子β及滲透率k[15-17]，為獲取紊流因子設(shè)計了此次驅(qū)替實驗。實驗主要采用智能驅(qū)替模擬系統(tǒng)，選取N 儲氣庫代表性巖心，分別設(shè)計采用0.50～0.95 MPa 共10組驅(qū)替壓力，測試不同驅(qū)替壓力下驅(qū)完2 mL 氣體的氣體流速（測試體積/測試時間）。實驗氣體為氮氣，測試氮氣體積為2 mL，實驗巖心長度為5.12 cm，直徑為2.53 cm，驅(qū)替環(huán)壓為10 MPa，實驗結(jié)果見表1。

表1 紊流因子測試實驗數(shù)據(jù)Table 1 Experimental data for turbulence factor

式中：p為壓力，單位MPa；X為流體流動方向；μ為氣體黏度，單位mPa·s；v為流體流速，單位cm/s；k為滲透率，單位10-3μm2；β為紊流因子，單位108/m；ρ為流體密度，單位g/cm3。

在此次驅(qū)替實驗條件下，F(xiàn)orchheimer 方程可以表達(dá)為式（17）。以MA（p12-p22）/(2zRTμlρpQp)為縱坐標(biāo)，ρpQp/（μA）為橫坐標(biāo)，可以擬合出一條直線，該直線截距為1/k，斜率為β，即紊流因子。

式中：M為氣體分子質(zhì)量，單位g/mol；A為實驗樣品橫截面積，單位cm2；p1為實驗樣品入口壓力，單位MPa；p2為實驗樣品出口壓力，單位MPa；z為氣體壓縮因子；R為通用氣體常數(shù)，R=8.314 472 m3·Pa/（K·mol）；T為溫度，單位K；l為實驗樣品長度，單位cm；ρp為泵中流體密度，單位g/cm3；Qp為泵中流體流量，單位cm3/h。

令式（17）中MA（p12-p22）/(2zRTμlρpQp)為y，ρpQp/（μA）為x，式（17）可以表達(dá)為式（18）。

通過紊流因子巖心實驗，以MA（p12-p22）/（2zRTμlρpQp）為縱坐標(biāo)，ρpQp/μA為橫坐標(biāo)，擬合繪制直線（圖2）的斜率（即紊流因子β）為1.63×108/m。

圖2 紊流因子巖心實驗擬合Fig.2 Turbulence factor regression curve base on experiment data

3.2 流體黏溫關(guān)系

實驗室一般使用黏度計或旋轉(zhuǎn)式流變儀來測量原油的黏度[18]。從N 油藏儲氣庫儲層中獲取原油樣品，采用自動密度黏度測定儀，測定剪切速率為60 s-1時，不同溫度（10.13～90.04 ℃）條件下原油黏度，實驗結(jié)果見圖3。氣水黏度計算通常采用4 種方法：Lohrentz-Brey-Clark（LBC 方法）、PFCT 方法、SUPERTRAPP 方法、Vesovic-Wakeham（VW 方法）。研究表明SUPERTRAPP 方法誤差較小[19-21]。因此，儲層壓力（22.6 MPa）條件下的油氣水黏度與溫度關(guān)系基于SuperTrapp軟件數(shù)據(jù)包，結(jié)合SuperTrapp模型計算獲得（圖3）。

圖3 油氣水黏溫關(guān)系（壓力：22.6MPa）Fig.3 Viscosity-temperature relationship curves of oil,gas and water

4 復(fù)雜斷塊油藏儲氣庫數(shù)值模擬

結(jié)合前期油氣藏數(shù)值模擬、油氣藏改建儲氣庫數(shù)值模擬方法[22]，建立了考慮儲氣庫多輪注采過程冷氣注入擾動儲層溫度場及高速非達(dá)西效應(yīng)等機(jī)理模擬的油藏型儲氣庫數(shù)值模擬方法和技術(shù)流程（圖4）。依據(jù)流體高壓物性隨溫度變化實驗及高速紊流實驗，獲取流體高壓物性參數(shù)隨溫度的函數(shù)關(guān)系及高速非達(dá)西紊流因子，分別用于儲氣庫注氣周期冷氣注入過程溫度場擾動數(shù)值模擬以及對流動方程（式2）中慣性力產(chǎn)生的附加壓降進(jìn)行校正。數(shù)值模擬采用笛卡爾網(wǎng)格剖分，網(wǎng)格數(shù)約為94萬，X、Y、Z方向平均步長分別為40 m×40 m×4 m。

圖4 油藏型儲氣庫數(shù)值模擬技術(shù)流程Fig.4 Workflow of numerical simulation for UGS rebuilt from oil reservoir

4.1 生產(chǎn)動態(tài)歷史擬合

建模階段屬性模型是通過地質(zhì)統(tǒng)計學(xué)數(shù)據(jù)分析、插值生成的，存在較強(qiáng)的不確定性。因此，需結(jié)合衰竭開發(fā)階段生產(chǎn)動態(tài)監(jiān)測資料來反演井間儲層物性參數(shù)。首先對區(qū)塊油、氣等進(jìn)行擬合（圖5a、圖5b），從而保證區(qū)塊物質(zhì)平衡，然后檢查單井瞬時產(chǎn)量擬合情況，在此基礎(chǔ)上依次開展單井靜壓、流壓擬合（圖5c、圖5d）。

4.2 冷氣注入對儲氣庫生產(chǎn)動態(tài)的影響

基于歷史擬合的模型，選取主力斷塊5口儲氣庫井，單井平均以11.2×104m3/d 注氣200 d，20×104m3/d 采氣120 d，注采平衡期均為15 d，設(shè)計考慮溫度場及高速非達(dá)西效應(yīng)影響的方案。其中模型的注入氣溫度為25 ℃，儲層中深溫度為87.8 ℃，儲層初始溫度梯度為3 ℃/hm，流體黏度隨溫度變化規(guī)律依據(jù)黏溫實驗數(shù)據(jù)（圖3）賦值，紊流因子依據(jù)非達(dá)西實驗結(jié)果賦值為1.63×108/m。研究將流體黏溫、高速紊流實驗結(jié)果敷設(shè)于模型里，把不考慮溫度場及非達(dá)西效應(yīng)的方案設(shè)為基礎(chǔ)方案。

模擬結(jié)果表明，相較于基礎(chǔ)方案，注冷氣擾動溫度場方案由于井控范圍內(nèi)儲層溫度下降，導(dǎo)致累產(chǎn)油量（圖6b）下降，累產(chǎn)水量增加（圖6c），而累產(chǎn)油下降幅度小于累產(chǎn)水增加幅度，使地層采出液量增多，地層壓力下降（圖6d）?？紤]非達(dá)西效應(yīng)影響時，一方面，除了黏滯阻力產(chǎn)生壓降外，高速慣性力還多產(chǎn)生了一部分壓降，因此，相同配產(chǎn)配注條件下，天然氣注入后部分采不出，隨著儲氣庫多周期運行，注氣末天然氣儲量及壓力逐漸增加；另一方面，在定產(chǎn)氣量生產(chǎn)條件下考慮非達(dá)西效應(yīng)，需增大生產(chǎn)壓差才能產(chǎn)出相同的天然氣量，因此，多周期運行后油、水累產(chǎn)量增加。總的來說，不考慮冷氣注入溫度場擾動和高速非達(dá)西效應(yīng)，將造成第三采氣期末累產(chǎn)油、氣量分別偏小3.29%、10.52%，累產(chǎn)水量偏大21.07%。

此外，模擬了4 組不同注入氣溫度（10、30、50、70℃）的儲氣庫多周期運行方案。結(jié)果表明，溫度場擾動對氣體滲流的影響微弱，遠(yuǎn)小于高速非達(dá)西效應(yīng)對氣體滲流的影響。但儲層溫度降低后對原油滲流能力降低幅度影響較大（圖7），不考慮溫度場擾動的基礎(chǔ)方案3周期內(nèi)累產(chǎn)油26 201 m3，考慮注入氣溫度為10 ℃時，原油采出量降低4 399 m3，降幅為16.79%。

圖7 注入氣溫度對原油累產(chǎn)的影響Fig.7 Effect of injection gas temperature on oil production cumulative

選取構(gòu)造高部位QK4 井，設(shè)計考慮高速非達(dá)西與不考慮高速非達(dá)西效應(yīng)各6項注采方案（注200 d，采120 d），方案中日采氣速度分別為（10、20、40、60、80、100）×104m3，對應(yīng)的日注氣速度分別為（5.6、11.2、22.4、33.6、44.9、56.1）×104m3。結(jié)果表明，在當(dāng)前物性及配產(chǎn)配注條件下，注氣末井控溫度范圍隨注氣速度呈對數(shù)上升，低速階段上升較快，高速階段由于受儲層物性及有限時率強(qiáng)注強(qiáng)采限制，上升速度變緩；考慮高速非達(dá)西效應(yīng)后，低速注采階段流體慣性力產(chǎn)生的附加壓降較小，對流體滲流影響較小，注氣末井控溫度范圍幾乎重合；隨著注采速度上升，井控溫度范圍逐漸增大，但注氣速度上升到44.9×104m3/d 后，考慮高速非達(dá)西效應(yīng)的方案注氣末井控溫度范圍幾乎恒定，未考慮高速非達(dá)西效應(yīng)的方案仍在增加，這表明注采速度較高時，受非達(dá)西效應(yīng)附加壓降影響，該物性下單井注采能力已達(dá)到極限，而不考慮非達(dá)西效應(yīng)影響的單井注采能力還有較大提升空間（圖8）。因此，基于該極限，可以確定不同儲層物性下受溫度場擾動及高速慣性力影響的單井合理注采氣能力。

圖8 注氣速度對注氣末期井控溫度范圍的影響Fig.8 Effect of gas injection rate on well control temperature range at the end of gas injection cycle

5 結(jié)論

1）建立了滲流-溫度雙場耦合數(shù)學(xué)模型，并基于有限體積法（FVM），在空間上采用兩點通量近似方案（TPFA），在時間上采用后向（隱式）歐拉格式對模型進(jìn)行耦合離散求解。該模型考慮了儲氣庫注氣周期冷氣注入擾動溫度場以及高速注采過程中高速非達(dá)西效應(yīng)對流體滲流的影響，更符合儲氣庫的特殊工況。

2）基于實驗回歸了紊流因子與流體黏溫關(guān)系，明確了油氣黏度大幅上升時水黏度幾乎不變，這一特性使得油氣相對滲流能力下降時水相對滲流能力反而上升?；趯嶒灁?shù)據(jù)及建立的模型，開展了儲氣庫多周期運行數(shù)值模擬。模擬實例表明，溫度場擾動、高速非達(dá)西效應(yīng)分別是累產(chǎn)油、氣量誤差的主控因素。儲層溫度場擾動減少的產(chǎn)油量比高速非達(dá)西效應(yīng)增加的產(chǎn)油量多。不考慮冷氣注入溫度場擾動和高速非達(dá)西效應(yīng)，將造成第三采氣期末的累產(chǎn)油、氣量分別偏小3.29 %、10.52 %，累產(chǎn)水量偏大21.07%。

3）生產(chǎn)動態(tài)敏感性分析結(jié)果表明，冷氣注入造成井控范圍內(nèi)儲層溫度下降，該范圍隨注氣速度增加呈對數(shù)上升，低速階段上升較快，高速階段由于受儲層物性及有限時率強(qiáng)注強(qiáng)采限制，上升速度變緩；高速慣性力附加壓力損失使相同配產(chǎn)配注條件下天然氣注入后部分采不出，隨著儲氣庫多周期運行，注氣末天然氣儲量及壓力逐漸增加，定天然氣量生產(chǎn)所需生產(chǎn)壓差增大，油水累產(chǎn)量隨之增加。