張東旭，張烈輝，唐慧瑩，趙玉龍

（西南石油大學油氣藏地質(zhì)及開發(fā)工程國家重點實驗室，成都 610500）

0 引言

致密油資源豐富、開發(fā)潛力大，在中國油氣資源中占據(jù)十分重要的地位[1-4]。與常規(guī)油氣儲集層相比，致密油儲集層表現(xiàn)為多尺度性、強非均質(zhì)性，儲集層有效應力狀態(tài)改變時，具有極強的流-固耦合效應[5-6]。傳統(tǒng)滲流理論一般將多孔介質(zhì)骨架視作剛性模型，但致密油儲集層多孔介質(zhì)變形與孔隙流體流動的耦合作用不可忽略，建立在傳統(tǒng)滲流理論基礎上的油藏數(shù)值模擬技術難以準確預測致密油儲集層開發(fā)動態(tài)[7-8]。

Terzaghi有效應力原理和 Biot三維固結理論共同奠定了流-固耦合研究的基礎[9]，許多學者在此基礎上建立了多相飽和流體滲流和多孔介質(zhì)變形耦合作用模型[10-12]。致密油儲集層作為一種典型的雙重介質(zhì)，包含基質(zhì)孔隙系統(tǒng)和天然裂縫系統(tǒng)，且二者通過竄流的方式進行流體交換，具有孔喉狹小、滲透率低的特點[5,13]。通常采用水平井水力壓裂的方式改造致密油儲集層，除原有小尺度的天然裂縫外，還出現(xiàn)了大尺度的導流裂縫[14-16]。隨著致密油的采出，巖石骨架有效應力重新分布，骨架發(fā)生變形，儲集層物性發(fā)生動態(tài)變化，這些變化又會影響流體在儲集空間的流動和壓力分布[17]。因此，致密油開發(fā)過程屬于多尺度、多物理過程，且存在強非線性的力學問題，如何建立考慮多尺度流動、多物理場耦合的模型并開展壓裂井生產(chǎn)動態(tài)預測，成為致密油研究的難點。本文提出一種有效的數(shù)學模型來模擬致密油儲集層基質(zhì)孔隙、天然裂縫及人工裂縫中的流-固耦合，在一定程度上完善了致密油儲集層多級壓裂水平井多尺度多場耦合滲流理論。

油藏數(shù)值模擬技術是預測壓裂水平井生產(chǎn)動態(tài)及其增產(chǎn)效果的有效方法，目前采用的數(shù)值模型主要包括連續(xù)介質(zhì)模型、離散裂縫模型和連續(xù)介質(zhì)-離散裂縫混合模型，其中連續(xù)介質(zhì)-離散裂縫混合模型具有雙重介質(zhì)和離散裂縫模型的優(yōu)點，在裂縫儲集層模擬中具有廣闊的應用前景[18-19]。流-固耦合問題的系統(tǒng)方程本身具有不可獨立求解和無法顯式消去兩類獨立變量（壓力、位移）的特征，根據(jù)求解方法的不同，可分為間接耦合方法（顯式迭代耦合、隱式迭代耦合）和直接求解的全耦合方法。間接耦合方法的方程積分交錯進行，耦合參數(shù)在兩場間交叉迭代。Liu等[20]基于嵌入式離散裂縫模型和有限元方法建立的適用于致密油儲集層多級壓裂水平井流-固耦合的數(shù)值模擬，采用迭代耦合求解滲流力學和地質(zhì)力學控制方程；Zhang等[8]基于離散裂縫模型和有限元方法建立的適用于裂縫性油藏注水開發(fā)過程的流-固耦合數(shù)值模擬，采用顯式耦合求解。此類方法運算速度快，但精度較低，且對時間步的要求較高，一般只能以較短的時步進行求解，同時可能存在不收斂的情況。全耦合法可實現(xiàn)滲流場和應力場的同步求解，不存在時間滯后現(xiàn)象，因其計算結果穩(wěn)定、誤差較小，計算結果與實際物理過程更為一致，逐漸被國內(nèi)外學者采用[21-22]。流-固耦合問題常用的數(shù)值模擬方法包括有限差分方法、有限體積法、有限元方法等，其中基于非結構網(wǎng)格的有限元方法因其能靈活處理任意邊界條件和復雜的幾何形狀，并且相對于其他數(shù)值方法在多物理場耦合分析方面具有優(yōu)勢，被廣泛應用于致密儲集層生產(chǎn)動態(tài)模擬[5,8,23-24]。

以往的研究大多集中于應力對儲集層滲流能力的影響，但對于致密油儲集層壓裂水平井開采的動態(tài)耦合及產(chǎn)能研究相對較少。因此，本文針對致密油儲集層開采的多尺度流體流動機理、多物理場耦合作用機理和多級壓裂水平井滲流理論等關鍵問題，建立了表征雙重介質(zhì)致密油儲集層流體流動和儲集層變形的數(shù)學模型，構建了致密油儲集層物性演化模型；采用有限元方法全耦合求解模型，分析多重介質(zhì)變形規(guī)律以及儲集層流體流動的復雜機制，并對不同參數(shù)條件下的產(chǎn)能進行了研究。

1 致密油儲集層流-固耦合模型

本文應用多孔介質(zhì)彈性理論與流-固耦合作用機理，描述致密油儲集層人工裂縫、天然裂縫、基質(zhì)系統(tǒng)的變形，以及流體與固體之間的相互耦合作用。致密油儲集層多尺度空間如圖1所示，模型的基本假設條件如下：①儲集層由基質(zhì)系統(tǒng)與天然裂縫系統(tǒng)組成，二者之間存在擬穩(wěn)態(tài)竄流；②多孔介質(zhì)視為完全飽和且各向同性的線彈性體；③儲集層為彈性變形，滿足小變形假設，遵從Terzaghi有效應力原理；④單相流體等溫滲流，其流動遵循達西定律，忽略重力影響；⑤原油從基質(zhì)系統(tǒng)經(jīng)天然裂縫流向壓裂裂縫和井筒。

圖1 致密油儲集層多尺度空間示意圖（根據(jù)文獻[16]修改）

1.1 儲集層巖體變形模型

1.2 流體滲流模型

1.3 流-固動態(tài)耦合模型

流-固耦合的過程包括：流體流動引起的壓力變化會影響作用在巖石骨架上的有效應力，導致巖體發(fā)生形變；巖體變形使巖石孔隙空間、孔喉半徑、裂縫開度等發(fā)生改變，從而改變巖體的滲透率、孔隙度和裂縫導流能力，同時又會影響流體的流動（見圖2）。

圖2 致密油儲集層多物理場耦合關系示意圖

1.4 初始和邊界條件

2 模型求解與驗證

2.1 離散與求解

（6）式、（12）式和（13）式構成了描述巖體變形和流體流動的控制方程，其中每個方程均是復雜的偏微分方程，而且多物理場耦合問題又是高度非線性化的，只能求解其數(shù)值解。本文采用有限元法對控制方程進行離散，將控制方程中位移和壓力表示為結點變量的插值函數(shù)形式：

基于有限元方法與離散裂縫模型，考慮人工裂縫內(nèi)滲流以及匯-源項影響，系統(tǒng)力學平衡控制方程（6）式和系統(tǒng)內(nèi)質(zhì)量守恒方程（12）式、（13）式可在空間上分別離散成如下矩陣形式[32]：

流-固全耦合有限元平衡方程在時域進行離散，控制方程的 3個主變量（儲集層節(jié)點位移、基質(zhì)系統(tǒng)壓力、天然裂縫系統(tǒng)壓力）均以隱式方式在同一時間步內(nèi)聯(lián)立解出；與主變量相關的次要變量及其耦合關系在每個時間步進行迭代更新，即可求得系統(tǒng)位移和壓力的值。流-固耦合具體求解流程如圖3所示：①模型初始化，包括流動初始化和地質(zhì)力學初始化，得到t=0時刻空間不同位置流壓場和巖體應力狀態(tài)；②流-固全耦合，將求解過程劃分成若干個時間步增量，采用時間步推進的方法聯(lián)立耦合求解壓力場和位移場，并檢查其收斂性，若不收斂，返回重復迭代；③利用第②步得到的位移、壓力（主要變量）計算該時間步下孔隙度、滲透率、裂縫導流能力，更新參數(shù)（次要變量）；④重復②—③步驟，直至模擬時間結束。

圖3 流-固全耦合求解流程圖

2.2 模型驗證

為驗證上述流-固全耦合模型的正確性，建立了致密油儲集層多級壓裂水平井模型，模型基本參數(shù)如表1所示，數(shù)值模擬參數(shù)取值主要參考文獻[5]、[8]、[33]。模型采用三角形網(wǎng)格進行網(wǎng)格剖分，其中基質(zhì)系統(tǒng)與天然裂縫系統(tǒng)共用一套網(wǎng)格，而壓裂裂縫在網(wǎng)格劃分過程中進行降維處理；由于裂縫周圍壓力梯度較大，所以在壓裂裂縫周邊進行網(wǎng)格加密，由此得到整個模型儲集層區(qū)域網(wǎng)格劃分圖（見圖4）。由數(shù)值計算程序得到全耦合剛性模型的模擬結果，并與商業(yè)軟件Eclipse采用有限差分模型得到的日產(chǎn)油量和累計產(chǎn)油量進行對比（見圖5），本文全耦合剛性模型數(shù)值模擬結果與Eclipse的模擬結果吻合性較好。同時，在二維均質(zhì)、線性彈性、不滲透的巖石中，均勻受壓裂縫沿裂縫長度方向的開度分布滿足解析解，則將本文模型和 Jiang等[34]模型計算的裂縫開度與解析解進行對比（見圖6）。可以看出，3個模型計算的裂縫開度相近，而本文模型與解析解更接近。從以上兩個方面驗證了本文流-固耦合模型的合理性與正確性。

表1 模型基本參數(shù)表

圖4 儲集層區(qū)域網(wǎng)格剖分示意圖

圖5 全耦合剛性模型與Eclipse模型模擬結果對比

圖6 不同模型計算的裂縫開度分布與解析解對比

2.3 模型對比

本文對比了全耦合模型（模型1）、基質(zhì)+天然裂縫剛性模型（模型2）、人工裂縫剛性模型（模型3）、全剛性模型（模型 4）模擬計算的多級壓裂水平井產(chǎn)量（見圖7），并以全耦合模型累計產(chǎn)油量為基數(shù)，計算其他模型累計產(chǎn)油量與基數(shù)的相對誤差，結果見表2?？梢钥闯?，全耦合與全剛性模型產(chǎn)量預測差異大；開發(fā)前期，全耦合模型計算日產(chǎn)油量遞減快，生產(chǎn)3 000 d后全耦合模型計算累計產(chǎn)量最小，全剛性模型計算累計產(chǎn)量最高。忽略基質(zhì)與天然裂縫系統(tǒng)的流-固耦合作用，壓裂水平井生產(chǎn)3 000 d產(chǎn)量預測誤差為0.92%，主要原因在于實際開采過程中基質(zhì)與天然裂縫系統(tǒng)物性逐漸變差，其向水平井筒的供給能力逐步減弱，而模型2不考慮該因素，因此模型2計算累計產(chǎn)量較模型1偏大。忽略人工裂縫系統(tǒng)的流-固耦合作用，壓裂水平井生產(chǎn)3 000 d產(chǎn)量預測誤差為36.46%，主要由于實際開采過程中人工裂縫導流能力強，裂縫內(nèi)壓力逐漸下降、開度不斷減小，而模型 3不考慮該因素，因此模型3計算累計產(chǎn)量較模型1更大。對比模型2—模型4產(chǎn)量預測相對誤差值，可以看出模型 4產(chǎn)量預測相對誤差（38.30%）并非模型2和模型3產(chǎn)量預測相對誤差簡單的線性疊加（37.38%），進一步說明基質(zhì)、天然裂縫、人工裂縫之間的影響具有相互耦合特性。

圖7 不同模型計算的多級壓裂水平井產(chǎn)量

表2 不同模型計算的累計產(chǎn)油量對比

綜上所述，流-固耦合作用在致密油產(chǎn)能預測中的影響不可忽略，未耦合模型預測的多級壓裂水平井產(chǎn)能偏高，流-固全耦合模型預測的產(chǎn)能更為精確；人工裂縫物性是決定壓裂水平井產(chǎn)能的重要因素，有必要優(yōu)化人工裂縫參數(shù)設計。

3 儲集層物性演化規(guī)律

采用模型驗證部分所用的數(shù)值模型，計算致密油儲集層多級壓裂水平井投入生產(chǎn)后不同時刻基質(zhì)系統(tǒng)的滲透率分布和不同時刻天然裂縫系統(tǒng)的滲透率分布（見圖8—圖11）。流體經(jīng)基質(zhì)系統(tǒng)和天然裂縫系統(tǒng)流入人工裂縫及井筒時，由于人工裂縫內(nèi)導流能力大、流動阻力小，近壓裂裂縫區(qū)域存在較大的流體壓力梯度，儲集層滲透率發(fā)生顯著變化。水平井生產(chǎn)2 000 d后，基質(zhì)系統(tǒng)滲透率由 0.120×10-3μm2最小減至約0.117×10-3μm2，天然裂縫系統(tǒng)滲透率由 1.500×10-3μm2最小減至約1.070×10-3μm2，由于基質(zhì)系統(tǒng)更為致密，其滲透率變化幅度小于天然裂縫系統(tǒng)。儲集層滲透率的最小值分布于人工裂縫周圍，隨著時間的延續(xù)，滲透率降低的范圍從人工裂縫附近逐漸向儲集層四周擴大。

圖8 生產(chǎn)100 d后基質(zhì)系統(tǒng)滲透率分布

圖9 生產(chǎn)2 000 d后基質(zhì)系統(tǒng)滲透率分布

圖10 生產(chǎn)100 d后天然裂縫滲透率分布

圖11 生產(chǎn)2 000 d后天然裂縫滲透率分布

致密油儲集層多級壓裂水平井生產(chǎn)2 000 d后，人工裂縫開度由2.884×10-4m最小減至1.440×10-4m左右（見圖12、圖13）。圖14為人工裂縫開度損失幅度和導流能力損失幅度曲線。開發(fā)前期，由于人工裂縫內(nèi)流體迅速采出，其附近壓力降低較快，導致人工裂縫開度與導流能力損失幅度大；開發(fā)500 d后，人工裂縫開度損失幅度為 40.37%，導流能力損失幅度為78.80%，隨后人工裂縫物性損失幅度緩慢上升；3 000 d后，人工裂縫開度損失幅度達52.12%，導流能力損失幅度達89.02%。

圖12 生產(chǎn)100 d后壓裂裂縫開度分布

圖13 生產(chǎn)2 000 d后壓裂裂縫開度分布

圖14 不同時刻人工裂縫開度和導流能力損失幅度曲線

4 壓裂參數(shù)優(yōu)化

4.1 參數(shù)敏感性

設計案例 1—案例 9，人工裂縫壓縮系數(shù)分別為0.060，0.062，0.064 MPa-1，人工裂縫開度分別為2.884×10-4，2.978×10-4，3.065×10-4m，天然裂縫壓縮系數(shù)分別為0.15，0.25，0.35 MPa-1，啟動壓力梯度分別為0，0.05，0.10 MPa/m，模擬計算9個案例多級壓裂水平井產(chǎn)量，并以案例 1的累計產(chǎn)油量為基數(shù)，計算其他 8個案例累計產(chǎn)油量與基數(shù)的相對誤差，結果見表3。圖15—圖18對比了不同參數(shù)下生產(chǎn)3 000 d的累計產(chǎn)油量，其中案例 5累計產(chǎn)油量最大，案例 9累計產(chǎn)油量最小，可見水平井產(chǎn)量對啟動壓力梯度最敏感，人工裂縫開度次之。啟動壓力梯度減緩了儲集層壓力波的傳播，對致密油井的產(chǎn)量有較大的抑制作用，提高人工裂縫初始導流能力將有助于降低產(chǎn)量損失，同時水平井累計產(chǎn)量隨人工裂縫和天然裂縫壓縮系數(shù)的增加而逐漸降低。

表3 不同敏感性參數(shù)對產(chǎn)量的影響

圖15 不同人工裂縫壓縮系數(shù)下的累計產(chǎn)油量

圖16 不同人工裂縫開度下的累計產(chǎn)油量

圖17 不同天然裂縫壓縮系數(shù)下的累計產(chǎn)油量

圖18 不同啟動壓力梯度下的累計產(chǎn)油量

4.2 裂縫參數(shù)優(yōu)化

人工裂縫幾何分布是影響致密油儲集層產(chǎn)能的重要因素之一，本文在保持壓裂裂縫數(shù)量與壓裂裂縫長度乘積不變的情況下，設計了11個案例（案例10—案例20），分別計算了不同壓裂裂縫參數(shù)下的多級壓裂水平井產(chǎn)量，并以案例10的累計產(chǎn)油量為基數(shù)，計算其他案例累計產(chǎn)油量與基數(shù)的相對誤差，案例設計與模擬結果見表4，圖19對比了不同案例生產(chǎn)3 000 d的累計產(chǎn)油量?？梢钥闯觯咐?0累計產(chǎn)油量最大，案例13累計產(chǎn)油量最小，隨著人工裂縫數(shù)量的增多，水平井產(chǎn)量不斷增加；各案例累計產(chǎn)油量曲線在生產(chǎn)后期接近平行，因此水平井人工裂縫加密對產(chǎn)量的影響主要集中于生產(chǎn)早期。圖20展示了不同人工裂縫數(shù)量與產(chǎn)量預測相對誤差的關系。分別取產(chǎn)量預測相對誤差曲線前5個點和后6個點，做擬合曲線，形成交會圖，裂縫條數(shù)為8～9條時，產(chǎn)量預測相對誤差曲線出現(xiàn)了拐點；當裂縫數(shù)量超過10條之后，水平井產(chǎn)量增幅減弱，裂縫條數(shù)對產(chǎn)量的影響程度也逐漸減小。因此壓裂施工優(yōu)化設計需要考慮人工裂縫間距、數(shù)量、長度的綜合影響，只增加裂縫條數(shù)無法取得預期的增產(chǎn)效果；在本文模擬條件下，最佳裂縫條數(shù)為9條。

表4 不同裂縫參數(shù)對產(chǎn)量的影響

圖19 不同案例多級壓裂水平井累計產(chǎn)油量對比

5 結論

建立了致密油儲集層多級壓裂水平井多重孔隙介質(zhì)變形與流體流動的全耦合數(shù)值模型，模型綜合考慮了儲集層基質(zhì)系統(tǒng)、天然裂縫和人工裂縫的變形特征與分區(qū)流體滲流規(guī)律，能夠準確模擬預測致密油儲集層多級壓裂水平井產(chǎn)能。

模型計算結果表明，致密油儲集層水平井生產(chǎn)過程中近人工裂縫區(qū)域儲集層物性變差，特別是人工裂縫開度和裂縫導流能力損失幅度分別達到 52.12%和89.02%；流-固耦合作用在致密油儲集層產(chǎn)能預測中的影響不可忽略，模擬致密油儲集層水平井生產(chǎn)3 000 d后，全耦合模型與未耦合模型產(chǎn)能預測的誤差達38.30%；參數(shù)敏感性分析表明，致密油儲集層水平井產(chǎn)能對啟動壓力梯度最敏感，人工裂縫開度次之，提高人工裂縫初始導流能力有助于提高致密油井產(chǎn)能；致密油儲集層壓裂施工設計需考慮人工裂縫導流能力、間距、數(shù)量、長度的綜合影響，片面追求增加裂縫條數(shù)無法取得預期的增產(chǎn)效果。

符號注釋：

b——比奧校正系數(shù)；B——應變矩陣；c——壓縮系數(shù)，Pa-1；C——柔度矩陣；Cklrt,Ⅰ——基質(zhì)系統(tǒng)柔度張量，Pa-1；Cklrt,Ⅱ——天然裂縫系統(tǒng)柔度張量，Pa-1；D——彈性矩陣；Dijkl,Ⅰ,Ⅱ——雙重介質(zhì)系統(tǒng)彈性張量，Pa；fi——重力項矢量，N/m3；——邊界處已知應力，Pa；G——啟動壓力梯度，Pa/m；h——人工裂縫面高度，m；H——流體通量，kg/（m2·s）；k——迭代步；K——滲透率，m2；K0——初始滲透率，m2；m——單位張量；Mb——比奧模量，Pa；Ms,I——基質(zhì)系統(tǒng)巖石骨架體積模量，Pa；n——邊界法線方向；n——邊界單位法向量；Np——壓力場的形函數(shù)；Nu——位移場的形函數(shù)；p——孔隙壓力，Pa；p0——初始孔隙壓力，Pa；p——節(jié)點壓力，Pa；Q——源匯項，kg/（m3·s）；Q——流動方程；sF——人工裂縫存儲系數(shù)，Pa-1；t——時間，s；ts——時間步；u0——系統(tǒng)初始位移，m；u——位移矢量，m；ui,j——ui對j方向求偏導數(shù)；u——節(jié)點位移，m；us——邊界處已知位移，m；V——流體速度，m/s；wF——人工裂縫寬度，m；wF,z——z單元壓裂裂縫寬度，m；wF,0——初始人工裂縫寬度，m；W——單位體積流體質(zhì)量，kg/m3；x，y——儲集層坐標位置，m；z——裂縫單元；α——竄流系數(shù)，m-2；?——哈密頓算子；μ——流體黏度，Pa·s；υ——巖石泊松比；δ——Kronecker符號；ε——應變張量；εij——應變分量；εv——巖石體應變；ρ——流體密度，kg/m3；σβn——作用于裂縫表面的法向應力，Pa；σ——總應力張量；σ0——系統(tǒng)初始應力，Pa；σij——應力分量，Pa；σij,j——σij對j方向求偏導數(shù)，Pa/m；σi′j——有效應力張量的分量，Pa；Γm——基質(zhì)系統(tǒng)邊界；Γf——天然裂縫系統(tǒng)邊界；ΓF——人工裂縫系統(tǒng)邊界；Ωe——整個儲集層系統(tǒng)；ΩF——人工裂縫系統(tǒng)；Ωs——儲集層邊界；φ——孔隙度，%；φ0——初始孔隙度，%。下標：F——人工裂縫；i，j，k，l，r，s——Voigt 標記方法中的指標；η——基質(zhì)系統(tǒng)或天然裂縫系統(tǒng)；τ——竄流；Ⅰ——基質(zhì)系統(tǒng)；Ⅱ——天然裂縫系統(tǒng)。