熱沖擊作用下花崗巖溫度場(chǎng)分布規(guī)律數(shù)值模擬研究

熊貴明，郤保平,b，吳陽(yáng)春，陳路海

(太原理工大學(xué) a.礦業(yè)工程學(xué)院，b.原位改性采礦教育部重點(diǎn)實(shí)驗(yàn)室，太原 030024)

近年來(lái)，為扎實(shí)推進(jìn)經(jīng)濟(jì)持續(xù)健康發(fā)展，我國(guó)大力發(fā)展高效清潔能源產(chǎn)業(yè)，堅(jiān)持走綠色低碳循環(huán)發(fā)展的新路子。干熱巖是一種清潔安全高效的能源資源，對(duì)干熱巖的開(kāi)發(fā)利用有利于增加清潔能源供給，優(yōu)化能源結(jié)構(gòu)，提高能源綜合利用效率，建立多元供應(yīng)體系，對(duì)保障國(guó)家能源安全有著重要意義。我國(guó)干熱巖資源儲(chǔ)量巨大，初步估算3～10 km深處干熱巖資源總計(jì)為2.5×1025J，相當(dāng)于8.56×1014t標(biāo)煤[1]。特別是干熱巖發(fā)電技術(shù)可大幅度降低溫室效應(yīng)和酸雨對(duì)環(huán)境的影響，且不受季節(jié)、氣候制約，可以有效取代煤炭、石化能源消耗，有效保護(hù)生態(tài)環(huán)境。更為重要的是干熱巖開(kāi)發(fā)成本低，利用干熱巖發(fā)電的成本僅為風(fēng)力發(fā)電的1/2，或太陽(yáng)能發(fā)電的1/10.正因?yàn)楦蔁釒r具有這么多優(yōu)點(diǎn)，干熱巖的開(kāi)發(fā)備受關(guān)注。

自20世紀(jì)70年代初學(xué)者們提出干熱巖開(kāi)發(fā)構(gòu)想以來(lái)，很多國(guó)家都投入了大量人力財(cái)力進(jìn)行研究[2-4]。到目前為止，國(guó)內(nèi)外對(duì)干熱巖的生熱、導(dǎo)熱、控?zé)?、?chǔ)熱、釋熱構(gòu)造和規(guī)律性研究仍然很少[5]。我國(guó)對(duì)干熱巖開(kāi)發(fā)研究較晚，且研究主要集中在勘探、鉆井、壓裂和監(jiān)測(cè)等方面，目前主要存在以下問(wèn)題[6-9]：

1) 在干熱巖鉆井過(guò)程中，如何控制鉆井液導(dǎo)致的井壁圍巖坍塌崩落、縮頸或地層破裂等失穩(wěn)問(wèn)題。

2) 在熱儲(chǔ)層建造過(guò)程中，怎樣通過(guò)水力致裂、熱激發(fā)或化學(xué)激發(fā)產(chǎn)生雜亂無(wú)章、分布均勻、沒(méi)有主流方向的裂縫網(wǎng)絡(luò)等裂縫擴(kuò)展問(wèn)題。

3) 在工質(zhì)與巖石進(jìn)行熱交換提取熱量的過(guò)程中，如何選擇工質(zhì)(水、液氮等)、流速、注液壓力和開(kāi)采年限等參數(shù)，確保經(jīng)濟(jì)利益最大化。

由此可以看出，在干熱巖開(kāi)發(fā)的整個(gè)過(guò)程中，熱沖擊作用會(huì)引起巖體內(nèi)的溫度場(chǎng)重新分布，進(jìn)而產(chǎn)生熱應(yīng)力，促使巖石熱力學(xué)性質(zhì)發(fā)生變化[10-14]。因此研究干熱巖在熱沖擊作用下的溫度場(chǎng)分布規(guī)律是干熱巖開(kāi)發(fā)的基礎(chǔ)。

研究干熱巖開(kāi)發(fā)整個(gè)過(guò)程，必須圍繞巖石在熱沖擊作用下的溫度場(chǎng)分布規(guī)律這一主線進(jìn)行。本文采用Comsol Multiphysics 5.3a軟件，依據(jù)某干熱巖項(xiàng)目人工熱儲(chǔ)層的花崗巖參數(shù)，建立熱沖擊作用下花崗巖的溫度場(chǎng)數(shù)學(xué)物理模型；通過(guò)對(duì)其進(jìn)行數(shù)值模擬，分析熱沖擊作用下花崗巖在不同工質(zhì)(H2O、液態(tài)CO2和液N2)和初始溫度時(shí)的溫度場(chǎng)變化規(guī)律。

1 花崗巖溫度場(chǎng)數(shù)學(xué)模型

熱沖擊和熱破裂兩個(gè)概念的內(nèi)涵如下：

1) 巖石受到急劇的加熱和冷卻時(shí)，其內(nèi)部產(chǎn)生很大的溫差，從而引起很大的沖擊熱應(yīng)力，這種現(xiàn)象稱為熱沖擊。一次大的熱沖擊所產(chǎn)生的熱應(yīng)力能超過(guò)巖石的屈服極限，從而導(dǎo)致巖石失穩(wěn)。

2) 巖石受熱后，由于其內(nèi)部各種造巖礦物顆粒的熱膨脹系數(shù)不同，巖石內(nèi)部產(chǎn)生熱應(yīng)力并產(chǎn)生裂縫；在熱應(yīng)力作用下巖石的內(nèi)部裂縫逐漸擴(kuò)展，裂縫相互連通形成裂縫網(wǎng)絡(luò)，從而引起儲(chǔ)層巖石物理性質(zhì)(孔隙度、滲透率等)的改變。這種由熱應(yīng)力引起的巖石破裂稱為巖石熱破裂。

總之，熱沖擊強(qiáng)調(diào)過(guò)程，而熱破裂強(qiáng)調(diào)結(jié)果。本文采用熱沖擊速度、溫度梯度等概念分別去闡釋巖石表面和內(nèi)部的溫度場(chǎng)分布規(guī)律。其中，熱沖擊速度是指單位時(shí)間內(nèi)溫度變化量，溫度梯度是指單位長(zhǎng)度內(nèi)溫度變化量。

1.1 基本假設(shè)

干熱巖的主要成分是花崗巖，花崗巖由多種成礦物質(zhì)組成。熱沖擊作用時(shí)，巖石內(nèi)部多種物質(zhì)因具有不同的熱效應(yīng)，相互制約、相互影響，因此熱沖擊過(guò)程是一個(gè)相當(dāng)復(fù)雜的過(guò)程。本文假設(shè)花崗巖僅受溫度場(chǎng)的影響。此外，還對(duì)模型作如下假設(shè)：

1) 花崗巖各向同性，不發(fā)生相變。

2) 花崗巖的初始溫度是均勻的。

3) 花崗巖熱沖擊過(guò)程中忽略熱對(duì)流和熱輻射，只考慮熱傳導(dǎo)。

4) 花崗巖在熱沖擊作用下的溫度變化過(guò)程滿足局部熱平衡。

1.2 巖石瞬態(tài)熱傳導(dǎo)模型

根據(jù)傳熱學(xué)理論，由能量守恒定律，得三維圓柱體瞬態(tài)溫度場(chǎng)方程[15]：

(1)

巖石內(nèi)部沒(méi)有熱源，即熱流量Φ=0，且圓柱體的幾何形狀和溫度作用都對(duì)稱于圓柱體的中心軸線。因此這一問(wèn)題屬于空間軸對(duì)稱問(wèn)題，式(1)可簡(jiǎn)化為：

(2)

式中：ρ為巖石密度；c為巖石比熱容；k為巖石導(dǎo)熱系數(shù)；T為巖石瞬時(shí)溫度；τ為時(shí)間；r為巖石試件半徑；z為巖石試件高度。為了簡(jiǎn)化模型，將巖石的熱力學(xué)參數(shù)(彈性模量、熱膨脹系數(shù)和泊松比)取為常數(shù)。

為了得到花崗巖溫度值的唯一結(jié)果，須給出一定的邊界條件，具體如下。

第一類邊界條件：取溫度為定值，即

T1=C.

(3)

第二類邊界條件：給定邊界上為絕熱邊界，即

(4)

第三類邊界條件：給定邊界上物體表面與周圍流體間的表面?zhèn)鳠嵯禂?shù)h及流體的溫度Tf，但不發(fā)生對(duì)流傳熱，即

(5)

綜合式(2)—(5)可得出，熱沖擊作用下花崗巖溫度場(chǎng)的物理模型如下:

(6)

式中：T為巖石邊界溫度；C為常數(shù)；n為巖石換熱表面的外法線的長(zhǎng)度；h為巖石表面?zhèn)鳠嵯禂?shù)；Tw為巖石表面溫度；Tf為流體溫度；r0為巖石的最大半徑；b為巖石最大高度；T0為巖石初始溫度。

2 熱沖擊作用下花崗巖有限元模型的建立

2.1 模型建立

圖1 花崗巖試件數(shù)值計(jì)算物理模型(單位：mm)Fig.1 The granite physical model of numerical calculation

選取某干熱巖項(xiàng)目人工熱儲(chǔ)層花崗巖為研究對(duì)象，對(duì)室內(nèi)實(shí)驗(yàn)進(jìn)行數(shù)值模擬。模擬對(duì)象采用φ50 mm×100 mm的圓柱體花崗巖試件，見(jiàn)圖1。由于該模型所加載的邊界條件均滿足軸對(duì)稱條件，在三維模型中取一截面進(jìn)行研究。在這里規(guī)定，沿圓柱體底面圓的半徑方向?yàn)閺较?，沿圓柱體的高度方向?yàn)檩S向。為了便于研究花崗巖的溫度場(chǎng)分布規(guī)律，在花崗巖試件上布置監(jiān)測(cè)點(diǎn)1—10(見(jiàn)圖1)。監(jiān)測(cè)點(diǎn)1—6相鄰兩點(diǎn)之間距離為5 mm，監(jiān)測(cè)點(diǎn)6—9相鄰兩點(diǎn)之間距離為25 mm.

2.2 參數(shù)選取

模擬中選取的花崗巖的比熱容c=850 J/(kg·K)，密度ρ=2 600 kg/m3，導(dǎo)熱系數(shù)λ=2.9 W/(m·K)，泊松比μ=0.25，彈性模量E=60 GPa，熱膨脹系數(shù)α=7×10-6K-1.

2.3 花崗巖熱沖擊模擬實(shí)驗(yàn)方案

花崗巖試件的初始溫度分別為373.15，473.15，573.15，673.15，773.15，873.15 K；工質(zhì)分別為常溫水(293.15 K)、液態(tài)二氧化碳(236.15 K)和液氮(77.15 K).本模擬方案為：將花崗巖試件加熱到初始溫度后，將沿軸向1/4部分的花崗巖試件浸入不同工質(zhì)；按照表1方案模擬花崗巖試件熱沖擊下溫度的變化情況。其中，底部和其他部分均為絕熱邊界。以圖1截面圖中左下角O為原點(diǎn)，則監(jiān)測(cè)點(diǎn)1—10的坐標(biāo)分別為(0,25)，(5,25)，(10,25)，(15,25)，(20,25)，(25,25)，(25,50)，(25,75)，(25,100)，(0,100).

表1 數(shù)值模擬實(shí)驗(yàn)方案Table 1 Numerical simulation experiment scheme

2.4 網(wǎng)格劃分與求解設(shè)置

網(wǎng)格劃分是數(shù)值計(jì)算的重要一環(huán)，網(wǎng)格劃分的質(zhì)量直接影響整個(gè)模型計(jì)算過(guò)程的精細(xì)程度。本模型為二維軸對(duì)稱，對(duì)稱軸為r=0，采用Comsol Multiphysics 5.3a軟件網(wǎng)格繪制功能中的三角形網(wǎng)格進(jìn)行劃分。

在Comsol Multiphysics 5.3a軟件中添加瞬態(tài)研究，熱沖擊時(shí)間設(shè)置為500 min，時(shí)間間隔設(shè)為0.1 min，使用參數(shù)化掃描，求解花崗巖在熱沖擊作用下的溫度場(chǎng)分布規(guī)律。

3 結(jié)果與討論

3.1 熱沖擊作用下花崗巖溫度場(chǎng)總體變化規(guī)律

在干熱巖開(kāi)發(fā)中，從鉆井、人工熱儲(chǔ)建造到熱能提取整個(gè)過(guò)程均涉及到巖石的溫度場(chǎng) ，因此，研究花崗巖溫度場(chǎng)的變化規(guī)律是干熱巖開(kāi)發(fā)最關(guān)注的信息點(diǎn)。

由圖2可知，伴隨著熱沖擊作用，無(wú)論是哪種工質(zhì)，熱沖擊作用下花崗巖的溫度場(chǎng)均以波的形式呈倒梯形向內(nèi)部以一定的速度傳播；工質(zhì)接觸的面積越大，溫度場(chǎng)傳播速度越快；熱沖擊波傳播所經(jīng)過(guò)的點(diǎn)溫度會(huì)突然降低，直至最后趨于穩(wěn)定。

圖2(a)、(b)和(c)中的溫度場(chǎng)云圖僅能看到溫度場(chǎng)特定時(shí)間的分布云圖，不能明顯地看出溫度隨時(shí)間的變化規(guī)律。

花崗巖初始溫度：873.15 K圖2 不同工質(zhì)下花崗巖熱沖擊作用1 min時(shí)剖面溫度場(chǎng)云圖Fig.2 Temperature field cloud map of 1min time profile of hot shock of granite under different working conditions

在不同工質(zhì)與初始溫度的條件下，熱沖擊作用時(shí)花崗巖的溫度隨時(shí)間的增加而降低。取離沖擊波源最遠(yuǎn)的監(jiān)測(cè)點(diǎn)10為研究對(duì)象，如圖3所示。在所研究的溫度范圍內(nèi)，對(duì)于工質(zhì)一(H2O)而言，當(dāng)初始溫度小于473.15 K時(shí)，花崗巖在1.1 min波及到整塊巖石；當(dāng)初始溫度大于473.15 K時(shí)，花崗巖在1.4 min波及到整塊巖石。此外，初始溫度小于673.15 K時(shí)，花崗巖達(dá)到與工質(zhì)相同溫度的時(shí)間隨著初始溫度的升高而增加；而初始溫度大于673.15 K時(shí)，花崗巖隨初始溫度的升高所需時(shí)間減少；當(dāng)初始溫度大于773.15 K時(shí)，花崗巖不再發(fā)生變化，且保持在236 min.

對(duì)于工質(zhì)二(液態(tài)CO2)而言，當(dāng)初始溫度小于473.15 K時(shí)，花崗巖在1.2 min波及到整塊巖石；當(dāng)初始溫度大于473.15 K時(shí)，花崗巖在1.5 min波及到整塊巖石；同時(shí)，達(dá)到與介質(zhì)相同的溫度的變化規(guī)律同工質(zhì)一；但在初始溫度大于773.15 K時(shí)，花崗巖不再發(fā)生變化且保持在243 min.

圖3 熱沖擊作用下花崗巖溫度隨時(shí)間的變化Fig.3 Temperature change with time under hot shock

對(duì)于工質(zhì)三(液N2)而言，花崗巖在1.5 min波及到整塊巖石，與巖石的初始溫度無(wú)關(guān)；同時(shí)，達(dá)到與介質(zhì)相同的溫度變化規(guī)律同工質(zhì)一；但在初始溫度大于773.15 K時(shí)，花崗巖不再發(fā)生變化且保持在294 min.

由此可見(jiàn)，在不同的工質(zhì)中，花崗巖溫度場(chǎng)的熱沖擊作用在100 s左右波及到整塊花崗巖且變化規(guī)律一致，可見(jiàn)熱沖擊具有時(shí)間效應(yīng)；隨著工質(zhì)溫度的降低，達(dá)到與介質(zhì)相同的溫度的時(shí)間變長(zhǎng)；在初始溫度673.15～773.15 K之間存在一個(gè)溫度閾值，可能會(huì)引起花崗巖的相態(tài)轉(zhuǎn)變，進(jìn)而影響花崗巖溫度場(chǎng)的變化規(guī)律。

由于熱沖擊作用下花崗巖的溫度場(chǎng)分布規(guī)律一致，下面僅選取工況一對(duì)熱沖擊作用下花崗巖表面及內(nèi)部溫度場(chǎng)進(jìn)行討論與分析。

3.2 熱沖擊作用下花崗巖表面溫度場(chǎng)變化規(guī)律

花崗巖表面監(jiān)測(cè)點(diǎn)的熱沖擊速度隨時(shí)間變化的曲線如圖4.可見(jiàn)，熱沖擊作用下花崗巖的熱沖擊速度隨時(shí)間的增加而降低；花崗巖表面在同一時(shí)間任一點(diǎn)的熱沖擊速度相等，表明熱沖擊作用下的花崗巖溫度是以沖擊波的形式向前傳遞，進(jìn)而影響花崗巖的溫度場(chǎng)分布。

花崗巖初始溫度：473.15 K圖4 花崗巖表面熱沖擊速度隨時(shí)間的變化Fig.4 Hot shock velocity of granite surface changes with time

花崗巖表面某點(diǎn)在不同初始溫度下的熱沖擊速度隨時(shí)間變化曲線如圖5所示。由圖5可知，熱沖擊作用下花崗巖溫度在15 s之內(nèi)迅速下降；隨著時(shí)間的增加，熱沖擊速度迅速降低且趨于零?？梢?jiàn)熱沖擊作用具有瞬時(shí)性，且熱沖擊速度隨著花崗巖初始溫度的升高而增大。

圖5 花崗巖表面點(diǎn)7熱沖擊速度隨時(shí)間的變化Fig.5 Hot shock velocity of granite surface point 7 with time

3.3 熱沖擊作用下花崗巖內(nèi)部溫度場(chǎng)變化規(guī)律

取花崗巖內(nèi)部監(jiān)測(cè)點(diǎn)1—5，各監(jiān)測(cè)點(diǎn)溫度梯度隨時(shí)間變化的曲線見(jiàn)圖6.可以看出，熱沖擊作用下花崗巖沿徑向的溫度梯度隨時(shí)間增加呈現(xiàn)先增加后減小，最后趨于零；花崗巖內(nèi)部監(jiān)測(cè)點(diǎn)5、4、3、2、1達(dá)到溫度梯度峰值所用的時(shí)間分別是0，0.8，1.2，3.3，3.6 min.在熱沖擊作用邊界上，花崗巖瞬間溫度梯度達(dá)到相當(dāng)大的值；離熱沖擊作用邊界距離越遠(yuǎn)，沿徑向熱傳導(dǎo)溫度梯度峰值越低且所用的時(shí)間越長(zhǎng)，說(shuō)明熱沖擊作用波沿徑向的傳播過(guò)程會(huì)沿程損失，熱沖擊作用程度隨熱邊界距離增大而減弱。

初始溫度：473.15 K圖6 花崗巖沿徑向溫度梯度隨時(shí)間的變化Fig.6 Temperature gradient of the granite in the radial direction changes with time

由于各監(jiān)測(cè)點(diǎn)溫度場(chǎng)變化規(guī)律一致，選取點(diǎn)1研究在不同初始溫度下花崗巖內(nèi)部溫度場(chǎng)的變化，見(jiàn)圖7.圖7顯示：熱沖擊作用下的花崗巖隨著初始溫度的升高，溫度梯度增大；不同初始溫度的花崗巖熱沖擊溫度梯度隨時(shí)間的增加呈現(xiàn)先增加后減小，最后趨于零。經(jīng)監(jiān)測(cè)，當(dāng)花崗巖的初始溫度為373.15 K時(shí)，其溫度梯度峰值是1 395 K/m；當(dāng)花崗巖初始溫度分別為473.15，573.15，673.15，773.15，873.15 K時(shí)，其溫度梯度峰值分別是373.15 K時(shí)溫度梯度峰值的2.30，1.61，1.35，1.26，1.20倍。這說(shuō)明溫度差越大，熱沖擊作用越強(qiáng)烈，巖石中熱量損失自然越大。

圖7 不同初始溫度下花崗巖點(diǎn)1溫度梯度隨時(shí)間的變化Fig.7 Variation of temperature gradient of granite point 1 with time at different temperatures

綜上所述，花崗巖的熱沖擊溫度變化曲線可分為三個(gè)階段：第一階段為溫度加速變化階段，溫度快速增加，溫度梯度隨時(shí)間的增加呈直線增大到峰值；第二階段為溫度減速變化階段，溫度減速增加，溫度梯度隨時(shí)間的增加由峰值減小到零；第三階段為溫度穩(wěn)定階段，溫度不再變化，溫度梯度為零，達(dá)到熱平衡。

3.4 熱沖擊作用下花崗巖表面溫度實(shí)驗(yàn)與數(shù)值模擬結(jié)果對(duì)比

由于目前實(shí)驗(yàn)設(shè)備的限制，熱沖擊作用下花崗巖的物理實(shí)驗(yàn)只能觀測(cè)花崗巖表面的溫度場(chǎng)規(guī)律。因此，選取工質(zhì)一(H2O)，將數(shù)值模擬所得花崗巖表面在熱沖擊作用下5 s時(shí)熱沖擊速度隨溫度的變化曲線與花崗巖試件在相同條件下實(shí)驗(yàn)所得的曲線進(jìn)行對(duì)比，結(jié)果如圖8所示。從熱沖擊速度隨溫度的變化趨勢(shì)來(lái)看，數(shù)值模擬與物理實(shí)驗(yàn)得到的結(jié)果相差不大。但數(shù)值模擬建立在簡(jiǎn)化的基礎(chǔ)上，未能真實(shí)地反映花崗巖本身性質(zhì)。

圖8 花崗巖表面熱沖擊速度隨初始溫度變化的實(shí)驗(yàn)與模擬結(jié)果對(duì)比Fig.8 Comparison of experimental and simulation results of thermal shock velocity of granite surface with initial temperature

4 結(jié)論

對(duì)某干熱巖項(xiàng)目人工熱儲(chǔ)層的花崗巖進(jìn)行了熱沖擊作用下的溫度場(chǎng)分布規(guī)律研究，得到如下結(jié)論:

1) 花崗巖在熱沖擊作用下溫度場(chǎng)以沖擊波形式傳播，工質(zhì)接觸的面積越大，溫度場(chǎng)傳播速度越快；熱平衡時(shí)間隨著工質(zhì)溫度降低而延長(zhǎng)，且在花崗巖初始溫度673.15～773.15 K之間存在著一個(gè)可能會(huì)引起巖石相態(tài)轉(zhuǎn)變的溫度閾值。

2) 花崗巖初始溫度越高，表面溫度變化越劇烈；而溫度梯度越大，內(nèi)部溫度變化越劇烈，巖石中熱量損失越大。

3) 花崗巖的熱沖擊溫度變化曲線可分為三個(gè)階段：第一階段為溫度加速變化階段，溫度梯度隨時(shí)間的增加呈直線增大到峰值；第二階段為溫度減速變化階段，溫度梯度隨時(shí)間的增加由峰值減小到零；第三階段為溫度穩(wěn)定階段，溫度梯度為零，達(dá)到熱平衡。

4) 巖石是多種物理場(chǎng)共同制約形成的產(chǎn)物，其自身成分也十分復(fù)雜。目前的數(shù)值模擬均是不同程度的簡(jiǎn)化處理，比如，假定巖石為各向同性、不發(fā)生相變的均質(zhì)體。因此，利用數(shù)值模擬研究巖石熱沖擊作用下的溫度場(chǎng)規(guī)律是可行的，但仍然無(wú)法完全替代室內(nèi)及現(xiàn)場(chǎng)試驗(yàn)，只能起到補(bǔ)充作用。

未來(lái)對(duì)花崗巖溫度場(chǎng)的研究，應(yīng)該考慮過(guò)高的溫度對(duì)巖石內(nèi)部形態(tài)和結(jié)構(gòu)的影響。