三維三相達西流法及其在庫車坳陷的應用

石廣仁，馬進山，常軍華

（中國石油天然氣股份有限公司 石油勘探開發(fā)研究院，北京 100083）

三維三相達西流法及其在庫車坳陷的應用

石廣仁，馬進山，常軍華

（中國石油天然氣股份有限公司 石油勘探開發(fā)研究院，北京 100083）

多相達西流法是油氣二次運移模擬技術中最復雜的方法，存在三大技術難點：①模型十分復雜，求解比油藏模擬的難度大得多；②收斂性、穩(wěn)定性不易保證，導致計算結果精度低、計算中途異常停止；③計算機耗時量巨大，在微機上算完一個大型盆地要花幾天時間。為此，采用了基于PEBI（直交平分）網的有限體積法、基于牛頓迭代的全隱式法及ORTHMIN（正交極小化）法，摸索出了“正確地選取一個模擬工區(qū)的最大允許時間步長”是解決上述三大技術難點的最大關鍵點，從而攻克了這三大技術難點。研發(fā)了三維三相達西流軟件，并應用于庫車坳陷，獲得了較好的油氣運聚史模擬結果。

有限體積法；PEBI網；最大允許時間步長；多相達西流法；盆地模擬；二次運移；庫車坳陷

烴類運移聚集史模擬是盆地模擬系統(tǒng)中最重要的部分，也是迄今為止技術上最薄弱的環(huán)節(jié)。當今世界上有4種烴類運聚定量模擬技術［1，2］：多相達西流法、流徑法、混合法及侵入逾滲法。其中，多相達西流法是技術最復雜也是采用最早的方法，曾被譽為“主力技術”，但其存在三大技術難點：①模型十分復雜，求解比油藏模擬的難度大得多；②收斂性、穩(wěn)定性不易保證，導致計算結果精度低、計算中途異常停止；③計算機耗時量巨大，在微機上算完一個大型盆地要花幾天時間。為此，筆者攻克了這三大技術難點，研發(fā)了三維三相達西流軟件，并將此應用于庫車坳陷，獲得了較好的油氣運聚史模擬結果。

1 三維三相達西流軟件

目前，世界上有5個多相達西流軟件（表1），其中PetroFlow 3D和Temis 3D應用較多。多相達西流軟件的模擬器多數(shù)為改造的三維三相黑油模型，只是采用的網格及數(shù)值解法不同而已。

表1 本三維三相達西流軟件與世界上同類軟件的主要特性比較Table 1 Com parison ofmajor characteristics of this 3-D 3-phase Darcy flow software w ith other sim ilar softwares in the world

表1中提到了3種類型的網格：①矩形網，是指固定矩形面組成的平行六面體；②三角網，是指可變三角面組成的四面體；③PEBI網，是指具有PEBI（perpendicular and bisection，譯為直交平分或垂直平分）性質的可變多角面組成的棱柱體。三維地質體的剖分按棱柱體進行，能確保每個網格塊的中心點與任一相鄰網格塊的中心點之間的連線垂直于這兩個網格塊的交面且被交面平分。

黑油模型在盆地模擬與油藏模擬中的模擬條件是不同的（表2），它在盆地模擬中的求解遠比它在油藏模擬中的求解困難，尤其是在求解的穩(wěn)定性和收斂性方面，這主要是由盆地模擬可變的模擬范圍所引起的。

1.1 三維地質體的網格化

在盆地演化期間，盆地的幾何體積因沉降及沉積物填充而大量地全局性地增大，因沉積壓實而小量地全局性地變小，因剝蝕而大量或小量地局部性地變小，因斷層及裂縫而大量或小量地局部性地增大，甚至因沉積間斷而基本不變。因此，地質家所研究的對象實為一個隨時間變化的地質體。為了定量地模擬一個盆地，其隨時間變化的地質體必須剖分為相當數(shù)量的小塊（即稱網格塊）。本軟件選用了PEBI網。

1.2 基于PEBI網的有限體積法

三維三相黑油模型的流動方程組［6］太復雜，以致無法對一個實際的地質模型給出解析解，應對之進行離散化而變成線性代數(shù)方程組。離散化的方法有有限差分法、有限元法和有限體積法。本軟件采用基于PEBI網的有限體積法，從而獲得：

1）水流動方程

表2 盆地模擬與油藏模擬的不同模擬條件Table 2 Different conditions between basin modeling and reservoir simulation

2）油流動方程

3）氣流動方程

上三式中：v——任一PEBI網格塊的符號；

Vi——任一PEBI網格塊（由x，y，z表達）；

?——關于v的三重積分；

div——向量場的散度；

al——等于［K］Krl/（μlBl）（l=w，o，g），因為［K］是絕對滲透率張量（即有方向性的絕對滲透率），故al也是張量；

Krw，Kro，Krg——水、油、氣的相對滲透率，可由油層物理試驗及改進的Stone公式確定；

μw，μo，μg——水、油、氣的粘度，可由PVT函數(shù)確定；

Bw，Bo，Bg——水、油、氣的地層體積因子，可由PVT函數(shù)確定；

grad——兩個相鄰PEBI網格塊之間的梯度；

Φw，Φo，Φg——水、油、氣的勢，分別是水壓力Pw、油壓力Po、氣壓力Pg的函數(shù)；

Qw，Qo，Qg——水、油、氣的源或匯，若源取正值，若匯取負值；

Rs——溶解氣油比，可由PVT函數(shù)確定；

t——地質時間；

φ——地層孔隙度；

Sw，So，Sg——水、油、氣的飽和度。

由于PEBI網的優(yōu)點（相鄰兩網格塊中心點的連線恰好就是這兩網格塊的公共面的法線方向），使方程組（1），（2），（3）的離散化變得比較容易，最終化成為水、油、氣三相離散方程組。然而，這個離散方程組只有3個方程，而未知數(shù)卻有6個（水、油、氣壓力Pw，Po，Pg和水、油、氣飽和度Sw，So，Sg）。所以，還需要其他3個獨立的方程。這3個獨立的方程是：飽和度平衡方程、油水系統(tǒng)中的毛細管力方程及油氣系統(tǒng)中的毛細管力方程。它們分別表示如下：

上三式中：Pcow——油水系統(tǒng)中的毛細管力；

Pcog——油氣系統(tǒng)中的毛細管力。

Pcow和Pcog可由油層物理試驗確定。

當給出了一個三維地質體的邊界條件、初始條件［6］時，就可以通過式（1）—（6）這6個方程，求解出該三維地質體的Pw，Po，Pg及Sw，So，Sg6個未知數(shù)。這6個未知數(shù)表示了該三維地質體油氣二次運移的基本狀況。

1.3 基于牛頓迭代的全隱式法

由于式（1）—（3）及其離散方程組都是非線性方程組，因此，①采用全隱式法，而不采用半隱式法，因為半隱式法可能因采用過時的飽和度而得到不合理的計算結果；②采用牛頓迭代法，而不采用其他迭代法，因為牛頓迭代法既有超線性收斂速度，又能保持矩陣的稀疏性。

牛頓算法分兩步進行：①生成Jacobi矩陣；②求解問題變?yōu)榇笮途€性方程組的Jacobi方程。

Jacobi方程是一個大型、稀疏、非對稱的線性代數(shù)方程組，即：

式中：A——二維（3I，3I）的系數(shù)矩陣，即Jacobi矩陣，其中I為網格塊個數(shù)；

X——一維（3I）的未知數(shù)矩陣，含有3I個未知數(shù)（δSw1，δSg1，δPo1，δSw2，δSg2，δPo2，…，δSwI，δSgI，δPoI）；

B——一維（3I）的常數(shù)矩陣，含有3I個余量。

1.4 大型線性方程組的求解

求解Jacobi方程，選取存儲空間盡可能小的存儲方式以及速度盡可能快的迭代算法是十分重要的。對于這類線性代數(shù)方程組，迄今尚缺少一致公認最優(yōu)的迭代算法。正交極小化（ORTHMIN，即orthogonalminimum）法經油藏模擬的實算表明是行之有效的方法，故選擇了正交極小化法。正交極小化法的具體步驟為：首先對矩陣A進行近似的不完全因子分解，這比其他方法更接近真解；再經過正交極小化法對AX=B進行求解。正交極小化法的優(yōu)點是：快速收斂，對迭代參數(shù)不敏感，對矩陣帶寬不敏感，對矩陣是否對稱不敏感。

1.5 解決三大技術難點的關鍵點

如何解決前述的三大技術難點？除了采用上述的基于PEBI網的有限體積法、基于牛頓迭代的全隱式法及正交極小化法外，還摸索出“正確地選取一個模擬工區(qū)的最大允許時間步長”是解決上述三大技術難點的最大關鍵點。

在理論上，全隱式法具有無條件的收斂性和穩(wěn)定性。然而，實際應用時發(fā)現(xiàn)并非如此。為了確保收斂性和穩(wěn)定性，合適地選擇正交極小化法的時間步長是至關重要的。一般來說，時間步長越大，計算機耗時量就越小，但收斂性和穩(wěn)定性就越差；時間步長越小，收斂性和穩(wěn)定性就越好，但計算機耗時量就越大。

在下面的庫車坳陷實例中，對時間步長的選取進行了試驗，結果發(fā)現(xiàn)：①若選取的時間步長太大，會影響模擬的穩(wěn)定性（即正交極小化發(fā)生浮點溢出而導致運行失?。?；②若選取的時間步長不足夠小，雖然沒有出現(xiàn)穩(wěn)定性問題，但再采用比它小一點的時間步長，會出現(xiàn)兩者的模擬結果存在不允許的差異（即收斂性問題）；③若選取的時間步長太小，會造成過長的計算時間而浪費機時。因此，應選取在確保穩(wěn)定性和收斂性條件下的最大時間步長，這個步長就稱之謂最大允許時間步長（0.025 Ma）。這里說明一下如何獲得0.025 Ma。令時間步長為Δt。試算表明：當Δt≥0.5 Ma時（例如Δt=5，2.5，1，0.5 Ma），計算發(fā)生浮點溢出而導致運行失??；當0.5 Ma＞Δt≥0.025 Ma時（例如Δt=0.25，0.1，0.05，0.025 Ma），雖然計算沒有發(fā)生浮點溢出而導致運行失敗，但使用這4個時間步長算出的計算結果的差異超過了允許的誤差；當Δt=0.01 Ma時，計算結果與Δt= 0.025 Ma時的差異沒有超過允許的誤差；當Δt= 0.005 Ma時，計算結果與Δt=0.025，0.01 Ma時的差異都沒有超過允許的誤差。這表示Δt=0.025，0.01，0.005Ma都是允許的。但Δt=0.025Ma時的計算機耗時量為9 h 38 min；而Δt=0.01，0.005 Ma時的計算機耗時量分別為19 h 39 min，37 h 8 min。顯然，Δt應選取0.025 Ma，它正是這3個Δt中的最大者。這是在平面網格塊最大個數(shù)取為900時的試算結果。實際上，由于最大允許時間步長基本上不依賴于平面網格塊最大個數(shù)，所以這種試算可在平面網格塊最大個數(shù)較小的條件下進行。例如，平面網格塊最大個數(shù)取為100，也可以獲得最大允許時間步長為0.025 Ma的結論，這樣可大大節(jié)省試算的機時（計算機耗時量為59 min）。

2 在庫車坳陷的應用

2.1 地質背景

庫車坳陷位于塔里木盆地的北部，北倚南天山造山帶，南界塔北隆起，東西長約380 km，南北寬約30～145 km，西寬東窄，面積約30 000 km2。截至2003底，該坳陷內已鉆探井102口，其中獲工業(yè)油流井2口，獲工業(yè)氣流井56口。油氣藏類型有凝析油氣藏、干氣氣藏和油藏3種。在該坳陷中已經發(fā)現(xiàn)2個油田及13個氣田。庫車坳陷已成為中國最富的天然氣聚集區(qū)之一（圖1；表3）。

圖1 庫車坳陷構造及油氣田分布［8］Fig.1 Structures and oil/gas fields in the Kuqa Depression［8］

表3 庫車坳陷地層年齡Table 3 Geological ages of formations in the Kuqa Depression

庫車坳陷是在海西運動基礎上于晚二疊世開始發(fā)育起來的中、新生代疊合坳陷。坳陷演化具有多階段性和復雜性，中、新生代陸相地層發(fā)育齊全，后期褶皺-沖斷變形改造強烈。主力烴源巖為中、上三疊統(tǒng)湖相泥巖和中、下侏羅統(tǒng)湖泊沼澤相煤系地層，平均厚度達600 m。暗色泥巖有機碳含量多在1.5％～5.0％，煤層的殘余有機碳均值為58％。源巖干酪根類型以腐殖型（Ⅲ型）為主，約占90％；其次是腐泥-腐殖型（Ⅱ）型，約為10％。源巖現(xiàn)今主體成熟度Ro在0.7％～1.8％，處于成熟-高成熟階段，以生氣為主。第三系巨厚（100～1 000 m）的鹽-膏-泥優(yōu)質區(qū)域蓋層封蓋了整個坳陷，而這個區(qū)域蓋層之下分布著古近系和白堊系的良好砂巖儲集層。優(yōu)越的源、蓋條件形成了克拉2、大北等大型氣藏以及牙哈、英買7、羊塔克等中型凝析油氣田［7～9］。

2.2 模擬條件

石廣仁和張慶春［10］曾使用簡易的方法對庫車坳陷做過油氣運聚史模擬，介紹并使用了兩類輸入參數(shù)：①坳陷內71口已鉆井的數(shù)據(jù)以及從地震剖面上所取的577口人工模擬井的數(shù)據(jù)，共計90 000余個；②地史、熱史、成巖史、生烴史及排烴史的模擬結果。該簡易方法視地層為一個具有厚度的平面，故為擬三維，而不是真三維；只算油、氣二相，不算油、氣、水三相；使用簡易的流體流動方法，故不屬于多相達西流法。本文討論的才是三維三相達西流法，但也同樣使用了上述的兩類輸入參數(shù)。

運移聚集史模型是正演模型，故歷史模擬是從古到今進行的。一般來說（假設盆地的最大埋深是在現(xiàn)今），三維地質體的網格塊個數(shù)隨之由小到大變化著，則現(xiàn)今網格塊個數(shù)為最大，這一個數(shù)是很大的。被模擬的庫車坳陷的面積約為30 000 km2，東西向最長為383 km，南北向最長為145 km，最大歷史埋深為10.4 km。在三維地質體的剖分中，網格塊垂向長度取為100 m，平面網格塊最大個數(shù)取為900，則實際的最大網格塊個數(shù)I=31 660（節(jié)點）。因此，式（7）中的二維系數(shù)矩陣A（3I，3I）最大時為A（94 980，94 980）。

計算機耗時量與最大網格塊個數(shù)I有關，還與地層個數(shù)、烴源層個數(shù)，尤其是模擬歷史、模擬時間步長有關。被模擬的庫車坳陷有8個地層（表3）、2個烴源層（T，J），模擬歷史為208～0 Ma，模擬時間步長取最大允許時間步長（0.025 Ma），計算機耗時量為9 h 38 min（這是用戶可以接受的）。該實例表明，目前微機可承受30 000節(jié)點的三維三相達西流模擬。

2.3 模擬結果

使用模擬結果，可以用三維立體圖及其他形式的圖件來顯示各地層的油、氣聚集量史。但是這里由于篇幅所限，①為了簡明起見，改用曲線形式表示各地層的油、氣聚集量史（圖2）；②因為白堊系是最大的油氣儲集層（圖2），用三維顯示其現(xiàn)今（0 Ma）的油聚集量（圖3a）以及氣聚集量（圖3b）；③為了分析古近系和白堊系的儲氣層，給出了平面等值線圖（圖4）。

圖2 庫車坳陷各相關地層油、氣聚集量史Fig.2 Oil and gas accumulation histories of each related formation in the Kuqa Depression

2.4 地質分析

2.4.1 坳陷的模擬總量與勘探成果吻合

表4列出了庫車坳陷生、排烴的模擬總量，其中排油系數(shù)=排油量/生油量，排氣系數(shù)=排氣量/生氣量［10］。在排烴史模擬結果的基礎上，這次采用三維三相達西流軟件繼續(xù)進行運聚史的計算，獲得了聚烴的模擬總量，其中聚油系數(shù)=聚油量/排油量，聚氣系數(shù)=聚氣量/排氣量（圖2；表4）。在以前庫車坳陷的圈閉評價中，曾獲得了該坳陷的油氣地質儲量（含預測、控制、探明儲量）（表4）［10］。從表4可見，模擬的聚油量（409.5× 106t）稍大于油的地質儲量（380×106t），相對誤差為7.8％；模擬的聚氣量（2 277×109m3）也稍大于氣的地質儲量（2 100×109m3），相對誤差為8.4％。這表明了模擬結果的合理性，也表明庫車坳陷還有深入勘探的前景。此外，模擬的聚氣量為模擬的聚油量的5.56倍，氣的地質儲量為油的地質儲量的5.53倍，兩者十分接近，均說明庫車坳陷是富氣的。

圖3 白堊系在0 Ma時的累計油、氣聚集量Fig.3 Cumulative oil and gas accumulations in the Cretaceous（K）at0 Ma

圖4 在0 Ma時a）古近系累計油聚集量和b）白堊系累計氣聚集量等值線Fig.4 Cumulative oil accumulation in the Paleogene（a）and cumulative gas accumulation in the Cretaceous（b）at 0 Ma

表4 庫車坳陷烴類生、排、聚總量（模擬）及地質儲量Table 4 Simulated total amount of hydrocarbons generated，expelled and accumulated in the Kuqa Depression and its reserves in p lace

由圖2可見，坳陷的油、氣聚集時期主要是新近紀（24～2 Ma）。烴類排、聚時期與主要構造的形成時期基本是吻合的。這個模擬結果與王庭斌教授通過多方面的分析所得的結論（即“庫車坳陷的氣藏主要形成、定型于新近紀”）［11］基本符合。

2.4.2 各相關地層的模擬總量與勘探成果吻合

從圖2a可見，新近系（N1j，N1-2k，N2k）無油；古近系（E）只有極少量的油；而白堊系（K）才是最大的儲油層（圖3a），約占坳陷儲油量的78％。從圖2b可見，新近系（N1j，N1-2k）只有極少量的氣；而白堊系（K）又是最大的儲氣層（圖3b），儲氣量約占坳陷的74％；其次是古近系（E），儲氣量約占坳陷的15％。這些分層的模擬結果與勘探成果吻合。

2.4.3 模擬的油氣田分布與勘探成果吻合

觀察一下運聚模擬結果（圖4）與已發(fā)現(xiàn)的油氣田（圖1）的符合程度。對比表明：①圖4a上最大的油聚集帶恰好是牙哈（YH）油田（圖1），其余較小的油聚集帶也有深入勘探的前景；②圖4b上較大的氣聚集帶恰好是大北（DB）、克拉2（KL2）、羊塔克（YT）、玉東（YD）、英買7（YM7）、紅旗（HQ）、提爾根（T）及牙哈（YH）等氣田（圖1）。十分明顯，模擬結果與目前勘探成果吻合。

2.5 結論和討論

由庫車坳陷的應用實例可以得出如下結論：①采用本三維三相達西流來研究油氣二次運移是有效可行的；②庫車坳陷油氣運移的模擬結果不僅與實際情況符合，而且表明還有深入勘探的前景。

然而，由于油氣運移機理的復雜性及地質因素的不確定性，要想建立一個較完美的全定量模型并得到較好的應用效果，對地質家來說仍是一個嚴重的挑戰(zhàn)。要想使用上述方法得到一個較好的應用，應注意兩點：①盡可能把地史、熱史、成巖史、生烴史和排烴史模擬得精確，這是油氣運移定量模擬的基礎；②盡量取準油氣運移的敏感性參數(shù)。

1 Hantschel T，Kauerauf A I.Fundamentals of basin modeling and petroleum systemsmodeling［M］.Berlin：Springer-Verlag，2009

2 石廣仁.油氣運聚定量模擬技術現(xiàn)狀、問題及設想［J］.石油與天然氣地質，2009，30（1）：1～10

3 袁益讓，韓玉笈，趙衛(wèi)東，等.多層油資源運移聚集的數(shù)值模擬和實際應用［J］.應用數(shù)學和力學，2002，23（8）：827～836

4 李長峰，袁益讓.三維兩相滲流驅動問題迎風區(qū)域分裂顯隱差分法［J］.計算數(shù)學，2007，29（2）：113～136

5 Mello U T，Rodrigues JR P，Rossa A L.A control-volume finiteelementmethod for three-dimensional multiphase basin modeling［J］.Marine＆amp；Petroleum Geology，2009，26（4）：504-518

6 石廣仁.油氣盆地數(shù)值模擬方法（第三版）［M］.北京：石油工業(yè)出版社，2004

7 周興熙.庫車油氣系統(tǒng)成藏作用與成藏模式［J］.石油勘探與開發(fā)，2001，28（2）：8～10

8 趙靖舟，戴金星.庫車前陸逆沖帶天然氣成藏期與成藏史［J］.石油學報，2002，23（2）：6～10

9 王庭斌.中國氣田的成藏特征分析［J］.石油與天然氣地質，2003，24（2）：103～110

10 石廣仁，張慶春.庫車坳陷的油氣運移全定量模擬［J］.地球科學——中國地質大學學報，2004，29（4）：391～396

11 王庭斌.中國氣藏主要形成、定型于新近紀以來的構造運動［J］.石油與天然氣地質，2004，25（2）：126～132

（編輯 李 軍）

3-D 3-phase Darcy flow method and its application to the Kuqa Depression

Shi Guangren，Ma Jinshan and Chang Junhua
（PetroChina Research Institute of Petroleum Exploration＆amp；Development，Beijing 100083，China）

Themulti-phase Darcy flow method is the most complicated among all the modeling methods for secondary hydrocarbonmigration.It has threemajor technical problems：（1）itsmodel is too complicated to easily generate a solution，at leastmuch difficult than that in reservoir simulation；（2）it is not easy to ensure the convergence and stability of themodel and thus rendering a calculation of low accuracy or causing abnormal shutdown；and（3）it is very time-consuming as the calculation of a larger basinmay take several days to complete.To tackle these problems，we tried the finite volumemethod based on PEBI（perpendicular and bisection）cells，the fully implicit solution based on the Newton iteration，and the ORTHMIN（orthogonalminimum）algorithm，and came to a conclusion that the key to deal with the problems is to correctly select the allowable maximum time step for a simulated block.Finally，we developed a program of3D 3-phase Darcy Flow.An application of the software to the Kuqa Depression had gained positive results.

finite volume method，PEBI gridding，allowable maximum time-step，multi-phase Darcy flow，basin modeling，secondarymigration，Kuqa Depression

TE19 < class="emphasis_bold">文獻標識碼：A

A

0253-9985（2010）04-0403-07

2010-05-28。

石廣仁（1940—），男，教授級高級工程師，地學定量。

中石油科技部項目（2008A-0602）。