梁景偉，王立群，郎曉彬，周威皓，于萌，蘇坤偉，程超，張虎臣

（中國石油大學(xué)（北京）理學(xué)院，北京 102249）

0 引言

油藏數(shù)值模擬即用數(shù)值方法模擬油、氣、水在油藏中的運動規(guī)律，依次建立的數(shù)學(xué)模型［1］可為油藏科學(xué)合理開發(fā)、提高原油采收率等提供依據(jù)。隨著計算機技術(shù)的高速發(fā)展，油藏數(shù)值模擬也發(fā)展迅速，傳統(tǒng)規(guī)則的塊中心網(wǎng)格已不能滿足日益精細(xì)的三維地質(zhì)模型描述的需要，特別是油藏數(shù)值模擬中存在斷層、裂縫、尖滅等復(fù)雜地質(zhì)結(jié)構(gòu)以及斜井、水平井等布井方式的情況下，塊中心網(wǎng)格存在嚴(yán)重的局限性［2］。角點網(wǎng)格的出現(xiàn)為克服以上問題提供了行之有效的方法［3-4］。1992年P(guān)onting首次在油藏數(shù)值模擬中引入角點網(wǎng)格［5］，它以8個頂點確定1個不規(guī)則網(wǎng)格塊，對幾何結(jié)構(gòu)復(fù)雜的地層描述更加靈活有效［6-8］，因此，被廣泛應(yīng)用于石油地質(zhì)描述和油藏數(shù)值模擬中。

計算2個相鄰網(wǎng)格塊間的傳導(dǎo)率是油藏數(shù)值模擬的關(guān)鍵步驟之一。相鄰網(wǎng)格塊間的傳導(dǎo)率和流體性質(zhì)、巖石性質(zhì)、流體巖石之間相互作用、兩網(wǎng)格塊相連方向及網(wǎng)格幾何尺寸有關(guān)。α相某一方向的傳導(dǎo)率（Τα）可以表示為［1］

其中：α＝w，ο，g（即油、水、氣）。

式中：K為該方向的絕對滲透率；Α為該方向上相鄰網(wǎng)格塊的交接面積；h為相鄰網(wǎng)格塊中心在該方向上的距離；Krα為 α 相的相對滲透率；μα為 α 相的黏度；Βα為α相的體積系數(shù)；下標(biāo)av表示取平均值。

由于地層的復(fù)雜性和網(wǎng)格塊形狀的不規(guī)則性，相鄰網(wǎng)格間的交接面形狀（包括三角形、四邊形、五邊形、六邊形）各異，也存在交接面面積為0（尖滅）或網(wǎng)格塊體積為0等特殊情形［9-10］，因此，很大程度上加大了角點網(wǎng)格傳導(dǎo)系數(shù)的計算難度。油藏數(shù)值模擬軟件Eclipse成功地運用了角點網(wǎng)格，有效地解決了復(fù)雜幾何地層結(jié)構(gòu)的數(shù)值模擬問題，形成了該軟件的一大特色，但其技術(shù)手冊中缺乏比較詳細(xì)的計算思路與步驟。文獻(xiàn)［11］，［12］比較系統(tǒng)地推導(dǎo)了角點網(wǎng)格傳導(dǎo)系數(shù)的計算方法，但是沒有給出確定實際相鄰網(wǎng)格間交接面形狀的方法，而確定交接面形狀是計算交接面面積的前提，故該方法在油藏數(shù)值模擬中難以實現(xiàn)。本文給出了確定復(fù)雜網(wǎng)格塊交接面形狀的方法，并將其應(yīng)用到角點網(wǎng)格傳導(dǎo)系數(shù)的計算中。

傳導(dǎo)系數(shù)只與角點網(wǎng)格的幾何結(jié)構(gòu)及其內(nèi)部的絕對滲透率有關(guān)，即它在油藏數(shù)值模擬過程中是1個靜態(tài)變量。由于在不同模擬時間步中都要調(diào)用傳導(dǎo)系數(shù)，其準(zhǔn)確性直接影響數(shù)值模擬的計算結(jié)果，因此，如何準(zhǔn)確計算角點網(wǎng)格中的傳導(dǎo)系數(shù)是油藏數(shù)值模擬的關(guān)鍵技術(shù)之一，而準(zhǔn)確求出相鄰網(wǎng)格塊間交接面更是角點網(wǎng)格中傳導(dǎo)系數(shù)計算的核心所在。本文詳細(xì)推導(dǎo)了角點網(wǎng)格傳導(dǎo)系數(shù)的計算公式，充分考慮了油藏局部各種復(fù)雜幾何結(jié)構(gòu)對傳導(dǎo)系數(shù)的影響，設(shè)計出簡單有效的算法，成功地計算出相鄰網(wǎng)格塊交接面的具體形狀和在3個坐標(biāo)軸上的投影面積；通過編寫高效的程序，針對1個有代表性的油藏實例（包括斷層、尖滅等復(fù)雜情況），計算出相應(yīng)角點網(wǎng)格的傳導(dǎo)系數(shù)，并與通用商業(yè)軟件Eclipse的結(jié)果進(jìn)行對比，其結(jié)果表明本文的算法具有較高精度。

1 角點網(wǎng)格數(shù)據(jù)體導(dǎo)入與處理

考慮網(wǎng)格數(shù)為Nx×Ny×Nz的地層結(jié)構(gòu)。在大多數(shù)商業(yè)軟件中，角點網(wǎng)格包括2部分?jǐn)?shù)據(jù)體。第1部分?jǐn)?shù)據(jù)體是網(wǎng)格豎線兩端點的坐標(biāo)組成的矩陣，矩陣的列數(shù)為 6，前 3 列是豎線起點坐標(biāo)（x1，y1，z1），后 3 列為豎線終點坐標(biāo)（x2，y2，z2）；矩陣的行數(shù)代表豎線的條數(shù)，共計

在計算與網(wǎng)格塊（i，j，k ） 所在4條豎線的編號時，定義中間變量Lindex：

則與網(wǎng)格塊（i，j，k） 所在4條豎線的編號（見圖1）為

第2部分?jǐn)?shù)據(jù)體表示每1個網(wǎng)格對應(yīng)8個頂點的深度值組成的三維數(shù)據(jù)體，數(shù)據(jù)體的維數(shù)為2Nx×2Ny×2Nz。

在計算與網(wǎng)格塊（i，j，k） 對應(yīng)8個頂點的編號前，需要把以上三維數(shù)據(jù)體轉(zhuǎn)成一維數(shù)據(jù)體，然后再進(jìn)行計算。定義中間變量Pindex：

則網(wǎng)格塊（i，j，k） 對應(yīng)8個頂點的編號（見圖1）為

圖1角點網(wǎng)格（i，j，k） 相鄰4條線與8個頂點的編號示意

2 傳導(dǎo)系數(shù)的計算

2.1 x方向傳導(dǎo)系數(shù)

令I(lǐng)，J分別代表2個相鄰網(wǎng)格塊的一維坐標(biāo)。網(wǎng)格塊I的中心點到網(wǎng)格塊I與網(wǎng)格塊J相鄰面（網(wǎng)格塊I的右側(cè)面）中心點的傳導(dǎo)系數(shù)為

式中：K11為x方向的絕對滲透率；dI為網(wǎng)格塊I的中心點到網(wǎng)格塊I與網(wǎng)格塊J相鄰面 （網(wǎng)格塊 I的右側(cè)面）中心點的向量（見圖 2）為面ABCD沿著dI方向的投影面積，m2；為向量 dI的長度，m；d 為長度，m。

如果需要考慮網(wǎng)格塊I的凈毛比RNTGI，那么該傳導(dǎo)系數(shù)為

圖2 角點網(wǎng)格示意

于是網(wǎng)格塊I的中心點到網(wǎng)格塊I與網(wǎng)格塊J相鄰面（網(wǎng)格塊I的右側(cè)面）中心點的傳導(dǎo)系數(shù)為

同理可得網(wǎng)格塊J中心到網(wǎng)格塊J與網(wǎng)格塊J相鄰面（網(wǎng)格塊J的左側(cè)面）中心點的傳導(dǎo)系數(shù)TJ。TJ的計算只需將計算TI的式（12）中I標(biāo)號換為J即可：

由于AI面與AJ面的單位外法線方向正好相反，即

AIJ為網(wǎng)格塊I與網(wǎng)格塊J交接面在x，y和z方向的投影面積的絕對值。于是式（14）可表示為

如果需要考慮單位換算，則需要引入達(dá)西常數(shù)Cdarcy；如果網(wǎng)格塊I與網(wǎng)格塊J之間存在斷層，則需要引入傳導(dǎo)系數(shù)乘子TMLTx，IJ，此時該傳導(dǎo)系數(shù)為

2.2 y方向傳導(dǎo)系數(shù)

y方向傳導(dǎo)系數(shù)的計算與x方向傳導(dǎo)系數(shù)的計算方法相同，只需要把式（32）中的下標(biāo)x改成y，并把x方向的絕對滲透率改成y方向的絕對滲透率，即

其中，AIJ，DI和DJ的計算公式分別與x方向傳導(dǎo)系數(shù)中AIJ，DI和DJ的計算公式相同。

2.3 z方向傳導(dǎo)系數(shù)

z方向傳導(dǎo)系數(shù)的計算方法與x方向傳導(dǎo)系數(shù)的計算方法相同，只需要把公式（32）中的下標(biāo)x改成z，并把x方向的絕對滲透率改成z方向的絕對滲透率，但在z方向不需要考慮凈毛比RNTG，即

其中，AIJ，DI和DJ的計算公式分別與x方向傳導(dǎo)系數(shù)中AIJ，DI和DJ的計算公式相同。

由于尖滅等特殊情況的出現(xiàn)，某些網(wǎng)格塊出現(xiàn)了體積為0的特殊情況，因此，在z方向上傳導(dǎo)系數(shù)的計算必須對這些特殊情況加以處理：

1）若I網(wǎng)格塊與J網(wǎng)格塊厚度均不為0，則按上述z方向計算式（36）—（38）計算；

2）若I網(wǎng)格塊厚度為0，J網(wǎng)格塊厚度不為0，則z方向傳導(dǎo)率定義為0；

3）若I網(wǎng)格塊與J網(wǎng)格塊厚度均為0，則z方向傳導(dǎo)率定義為0；

4）若I網(wǎng)格塊厚度不為0，J網(wǎng)格塊厚度為0，則沿z方向依次向下尋找第1個厚度不為0的網(wǎng)格塊M，然后計算I網(wǎng)格塊與M網(wǎng)格塊之間的傳導(dǎo)系數(shù)，并作為I網(wǎng)格塊與J網(wǎng)格塊之間的傳導(dǎo)系數(shù)。

3 相鄰網(wǎng)格塊間交接面的計算

為了計算角點網(wǎng)格的傳導(dǎo)系數(shù)，由式（32）、式（33）、式（34）可知，其關(guān)鍵需要計算2個網(wǎng)格的交接面分別沿著x，y和z方向的投影面積。在圖2中四邊形aBCd即為2個網(wǎng)格塊的交界面，由于4個點均位于生成角點網(wǎng)格的豎線上，而豎線和各點的深度值都是已知的，因此，可以通過式（3）、式（5）求得 4個點的編號及坐標(biāo)。

盡管I與J之間交接面的幾何結(jié)構(gòu)非常復(fù)雜，但可將其分為以下3種情況處理。

3.1 相鄰側(cè)面完全重合

相鄰側(cè)面完全重合是指點A與點a，點B與點b，點C與點c，點D與點d完全重合，即網(wǎng)格塊I與網(wǎng)格塊J相鄰的側(cè)面完全重合。

由空間解釋幾何理論可知，如果空間四邊形ABCD（abcd）4 個頂點在同一平面上，那么，其在 x，y，z方向投影的有向面積Sx，Sy，Sz組成的向量S，可以用2個對角線向量的叉乘的一半表示，即

為了簡單起見，即使空間四邊形ABCD（abcd）4個頂點不在同一平面上，那么，其在x，y，z方向投影的有向面積組成的向量仍用式（38）進(jìn)行近似處理。

由于

而四邊形 ABCD（abcd）在 x，y，z方向的投影絕對面積分別為

3.2 無交接面

無交接面包括網(wǎng)格塊I與網(wǎng)格塊J相鄰的側(cè)面完全錯開，或者其中1個網(wǎng)格塊相鄰的側(cè)面面積為0，此時交接面面積為0。

3.3 交接面不完全重合

第1步：判斷線段AD與ad、AD與bc、BC與ad、BC與bc是否相交；如果相交，求出相應(yīng)交點坐標(biāo)。

由圖2可知，為了判斷每對線段是否相交，把這4對線段所在的面投影到 yoz平面上，即點 A，B，C，D，a，b，c，d 在 yoz面上的投影 A′，B′，C′，D′，a′，b′，c′，d′的 x 坐標(biāo)值為0。假設(shè)給定4個點坐標(biāo)分別為上的投影為

如果 L=0，則線段 A′D′與 a′d′平行或重合，即線段AD與ad無交點。否則，求交點：

如果L≠0，且滿足：

則線段 A′D′與 a′d′相交于（0，y，z）點。 否則，2 線段不相交。

由于線段 A′D′與 a′d′是線段 AD 與 ad 在 yoz面上的投影并相交，若線段AD與ad共面，則可以推出AD與ad相交，交點的x坐標(biāo)為

由此得到空間線段AD與ad的交點坐標(biāo)P x，y，（）z。

若線段AD與ad不共面，仍用公式（49）計算出P點坐標(biāo)，并將AP與PD連接成的折線近似地代替空間線段AD，同時將aP與Pd連接成的折線近似地代替空間線段ad。

第2步：判斷內(nèi)點。

如圖2所示，判斷點A，B和點a，b之間的位置關(guān)系時，若點A或點B位于點a，b之間，則把點A或點B稱為內(nèi)點；同理，若點a或點b位于點A，B之間，則把點a或點b稱為內(nèi)點。同時應(yīng)注意，點A，B，a，b中最多有2個內(nèi)點。用相同的方法判斷點C，D和點c，d之間的位置關(guān)系，確定這4個點中的內(nèi)點，如圖3所示。

圖3 相鄰網(wǎng)格不同交接面示意

圖3中，在AB邊上，點A在點a和點b中間，因此，A為內(nèi)點（見圖3a），同理，點b在點A和點B中間，因此，點b也為內(nèi)點（見圖3b）。

第3步：排序。

所有內(nèi)點和交點所圍成的面即為網(wǎng)格塊I與網(wǎng)格塊J的交接面。對所有內(nèi)點和線段交點，按逆時針位置關(guān)系，從A，B點所在線段開始進(jìn)行排序，確定交接面頂點位置關(guān)系。如圖2，交接面中只有內(nèi)點沒有交點，所有內(nèi)點形成四邊形C，d，a，B即為交接面。例如：圖3d中，點 A，b，C，d為內(nèi)點，點 3和點 6為交點，因此,排列順序為 A，b，3，C，d，6。

第4步：計算交接面面積。

不完全重合的交接面可分為三角形、四邊形、五邊形、六邊形4種，對于三角形和四邊形的面積計算可按上述相鄰側(cè)面完全重合情況中介紹的方法進(jìn)行計算。對于五邊形和六邊形，分別切為1個四邊形、1個三角形以及2個四邊形的組合。如：圖3c為五邊形，切為四邊形1234和三角形451，分別計算面積；圖3d為六邊形，切為四邊形1234和四邊形4561，分別計算面積。

4 實例驗證

為檢驗以上數(shù)值模型和程序的正確性，選擇某油藏實測地層數(shù)據(jù)進(jìn)行建模驗證。考慮網(wǎng)格數(shù)為24×25×12的一個真實地層，具體地層結(jié)構(gòu)見圖4。

圖4 地層結(jié)構(gòu)示意

從圖4可見，該地層幾何結(jié)構(gòu)復(fù)雜，包括斷層、尖滅等特殊情況，比較有代表性。由于該地層結(jié)構(gòu)包含了相鄰網(wǎng)格塊間交接面的各種不同幾何形態(tài)，因此傳導(dǎo)系數(shù)的計算復(fù)雜度比較高。

該地層在x，y，z 3個方向的滲透率如圖5所示。

圖5 某地層3個方向滲透率示意

圖6顯示了本文算法與Eclipse軟件就該地層所計算的傳導(dǎo)系數(shù)的絕對誤差。其中，橫軸表示z方向的層序，縱軸表示該層最大絕對誤差值。圖中數(shù)值表明，所有方向傳導(dǎo)系數(shù)最大絕對誤差為0.093，而z方向傳導(dǎo)系數(shù)的絕對誤差幾乎為0。

圖6 本文算法與Eclipse傳導(dǎo)系數(shù)絕對誤差示意

圖7顯示了本文算法與Eclipse軟件就該地層所計算的傳導(dǎo)系數(shù)的相對誤差。其中，橫軸表示z方向的層序，縱軸表示該層最大相對誤差值。圖中數(shù)值表明，所有方向傳導(dǎo)系數(shù)最大相對誤差是0.010，而z方向的傳導(dǎo)系數(shù)的相對誤差幾乎為0。

圖7 本文算法與Eclipse傳導(dǎo)系數(shù)相對誤差示意

由圖6、圖7可見，z方向的絕對誤差和相對誤差幾乎為0，x方向和y方向的2種誤差都略大于0。產(chǎn)生以上結(jié)果的原因是交接面計算的精度問題。z方向上的交接面永遠(yuǎn)是四邊形，而且上網(wǎng)格塊I的下表面永遠(yuǎn)等于下網(wǎng)格塊J的上表面，也就是說上下網(wǎng)格塊不產(chǎn)生偏移；因此，z方向上的傳導(dǎo)系數(shù)計算遠(yuǎn)比x，y方向上的傳導(dǎo)系數(shù)簡單，不同算法的數(shù)值結(jié)果差別幾乎為0。由于斷層、尖滅等特殊地質(zhì)條件的出現(xiàn)，在x，y方向上，不同相鄰網(wǎng)格塊間的交接面形狀差異性很大，不同的算法對該形狀在各個坐標(biāo)軸上投影面積的計算存在一定的差異性。但從整體上看，2種算法間的誤差比較小，從工程意義上講是可接受的。

5 結(jié)論

1）建立了基于角點網(wǎng)格傳導(dǎo)系數(shù)數(shù)值計算模型。該模型充分考慮了角點網(wǎng)格塊間各種不同交接面對傳導(dǎo)系數(shù)計算的影響，在此基礎(chǔ)上詳細(xì)推導(dǎo)了傳導(dǎo)系數(shù)的計算公式，并通過具體實例驗證了算法的準(zhǔn)確性。

2）將以上算法應(yīng)用于其他實際油藏，并與Eclipse軟件計算結(jié)果進(jìn)行比較，其差異性均較小，說明本文算法具有廣泛的適用性。

［1］Chen Z.Reservoir simulation:Mathematical techniques in oil recovery［M］.Philadelphia：Siam，2007：103-129.

［2］Helnemann Z E，Brand C W，Munka M，et al.Modeling reservoir geometry with irregular grids［R］.SPE 18412，1991.

［3］Ding Y，Lemonnier P.Use of corner point geometry in reservoir simulation［R］.SPE 29933，1995.

［4］Ding Y.Permeability upscaling on corner-point geometry in the nearwell region［R］.SPE 81431，2003.

［5］Ponting D K.Corner point geometry in reservoir simulation.The mathematics of oil recovery ［M］.Oxford：Oxford University Press，1992：45-65.

［6］葉繼根，吳向紅，朱怡翔，等.大規(guī)模角點網(wǎng)格計算機輔助油藏模擬歷史擬合方法研究［J］.石油學(xué)報，2007，28（2）：83-86.

［7］毛小平，張志庭，錢真.用角點網(wǎng)格模型表達(dá)地質(zhì)模型的剖析及在油氣成藏過程模擬中的應(yīng)用［J］.地質(zhì)學(xué)刊，2012，36（3）：265-273.

［8］Juanes R，Matringe S F，Thomas L K.Implementation and application of a hybrid multipoint flux approximation for reservoir simulation on corner-point grids［R］.SPE 95928，2005.

［9］Aarnes J，Krogstad S，Lie K A.Multiscale mixed/mimetic methods on corner-point grids［J］.Computational Geosciences，2008，12（3）：297-315.

［10］Sammon P H.Calculation of convective and dispersive flows for complex corner point grids［R］.SPE 62929，2000.

［11］劉紅衛(wèi)，李愛華，趙國忠.角點網(wǎng)格傳導(dǎo)率計算技術(shù)研究［J］.大慶石油地質(zhì)與開發(fā)，2005，24（6）：35-36.

［12］計秉玉，趙國忠，李潔.多學(xué)科集成化油藏研究方法與應(yīng)用［M］.石油工業(yè)出版社，2009：85-96.