考慮熱流固耦合作用的多孔介質孔隙尺度兩相流動模擬1)

蔡少斌 楊永飛 劉 杰

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

引言

經過多年的勘探與開發,我國淺層和中淺層的油氣資源已經進入了較為成熟的開發階段.伴隨著勘探開發的深入,人們認識到在深層地層中同樣孕育著豐富的可采能源,例如深層油氣藏、地熱資源等.深層能源儲藏具有可觀的挖潛能力,是未來助力我國增儲上產的重要來源[1].

深部地層內的流體流動同時受到了溫度、壓力及應力場的綜合影響[2].伴隨地層深度加深,地層溫度不斷地升高.不同地區的地溫梯度不同,一般處于20～30 °C/km,部分異常區域可達40～80 °C/km.當儲層深度到達深層時,地層溫度可以達到100 °C 以上.綜合考慮以上因素,流體在巖石多孔介質內的流動受到了多物理場的綜合影響[3].孔隙尺度下考慮兩相流體[4]在熱場、流場、應力場多場耦合作用下的流動模擬研究,為解釋油氣資源在深部地層的運移規律提供了一種解決方案.同時,多場耦合的研究也有益于其他地下工程的研究,如天然氣水合物開采[5]、二氧化碳地質埋存工程[6]等.

最早提出流固耦合理論研究模型的是Terzaghi,他提出了飽和多孔介質的有效應力公式和一維有效固結模型[7].李錫夔等[8-9]提出了一種混合有限元方法計算非飽和多孔介質中的熱流固過程.宋睿等[10-11]通過使用ANSYS 和CFX 求解熱?流?固耦合模型,結果表明隨著巖石有效壓力和溫度的增加,孔隙度和滲透率都會下降.Wu 等[12]建立了具有周期性孔隙形態的分形理論模型,用于模擬復雜多孔介質中的輸運,并基于分形模型,進行了熱?流?固耦合流動模擬.鄧佳等[13]通過考慮頁巖儲層的應力敏感特征,分析了儲層的應力敏感性對頁巖氣儲層產量的影響.王沫然和王梓巖等[14]采用了跨尺度混合模擬算法,研究了深層頁巖氣藏的流固耦合問題.樊冬艷等[15]研究了基于離散裂縫網絡模型的干熱巖儲層熱?流?固耦合問題,結果表明出口端流量及溫度會對流量產生影響.孫致學等[16]將裂隙巖體視為離散裂隙網絡和基質巖體的雙重介質模型,研究了基于熱?流?固耦合模型的增強型地熱系統的生產過程.郭穎等[17]基于Biot 方程、Darcy 定律及Lord-Schulman 準則研究了土壤在載荷作用下物理量的變化規律,結果表明載荷對孔隙度和滲透率皆有明顯影響.

Darcy-Brinkman 方法[18]被廣泛地應用于多孔介質多尺度多場耦合流動模擬研究中,如模擬酸巖反應過程[19]、流固耦合過程[20]等.為了能夠定量表征溫度場、應力場及流場對發生在巖石多孔介質內的多相滲流的影響,本研究使用了Darcy-Brinkman-Biot[21]耦合Duhamel-Neumann[22]熱彈性應力準則的模型開展模擬工作.采用體積平均方法(volume averaging method,VAM)[23]引入系列參數用以表征控制單元內的油相、水相及巖石固體顆粒的分布.如圖1 所示,模型中水相、油相、固體顆粒相在控制單元內共存,水相與油相占據孔隙空間.孔隙空間內的油水兩相流動采用流體體積模型(volume of fluids,VOF)進行求解.巖石固體顆粒由于流固耦合作用產生的彈性變形采用Biot 彈性變形方程求解.溫度場的變化引起的熱彈性應力則使用了Duhamel-Neumann 熱彈性應力準則進行求解,從而實現了在孔隙尺度上的熱流固耦合流動模擬.

圖1 模型示意圖Fig.1 Illustration of model

綜上,本工作采用Darcy-Brinkman-Biot 方法,耦合了Duhamel-Neumann 熱彈性應力準則,推導了考慮熱流固三場耦合作用下的多孔介質兩相流動模擬模型,研究了多場耦合作用對油水兩相滲流的影響,以期為油氣田開發工程提供參考數據.

1 數學模型

假設模型為飽和流體的多孔介質,其中油相Vo、水相Vw及固體相Vs分別占據模型的不同空間位置,它們之間的關系為

其中,V為模型的總體積,m3;Vo,Vw及Vs分別代表油相、水相及固體相占據的體積,m3.

油相與水相統稱為流體相,以相分數 φf表示;固體相以固體體積分數 φs表示(如圖1 所示),二者的計算式分別為

顯然, φf與 φs兩者之間具有以下關系

通過 φf及 φs可以區分流體及固體控制的區域.在流體控制區域中,分布著水相與油相,分別定義油相飽和度為 αo及水相飽和度 αw,二者的計算公式分別為

顯然, αw與 αo之間具有如下關系

因此,通過圖1 中所定義的參數,可以判斷油相、水相及固體相在模型中的分布,其判定標準為

1.1 不可壓縮兩相流體控制方程

對于不可壓縮、不可混溶的兩相牛頓流體,可以采用流體體積方法[24]描述多相流動過程.采用Carrillo 等[21]基于體積平均方法推廣的流體體積方法公式來描述模型內的兩相流動過程,即

式中,t為時間,s;Uf為流體流速,m/s;Ur為相對流速,m/s;ρf為流體密度, m3/kg;p為流體壓力,Pa;g為重力加速度,9.8 m/s2;為剪切應力張量,N/m2;Dm,n表示m相對n相的動能交換,m,n= w,o,s.將動能交換項Dm,n展開,方程(11)可以寫成如下形式

式中, γ 表示界面張力,N/m;nw,o表示油水界面單位法向量.通過計算油相、水相及固體相之間的動能交換,流體與固體之間的耦合效應便可以通過模擬獲得.

1.2 線彈性固體變形控制方程

本研究采用的模型符合連續介質假設,為了研究變形多孔介質對兩相滲流過程的影響,研究假設巖石骨架在彈性變形階段為線彈性變形.采用Carriillo和Bourg[25]體積平均方法推廣的Biot 變形方程

式中,Us為固體顆粒運動速度,m/s;ρs為固體顆粒密度, k g/m3;σ 為體積平均后的彈性應力張量,N/m;τs(N/m2) 為Terzaghi 有效應力張量,τs=Pconf?Ip(Pconf為圍壓,Pa;I為Biot 系數,表示壓力在流體與固體之間的傳遞關系,在本研究中假定I=1,即壓力在流體與固體之間完全傳遞);Bs,i為固體相與油相或水相之間的動能交換項,i= w,o,假設流體的所有剪切動能損失都轉移到固體上,根據式(12)中流體與固體之間的動能交換項,可以獲得如式(16)所示的平衡方程

1.3 線性熱彈性應力控制方程

在流固耦合的基礎上,進一步考慮巖石骨架在熱場作用下產生的熱彈性應力對兩相流動過程的影響.

傳熱過程中的動能方程為

流體換熱過程中的熱傳導方程為

其中,c為比熱容,J/(kg·K?1);T為溫度,K;K為導熱系數W/(m·K).

溫度對骨架產生的熱應力為

其中, σt為熱應力,N/m2;T為溫度,K;β 為熱應力系數.

當巖石的溫度發生改變時,熱應力作用于于巖石固體顆粒上,發生形變.在考慮流固耦合作用的基礎上,固體受到的應力可以分為兩部分:一部分是由于流固耦合作用發生的應力, σε;另一部分是熱應力, σt.模型考慮了Duhamel-Neumann[22]熱彈性應力準則,用以更新熱應力影響下的總應力

其中, σi,i=x,y表示總應力在x和y方向上的分量,N/m2.

2 模型與參數

2.1 模型參數

為了模擬油水兩相在多孔介質內受熱?流?固耦合作用下的流動過程,采用二維Berea[26]砂巖巖心切片作為模擬模型,開展流動模擬工作.圖2 展示的是本研究所采用的巖心切片示意圖,切片體素大小為500×500,分辨率為3.003 5 μm,孔隙度 φf=0.23.模型中,藍色部分代表水相,紅色部分代表油相,白色部分代表巖石骨架.為了使得模型能夠穩定計算,在模型的進出口部分,分別沿著y軸方向設置了寬度為100 體素大小的入流區域與出流區域(如圖2中右上藍色區域及左下紅色區域所示).入流區域的設置使得模型在初始時刻始終保持模型內油水界面的存在,從而起到了穩定計算的作用.

圖2 模型示意圖Fig.2 Illustration of simulation model

模型模擬的過程為水相驅替油相的過程,以速度邊界條件作為控制條件,在入口處以0.001 m/s 的注入流速注入水相.模型入口溫度設置為350 K,出口溫度設置為300 K,模型內部初始溫度場設置為300 K.模型網格采用六面體結構化網格,網格數為10000 個.

2.2 模擬參數

模型采用的模擬參數如表 1 所示.模型設置水相為潤濕相,其潤濕角設置為45°.流體物性參數假設在溫度場的作用下為恒定值.巖石的物理性質設置參考了砂巖實驗的數據[27-28],其中熱擴散系數通過巖石密度、比熱容及導熱系數確定,計算公式為αt=K/(ρsc).模擬直至油水兩相飽和度不變停止.

表1 模擬參數設置Table 1 Simulation parameters

2.3 模型求解

模型的求解流程如圖3 所示.模型求解借助了著名的基于有限體積方法(finite volume method,FVM)的開源CFD 軟件OpenFOAM[29]完成.在每個求解時間步內,分別依次對流場、熱場及固體場進行求解.

圖3 模型求解過程示意圖Fig.3 Illustration of solution algorithm

對于流體場的求解,首先使用MULES (multidimensional universal limiter of explicit solution algorithm)算法對流體體積分數進行顯式求解.MULES方法能夠很好的保證流體體積分數的有界性,從而減少計算越界的問題.對于整個流體控制方程的求解采用了SIMPLE (semi-implicit method for pressure linked equations)算法.SIMPLE 算法引入了壓力修正方程,使得流體速度場可以不斷地進行修正,從而使計算得到的速度滿足控制方程收斂條件.將計算得到的速度場傳入溫度場進行求解,最后進行固體場的求解.求解固體變形的算法采用了Jasak 和Weller[30]在構建fsiFoam 時所用的算法.求解時間步長滿足柯朗?弗里德里希斯?列維數(Courant-Friedrichs-Lewy number).對于二維模型,CFL 數計算公式為

式中,Ux和Uy分別代表了速度在x,y方向上的分量,m/s;Δt表示求解時間步長,s;Δx和 Δy分別代表了x,y方向上的求解空間長度,m;C 代表判定標準,為常數,在本研究中,經過反復測試,取C = 0.15 模型能夠穩定計算.

3 結果與討論

3.1 模擬結果

圖4 展示的是恒定流速0.001 m/s 為注入流速,水相潤濕角設置為45°,入口溫度為350 K,初始內部溫度與出口溫度為300 K 的模型的模擬結果.從圖中可以看出,水相在巖心切片的右側形成了優勢主力流動通道,驅替前端在注入速度的作用下發生了潤濕滯后.隨著驅替過程的進行,水相不斷地將油相驅替出巖心切片,最終當水相到達出口時,主力流動通道達到穩定狀態,模型內的兩相飽和度趨于不變.

圖4 不同模擬時刻的相分布Fig.4 Simulation results of phase distribution at different time steps

圖5 展示的是圖4 中的模型在0.15 s 即兩相飽和度達到穩定狀態時的模型內的應力、應變分布狀態.兩圖中不同顏色的線段代表應力、應變數值大小的等值線分布.從應力分布圖中可以看出,在驅替的進行過程中,流體對于巖石骨架產生了擠壓應力.在流動通道的兩側,巖石骨架由于受擠而產生了相應的壓力梯度.巖石骨架的應力主要由流固作用與熱固作用產生,流體擠壓巖石產生了彈性應力,溫度的增加使得巖石產生了向外膨脹的熱應力,兩者方向相反,互相抵消一部分影響.在本例中,巖石的應力為受擠應力,即流體對巖石骨架的影響大于溫度對巖石骨架的影響.從應變分布圖中可以看出,巖石的彈性應力使得巖石發生了彈性形變,發生形變的數量級分布在0.1～1 μm 之間.

圖5 0.15 s 時應力及應變模擬結果Fig.5 Simulation results of (a) stress and (b) strain profile at 0.15 s

相對滲透率是反應兩相流體在巖石多孔介質內相對流動能力的重要參數.圖6 展示的是圖4 模型的歸一化相滲曲線.首先引入有效含水飽和度Se用以計算歸一化相滲曲線,Se的計算公式如式(22)所示

圖6 歸一化相滲曲線與實驗結果對比Fig.6 Normalized relative permeability curve vs.experimental result

式中,Sw表示含水飽和度;Swc表示束縛水飽和度,本例中,由于模型初始狀態為全部飽和油相,因此Swc= 0;Sor表示殘余油飽和度,取驅替至油相飽和度不變時的飽和度.

相滲曲線的計算采用的是van Genuchten[31]模型,計算公式如式(23)所示

式中,Krw與Kro分別代表水相與油相的相對滲透率;m為衡量流體潤濕能力的參數,m數值越高代表流體的潤濕性越好,相對滲流能力越高.

圖6 展示模型歸一化后的相滲曲線與Oak[32]的Berea 砂巖實驗結果的對比,為了方便對比,實驗數據使用式(22)進行了歸一化.模型使用式(23)計算相滲曲線,發現當m= 0.75 時,模擬結果與實驗所獲得的相滲曲線吻合較好.模擬所獲得的相滲曲線等滲點Se= 0.58,實驗獲得的等滲點Se約等于0.62.在低含水飽和度下,水相對流動能力較小,油相相對流動能力較大.隨著驅替過程的不斷地進行,水相的相對滲流能力不斷增加,油相的相對滲流能力快速下降.在驅替過程達到末期時,水相已經形成了優勢流動通道,相對滲流能力遠大于油相,油相基本不再流動.

對比實驗結果,模擬結果對油相的整體相對滲流能力的預測偏低,對水相相對滲流能力的預測在等滲點的左側偏高,在等滲點的右側偏低.在實際巖心驅替的過程中,在水濕巖心條件下,水相在驅替早期含水飽和度較低的情況下,會首先潤濕壁面,在壁面處形成潤濕層,從而造成了水相在貼近壁面處流動,而油相占據孔隙中央區域流動的現象.這樣的現象便造成了在驅替實驗的早期,水相相對滲透率較低的情況.隨著驅替過程的進行,水相逐漸填充了孔隙空間,油相被驅替出巖心,流動阻力逐漸減小,使得水相相對滲透率大幅增加,導致驅替階段后期,水相相對滲透率較高.而在考慮了熱流固耦合的綜合影響下,在驅替的早期,由于固體發生了彈性形變(圖5),使得水相的流動通道變大,因此相對滲流能力略有增加.在達到等滲點后,隨著溫度在流場內的傳遞,巖石受熱應力的影響,形變大小減小,水相流動通道減小,使得相對滲流能力略有減小.

3.2 注入PV 數的影響

注入PV 數(pore volume number)指的是累計注入的流體體積占據孔隙空間的倍數.PV 數是衡量流體注入能力及驅替階段的重要參數.PV 數越高,代表越多的注入流體進入了巖心,使得巖心內預先飽和的流體飽和度下降.并且隨著注入流體的不斷增加,注入流體在巖心內的分布面積變大.

保持圖4 中的模擬條件不變,圖7 展示的是PV 數為5～160 之間的模擬結果.圖7(a)～圖7(c)分別展示的是模型沿著驅替方向中心軸線的溫度、應力及應變分布(圖中應力、應變為零的部分為孔隙空間).

圖7 不同注入PV 數的模擬結果Fig.7 Simulation results under different PV numbers

從圖7(a)可以看出,隨著注入PV 數的不斷增加,熱量不斷地沿著驅替方向傳播,巖心固體不斷地被加熱.如圖7(b)所示,隨著PV 數增加,固體溫度升高,受到的熱膨脹應力不斷上升,總應力不斷下降,并且沿著驅替方向,總應力也在不斷地下降.根據式(19)和式(20),熱膨脹應力與固體所受的流體所施加的應力方向相反,因此兩者會相互抵消,造成隨著溫度的上升,而總應力下降的情況.由于總應力隨著PV 數的增加而下降,如圖7(c)所示,發生的彈性形變也在不斷地下降,與應力的變化呈現出相同的規律.

3.3 注入溫差的影響

為了探究溫度變化對兩相流體在多孔介質內流動過程的影響,本研究設置了五組模擬,保持表 1 中的參數設置不變,分別設置入口溫度為350 K,400 K,450 K,500 K,550 K 開展模擬.水相潤濕角為45°,入口端注入流速為0.001 m/s.

圖8 展示的是五組不同注入溫差條件模型在100 PV 的模擬結果,五組模型中溫度場沿著驅替方向中心軸線的分布情況如圖8(a) 所示;圖8(b) 和圖8(c)分別展示了應力、應變沿著驅替方向的分布情況;圖8(d)是孔隙度與注入溫差的關系.

圖8 不同注入溫度下的模擬結果Fig.8 Simulation results under different temperature difference

隨著注入溫差的不斷增加,模型內的溫度梯度不斷增加(如圖8(a)所示).此時,熱膨脹應力成為了造成巖石彈性變形的主導力.如圖8(b)和圖8(c)所示,隨著溫度的增加,巖石固體所受的應力、應變隨之增加.與圖7(b)和圖7(c)相比,當模型的注入溫度達到400 K時,熱膨脹應力已經抵消了流體與骨架之間的流固耦合相應,此時流體對骨架產生的應力已經不足以抵消熱膨脹應力,因此巖石的應力、應變隨著溫度的升高而增加.同時,對比圖7(b)、圖7(c)與圖8(b)、圖8(c)可以發現,圖7 中原本孔隙占據的空間(應力、應變為零處)在圖8 中觀察到了應力、應變的增加.這是由于熱膨脹應力導致了巖石固體發生向外擴張的形變(如圖8(c),當T>400 K,應變量在5～9 μm),導致原本由孔隙占據的空間,變成了固體占據的空間.

固體顆粒受熱變形導致的直接結果便是巖心孔隙度的變化.圖8(d)展示的孔隙度隨著巖心兩端驅替溫差的變化(孔隙度的計算方法如式(2)所示).整體而言,孔隙度隨著溫度的增加而增加,當DT>150 K(入口溫度為450 K 時),孔隙度變化幅度不大.從DT=0～250 K,孔隙度的變化幅度為8%.Sengun[33]認為溫度與孔隙度的變化是正相關的關系,在他的研究中觀察到了當巖心從105 °C 加熱到600 °C 后,孔隙度增加的現象.Zhang 等[34]也在實驗中觀察到了孔隙度隨溫度增加的現象.

溫度導致的微小孔隙結構變化改變了兩相流體的相對滲流特征.如圖9 所示,溫度的改變使得水相的相對滲流能力提高,而油相的相對滲流能力幾乎沒有改變,導致等滲點向原等滲點的左邊移動[35-36].

圖9 不同注入溫度下的歸一化相滲曲線Fig.9 Normalized relative permeability curves at different injecting temperature

4 結論

(1)基于Darcy-Brinkman-Biot 的流固耦合數值方法,使用了Duhamel-Neumann 熱彈性應力準則計算熱應力,實現了一個多孔介質內考慮熱流固耦合作用的兩相流動模型;

(2)熱流固耦合的綜合作用改變孔隙結構參數,影響了兩相流體在多孔介質內的滲流過程.主要體現為孔隙的擴大提高了相對滲流能力,孔隙的減小降低了相對滲流能力;

(3)在注入溫度不變的情況下,隨著注入PV 數的增加,巖石所受熱應力與流固耦合作用產生的應力達到平衡.熱應力與流固耦合作用產生的應力方向相反,使得總應力比單獨考慮流固耦合作用下的應力小;

(4)注入溫度的增加使得模型孔隙度增加,增加幅度為8%.但當注入溫差達到150 K 后,孔隙度不再有明顯增加.溫度的增加改變了孔隙結構,間接導致了兩相滲流規律的變化,使得水相的相對滲流能力隨著溫度的增加而增加,等滲點左移.