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

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

0 引言

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

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

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

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

1 控制方程

為模擬帶裂縫儲(chǔ)層中的孔隙應(yīng)力問(wèn)題，孔隙介質(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，表示在邊界上法向速度的約束，即邊界沒(méi)有流體流出。由(??)，可得到問(wèn)題的弱形式為

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

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

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

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

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

圖1 正方形無(wú)裂縫巖石區(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ì)，并且沒(méi)有流體流出。圖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ū)油模型

接下來(lái)，我們考慮如圖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)速度，因此通過(guò)在左下角的持續(xù)注水，流體迅速沿對(duì)角線x=y上的裂縫流動(dòng)。而在巖石中水的流動(dòng)也受到影響并有向?qū)蔷€靠攏的趨勢(shì)。而在另外一條裂縫x+y=10 處，水的飽和度與無(wú)裂縫區(qū)域比較也明顯增加。圖4(b)所示的速度分布明顯顯示出裂縫區(qū)域的速度大于無(wú)裂縫區(qū)域的速度。由于裂縫分布的對(duì)稱性，壓力、位移與地應(yīng)力也呈現(xiàn)出對(duì)稱性，并且由于區(qū)域的非均勻性，計(jì)算得到的解明顯沒(méi)有均勻巖石區(qū)域的解光滑。

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

4 結(jié)論與展望

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

本文主要關(guān)注于有限元方法的工程應(yīng)用，所考慮的模型相對(duì)較復(fù)雜，難以構(gòu)造比較合適的真解實(shí)例，我們也沒(méi)有對(duì)算法進(jìn)行理論分析，主要通過(guò)檢查計(jì)算結(jié)果是否符合物理實(shí)際做為評(píng)價(jià)結(jié)果正確性的標(biāo)準(zhǔn)。當(dāng)然，在未來(lái)的工作中，我們不但要考慮更實(shí)際的三維問(wèn)題，還要從理論層面對(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)用到該非線性問(wèn)題的求解當(dāng)中，并設(shè)計(jì)更高效穩(wěn)健的算法和程序。