能源樁全生命周期熱響應(yīng)半徑簡化計算方法

王哲，劉耶軍，張正威，翁凱文，鄭秀玲，許四法

(1.浙江工業(yè)大學(xué)建筑工程學(xué)院，浙江杭州，310014；2.浙江農(nóng)林大學(xué)風(fēng)景園林與建筑學(xué)院，浙江杭州，311300)

隨著能源壓力的增加，能源樁技術(shù)迅速發(fā)展。能源樁是將熱交換管預(yù)先埋置在建筑樁基礎(chǔ)中，并通過管內(nèi)循環(huán)流體與淺層地?zé)崮苓M(jìn)行熱交換。與傳統(tǒng)的地源熱泵技術(shù)相比，具有占地空間小、不需要額外的鉆孔費用等優(yōu)點，因而，在實際工程中逐漸得到廣泛應(yīng)用。鑒于傳統(tǒng)能源樁傳熱特性問題，國內(nèi)外學(xué)者提出了多種地下?lián)Q熱器傳熱解析模型[1-6]，然而，如何根據(jù)現(xiàn)場地質(zhì)條件確定能源樁埋管換熱器間距一直是困擾工程師的主要問題之一。目前，有關(guān)能源樁埋管換熱器熱響應(yīng)半徑的研究較少，且考慮到我國幅員遼闊，不同熱工分區(qū)的巖土體熱物性參數(shù)、系統(tǒng)規(guī)模程度和使用方式也會有較大差異，因此，常常給工程設(shè)計帶來非常大的不確定性。國內(nèi)外關(guān)于地下?lián)Q熱器布置間距取值的相關(guān)規(guī)范也很不相同，例如瑞典為10 m，奧地利為2.5 m[7]，中國“地源熱泵系統(tǒng)工程技術(shù)規(guī)范”[8]中布置間距取值為3～6 m，但都沒有給出相應(yīng)的計算方法。HART 等[9]在研究單孔埋管地下?lián)Q熱器傳熱特性時，根據(jù)遠(yuǎn)場邊界處巖土體溫度等于初始溫度這一條件計算出遠(yuǎn)場邊界熱響應(yīng)半徑r∞。對于工程中常用的群管換熱器，因為巖土體溫度隨著遠(yuǎn)離熱源呈指數(shù)衰減，所以，根據(jù)r∞布置來安排鉆孔將會浪費大量的場地。SIGNORELLI等[10]通過算例計算得出豎直埋管地下?lián)Q熱器熱響應(yīng)半徑約為15 m。王艷等[11]采用數(shù)值仿真軟件對豎直埋管換熱器熱響應(yīng)半徑的影響因素進(jìn)行了模擬分析，但也沒給出相應(yīng)的結(jié)論或計算方法。王美燕等[12]提出了一種運(yùn)用量綱—參數(shù)曲線的豎直埋管換熱器熱響應(yīng)半徑計算方法，但該方法要對建筑物冷熱負(fù)荷平衡問題進(jìn)行研究，對工程應(yīng)用有一定的難度。李曉星等[13]提出了一種地源熱泵地下?lián)Q熱器熱響應(yīng)半徑計算方法，但所研究的鉆孔孔徑在100～200 mm 之間，不太適合計算大孔徑能源樁埋管換熱器的熱響應(yīng)半徑。張正威等[14]在不考慮鉆孔孔徑影響的基礎(chǔ)上提出了地源熱泵埋管換熱器熱響應(yīng)半徑簡化計算方法，該方法考慮了巖土體的熱物性參數(shù)、地源熱泵系統(tǒng)的運(yùn)行時間以及現(xiàn)場鉆孔布置方式等的影響，但也僅適用于傳統(tǒng)地源熱泵系統(tǒng)的運(yùn)行方式。綜上所述，人們對現(xiàn)有的能源樁埋管換熱器熱響應(yīng)半徑研究較少，缺乏實際的工程應(yīng)用依據(jù)，大多局限在地源熱泵豎直埋管換熱器熱響應(yīng)半徑的研究。趙石嬈等[15]的研究表明，豎直埋管換熱器的軸向長度比徑向長度差幾個數(shù)量級，現(xiàn)有傳熱解析模型在能源樁和地源熱泵系統(tǒng)使用時間內(nèi)溫度場計算結(jié)果沒有明顯差異。為此，本文作者采用比較簡單的無限長線熱源傳熱解析模型，在不考慮熱—力相互耦合作用下[16-18]，分析樁周巖土體過余溫度場的分布特征，在大量計算和分析的基礎(chǔ)上，利用最小二乘法擬合出單工況荷載作用下能源樁熱響應(yīng)半徑計算公式。通過數(shù)值模擬方法研究變負(fù)荷作用對能源樁熱響應(yīng)半徑的影響，并提出能源樁全生命周期熱響應(yīng)半徑簡化計算方法。

1 長時間作用下能源樁傳熱特性

國內(nèi)外許多學(xué)者對豎直埋管地下?lián)Q熱器的傳熱特性進(jìn)行了研究，并提出了多種傳熱解析模型，一般可分為線熱源模型、面熱源模型和體熱源模型。假設(shè)地下巖土的導(dǎo)熱系數(shù)為1.2 W/(m·K)，熱擴(kuò)散系數(shù)為10-6m2/s，能源樁每延米放熱功率為50 W/m，樁徑為0.6 m，樁長為30 m。圖1所示為不同解析模型樁壁處過余溫度隨時間的變化。從圖1可以看出：長時間作用下幾種代表性解析模型的計算結(jié)果差異很小。為此，本文采用相對簡單的無限長線熱源模型[1]：

式中：θ(r,τ)為過余溫度，oC；r為圓柱坐標(biāo)系中的徑向坐標(biāo)，m；a為地下巖土體的熱擴(kuò)散系數(shù)，m2/s；τ為時間，s；ql為每延米換熱功率，W/m；k為巖土體的導(dǎo)熱系數(shù)，W/(m·K)；Ei(x)為指數(shù)積分函數(shù)。

圖1 不同解析模型中過余溫度計算結(jié)果對比Fig.1 Comparison of excess temperature calculation results of different analytical models

圖2所示為不同計算時間下能源樁樁周巖土體過余溫度的徑向分布圖。從圖2可以看出：樁壁的過余溫度是最大的，并且遠(yuǎn)離樁壁呈指數(shù)衰減；隨著時間的增加，過余溫度持續(xù)傳遞到周圍的巖土體中，過余溫度逐漸增加，熱響應(yīng)半徑也隨之增大。

2 單工況荷載作用下熱響應(yīng)半徑

2.1 群樁埋管換熱器樁壁過余溫度計算

圖2 過余溫度沿徑向分布Fig.2 Excess temperature distribution in radial direction

在實際能源樁系統(tǒng)工程中，能源樁埋管換熱器通常由多個樁組成，布置方式主要有3 種類型：單排樁、雙排樁、多排樁(3排及3排以上)等間距布置。在多排布置方式下，樁編號為1，2，3，…，i，…，N。

地下巖土體假定為均勻、各向同性材料，熱傳導(dǎo)過程在均勻介質(zhì)中是線性的，根據(jù)疊加原理，編號為i的樁樁壁過余溫度可按下式計算[13]：

式中：θbi(τ)為第i個樁樁壁過余溫度，℃；θbi0(rb,τ)為第i個樁自身引起的樁壁過余溫度，℃；θij(rij,τ)為第j個樁引起的第i個樁樁壁附加過余溫度，℃；rb為樁半徑，m；τ為時間，s；rij為第j個樁至第i個樁的距離，m；li為第i個樁周圍其他樁對第i個樁的過余溫度影響系數(shù)。

2.2 熱響應(yīng)半徑定義

埋在地下巖土層中的能源樁埋管換熱器會對整個巖土層的溫度產(chǎn)生影響。對于單個能源樁，所產(chǎn)生的巖土體過余溫度隨離樁壁距離的增加而呈指數(shù)衰減，當(dāng)衰減至足夠小且在能源樁設(shè)計中可以忽略不計時，其對應(yīng)的距離便是能源樁埋管換熱器的熱響應(yīng)半徑。對于群樁埋管換熱器，還需要考慮其他能源樁所引起的溫度場的影響。在能源樁的整個系統(tǒng)運(yùn)行時間內(nèi)，當(dāng)周圍其他樁引起的過余溫度影響系數(shù)l≤5%時，工程上相鄰樁之間的熱干擾可以忽略不計，此時，相鄰能源樁中心線之間的垂直距離便可作為能源樁埋管換熱器的熱響應(yīng)半徑。一般在距離計算樁2倍樁距外的其他樁對計算樁的附加過余溫度影響系數(shù)很小，因此，對3 排以上多排樁可按3 排計算其熱響應(yīng)半徑[13]。

2.3 熱響應(yīng)半徑計算方法

參考文獻(xiàn)[8]，巖土熱擴(kuò)散系取0.5×10-6~1.6×10-6m2/s，假設(shè)能源樁系統(tǒng)處于冷熱負(fù)荷平衡狀態(tài)，計算時間可按單一制冷或者制熱工況取1~5月，樁基直徑按工程常規(guī)取0.4~1.2 m。確定能源樁的熱響應(yīng)半徑應(yīng)該考慮巖土體內(nèi)垂直于能源樁軸線方向的熱響應(yīng)最大區(qū)域，此時可忽略地表的影響。根據(jù)無限長線熱源模型(見式(1))，分別計算不同樁徑、不同巖土熱擴(kuò)散系數(shù)和不同運(yùn)行時間下能源樁埋管換熱器的熱響應(yīng)半徑，見圖3。經(jīng)擬合得：

式中：rc為熱響應(yīng)半徑，m。Ks和Kb分別為巖土介質(zhì)和樁半徑修正系數(shù)，當(dāng)樁為單樁時，Ks=1.68，Kb=2.13；當(dāng)樁單排布置時，Ks=2.10，Kb=1.98；當(dāng)樁雙排布置時，Ks=2.39，Kb=1.88；當(dāng)樁多排布置時，Ks=2.56，Kb=1.82。

圖3 不同條件下熱響應(yīng)半徑Fig.3 Thermal response radii at different conditions

2.4 計算結(jié)果與分析

表1所示為單樁，單排樁、雙排樁和多排樁布置狀況下，擬合公式(3)中rc計算結(jié)果與式(1)解析解的對比。從表1可以看出，擬合公式的計算結(jié)果與解析解之間最大相對誤差分別為4.62%，4.45%，3.77%和3.32%，結(jié)果表明擬合公式(3)具有較高的計算精度，符合工程要求。

3 熱響應(yīng)半徑影響因素分析

3.1 建筑負(fù)荷模型

本文利用浙江農(nóng)林大學(xué)衣錦校區(qū)學(xué)術(shù)交流中心作為建筑模擬對象，該建筑南北朝向，建筑面積為418 m2，運(yùn)用到能源樁系統(tǒng)的客房首層8 間，二層5 間，空調(diào)面積為162.79 m2，層高3.6 m。用DeST負(fù)荷軟件來模擬夏熱冬冷地區(qū)典型城市杭州1 a 的逐時冷熱負(fù)荷，建筑圍護(hù)結(jié)構(gòu)參數(shù)按照規(guī)范[19]所給出的不同熱工分區(qū)限值來設(shè)置。杭州采暖運(yùn)行期定于每年11月15日至次年3月15日，空調(diào)運(yùn)行期為每年05-15—09-15。圖4所示為杭州某建筑1 a逐時冷熱負(fù)荷的模擬結(jié)果。

在進(jìn)行動態(tài)負(fù)荷作用下能源樁系統(tǒng)設(shè)計時，根據(jù)建筑物冷熱負(fù)荷峰值與所在地建議的能源樁每延米取放熱量計算最大樁數(shù)；由DeST負(fù)荷軟件模擬得出的建筑全年動態(tài)逐時冷熱負(fù)荷和熱泵機(jī)組的Ccop(制冷性能系數(shù))和Hcop(制熱性能系數(shù))可計算出地下巖土體需要承受的逐時冷熱負(fù)荷；將上述二者結(jié)合可以獲得每個能源樁需要承擔(dān)的逐時冷熱負(fù)荷。

3.2 數(shù)值模型

由于現(xiàn)有解析模型假設(shè)熱源每延米換熱功率是恒定的，而實際工程項目中能源樁每延米換熱功率是關(guān)于時間的函數(shù)，為了能準(zhǔn)確研究動態(tài)負(fù)荷作用下的熱響應(yīng)半徑，本文利用COMSOL Multiphysics 有限元軟件建立數(shù)值模型進(jìn)行計算分析。

3.2.1 幾何模型建立與材料屬性賦予

單樁內(nèi)建立1根U型管。U型管在樁中，樁和模擬范圍內(nèi)的巖土體可看成是圓柱體。模型的幾何參數(shù)參考規(guī)范[8]給出的范圍，U 型管內(nèi)徑取25 mm，外徑為32 mm，埋深為30 m，支管間距為180 mm，樁徑為600 mm；巖土體模擬半徑根據(jù)文獻(xiàn)[13]中的模型的計算時間和巖土體的熱擴(kuò)散系數(shù)取6 m。U 型管為聚乙烯管，管內(nèi)的換熱液為水，熱物性參數(shù)如表2所示。

3.2.2 邊界條件和初始條件

一般來說，地表以下5 m的巖土體溫度全年基本不變，雖有波動，但這段距離占樁長的比例很小，因此，將巖土體最遠(yuǎn)邊界處和底部設(shè)定為巖土體遠(yuǎn)邊界恒溫條件；巖土體表面暫不考慮外界四季溫度變化所帶來的影響，也設(shè)置為恒溫邊界條件；根據(jù)埋管內(nèi)徑，U 型管入口流速設(shè)定為0.6 m/s。假定能源樁以每延米50 W/m 的恒定功率放熱；巖土體和管內(nèi)水的初始溫度為17.3 ℃。

表1 熱響應(yīng)半徑rc計算結(jié)果與式(1)解析解對比Table1 Comparison ofrcobtained from of the fitting formula and analytical solution of equation(1)

圖4 杭州某建筑1 a逐時冷熱負(fù)荷模擬結(jié)果Fig.4 Simulation results of hourly cooling and heating load of a building in Hangzhou during one year

表2 熱物性參數(shù)Table2 Thermophysical parameters

3.2.3 網(wǎng)格劃分與求解計算

網(wǎng)格劃分時對U 型管周圍進(jìn)行了詳細(xì)劃分，特別是在入口和出口處，巖土體由里到外網(wǎng)格密度由密到疏，網(wǎng)格劃分如圖5所示。劃分網(wǎng)格時，選取3 種不同網(wǎng)格單元數(shù)(138 855，252 751 和572 249 個)計算的過余溫度相對誤差在5%之內(nèi)，滿足精度要求，因此，選擇第1種網(wǎng)格單元數(shù)進(jìn)行劃分以節(jié)約計算成本。使用時間步長為1 h的瞬態(tài)求解器求解模型。

3.2.4 模型驗證

為了驗證所建三維數(shù)值模型的正確性，對比無限長線熱源傳熱解析模型，每延米放熱功率統(tǒng)一取50 W/m，計算時間為3月，巖土體的熱物性參數(shù)與表2中的保持一致。過余溫度數(shù)值模型與解析模型計算結(jié)果對比如圖6所示。由圖6可知：過余溫度數(shù)值模型的計算結(jié)果和解析模型趨勢基本一致，驗證了本文三維數(shù)值模型建立的正確性。

圖5 數(shù)值模型網(wǎng)格與幾何示意圖Fig.5 Mesh and geometry diagrams of numerical model

圖6 數(shù)值模型與解析模型中過余溫度計算結(jié)果對比Fig.6 Comparison of calculation results of excess temperature between the numerical model and the analytical model

3.3 影響因素分析

能源樁的運(yùn)行工況往往是周期性的，在1個周期內(nèi)存在制熱運(yùn)行期、制冷運(yùn)行期以及間歇期，運(yùn)行期內(nèi)負(fù)荷隨時間變化且存在1個峰值負(fù)荷，間歇期沒有負(fù)荷作用。為此，本文分別研究不同荷載分布、不同峰值負(fù)荷、不同運(yùn)行期以及全生命周期對熱響應(yīng)半徑的影響。

3.3.1 荷載分布的影響

負(fù)荷采用3.1 節(jié)中的計算結(jié)果。采用恒溫法[20]計算得到杭州地區(qū)的能源樁每延米換熱功率約為50 W/m；熱泵機(jī)組的額定工況Hcop取4.5；制熱工況運(yùn)行時間為4月，從每年11月15日到次年的3月15日。隨時間變化的負(fù)荷也可簡化成矩形負(fù)荷，負(fù)荷峰一致，見圖7。

圖7 不同荷載分布形式Fig.7 Different load distributions

表3所示為運(yùn)行時間為4月時熱響應(yīng)半徑rc計算結(jié)果(見式(3))與2 種荷載分布下rc數(shù)值解對比。從表3可以看出：由本文擬合公式所得rc計算結(jié)果與2 種荷載下rc數(shù)值解的最大相對誤差分別為3.07%和2.01%，表明相同作用時間內(nèi)不同荷載分布的變化對熱響應(yīng)半徑的影響較小，在實際工程中可將逐時負(fù)荷簡化成矩形負(fù)荷來計算。

3.3.2 峰值負(fù)荷的影響

取峰值負(fù)荷分別為900，1 200 和1 500 W，當(dāng)運(yùn)行時間為4月時本文rc計算結(jié)果(見式(3))與不同峰值負(fù)荷下數(shù)值解對比如表4所示。從表4可以看出：本文擬合公式rc計算結(jié)果與3種峰值負(fù)荷數(shù)值解的最大相對誤差均為2.01%，表明不同峰值負(fù)荷對能源樁熱響應(yīng)半徑影響很小。

表3 熱響應(yīng)半徑計算結(jié)果與不同荷載分布下數(shù)值解對比Table3 Comparison of calculation results and numerical solutions of thermal response radius with different load distributions

3.3.3 運(yùn)行期的影響

取運(yùn)行期分別為3月、4月和5月，當(dāng)峰值負(fù)荷為1 200 W 時本文rc計算結(jié)果(見式(3))與不同運(yùn)行期數(shù)值解對比如表5所示。從表5可以看出：擬合公式中rc計算結(jié)果與3種不同運(yùn)行期數(shù)值解的最大相對誤差分別為2.67%，2.01%和2.02%，表明運(yùn)行期時間越長，能源樁熱響應(yīng)半徑越大。

3.3.4 全生命周期的影響

在能源樁埋管換熱器實際運(yùn)行中，通常是由加熱和冷卻2種過程組成的。假定荷載分布為矩形負(fù)荷，峰值負(fù)荷為1 500 W，制冷制熱運(yùn)行期為4月，間歇期為2月，運(yùn)行30 a。

能源樁全生命周期內(nèi)冷熱負(fù)荷的總量可以相互抵消，計算時可取冷熱負(fù)荷運(yùn)行時最大熱響應(yīng)半徑。在能源樁30 a 運(yùn)行周期內(nèi)，樁壁處最大溫升為11.62 ℃，最大熱響應(yīng)半徑為4.76 m。在相同工況下，式(3)熱響應(yīng)半徑計算結(jié)果為4.59 m，兩者相對誤差為3.57 %，表明在計算整個生命周期內(nèi)的熱響應(yīng)半徑時，擬合計算公式(3)仍具有較高的計算精度。

4 能源樁全生命周期熱響應(yīng)半徑簡化計算方法

為便于工程應(yīng)用，取實際工程中應(yīng)用最多的多排樁布置方式，計算時間為5月。分別計算不同樁徑和不同巖土體熱擴(kuò)散系數(shù)條件下能源樁熱響應(yīng)半徑隨時間的變化關(guān)系，見圖8。圖中虛線表示的是“地源熱泵系統(tǒng)工程技術(shù)規(guī)范”[8]中所推薦的鉆孔間距范圍。從圖8可以看出：能源樁埋管換熱器的熱響應(yīng)半徑隨著樁周巖土體熱擴(kuò)散系數(shù)的增大而增大，隨系統(tǒng)運(yùn)行時間的增加而增大，隨能源樁樁徑的增大而增大。從圖8還可以看出：若能源樁在熱擴(kuò)散系數(shù)為0.5×10-6m2/s 的巖土體中運(yùn)行時，按單一工況運(yùn)行4月，能源樁熱響應(yīng)半徑便超過了地源熱泵規(guī)范所推薦的上限值。因此，在設(shè)計能源樁間距時，不能簡單采取地源熱泵規(guī)范所給的值，還需根據(jù)實際現(xiàn)場地質(zhì)條件、運(yùn)行時間以及樁徑來選取。

表5 不同運(yùn)行期熱響應(yīng)半徑計算結(jié)果與數(shù)值解對比Table5 Comparison of calculation results and numerical solutions of thermal response radius at different operating periods

圖8 能源樁不同樁徑、不同熱擴(kuò)散系數(shù)下的熱響應(yīng)半徑Fig.8 Thermal response radius of energy piles with different pile diameters and thermal diffusivities

5 結(jié)論

1)巖土體中的過余溫度在樁壁處最大，隨著與樁壁距離的增加而迅速衰減，隨時間的增加而增大。

2)在單樁、單排樁、雙排樁和多排樁布置情況下，本文擬合公式的熱響應(yīng)半徑計算結(jié)果與無限長線熱源模型解析解之間的最大相對誤差分別為4.62%，4.45%，3.77%和3.32%，表明該方法具有較高的計算精度，符合工程要求。

3)荷載分布形式以及峰值負(fù)荷只影響巖土體內(nèi)過余溫度及其梯度；巖土體熱擴(kuò)散系數(shù)及運(yùn)行時間影響溫度傳遞區(qū)域；全生命周期內(nèi)冷熱負(fù)荷的總量可以相互抵消，最大熱響應(yīng)半徑計算可按單工況荷載作用下熱響應(yīng)半徑計算方法。

4)選取工程中應(yīng)用最多的多排樁布置形式，繪制了不同樁徑和不同巖土體熱擴(kuò)散系數(shù)條件下能源樁熱響應(yīng)半徑隨時間的變化關(guān)系，為工程應(yīng)用提供參考。