能源樁傳熱數學新模型及其在熱響應測試的適用性

閆振國 ，曾姝，楊軍

(1.清華大學 土木工程系，北京，100080；2.海軍研究院，北京，100071)

以能源樁為代表的地下能源結構可以充分利用結構內部空間布設循環回路，無需支出額外的鉆孔費用，其使用量逐年增加、應用場景不斷增多[1]。在此過程中，土為地熱能儲存、傳遞的載體，準確掌握其傳熱能力是能源樁設計的關鍵，也是提高能源樁傳熱能力[2-3]的基礎。

目前，獲取土體熱物性參數的方法主要包括查表，探針、平板測試以及熱響應測試[4]。

1) 查表法主要根據勘察報告提供的土體類別查詢得到熱物性參數的取值范圍。

2) 探針、平板測試通常在試驗室內進行，通過測試儀器對試樣的檢測，得到熱物性參數結果。試樣在取樣、運輸過程中可能產生擾動，溫度、含水率可能發生變化，對結果的準確性產生影響。

3) 熱響應測試通過工質在地熱交換器內部管道循環的過程與周圍土發生熱交換，根據進出口水溫及/或土體溫度變化反演熱物性參數[5]。

熱物性參數反演時，通過反演方法不斷調整識別模型中參數的取值并進行計算，當計算結果與記錄值之間的誤差滿足設定，即為反演結果。模型與實際傳熱情況的匹配程度、參數反演方法決定了反演參數的準確程度和計算時間[6]。其中，識別模型主要分為數學模型和數值模型。

1) 數值模型可精細地還原整個傳熱過程[7]，隨著計算機科技的進步，模型精確度不斷提高、計算時間不斷減少，但模型建立過程中所需參數多、普適性較差；

2) 數學模型所需參數少，隨著模型的發展，計算結果的準確度不斷提高。

建立數學模型時，假定能源樁傳熱過程中僅發生熱傳導[8]，該形式是能源樁由樁內回路向樁周土傳熱過程中的主要傳熱形式[9-12]，基于傅里葉熱傳導定律，人們提出線、圓柱面等不同幾何構型熱源的數學模型[13]。

數學模型用于土體熱物性參數的反演，始于鉆孔熱交換器的熱響應測試。當鉆孔熱交換器傳熱進入準穩態[13]，其自身傳熱可用熱阻表征，外側傳熱可使用線熱源模型(LSM)描述[14]。利用線熱源模型過余溫度與傳熱時間的對數呈線性相關的特點，便捷地反演土的綜合導熱系數[15-17]。當上述方法用于能源樁時，LSM 無法描述熱源半徑，樁身達到準穩態所需時間較長。空心圓柱面熱源模型(HCSM)能夠考慮熱源半徑，但模型的等功率傳熱假定無法準確描述達到準穩態傳熱前的動態傳熱過程[18]。

為克服上述問題，MAN等[19]提出熱量向熱源兩側同時傳遞的無限長實心圓柱面熱源模型(SCSM)，適用于描述能源樁的傳熱[20]。但是，該模型無法考慮樁土的熱物性差異，反演的熱物性參數是能源樁與樁周土的等效值，當兩者的熱物性差異較大時，結果與實際值相差較大。現有數學模型中，LI等[21-22]分別根據復合介質的線熱源模型和傳熱過程中的能量守恒提出復合介質的圓柱面熱源模型(ICCSM)，該模型充分考慮熱源兩側的熱物性差異，已成功運用于能源樁的熱響應測試[23-24]。實際工程中，循環回路通常綁扎在鋼筋籠的內側，根據CECINATO等[9]的研究結論，當樁周土為砂土時，增加能源樁保護層厚度，即增加熱源外側部分區域的導熱系數，可有效提高傳熱效率，現有數學模型對此鮮有涉及。近年來，能源樁用于大型LNG 儲罐、大型地下空間等結構時，保護層厚度、鋼筋直徑往往較大，有必要建立增加熱源外側熱物性差異的數學模型。

基于上述現狀，本文充分考慮熱源外側混凝土(包括但不限于保護層)與土的熱物性差異，首先，提出圓柱面熱源位于復合介質內層的傳熱模型(ILCCSM)，并求得顯式解析解；其次，以擬建能源樁作為分析對象，使用COMSOL 多場耦合有限元軟件模擬能源樁的傳熱過程，并將樁表面沿長度方向的過余溫度平均值作為記錄數據；最后，分別采用本文模型及LSM，SCSM 和ICCSM，通過模式搜索法反演土的導熱系數和熱擴散系數，并與模型輸入值比較，分析各模型反演結果的準確度，評價模型的適用性。

1 圓柱面熱源位于復合介質內層的傳熱模型

圓柱面熱源位于復合介質內層的傳熱模型(ILCCSM)如圖1 所示，整個傳熱空間由熱源位置r1和復合介質界面位置r2劃分為1，2和3區域，區域內的溫度和分配功率分別為Ti和qi，i=1，2和3。復合介質界面內側(區域1 和2)的熱物性參數以下角標i表示，復合介質界面外側(區域3)的熱物性參數以下角標o表示。

圖1 復合介質內層的圓柱面熱源模型示意圖Fig.1 Schematic diagram of layout of a cylindrical heat source at inner layer of composite medium

模型建立過程中遵從以下假定[22]：

1) 雙層復合介質均為均勻常物性材料；

2) 溫度場的初始溫度T0；

3) 熱源與兩側介質緊密接觸。

傳熱空間各處的過余溫度θi=Ti(t)-T0，滿足如下熱傳導方程：

式中：θi為過余溫度，℃；T0為初始溫度，℃；qi為傳熱功率，W/m；aj為熱擴散系數，m2/s；t為計算傳熱時間，s；r為計算傳熱位置，m；i=1，2，3；j=i，o。空間各處初始溫度為T0，

基于單樁傳熱假定，傳熱過程中圓心位置處為絕熱狀態，見式(3)；傳熱影響有限，見式(4)。

在r1(r1≤r2)處，熱源以恒定功率ql向兩側傳熱，根據能量守恒原理可得式(5)；且在該位置處，2個區域溫度相等，見式(6)。

式中：kj為導熱系數，W·m-1·℃-1。

在r2界面位置處，兩側空間的傳熱量和溫度均相等，則有：

式(1)～(8)完整地描述了傳熱空間的溫度場分布。

采用如圖2 所示的Laplace 變換及逆變換流程對上述方程組進行求解，具體r2外側樁周土的過余溫度如式(9)～(11)所示。

圖2 控制方程求解流程Fig.2 Solving flow of governing equation

式中：J0和J1分別為0 階和1 階第一類貝塞爾函數；Y0和Y1分別為0階和1階第二類貝塞爾函數。

當r1=r2時，式(9)～(11)的解與文獻[22]中熱源外側過余溫度解析解相同，證明本模型解析解的正確性。

2 能源樁熱響應測試數值模型

能源樁進行熱響應測試時，恒定功率的熱量經由回路內的流動工質、傳熱性能較佳的螺旋形循環回路[25]以及樁身混凝土與樁周土發生熱交換，數值模型需對上述傳熱過程進行準確描述。與此同時，合理的地溫分布、符合實際的邊界條件也必不可少。

2.1 模型控制方程

數值模型通過COMSOL 多物理場耦合有限元軟件建立。根據能源樁的傳熱特性，分別使用管道傳熱溫度場和固體傳熱溫度場[24]。其中，根據循環回路管道長徑比較大的特征，使用線單元描述循環工質的瞬態流動和傳熱，包括循環工質的連續性方程、動量方程以及循環工質的熱對流、與管壁發生熱傳導過程中的能量守恒方程，如式(12)～(14)所示。

式中：A為回路截面面積，m2；ρw為循環工質的密度，kg/m3；u為循環工質流速，m/s；p為水壓，Pa；fD為達西摩擦因數；dh為回路水力內徑，m；cw為循環工質的比熱容，J·m-1·℃-1；T為計算溫度，℃；kw為循環工質的導熱系數，W·m-1·℃-1；Qwall為通過管壁的傳熱功率，W。

固體傳熱溫度場包括樁身和周圍土體的熱傳導方程，如式(15)所示。

式中：ρs為固體的密度，kg/m3；cs為固體的比熱容，J·m-1·℃-1；Ts為固體溫度，℃；ks為固體的導熱系數，W·m-1·℃-1。

2.2 模型建立

基于上述控制方程建立數值模型，如圖3 所示。其中，能源樁樁長為20 m，樁徑為1 m；樁身內螺旋形回路直徑為0.8 m，匝距為0.3 m；回路管道內徑為25 mm，壁厚為2.3 mm。樁內循環工質、回路和樁身的材質分別為水、聚乙烯和混凝土。樁周土采用黏土，其熱物性參數隨飽和度基本呈線性增加，通過對能源樁在干黏土和飽和黏土中傳熱的結果進行內插，即可得到能源樁在不同飽和度的黏土中傳熱的結果[21-22]。上述傳熱介質的熱物性參數詳見表1。

圖3 有限元模型示意圖Fig.3 Schematic diagram of layout of finite element model

表1 傳熱介質的熱物性參數Table 1 Thermal properties of media

熱響應測試的傳熱功率可表示為回路進出口水溫的函數，如式(16)所示。

式中：Tin和Tout分別為能源樁的進口溫度和出口溫度，℃。

模擬過程中，選取能源樁的傳熱功率為 100 W/m[26]，則循環工質的流速為0.5 L/s，進出口溫差為1 ℃。

試驗驗證結果表明，在適宜的初始條件和邊界條件下，上述模型能夠準確模擬包括熱響應測試在內的能源樁傳熱問題。

2.2.1 地溫分布

地溫隨深度發生改變，在地表以下數米范圍內受到以年為周期的大氣溫度影響，并體現出一定的滯后性；達一定深度后地溫基本保持恒定。具體可表示為式(17)所示的函數[1]。

式中：Tm為年平均氣溫，℃；Ta為年溫度變化幅值，℃；z為地層深度，m；ω為溫度變化頻率，為2π/365 d-1；as為土壤熱擴散系數，m2/s；t′為時間，d；t0為最高氣溫出現的時間，d。

根據雄安地區的氣象統計數據，土體各月的溫度梯度計算結果如圖4所示。選擇6月作為熱響應測試的時間，模型初始溫度使用該月的地熱溫度曲線。

圖4 地溫分布Fig.4 Distribution of ground temperature

2.2.2 邊界條件

數值模型的邊界條件是制約和影響計算準確度和計算時長的因素之一，在與傳熱實際相符的前提下，既要滿足計算精度，又要考慮計算時長。實際傳熱時，土體受影響區域面積隨傳熱時間增加，但影響范圍有限，定義該影響范圍為模型的計算空間。空間邊界除上表面外，其他邊界均定義為絕熱。對于上表面，日內環境溫度變化對地溫的影響僅限于地表以下1 m范圍內[27]；能源樁外部與熱響應測試儀連接的回路由隔熱棉包裹，受環境溫度影響較小，因此定義上表面為等溫面，溫度為月平均氣溫。使用上述邊界條件，根據熱響應測試時間及相應計算結果對有限元網格進行敏感性分析，確定模型的幾何尺寸，最終確定數值模型。

2.2.3 敏感性分析

模型計算時，以出口水溫作為監測量進行網格敏感性分析。根據不同的網格劃分，單元數量為59 662～700 716 個。出口水溫隨網格尺寸減小(單元數增加)而發生變化，當相對誤差在1%以內時，可忽略網格對結果的影響[24]。經過1 h傳熱的計算，出口水溫的相對誤差均在1%以內，但計算時間從112 s 增加至8 988 s。其中，單元數量為 132 113個時，計算時間213 s，此時整個模型網格非均勻劃分，溫度梯度較大區域較密集，溫度梯度較小區域較粗糙，如圖5所示。后續使用該網格進行計算。

圖5 模型網格劃分Fig.5 Model meshing

確定模型網格劃分形式后，對計算區域進行依賴性測試。分別選用10，12和14 m作為模型正方形橫截面的邊長。經計算，7 d 傳熱期內，采用上述3種邊長均滿足傳熱要求，后續研究模型橫截面邊長選用10 m。

2.3 計算結果

使用上述模型，分別模擬能源樁在飽和黏土、干黏土中為期7 d的熱響應測試，樁表面沿深度方向的平均溫度如圖6所示。相比飽和黏土，能源樁在干黏土中傳熱時樁表面的平均溫度明顯增大。

圖6 樁表面沿深度方向平均溫度Fig.6 Mean temperature of pile's surface along depth direction

3 參數優化方法

熱物性參數反演過程中，使用數學模型試算熱物性參數不同時的溫度，當與記錄數據的相對誤差滿足精度要求時，即作為反演結果。除使用數學模型外，適宜的優化算法是形成準確、高效反演熱物性參數優化方法的有效保障。

3.1 優化算法

優化算法決定了參數反演過程中的迭代次數和結果精度。當目標函數可導時，根據目標函數的一階導數(梯度)、二階導數(Hesse矩陣)或一階偏導數矩陣(Jacob矩陣)，優化收斂方向和收斂步長，使迭代得到最優值。當不確定目標函數是否可導時，可采用坐標輪換法/全局搜索法直接搜索最優方向和最優步長，最終迭代得到反演結果。為提高坐標輪換法/全局搜索法的優化速度，沿相鄰幾點連線的方向進行迭代搜索，直到滿足誤差要求，即本文使用的模式搜索法。該方法通過設定模式搜索法的初始網格尺寸、最大迭代數等優化準則，明確參數的初始值和取值范圍，最終得到取值范圍內的最優反演結果。

3.2 參數靈敏度分析

參數的靈敏度是表征其可識別性和結果不確定性的指標。本文借鑒相關研究結論，采用相對靈敏度系數的概念，其表達式[23]如下：

其中：β為模型中的任意變量。

能源樁傳熱過程中，熱物性參數的靈敏度較高，表明熱物性參數的微小變化導致較大的溫度變化；而靈敏度較小則說明溫度對該參數的變化不敏感，應避免在靈敏度較低階段進行參數反演。

混凝土能源樁分別在飽和黏土以及干黏土中傳熱時，不同數學模型計算樁表面過余溫度的相對靈敏度如圖7所示。由圖7可見：各模型計算的相對靈敏度變化規律基本相同，k2和a2的相對靈敏度分別為負值和正值，絕對值隨傳熱時間增加，LSM 由于無法體現熱源半徑對傳熱的影響，在傳熱初期相對靈敏度為0。對于飽和黏土，SCSM 計算k2的相對靈敏度最大，LSM 在48 h 后，計算k2的相對靈敏度僅比SCSM 的小，而ICCSM 和ILCCSM計算k2的相對靈敏度始終相近；相同傳熱時間下，由SCSM、ICCSM 和ILCCSM 計算a2的相對靈敏度結果依次減小，差值逐漸增加，并趨于穩定，LSM 在傳熱24 h 后數值最大。對于干黏土，各模型計算的相對靈敏度變化規律與飽和黏土相似，計算值較大。根據模型能否考慮復合介質的熱物性差異，計算結果明顯分為2類，忽略熱物性差異的模型各參數相對靈敏度較高。研究過程中，為保證熱物性參數反演值的準確性，選擇1～7 d 傳熱時間的溫度變化作為熱響應測試數據。將計算過余溫度與圖6所示過余溫度的均方根誤差作為目標函數，函數值為0.2 ℃，此時，對于飽和黏土和干黏土，相對誤差分別在10%和5%以內。

圖7 不同模型計算樁表面過余溫度時土體熱物性參數的相對靈敏度Fig.7 Relative sensitivity of soil thermal parameters when calculating pile surface excess temperature by different models

3.3 參數反演過程

利用優化算法不斷改變數學模型中的參數，完成土體熱物性參數反演的過程，如圖8所示。首先，為待反演參數k2和a2賦初始值，分別使用LSM，SCSM，ICCSM和ILCCSM計算樁身過余溫度，通過模式搜索法迭代，當目標函數滿足誤差要求時，數學模型所用熱物性參數即為反演結果。

圖8 熱物性參數反演流程Fig.8 Process of thermal parameters inversion

4 參數反演結果

混凝土能源樁在飽和黏土中傳熱時，模型輸入值和不同數學模型的反演值如圖9(a)所示，兩者的相對誤差如圖10(a)所示。根據與模型輸入值的相對位置，不同模型的反演值可分為兩類。其中，SCSM，ICCSM和ILCCSM反演的a2與模型輸入值基本相等，相對誤差在1%以內；k2與模型輸入值的差值逐漸減小，SCSM 的相對誤差較大，為6.2%，ILCCSM 的相對誤差最小，為4.4%，均在10%以內。LSM 的反演結果與模型輸入值相差較大，k2和a2的相對誤差分別為10.4%和56.0%。經計算，LSM 反演的熱物性參數相對誤差達到10%以內需要25 d。

圖10 熱物性參數反演值與模型輸入值的相對誤差Fig.10 Relative error between calculated value and model input value

與線熱源模型相比，采用向兩側同時傳熱的圓柱面熱源模型均可準確地反演飽和黏土的熱物性參數。其中，考慮熱源兩側的熱物性差異及保護層厚度，可進一步提高熱物性參數反演的準確度，但提高程度有限。

混凝土能源樁在干黏土中傳熱時，模型輸入值和不同數學模型的反演值及其相對誤差分別如圖9(b)和圖10(b)所示。不同模型對干黏土熱物性的反演結果分布較分散。對于a2，LSM 和ICCSM的反演結果相對誤差較大，均大于20%；ILCCSM和SCSM 的計算結果相對誤差均在10%以內，其中ILCCSM的相對誤差僅為2.7%。對于k2，SCSM的反演結果的相對誤差近200%；LSM反演結果的相對誤差為44.0%；ILCCSM和ICCSM反演結果的相對誤差較小，均在10% 以內。綜合考慮，ILCCSM的反演結果相較于其他模型更準確。

各數學模型反演干黏土熱物性參數的結果表明：熱源幾何構型仍影響反演結果的準確性，但熱源兩側的熱物性差異以及在此基礎上保護層厚度的存在對反演結果的準確性影響較大。根據表1，混凝土的導熱系數較飽和黏土和干黏土分別提高17.65%和300.00%；熱擴散系數較飽和黏土提高28.56%，而較干黏土減小25.18%，混凝土與干黏土的熱物性差異遠大于飽和黏土，對熱物性參數反演的準確性影響也更大。

綜合上述分析，在利用熱響應測試反演熱物性參數的過程中，能源樁與樁周土的熱物性差異較小，采用熱量向圓柱面熱源兩側同時傳遞的數學模型即可得到較準確的結果；而熱物性差異較大時，ILCCSM 反演的熱物性參數較LSM、SCSM、ICCSM準確。實際測試過程中，由于樁周土的熱物性參數未知、熱源外混凝土層可能較厚，統一采用ILCCSM反演其熱物性參數更為適宜。

5 結論

1) 提出了圓柱面熱源位于復合介質內層的傳熱模型，在充分考慮熱源兩側熱物性差異的基礎上，進一步增設熱源外側傳熱介質的差異，更貼近能源樁需通過保護層與樁周土發生熱交換的實際情況。

2) 能源樁樁徑對熱物性參數反演的影響較大，在為期7 d的熱響應測試中，線熱源模型反演的熱物性參數相對誤差始終較大。樁-土的熱物性差異影響樁周土熱物性參數的反演結果，當差異較大時，使用本文數學模型反演得到的熱物性參數較線熱源及其他圓柱面熱源模型準確。

3) 使用本文數學模型反演土的熱物性參數，受樁徑、樁-土熱物性差異影響小，所需熱響應測試時間短，適用范圍較其他模型廣。

4) 本研究主要圍繞能源樁傳熱過程的熱傳導提出數學模型，未能考慮土體中可能存在的滲流對傳熱的影響；反演過程中，模式搜索法作為優化方法，反演結果的準確程度與定義的終止條件有關。上述2個方面均制約反演樁周土熱物性參數的準確性和本文方法的繼續推廣，擬在后續工作中重點研究。