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

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

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

國家能源局再次強調要持續大力推進天然氣“產供儲銷”體系建設，儲氣庫作為其中重要一環備受關注?！笆奈濉币巹濐A計2025 年建成350×108m3調峰及儲備能力，但目前中國儲氣庫建設優質資源缺乏，因此，油藏被納入建庫選址范疇。

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

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

1 工區概況

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

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

2 儲氣庫滲流-溫度耦合數學模型

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

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

多孔介質中流體質量守恒方程：

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

考慮高速非達西效應的流動方程采用經典的Forchheimer方程：

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

地層流體密度是壓力與溫度的函數，表達式為：

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

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

流體飽和度方程：

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

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

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

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

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

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

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

飽和度和毛細管壓力之間關系：

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

通常使用LEVERETT[12]提出的J函數：

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

初始條件：

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

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

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

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

式中：Λm為流體導熱系數，單位W/（m·K）；Λr為巖石導熱系數，單位W/（m·K）。

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

3 流體高速非達西及高壓物性實驗

3.1 高速非達西實驗

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

表1 紊流因子測試實驗數據Table 1 Experimental data for turbulence factor

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

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

式中：M為氣體分子質量，單位g/mol；A為實驗樣品橫截面積，單位cm2；p1為實驗樣品入口壓力，單位MPa；p2為實驗樣品出口壓力，單位MPa；z為氣體壓縮因子；R為通用氣體常數，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）可以表達為式（18）。

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

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

3.2 流體黏溫關系

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

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

4 復雜斷塊油藏儲氣庫數值模擬

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

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

4.1 生產動態歷史擬合

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

4.2 冷氣注入對儲氣庫生產動態的影響

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

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

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

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

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

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

5 結論

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

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

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