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

熊貴明，郤保平,b，吳陽春，陳路海

(太原理工大學 a.礦業(yè)工程學院，b.原位改性采礦教育部重點實驗室，太原 030024)

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

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

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

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

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

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

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

1 花崗巖溫度場數(shù)學模型

熱沖擊和熱破裂兩個概念的內涵如下：

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

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

總之，熱沖擊強調過程，而熱破裂強調結果。本文采用熱沖擊速度、溫度梯度等概念分別去闡釋巖石表面和內部的溫度場分布規(guī)律。其中，熱沖擊速度是指單位時間內溫度變化量，溫度梯度是指單位長度內溫度變化量。

1.1 基本假設

干熱巖的主要成分是花崗巖，花崗巖由多種成礦物質組成。熱沖擊作用時，巖石內部多種物質因具有不同的熱效應，相互制約、相互影響，因此熱沖擊過程是一個相當復雜的過程。本文假設花崗巖僅受溫度場的影響。此外，還對模型作如下假設：

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

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

3) 花崗巖熱沖擊過程中忽略熱對流和熱輻射，只考慮熱傳導。

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

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

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

(1)

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

(2)

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

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

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

T1=C.

(3)

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

(4)

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

(5)

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

(6)

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

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

2.1 模型建立

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

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

2.2 參數(shù)選取

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

2.3 花崗巖熱沖擊模擬實驗方案

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

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

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

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

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

3 結果與討論

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.3 熱沖擊作用下花崗巖內部溫度場變化規(guī)律

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

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

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

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

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

3.4 熱沖擊作用下花崗巖表面溫度實驗與數(shù)值模擬結果對比

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

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

4 結論

對某干熱巖項目人工熱儲層的花崗巖進行了熱沖擊作用下的溫度場分布規(guī)律研究，得到如下結論:

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

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

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

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

未來對花崗巖溫度場的研究，應該考慮過高的溫度對巖石內部形態(tài)和結構的影響。