基于有限元法的裂縫頁巖計算機建模仿真研究

郭 晨，劉婉萍

(長安大學 信息工程學院，陜西 西安 710064)

0 引 言

隨著世界能源需求量的不斷增加，埋藏在淺層的常規石油資源已日漸枯竭[1]，尋找有開采價值的非常規資源成為重中之重的事情。頁巖油頁巖氣是重要的非常規資源，其主要是指以頁巖為主的頁巖層系中所含的石油資源[2]。中國在四川盆地、蘇里格等油田發現了儲藏豐富的致密油氣田，與此同時，隨著勘探技術[3]的發展，頁巖油氣的開采成本也在不斷降低。面對儲藏豐富的頁巖油氣田，中國對于其開采和探測技術仍然處于摸索和依賴國外技術的階段。

頁巖作為一種致密沉積巖，是通過若干年的細微顆粒及薄層沉積形成的。大多數情況下，頁巖在其水平方向可看作同一時期的沉淀物，在垂直方向可看作不同的沉淀物。因此可以將頁巖看成不同介質組成的層狀混合物。對于頁巖混合物來說，在每層介質的形成過程中，由于各種地質結構以及氧化等原因導致在層與層之間形成很薄的薄層。此外，頁巖在形成過程中，由于擠壓等原因會形成裂縫[4]，這些裂縫會以不同形態存在于頁巖的不同位置。這些薄層和裂縫的存在對頁巖油氣田的開采和探測具有一定的影響。因此國內外許多學者對此做了研究：美國Brace和加拿大Stesky等人對電測井響應進行了數值模擬，研究表明裂縫的存在會使地層的電特性出現顯著變化[5-6]。2007年，中國石油大學的沈金松對水平裂縫進行了建模分析，分析了幾種裂縫模型在不同圍壓條件下的電阻率變化規律[7]。2017年，唐巨鵬等人研究了預制裂縫角度對頁巖水力壓裂效果影響[8]。2018年，Matthew Josh和Han Tongcheng研究了定向裂縫對合成多孔砂巖介電特性的影響[9]。

在之前的研究中，對于電參數的求解問題，學者一方面通過目前已有的成熟的解析公式進行求解，另一方面利用實驗室的測量手段進行求解。文中基于有限元法的計算機建模仿真，對不同幾何結構的頁巖模型進行仿真分析。對于裂縫頁巖電特性的研究分析，學者只是對定向裂縫的電特性進行了研究，主要對層狀結構、層狀加薄膜結構以及含有不同形狀、方向裂縫頁巖結構進行了建模與仿真，分析不同結構對頁巖電特性的影響并得出結論，結論為實際巖石的反演提供了理論依據，具有一定的應用價值。

1 有限元法

對于數值化分析方法，目前用于分析巖石電特性的主要有有限差分方法(FDM)[10-12]、有限元法(FEM)[13]、邊界積分法(BIEM)[14]以及邊界元法(BEM)[15]。文中主要用有限元法(finite element method，FEM)對頁巖模型進行建模、仿真與計算。有限元法的核心內容是將一個連續三維幾何體剖分成有限個面積不同、形狀不同但是相連的幾何體，即對求解的連續區域進行離散化，然后根據初始條件求解出每個剖分后區域電參數的近似解。最終，將每個離散區域進行總裝，并利用隨機法、迭代法、直接法等求解方法求解總裝好的方程組，得出的解就是真解的近似解，若真解的近似解在誤差范圍內，則得出的近似解就是真解，若不滿足，則需要重新計算，直至滿足設計準則范圍內的允許值為止。

有限元法的實現需要三個階段：前期處理階段、計算求解階段和后期處理階段。在前期處理階段主要是對實際問題進行建模并剖分的工作，在后期處理階段主要對結果的準確性進行評估分析，進而得到精確度高的解。

基于有限元的分析方法，國內外研究開發了許多功能強大的算法軟件，文中主要利用COMSOL Multiphysics仿真軟件進行建模與仿真。

2 仿真方法

有限元法的核心步驟分別為：確定物理場及幾何區域；建立幾何模型；有限元網格劃分;確定初始條件;求解每個子域的近似解;聯合離散域方程;求解等效電參數。

2.1 物理場選取

基于有限元法的COMSOL Multiphysics是一款高級多物理場的仿真軟件，包含結構力學、聲學、光學、電化學以及電磁分析等物理場。根據不同的實際問題選擇合適的物理場，文中主要采用電磁方法對裂縫頁巖模型的電學特性進行了建模及分析。因此，在仿真中，物理場都選取電磁波場。

2.2 仿真模型

由于頁巖具有明顯的層理結構，并且在每層形成過程中，介質表面會形成極薄的薄膜層或者由于擠壓等原因會產生裂縫層。因此在文中主要有三種頁巖結構模型：層狀結構、層狀加薄膜結構和橢球形內含物結構。圖1是結構模型示意圖。

圖1 三種頁巖結構模型

2.3 網格劃分

圖2是三種不同的網格劃分示意圖。

不同的網格劃分影響最終求得的解的精度，網格劃分越細，解的精度越高，所需的計算機內存越大，耗費的時間越長。為了保證計算解準確的情況下能快速計算頁巖電參數，經過大量的仿真計算，利用標準網格進行剖分能快速有效地計算電參數。因此，文中在仿真中，都是用標準網格進行剖分。

圖2 三種網格剖分圖

2.4 等效電參數求解

最終仿真計算真解時，采用牛頓迭代法進行方程組計算，當相對容差小于10-3時終止計算。仿真計算結束后，通過COMSOL的后期處理功能，計算不同區域、不用方向的電場模量，根據以下公式求得不同頁巖結構的等效介電常數。

(1)

其中，D是電通量密度，E是電場強度。

3 仿真結果

由于頁巖在地形中的分布特點，主要對三種頁巖結構進行建模仿真，其中包括層狀結構、層狀加薄膜結構和橢球形內含物結構。

3.1 層狀結構

在地層中，由于頁巖具有很明顯的層理結構，在仿真建模中可以等效為層狀模型，為了驗證仿真結果的正確性，將COMSOL仿真結果與復折射率模型解析公式進行了對比，復折射率公式如下：

(2)

其中，ε1、ε2分別是上下兩層介質的相對介電常數；d1、d2分別是上下兩層介質的厚度，d=d1+d2。

模型尺寸及材料設置：兩層立方體的長寬高分別為1 mm、1 mm、0.5 mm，下層介質的介電常數為20，改變上層介質的介電常數，結果如圖3所示。

圖3 COMSOL仿真結果與解析解對比

從圖中可看出，COMSOL仿真結果與解析解的結果相吻合，由此證明了該模型與方法的正確性與有效性。

3.2 層狀加薄膜結構

頁巖具有明顯的層理結構，但在每一層介質的交界面上會出現一層類似于薄膜的結構。對于含有薄膜結構的頁巖來說，薄膜的厚度僅僅是整個模型尺寸的0.1%，或者更少。但是這種薄膜層對模型的等效介電常數影響極大。因此，文中對含有薄膜層結構進行了建模仿真。模型的尺寸及材料設置為：兩層立方體的長寬高分別為1 mm、1 mm、0.5 mm，薄膜厚度為10-7mm，上層介質的介電常數為20，電導率為0.001 s/m，下層介質的介電常數為20，電導率為0.01 s/m，薄膜的介電常數為80，電導率為10-5s/m。結果對比如圖4所示。

圖4 含薄膜層與不含薄膜層仿真結果對比

由圖中可看出，含有薄膜層結構模型的仿真結果在低頻處出現了巨大增強現象，而在高頻處和不含薄膜層結構的仿真結果基本一致。

3.3 橢球形內含物結構

頁巖在實際形成過程中，由于受到擠壓等原因產生一系列裂縫，這些裂縫以不同的形狀、大小分布在頁巖中，這對頁巖的等效電參數帶來了一定的影響，進而加大了對頁巖的開采難度。在眾多的裂縫結構中，橢球形裂縫結構最為常見，其中針狀裂縫和圓盤狀裂縫最具代表性。因此，文中對不同孔隙度下，不同角度分布的含裂縫頁巖進行了建模仿真。模型如圖5所示。

模型設材料置：立方體內包裹著橢球體，立方體的介電常數為4，電導率為0.01 s/m，橢球的介電常數為80，電導率為0.1 s/m。不同裂縫三方向上等效介電常數仿真圖如圖6所示(其中，圖(a)和(b)的孔隙度為2.09%，圖(c)和(d)的孔隙度為1.68%，圖(e)的孔隙度為0.523%，圖(f)的孔隙度為8.38%)。

圖5 橢球不同角度(0°30°45°60°90°)分布示意圖

圖6 不同裂縫三方向上等效介電常數仿真圖

由圖6可以看出，含橢球形裂縫結構三方向上的等效介電常數值不相同，這主要是由于橢球在三個方向上的尺寸不同，從而對電場的阻礙程度不同所導致的。針狀裂縫Y方向上的等效介電常數小于X、Z方向上的等效介電常數值，而圓盤狀裂縫Y方向上的等效介電常數則大于X、Z方向上的等效介電常數值。針狀裂縫和圓盤狀裂縫在與水平方向的角度呈45°時，X、Z方向上的等效介電常數值相等。

綜上所述，在寬頻帶內，若低頻處的等效介電常數出現了巨大的增強現象，則在具有層理結構的頁巖中夾雜著極薄的薄膜層；若三方向上的等效介電常數值不相等，則具有層理結構的頁巖內部具有不規則的內含物，大多數情況下，等效其為橢球內含物；當Y方向上的等效介電常數小于X、Z方向上的等效介電常數值時，則在具有層理結構的頁巖中包含著針狀裂縫；當Y方向上的等效介電常數大于X、Z方向上的等效介電常數值時，則在具有層理結構的頁巖中包含著圓盤狀裂縫。研究結果可為實際巖石的反演提供理論依據。

4 結束語

主要運用有限元法對不同頁巖內部幾何結構進行建模與仿真，不僅解決了理論公式和實驗測量無法進行內部幾何結構的建模，而且形象準確地計算了不同幾何結構頁巖的電參數，根據其電特性規律總結出不同幾何分布規律。仿真結果表明，含有薄膜層結構的頁巖其電特性在低頻處的等效介電常數出現了巨大的增強現象，含有橢球型裂縫結構的頁巖其電特性在三方向上的等效介電常數值不相等。這些結論為實際頁巖的反演提供了理論依據，具有一定的應用價值。