陳小凡 唐 潮 杜志敏 湯連東 魏嘉寶 馬 旭
1.“油氣藏地質(zhì)及開發(fā)工程”國家重點(diǎn)實(shí)驗(yàn)室·西南石油大學(xué) 2.中國石化勝利油田分公司河口采油廠3.中國石油長慶油田分公司第五采氣廠
合理考慮不同類型介質(zhì)中頁巖氣的流動(dòng)機(jī)理對(duì)于頁巖氣數(shù)值模擬具有重要意義[1]?,F(xiàn)有的頁巖氣數(shù)值模擬模型主要包括雙重介質(zhì)、多重離散介質(zhì)和等效介質(zhì)模型,其中雙重介質(zhì)模型使用較為廣泛。Kucuk和Sawyer[2]首次研究了基于雙孔模型的頁巖氣儲(chǔ)層壓力變化。隨后,Bumb和McKee[3]通過對(duì)Langmuir等溫吸附方程添加額外的吸附系數(shù)來研究吸附—解吸行為對(duì)不穩(wěn)定流動(dòng)瞬態(tài)行為的影響。然而,以上研究工作都忽略了納米及微米尺度的擴(kuò)散過程。Carlson和Mercer[4]通過在雙孔模型中引入解吸擴(kuò)散項(xiàng)研究頁巖氣直井壓力的變化,由于該模型沒有考慮滑脫效應(yīng),在生產(chǎn)初期對(duì)頁巖氣產(chǎn)能預(yù)測結(jié)果偏小。Swami等[5]建立了一個(gè)考慮Knudsen擴(kuò)散、滑脫以及吸附—解吸過程的雙重介質(zhì)模型,并通過實(shí)驗(yàn)室數(shù)據(jù)進(jìn)行了驗(yàn)證。
雖然雙重介質(zhì)模型在各類商業(yè)軟件中廣泛使用,由于其固有缺陷,即使進(jìn)行全網(wǎng)格或局部加密,仍然對(duì)多尺度以及大面積分布的縫網(wǎng)系統(tǒng)適應(yīng)性較差。首先,傳統(tǒng)雙重介質(zhì)模型均沒有考慮頁巖中多尺度流動(dòng)機(jī)理,不能真實(shí)反映頁巖氣井生產(chǎn)實(shí)際[6-8]。其次,由于頁巖儲(chǔ)層中存在復(fù)雜的多尺度縫網(wǎng),造成數(shù)值模擬工作中滲透率場非均質(zhì)性嚴(yán)重,如果要對(duì)局部裂縫區(qū)域進(jìn)行網(wǎng)格剖分,勢必產(chǎn)生大量的極小化網(wǎng)格,導(dǎo)致計(jì)算收斂性差甚至無法收斂[9-10]。同時(shí),基質(zhì)—裂縫界面處流度的劇烈變化也產(chǎn)生了階躍現(xiàn)象,造成的非物理振蕩會(huì)對(duì)計(jì)算精度產(chǎn)生影響[11-12]。因此,雙重介質(zhì)模型應(yīng)用在頁巖氣模擬中最大的限制在于難以選擇合適的特征體積單元,由于裂縫的多尺度性,裂縫性儲(chǔ)層是否存在特征體積單元或特征體積單元大小的確定仍存在較大爭議。因此,部分學(xué)者提出采用多重介質(zhì)模型來對(duì)頁巖氣井產(chǎn)能進(jìn)行模擬研究[13-14]?;诙嘀亟橘|(zhì)的概念,Schepers[15]、Dehghanpour等[16]分別建立了將吸附—解吸過程與基質(zhì)中流動(dòng)耦合的達(dá)西流動(dòng)模型。Wu等[17]和張烈輝等[18]將裂縫劃分為天然微裂縫和人工壓裂縫,建立了考慮裂縫應(yīng)力敏感和滑脫效應(yīng)的致密裂縫性油藏多重介質(zhì)模型,氣體在基質(zhì)中的流動(dòng)考慮滑脫效應(yīng),在裂縫中的流動(dòng)為高速非達(dá)西流動(dòng),并比較了多重介質(zhì)模型與雙重介質(zhì)模型的差別。Aboaba和Cheng[19]利用線性流模型描述了頁巖氣儲(chǔ)層中壓裂水平井的典型產(chǎn)能曲線及產(chǎn)量變化規(guī)律,但沒有考慮吸附和擴(kuò)散的影響。Wang等[20]采用嵌入式離散裂縫描述大尺度裂縫、連續(xù)模型描述中小尺度裂縫的方法建立了一種模擬頁巖多尺度裂縫中的流體流動(dòng)模型。方文超等[21]考慮致密儲(chǔ)集層的可壓縮性和基質(zhì)流體的非線性滲流,建立了基于二維體積壓裂的多尺度滲流裂縫模型。
在前人研究的基礎(chǔ)上,筆者將基質(zhì)—天然微裂縫視為雙重介質(zhì)、人工壓裂縫視為離散介質(zhì),考慮基質(zhì)孔隙中頁巖氣的吸附—解吸效應(yīng),建立了描述頁巖氣復(fù)雜流動(dòng)過程的數(shù)學(xué)模型,基于有限體積方法離散得到頁巖氣三維滲流數(shù)值模型,利用降維方法將壓裂縫處理成二維平面實(shí)體,井筒處理為一維線元實(shí)體,再利用順序求解方法進(jìn)行求解,進(jìn)而模擬頁巖氣多段壓裂水平井生產(chǎn)動(dòng)態(tài)和儲(chǔ)層壓力分布變化。所取得的研究成果可用于頁巖氣儲(chǔ)層體積壓裂設(shè)計(jì)以及多段壓裂水平井的生產(chǎn)動(dòng)態(tài)預(yù)測。
圖1為頁巖儲(chǔ)層中一口多段壓裂水平井網(wǎng)格模型示意圖,整個(gè)模擬區(qū)域沿x軸方向長度為1 200 m、y軸方向?qū)挾葹?00 m、z軸方向高度為100 m。水平井段長度為700 m,多段壓裂水平井空間配置如圖1-a所示,人工壓裂縫主縫半長100 m,縫高60 m;Ⅰ級(jí)分支裂縫縫長70 m,縫高60 ;Ⅱ級(jí)分支裂縫縫長50 m,縫高60 m。人工壓裂縫用二維面元實(shí)體表示,水平井段用一維線元實(shí)體表示。為簡化模型,作出如下假設(shè):①氣藏為矩形封閉頁巖氣藏,流動(dòng)為單相氣體等溫滲流,忽略重力、氣體滑脫、Knudsen擴(kuò)散及應(yīng)力敏感影響;②水平井位于氣藏中心位置,人工壓裂主裂縫垂直于井筒且關(guān)于井筒對(duì)稱,氣體只能通過人工壓裂縫流入井筒;③氣體在人工壓裂縫、天然微裂縫中的滲流服從達(dá)西定律,基質(zhì)中頁巖氣吸附—解吸過程采用Langmuir等溫吸附方程描述;④氣井定壓生產(chǎn),忽略水平井水平段壓降。
1.2.1 控制方程推導(dǎo)
根據(jù)真實(shí)氣體狀態(tài)方程,頁巖氣密度計(jì)算式為:
式中ρgi(i=m、f、hf)表示基質(zhì)、天然微裂縫和人工壓裂縫中氣體密度,kg/m3;Mg表示氣體摩爾質(zhì)量,g/mol;pi(i=m、f、hf)表示基質(zhì)、天然微裂縫和人工壓裂縫壓力,MPa;Z表示氣體偏差因子,無量綱;
圖1 網(wǎng)格模型示意圖
注:點(diǎn)C、F為兩個(gè)相鄰控制體的形心;f為界面(綠色背景);b為界面f的邊界(紅色線);dhf為人工裂縫無因次開度;dCf為控制體形心指向界面形心的位置向量,m;Sf為界面f上指向外側(cè)的法向向量,m2R表示氣體普適常數(shù),8 314 J/(kmol·K);T表示氣藏溫度,K。
依據(jù)物質(zhì)守恒定律,基質(zhì)中頁巖氣連續(xù)性方程為:
式中φm表示基質(zhì)孔隙度,無因次;t表示時(shí)間,s;Km表示基質(zhì)滲透率,mD;μg表示氣體黏度,mPa·s;qm表示基質(zhì)解吸氣率,kg/(m3·s);qmf表示基質(zhì)向微裂縫竄流率,kg/(m3·s)。
式(2)左邊第1項(xiàng)為單位體積基質(zhì)微元內(nèi)流體質(zhì)量變化,第2項(xiàng)為通過微元表面流出的流體質(zhì)量通量,第3項(xiàng)為基質(zhì)解吸氣量,第4項(xiàng)為基質(zhì)向裂縫竄流量。在此假設(shè)微孔隙內(nèi)為單相氣,由頁巖基質(zhì)表面的吸附氣和孔隙內(nèi)游離氣組成,根據(jù)Langmuir等溫吸附方程得到基質(zhì)中解吸氣量計(jì)算式為:
式中Vm表示單位體積基質(zhì)吸附氣量,m3;Vstd表示標(biāo)準(zhǔn)狀況下氣體摩爾體積,m3/mol;VL表示Langmuir體積,m3/kg;pL表示Langmuir壓力,MPa。
而頁巖氣在微裂縫中的連續(xù)性方程為:式中φf表示天然微裂縫孔隙度,無因次;Kf表示天然微裂縫滲透率,mD;qhf表示微裂縫向人工壓裂縫竄流率,kg/(m3·s)。
同理,在水力壓裂縫中的連續(xù)性方程為:
式中φhf表示人工裂縫孔隙度,無因次;Khf表示人工壓裂縫滲透率,mD;qwell表示水平井產(chǎn)氣率,kg/(m3·s)。
1.2.2 初始條件
式中pinit表示原始地層壓力,MPa。
1.2.3 邊界條件
令Γout表示求解域外邊界,Γin表示內(nèi)邊界。假設(shè)模型外邊界為封閉邊界,生產(chǎn)井定壓生產(chǎn),則邊界條件為:
由于四面體網(wǎng)格幾何中心與形心重合,為簡化計(jì)算步驟,本文以四面體網(wǎng)格和Delaunay三角形網(wǎng)格分別對(duì)基質(zhì)—微裂縫和人工壓裂縫進(jìn)行網(wǎng)格剖分。其他類型網(wǎng)格需要重新計(jì)算控制體形心,在此不再詳述。首先將基質(zhì)—微裂縫系統(tǒng)用四面體網(wǎng)格表示,考慮為連續(xù)介質(zhì)(圖1-b);人工壓裂縫降維處理成二維平面,由四面體網(wǎng)格間Delaunay三角形界面表示(圖1-c),對(duì)裂縫降維是提高多尺度模擬計(jì)算收斂性的關(guān)鍵處理方法,若在網(wǎng)格系統(tǒng)內(nèi)將裂縫考慮為三維,必然會(huì)在開度很小的裂縫空間內(nèi)進(jìn)行四面體網(wǎng)格剖分,進(jìn)而形成大量的極小化網(wǎng)格,導(dǎo)致后續(xù)求解過程無法繼續(xù)[22]。Juanes等[23]研究表明,將裂縫考慮成二維比將其考慮成三維的計(jì)算收斂性顯著提高??刂企w單元由網(wǎng)格剖分所得到的四面體確定(圖1-c),該方法可保證控制體單元互不重疊地覆蓋整個(gè)研究區(qū)域。然后,將整個(gè)求解區(qū)域Γd劃分為兩個(gè)子域Γm-f和Γhf,分別代表基質(zhì)—微裂縫區(qū)域和人工壓裂縫區(qū)域,根據(jù)連續(xù)介質(zhì)場理論,整個(gè)模擬區(qū)域流動(dòng)控制方程(Flow Governing Equation,簡稱為FGE)的積分形式為:
本節(jié)以基質(zhì)系統(tǒng)流動(dòng)方程為例說明在四面體網(wǎng)格中的求解方法。微裂縫與人工壓裂縫求解方法與之類似。首先將式(1)、(3)代入式(2),并在控制體VC范圍內(nèi)進(jìn)行積分,即
式中VC表示以C為形心的控制體;dV表示控制體微元體積,m3。
式(9)說明,對(duì)于任意控制體,在某一時(shí)間段內(nèi)基質(zhì)系統(tǒng)氣體質(zhì)量流量的變化等于流入流出基質(zhì)氣量與基質(zhì)解吸氣量以及基質(zhì)—微裂縫竄流量之和,因此保證了該模型在任意局部網(wǎng)格仍然遵守質(zhì)量守恒定律。
式(9)中不穩(wěn)定流動(dòng)項(xiàng)可以近似表達(dá)為:
根據(jù)散度定理,將體積積分轉(zhuǎn)化為面積分,則式(9)中對(duì)流項(xiàng)可以寫為:
進(jìn)一步對(duì)控制體VC的面積分寫作對(duì)控制體各界面上通量面積分之和的形式為:
式中SC表示構(gòu)成控制體VC的所有界面,以圖1-c為例,控制體VC由4個(gè)三角形界面構(gòu)成。
進(jìn)一步由梯形積分公式可知,在界面f上的面積分可以寫為:
式中Sf表示界面f上指向外側(cè)的法向向量,有
因此,式(9)中對(duì)流項(xiàng)可以寫為:
式中λm表示基質(zhì)中氣體流度,m2/(mPa·s)。同理,得到源匯項(xiàng)有限體積計(jì)算式為:
式(15)等號(hào)右邊第1項(xiàng)為控制體內(nèi)解吸氣量,第2項(xiàng)為基質(zhì)—微裂縫竄流量。由于解吸氣量是時(shí)間的函數(shù),則根據(jù)式(3)可以將基質(zhì)系統(tǒng)流動(dòng)方程的有限體積離散格式改寫為:
同理可以得到微裂縫及人工壓裂縫中流動(dòng)方程離散格式為:
式中λf和λhf表示天然微裂縫和人工壓裂縫中氣體流度,m2/(mPa·s); SCF表示控制體VC與大尺度裂縫共用的界面;b表示構(gòu)成界面f的所有邊界;lb表示垂直邊界b的外法向向量,其中 l是邊界b的長度(二維界面厚度為1)。
所謂順序求解方法,就是在每個(gè)時(shí)間步上,首先求解出某個(gè)變量,再代入其他變量表達(dá)式進(jìn)行迭代求解。此方法保證了在每個(gè)時(shí)間步上計(jì)算量小于整體求解方法。假設(shè)當(dāng)前時(shí)間步為k,則所有與人工壓裂縫壓力及基質(zhì)壓力有關(guān)的變量都采用(k+1)時(shí)間步的值進(jìn)行隱式求解,則式(16)、(18)可以寫為:
注意到式(19)、(20)中pf使用的是k時(shí)間步的值,是一個(gè)已知數(shù)。因此兩式中各含有一個(gè)未知變量phf和pm。因此可以使用Newton-Raphson方法迭代求解。將求解得到的帶入微裂縫流動(dòng)方程對(duì)裂縫壓力進(jìn)行顯式求解,得
顯然,為求解式(19)、(20)、(21),需要計(jì)算控制體內(nèi)壓力梯度的變化,在此給出一種基于格林高斯定理的方法,該方法相對(duì)簡單且適合各種幾何結(jié)構(gòu)的網(wǎng)格(結(jié)構(gòu)/非結(jié)構(gòu),正交/非正交)[24-25]。體積為VC、形心為C的控制體內(nèi)平均壓力梯度計(jì)算式為:
根據(jù)散度定理將體積積分轉(zhuǎn)化為面積分,有
對(duì)于離散界面,式(23)可以寫為:
然后,將沿控制體表面的面積分用積分中值定理近似表示為表面形心處的插值乘以表面向量的形式,即
由式(25)可知,為計(jì)算控制體內(nèi)壓力梯度,首先需要知道控制體表面向量Sf,并通過控制體形心處壓力值pC和相鄰控制體形心處壓力值pF插值得到界面f上的壓力,其插值格式為:
其中g(shù)C和gF為考慮控制體體積VC和VF的權(quán)重值,計(jì)算公式為:
考慮到相鄰的兩個(gè)控制體(網(wǎng)格編號(hào)分別為j和k),表面向量Sf不能同時(shí)向外,因此對(duì)于特定網(wǎng)格定義表面向量方向由控制體網(wǎng)格編號(hào)確定,始終由編號(hào)較小的控制體指向編號(hào)較大的控制體,則式(25)改寫為:
由于本文中采用的非結(jié)構(gòu)網(wǎng)格,其網(wǎng)格形態(tài)往往是非正交的。即控制體表面向量與連接控制體形心的向量不共線,也就是說控制體表面法向向量與流體速度(壓力梯度)方向不一致。這種情況下,因?yàn)榇嬖谝粋€(gè)垂直于向量的分量,控制體的壓力梯度不能寫成f(pC, pF)形式[26-29]。
定義e為沿著形心C、F連線所確定方向的單位矢量,m,則沿著e方向的壓力梯度可以表示為:
式中rC、rF表示形心C和F的位置向量,m;dCF表示相鄰控制體形心C與F間的距離,m。
因此,為了實(shí)現(xiàn)非正交網(wǎng)格中通量的線性化,表面向量Sf應(yīng)該寫作兩個(gè)向量Ef、Tf之和,即
其中向量Ef沿著形心C、F連線方向,m2;Tf為Sf的非正交分量,m2,則通過界面f的通量可以表示為:
對(duì)于向量Ef、Tf的計(jì)算通常有以下3種方法(表1)。
上述3種分解方法所得到的結(jié)果均滿足式(7)中擴(kuò)散通量計(jì)算需要,其主要差別主要體現(xiàn)在非正交網(wǎng)格計(jì)算過程中的準(zhǔn)確性和穩(wěn)定性。根據(jù)前人研究發(fā)現(xiàn),即使在高度非正交的網(wǎng)格中超松弛修正方法仍能保證通量計(jì)算的穩(wěn)定性[24-25,30-31]。
表1 非正交網(wǎng)格中控制體表面向量分解方法
2.6.1 模型對(duì)比驗(yàn)證
為了驗(yàn)證有限體積方法的正確性,建立矩形封閉氣藏模型,在基質(zhì)和微裂縫中僅考慮黏性流的影響,忽略吸附—解吸機(jī)制,并認(rèn)為人工壓裂縫為無限導(dǎo)流縫,模型參數(shù)如表2所示。將模擬結(jié)果與油氣藏商業(yè)數(shù)值模擬軟件Eclipse(以下簡稱軟件Eclipse)計(jì)算的結(jié)果進(jìn)行對(duì)比,如圖2所示,本文模型與軟件Eclipse計(jì)算的水平井產(chǎn)氣量基本一致,說明本文采用的數(shù)值計(jì)算方法正確、可行。
如圖2所示,在早期段本文數(shù)值解與商業(yè)軟件計(jì)算結(jié)果存在一定誤差,原因在于壓裂縫附近網(wǎng)格精度不夠?qū)е拢ㄗ畲髥卧介L15 m,最小單元步長3 m,最大單元遞增率1.2,單元曲率因子0.5,狹窄區(qū)域分辨率0.8),可以通過對(duì)壓裂縫位置處的網(wǎng)格進(jìn)行加密處理以提高計(jì)算精度。
2.6.2 求解方法對(duì)比
基于前述模型,通過不同求解方法(順序求解方法和全隱式求解方法)進(jìn)行計(jì)算,結(jié)果顯示差距主要在生產(chǎn)初期,順序求解方法計(jì)算得到的開井瞬時(shí)產(chǎn)氣量為195×104m3/d,而全隱式求解方法計(jì)算得到的開井瞬時(shí)產(chǎn)氣量為167×104m3/d,但隨著計(jì)算推進(jìn),兩種方法計(jì)算的產(chǎn)氣量迅速趨于一致(圖3),進(jìn)一步驗(yàn)證了模型的正確性。
表2 多段壓裂水平井?dāng)?shù)值模擬參數(shù)表[32]
圖2 本文數(shù)值解與軟件Eclipse計(jì)算結(jié)果對(duì)比圖
根據(jù)本文參考文獻(xiàn)[1],建立頁巖氣多段壓裂水平井模型,原始地層壓力30 MPa,地層溫度為343.15 K,基質(zhì)孔隙度為8%,基質(zhì)滲透率為0.003 mD,天然微裂縫孔隙度為0.3%,天然微裂縫滲透 率 為 0.5 mD,Langmuir體 積 為 4×10-3m3/kg,Langmuir壓力為5 MPa,頁巖密度為2 600 kg/m3,甲烷分子摩爾質(zhì)量為16 g/mol,標(biāo)況下頁巖氣摩爾體積為 0.022 4 m3/mol,氣體黏度為 1.85×10-2mPa·s,井底流壓為5 MPa。通過軟件Eclipse和CMG計(jì)算發(fā)現(xiàn),在上述模型壓力變化過程中氣體偏差因子變化范圍介于0.908 9~0.982 8,利用多項(xiàng)式擬合將偏差因子改寫成壓力函數(shù)的形式。多段壓裂水平井空間配置如圖1-a所示,人工壓裂縫主縫導(dǎo)流能力為15 D·cm,Ⅰ級(jí)分支裂縫導(dǎo)流能力為10 D·cm,Ⅱ級(jí)分支裂縫導(dǎo)流能力為7 D·cm。
3.2.1 Langmuir體積的影響
基于物理模型參數(shù),模擬頁巖氣多段壓裂水平井生產(chǎn)情況,如圖4所示,在前5年儲(chǔ)層壓力下降區(qū)域主要集中在體積改造區(qū)域,而未壓裂區(qū)域內(nèi)儲(chǔ)層壓力變化幅度較小,表明該階段產(chǎn)出的氣體主要來自體積改造區(qū)域內(nèi)的游離氣及解吸氣。
圖5為不同Langmuir體積所對(duì)應(yīng)的平均地層壓力、解吸氣量、產(chǎn)氣量以及解吸氣貢獻(xiàn)比變化曲線??梢钥闯觯馕鼩馐沟貙訅毫ι仙?,但作用有限,對(duì)產(chǎn)氣量影響不大;隨著生產(chǎn)時(shí)間延長,解吸氣量絕對(duì)值逐漸減少,但是在產(chǎn)氣量中所占的比例逐漸上升。
圖3 不同求解方法計(jì)算結(jié)果對(duì)比圖
圖4 不同生產(chǎn)時(shí)間地層壓力分布圖
圖5 不同Langmuir體積下水平井生產(chǎn)動(dòng)態(tài)預(yù)測曲線圖
3.2.2 裂縫形態(tài)的影響
通過改變壓裂段數(shù)進(jìn)一步研究裂縫形態(tài)對(duì)頁巖氣水平井產(chǎn)氣量的影響,模擬結(jié)果顯示,生產(chǎn)5年時(shí),3段壓裂模式下的壓降區(qū)域范圍與4段壓裂模式下的壓降區(qū)域范圍較接近(圖6);4段壓裂模式下投產(chǎn)初期產(chǎn)氣量(64×104m3/d)為3段壓裂模式下初期產(chǎn)氣量(47×104m3/d)的1.37倍,但4段壓裂模式下產(chǎn)氣量遞減速度更快,5年后兩種壓裂模式下的水平井產(chǎn)氣量基本一致,且該階段的累產(chǎn)氣量曲線平行(圖7)。由此可見,壓裂段數(shù)繼續(xù)增加對(duì)于增加累產(chǎn)氣量效果已不顯著,需要通過增加裂縫半長來擴(kuò)大儲(chǔ)層改造體積。確定合理的壓裂段數(shù),同時(shí)獲得較長的壓裂縫長,是頁巖氣水平井增產(chǎn)改造措施的核心內(nèi)容。
圖6 不同壓裂段數(shù)水平井生產(chǎn)5年時(shí)地層壓力分布圖
圖7 不同壓裂段數(shù)下頁巖氣水平井生產(chǎn)曲線圖
1)采用所建立的模型與軟件Eclipse計(jì)算的多段壓裂水平井產(chǎn)氣量基本一致,說明該數(shù)值計(jì)算方法正確、可行;通過順序求解方法和全隱式求解方法的對(duì)比,進(jìn)一步驗(yàn)證模型的正確性。
2)解吸氣對(duì)地層壓力有補(bǔ)充作用,但作用有限,對(duì)產(chǎn)氣量影響不大,隨著生產(chǎn)時(shí)間的延長解吸氣量在產(chǎn)氣量中所占比例逐漸上升。
3)確定合理的壓裂段數(shù),同時(shí)獲得較長的壓裂縫長,是頁巖氣水平井增產(chǎn)改造措施的核心內(nèi)容。