李澤錕， 杜震宇

(太原理工大學 土木工程學院， 山西 太原 030024)

隨著建筑能耗的增加和環(huán)境問題的日益嚴峻，開發(fā)和利用可再生能源受到越來越多的關注[1].土壤源熱泵因充分利用淺層地熱、穩(wěn)定性高、地上空間占用較小等優(yōu)點受到越來越多的關注.地埋管換熱器作為土壤源熱泵的核心組成部分，其換熱能力受周圍巖土體的熱物性和地下水的影響較大.為了分析地埋管換熱器的換熱量和周圍巖土體的溫度變化，國內外學者構建了許多地埋管換熱器的傳熱模型.黃雪婷[2]通過分層滲流模型，研究地下水位線高度和地下水流速對地埋管換熱器換熱的影響.Abdelaziz等[3]根據(jù)有限長線熱源模型，得到巖土體內任意一點的溫度響應，并根據(jù)各層的導熱系數(shù)得到分層換熱量的經(jīng)驗公式.文獻[4-5]通過有限長線熱源模型和移動線熱源模型，建立巖土體分層滲流解析模型.Sel?uk等[6]對文獻[3]的模型進行優(yōu)化，并考慮滲流的作用.李永等[7]通過熱物性測試實驗的數(shù)據(jù)，建立三維分層巖土體模型，并考慮了巖土體的初始溫度.張山等[8]建立三維非穩(wěn)態(tài)傳熱模型，分析地下水滲流、流速和含水層厚度對地埋管換熱器周圍巖土體溫度場和換熱量的影響.文獻[9-10]建立三維數(shù)值模型，將含水層處的導熱系數(shù)等效為無滲流的導熱系數(shù)，并研究含水層厚度對熱物性檢測的影響.

目前，關于地下水對地埋管換熱器換熱的研究還相對較少.因此，本文建立地埋管換熱器的三維非穩(wěn)態(tài)分層滲流傳熱模型，分析含水層厚度、地下水流速和巖土體孔隙率對地埋管換熱器換熱的影響.

1 傳熱模型的建立

1.1 模型的假設

地埋管換熱器的換熱過程是復雜的非穩(wěn)態(tài)換熱，需對模型進行簡化，有以下5個假設.1) 對計算區(qū)域的巖土體進行分層，各層巖土體視為均質介質，且?guī)r土體熱物性不隨時間和空間而變化.2) 假設循環(huán)流體、地埋管換熱器、回填材料和巖土體間不存在接觸熱阻.3) 假設地下水滲流僅存在于水平方向，且速度恒定.4) 假設固體和液體在接觸的瞬間達到熱平衡狀態(tài).5) 忽略地埋管換熱器底部彎管的影響.

1.2 幾何模型

將某項目土壤源熱泵系統(tǒng)作為建立豎直地埋管換熱器傳熱模型的基礎[11]，該項目處于典型的黃土高原地區(qū)，根據(jù)該地區(qū)的地質條件，沿軸向將巖土體分為5層，自上而下分別為濕陷性黃土、黃土、卵石、強風化砂巖和中風化砂巖.該項目采用單U型地埋管換熱器，管內換熱介質為水，回填材料為砂土，U型管為聚乙烯(PE)管，傳熱模型的熱物性參數(shù)，如表1所示.表1中：ρ為密度；α為熱擴散系數(shù)；λ為導熱系數(shù)；φ為孔隙率；K為滲透率.

表1 傳熱模型的熱物性參數(shù)Tab.1 Thermal physical property parameters of heat transfer model

圖1 地埋管換熱器及巖土體模型圖(單位：m)Fig.1 Model of ground heat exchanger and rock-soil body (unit： m)

為了防止邊界對模擬的干擾，將計算區(qū)域設置為6 m×6 m×80 m(長×寬×高)的六面體巖土體，鉆井(圓柱體)設置于計算區(qū)域的幾何中心，其埋深為75 m.地埋管換熱器及巖土體模型圖，如圖1所示.

鉆井與鉆井內部的相關參數(shù)，如表2所示.表2中：d為鉆井直徑；b0為支管間距；d1為U型管外徑；d2為U型管內徑；L為地埋管換熱器的埋深.

通過Fluent流體仿真軟件，建立地埋管換熱器的三維非穩(wěn)態(tài)分層滲流傳熱模型.對幾何模型進行網(wǎng)格劃分，由于計算區(qū)域的溫度軸向變化相對較小，故豎直方向的網(wǎng)格相對稀疏，而水平方向的網(wǎng)格則需進行加密處理.忽略底部彎管的影響，將U型管簡化為1根進水管和1根出水管，通過用戶自定義函數(shù)(UDF)，將進水支管與回水支管的底部相連接.

表2 鉆井與鉆井內部的相關參數(shù)Tab.2 Related parameters of borehole and borehole internal (mm)

1.3 數(shù)學模型

1.3.1 多孔介質的控制方程 多孔介質的質量守恒方程為

(1)

式(1)中：ρs，w為多孔介質中流體密度；u為多孔介質中流體流速；Q為控制體內熱源的強度；t為時間.

多孔介質的動量守恒方程為

(2)

多孔介質的能量守恒方程為

(3)

式(3)中：cp，s為多空介質中固體的體積比熱容；θs為巖土體溫度；(ρcp，s)tot為多孔介質總的體積比熱容，(ρcp，s)tot=φ(ρcp，s)s，w+(1-φ)(ρcp，s)s；λtot為多孔介質總的導熱系數(shù)，λtot=φλs，w+(1-φ)λs，λs，w，λs分別為多孔介質中流體和固體的導熱系數(shù).

1.3.2 管內循環(huán)流體的控制方程 管內循環(huán)流體的質量守恒方程為

(4)

式(4)中：ρw為管內流體的密度，當流體為不可壓縮流體時，密度為定值；U為管內流體的流速.

管內循環(huán)流體的動量守恒方程為

(5)

管內循環(huán)流體的能量守恒方程為

(6)

式(6)中：θw，λw，cp，w分別為管內流體的溫度、導熱系數(shù)和體積比熱容；Sθw為該方程的源項.

1.3.3 數(shù)值計算方法 采用Realizablek-ε雙方程湍流模型，通過壓力隱式算子分割(PISO)算法對狀態(tài)方程進行求解.離散格式為二階迎風格式.進行非穩(wěn)態(tài)計算時，時間步長設置為60 s，收斂條件的能量殘差設置為1×10-8，其余殘差設置為1×10-6.

1.4 定解條件

1.4.1 初始條件 將巖土體設置為多孔介質，巖土體周圍的邊界設置為恒溫邊界，其溫度分布與巖土體初始溫度分布相同.進行地埋管換熱器的模擬計算時，巖土體初始溫度對結果影響較大[12]，故需考慮巖土體軸向的初始溫度變化.經(jīng)實際測量，巖土體并未出現(xiàn)明顯的恒溫層.深度為0～20 m的巖土體溫度受外界溫度的影響較大，深度為20 m以下的巖土體溫度因地熱作用隨深度的增加而升高.考慮地熱作用，對實測數(shù)據(jù)進行擬合，可得巖土體初始溫度軸向變化的公式為

θs，0=-0.012z+11.04， 0≤z≤20，

(7)

θs，0=0.016z+10.46， 20

(8)

式(7)，(8)中：θs，0為巖土體初始溫度，℃；z為巖土體的深度，m.

1.4.2 邊界條件 當埋深為0.005，0.045 m時，地埋管換熱器換熱受太陽輻射的影響較大；當埋深為0.300 m以下時，地埋管換熱器換熱幾乎未受太陽輻射的影響[13].由于文中地埋管換熱器的埋深為70 m，地埋管換熱器的頂端距離地面為5 m，因此，可忽略太陽輻射的影響，僅考慮外界空氣溫度變化對地埋管換熱器換熱的影響.將計算區(qū)域頂部設置為第3類邊界條件，其中，流體溫度設定為夏季典型日隨時間變化的空氣溫度，則外界空氣溫度θup為

(9)

利用季節(jié)平均溫度，經(jīng)計算，可得地面與周圍空氣的對流換熱系數(shù)為0.48 W·(m·K)-1.在含水層設置一定流速的地下水滲流，地下水溫度與巖土體初始溫度相同.U型管入口設置為速度入口，出口設置為壓力出口.

1.5 網(wǎng)格獨立性檢驗

作為數(shù)值解的網(wǎng)格應足夠細密，使進一步加密網(wǎng)格對數(shù)值解幾乎沒有影響.這種數(shù)值解稱為網(wǎng)格獨立的解[14].對網(wǎng)格數(shù)量分別為223 327，497 680，791 652，1 062 233，1 267 854的模型進行模擬，設置進口溫度為14 ℃，進口流速為0.23 m·s-1.

待巖土體的溫度場穩(wěn)定后，進行網(wǎng)格獨立性檢驗，其結果如圖2所示.圖2中：θs，h為回水側巖土體軸向溫度；θout為出水溫度.由圖2可知：當網(wǎng)格數(shù)量不低于497 680時，回水側巖土體軸向溫度分布基本不變.綜合考慮計算效率和模型的精確度，選擇497 680個網(wǎng)格的模型.

(a) 回水側巖土體軸向溫度 (b) 出水溫度 圖2 網(wǎng)格獨立性檢驗Fig.2 Grid independence verification

2 模擬與結果分析

2.1 評價指標

單位井深換熱量是衡量地埋管換熱器換熱能力的參數(shù)，表達式為

(10)

根據(jù)文獻[15]，分層單位井深換熱量可直觀地反映各層巖土體的換熱能力，表達式為

(11)

式(11)中：qi為第i層巖土體的單位井深換熱量，W·m-1；θin，i，θout，i分別為第i層巖土體頂面的進水溫度和出水溫度，℃；Li為第i層巖土體的厚度，m.

為了更加簡單、有效地觀察各層巖土體的換熱量占總換熱量的份額，確定最優(yōu)的埋深，提出分層單位井深換熱比率，即

(12)

式(12)中：Ni為第i層的分層單位井深換熱比率.

地下水對第i層巖土體位置處地埋管換熱器換熱的影響可通過該層單位井深換熱量增加率進行表示，即

(13)

2.2 模型的驗證

為了驗證模型的準確性，選取文獻[11]中2011年7月24日在不同埋深的供、回水側溫度測量數(shù)據(jù)(實驗值)與模擬結果(模擬值)進行對比，入口溫度隨該日逐時冷負荷的變化而變化.

夏季制冷工況下，供、回水側溫度模擬值與實驗值的對比，如表3所示.表3中：θexp為供、回水側溫度的實驗值；θc為供、回水側溫度的模擬值；η為供、回水側溫度的實驗值和模擬值的相對誤差.

由表3可知：模擬值與實驗值吻合較好，兩者之間的相對誤差不超過3%，誤差產(chǎn)生的原因可能是實驗中的埋管以管群形式存在，各埋管相互之間有一定的干擾；模擬值與實驗值的誤差相對較小，且在允許范圍內，故該模型具有準確性.

表3 供、回水側溫度模擬值與實驗值的對比Tab.3 Comparison between simulation values and experimental values of water supply and return side temperatures

2.3 地下水對地埋管換熱器換熱的影響

無量綱貝克萊數(shù)是滲流產(chǎn)生的對流和熱傳導的影響程度比值，可用于判斷地埋管換熱器換熱受地下水影響的大小，其表達式為

Pe=ρucuvlm/λtot.

(14)

式(14)中：Pe為貝克萊數(shù)；ρu為地下水的密度，kg·m-3；cu為地下水的比熱容，J·(kg·K)-1；v為地下水的流速，m·s-1；lm為特征長度，文中為含水層厚度，m.

由貝克萊數(shù)的表達式可知，其大小受lm，v，λtot的影響，而λtot與λs，φ有關.為了研究地下水對地埋

圖3 Δqi隨Pe的變化情況(lm=20 m)Fig.3 Variation of Pe with Δqi (lm=20 m)

管換熱器的換熱作用，分別模擬不同的含水層厚度、地下水流速和巖土體孔隙率對Δqi的影響.

假設含水層的頂部位于深度27.6 m處，當含水層厚度為20 m時，模擬地下水流速分別為3.0×10-5，2.5×10-5，1.5×10-5，1.0×10-5，7.0×10-6，5.0×10-6，3.0×10-6，1.0×10-6m·s-1的情況，可得到該層單位井深換熱量增加率.同時，在上述地下水的流速下，模擬巖土體孔隙率分別為20%，30%，40%，50%，60%，70%，80%的情況，可得到該層單位井深換熱量增加率.Pe的范圍為46.842～2 810.520，Δqi隨Pe的變化情況，如圖3所示.

由圖3可知：地下水對地埋管換熱器換熱具有促進作用，但當?shù)叵滤魉佥^低或巖土體孔隙率較大時，地下水反而會對換熱產(chǎn)生抑制作用，即Pe為46.842～93.684時，單位井深換熱量增加率Δqi小于0；在含水層厚度為20 m的情況下，當Pe大于140.52時，地下水對地埋管換熱器換熱產(chǎn)生顯著影響，當Pe小于140.52時，單位井深換熱量增加率Δqi小于0.001，地下水對換熱的影響遠小于導熱，地下水的作用可以忽略.

將Pe與Δqi進行擬合，得到當含水層厚度為20 m時，Δqi隨Pe的變化方程為

(15)

式(15)中：R2為擬合度.

不同含水層厚度的情況下，Δqi隨Pe的變化情況，如圖4所示.圖4中：含水層厚度為5，10，15，25，30 m時，對應的Pe取值范圍為[11.71～702.63]， [23.42～1 405.26]， [35.13～2 107.89]，[58.55～3 513.15]，[70.26～4 215.78].

不同含水層厚度對應的Pe臨界值，如表4所示.表4中：Pecr為Pe的臨界值.

圖4 Δqi隨Pe的變化情況 Fig.4 Variation of Pe with Δqi

由圖4和表4可得以下3個結論.

1)Pe的臨界值隨著含水層厚度的增加而增大；當Pe大于臨界值時，地下水對地埋管換熱器換熱的影響比較顯著；當Pe小于臨界值時，地下水對地埋管換熱器換熱的影響可以忽略.

2) 隨著含水層厚度的增加，含水層單位井深換熱量增加率逐漸減小；當含水層厚度小于15 m時，單位井深換熱量增加率最大時接近0.5；當含水層厚度大于15 m時，單位井深換熱量增加率最大值逐漸減??；當含水層厚度為30 m時，單位井深換熱量增加率僅為0.32；當含水層厚度一定時，單位井深換熱量增加率的曲線斜率隨著Pe的增大而減小.

3) 并不是含水層厚度越大，其換熱效率就越高，地下水的作用具有一定限度.在Pe較小的情況下，地下水作用增大，導致?lián)Q熱量變化的幅度更加明顯.

當含水層厚度分別為5，10，15，25，30 m時，單位井深換熱量增加率Δqi和Pe成對數(shù)關系.Δqi和Pe之間變化關系的擬合曲線方程分別為

表4 不同含水層厚度對應的Pe臨界值Tab.4 Critical values of Pe with different aquifer thicknesses

Δqi=0.222 5ln(Pe)-0.966 7，R2=0.996 5，

(16)

Δqi=0.194 0ln(Pe)-0.859 5，R2=0.997 3，

(17)

Δqi=0.162 2ln(Pe)-0.754 1，R2=0.995 2，

(18)

Δqi=0.132 6ln(Pe)-0.704 0，R2=0.993 8，

(19)

Δqi=0.110 9ln(Pe)-0.587 2，R2=0.997 7.

(20)

在含水層厚度不同的情況下，地埋管換熱器在地下水作用下的單位井深換熱量增加率隨Pe變化的曲線方程為Δqi=aln(Pe)+b.其中，系數(shù)a，b隨lm的變化而變化.

通過擬合，可以得到a，b關于含水層厚度lm的曲線方程，分別為

(21)

(22)

3 新模型的建立與驗證

3.1 修正導熱系數(shù)

由于每層巖土體的單位井深換熱量不同，可以通過各層巖土體的厚度和導熱系數(shù)的比值，得到各層巖土體的分層單位井深換熱量[2]，即

(23)

式(23)中：λi為第i層巖土體的導熱系數(shù)，W·(m·K)-1.

式(23)將地埋管換熱器換熱的熱流密度分為平均部分和加權部分，平均部分反映鉆井內部均勻的換熱量，加權部分則反映不同外部巖土體產(chǎn)生的差異.該公式在不考慮地下水作用時比較準確，但當?shù)叵滤淖饔幂^為明顯時，加權部分會產(chǎn)生較大的誤差.如果將滲流作用考慮在內，則需對含水層處的導熱系數(shù)進行修正.假設計算區(qū)域內的巖土體軸向分為n層，其中，第m層為含水層，通過Δqi和式(23)，可得修正導熱系數(shù)公式為

.

(24)

3.2 新模型的驗證

通過COMSOL仿真軟件，建立地埋管換熱器的軸向-周向二維模型，并將鉆井等效為長度為70 m的線熱源，其頂部距離地表5 m.根據(jù)文獻[11]的方法，對計算區(qū)域進行豎向分層.其中，含水層深度為27.6～37.6 m，地下水流速為1×10-5m·s-1，流動僅存在于水平方向.通過式(14)進行計算，可得Pe為409.3.通過式(12)，(13)，(24)進行計算，可得修正后的各層導熱系數(shù)及巖土體的分層單位井深換熱比率，結果如表5所示.

將各層巖土體的單位井深換熱量代入對應的線熱源模型，可建立新模型.

表5 各層巖土體的相關參數(shù)Tab.5 Related parameters of rock-soil in each layer

對新模型分別進行短期模擬和長期模擬，通過地埋管換熱器周圍巖土體的溫度分布對模型進行驗證.短期模擬是將新模型與Fluent模型分別模擬4 d，對比模擬結束后距鉆井壁0.5 m處巖土體溫度的軸向變化，以及地下水下游的鉆井壁溫度在埋深35 m處隨時間的變化，以此驗證新模型的準確性.

新模型與Fluent模型的模擬對比，如圖5所示.圖5中：θb為鉆井壁溫度.

(a) 巖土體溫度 (b) 鉆井壁溫度 圖5 新模型與Fluent模型的模擬對比Fig.5 Simulation comparison between new model and Fluent model

由圖5可得以下4個結論.

1) 兩個模型的誤差在5%以內，變化趨勢幾乎完全相同，新模型的溫度略高于Fluent模型，這可能是因為新模型忽略了回填材料對換熱的影響.

2) 在含水層區(qū)域，由于受到地下水的影響，巖土體的熱擴散能力有較大的提升，地埋管換熱器供水側巖土體溫度幾乎與初始溫度相同，溫差不超過0.5 ℃.

3) 在地下水流動的下游，在4 d的運行后，兩個模型在埋深35 m處鉆井壁溫度隨時間變化的曲線趨勢基本相同，誤差很小，最大誤差僅為4.2%.

4) 新模型對短時間內巖土體溫度場分布的模擬較為準確，而且新模型的計算時間遠小于Fluent模型，對于整個供冷期的模擬時間僅為286 s左右.在研究地埋管換熱對巖土體溫度分布的影響時，新模型具有較大優(yōu)勢.

圖6 巖土體溫度的軸向變化Fig.6 Axial variation of rock-soil temperature

長期模擬是對整個供冷期進行模擬，可得地下水下游距鉆井壁2.5 m處巖土體溫度的軸向變化，如圖6所示.

由圖6可得以下2個結論.

1) 新模型與實驗數(shù)據(jù)的誤差在10%以內，誤差在工程允許范圍內，故新模型對長期溫度場分布的模擬具有較高的準確性.

2) 相較于實驗數(shù)據(jù)，新模型產(chǎn)生的誤差主要集中于地埋管換熱器的后半段，這是因為新模型中分層換熱量的大小與各層巖土體的厚度具有相關性，在導熱系數(shù)和各巖土層厚度確定時，各層單位井深換熱比率不變，故分層單位井深換熱量僅隨輸入負荷的變化而變化.

實驗數(shù)據(jù)取自黃土高原地區(qū)某處土壤源熱泵系統(tǒng)，在此系統(tǒng)中，雖然50～70 m處巖土體的導熱系數(shù)較高，但由于熱回流的存在，導致地埋管換熱器在57～70 m處的換熱效率很低，因此，新模型在后半段產(chǎn)生較大的誤差，新模型對埋深合理的地埋管換熱器的準確性更高.

4 結論

基于黃土高原地區(qū)某土壤源熱泵的實驗數(shù)據(jù)，通過Fluent流體仿真軟件建立地埋管換熱器三維非穩(wěn)態(tài)分層換熱模型，并借助實驗數(shù)據(jù)對模型的準確性進行驗證.

在不同的含水層厚度、地下水流速和巖土體孔隙率的條件下，利用該模型分析貝克萊數(shù)與地埋管換熱器換熱之間的關系.通過含水層處換熱量增加率，對導熱系數(shù)進行修正，進而可得到計算速度更快的換熱模型.由此得到以下3個結論.

1) 當含水層厚度為20 m時，對不同地下水流速、巖土體孔隙率的情況下，地下水對地埋管換熱器的單位井深換熱量增加率進行模擬.由模擬結果可知，在以含水層厚度為特征長度的情況下，貝克萊數(shù)存在臨界值，當貝克萊數(shù)大于臨界值時，地下水對地埋管換熱器換熱具有促進作用.在流速較小或巖土體孔隙率較大的情況下，地下水反而對地埋管換熱器換熱起到抑制的作用.

2) 對含水層厚度分別為5，10，15，25，30 m的情況進行模擬，得到不同含水層厚度下貝克萊數(shù)的臨界值.同時，通過模擬數(shù)據(jù)可知，單位井深換熱量增加率與貝克萊數(shù)成對數(shù)關系.建立不同含水層厚度的貝克萊數(shù)與單位井深換熱量增加率Δqi之間的方程，方程的形式為Δqi=aln(Pe)+b，并通過數(shù)據(jù)擬合得到a，b關于含水層厚度的方程.

3) 通過單位井深換熱量增加率對含水層處的導熱系數(shù)進行修正，利用修正后的導熱系數(shù)，可得到分層單位井深換熱量，進而建立新的地埋管換熱器分層滲流傳熱模型.將新模型與Fluent模型、實驗數(shù)據(jù)進行對比，可證明無論是短期模擬，還是長期模擬，新模型用于地埋管換熱器周圍巖土體的溫度場分布研究時具有較高的準確性，且計算耗時較短.