基于嵌入式離散裂縫和擴展有限元的裂縫性頁巖油藏流固耦合高效數(shù)值模擬方法

牛 駿, 蘇建政, 嚴 俠, 汪友平, 孫 海

(1.國家能源頁巖油研發(fā)中心， 中國石化石油勘探開發(fā)研究院，北京 100083；2.中國石油大學(華東)石油工程學院，青島 266580)

頁巖油氣作為一種不可或缺的非常規(guī)油氣資源，是目前研究和開發(fā)的熱點[1]。頁巖油氣藏儲層致密，其滲透率非常低，一般必須經(jīng)過大規(guī)模水力壓裂才可以進行商業(yè)化開采[2]。壓裂后的頁巖儲層，由于水力裂縫的存在使其在滲流機理和流動特征上均不同于常規(guī)儲層，表現(xiàn)出明顯的應力敏感特征，且具有較強的非均質各向異性，常規(guī)基于修正絕對滲透率的應力敏感模型已不適用，亟需建立滲流場和應力場耦合的流固耦合模型[3]。然而，相比于常規(guī)連續(xù)性儲層，頁巖油藏由于水力裂縫的存在，導致流動和力學性質發(fā)生較大變化，流固耦合效應更明顯，處理難度更大，常用的裂縫性儲層流固耦合數(shù)學模型有：等效連續(xù)介質模型[4-5]、雙重介質模型[6-7]和離散介質模型[8-9]。

等效連續(xù)介質模型將裂縫性介質看作一個連續(xù)體，首先需要判斷表征單元體的存在性，并確定其大小，再計算等效滲流和力學參數(shù)表征裂縫對流固兩場的影響，在此基礎上根據(jù)成熟的連續(xù)介質力學理論進行耦合流動分析，模型數(shù)據(jù)需求簡單而且計算效率高，一直以來在裂隙介質的水力學研究中占據(jù)重要地位。然而，由于等效連續(xù)介質模型過于簡化、等效參數(shù)計算不夠準確(尤其對于多相流問題)，無法準確刻畫裂縫性介質中的多種流動模式共存特征，因此，在油藏實際應用中其分析結果往往誤差較大；雙重介質模型最早由Barenblatt等[6]建立，Warren等[10]對其進行了簡化和推廣，此后眾多學者將多孔彈性力學與雙重介質模型相結合用于裂縫性介質流固耦合研究，這類模型將裂縫系統(tǒng)和基巖系統(tǒng)之間的物質交換考慮在內，可以實際地模擬裂縫內的優(yōu)先流現(xiàn)象。然而被裂縫分割成的基質巖塊系統(tǒng)被假定具有規(guī)則的大小和形狀，過于簡化卻不能充分表現(xiàn)出裂縫性介質的不連續(xù)性、各向異性等特點，同時基巖和裂縫間的竄流系數(shù)也難以準確獲取(尤其對于多相、非穩(wěn)態(tài)問題)，此外該模型僅適用于高度發(fā)育、空間尺度跨越小的裂縫性介質；離散介質模型對裂縫顯示表征，能夠反應裂縫對滲流場與應力場的影響。由于裂縫開度和油藏尺度相差極大，需要對裂縫進行降維處理，以防裂縫周圍局部加密而增加網(wǎng)格量。然而裂縫兩側的位移場具有間斷性，裂縫降維后采用常規(guī)數(shù)值方法難以有效解決，對此學者提出了相應的解決方法，主要分為以下三類：界面雙節(jié)點單元法[11]、位移不連續(xù)法[12]和擴展有限元法[13]。界面雙節(jié)點單元法在常規(guī)有限單元法的基礎上，在裂縫單元節(jié)點處新增一個額外節(jié)點，反映裂縫兩側變量的不連續(xù)性。該方法求解較簡單，但其網(wǎng)格劃分和編號存在一定難度，并且數(shù)值求解穩(wěn)定性和收斂性較差。位移不連續(xù)法作為一種邊界元方法，求解簡單、網(wǎng)格量小、精度較高，但其僅適用于均質儲層，對于實際強非均質各向異性儲層難以推導解析解建立離散格式，并且該方法最終形成的系數(shù)矩陣為滿秩矩陣，對于大規(guī)模問題的求解難度大。擴展有限元法在標準有限元框架下，在連續(xù)區(qū)域內采取標準有限元方法，在包含不連續(xù)邊界的區(qū)域內，增加間斷形函數(shù)以修正有限元的位移近似函數(shù)，采取水平集方法來描述間斷界面，使間斷的描述獨立于有限元網(wǎng)格，避免了以裂縫為約束的復雜非結構網(wǎng)格劃分，但擴展有限元作為一種改進的有限元方法，其方程離散為全局離散，不滿足單元的局部物質守恒，因此對于滲流方程(尤其多相流)，其計算精度較差。

近些年來，Yan等[14]基于等效連續(xù)介質模型和嵌套離散裂縫模型，提出了一種有效的數(shù)學模型來模擬含多尺度裂縫和溶洞的裂隙型多孔介質中流體力學耦合問題；張世明等[15]基于離散介質模型，提出了一種求解裂縫性介質等效滲透率的新方法；李友全等[16]基于離散介質模型，對致密油藏壓裂直井進行了動態(tài)分析；孫若凡等[17]基于雙重介質模型，對致密油藏壓裂水平井進行了模擬分析；盛茂等[18]基于雙重介質模型，建立了一種裂縫性頁巖氣藏的流固耦合模型并采用有限元法進行數(shù)值求解。然而目前并沒有基于離散介質模型并使用擴展有限元法求解的研究。于是，將擴展有限元法與嵌入式離散裂縫模型[19]相結合，建立一種適用于裂縫性頁巖油藏流固耦合數(shù)值模擬的新方法，該方法可以考慮裂縫開度的變化以及裂縫內流壓對開度變化的影響，分別采用擴展有限元和模擬有限差分[19-20]求解力學和滲流方程，模擬有限差方法適用于強非均質各向異性儲層，且具有良好的局部守恒性，可以準確求解壓力場和速度場。嵌入式離散裂縫模型網(wǎng)格劃分與擴展有限元類似，無需考慮儲層內的裂縫形態(tài)，大大降低了網(wǎng)格劃分難度，其計算精度與基于匹配性非結構網(wǎng)格的常規(guī)離散裂縫模型基本一致。

1 流固耦合數(shù)學模型

基于二維單相流固耦合問題來闡述本文數(shù)值模擬方法的原理與思想。假設在耦合流動過程中溫度恒定，考慮基巖和流體的壓縮性，基巖與裂縫中的流體流動均滿足達西定律，巖石變形為微小線彈性變形，不考慮慣性力。采用常規(guī)多孔介質彈性力學符號約定，即拉伸為正，壓縮為負。

1.1 巖石變形數(shù)學模型

(1)

σ=σeff-αpmI

(2)

(3)

(4)

σ·nf=-pf·nfon Γf

(5)

圖1 裂縫性介質及其內、外邊界示意圖Fig.1 The sketch of a fractured porous medium and its boundaries

1.2 嵌入式離散裂縫模型

基巖系統(tǒng)流動數(shù)學模型：

(6)

(7)

(8)

裂縫系統(tǒng)流動數(shù)學模型：

(9)

(10)

式中:vm和vf分別為基巖和裂縫內滲流速度，m/s；km和Kf分別為基巖和裂縫滲透率，m2；μ為流體黏度，mPa·s；ρl為流體密度，kg/m3；φ為基巖孔隙度；Ks和Kl分別為巖石顆粒和流體的體積模量，Pa；εv為基巖體積應變；qm和qf分別為基巖和裂縫的源匯項，s-1；Vm和Vf分別為基巖和裂縫的單元體積，m3；bf0和bf分別為裂縫初始和目前的開度，m；qmf為基巖與裂縫間的竄流量，m3/s；qff為相交裂縫間的竄流量，m3/s；δ為狄克拉函數(shù)；qmf、qff、δmf和δff的詳細計算格式和取值方法見參考文獻[19]；Γq為流量外邊界，Γp為定壓外邊界，如圖1所示。整個系統(tǒng)的初始壓力為p0。

2 模型離散求解

圖2 網(wǎng)格單元數(shù)值離散及自變量分布示意圖Fig.2 The schematic of discretization and locations of the solution variables

2.1 擴展有限元離散力學方程

擴展有限元根據(jù)單位分解法[21](PUM)的思想，在常規(guī)有限元連續(xù)位移場的基礎上，增加不連續(xù)位移場(裂縫段采用Heaviside函數(shù)，裂縫尖端采用尖端漸進函數(shù))來表征計算區(qū)域內的位移間斷特性，因此，其單元位移的近似解可以表示為

H[φ(xi)]}ai+

(11)

(12)

H(ξ)=1, ?ξ>0 andH(ξ)=-1, ?ξ<0

(13)

F(r,θ)=

(14)

將式(2)代入式(1)，考慮邊界條件[式(4)和式(5)]，得到方程(1)的積分弱形式為

(15)

式(15)中：δu為位移試函數(shù)；δu+和δu-分別表示裂縫正負兩側的位移試函數(shù)，如圖1所示。

裂縫壁面發(fā)生相對位移后，裂縫開度及滲透率(立方定律)表達式為

df=df0+(u+-u-)·nf

(16)

(17)

在擴展有限元中，采用水平集方法[22]確定裂縫的幾何信息，同時采用高斯數(shù)值積分方法計算單元和邊界上的積分，而對于裂縫穿過的單元，常規(guī)的高斯積分方法不適用，需要根據(jù)裂縫形狀對這些單元進行單元分解，再在每個分解子單元內進行數(shù)值積分。

2.2 模擬有限差分離散滲流方程

(18)

對于方程(7)，在基質網(wǎng)格單元上積分且應用散度定理得

(19)

考慮裂縫滲透率為標量且為一維流動，因此采用有限差分進行離散。將方程(9)代入方程(10)，方程左右兩端同乘裂縫網(wǎng)格單元體積得

(20)

2.3 整體計算格式

以存在一條裂縫為例，考慮基巖網(wǎng)格單元邊界面上的速度連續(xù)，將力學方程和滲流方程的數(shù)值離散格式組裝到一起，得到嵌入式離散裂縫流固耦合模型數(shù)值計算格式為

(21)

(22)

(23)

式中:Ku為擴展有限元位移剛度矩陣；Qum和Quf分別為基巖和裂縫的流固耦合系數(shù)矩陣；Km為基巖模擬有限差分系數(shù)矩陣；Tmf為基巖與裂縫竄流系數(shù)矩陣；Kf為嵌入式裂縫有限差分系數(shù)矩陣；Sm和Sf分別為基巖和裂縫的壓縮項系數(shù)矩陣；u、a和b分別為節(jié)點位移和額外未知量，整體對應于系數(shù)矩陣A的第一行；vm、pm和πm分別為網(wǎng)格邊界流速、網(wǎng)格中心壓力和網(wǎng)格邊界壓力，整體對應于系數(shù)矩陣A的第二行；pf為裂縫單元壓力；C為體積力載荷、邊界應力載荷以及基巖和裂縫源匯項列陣。

對于方程(21)中的時間導數(shù)項，采用全隱式差分格式進行離散。

(Adt+B)Xn+1=BXn+Cdt

(24)

式(24)中：n表示第n個時間步；dt為時間步長。

3 算例分析

首先分別通過與一維多孔介質和二維裂縫性介質流固耦合問題的解析解和精細網(wǎng)格數(shù)值解進行對比，驗證所提方法的正確性；然后分析不同條件下裂縫開度的變化情況以及裂縫開度變化對耦合流動的影響；最后對于某一二維頁巖油藏彈性開采進行數(shù)值模擬，對比基巖滲透率為對角張量和全張量時的差異，說明所提方法對于全張量滲透率的適用性。本文算例均采用平面應變假設。

3.1 模型驗證

考慮如圖3(a)所示的深度200 m的一維流固耦合模型，初始流壓和位移均為0，在上、下邊界同時施加100 kN/m2的恒定邊界載荷，左、右邊界封閉，上、下為邊界定壓(壓力為0)，模型參數(shù)如表1所示。采用所提方法對該一維模型進行數(shù)值求解并與其解析解進行對比，圖3(b)給出了不同時刻流體壓力的數(shù)值解與解析解對比結果(由于模型上、下對稱，圖中只給出了上半部分的壓力對比)，能夠看出數(shù)值解與解析解基本一致。

表1 模型基本參數(shù)Table 1 The model basic parameters

圖3 一維模型幾何示意圖及不同時間流體壓力數(shù)值解與解析解對比Fig.3 The sketch of a 1D model and the comparison of numerical solution and analytical solution of fluid pressure at different times

如圖4(a)所示，為一10 m×20 m 的裂縫性介質模型，其中裂縫長10 m、開度0. 001 m，垂直位于模型區(qū)域中心。初始流壓和位移均為0，模型的力學邊界條件如圖4(a)所示，滲流邊界條件為：上邊界定壓(壓力為0)，其余邊界封閉。裂縫孔隙為1.0，初始滲透率為8.33×10-8m2，其他參數(shù)如表1所示。由于裂縫的存在，該模型難以推導解析解，因此，將裂縫顯示處理，通過在裂縫周圍進行局部網(wǎng)格加密構建精細網(wǎng)格，如圖4(b)所示，采用有限單元法求解作為參考解。圖4(c)為所提方法采用的計算網(wǎng)格，可以看出其網(wǎng)格劃分難度及網(wǎng)格數(shù)量都大大降低。

不考慮裂縫開度變化，分別采用精細網(wǎng)格有限元方法和所提方法對該裂縫性介質流固耦合模型進行模擬，圖4(d)和圖4(e)分別為模擬10 000 s后兩種方法得到的壓力場圖，圖5(a)和圖5(b)分別給出了兩種方法計算得到的測試點A和B處壓力和y方向位移隨時間的變化曲線對比。通過對比可以看出，所提方法的計算結果與精細網(wǎng)格有限元方法參考解基本一致。因此，本節(jié)的兩個算例驗證了所提方法的正確性。

圖4 裂縫性介質模型示意圖、不同方法網(wǎng)格劃分結果及模擬10 000 s后壓力場圖Fig.4 The model and grids of a fractured medium, pressure distribution after 10 000 seconds for different methods

圖5 兩種方法計算得到的測試點壓力和位移隨時間變化曲線Fig.5 Pressure and displacement (y) the test point various time obtained from two methods

3.2 裂縫開度隨時間變化

由于本文模型分別采用嵌入式離散裂縫和擴展有限元中將裂縫的滲流和力學特征顯示表征，因此可以有效刻畫裂縫開度隨時間的變化情況以及裂縫開度變化對耦合流動的影響?；趫D4(a)所示的裂縫多孔介質模型，模型參數(shù)與上一個算例一致，不考慮裂縫內流體壓力的影響，分別采用所提方法和精細網(wǎng)格有限元方法計算100、2 000、10 000 s時的裂縫開度分布，結果如圖6(a)所示，能夠看出兩者結果基本一致，驗證了所提方法用于計算裂縫開度的正確性。圖6(b)給出了不考慮與考慮裂縫內流體壓力影響兩種情況下所提方法計算得到的不同時間裂縫開度分布，可以看出模擬100 s時，兩者裂縫開度相差很大，隨著模擬時間增長，兩者之間的差距逐漸減小，到模擬10 000 s時，兩者裂縫開度基本一致。這是由于早期裂縫內流體壓力較高，考慮其影響后對裂縫起到支撐作用，可以有效阻止裂縫閉合，隨著模擬時間增長，裂縫內的流體壓力逐漸減小，其支撐作用逐漸減弱。

圖6 不同時間裂縫開度分布Fig.6 Fracture aperture distribution at different time

圖7 模擬10 000 s后的壓力場分布圖Fig.7 Pressure distribution after 10 000 seconds

在以上模型的基礎上，將上邊界載荷增加到200 kN/m2，其他模型參數(shù)和條件均不變。采用所提方法分別對裂縫開度不變和變化兩種情況下的流固耦合進行模擬，圖7(a)和7(b)分別給出了模擬10 000 s后兩種情況下壓力場分布，圖8給出了兩種情況下A點處壓力及兩者誤差隨時間的變化曲線，可以看出：考慮裂縫開度變化后，隨著模擬時間增加，裂縫開度和導流能力逐漸減小，流體流出速度減小，區(qū)域內的流體壓力增加，兩者相對誤差逐漸增大。

圖8 兩種情況下A點處壓力及兩者相對誤差隨時間變化曲線Fig.8 Pressure and relative error at point A for various time in the two cases

3.3 裂縫性頁巖油藏彈性開采

圖9 裂縫性頁巖油藏模型示意圖及網(wǎng)格劃分結果Fig.9 The sketch of fractured shale reservoir model and meshing results

表2 頁巖油藏基本參數(shù)Table 2 The reservoir parameters

圖10～圖12分別給出了基巖滲透率為對角張量和全張量兩種情況下，模擬10 d后的壓力場和x方向位移場圖、不同裂縫的開度變化情況(裂縫1為最左側裂縫，裂縫3為從左側數(shù)第三條裂縫)以及累產(chǎn)隨時間變化情況。可以看出，考慮全張量滲透率后，壓力場、位移場、裂縫開度和累積產(chǎn)量相比于對角張量滲透率條件下的結果有一定變化，說明對于基巖滲透率為全張量的實際強各向異性儲層，不能采用有限差分或有限體積等忽略非主對角線值的常規(guī)數(shù)值方法。

圖10 對角張量滲透率下彈性開采10 d時壓力場和x方向位移Fig.10 Pressure and displacement distribution after 10 days for diagonal-tensorpermeability

圖11 全張量滲透率下彈性開采10 d時壓力場和x方向位移Fig.11 Pressure and displacement (x) distribution after 10 days for full-tensorpermeability

圖12 不同基巖滲透率條件下，模擬10 d裂縫開度和累產(chǎn)對比Fig.12 Fracture apertures and cumulative production for different matrix permeability on different bedrock permeability conditions after 10 days

4 結論

通過耦合嵌入式離散裂縫模型和擴展有限元提出一種適用于裂縫性頁巖油藏的流固耦合高效數(shù)值模擬方法，該方法的創(chuàng)新性以及所取得的結論如下。

(1)該方法對裂縫使用降維顯式處理，能夠準確刻畫裂縫對滲流場和應力場的互相影響，同時可以模擬裂縫開度的動態(tài)變化過程以及裂縫內流體壓力對裂縫變形的影響。該模型劃分網(wǎng)格時不必考慮裂縫幾何形態(tài)，只需對基巖系統(tǒng)進行正交網(wǎng)格剖分，使得網(wǎng)格劃分時的復雜度大幅降低。

(2)對于模型中流動方程求解采取模擬有限差分方法，滿足局部質量守恒并且可以有效處理全張量形式滲透率；力學方程求解采取擴展有限元方法，可以準確處理裂縫兩側的位移間斷性。

(3)分別通過與一維多孔介質和二維裂縫性介質流固耦合問題的解析解和精細網(wǎng)格數(shù)值解進行對比，驗證所提方法的正確性，分析了不同條件下裂縫開度的變化情況以及裂縫開度變化對耦合流動的影響。

(4)針對一、二維裂縫性頁巖油藏彈性開采進行了數(shù)值模擬，對比基巖滲透率為對角張量和全張量時的差異，說明所提方法對于全張量滲透率的適用性。