三維與二維結(jié)合地源熱泵地下?lián)Q熱器復(fù)合數(shù)值模型

T.M.Eikevik

(1 上海交通大學(xué)機(jī)械與動(dòng)力工程學(xué)院 上海 200240； 2 挪威科技大學(xué) 特隆赫姆 7491)

鑒于我國(guó)面臨的嚴(yán)峻能源和碳排放形勢(shì)[1]，地源熱泵技術(shù)以低能耗特點(diǎn)如今成為傳統(tǒng)空調(diào)系統(tǒng)的替代方案之一。地埋管換熱器作為地源熱泵系統(tǒng)最主要的部件之一備受研究者關(guān)注并由此演化出諸多模型。L. R.Ingersoll等[2]提出的無(wú)限線熱源模型屬于最早期的地埋管模型。該模型計(jì)算簡(jiǎn)便但只考慮了土壤中徑向的熱交換而忽略了軸向的熱交換。而這一點(diǎn)被D. Marcotte 等[3-5]證明了重要性。此外，無(wú)限線源模型鑒于其積分式特點(diǎn)只適用于相對(duì)管徑較小的場(chǎng)合[6]。另一種廣泛使用的解析模型是L.Ingersioll等[7-9]提出的無(wú)限圓柱熱源模型。該模型由無(wú)限線熱源模型發(fā)展而來(lái)，并繼承了無(wú)限線熱源模型純熱傳導(dǎo)假設(shè)和不考慮軸向熱交換的特點(diǎn)，因而同樣不適用于長(zhǎng)期模擬[10]。為了彌補(bǔ)這一缺陷，P. Eskilson[11]引入了有限線熱源模型。P. Eskilson的模型考慮了軸向換熱及多管疊加的熱效應(yīng)，但其缺陷在于計(jì)算相對(duì)耗時(shí)，且和其他解析模型一樣設(shè)定沿管長(zhǎng)的熱流密度為一定值。

相較于解析模型，數(shù)值模型通常有高精度和對(duì)應(yīng)實(shí)際工程場(chǎng)景可定制的特點(diǎn)。然而，數(shù)值模型尤其是三維數(shù)值模型對(duì)計(jì)算量的要求極高[12-13]。因此多篇論文采用了不同方法來(lái)降低計(jì)算耗時(shí)，這一嘗試主要分成普通三維模型和準(zhǔn)三維模型兩類(lèi)：一類(lèi)是以T.Y.Ozudogru 等[14]為代表的普通三維模型，該模型通過(guò)一維線性單元來(lái)簡(jiǎn)化模擬過(guò)程卻刻意忽略了管內(nèi)液體的流動(dòng)。不考慮流動(dòng)的影響在流速較大、管長(zhǎng)較長(zhǎng)的情況下尤為明顯。另一類(lèi)是以Li Zhongjian[15]為代表的準(zhǔn)三維模型，該模型考慮到模型的對(duì)稱(chēng)性因而只計(jì)算了其中一邊。這一類(lèi)準(zhǔn)三維模型[16-18]的主要特點(diǎn)是不考慮地表和埋管底部的影響直接把模型簡(jiǎn)化為水平二維模型，同樣不能反映軸向熱交換，因而不適用于長(zhǎng)期模擬。

綜上所述，這些為了追求減少計(jì)算耗時(shí)的簡(jiǎn)化削弱了原三維模型的準(zhǔn)確性。很多簡(jiǎn)化只是局部性的，不能從根本上改變?nèi)S模型耗時(shí)的特點(diǎn)，無(wú)法滿足不規(guī)則管群溫度分布的模擬。由此可見(jiàn)，單獨(dú)一個(gè)的三維模型或二維模型不能很好地模擬真實(shí)氣候條件下的實(shí)際問(wèn)題。

本文提出了一個(gè)三維與二維結(jié)合地源熱泵地下?lián)Q熱器復(fù)合數(shù)值模型，能引入實(shí)際氣候數(shù)據(jù)的影響，獲得不規(guī)則管群周邊溫度場(chǎng)分布情況，且提高計(jì)算效率。在新建的復(fù)合模型中，三維模型導(dǎo)出的模擬結(jié)果被用作二維模型輸入的邊界條件。結(jié)合兩者的復(fù)合模型既發(fā)揮了三維模型準(zhǔn)確性高的優(yōu)點(diǎn)，又通過(guò)二維模型顯著降低了管群模擬的運(yùn)算量。

1 換熱器模型的建立

復(fù)合模型由一個(gè)三維模型和一個(gè)二維模型組成。三維模型由兩組U形管、管內(nèi)流體、回填材料、管外及表層土壤組成。U形管外土壤的計(jì)算區(qū)域半徑為2 m，計(jì)算區(qū)域深度距U形管底部2 m，具體結(jié)構(gòu)如圖1所示，換熱模型參數(shù)如表1所示。二維模型可根據(jù)實(shí)際管群分布建立。實(shí)際管群分布為呈直角分布的兩組矩形埋管陣列，管外土壤的計(jì)算半徑為4 m，如圖2所示。選擇陣列中最具代表性的共27根地埋管組成的直角區(qū)域進(jìn)行建模。為進(jìn)一步簡(jiǎn)化，二維模型采用單位管化，相應(yīng)的參數(shù)也做了適當(dāng)調(diào)整[19]，相應(yīng)的網(wǎng)格如圖3所示，其中離U形管較近的網(wǎng)格較為緊湊，離U形管較遠(yuǎn)的網(wǎng)格較為稀疏。復(fù)合模型中，三維模型在運(yùn)行狀態(tài)下不同深度的鉆孔壁溫度模擬結(jié)果可視為二維模型的鉆孔壁溫度；其導(dǎo)出的鉆孔壁暫態(tài)溫度序列可被植入二維模型的邊界條件函數(shù)。因此，在計(jì)算管群埋管區(qū)域某一個(gè)深度范圍內(nèi)的溫度分布時(shí)可以針對(duì)性的選取三維模型該范圍內(nèi)等間距的輸出結(jié)果作為輸入邊界條件帶入二維模擬，可避免運(yùn)行三維管群模型需要計(jì)算每一層網(wǎng)格所負(fù)擔(dān)的龐大計(jì)算量。換熱器模型建立過(guò)程中所作假設(shè)如下：1)地表溫度隨大氣溫度變化；2)管內(nèi)換熱流體溫度僅隨流動(dòng)方向一維變化；3)土壤及回填材料為物性均勻；4)材料物性參數(shù)在模擬過(guò)程中保持不變；5)管間的相互影響僅限于垂直于U形方向。

圖1 三維模型示意圖Fig.1 Schematic diagram for the 3D model

系統(tǒng)在運(yùn)行狀態(tài)下流體暫態(tài)的總換熱量等于因軸向流動(dòng)引起的熱量變化與徑向?qū)嵋鸬膿Q熱量之和，管壁與流體邊界上的能量平衡方程為：

(1)

圖2 二維管群模型示意圖Fig.2 Schematic diagram of the boreholes distribution of the site in Hangzhou

圖3 三維及二維網(wǎng)格模型Fig.3 The numerical mesh used to represent 3D and 2D model

系統(tǒng)停機(jī)狀態(tài)下流體暫態(tài)的總換熱量?jī)H受到徑向的導(dǎo)熱的影響，其平衡方程為：

(2)

在U形管，回填材料和土壤內(nèi)部的能量平衡方程為：

(3)

土壤的熱力學(xué)能初值設(shè)為零，變化方程為：

(4)

式中：T為溫度，K；τ為時(shí)間，s；m為質(zhì)量，kg；r為半徑，m；z為深度，m；u為流速，m/s；α為熱擴(kuò)散系數(shù)，m2/s；ρ為流體的密度，kg/m3；c為比熱容，J/(kg·K)；h為表面?zhèn)鳠嵯禂?shù)，W/(m2·K)；下標(biāo)w，p，f，s分別指代流體、U形管、回填材料以及土壤；下標(biāo)pi，po，b，in分別指代U形管內(nèi)、外管壁、鉆孔壁及入口坐標(biāo)位置。

假定各接觸面緊密結(jié)合，土壤表面的溫度，進(jìn)水口的溫度和速度均分別由溫度函數(shù)序列Tudf(τ)和速度函數(shù)序列uudf(τ)給出。邊界上的溫度關(guān)系及輸入條件為：

(5)

(6)

歷史氣候數(shù)據(jù)的獲取有兩種途徑，國(guó)家氣象局的網(wǎng)站或使用氣象站數(shù)據(jù)。如果都沒(méi)有，可以考慮建立氣溫的理想化周期波動(dòng)曲線：

T(τ)=0.5Tmcos(2×10-7τ+φm)+

0.5Tdcos(2×10-7τ+φd)+Tavg

(7)

式中：τ為時(shí)間，s；Tm，Td分別為年周期及日周期溫差,K；φm，φd分別為年周期及日周期的相位；Tavg為當(dāng)?shù)啬昶骄鶜鉁?，K。

二維模型傳熱方程的極坐標(biāo)表達(dá)式為：

(8)

式中：θ為角度，(°)；k為導(dǎo)熱系數(shù)，W/(m·K)。三維模型在運(yùn)行狀態(tài)下不同深度的孔壁溫度模擬結(jié)果可視為二維模型的鉆孔壁溫度，因而其導(dǎo)出的鉆孔壁暫態(tài)溫度序列Tudf′(τ)可構(gòu)成二維模型的鉆孔壁輸入邊界條件：

Ts(τ)|r=rb=Tudf′(τ)

(9)

同時(shí)三維模型導(dǎo)出的不同深度穩(wěn)定的土壤初始溫度構(gòu)成不同深度二維模型土壤初始溫度。

三維及二維模型模擬流程如圖4所示。

圖4 三維及二維模型模擬流程Fig.4 The simulation process of 3D and 2D model

2 換熱器模型的驗(yàn)證

為了驗(yàn)證本文采納的三維模型暫態(tài)模擬方法的準(zhǔn)確性，此處比照B.A.Beier等[20]提出的沙盒實(shí)驗(yàn)建立文獻(xiàn)中同尺寸的地下?lián)Q熱器暫態(tài)數(shù)值模型。模型采用的具體幾何與熱物性參數(shù)見(jiàn)表2，同時(shí)模型采用了與實(shí)驗(yàn)相同的入口流速和溫度作為輸入邊界條件。模擬預(yù)測(cè)的出水口暫態(tài)溫度文與文獻(xiàn)中傳感器實(shí)測(cè)的出水口溫度對(duì)比如圖5所示。其平均絕對(duì)比例誤差(MAPE)為3.13%，平均絕對(duì)誤差(MAE)為0.416，相較傳感器在測(cè)量時(shí)隨機(jī)產(chǎn)生的標(biāo)準(zhǔn)差0.391 K，可認(rèn)為所建立的三維模型能很好地反映實(shí)際情況。

圖5 模型預(yù)測(cè)出水溫度與實(shí)測(cè)出水溫度對(duì)比Fig.5 Comparison for predicted and measured outlet temperature

3 模擬結(jié)果與分析

論文采納了杭州市內(nèi)的氣象站數(shù)據(jù)，跨度為一個(gè)完整的供冷季、過(guò)渡季和供暖季。模型的模擬時(shí)間從2015年6月12日—2016年3月17日，期間分為供冷季、過(guò)渡季及供暖季。地源熱泵在供冷季和供暖季中分別向地下釋放熱量和冷量，在過(guò)渡季長(zhǎng)時(shí)間處于停機(jī)狀態(tài)。從15年6月12日開(kāi)始模擬直至16年3月17日的熱力學(xué)能與氣溫?cái)?shù)據(jù)如圖6所示。地下?lián)Q熱器模型采納的初始化條件為流體、U形管、回填材料和土壤各節(jié)點(diǎn)采用統(tǒng)一溫度設(shè)定為291.8 K，通過(guò)將第一步的邊界條件作為定值重復(fù)模擬50個(gè)假想日，可以有效緩沖初始化帶來(lái)的誤差影響。在模擬中第一個(gè)假想日的熱力學(xué)能增量為0.039 4 MJ，供冷季開(kāi)始前最后一個(gè)假想日的熱力學(xué)能增量?jī)H有0.001 9 MJ，熱力學(xué)能隨時(shí)間的斜率之比降至初始值的4.9%，近似于6月大氣條件下的土壤熱力學(xué)能自然增量的速度。表明經(jīng)過(guò)50個(gè)假想日的前置穩(wěn)定期后，土壤已基本達(dá)到適宜后續(xù)模擬運(yùn)算的溫度分布。

表2 沙盒實(shí)驗(yàn)參數(shù)

圖6 土壤熱力學(xué)能與氣溫對(duì)照數(shù)據(jù)Fig.6 Data of soil internal energy and air temperature

圖6中大氣溫度曲線在8月上旬(8月4日14∶00)達(dá)到峰值后開(kāi)始下行并在一月下旬(1月16日6∶00)達(dá)到谷值。地埋管在供冷季向地下釋放熱量并在供暖季從地下抽取熱量。對(duì)應(yīng)的土壤熱力學(xué)能曲線在供冷季上升，在過(guò)渡季中基本保持不變，并在供暖季開(kāi)始后逐步下降。與僅受氣溫影響的熱力學(xué)能對(duì)比可見(jiàn)系統(tǒng)運(yùn)轉(zhuǎn)下的內(nèi)能曲線變化在供冷季和供暖季明顯更加陡峭。在供冷季結(jié)束后整個(gè)三維計(jì)算區(qū)域的熱力學(xué)能上升了3.96 MJ，平均溫升為1.09 K。一個(gè)供冷季、過(guò)渡季、供暖季結(jié)束后，整個(gè)三維計(jì)算區(qū)域的熱力學(xué)能上升了0.292 MJ，與僅受氣溫影響的情況相比，土壤的相對(duì)熱力學(xué)能增量為1.037 MJ，相對(duì)溫度增量為0.28 K，如圖7所示。這反映了長(zhǎng)江三角洲地區(qū)每年對(duì)冷量的需求高于對(duì)熱量的需求，因而在年復(fù)一年的換熱過(guò)程后，將有大量的熱量被堆積在土壤中，進(jìn)而影響地埋管換熱器的換熱效率。

圖7 不同深度觀測(cè)點(diǎn)的土壤溫度Fig.7 The soil temperature at at different depth observation points

指定-0.5 m及-1.5 m的進(jìn)出水管連線中心軸上的觀察點(diǎn)，考察其全年的溫度變化情況。-1.5 m的觀測(cè)點(diǎn)的平均溫度及和僅受氣溫影響情況下的相對(duì)溫升峰值滯后于-0.5 m的觀測(cè)點(diǎn)約2周的時(shí)間。計(jì)算區(qū)域內(nèi)土壤的平均溫度和僅受氣溫影響情況下的相對(duì)溫升峰值滯后于-0.5 m的觀測(cè)點(diǎn)約9周的時(shí)間。說(shuō)明沿埋管深度深的絕對(duì)溫度和相對(duì)僅受氣溫影響情況下的溫升曲線均在相位上滯后于沿埋管深度較淺的點(diǎn)。

圖8 進(jìn)出水口溫度對(duì)比Fig.8 Comparison of the inlet fluid temperature and outlet fluid temperature

杭州市內(nèi)的地源熱泵系統(tǒng)在僅在需要供冷或供熱季節(jié)的工作時(shí)間開(kāi)機(jī)，在多數(shù)周末、節(jié)假日停機(jī)。尤其在9月末到11月初這段氣溫適宜的過(guò)渡季期間，系統(tǒng)有較長(zhǎng)的停機(jī)時(shí)段。在系統(tǒng)停機(jī)期間，土壤各點(diǎn)的溫度及進(jìn)出口處的水溫僅受大氣及環(huán)境溫度的影響，顯示出較為平滑的曲線。由圖8整理得知，當(dāng)系統(tǒng)運(yùn)行時(shí)，整個(gè)供冷季地埋管的出水口水溫與同時(shí)刻氣溫溫差均值為-8.313 K，整個(gè)供暖季地埋管的出水口水溫與同時(shí)刻氣溫溫差均值為9.077 K。

圖9所示為在完整的供冷季、過(guò)渡季、供暖季結(jié)束后埋管周?chē)臏囟确植记闆r，埋管區(qū)域一側(cè)及轉(zhuǎn)角處受到熱堆積的最大影響半徑見(jiàn)表3。由表3可知，沿埋管深度越深受熱堆積影響的范圍越小，埋管區(qū)域在轉(zhuǎn)角處的熱堆積影響半徑顯著大于埋管區(qū)域一側(cè)的影響半徑，其影響半徑的比值隨深度減小，反映出不規(guī)則排列地埋管群的熱堆積在較深的地層的影響范圍更局限于其自身的排列形狀。通過(guò)用系統(tǒng)運(yùn)轉(zhuǎn)條件下的截面平均溫度減去僅受大氣影響條件下的平均溫度可得不同深度截面的平均溫升，如圖10所示。-1 m深和-10 m以下的截面有不同的均溫變化規(guī)律，并且-1 m深度在整個(gè)供冷季的均溫峰值低于-10 m的均溫峰值。這一現(xiàn)象表現(xiàn)出與大氣的換熱在帶走土壤熱堆積尤其是表層土壤熱堆積方面的顯著影響。結(jié)合三維模型的模擬結(jié)果可知，系統(tǒng)運(yùn)轉(zhuǎn)所造成的熱堆積效應(yīng)主要集中在7 m以下稍深的土壤層而非與大氣頻繁換熱的表層土壤(最大溫差0.5 K)。

表3 不同深度處的最大影響半徑

圖9 供暖季結(jié)束時(shí)不同深度的二維溫度分布Fig.9 Temperature distribution at different depth at March 17th, 2016

圖10 截面平均溫升隨不同深度的變化Fig.10 Cross section average increment of temperatureat changes with different depth

3 結(jié)論

本文建立了基于實(shí)際氣候條件下的三維與二維結(jié)合地源熱泵地下?lián)Q熱器復(fù)合數(shù)值模型。實(shí)現(xiàn)了用較小的計(jì)算耗時(shí)處理模擬管群換熱的問(wèn)題。相較于現(xiàn)有模型，新建立的復(fù)合模型利用三維模型的暫態(tài)模擬結(jié)果作為二維模型輸入的邊界條件，進(jìn)而在降低管群運(yùn)算復(fù)雜度的同時(shí)保證了精確性。案例中，采用了杭州市內(nèi)的實(shí)測(cè)數(shù)據(jù)，模擬時(shí)間為2015年6月12日—2016年3月17日，跨度為一個(gè)完整的供冷季、過(guò)渡季和供暖季，主要結(jié)論如下：

1)當(dāng)計(jì)算三維換熱時(shí)，不需要采用設(shè)置初始溫度的方法，僅需設(shè)定初始?xì)鉁?，將第一步的邊界條件作為定值重復(fù)模擬50個(gè)假想日的時(shí)長(zhǎng)作為緩沖，可將熱力學(xué)能隨時(shí)間的日均斜率降至初始值的4.9%，有效沖抵了初始化時(shí)全部節(jié)點(diǎn)設(shè)置同一個(gè)初始溫度帶來(lái)的誤差干擾。

2)影響半徑為2 m以?xún)?nèi)的土壤在模擬一個(gè)完整的供冷季、過(guò)渡季和供暖季后的熱力學(xué)能增量為0.292 MJ，相較僅受氣溫影響的土壤的相對(duì)熱力學(xué)能增量為1.037 MJ，相對(duì)溫度增量為0.28 K。該增量反映出長(zhǎng)三角地區(qū)用戶對(duì)冷量的需求高于對(duì)熱量的需求，因此釋放到土壤中多余的熱量將隨時(shí)間逐步堆積，對(duì)換熱效率產(chǎn)生不利影響。

3)沿埋管深度深的絕對(duì)溫度和相對(duì)僅受氣溫影響情況下的溫升曲線均在相位上滯后于沿埋管深度較淺的點(diǎn)。在進(jìn)出水口連線中心軸上-1.5 m的觀測(cè)點(diǎn)的平均溫度及和僅受氣溫影響情況下的相對(duì)溫升峰值滯后于-0.5 m的觀測(cè)點(diǎn)約2周的時(shí)間。計(jì)算區(qū)域內(nèi)土壤的平均溫度和僅受氣溫影響情況下的相對(duì)溫升峰值滯后于-0.5 m的觀測(cè)點(diǎn)約9周的時(shí)間。

4)系統(tǒng)運(yùn)行條件下，整個(gè)供冷季地埋管的出水口水溫與同時(shí)刻氣溫溫差均值為-8.313 K，整個(gè)供暖季地源熱泵的出水口水溫與同時(shí)刻氣溫溫差均值為9.077 K。

5)系統(tǒng)運(yùn)轉(zhuǎn)造成的熱堆積效應(yīng)主要集中在7 m以下稍深的土壤層而非與大氣頻繁換熱的表層土壤。埋管區(qū)域在轉(zhuǎn)角處的熱堆積影響半徑大于其一側(cè)的影響半徑，且其絕對(duì)值與比值均隨深度的增大而減少，說(shuō)明不規(guī)則排列地埋管群的熱堆積的影響范圍在越深的地層越局限于自身的排列形狀。

[1] 中華人民共和國(guó)國(guó)家統(tǒng)計(jì)局. 2014年中國(guó)統(tǒng)計(jì)年鑒[M].北京: 中國(guó)統(tǒng)計(jì)出版社,2014：54-56. (National Bureau of Statics of China. China statistical yearbook in 2014[M]. Beijing: China Statistic Press, 2014：54-56.)

[2] INGERSOLL L R, PLASS H J. Theory of the ground pipe heat source for the heat pump[J].Heating Piping & Air Conditioning, 1948, 54(7): 339-348.

[3] MARCOTTE D, PASQUIER P, SHERIFF F, et al. The importance of axial effects for borehole design of geothermal heat-pump systems[J]. Renewable Energy, 2010, 35(4): 763-770.

[4] ZENG H Y, DIAO N R, FANG Z H. A finite line-source model for boreholes in geothermal heat exchangers[J]. Heat Transfer—Asian Research, 2002, 31(7): 558-567.

[5] MOLINA-GIRALDO N, BLUM P, ZHU K, et al. A moving finite line source model to simulate borehole heat exchangers with groundwater advection[J]. International Journal of Thermal Sciences, 2011, 50(12): 2506-2513.

[6] YANG H, CUI P, FANG Z. Vertical-borehole ground-coupled heat pumps: A review of models and systems[J]. Applied Energy, 2010, 87(1): 16-27.

[7] INGERSIOLL L, ZOBEL O J, INGERSOLL A C. Heat conduction: with engineering geological and other applications[M]. US: University of Wisconsin Press, 1954.

[8] KAVANAUGH S P. Simulation and experimental verification of vertical ground-coupled heat pump systems[D].US: Oklahoma State University Stillwater, 1985.

[9] CARSLAW H S, JAEGER J C. Conduction of heat in solids[M]. Oxford: Clarendon Press, 1960.

[10] REES S J, HE Miaomiao. A three-dimensional numerical model of borehole heat exchanger heat transfer and fluid flow[J]. Geothermics, 2013, 46: 1-13.

[11] ESKILSON P. Thermal analysis of heat extraction boreholes[M]. 1987.

[12] HE Miaomiao, REES S J, SHAO Li. Simulation of a domestic ground source heat pump system using a three-dimensional numerical borehole heat exchanger model[J]. Journal of Building Performance Simulation, 2011, 4(2): 141-155.

[13] HE Miaomiao. Numerical modelling of geothermal borehole heat exchanger systems[D].Leicester: De Montfort University, 2012.

[14] OZUDOGRU T Y, OLGUN C G, SENOL A. 3D numerical modeling of vertical geothermal heat exchangers[J]. Geothermics, 2014, 51: 312-324.

[15] LI Zhongjian. A new constant heat flux model for vertical U-tube ground heat exchangers[J]. Energy and Buildings, 2012, 45: 311-316.

[16] ZHANG Linfeng, ZHANG Quan, HUANG Gongsheng. A transient quasi-3D entire time scale line source model for the fluid and ground temperature prediction of vertical ground heat exchangers (GHEs)[J]. Applied Energy, 2016, 170: 65-75.

[17] MA Weiwu, LI Min, LI Ping, et al.New quasi-3D model for heat transfer in U-shaped GHEs (ground heat exchangers): Effective overall thermal resistance[J]. Energy, 2015, 90: 578-587.

[18] WETTER M, HUBER A. Vertical borehole heat exchanger EWS Model[J]. Trnsys Type, 1997, 451.

[19] GU Y, O′N(xiāo)EAL D L. Development of an equivalent diameter expression for vertical U-tubes used in ground-coupled heat pumps[J]. Ashrae Transactions, 1998, 104(2): 347-355.

[20] BEIER R A, SMITH M D, SPITLER J D. Reference data sets for vertical borehole ground heat exchanger models and thermal response test analysis[J]. Geothermics, 2011, 40(1): 79-85.