縫洞型介質流動模擬的多尺度分解法*

張慶福 黃朝琴? 姚軍 李陽2) 嚴俠

1) (中國石油大學(華東)石油工程學院,青島 266580)

2) (中石化油田勘探開發部,北京 100728)

(2018 年8 月31 日收到; 2019 年1 月21 日收到修改稿)

縫洞型介質通常具有非均質性強、結構多尺度的特征. 傳統數值方法在解決此類多尺度流動問題時,難以兼顧計算精度與計算效率,無法實際應用. 對此,本文提出了多孔介質流體流動的多尺度分解法,并應用于縫洞介質流動模擬,能夠大幅減小計算的復雜度,同時可以通過控制均化程度控制計算精度. 該方法將求解空間分為若干個子空間的正交直和,從而獲得一個近線性的計算復雜度; 以分層計算的方式實現了快速計算,另外這種方法是一種無網格方法,具有較好的地層適應性. 同時,采用離散縫洞模型簡化縫洞結構,進一步提高了計算效率. 詳細闡述了基于多尺度分解法的多孔介質流體流動數值計算格式的建立,重點介紹了如何在不同的層次上計算基函數. 數值結果表明,本文提出的計算方法不僅能夠準確捕捉多孔介質中的精細流動特征,而且具有很高的計算效率,是一種有效的流動模擬方法.

1 引 言

縫洞型介質廣泛存在于自然界中,過去幾年間,針對其開展的數值模擬技術引起了廣泛關注,但是縫洞型介質流動模擬給常規數值方法帶來很大的挑戰. 縫洞型介質除了本身具有的非均質性,其中所包含的裂縫和溶洞橫跨多個尺度. 如果要描述縫洞介質所有尺度上的流動特征,則需要很精細的網格劃分,使用常規數值方法求解需要巨大的計算量,即便使用超級計算機也難以滿足需求. 因此,常規數值方法難以描述流體所有尺度上的流動特征. 尺度升級方法被廣泛視為能減小計算量的方法[1,2],但是由于其忽視了小尺度信息導致不能反映精細的介質非均質性. 近年來,多尺度方法被應用于流動模擬[3?7],這種方法能夠在粗網格上通過計算多尺度基函數捕捉小尺度特征,從而在保證計算精度的同時減少計算量. 然而,當粗網格和細網格尺度相差較大時,多尺度方法在差值算子及約束算子的構造方面存在著困難. 多重網格方法是公認的求解偏微分方程的最快的方法之一[8?11],因此,本文首先將多重網格法應用于非均質介質流動問題,并進一步構建了一種新的縫洞型介質流動模擬的多重網格法(多尺度分解法).

雖然多重網格方法被拓展至求解多種偏微分方程,但是當方程系數粗糙時,會導致該方法收斂性差[12,13]. 因此,對于強非均質性油藏,傳統多重網格法也存在收斂性差的問題. 近年來,雖然改進的多重網格方法在一定程度上能夠緩和粗糙系數對收斂性的影響[14?17],但能夠精確適應粗糙系數的多重網格法一直都是數學界的難題[18].

目前,主要有三種方法用于縫洞型碳酸鹽巖數值模擬: 1)等效介質模型[19]; 2)多重介質模型[20];3)離散縫洞網絡(DFVN)模型[21,22]. DFVN 模型能精細描述縫洞,基巖和裂縫系統視為滲流區域,溶洞系統視為自由流區域,是離散裂縫模型的有效延伸. DFVN 模型屬于宏觀尺度的精細模型,能克服三重介質模型和等效介質模型的瓶頸,但受目前計算機硬件和數值模擬技術的制約,只能處理小規模的縫洞型油藏. 針對此問題,對于小尺度裂縫和溶洞,對其進行等效處理,用等效后的滲透率場描述其造成的非均質性; 使用離散裂縫模型描述大尺度裂縫,準確刻畫裂縫流動狀態.

本文基于多尺度分解思想[23,24],構建了縫洞型介質多重網格計算格式,能夠準確描述介質的非均質性,并能準確模擬大尺度裂縫存在時的壓力差分布. 文中首先介紹了離散縫洞的數學模型和多尺度分解法的基本原理,然后構建了縫洞介質的多尺度分解計算格式,并給出了幾種不同的算例. 數值結果表明,多尺度分解法可以在充分反映裂縫溶洞影響的同時大幅減少計算復雜度,是一種潛力很大的流動模擬數值方法.

2 縫洞混合流動數學模型

2.1 縫洞介質數學模型

本文通過計算等效滲透率來描述介質中的小尺度裂縫和溶洞導致的非均質性; 采用離散裂縫模型表征大尺度裂縫,兩者結合形成多尺度縫洞混合模型(圖1).

圖1 縫洞型介質示意圖Fig. 1. Schematic of fractured-vuggy porous medium.

將基巖和裂縫系統視為滲流區域,其運動方程符合Darcy 定律,即

式中μ為流體黏度(mPa·s);Ki為滲透率( μ m2);vi為滲流速度(m/s); 對于巖塊l=m,對于裂縫l=f;p為多孔介質中的平均壓力(Pa);ρ是流體密度(kg/m3);f為單位質量力(m/s2).

溶洞的空間尺度較大,其流動模式視為自由流,為簡單起見本文僅考慮牛頓流體的低雷諾數運動. 為了避免引入復雜的耦合邊界條件,降低自由流-滲流耦合問題的復雜性,采用Brinkman 方程描述溶洞內的自由流:

式中vv是流體真實速度(m/s);Kv是溶洞內滲透率( μ m2);?為溶洞內孔隙度(無充填時為1);pv為流體壓力(Pa).

假設流體流動過程溫度不變,綜合以上方程可以推導出DFVN 模型的控制方程組:

式中

v為速度場. 結合合適的邊界條件可得到完整的離散縫洞網絡單相穩態滲流數學模型.

2.2 等效滲透率計算

基于流量等效計算等效滲透率的方法,往往僅考慮目標單元內的非均質性,這樣得到的等效滲透會存在誤差. 對此,本文采用超樣本方法,在包含目標單元的更大區域內建立離散縫洞網絡模型并進行有限元數值求解,從而獲取超樣本單元內的壓力和速度場; 最后,基于體積平均法計算目標網格內的等效滲透率.

在超樣本單元內結合封閉定壓邊界條件建立離散縫洞網絡模型,采用有限單元法進行求解[21],根據求得的壓力及速度場,采用體積平均法計算目標網格上的速度和壓力梯度的體積平均值[25]:

式中j=x,y(表示定壓邊界施加的坐標軸方向);Vb表示目標原始網格的體積;ul和 (?p)l分別表示目標原始網格中第l個單元的速度和壓力梯度;Vl表示第l個單元的體積;N表示目標原始網格中單元的總數. 在計算得到后,根據達西定律,可以求得等效滲透率,詳細計算過程可以參考文獻[26].

3 多尺度分解法

多尺度分解法的目的是盡可能快的、以近線性的計算復雜度求解偏微分方程. 令q為網格層數,定義I(q)為每一層的細網格單元索引i= (i1,· · · ,iq) .令i(k)=(i1,··· ,ik),I(k)={i(k):i ∈I(q)} ,其 中1 ≤k≤q,i=(i1,· · · ,iq)∈I(q);i(k,k′) 為滿足j ∈I(k′)和j(k)=i的 單 元,其 中k′ ∈{k,· · · ,q}.為I(k)×I(k+1)矩陣,滿足當j∈/i(k,k+1)時同 時其 中I(k)為I(k)×I(k)的單位矩陣.

不同網格層上的測試函數定義為

如 圖2 所 示,?被 分 為 若 干 2?k×2?k的 子 單 元其 中為的體積,為的指示函數. 對于邊界條件問題,使用非零邊界的測試函數. 通過假設A和B在進行一個博弈游戲來定義基函數,具體理論證明見文獻[24]. 根據博弈理論[27,28]并結合流動方程,B的最佳選擇為

圖2 區域?的網格剖分示意圖Fig. 2. Schematic of grid partition of solution space.

其中為基函數,定義為

在實際計算中,可以通過求解以下方程獲得

可以證明以下定理成立:

1) 最優問題(9)式可以確定唯一一個由方程(8)定義的ψi;

以均勻介質為例,首次給出了流動問題的多重網格法的基函數,如圖3 所示,每一層的基函數的共同特點是這些基函數都是局部化的. 本文沒有對基函數進行截取,當基函數包含所有信息時可視為精確的基函數,單項流動中也存在其他高精度基函數[29]. 局部化的思想是目前的熱點問題之一,該思想有利于信息的集中處理,目前大數據等研究中均有涉及.

定義空間

其中k∈{1,··· ,q}. 從(6)式可以推出此空間是嵌套的. 因此,對于k∈{1,·· · ,q ?1},

圖3 分層基函數示意圖Fig. 3. Illustration of gamblets .

令A(k)是剛度矩陣那么至此,以上嵌套公式足以進行多重網格計算,但是本文進一步對其進行多尺度分解. 定義J(k)為k-元組,其中令W(k)為J(k)×I(k)矩 陣,滿 足 當j(k?1)=i(k?1)時對于令

其中i∈J(k),定義空間

III(k)為II(k?1)關于II(k)的正交補集. 令⊕λ表示正交直和,那么可以得到空間的多尺度分解為

進一步可以推導出

此為壓力p的正交分解,p(k+1)?p(k)為流動方程在空間III(k+1)里的解. 流動方程的及多尺度分解如圖4 和圖5 所示,可以理解為在每一層上的基函數.

當介質內存在大尺度裂縫時,則采用離散裂縫模型表征大裂縫,即對裂縫進行降維處理[8]. 假設裂縫介質流動方程為FEQ,離散裂縫模型的計算格式為

其中di為裂縫開度,?m為基巖區域,?f為裂縫區域. 在每一個子空間內求解離散裂縫模型,捕捉裂縫對滲流的影響,然后將各個子空間的解加起來獲取最終縫洞介質數值解,即

4 數值算例

本節首先通過一個數值算例與參考算例的對比驗證了程序的正確性,把使用傳統模擬有限差分法求得的數值解作為參考解; 然后通過復雜離散裂縫模型數值算例進一步驗證了方法的正確性和程序的魯棒性. 對于單相流動而言,著重對比壓力差分布,為了實現壓力差的對比,設介質中一點的壓力為已知,令壓力分布為正值.

4.1 小尺度縫洞模型算例

考慮圖6(a)所示的縫洞型介質物理模型,模型大小為100 m × 100 m. 模型包含小裂縫和溶洞,不存在大裂縫,模型中的參數如表1 所列. 溶洞簡化為規則球形,溶洞內不考慮充填(滲透率無充填時為無窮大),其中的流動視為自由流. 圖6(b)是通過均化理論求得的等效滲透率場圖. 令流體自上而下流動,邊界為不流動邊界.

圖4 基函數示意圖Fig. 4. Illustration of basis functions

圖5 區域?的多尺度解示意圖Fig. 5. Illustration of solution of the multiple spaces.

令q=4,即網格系統包含4 層網格,如圖6(c)所示. 采用多尺度分解法求解得到每一層上的解,并根據(17)式獲得最終的壓力分布. 小尺度裂縫與溶洞擁有較高的滲透率,會造成多孔介質的非均質性,進而影響整體壓力分布. 在數值算例中采用L2范數計算壓力誤差

表1 裂縫性介質模型基本參數Table 1. Parameters of fractured porous medium.

式中pf為參考壓力解,pmg為多尺度分解法所計算的壓力解.

圖7 給出了參考解和多重網格解的對比. 可以看出,當k=3 時多重網格法所得壓力分布與參考解幾乎一致,因此多重網格解能夠精確地反映非均質性,從而捕捉到小尺度裂縫和溶洞的影響. 表2 列出了不同k時的計算誤差,其中當k=4時為多層分解法的精確解,當k=3 和k=2 時為均化計算,計算誤差均很小,體現了該方法的精確性. 小尺度縫洞模型的計算結果表明,此方法可以有效地進行非均質地層以及小尺度縫洞介質的流動模擬.

表2 對于小尺度縫洞模型,不同k時的計算誤差Table 2. Relative error in differentkfor a smallscale-fractured vuggy porous medium.

圖6 (a)小尺度縫洞型油藏幾何模型; (b)等效滲透率場圖; (c) 多層網格系統Fig. 6. (a) Geometrical model of a small-scale-fractured vuggy porous oil reservoir; (b) equivalent permeability of fractured-vuggy medium; (c) nested grid system.

4.2 長裂縫模型算例

算例2 考慮存在長裂縫的情況. 考慮如圖8 所示的裂縫型介質,在小裂縫和溶洞的背景下包含一條長裂縫. 裂縫開度df=1 × 10–3m,裂縫滲透率為/12 μ m2,基巖孔隙度?=0.2 . 邊界均為不流動邊界,流體從上往下流動. 采用離散裂縫模型表征裂縫.

本算例考慮長裂縫與小裂縫溶洞共存的情況.壓力分布圖9 展示了不同k時多尺度分解法與參考解的對比. 當k=3 時計算誤差為0.0131,當k=2 時計算誤差為0.0741,在進行均化的前提下依然能夠反映裂縫的存在,保證較高的計算精度.參考解與多重網格解的對比結果說明,本文構建的多尺度分解法能夠有效地處理長裂縫的情況,并能通過均化大幅減少計算量,同時保持較高的計算精度.

圖7 對于小尺度縫洞模型,參考解和多重網格解對比 (a) 參考解; (b)k=3 時的多重網格解Fig. 7. Comparison of reference solution and gamblets solution for a small-scale-fractured vuggy porous medium: (a) Reference solution; (b) gamblets solution withk=3.

圖8 長裂縫介質模型Fig. 8. Geometrical model of a fractured-vuggy porous medium with a long fracture.

4.3 大尺度縫洞模型算例

自然介質中,除了小裂縫,往往伴隨著壓裂等增產措施產生大裂縫,該算例檢驗本文提出的方法對大裂縫的模擬能力. 圖10 所示為10 m × 10 m縫洞型介質,其中包含的小尺度裂縫和溶洞的參數與算例1 相同. 在包含小尺度溶洞的同時,還包含著尺度較大的裂縫網絡系統. 裂縫網絡由6 條長裂縫組成,長裂縫開度df=1 × 10–3m,裂縫滲透率為/12 μ m2,基巖孔隙度?=0.2 . 流體從上至下流動.

多孔介質中的長裂縫作為導流通道,對壓力分布有顯著的影響,并且裂縫之間也會相互影響. 該方法在每一個子空間內求解離散裂縫模型,圖11為離散裂縫模型在各層上的解,可以很明顯地看出,各個層次上的解都能反映出裂縫的存在. 然后將k個層次上的解加起來,即得到最終的解.

圖9 對于長裂縫模型,參考解和多重網格解對比 (a)參考解; (b)k=3 時的多重網絡解; (c)k=2 時的多重網格解Fig. 9. Comparison of reference solution and gamblets solution for a fractured-vuggy porous medium with a long fracture: (a) Reference solution; (b) gamblets solution withk=3; (c) gamblets solution withk=2.

圖10 大尺度縫洞介質幾何模型Fig. 10. Geometrical model of a large-scale-fractured vuggy porous medium.

圖12 給出了參考解和多重網格解的對比,可以看出,本文構建的數值算法在捕捉小尺度裂縫和溶洞的同時,能夠精確模擬長裂縫對壓力場的影響,并反映裂縫間的相互作用. 結合表3 可知,當均化程度較大(k=1)時存在較大誤差. 圖13 給出了沿x=5 m 的參考解與不同k時多重網格解的對比. 裂縫在流動過程中可視為一個等勢體,對應于圖11 中較平緩的部分,與實際符合,隨著均化程度加大誤差變大. 該算例展示了本文構建的數值方法對長裂縫的精確模擬能力.

實際上,對于任意形狀的一個裂縫性介質,一旦求出剛度矩陣及對應的W和 π ,就可以有效地進行任意形狀的裂縫性介質數值模擬. 因此,該方法能夠處理形狀復雜的裂縫系統. 但是,對于裂縫分布極其復雜的裂縫系統,尤其裂縫之間距離非常小的情況,多層網格剖分會存在一定困難,造成分層的基函數不適應高滲透的裂縫,進而影響該方法的計算精度.

圖11 大尺度縫洞介質模型在各層上的解Fig. 11. Solutions in different levels for a large-scale-fractured vuggy porous medium.

圖12 對于大尺度縫洞介質模型,參考解和多重網格解對比 (a) 參考解; (b)k=3 時的多重網格解; (c)k=2 時的多重網格解;(d)k=1 時的多重網格解Fig. 12. Comparison of reference solution and gamblets solution for a large-scale-fractured vuggy porous medium: (a) Reference solution; (b) gamblets solution withk=3; (c) gamblets solution withk=2; (d) gamblets solution withk=1.

表3 對于大尺度縫洞介質模型,不同k時的計算誤差Table 3. Relative error with differentkfor a large-scalefractured vuggy porous medium.

圖13 沿x=5 m 的參考解和不同均化程度下多重網格解對比曲線圖Fig. 13. Comparison of reference solution and gamblets solution with differentkwhenx=5 m.

5 結 論

1)隨著現代地質建模技術的發展,縫洞型地質模型可能包含數百萬甚至數億個網格,采用傳統的數值方法對其進行求解,計算量巨大,甚至超出當今計算機的計算能力. 本文首次提出縫洞介質的多尺度分解計算模式,構建局部化的基函數,該方法能夠快速求解流動方程,非常適合對強非均質型油藏進行精細流動模擬;

2)結合離散裂縫模型和等效介質模型,既可以精確反映大尺度裂縫的流體流動,又可以捕捉小尺度溶洞和裂縫的影響,從而可以對縫洞型介質進行高效流動模擬;

3)每一層網格的基函數是獨立計算的,可以采用并行計算得到,進一步減少了計算量. 因此,本文提出的方法對于縫洞型油藏數值模擬有很高的潛在價值. 滲透率場對計算精度有一定影響,因此提高算法的魯棒性是研究重點. 下一步可采用壓力梯度對比的形式驗證方法對多相流動計算的精確性[30].