劉子平, 李俊翔, 劉 琦, 馬 強(qiáng), 魏華祎

(1. 中國石油川慶鉆探工程有限公司，成都 610051; 2. 四川大學(xué)數(shù)學(xué)學(xué)院，成都 610064;3. 湘潭大學(xué)數(shù)學(xué)與計(jì)算科學(xué)學(xué)院，湘潭 411105)

0 引言

在石油工程中，油藏模擬器是描述油藏中復(fù)雜流動(dòng)機(jī)理的工具，大多數(shù)石油與天然氣公司都使用油藏模擬器來設(shè)計(jì)與預(yù)測石油或者天然氣的產(chǎn)出與效果。我國蘊(yùn)藏著豐富的頁巖氣資源，其大規(guī)模開發(fā)將是我國滿足天然氣需求的重要途徑和保障，但相對(duì)于傳統(tǒng)油氣資源，人們對(duì)于頁巖氣開采中流體流動(dòng)機(jī)理認(rèn)識(shí)還不是很成熟，而數(shù)值模擬是研究頁巖氣藏流動(dòng)機(jī)理的重要手段，目前相關(guān)的數(shù)值模擬算法和軟件與國外還有很大差距。

大多數(shù)模擬器都關(guān)注于孔隙介質(zhì)中的流動(dòng)過程，而結(jié)構(gòu)本身的地應(yīng)力則沒有得到相對(duì)應(yīng)的關(guān)注。事實(shí)上，多孔介質(zhì)是可變形的，需要我們考慮流體流動(dòng)與巖石變形之間的相互作用。在常規(guī)的固化巖石儲(chǔ)層中，我們可以調(diào)整流動(dòng)方程來充分描述地質(zhì)結(jié)構(gòu)中的低孔隙壓縮性，剛度過載以及輕微的壓降。但在非常規(guī)儲(chǔ)層比如頁巖氣的開發(fā)中，非固化的介質(zhì)、大的壓降以及水力裂縫使得流動(dòng)與結(jié)構(gòu)變形之間的相互依賴與相互作用變得顯著[?,?,?,?]，從而必須對(duì)此現(xiàn)象進(jìn)行顯式地描述。這就需要我們建立流動(dòng)與地應(yīng)力相互耦合的數(shù)學(xué)模型并發(fā)展相應(yīng)的數(shù)值計(jì)算方法來分析此類復(fù)雜問題。

在進(jìn)行數(shù)值模擬中，大多數(shù)油藏系統(tǒng)都是非均勻的，它們由其流動(dòng)性質(zhì)上的空間變化來刻畫，比如滲透率需要使用二階張量進(jìn)行描述。為此，在模擬中，我們有必要獲得非常精確的速度分布，因?yàn)榱黧w的流動(dòng)路徑將由速度方向決定。此外，得到關(guān)于飽和度解的準(zhǔn)確描述也總是需要的。在石油工程數(shù)值模擬中，常常使用有限體積法或者單元中心有限差分離散來描述非變形的多孔介質(zhì)，因?yàn)檫@些方法可以保證局部的質(zhì)量守恒[?]，而對(duì)于描述固體的變形，有限元方法則是更好的選擇[6–7]?；旌嫌邢拊椒ㄍㄟ^將質(zhì)量守恒方程與達(dá)西定律方程同時(shí)進(jìn)行求解，既能滿足局部的質(zhì)量守恒也可以給出速度場更加精確與連續(xù)的描述[?,?,?]。由此，混合有限元方法也在描述多孔介質(zhì)中的流動(dòng)與巖體形變的耦合中得到了大量的使用[?,?,?,?,?,?,?,?,?,?]。

本文針對(duì)頁巖氣勘探開發(fā)過程中水氣或水油兩相流耦合的地應(yīng)力計(jì)算與模擬問題展開研究，構(gòu)建地質(zhì)孔隙結(jié)構(gòu)非線性兩相流耦合地應(yīng)力數(shù)學(xué)模型，提出合理的離散數(shù)值計(jì)算格式，構(gòu)建相應(yīng)的非線性混合有限元算法，并針對(duì)典型的兩相流地應(yīng)力問題進(jìn)行計(jì)算模擬。

1 控制方程

為模擬帶裂縫儲(chǔ)層中的孔隙應(yīng)力問題，孔隙介質(zhì)中的流體流動(dòng)與地應(yīng)力產(chǎn)生了耦合，而地應(yīng)力導(dǎo)致了巖石結(jié)構(gòu)的線彈性變形。流體的流動(dòng)與結(jié)構(gòu)的應(yīng)力應(yīng)變充分耦合，需要在統(tǒng)一的系統(tǒng)下進(jìn)行求解以保證精度與數(shù)值穩(wěn)定性。一般地，我們首先推導(dǎo)兩相流作用下的地應(yīng)力耦合模型，基于兩相流的質(zhì)量守恒方程為

2 混合有限元分析與算法

下面對(duì)耦合地應(yīng)力模型構(gòu)建有限元算法，令未知量S0、vt、p、u分別屬于函數(shù)空間

注意這里vini= 0 下標(biāo)i為坐標(biāo)分量求和，并非流體標(biāo)識(shí)0 與1，表示在邊界上法向速度的約束，即邊界沒有流體流出。由(??)，可得到問題的弱形式為

其中σt=bp0I+σ0+σeff，(·,·)Ω表示區(qū)域Ω上內(nèi)積，〈·,·〉Γ2表示邊界Γ2上內(nèi)積，測試函數(shù)

注意到左端整體矩陣中的子矩陣包含了未知量，并且右端項(xiàng)也包含了非知量。為簡單起見，我們采用固定點(diǎn)迭代計(jì)算求解該非線性問題。由此建立有限元計(jì)算流程如下：

3 數(shù)值算例與討論

3.1 無裂縫水驅(qū)油模型

下面給出數(shù)值算例，首先我們考慮一個(gè)均勻無裂縫區(qū)域下的水驅(qū)油兩相流地應(yīng)力模型。設(shè)模型中0 表示水相，1 表示油相，區(qū)域Ω=[0,10]2m2由32×32×2 的一致三角形單元離散，如圖1 所示。

圖1 正方形無裂縫巖石區(qū)域，使用三角形對(duì)區(qū)域作離散

表1 為巖石，水與油的材料參數(shù)，其中λ表示彈性拉梅第一常數(shù)，μ表示巖石剪切模量，由此可以得到四階彈性系數(shù)張量為

表1 巖石與液體參數(shù)

C=Cijkl=λδijδkl+μδikδjl+μδilδjk,

圖2 給出了在模擬1 000 天注水后區(qū)域中的水飽和度、速度、壓力、速度模量以及地應(yīng)力模量的分布。圖2(a)中飽和度的分布給出了注入水前端的均勻一致傳播，我們可以看到與注入井相距相同的徑向距離有相同的飽度值。這是因?yàn)閮?chǔ)層具有均勻的性質(zhì)，并且沒有流體流出。圖2(c)中的壓力分布給出了一個(gè)在注入井處與產(chǎn)出井處的壓力梯度，速度最高處位于注入井與產(chǎn)出井附近。計(jì)算得到的飽和度，壓力，速度解顯示了我們的算法能給出孔隙介質(zhì)多相流動(dòng)的準(zhǔn)確物理行為。

圖2 均勻巖石結(jié)構(gòu)注水1 000 天后的計(jì)算解

由于整個(gè)結(jié)構(gòu)具有各向同性的彈性剛度，位移解與地應(yīng)力非常光滑并且依賴于壓力梯度，并且如我們所期望的，位移顯示出關(guān)于對(duì)角線x=y的對(duì)稱性。該模型的計(jì)算結(jié)果與文獻(xiàn)[11]的結(jié)果相符，顯示了我們構(gòu)造的算法的正確性。

3.2 交叉裂縫水驅(qū)油模型

接下來，我們考慮如圖3 所示的非均勻巖石儲(chǔ)層區(qū)域。在區(qū)域中存在一個(gè)交叉裂縫，流體在該裂縫中具有更大的滲透率，設(shè)其為其它區(qū)域巖石滲透率的3 倍，模型的其它參數(shù)仍如表1 保持不變。由于裂縫尺寸較小，我們對(duì)區(qū)域進(jìn)行如圖3 所示的非結(jié)構(gòu)化有限元網(wǎng)格剖分。與上節(jié)中的條件一樣，我們?nèi)匀辉谧笙陆沁M(jìn)行注水，油在右上角流出。模擬時(shí)間設(shè)置為152 天。

圖3 交叉裂縫巖石區(qū)域以及非結(jié)構(gòu)化有限元網(wǎng)格離散

圖4 給出了在152 天時(shí)計(jì)算得到的水飽和度、速度、壓力、位移以及地應(yīng)力解?？梢钥吹?，裂縫區(qū)域中水的流動(dòng)速度遠(yuǎn)大于在巖石中的流動(dòng)速度，因此通過在左下角的持續(xù)注水，流體迅速沿對(duì)角線x=y上的裂縫流動(dòng)。而在巖石中水的流動(dòng)也受到影響并有向?qū)蔷€靠攏的趨勢。而在另外一條裂縫x+y=10 處，水的飽和度與無裂縫區(qū)域比較也明顯增加。圖4(b)所示的速度分布明顯顯示出裂縫區(qū)域的速度大于無裂縫區(qū)域的速度。由于裂縫分布的對(duì)稱性，壓力、位移與地應(yīng)力也呈現(xiàn)出對(duì)稱性，并且由于區(qū)域的非均勻性，計(jì)算得到的解明顯沒有均勻巖石區(qū)域的解光滑。

圖4 非均勻交叉裂縫區(qū)域注水152 天后的計(jì)算解

4 結(jié)論與展望

本文建立了考慮裂縫儲(chǔ)層的地應(yīng)力兩相流耦合模型，推導(dǎo)了地應(yīng)力模型的非線性混合有限元離散格式，構(gòu)建了二維非線性有限元計(jì)算流程，針對(duì)典型的帶裂縫的水驅(qū)油模型問題，實(shí)現(xiàn)了非線性有限元計(jì)算模擬，驗(yàn)證了算法的正確性與有效性。

本文主要關(guān)注于有限元方法的工程應(yīng)用，所考慮的模型相對(duì)較復(fù)雜，難以構(gòu)造比較合適的真解實(shí)例，我們也沒有對(duì)算法進(jìn)行理論分析，主要通過檢查計(jì)算結(jié)果是否符合物理實(shí)際做為評(píng)價(jià)結(jié)果正確性的標(biāo)準(zhǔn)。當(dāng)然，在未來的工作中，我們不但要考慮更實(shí)際的三維問題，還要從理論層面對(duì)算法進(jìn)行適當(dāng)?shù)姆治?，以期建立本文算法收斂性理論?/p>

關(guān)于下一步的工作，一方面針對(duì)三維帶裂縫的非均勻儲(chǔ)層區(qū)域進(jìn)行有效的建模，使用自適應(yīng)有限元網(wǎng)格對(duì)區(qū)域進(jìn)行更精細(xì)的離散，另一方面希望把適用于多邊形和多面體的非標(biāo)準(zhǔn)有限元應(yīng)用到該非線性問題的求解當(dāng)中，并設(shè)計(jì)更高效穩(wěn)健的算法和程序。