蔡少斌 楊永飛 劉 杰

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

引言

經(jīng)過(guò)多年的勘探與開(kāi)發(fā),我國(guó)淺層和中淺層的油氣資源已經(jīng)進(jìn)入了較為成熟的開(kāi)發(fā)階段.伴隨著勘探開(kāi)發(fā)的深入,人們認(rèn)識(shí)到在深層地層中同樣孕育著豐富的可采能源,例如深層油氣藏、地?zé)豳Y源等.深層能源儲(chǔ)藏具有可觀的挖潛能力,是未來(lái)助力我國(guó)增儲(chǔ)上產(chǎn)的重要來(lái)源[1].

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

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

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

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

綜上,本工作采用Darcy-Brinkman-Biot 方法,耦合了Duhamel-Neumann 熱彈性應(yīng)力準(zhǔn)則,推導(dǎo)了考慮熱流固三場(chǎng)耦合作用下的多孔介質(zhì)兩相流動(dòng)模擬模型,研究了多場(chǎng)耦合作用對(duì)油水兩相滲流的影響,以期為油氣田開(kāi)發(fā)工程提供參考數(shù)據(jù).

1 數(shù)學(xué)模型

假設(shè)模型為飽和流體的多孔介質(zhì),其中油相Vo、水相Vw及固體相Vs分別占據(jù)模型的不同空間位置,它們之間的關(guān)系為

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

油相與水相統(tǒng)稱為流體相,以相分?jǐn)?shù) φf(shuō)表示;固體相以固體體積分?jǐn)?shù) φs表示(如圖1 所示),二者的計(jì)算式分別為

顯然, φf(shuō)與 φs兩者之間具有以下關(guān)系

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

顯然, αw與 αo之間具有如下關(guān)系

因此,通過(guò)圖1 中所定義的參數(shù),可以判斷油相、水相及固體相在模型中的分布,其判定標(biāo)準(zhǔn)為

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

對(duì)于不可壓縮、不可混溶的兩相牛頓流體,可以采用流體體積方法[24]描述多相流動(dòng)過(guò)程.采用Carrillo 等[21]基于體積平均方法推廣的流體體積方法公式來(lái)描述模型內(nèi)的兩相流動(dòng)過(guò)程,即

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

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

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

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

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

1.3 線性熱彈性應(yīng)力控制方程

在流固耦合的基礎(chǔ)上,進(jìn)一步考慮巖石骨架在熱場(chǎng)作用下產(chǎn)生的熱彈性應(yīng)力對(duì)兩相流動(dòng)過(guò)程的影響.

傳熱過(guò)程中的動(dòng)能方程為

流體換熱過(guò)程中的熱傳導(dǎo)方程為

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

溫度對(duì)骨架產(chǎn)生的熱應(yīng)力為

其中, σt為熱應(yīng)力,N/m2;T為溫度,K;β 為熱應(yīng)力系數(shù).

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

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

2 模型與參數(shù)

2.1 模型參數(shù)

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

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

模型模擬的過(guò)程為水相驅(qū)替油相的過(guò)程,以速度邊界條件作為控制條件,在入口處以0.001 m/s 的注入流速注入水相.模型入口溫度設(shè)置為350 K,出口溫度設(shè)置為300 K,模型內(nèi)部初始溫度場(chǎng)設(shè)置為300 K.模型網(wǎng)格采用六面體結(jié)構(gòu)化網(wǎng)格,網(wǎng)格數(shù)為10000 個(gè).

2.2 模擬參數(shù)

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

表1 模擬參數(shù)設(shè)置Table 1 Simulation parameters

2.3 模型求解

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

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

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

式中,Ux和Uy分別代表了速度在x,y方向上的分量,m/s;Δt表示求解時(shí)間步長(zhǎng),s;Δx和 Δy分別代表了x,y方向上的求解空間長(zhǎng)度,m;C 代表判定標(biāo)準(zhǔn),為常數(shù),在本研究中,經(jīng)過(guò)反復(fù)測(cè)試,取C = 0.15 模型能夠穩(wěn)定計(jì)算.

3 結(jié)果與討論

3.1 模擬結(jié)果

圖4 展示的是恒定流速0.001 m/s 為注入流速,水相潤(rùn)濕角設(shè)置為45°,入口溫度為350 K,初始內(nèi)部溫度與出口溫度為300 K 的模型的模擬結(jié)果.從圖中可以看出,水相在巖心切片的右側(cè)形成了優(yōu)勢(shì)主力流動(dòng)通道,驅(qū)替前端在注入速度的作用下發(fā)生了潤(rùn)濕滯后.隨著驅(qū)替過(guò)程的進(jìn)行,水相不斷地將油相驅(qū)替出巖心切片,最終當(dāng)水相到達(dá)出口時(shí),主力流動(dòng)通道達(dá)到穩(wěn)定狀態(tài),模型內(nèi)的兩相飽和度趨于不變.

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

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

圖5 0.15 s 時(shí)應(yīng)力及應(yīng)變模擬結(jié)果Fig.5 Simulation results of (a) stress and (b) strain profile at 0.15 s

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

圖6 歸一化相滲曲線與實(shí)驗(yàn)結(jié)果對(duì)比Fig.6 Normalized relative permeability curve vs.experimental result

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

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

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

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

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

3.2 注入PV 數(shù)的影響

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

保持圖4 中的模擬條件不變,圖7 展示的是PV 數(shù)為5～160 之間的模擬結(jié)果.圖7(a)～圖7(c)分別展示的是模型沿著驅(qū)替方向中心軸線的溫度、應(yīng)力及應(yīng)變分布(圖中應(yīng)力、應(yīng)變?yōu)榱愕牟糠譃榭紫犊臻g).

圖7 不同注入PV 數(shù)的模擬結(jié)果Fig.7 Simulation results under different PV numbers

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

3.3 注入溫差的影響

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

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

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

隨著注入溫差的不斷增加,模型內(nèi)的溫度梯度不斷增加(如圖8(a)所示).此時(shí),熱膨脹應(yīng)力成為了造成巖石彈性變形的主導(dǎo)力.如圖8(b)和圖8(c)所示,隨著溫度的增加,巖石固體所受的應(yīng)力、應(yīng)變隨之增加.與圖7(b)和圖7(c)相比,當(dāng)模型的注入溫度達(dá)到400 K時(shí),熱膨脹應(yīng)力已經(jīng)抵消了流體與骨架之間的流固耦合相應(yīng),此時(shí)流體對(duì)骨架產(chǎn)生的應(yīng)力已經(jīng)不足以抵消熱膨脹應(yīng)力,因此巖石的應(yīng)力、應(yīng)變隨著溫度的升高而增加.同時(shí),對(duì)比圖7(b)、圖7(c)與圖8(b)、圖8(c)可以發(fā)現(xiàn),圖7 中原本孔隙占據(jù)的空間(應(yīng)力、應(yīng)變?yōu)榱闾?在圖8 中觀察到了應(yīng)力、應(yīng)變的增加.這是由于熱膨脹應(yīng)力導(dǎo)致了巖石固體發(fā)生向外擴(kuò)張的形變(如圖8(c),當(dāng)T>400 K,應(yīng)變量在5～9 μm),導(dǎo)致原本由孔隙占據(jù)的空間,變成了固體占據(jù)的空間.

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

溫度導(dǎo)致的微小孔隙結(jié)構(gòu)變化改變了兩相流體的相對(duì)滲流特征.如圖9 所示,溫度的改變使得水相的相對(duì)滲流能力提高,而油相的相對(duì)滲流能力幾乎沒(méi)有改變,導(dǎo)致等滲點(diǎn)向原等滲點(diǎn)的左邊移動(dòng)[35-36].

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

4 結(jié)論

(1)基于Darcy-Brinkman-Biot 的流固耦合數(shù)值方法,使用了Duhamel-Neumann 熱彈性應(yīng)力準(zhǔn)則計(jì)算熱應(yīng)力,實(shí)現(xiàn)了一個(gè)多孔介質(zhì)內(nèi)考慮熱流固耦合作用的兩相流動(dòng)模型;

(2)熱流固耦合的綜合作用改變孔隙結(jié)構(gòu)參數(shù),影響了兩相流體在多孔介質(zhì)內(nèi)的滲流過(guò)程.主要體現(xiàn)為孔隙的擴(kuò)大提高了相對(duì)滲流能力,孔隙的減小降低了相對(duì)滲流能力;

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

(4)注入溫度的增加使得模型孔隙度增加,增加幅度為8%.但當(dāng)注入溫差達(dá)到150 K 后,孔隙度不再有明顯增加.溫度的增加改變了孔隙結(jié)構(gòu),間接導(dǎo)致了兩相滲流規(guī)律的變化,使得水相的相對(duì)滲流能力隨著溫度的增加而增加,等滲點(diǎn)左移.