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

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

基于有限體積法的三維油氣運(yùn)聚模擬技術(shù)

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

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

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

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

0 引言

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

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

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

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

1 網(wǎng)格劃分與關(guān)鍵地質(zhì)問題處理

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

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

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

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

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

1.2 關(guān)鍵地質(zhì)問題的處理

1.2.1 地層物性非均質(zhì)性處理方法

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

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

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

1.2.3 混合巖性的處理方法

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

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

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

假設(shè)網(wǎng)格體運(yùn)動(dòng)速度與流體流動(dòng)速度相比可以忽略不計(jì)，建立變網(wǎng)格滲流方程。

2.1.1 質(zhì)量守恒方程

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

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

（1）式對(duì)時(shí)間積分得：

2.1.2 滲流運(yùn)動(dòng)方程

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

2.1.3 流動(dòng)方程

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

油相流動(dòng)方程為：

氣相流動(dòng)方程為：

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

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

輔助方程為：

2.2 初始條件和邊界條件

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

2.3 全張量滲透率的分解與計(jì)算

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

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

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

2.4 傳導(dǎo)率的計(jì)算

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

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

要求Sij與nij、Sji與nji位于界面的同側(cè)且點(diǎn)積大于0。

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

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

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

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

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

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

3.2 自動(dòng)調(diào)整時(shí)間步長(zhǎng)以提高計(jì)算速度

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

圖1 自動(dòng)調(diào)整時(shí)間步長(zhǎng)實(shí)例

3.3 多核并行計(jì)算以提高計(jì)算速度

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

圖2 多核并行計(jì)算示意圖

4 實(shí)例分析

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

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

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

4.2 目標(biāo)層段特征

選擇Es—Ed作為模擬目標(biāo)層。其中，烴源巖層為沙河街組和東三段，儲(chǔ)集層為東營(yíng)組一段和二段，蓋層為館陶組的火山巖。

4.2.1 探明儲(chǔ)量分布

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

4.2.2 生烴潛力

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

圖3 南堡凹陷構(gòu)造簡(jiǎn)圖

表1 東營(yíng)組一段石油探明儲(chǔ)量在各構(gòu)造的分布比例

4.2.3 沉積特征與儲(chǔ)集層物性

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

4.2.4 蓋層分布

直接覆蓋在東營(yíng)組之上的館陶組火山巖是主要蓋

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

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

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

4.2.5 斷層作用

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

表4 斷層帶滲透率特征

4.3 三維地質(zhì)體構(gòu)建及主要參數(shù)

4.3.1 三維地質(zhì)體

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

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

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

表5 模擬層段地層特征

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

4.4 模擬結(jié)果

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

4.4.1 含油飽和度與石油聚集量模擬結(jié)果

模擬時(shí)間從距今25 Ma開始至今，按含油飽和度大于10%統(tǒng)計(jì)，不同地質(zhì)時(shí)期石油聚集量和聚集區(qū)面積見表7。由表7可見，在距今20 Ma時(shí)，石油開始聚集；在距今15 Ma時(shí)，已聚集石油7 074×104t，占總聚集量的17.3%，此時(shí)石油主要聚集在北部的M10井和GC1井附近及中部的NP2-16井附近（見圖4a）；在距今8 Ma時(shí)，已聚集石油28 626×104t，占總聚集量的69.8%，此時(shí)在西南部的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 不同地質(zhì)時(shí)期石油聚集量和聚集區(qū)面積

4.4.2 油氣運(yùn)移主要路徑追蹤

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

圖4 東一段含油飽和度模擬結(jié)果

圖5 東一段石油資源豐度模擬結(jié)果

4.4.3 模擬結(jié)果分析

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

方面進(jìn)行分析。①含油飽和度模擬結(jié)果與實(shí)測(cè)數(shù)據(jù)對(duì)比：聚集區(qū)（見圖4c）含油飽和度主要分布在40%～70%，加權(quán)平均值約在60%左右；70個(gè)油藏的實(shí)測(cè)值分布在60%～63%，平均值為61.17%，兩者比較接近。②聚集量分析。按含油飽和度大于10%的聚集量計(jì)算，模擬聚集量為41 006×104t。目前，目標(biāo)層探明量為34 613.53×104t（見表1），探明率為84.4%，說明該目標(biāo)層的勘探程度較高，比較符合目前勘探進(jìn)程。③聚集區(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個(gè)分區(qū)石油聚集位置與勘探結(jié)果較為一致（見圖3），各區(qū)聚集量比例與實(shí)際勘探數(shù)據(jù)（見表1）雖有所差異，但總體偏差不大。對(duì)比目標(biāo)層已發(fā)現(xiàn)的油田分布，模擬結(jié)果預(yù)測(cè)出在C區(qū)NP2-16井附近、A區(qū)NP5-4井北側(cè)還有勘探潛力，是下步應(yīng)該關(guān)注的目標(biāo)。

圖6 南堡凹陷石油運(yùn)移主要路徑追蹤結(jié)果

5 結(jié)論

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

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

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

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

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

三維油氣運(yùn)聚模擬技術(shù)在實(shí)際應(yīng)用中還存在地質(zhì)

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

符號(hào)注釋：

Bl——l相流體的體積系數(shù)，m3/m3；cK,m——m方向的滲透率壓縮系數(shù)，m=x,y,z，Pa-1；——網(wǎng)格i，j相對(duì)于基準(zhǔn)面的埋深，m；g——重力加速度，取9.8 m/s2；g——重力加速度向量，m/s2；i，j——網(wǎng)格編號(hào)；l——流體相，,分別代表水、油、氣相；Kx，Ky，Kz——x，y，z方向的絕對(duì)滲透率，m2；K——絕對(duì)滲透率張量，m2；Krl——l相流體的相對(duì)滲透率，m2；K0,m——參考?jí)毫0下m方向的絕對(duì)滲透率，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——參考?jí)毫Γ琍a；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——時(shí)間，s；——時(shí)間步長(zhǎng)，s；ul——通過任一面元的滲流速度向量，m/s；V——網(wǎng)格體積，m3；——l相流體密度，kg/m3；——l相流體標(biāo)準(zhǔn)狀態(tài)下的密度，kg/m3；f——孔隙度，%；lm——l相流體的黏度，Pa·s；——網(wǎng)格i，j間l相流體的勢(shì)差，Pa；——網(wǎng)格i，j間l相流體的門限勢(shì)差，Pa。下標(biāo)：sc——標(biāo)準(zhǔn)狀態(tài)。

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

[2] 郭秋麟,楊文靜,肖中堯,等.不整合面下縫洞巖體油氣運(yùn)聚模型[J].石油實(shí)驗(yàn)地質(zhì),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ìn)山,等.三維三相烴類二次運(yùn)移模型[J].石油學(xué)報(bào),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ìn)山,常軍華.三維三相達(dá)西流法及其在庫車坳陷的應(yīng)用[J].石油與天然氣地質(zhì),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] 石廣仁.油氣運(yùn)聚定量模擬技術(shù)現(xiàn)狀、問題及設(shè)想[J].石油與天然氣地質(zhì),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] 馮勇,石廣仁,米石云,等.有限體積法及其在盆地模擬中的應(yīng)用[J].西南石油學(xué)院學(xué)報(bào),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] 范天一,宋新民,吳淑紅,等.低滲透油藏水驅(qū)動(dòng)態(tài)裂縫數(shù)學(xué)模型及數(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ā),桂寶玲,等.東營(yíng)凹陷民豐洼陷邊界斷層三維幾何學(xué)及運(yùn)動(dòng)學(xué)特征[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)代地質(zhì),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] 龐雄奇,霍志鵬,范泊江,等.渤海灣盆地南堡凹陷源控油氣作用及成藏體系評(píng)價(jià)[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ī)?chǔ)集層特征: 以南堡凹陷南部古近系沙一段為例[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] 張林曄,包友書,李鉅源,等.湖相頁巖油可動(dòng)性: 以渤海灣盆地濟(jì)陽坳陷東營(yíng)凹陷為例[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

國(guó)家重大油氣專項(xiàng)“巖性地層區(qū)帶、圈閉評(píng)價(jià)與儲(chǔ)層預(yù)測(cè)技術(shù)”（2011ZX05001-003）；中國(guó)石油重大科技專項(xiàng)“中國(guó)石油第四次油氣資源評(píng)價(jià)”（2013E-0502）；中國(guó)石油天然氣集團(tuán)公司自主知識(shí)產(chǎn)權(quán)軟件推廣應(yīng)用（2014D-1809）

TE122.2

A

1000-0747(2015)06-0817-09

10.11698/PED.2015.06.17

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

2014-12-31

2015-10-10