饒 翔， 程林松， 曹仁義*， 安小平， 雷征東

(1.中國石油大學(xué)(北京)石油工程教育部重點(diǎn)實(shí)驗(yàn)室，北京 102249；2.中國石油大學(xué)(北京)石油工程學(xué)院，北京 102249；3.中國石油長慶油田分公司勘探開發(fā)研究院，西安 710018；4.中國石油勘探開發(fā)研究院，北京 100083)

在開展復(fù)雜裂縫網(wǎng)絡(luò)傳質(zhì)的數(shù)值模擬時(shí)，為了能夠準(zhǔn)確刻畫油藏中裂縫的形態(tài)，離散裂縫模型(DFM)需要生成高質(zhì)量的非結(jié)構(gòu)網(wǎng)格去匹配裂縫在油藏空間中的幾何形態(tài)。對于在三維油藏中的大量傾斜裂縫及其形成的復(fù)雜縫網(wǎng)，與之匹配的高質(zhì)量非結(jié)構(gòu)網(wǎng)格的生成是十分困難的。如果需要考慮裂縫的動(dòng)態(tài)行為，則需要在每一時(shí)間步重新生成非結(jié)構(gòu)網(wǎng)格去匹配該時(shí)間步的縫網(wǎng)形態(tài)，會(huì)顯著降低計(jì)算效率。而嵌入式離散裂縫模型(EDFM)僅需要利用結(jié)構(gòu)化網(wǎng)格剖分基質(zhì)系統(tǒng)，形成背景網(wǎng)格，將顯式表征的裂縫嵌入到背景網(wǎng)格之中形成裂縫網(wǎng)格，再獲取網(wǎng)格之間的連接及相應(yīng)的傳導(dǎo)系數(shù)，可以解決高質(zhì)量非結(jié)構(gòu)網(wǎng)格生成的問題，在考慮裂縫擴(kuò)展或增刪裂縫的情況下無需更新背景網(wǎng)格，模擬效率高，操作靈活。Lee等[1]首次提出將裂縫嵌入到基質(zhì)網(wǎng)格中，采用與水平井相似的方式將裂縫視為基質(zhì)網(wǎng)格中的源匯項(xiàng)來處理基質(zhì)與裂縫之間的傳質(zhì)。隨后，Li等[2]擴(kuò)展了前述的初步想法首次提出了EDFM。在EDFM中，由于結(jié)構(gòu)化網(wǎng)格的使用，一般采用既能滿足局部物質(zhì)守恒，又簡單易行的有限體積、有限差分或者模擬有限差分等方法去獲取滲流控制方程的離散格式。Hajibeygi等[3]利用循環(huán)多尺度有限體積方法建立了相應(yīng)的EDFM，可以有效提高模擬計(jì)算效率。Moinfar等[4-5]首次將EDFM從二維擴(kuò)展到三維，并簡單分析了天然裂縫和水力壓裂縫動(dòng)態(tài)行為對滲流場及水平井產(chǎn)能的影響。EDFM已被用于開展多種滲流問題中復(fù)雜縫網(wǎng)傳質(zhì)的數(shù)值模擬研究[6-12]。Zhang等[13]、Yan等[14]利用模擬有限差分方法離散滲流控制方程，達(dá)到同時(shí)求解速度場和壓力場的目的；Cao等[15]基于邊界積分方程得到基質(zhì)網(wǎng)格與裂縫網(wǎng)格之間傳質(zhì)量的離散格式；Rao等[16]基于格林元方法中對方程中含時(shí)項(xiàng)處理的思想提出了新的基質(zhì)網(wǎng)格與裂縫網(wǎng)格傳質(zhì)量的近似格式；Rao等[17]利用EDFM研究水侵層對注水吞吐過程的影響。目前，EDFM已被廣泛應(yīng)用于各種復(fù)雜滲流問題的模擬研究中，然而Tene等[18]研究發(fā)現(xiàn)EDFM無法有效處理裂縫網(wǎng)格滲透率低于所在基質(zhì)網(wǎng)格滲透率的情況，如油藏中的斷層、隔夾層等。Jiang等[19]的研究也表明EDFM在多相流橫穿裂縫時(shí)會(huì)出現(xiàn)明顯誤差。造成EDFM局限性的本質(zhì)原因在于僅將裂縫處理為其所在基質(zhì)網(wǎng)格中的源匯項(xiàng)，一般僅能適用于多段壓裂水平井生產(chǎn)或注入過程的數(shù)值模擬。為了解決EDFM的上述局限性，Tene等[18]提出了基于投影的嵌入式離散裂縫模型(pEDFM)；Jiang等[19]對pEDFM做出了一些改進(jìn)，改進(jìn)之處主要是提出了裂縫單元投影面的選取原則以及修正了傳導(dǎo)系數(shù)的計(jì)算公式。但pEDFM與EDFM的相似之處在于都是將裂縫嵌入到被結(jié)構(gòu)化網(wǎng)格剖分的基質(zhì)系統(tǒng)獲取裂縫網(wǎng)格的分布，建立網(wǎng)格之間的連接，不同的是pEDF較EDFM添加了一類新的連接。

EDFM和pEDFM的前處理主體部分基本相同，從二維模型升級(jí)到三維模型的本質(zhì)即是前處理算法的維度升級(jí)，由于三維油藏中裂縫面可以任意傾斜或者具有復(fù)雜的幾何形狀，這給前處理過程帶來了挑戰(zhàn)，并且在以往的三維EDFM中一般只有矩形平面縫，且均未給出相應(yīng)的前處理算法細(xì)節(jié)。為此，提出一種可適用于由復(fù)雜形狀傾斜縫構(gòu)成的復(fù)雜裂縫網(wǎng)絡(luò)的EDFM前處理算法，以期依據(jù)該前處理算法實(shí)現(xiàn)相應(yīng)的三維EDFM。

1 裂縫面的參數(shù)化

在EDFM中，對裂縫采用降維處理，即在三維油藏中，用一個(gè)二維平面或曲面刻畫一條裂縫(記作裂縫面)。一般情況下，三維油藏用三維歐式空間刻畫(油藏空間尺度不大，無需考慮地面曲率的影響)，裂縫面是一個(gè)二維曲面。根據(jù)經(jīng)典的曲面論，二維裂縫曲面可以用該曲面的第一、第二基本形式和該面上某點(diǎn)的向徑唯一確定。然而，如果裂縫面用二維曲面表征，則難以獲取裂縫面與基質(zhì)網(wǎng)格之間的連接關(guān)系及相應(yīng)裂縫網(wǎng)格的幾何參數(shù)，并且在EDFM或者DFM中，裂縫一般是二維帶邊平面或者一維直線段，因此利用二維平面來刻畫三維油藏中的裂縫，即一條裂縫對應(yīng)著一個(gè)二維平面。而要確定一個(gè)三維歐式空間(油藏)中的二維平面，則需要該平面上某點(diǎn)A0的向徑(記作基準(zhǔn)向量β0)和該裂縫面上的兩個(gè)線性無關(guān)(或正交)的向量β1和β2，這兩個(gè)線性無關(guān)的向量實(shí)際上構(gòu)成了該裂縫所在平面以A0為原點(diǎn)的局部坐標(biāo)系的兩個(gè)基向量，即該裂縫所在平面上的點(diǎn)可以表示為

r=β0+uβ1+vβ2

(1)

式(1)中：r為裂縫所在平面上點(diǎn)的矢徑；u和v分別為β1和β2對應(yīng)的參數(shù)，即偶對 (u,v)在局部坐標(biāo)系{A0，β1,β2}中的坐標(biāo)，與該裂縫所在平面上的點(diǎn)一一對應(yīng)。

為了進(jìn)一步確定裂縫在該平面上的形狀，則需要給出u和v需要滿足的方程函數(shù)關(guān)系f(u,v)≤0，其中，f(u,v)=0確定裂縫的邊界形狀，f(u,v)<0即是裂縫內(nèi)部。

基于以上的裂縫面參數(shù)化表示，得出具體的前處理算法如下。

2 前處理具體算法

EDFM實(shí)質(zhì)上是將DFM生成非結(jié)構(gòu)網(wǎng)格匹配裂縫幾何形態(tài)的困難轉(zhuǎn)變?yōu)楂@取裂縫網(wǎng)格分布及網(wǎng)格之間連接關(guān)系的困難，這種轉(zhuǎn)變使得EDFM避免了在復(fù)雜縫網(wǎng)情況下高質(zhì)量非結(jié)構(gòu)網(wǎng)格生成的困難，并且只需要一套預(yù)先設(shè)定的固定的背景網(wǎng)格，使得在求解大量問題時(shí)的計(jì)算效率更高。因此獲取裂縫網(wǎng)格分布及網(wǎng)格之間連接關(guān)系是嵌入式離散裂縫模型前處理的關(guān)鍵，同時(shí)也說明從二維EDFM到三維EDFM維度升級(jí)的本質(zhì)困難在于前處理算法的維度升級(jí)。

在獲取裂縫網(wǎng)格分布及網(wǎng)格之間連接關(guān)系這個(gè)問題上，三維EDFM比二維EDFM要復(fù)雜得多，因?yàn)樵诙SEDFM中是二維的基質(zhì)網(wǎng)格和一維直線段裂縫，故只需要求解一維直線段裂縫與基質(zhì)網(wǎng)格各邊的交點(diǎn)以及不同直線段裂縫之間的交點(diǎn)即可確定裂縫網(wǎng)格的分布以及網(wǎng)格之間的連接關(guān)系，但在三維EDFM中是三維的基質(zhì)網(wǎng)格和二維的裂縫面，裂縫面與三維基質(zhì)網(wǎng)格之間以及不同裂縫面之間的幾何關(guān)系復(fù)雜。對此，三維EDFM前處理算法的基本思路分為六個(gè)步驟：①計(jì)算裂縫面(裂縫面指裂縫降維處理后的幾何實(shí)體；裂縫所在平面是指裂縫面實(shí)體所在的三維歐式空間中的無限大平面)與基質(zhì)網(wǎng)格線的交點(diǎn)；②計(jì)算裂縫面邊界與基質(zhì)網(wǎng)格面的交點(diǎn)；③確定預(yù)裂縫網(wǎng)格分布；④求解不同裂縫面之間的交線方程；⑤確定最終的裂縫網(wǎng)格分布；⑥確定網(wǎng)格之間連接關(guān)系。

按照上述的六個(gè)步驟，具體的算法細(xì)節(jié)如下。

2.1 計(jì)算裂縫面與基質(zhì)網(wǎng)格線的交點(diǎn)

假設(shè)該裂縫面的參數(shù)化形式如下：

r=β0+uβ1+vβ2，f(u,v)≤0

(2)

式(2)中：β0=(β01,β02,β03)，β1=(β11,β12,β13)，β2=(β21,β22,β23)。

則計(jì)算裂縫面與基質(zhì)網(wǎng)格線交點(diǎn)的具體算法步驟如下。

步驟1首先判斷該裂縫面是否垂直于x-y平面Γxy、y-z平面Γyz或x-z平面Γxz，判斷方法是：先計(jì)算該裂縫面法向量n=β1β2，若n·(0,0,1)=0，則該裂縫面垂直于Γxy，若n·(1,0,0)=0，則該裂縫面垂直于Γyz，若n·(0,1,0)=0，則該裂縫面垂直于Γxz。

步驟2計(jì)算網(wǎng)格線與裂縫所在平面的交點(diǎn)，以垂直于x-y平面的z方向網(wǎng)格線為例說明，對于某一條z方向網(wǎng)格線，其x和y是確定的，不妨假設(shè)為x0和y0，則將x0和y0代入式(2)，可以得到線性方程組：

(3)

若裂縫所在平面不垂直于x-y平面[圖1(a)]，根據(jù)幾何關(guān)系可知裂縫所在平面與該z方向網(wǎng)格線必有唯一交點(diǎn)，即上述線性方程組系數(shù)矩陣非奇異必有唯一解。解方程組得到(u0,v0)后，再判斷(u0,v0)是否滿足f(u0,v0)≤0，若不滿足，則說明該z方向網(wǎng)格線與裂縫面無交點(diǎn)，若滿足，則說明(u0,v0)對應(yīng)的點(diǎn)是該z方向網(wǎng)格線與裂縫面的交點(diǎn)，則計(jì)算z0=β03+u0β13+v0β23，存儲(chǔ)(x0,y0,z0)和(u0,v0)分別為該交點(diǎn)的三維歐式空間坐標(biāo)和二維裂縫面上局部坐標(biāo)。

圖1 裂縫面與基質(zhì)網(wǎng)格線之間的三種幾何關(guān)系

若裂縫所在平面垂直于x-y平面，根據(jù)幾何關(guān)系可知上述線性方程組系數(shù)矩陣奇異，或者無解,如圖1(b)所示，對應(yīng)于該網(wǎng)格線與裂縫所在平面沒有交點(diǎn)；或者無窮多組解，如圖1(c)所示，對應(yīng)于該網(wǎng)格線包含在裂縫所在平面有無窮多個(gè)交點(diǎn)。此時(shí)需要判斷是上述哪種情況，方法是任意取一個(gè)z0值，將y0和z0代入式(2)，可以得到一個(gè)線性方程組：

(4)

若式(4)中方程組的系數(shù)矩陣奇異，即說明原z方向網(wǎng)格線與裂縫所在平面沒有交點(diǎn)，即無解；反之，則說明該z方向網(wǎng)格線包含在該裂縫平面內(nèi)，是第二種情況。在無解情況下，由于沒有交點(diǎn)，因此無需進(jìn)行其他操作；在無窮多組解下，由于該z網(wǎng)格方向線包含在裂縫所在平面，因此通過二分法結(jié)合f(u0,v0)≤0的條件來計(jì)算z0的取值范圍，從而將z0取值范圍內(nèi)的同時(shí)也屬于基質(zhì)網(wǎng)格頂點(diǎn)(格點(diǎn))的(x0,y0,z0)以及(x0,y0,z0,max)、(x0,y0,z0,min)存儲(chǔ)為該z方向網(wǎng)格線與裂縫面的交點(diǎn)，同時(shí)存儲(chǔ)這些交點(diǎn)對應(yīng)二維裂縫面上的局部坐標(biāo)(u0,v0)[該局部坐標(biāo)通過方程(3)求解得到]。

同理，對于x方向坐標(biāo)線和y方向坐標(biāo)線，采用類似的方法計(jì)算得到坐標(biāo)線與裂縫面的交點(diǎn)。

2.2 計(jì)算裂縫面邊界與基質(zhì)網(wǎng)格面的交點(diǎn)

依然假設(shè)該裂縫面的參數(shù)化表示如式(2)所示，則計(jì)算裂縫面邊界與基質(zhì)網(wǎng)格面交點(diǎn)的具體算法步驟如下。

步驟1首先判斷裂縫面Γf是否平行于x-y平面Γxy、y-z平面Γyz或x-z平面Γxz，其判斷方法類似于2.1節(jié)的第一步；

步驟2接下來以某一x-y基質(zhì)網(wǎng)格面Γ1(其z值是一個(gè)常數(shù)，記為z0)為例介紹計(jì)算裂縫面Γf與Γ1交點(diǎn)的算法。

圖2 裂縫面與基質(zhì)網(wǎng)格面之間的四種幾何關(guān)系

若裂縫面Γf不與Γxy平行，若存在裂縫面某邊界Γ1Γf在基質(zhì)網(wǎng)格面Γ1內(nèi)，如圖2(a)所示，因?yàn)檫@兩個(gè)面均是平面，故該邊界是一條直線段，且該直線段一般與Γ1的網(wǎng)格線有交點(diǎn)，即這些交點(diǎn)在2.1節(jié)中已計(jì)算出來，此時(shí)無需任何操作；同時(shí)因?yàn)樵撨吔缡且粭l直線，故方程(5)是線性方程組，因此若該方程組的系數(shù)矩陣是否奇異則可以判斷?1Γf在基質(zhì)網(wǎng)格面Γ1內(nèi)，如圖2(b)～圖2(d)所示，在這種情況下由幾何關(guān)系可知該裂縫面邊界與Γ1至多有有限個(gè)交點(diǎn)，裂縫面Γf邊界方程是f(u,v)=0，因此交點(diǎn)對應(yīng)的裂縫面Γf上的局部坐標(biāo)(u0,v0)滿足下述方程組：

(5)

求解式(5)即可解出交點(diǎn)的二維裂縫面局部坐標(biāo)(u0,v0)，再將(u0,v0)代入到裂縫面Γf參數(shù)化形式中，計(jì)算交點(diǎn)對應(yīng)的三維歐式空間坐標(biāo)(x0,y0,z0)。特殊的，當(dāng)裂縫面Γf是矩形、平行四邊形或者其他多邊形時(shí)，裂縫面Γf邊界是由直線段構(gòu)成的，因此上述方程組實(shí)際上是線性方程組，易于求解。若裂縫面Γf邊界是橢圓或者更高次代數(shù)曲線時(shí)，則需要求解二次或更高次的二元代數(shù)方程組，此時(shí)有兩種方法，第一種是利用NR迭代計(jì)算，但需要注意若迭代不收斂，則說明裂縫面Γf邊界與Γxy沒有交點(diǎn)；第二種是先利用方程組中的第一個(gè)線性方程，用u0表示v0，或者用v0表示u0，再代入f(u0,v0)=0中得到一元二次或更高次方程求解實(shí)數(shù)根，若沒有實(shí)數(shù)根，則說明裂縫面Γf邊界與Γxy沒有交點(diǎn)。此處是三維EDFM前處理中矩形縫與橢圓形等復(fù)雜形狀縫在幾何參數(shù)計(jì)算中的唯一不同之處。

若裂縫面平行于Γxy，則可知裂縫面邊界與Γ1要么沒有交點(diǎn)，此時(shí)無需任何操作；要么無窮多個(gè)交點(diǎn)(即該裂縫面邊界全包含在Γ1內(nèi))，這種情況在EDFM中一般不會(huì)出現(xiàn)，因?yàn)榇藭r(shí)裂縫面與基質(zhì)網(wǎng)格面重合，即裂縫面在幾何形態(tài)上成為了基質(zhì)網(wǎng)格的邊界。

2.3 確定預(yù)裂縫網(wǎng)格分布

在前述基質(zhì)網(wǎng)格線與每個(gè)裂縫面之間以及基質(zhì)網(wǎng)格面與每個(gè)裂縫面邊界之間的所有交點(diǎn)計(jì)算得到后，判斷同一裂縫面上的所有交點(diǎn)所屬的基質(zhì)網(wǎng)格編號(hào)，在同一基質(zhì)網(wǎng)格中的所有交點(diǎn)構(gòu)成的一個(gè)凸多邊形，定義為該裂縫面上的一個(gè)預(yù)裂縫網(wǎng)格(此處是預(yù)裂縫網(wǎng)格，是因?yàn)樽罱K的裂縫網(wǎng)格還需要考慮不同裂縫面之間的交線對預(yù)裂縫網(wǎng)格的進(jìn)一步剖分)，一個(gè)裂縫面上的全體預(yù)裂縫網(wǎng)格構(gòu)成該裂縫面上的預(yù)裂縫網(wǎng)格分布。如圖3所示，點(diǎn)1、2、3、4、5、6、7是通過2.1節(jié)算法計(jì)算得到的該橢圓裂縫面與基質(zhì)網(wǎng)格線的交點(diǎn)，點(diǎn)8、9、10、11、12、13、14、15、16、17是通過2.2節(jié)算法計(jì)算得到的橢圓形裂縫面與基質(zhì)網(wǎng)格面的交點(diǎn)，然后屬于同一基質(zhì)網(wǎng)格的節(jié)點(diǎn)構(gòu)成一個(gè)預(yù)裂縫網(wǎng)格，從而得到12個(gè)預(yù)裂縫網(wǎng)格，分別是：1-2-8、1-2-9、2-3-10-8、2-3-11-9、3-4-12-10、3-4-13-11、4-5-14-12、4-5-15-13、5-6-16-14、5-6-17-1、6-7-16、6-7-17。

圖3 橢圓形裂縫面預(yù)裂縫網(wǎng)格示意圖

需要指出以下問題。

(1)其中凸多邊形的生成及相應(yīng)的交點(diǎn)排列順序是利用MATLAB中DelaunayTri函數(shù)和convexhull函數(shù)得到的，需要注意的是，根據(jù)這兩個(gè)函數(shù)的使用規(guī)則，并不能夠處理這些交點(diǎn)的三維坐標(biāo)(x,y,z)，因此需要將每個(gè)交點(diǎn)對應(yīng)的二維裂縫面上的局部坐標(biāo)(u,v)代入函數(shù)中處理得到凸多邊形，這也說明將二維裂縫面參數(shù)化從而只需要用兩個(gè)參數(shù)u、v去刻畫裂縫面的處理是十分必要的。

(2)如果一個(gè)基質(zhì)網(wǎng)格內(nèi)包含某裂縫面交點(diǎn)的數(shù)量小于3，則無法構(gòu)成一個(gè)凸多邊形，則認(rèn)為該裂縫面沒有屬于這個(gè)基質(zhì)網(wǎng)格中的預(yù)裂縫網(wǎng)格。

按照上述方法，可以求出所有裂縫面上的預(yù)裂縫網(wǎng)格分布，其中每個(gè)預(yù)裂縫網(wǎng)格是該凸多邊形頂點(diǎn)編號(hào)按照頂點(diǎn)順序組成的有序數(shù)組來表征和存儲(chǔ)的。

2.4 求解不同裂縫面之間的交線方程

由上述算法可知，預(yù)裂縫網(wǎng)格分布僅是基質(zhì)網(wǎng)格與每一個(gè)單獨(dú)裂縫面相交形成的裂縫網(wǎng)格分布，而實(shí)際情況中很可能還存在多個(gè)裂縫面相交，繼而將預(yù)裂縫網(wǎng)格分布繼續(xù)剖分的情況，因此接下來首先是要判斷各個(gè)裂縫面之間是否會(huì)相交，如果相交，則需要求解交線的方程。具體算法如下。

不妨假設(shè)是對裂縫面1和裂縫面2進(jìn)行處理，其裂縫面參數(shù)化形式分別如下。

裂縫面1：

r=α0+uα1+vα2，f1(u,v)≤0

(6)

式(6)中：α0=(α01,α02,α03)；α1=(α11,α12,α13)；α2=(α21,α22,α23)。

裂縫面2：

r=β0+aβ1+bβ2，f2(u,v)≤0

(7)

式(7)中：β0=(β01,β02,β03)；β1=(β11,β12,β13)；β2=(β21,β22,β23)。

通過計(jì)算裂縫面1和裂縫面2的法向量n1=α1×α2和n2=β1×β2，若n1×n2=0，則裂縫面1與裂縫面2平行，即這兩個(gè)裂縫面之間不會(huì)相交；若n1×n2≠0，則裂縫面1與裂縫面2不平行，即裂縫面1所在平面和裂縫面2所在平面(此處是裂縫面所在平面，而不是裂縫面)必會(huì)存在交線。然后計(jì)算這兩個(gè)裂縫面所在平面的交線方程，提出的方法如下：

首先證明一個(gè)引理。

引理如圖4所示，若裂縫面1和裂縫面2不平行，則直線r1=α0+uα1和r2=α0+vα2中至少存在一條直線與裂縫面2所在平面存在唯一交點(diǎn)。

證明若直線r1=α0+uα1和r2=α0+vα2均與裂縫面2所在平面沒有交點(diǎn)或者有無窮多個(gè)交點(diǎn)，則說明這兩條直線的方向向量α1和α2均平行于裂縫2所在平面，又因?yàn)棣?和α2線性無關(guān)，故推出裂縫面1所在平面與裂縫面2所在平面平行，產(chǎn)生矛盾，因此假設(shè)不成立，即直線r1=α0+uα1和r2=α0+vα2中至少存在一條直線與裂縫面2所在平面存在唯一交點(diǎn)。

圖4 引理證明過程示意圖

基于上述引理，如圖5所示，不妨設(shè)直線r1=α0+uα1與裂縫面2所在平面存在唯一交點(diǎn)，結(jié)合裂縫面2的參數(shù)化形式得到方程α0+uα1=β0+aβ1+bβ2存在唯一解，即下述線性方程組存在唯一解：

(8)

圖5 交線方程計(jì)算示意圖

求解上述線性方程組得到即可得到(u,a,b)，繼而計(jì)算出裂縫面1所在平面上直線r1=α0+uα1與裂縫面2所在平面的交點(diǎn)，記作點(diǎn)A。將直線r1=α0+uα1沿α2方向平移得到直線r3=(α0+α2)+uα1，可知直線r3仍是裂縫面1所在平面上的直線并與直線r1平行，因此直線r3與裂縫面2所在平面也存在唯一交點(diǎn)，即方程(α0+α2)+uα1=β0+aβ1+bβ2也存在唯一解，類似得解如下線性方程組即可計(jì)算出直線r3與裂縫面2所在平面的交點(diǎn)，記作點(diǎn)B。

(9)

點(diǎn)A和點(diǎn)B都是裂縫面1所在平面和裂縫面2所在平面之間的交點(diǎn)，又因?yàn)閮牲c(diǎn)確定一條直線，故AB連線所在直線即是裂縫面1所在平面與裂縫面2所在平面之間的交線，從而可由A、B這兩點(diǎn)的坐標(biāo)確定出交線的方程。

同理，通過上述算法可以求解出各裂縫面所在平面之間的交線方程。

2.5 確定最終的裂縫網(wǎng)格分布

以裂縫面1為例說明具體的過程，若某裂縫面(稱為裂縫面2)所在平面與裂縫面1所在平面存在交線且交線方程為r=γ0+mγ1。接下來就是要判斷該交線會(huì)不會(huì)繼續(xù)剖分某些預(yù)裂縫網(wǎng)格，具體算法如下。

以裂縫面1上的某一個(gè)預(yù)裂縫網(wǎng)格(稱為預(yù)裂縫網(wǎng)格1)為例說明，假設(shè)該預(yù)裂縫網(wǎng)格是一個(gè)凸四邊形，記作ABCD，并且ABCD是其構(gòu)成凸多邊形的連點(diǎn)順序。

步驟1首先判斷邊AB、BC、CD、DA是否與交線平行，若平行，則說明該邊與交線沒有交點(diǎn)或者該邊包含在交線內(nèi)，這種情況下無需任何操作，若不平行，則該邊所在直線與該交線必有唯一交點(diǎn)E，求解出該交點(diǎn)E，并判斷點(diǎn)E是否在該邊上(因?yàn)樵撨吺且粭l有限長線段)，若在該邊上，則存儲(chǔ)交線與該邊的交點(diǎn)E。

步驟2若存儲(chǔ)的交點(diǎn)數(shù)小于2，則說明該交線對該“預(yù)裂縫網(wǎng)格”沒有影響，無需繼續(xù)剖分；若存儲(chǔ)的交點(diǎn)數(shù)等于2(不可能大于2，因?yàn)椤邦A(yù)裂縫網(wǎng)格是凸多邊形”)，則需要根據(jù)存儲(chǔ)的交點(diǎn)所在的具體位置將該“預(yù)裂縫網(wǎng)格”繼續(xù)剖分成兩個(gè)子凸多邊形裂縫網(wǎng)格。以圖6為例，凸四邊形ABCD是某裂縫面上的一個(gè)預(yù)裂縫網(wǎng)格，交線與該網(wǎng)格相交于E、F兩點(diǎn)，故將ABCD剖分為AEFD和BCFE這兩個(gè)裂縫網(wǎng)格。

圖6 某預(yù)裂縫網(wǎng)格被交線剖分為兩個(gè)裂縫網(wǎng)格

步驟3對裂縫面1上的所有“預(yù)裂縫網(wǎng)格”按照步驟1、步驟2循環(huán)，可以得到裂縫面2所在平面與裂縫面1所在平面的交線r=γ0+mγ1對裂縫面1上“預(yù)裂縫網(wǎng)格”繼續(xù)剖分后的網(wǎng)格分布。

步驟4對與裂縫面1所在平面存在交線的裂縫面循環(huán)步驟3，則可以得到裂縫面1上最終的裂縫網(wǎng)格分布。

同理可以得到其他裂縫面上最終的裂縫網(wǎng)格分布。

如圖7(a)所示，兩條裂縫相交后，兩條裂縫面上的“預(yù)裂縫網(wǎng)格分布”被交線進(jìn)一步剖分，圖7(b)、圖7(c)分別是兩條裂縫面上最終裂縫網(wǎng)格分布的二維視角圖。

圖7 兩條裂縫相交后裂縫網(wǎng)格分布

2.6 確定網(wǎng)格之間的連接關(guān)系

在EDFM中，存在四類網(wǎng)格之間的連接，分別是：①兩相鄰基質(zhì)網(wǎng)格之間的連接；②同一裂縫面上且位于不同基質(zhì)網(wǎng)格內(nèi)的兩相鄰裂縫網(wǎng)格之間的連接；③同一基質(zhì)網(wǎng)格內(nèi)裂縫網(wǎng)格之間的連接；④基質(zhì)網(wǎng)格與其所包含的裂縫網(wǎng)格之間的連接。

2.6.1 第一類連接

僅涉及相鄰基質(zhì)網(wǎng)格之間的連接，是比較平凡的，剩余三類均與裂縫網(wǎng)格有關(guān)，2.6節(jié)針對這三類連接分別簡要給出確定相應(yīng)連接關(guān)系的方法。

2.6.2 第二類連接

基于每個(gè)裂縫面上的裂縫網(wǎng)格分布，其中，每個(gè)裂縫網(wǎng)格是由該凸多邊形網(wǎng)格頂點(diǎn)編號(hào)按照頂點(diǎn)順序組成的有序數(shù)組來表征的，記這個(gè)有序數(shù)組為該裂縫網(wǎng)格的表征數(shù)組。對于同一裂縫面上的兩個(gè)裂縫網(wǎng)格(記為e1、e2)，首先計(jì)算這兩個(gè)裂縫網(wǎng)格表征數(shù)組中相同頂點(diǎn)編號(hào)的數(shù)量(記為n)，若n<2，則e1和e2之間沒有公共邊，即e1和e2無連接；由于都是凸多邊形裂縫網(wǎng)格，因此不可能出現(xiàn)n>2的情況；若n=2，不妨設(shè)e1和e2公共的兩個(gè)頂點(diǎn)分別是A、B,則e1和e2之間具有公共邊AB，同時(shí)可以計(jì)算出該公共邊長度(記為lAB)，e1和e2中心到AB的距離(分別記為de1,AB，de2,AB)，并可以計(jì)算后續(xù)塊中心有限體積離散格式中傳導(dǎo)系數(shù)所需的一些幾何參數(shù)值。

2.6.3 第三類及第四類連接

在第三類連接中，有多個(gè)裂縫網(wǎng)格在同一個(gè)基質(zhì)網(wǎng)格中，這種情況是由于有不同的裂縫面在該基質(zhì)網(wǎng)格中相交造成的。在三維基質(zhì)網(wǎng)格和二維裂縫面情況下，該基質(zhì)網(wǎng)格與包含的裂縫網(wǎng)格之間的幾何關(guān)系遠(yuǎn)比二維基質(zhì)網(wǎng)格和一維裂縫線段的情況復(fù)雜得多，缺乏一個(gè)絕對精確的確定裂縫網(wǎng)格之間連接關(guān)系及相應(yīng)幾何因子計(jì)算的方法。對此，假設(shè)在同一基質(zhì)網(wǎng)格中的所有裂縫網(wǎng)格之間均存在連接(在實(shí)際地質(zhì)條件下，這種假設(shè)也是合理的，因?yàn)榛|(zhì)網(wǎng)格尺寸一般不大)。基于上述假設(shè)，由于2.3、2.5節(jié)中的算法可以獲取每個(gè)基質(zhì)網(wǎng)格內(nèi)包含的裂縫網(wǎng)格編號(hào)，因此便可以建立在同一基質(zhì)網(wǎng)格中的這些裂縫網(wǎng)格之間的連接，同時(shí)也可以建立該基質(zhì)網(wǎng)格與所包含的每個(gè)裂縫網(wǎng)格之間的連接。

3 數(shù)值實(shí)例

根據(jù)前處理算法，建立能夠適用于任意形狀傾斜縫的三維EDFM，并給出一個(gè)三維油藏復(fù)雜縫網(wǎng)情況下單相原油衰竭開發(fā)的數(shù)值實(shí)例。

如圖8所示，油藏工區(qū)大小為1 000 m×700 m×20 m，基質(zhì)系統(tǒng)網(wǎng)格的總數(shù)是100×70×2，每個(gè)基質(zhì)網(wǎng)格的尺寸為10 m×10 m×10 m。有一口多段壓裂水平井定井底流壓(10 MPa)生產(chǎn)，體積壓裂形成的有效縫網(wǎng)包含10條裂縫，其中5條人工主裂縫和5條次生裂縫，其中有兩條橢圓形主裂縫?？梢钥吹竭@10條裂縫具有不同程度的傾斜(紅色網(wǎng)格表示形成的裂縫網(wǎng)格)，且在縱向上貫穿的層位各有不同。油藏及流體的基本物性參數(shù)如表1所示。

圖8 油藏模型示意圖

表1 油藏及流體的基本物性參數(shù)

圖9所示分別為生產(chǎn)三年的日產(chǎn)油量、累產(chǎn)油量、平均地層壓力及采出程度曲線，圖10所示分別為生產(chǎn)20、500 d后油藏上下兩層的壓力分布。由圖10可知，依據(jù)提出的適用于任意形狀傾斜縫的前處理算法而建立的三維EDFM能夠有效處理三維油藏復(fù)雜縫網(wǎng)傳質(zhì)的數(shù)值模擬。

4 結(jié)論

(1)首先將裂縫面參數(shù)化，然后基于裂縫面參數(shù)化形式，將三維EDFM前處理過程分為六個(gè)步驟：①計(jì)算裂縫面與基質(zhì)網(wǎng)格線的交點(diǎn)；②計(jì)算裂縫面邊界與基質(zhì)網(wǎng)格面的交點(diǎn)；③確定預(yù)裂縫網(wǎng)格分布；④求解不同裂縫面之間的交線方程；⑤確定最終的裂縫網(wǎng)格分布；⑥確定網(wǎng)格之間連接關(guān)系。

(2)針對上述的六個(gè)步驟分別提出具體的算法，該套算法能夠?qū)崿F(xiàn)任意形狀傾斜縫形成的復(fù)雜裂縫網(wǎng)絡(luò)在EDFM中的前處理。

(3)基于提出的前處理算法建立了能夠適用于任意形狀傾斜縫的三維EDFM，并給出一個(gè)三維油藏復(fù)雜縫網(wǎng)情況下衰竭開發(fā)的數(shù)值實(shí)例。

圖9 模擬計(jì)算的生產(chǎn)數(shù)據(jù)曲線

圖10 模擬計(jì)算的壓力分布