基于有限體積法的三維油氣運聚模擬技術

郭秋麟，陳寧生，謝紅兵，吳曉智，劉繼豐，趙錫然，高日麗，胡俊文

（中國石油勘探開發(fā)研究院）

基于有限體積法的三維油氣運聚模擬技術

郭秋麟，陳寧生，謝紅兵，吳曉智，劉繼豐，趙錫然，高日麗，胡俊文

（中國石油勘探開發(fā)研究院）

針對三維油氣運聚模擬技術地質模型過于簡化所導致的應用效果難于達到實際需求的問題，建立順層柱狀PEBI網(wǎng)格動態(tài)三維地質模型，解決不同沉積相內儲集層非均質性問題，降低混合巖性對模擬精度的影響。構建基于有限體積法的三維數(shù)值模擬技術，包括變網(wǎng)格滲流方程構建、全張量滲透率的分解與計算、傳導率的計算等。通過對關鍵參數(shù)曲線的光滑、網(wǎng)格流動上下游穩(wěn)定性處理等，提高牛頓法迭代收斂性能。采用自動調整時間步長與多核并行計算的技術，提高軟件運行效率。模擬計算了渤海灣盆地南堡凹陷油氣運聚過程、不同時期地層含油飽和度與油氣資源豐度，模擬的含油飽和度、油氣聚集量和分布位置符合當前勘探現(xiàn)狀，揭示了目標層資源探明率為84.4%，主要資源分布在B區(qū)和C區(qū)，預測未發(fā)現(xiàn)資源主要分布在C區(qū)NP2-16井附近和A區(qū)NP5-4井北側的構造上。圖6表7參15

三維油氣運聚模擬技術；有限體積法；順層柱狀PEBI網(wǎng)格；變網(wǎng)格滲流方程；南堡凹陷

0 引言

目前，油氣運聚模擬主要有3種方法，即二維流線法（flowpath）[1]、侵入逾滲法（invasion percolation）[1-2]和三維多相達西流法（multi-phase Darcy flow）[1,3-6]。其中，流線模擬法適用于二維構造面上的油氣運聚模擬，僅能模擬構造型油氣藏的運聚；侵入逾滲法主要用于模擬油氣運聚路徑，既可在二維空間也可在三維空間使用。以上兩種方法都要求地質模型是靜態(tài)的，模擬網(wǎng)格不變。三維多相達西流法是各種運聚定量模擬技術中考慮因素最全面、技術較成熟的方法[7]，其有3種核心算法：有限元法（如PetroMod、3D SEMI）、有限體積法（如Temispack）和有限差分法（如BasinMod等）。每種方法采用的三維網(wǎng)格有所差異，有限差分法僅適用于規(guī)則的中心網(wǎng)格，如矩形網(wǎng)格；有限元法適用于規(guī)則或不規(guī)則的角點網(wǎng)格，如矩形網(wǎng)格、角點網(wǎng)

格、四面體網(wǎng)格等；有限體積法適用于規(guī)則或不規(guī)則中心網(wǎng)格，如矩形網(wǎng)格、PEBI網(wǎng)格（垂直正交網(wǎng)格）等。各種算法及相應的網(wǎng)格建模技術各有優(yōu)缺點。隨著地質認識的深入和油氣勘探的發(fā)展，對三維地質模型的要求越來越高，建模中較簡單且較常用的矩形網(wǎng)格已很難滿足復雜地區(qū)的建模需要。由于PEBI網(wǎng)格建模技術更加靈活，因此適用范圍也更寬。

基于有限體積法的三維油氣運聚模擬技術在國內外已開展了許多研究。2001年，馮勇等[8]研究了PEBI網(wǎng)格和有限體積法相結合的方法，但應用效果不明顯；2003年，石廣仁等[3]對該方法進行了改進；2009年，Hantschelh T和Kauerauf A I[1]對該技術進行了較深入的研究；IBM公司W(wǎng)atson實驗室[4]提出了一種三維控制體積有限元法；2010年，石廣仁等[5]發(fā)展了基于PEBI網(wǎng)格的有限體積法，并在庫車坳陷進行了應用，取得了初步應用實效。以上研究的地質建模網(wǎng)格除IBM公司外均為水平柱狀PEBI網(wǎng)格，即垂向上的網(wǎng)格面與水平面平行。這類柱狀網(wǎng)格的頂?shù)酌媾c地層面相交，在地質上稱為穿層或穿時。這樣劃分網(wǎng)格雖可以提高運算速度，但卻損失了建模精度，不利于復雜地區(qū)油氣運聚精細模擬。

本文從地質模型的建立、滲流方程的構建、傳導率的全張量計算、牛頓法迭代穩(wěn)定性與計算效率的提高等方面，研發(fā)了基于有限體積法的三維油氣運聚模擬技術，包括：①建立了順層柱狀PEBI網(wǎng)格三維動態(tài)地質模型，精細刻畫地層的演化，初步解決了地層非均質性、斷層等引起的滲流特定性及混合巖性等地質難題；②構建變網(wǎng)格條件下的滲流方程，代替定網(wǎng)格滲流方程，更有效地實現(xiàn)了質量守恒；③引入了矢量滲透率（即全張量滲透率），解決復雜的滲流問題。該技術成功應用于南堡凹陷，取得了良好效果。

1 網(wǎng)格劃分與關鍵地質問題處理

1.1 三維地質體網(wǎng)格劃分方法

1.1.1 順層柱狀PEBI網(wǎng)格

PEBI網(wǎng)格中心點與所有相鄰網(wǎng)格中心點的連線均垂直通過相應的網(wǎng)格面，這種網(wǎng)格劃分方法可以提高滲流模擬速度。目前一般采用水平柱狀PEBI網(wǎng)格，即平面為PEBI網(wǎng)格，垂向為水平網(wǎng)格，網(wǎng)格面與地層面相交（“穿時”），地質模型精度受影響。本文采用順層柱狀PEBI網(wǎng)格，即平面為PEBI網(wǎng)格，垂向為地層面網(wǎng)格，這種網(wǎng)格在平面上能夠根據(jù)已知數(shù)據(jù)點的分布構建最優(yōu)平面網(wǎng)格，最大化提高模擬運行效率；在垂向上按地層面劃分，保持網(wǎng)格面與地層面一致，避免了“穿時”，提高了地質模型精度。

1.1.2 動態(tài)網(wǎng)格

常規(guī)模擬網(wǎng)格為靜態(tài)網(wǎng)格[9]，本文在不同地質歷史時期采用動態(tài)網(wǎng)格，即在地層埋藏過程中網(wǎng)格與地層的變化一致，網(wǎng)格內的地層一直不變，僅屬性因壓實等地質作用發(fā)生變化（網(wǎng)格體積變內容不變）。動態(tài)網(wǎng)格的優(yōu)點是地質模型更接近實際，缺點是網(wǎng)格體積受地層厚度變化影響大，降低計算速度。

1.2 關鍵地質問題的處理

1.2.1 地層物性非均質性處理方法

筆者采用隨機抽樣解決非均質性問題。首先按沉積相類型統(tǒng)計孔隙度、滲透率和孔喉半徑等參數(shù)的最大值、最小值、均值和方差，并建立分布模型；然后按分布模型隨機抽樣獲得不同網(wǎng)格的孔隙度、滲透率和孔喉半徑等參數(shù)。

1.2.2 斷層、河道等特定方向的滲透率處理方法

首先，將滲透率分為3個方向：主方向滲透率（Kx，河道、斷層的走向）、副方向滲透率（Ky）和垂向滲透率（Kz）；然后，用矢量或方位角表示主方向，副方向與主方向垂直，垂向是指垂直地層面的方向（地層傾向）；最后，將河道、斷層帶等所在的網(wǎng)格分別賦值（可采用隨機抽樣方法），包括Kx、Ky、Kz和矢量（或方位角）。

1.2.3 混合巖性的處理方法

油氣運移模擬網(wǎng)格比較大，網(wǎng)格內一般包含多種巖性，在網(wǎng)格參數(shù)賦值時需要特殊處理。本文通過引入有效儲集層比例參數(shù)，設置有效儲集層存在的下限值等處理方法，降低混合巖性對模擬精度的影響程度。以砂泥巖層為例，處理方法如下：設定每個網(wǎng)格的有效儲集層比例（f），給定有效儲集層存在的下限值（fmin）。當f >fmin時，該網(wǎng)格為“儲集層或輸導層”，網(wǎng)格物性取本網(wǎng)格的砂巖物性，有效孔隙空間為砂巖的孔隙空間；反之，屬于“非滲透層”網(wǎng)格，物性取本網(wǎng)格的泥巖物性。

2 基于有限體積法的三維數(shù)值模擬技術

2.1 變網(wǎng)格滲流方程

假設網(wǎng)格體運動速度與流體流動速度相比可以忽略不計，建立變網(wǎng)格滲流方程。

2.1.1 質量守恒方程

根據(jù)質量守恒定律，任意控制體?中，源產生的質量減去流出的質量等于流體質量的增量。在多相滲流的情況下，控制體?內，l相流體滲流的質量守恒方程為：

上式即為變網(wǎng)格滲流方程，V=V（t）是與時間有關的變量，下文相關方程的理論推導均以（1）式為基礎。

（1）式對時間積分得：

2.1.2 滲流運動方程

為了更好地描述地層非均質性，筆者采用全張量的滲透率；為了解決滲透率張量、網(wǎng)格界面法向量、流體勢梯度向量3者不正交的問題，采用了達西定律的全張量形式：

2.1.3 流動方程

三相滲流數(shù)學模型如下。水相流動方程為：

油相流動方程為：

氣相流動方程為：

油水系統(tǒng)中的毛細管力方程為：

油氣系統(tǒng)中的毛細管力方程為：

輔助方程為：

2.2 初始條件和邊界條件

邊界條件分兩類，即封閉邊界和流動邊界。在封閉邊界外，網(wǎng)格的流體屬性對流動無影響，在實際計算中忽略該類邊界面即可，邊界外網(wǎng)格的變量及流體屬性無須更新。流動邊界按定壓邊界處理，對于邊界外的網(wǎng)格，壓力為靜水壓力，油氣飽和度為0。

2.3 全張量滲透率的分解與計算

在區(qū)域坐標系O-xyz中，設壓力p下地層的滲透率張量為：

則壓力p下的滲透率張量K可以變換為：

根據(jù)滲透率與壓力的關系：

2.4 傳導率的計算

假設網(wǎng)格i,j是相鄰網(wǎng)格，使用單點上游法時，l相流體的傳導率計算公式如下：網(wǎng)格i中心到網(wǎng)格i,j間界面的傳導率為：

網(wǎng)格j中心到網(wǎng)格i,j間界面的傳導率為：

要求Sij與nij、Sji與nji位于界面的同側且點積大于0。

網(wǎng)格i,j間l相流體的流量為：

網(wǎng)格i,j間l相流體的勢差為：

對網(wǎng)格間勢差的修正系數(shù)，只能為正值。

3 提高穩(wěn)定性與計算效率的處理方法

3.1 特殊處理以提高牛頓法的收斂性能

①參數(shù)曲線光滑性處理。為提高收斂性能，對相對滲透率曲線、PVT曲線以及存在啟動壓力梯度的滲流曲線等進行光滑處理。②上下游穩(wěn)定性處理。在1個時間步內，要求網(wǎng)格間上下游關系不變，以確保雅可比矩陣的穩(wěn)定，從而提高運行速度。③網(wǎng)格體積異常的處理。體積為0的網(wǎng)格，在模擬計算中需要進行網(wǎng)格鄰接關系處理，以消除0體積網(wǎng)格變量，避免導致雅可比矩陣奇異。小網(wǎng)格是指孔隙體積較小的網(wǎng)格，是影響收斂性的關鍵因素之一。采取飽和度上游化、壓力平均化等處理方法，確保矩陣方程的求解效率和精度。

3.2 自動調整時間步長以提高計算速度

由于收斂性也依賴于初始解的近似程度，要求時間步長不能過大；為了降低數(shù)值模擬過程總計算量，在保證計算精度的前提下，盡量采用較大時間步長。為解決這一矛盾，根據(jù)預先設定的目標迭代次數(shù)和最大迭代次數(shù)來自動調整時間步長。圖1為時間步長動態(tài)變化實例。

圖1 自動調整時間步長實例

3.3 多核并行計算以提高計算速度

網(wǎng)格內計算并行設計按CPU的緩存單元大小的整倍數(shù)來分配并行單元規(guī)模即可。網(wǎng)格間的并行設計分兩個步驟。第1步，網(wǎng)格層內計算并行化（見圖2a）：同層之內，橫向流動，只跟相鄰兩個網(wǎng)格單元有關，因此可分層并行計算，但考慮到偽共享問題，采用隔層并行計算，避免緩存沖突。首先并行層內線程-1-1、層內線程-1-2、層內線程-1-3等，其次并行層內線程-2-1、層內線程-2-2、層內線程-2-3等。第2步，網(wǎng)格層間計算并行化（見圖2b）：相鄰兩層之間垂向流動只跟上下層間相鄰兩個網(wǎng)格單元有關，因此可采用分層面并行。首先并行層間線程-1-1、層間線程-1-2、層間線程-1-3等，其次并行層間線程-2-1、層間線程-2-2等。

圖2 多核并行計算示意圖

4 實例分析

4.1 模擬區(qū)地質背景

渤海灣盆地南堡凹陷為發(fā)育在奧陶系、石炭-二疊系和中生界基底之上的第三紀沉積凹陷，第三系沉積巖厚度最大可達7 500 m，包括古近系沙河街組（Es）、東營組（Ed）以及新近系館陶組（Ng）、明化鎮(zhèn)組（Nm）

等，面積1 932 km2，其中灘海面積1 100 km2。南堡凹陷是一個具有北斷南超特征的典型箕狀凹陷（見圖3），北區(qū)發(fā)育柳贊、高尚堡披覆背斜構造帶、拾場次凹；南區(qū)發(fā)育老爺廟、北堡背斜構造帶；灘海區(qū)域發(fā)育南堡1號—5號5個構造帶；中央地帶發(fā)育林雀次凹，該次凹是最重要的生烴中心[10-11]。

4.2 目標層段特征

選擇Es—Ed作為模擬目標層。其中，烴源巖層為沙河街組和東三段，儲集層為東營組一段和二段，蓋層為館陶組的火山巖。

4.2.1 探明儲量分布

原油主要聚集在東營組一段，該段探明石油地質儲量3.46×108t，主要分布在凹陷西南部，包括南堡1號、北堡構造、南堡2號和老堡南1構造（見表1、圖3）；已發(fā)現(xiàn)約1.2×1011m3天然氣，以油溶氣為主，主要分布在南堡1號構造、南堡2號構造和北堡構造[11]。

4.2.2 生烴潛力

發(fā)育沙河街組三段（Es3）、沙一段（Es1）和東營組三段（Ed3）等3套烴源巖，其中Es3為主力烴源巖[12-13]。Es3巖性為暗色泥巖、頁巖，有機質類型為Ⅱ型，其中H/C和O/C的平均值分別為1.15和0.14，平均有機碳含量為1.87%、生烴潛量（S1+S2）為6.56 mg/g、氯仿瀝青“A”含量為0.15%、總烴含量7.33×10-4?？傮w評價為優(yōu)質烴源巖，處于成熟階段。

圖3 南堡凹陷構造簡圖

表1 東營組一段石油探明儲量在各構造的分布比例

4.2.3 沉積特征與儲集層物性

目標層段以河流—湖泊沉積體系為主。其中，沙河街組下部發(fā)育厚層較深湖—深湖相暗色泥巖、頁巖，為該區(qū)油氣的大量生成奠定了物質基礎；中上部沙一段是一套辮狀河三角洲前緣相帶的砂泥巖地層，砂層較為發(fā)育[13-15]。東營組以扇三角洲—河流沉積體系為主，發(fā)育優(yōu)質儲集層；館陶組和明化鎮(zhèn)組以辮狀河沉積為主。東營組砂巖孔隙度、滲透率與其所處的沉積相帶有關（見表2、表3）。

4.2.4 蓋層分布

直接覆蓋在東營組之上的館陶組火山巖是主要蓋

層?；鹕綆r巖性以玄武巖為主，凝灰?guī)r為輔，發(fā)育程度具明顯的不均勻性，呈現(xiàn)東北薄、西南厚的分布特點。在1號構造西南部厚度最大，累計厚度達250～500 m；在東南部2號構造附近，火山巖厚度明顯減薄，一般小于50 m；在北堡地區(qū)，火山巖厚度在50～300 m；在老爺廟地區(qū)，火山巖厚度在50～100 m；高尚堡、柳贊地區(qū)火山巖厚度明顯變薄，僅為數(shù)米。蓋層分布控制了東一段油藏的分布。

表2 南堡凹陷不同沉積相帶東營組砂巖孔隙度

表3 南堡凹陷不同沉積相帶東營組砂巖滲透率

4.2.5 斷層作用

南堡凹陷斷裂比較發(fā)育，斷層斷距較大，活動時間長，具有同生性。斷層的作用有兩點：一是可形成斷塊圈閉；二是作為通道，溝通油源。將目標層段按斷層切穿烴源層的程度分為3組，分別為切穿沙三段斷層（F-Es3）、切穿沙一段斷層（F-Es1）和切穿東三段斷層（F-Ed3），各組斷層的滲透率特征見表4。

表4 斷層帶滲透率特征

4.3 三維地質體構建及主要參數(shù)

4.3.1 三維地質體

在平面上，考慮到研究區(qū)最小的圈閉目標約0.4 km2，將研究區(qū)劃分為4 750個PEBI網(wǎng)格；在縱向上，將烴源巖層合并，確定5個模擬層，構成完整的生、儲、蓋組合。該組合由下至上分別為：①烴源層，由東營組三段和沙河街組組成；②輸導層，即東營組二段；③目標層，即東營組一段；④蓋層，即館陶組，所發(fā)育的火山巖具有良好的封蓋性；⑤出水口，即模擬上邊界（見表5）。

4.3.2 地質體參數(shù)和模型參數(shù)

東營組的孔隙度和滲透率參數(shù)來自沉積相帶及對應的物性數(shù)據(jù)（見表2、表3）。其他相關參數(shù)見表5和表6。

表5 模擬層段地層特征

表6 輸入模型的流體PVT數(shù)據(jù)

4.4 模擬結果

模擬結果包括各地質歷史時期目標層含油飽和度、油氣資源豐度和聚集量等關鍵數(shù)據(jù)和圖表。

4.4.1 含油飽和度與石油聚集量模擬結果

模擬時間從距今25 Ma開始至今，按含油飽和度大于10%統(tǒng)計，不同地質時期石油聚集量和聚集區(qū)面積見表7。由表7可見，在距今20 Ma時，石油開始聚集；在距今15 Ma時，已聚集石油7 074×104t，占總聚集量的17.3%，此時石油主要聚集在北部的M10井和GC1井附近及中部的NP2-16井附近（見圖4a）；在距今8 Ma時，已聚集石油28 626×104t，占總聚集量的69.8%，此時在西南部的NP1井、LPN1井附近也出現(xiàn)了石油聚集（見圖4b）；現(xiàn)今，石油總聚集量為41 006×104t，主要聚集在西南部和北部等（見圖4c）。資源豐度圖（見圖5）揭示，石油主要聚集在B區(qū)的NP1井、C區(qū)的LPN1井和NP2-16井附近，北部A區(qū)只有少量聚集，東部D區(qū)沒有聚集。

表7 不同地質時期石油聚集量和聚集區(qū)面積

4.4.2 油氣運移主要路徑追蹤

統(tǒng)計分析模擬網(wǎng)格體各面石油流量（流出為正，流入為負），確定最大流量對應的方向為主要運移路徑。有兩種追蹤方式：①正常追蹤，即從烴源巖層開始追蹤，當流量為0或遇到出水口時，不再追蹤。圖6a為研究區(qū)正常追蹤結果，該圖揭示在東部蓋層不好的位置綠色流線（運移路徑）直接向上進入出水口，說明石油不在目標層聚集，在蓋層與圈閉配合良好的西南部及北部地區(qū)發(fā)生較多的聚集；圖6c為研究區(qū)正常追蹤結果的另一種表示形式，該圖記錄了石油通過側向運移進入到各油藏的過程（藍色流線，即運移路徑）。②反向追蹤，即從聚集區(qū)反追蹤到烴源巖層。圖6b為石油聚集區(qū)反向追蹤結果，該圖揭示主要聚集區(qū)的石油源于下部烴源巖層的生烴中心（紫色部位）及附近。正常追蹤和反向追蹤都有利于油源跟蹤分析和油氣成藏研究，對分析油氣聚集與散失具有重要意義。

圖4 東一段含油飽和度模擬結果

圖5 東一段石油資源豐度模擬結果

4.4.3 模擬結果分析

以下從含油飽和度、模擬聚集量和油氣分布位置3

方面進行分析。①含油飽和度模擬結果與實測數(shù)據(jù)對比：聚集區(qū)（見圖4c）含油飽和度主要分布在40%～70%，加權平均值約在60%左右；70個油藏的實測值分布在60%～63%，平均值為61.17%，兩者比較接近。②聚集量分析。按含油飽和度大于10%的聚集量計算，模擬聚集量為41 006×104t。目前，目標層探明量為34 613.53×104t（見表1），探明率為84.4%，說明該目標層的勘探程度較高，比較符合目前勘探進程。③聚集區(qū)位置分析。高資源豐度（大于300×104t/km2）區(qū)與資源規(guī)模分布基本一致，主要分布在B區(qū)和C區(qū)（見圖5），其中，B區(qū)聚集18 417×104t，占總量的44.9%，C區(qū)聚集19 724×104t，占總量的48.1%，A區(qū)聚集2 865×104t，僅占總量的7%，D沒有聚集量。4個分區(qū)石油聚集位置與勘探結果較為一致（見圖3），各區(qū)聚集量比例與實際勘探數(shù)據(jù)（見表1）雖有所差異，但總體偏差不大。對比目標層已發(fā)現(xiàn)的油田分布，模擬結果預測出在C區(qū)NP2-16井附近、A區(qū)NP5-4井北側還有勘探潛力，是下步應該關注的目標。

圖6 南堡凹陷石油運移主要路徑追蹤結果

5 結論

建立的順層柱狀PEBI網(wǎng)格動態(tài)三維地質模型可以更好地刻畫地層的變化，采用隨機抽樣方法可以較好地解決不同沉積相內儲集層非均質性問題，引入有效儲集層比例參數(shù)，設置有效儲集層存在的下限值等方法可以降低混合巖性對模擬精度的影響程度。

基于順層柱狀PEBI網(wǎng)格動態(tài)三維地質模型建立的滲流方程，是一種改進的黑油模型方程，與原黑油模型方程的最大區(qū)別在于方程中的網(wǎng)格體積是與時間有關的變量；變網(wǎng)格模擬技術最大的優(yōu)點是能有效提高地層與流體變化的模擬精度。

采用全張量滲透率不僅可以滿足變網(wǎng)格滲流方程的要求，也能夠更好地刻畫類似河道輸導層、斷裂與裂縫輸導通道所引起的地層滲透率各向異性，提高復雜地區(qū)的模擬精度。

通過對關鍵參數(shù)曲線的光滑處理、網(wǎng)格流動上下游的穩(wěn)定性處理、小體積網(wǎng)格與零體積網(wǎng)格的特殊處理，可以有效提高牛頓法的收斂性能；采用自動調整時間步長方法，可以在確保收斂穩(wěn)定性的同時提高運算效率。

模擬實例與勘探現(xiàn)狀的對比表明，模擬的含油飽和度、油氣聚集量和分布位置等方面結果符合勘探現(xiàn)狀；實例揭示了目標層資源探明率為84.4%，主要資源分布在B區(qū)和C區(qū)，待發(fā)現(xiàn)資源主要分布在C區(qū)NP2-16井附近和A區(qū)NP5-4井北側的構造上。

三維油氣運聚模擬技術在實際應用中還存在地質

模型過于簡化、地質參數(shù)不能滿足數(shù)值模型要求等多方面的問題，本文在較精細化的地質模型與數(shù)值模型緊密結合方面做了嘗試，效果較好。

符號注釋：

Bl——l相流體的體積系數(shù)，m3/m3；cK,m——m方向的滲透率壓縮系數(shù)，m=x,y,z，Pa-1；——網(wǎng)格i，j相對于基準面的埋深，m；g——重力加速度，取9.8 m/s2；g——重力加速度向量，m/s2；i，j——網(wǎng)格編號；l——流體相，,分別代表水、油、氣相；Kx，Ky，Kz——x，y，z方向的絕對滲透率，m2；K——絕對滲透率張量，m2；Krl——l相流體的相對滲透率，m2；K0,m——參考壓力p0下m方向的絕對滲透率，m2；Lij——網(wǎng)格i中心到網(wǎng)格i，j界面中心的向量，m；——向量Lij的模，m；Lji——網(wǎng)格j中心到網(wǎng)格i， j界面中心的向量，m；——向量Lji的模，m；nij——向量Lij方向的單位向量；nji——向量Lji方向的單位向量；p0——參考壓力，Pa；pcow——油/水系統(tǒng)毛管壓力，Pa；——油/氣系統(tǒng)毛管壓力，Pa；pl——l相流體的壓力，Pa；?pl——l相流體的壓力梯度，Pa/m；ql——控制體中的l相源匯量，m3/s；Rs——溶解氣油比，m3/m3；Sij——網(wǎng)格i，j間界面的面積向量，m2；Sl——l相流體飽和度，%； t——時間，s；——時間步長，s；ul——通過任一面元的滲流速度向量，m/s；V——網(wǎng)格體積，m3；——l相流體密度，kg/m3；——l相流體標準狀態(tài)下的密度，kg/m3；f——孔隙度，%；lm——l相流體的黏度，Pa·s；——網(wǎng)格i，j間l相流體的勢差，Pa；——網(wǎng)格i，j間l相流體的門限勢差，Pa。下標：sc——標準狀態(tài)。

[1] Hantschel T,Kauerauf A I.Fundamentals of basin modeling and petroleum systems modeling[M].Berlin: Springer-Verlag,2009.

[2] 郭秋麟,楊文靜,肖中堯,等.不整合面下縫洞巖體油氣運聚模型[J].石油實驗地質,2013,35(5): 495-499.Guo Qiulin,Yang Wenjing,Xiao Zhongyao,et al.Hydrocarbon migration and accumulation model of fractured-vuggy reservoir under unconformity surface[J].Petroleum Geology &Experiment,2013,35(5): 495-499.

[3] 石廣仁,張慶春,馬進山,等.三維三相烴類二次運移模型[J].石油學報,2003,24(2): 38-42.Shi Guangren,Zhang Qingchun,Ma Jinshan,et al.3-D and 3-phase model for secondary migration of hydrocarbon[J].Acta Petrolei Sinica,2003,24(2): 38-42.

[4] Mello U T,Rodrigues J R P,Rossa A L.A control-volume finite element method for three-dimensional multiphase basin modeling[J].Marine ＆ Petroleum Geology,2009,26(4): 504-518.

[5] 石廣仁,馬進山,常軍華.三維三相達西流法及其在庫車坳陷的應用[J].石油與天然氣地質,2010,31(4): 403-409.Shi Guangren,Ma Jinshan,Chang Junhua.3-D 3-phasce Darcy flow method and its application to the Kuqa Depression[J].Oil &Gas Geology,2010,31(4): 403-409.

[6] Shi Guangren,Ma Jinshan,Yang Xinshe,et al.Finite volume method for solving a modified 3-d 3-phase black-oil hydrocarbon secondary migration model,and its application to the Kuqa Depression of the Tarim Basin in Western China[J].Advances in Petroleum Exploration and Development,2011,2(1): 1-12.

[7] 石廣仁.油氣運聚定量模擬技術現(xiàn)狀、問題及設想[J].石油與天然氣地質,2009,30(1): 1-10.Shi Guangren.Status,problems and proposals of the quantitative modeling techniques for hydrocarbon migration and accumulation[J].Oil &Gas Geology,2009,30(1): 1-10.

[8] 馮勇,石廣仁,米石云,等.有限體積法及其在盆地模擬中的應用[J].西南石油學院學報,2001,23(5): 12-15.Feng Yong,Shi Guangren,Mi Shiyun,et al.Finite volume method and its application in basin modeling[J].Journal of Southwest Petroleum Institute,2001,23(5): 12-15.

[9] 范天一,宋新民,吳淑紅,等.低滲透油藏水驅動態(tài)裂縫數(shù)學模型及數(shù)值模擬[J].石油勘探與開發(fā),2015,42(4): 496-501.Fan Tianyi,Song Xinmin,Wu Shuhong,et al.A mathematical model and numerical simulation of waterflood induced dynamic fractures of low permeability reservoirs[J].Petroleum Exploration and Development,2015,42(4): 496-501.

[10] 高麗明,何登發(fā),桂寶玲,等.東營凹陷民豐洼陷邊界斷層三維幾何學及運動學特征[J].石油勘探與開發(fā),2014,41(5): 546-553.Gao Liming,He Dengfa,Gui Baoling,et al.3D geometrical and kinematic characteristics of boundary faults in Minfeng subsag,Dongying Sag,Bohai Bay Basin[J].Petroleum Exploration and Development,2014,41(5): 546-553.

[11] 姜福杰,董月霞,龐雄奇,等.南堡凹陷油氣分布特征與主控因素分析[J].現(xiàn)代地質,2013,27(5): 1258-1264.Jiang Fujie,Dong Yuexia,Pang Xiongqi,et al.Hydrocarbon distribution features and main controlling factors in the Nanpu Sag[J].Geoscience,2013,27(5): 1258-1264.

[12] 龐雄奇,霍志鵬,范泊江,等.渤海灣盆地南堡凹陷源控油氣作用及成藏體系評價[J].天然氣工業(yè),2014,34(1): 28-36.Pang Xiongqi,Huo Zhipeng,Fan Bojiang,et al.Control of source rock on hydrocarbon accumulation and assessment of gas pools in the Nanpu Sag,Bohai Bay Bain[J].Natural Gas Industry,2014,34(1): 28-36.

[13] 蒲秀剛,周立宏,韓文中,等.歧口凹陷沙一下亞段斜坡區(qū)重力流沉積與致密油勘探[J].石油勘探與開發(fā),2014,41(2): 138-149.Pu Xiugang,Zhou Lihong,Han Wenzhong,et al.Gravity flow sedimentation and tight oil exploration in lower first member of Shahejie Formation in slope area of Qikou Sag,Bohai Bay Basin[J].Petroleum Exploration and Development,2014,41(2): 138-149.

[14] 董月霞,楊賞,陳蕾,等.渤海灣盆地辮狀河三角洲沉積與深部儲集層特征: 以南堡凹陷南部古近系沙一段為例[J].石油勘探與開發(fā),2014,41(4): 385-393.Dong Yuexia,Yang Shang,Chen Lei,et al.Braided river delta deposition and deep reservoir characteristics in Bohai Bay Basin: A case study of Paleogene Sha 1 Member in the south area of Nanpu Sag[J].Petroleum Exploration and Development,2014,41(4): 385-393.

[15] 張林曄,包友書,李鉅源,等.湖相頁巖油可動性: 以渤海灣盆地濟陽坳陷東營凹陷為例[J].石油勘探與開發(fā),2014,41(6): 641-649.Zhang Linye,Bao Youshu,Li Juyuan,et al.Movability of lacustrine shale oil: A case study of Dongying Sag,Jiyang Depression,Bohai Bay Basin[J].Petroleum Exploration and Development,2014,41(6): 641-649.

（編輯 黃昌武）

Three-dimensional hydrocarbon migration and accumulation modeling based on finite volume method

Guo Qiulin,Chen Ningsheng,Xie Hongbing,Wu Xiaozhi,Liu Jifeng,Zhao Xiran,Gao Rili,Hu Junwen
(PetroChina Research Institute of Petroleum Exploration &Development,Beijing 100083,China)

Aiming at the issue that the application effect of current 3-D hydrocarbon migration and accumulation modeling technique cannot meet the actual requirement because of over-simplified geologic model,a 3-D dynamic geologic model using PEBI grid in bedding and columnar shape was established,which can be used to solve the heterogeneity in reservoirs of various sedimentary facies,and to reduce the effect of mixed lithologies on modeling accuracy.A 3-D numerical modeling technique based on finite volume method was established,including the construction of variable grid flow equation,decomposition and calculation of full-tensor permeability,and calculation of conductivity etc.Some related processing techniques,such as smoothing of key parameter curves,stability processing of upstream and downstream of grid flow,etc.were used to improve the convergence of Newton’s iteration.Automatic adjustment of time step and multi-core parallel computation were taken to enhance the running efficiency of the software.In a case study of Nanpu sag in the Bohai Bay Basin,the hydrocarbon migration and accumulation process,oil saturation and hydrocarbon resources abundance in various periods were simulated,the modeled oil saturation,accumulation quantity and hydrocarbon distribution locations coincide with current exploration results.The results also showed that the proved ratio of hydrocarbon resources in target layers is 84.4%,the majority of the resources are distributed in Block B and Block C,and it is predicted that undiscovered resources are mainly near Well NP2-16 in Block C and in the structure north of Well NP5-4 in Block A.

3-D hydrocarbon migration and accumulation modeling technique;finite volume method;PEBI grid in bedding and columnar shape;variable grid flow equation;Nanpu sag

國家重大油氣專項“巖性地層區(qū)帶、圈閉評價與儲層預測技術”（2011ZX05001-003）；中國石油重大科技專項“中國石油第四次油氣資源評價”（2013E-0502）；中國石油天然氣集團公司自主知識產權軟件推廣應用（2014D-1809）

TE122.2

A

1000-0747(2015)06-0817-09

10.11698/PED.2015.06.17

郭秋麟（1963-），男，福建漳州人，博士，中國石油勘探開發(fā)研究院教授級高級工程師，主要從事盆地模擬與油氣資源評價研究。地址：北京市海淀區(qū)學院路20號，中國石油勘探開發(fā)研究院油氣資源規(guī)劃所，郵政編碼：100083。E-mail：qlguo@petrochina.com.cn

2014-12-31

2015-10-10