有滲流地埋管傳熱模型及快速算法＊

陳友明，張訓(xùn)水，杜立志，劉向偉

（湖南大學(xué) 土木工程學(xué)院，湖南 長沙 410082）

地源熱泵技術(shù)是一種高效清潔的空調(diào)方式，地埋管換熱器與周圍土壤間傳熱的研究是地源熱泵技術(shù)研究與應(yīng)用的關(guān)鍵.Eskilson[1]率先求得了單鉆孔地埋管換熱器在階躍熱流作用下的溫度響應(yīng)（即g－函數(shù)），對于多個(gè)鉆孔的地埋管換熱器所引起的溫度響應(yīng)可在單鉆孔傳熱分析的基礎(chǔ)上采用疊加原理進(jìn)行分析，這種方法已應(yīng)用在地源熱泵設(shè)計(jì)軟件EED[2]和GLHEPRO[3]中.Zeng等[4]繼承了這一方法，并根據(jù)線性疊加原理、格林函數(shù)法及虛擬熱源方法，進(jìn)一步求得有限長線熱源在半無限大介質(zhì)中溫度分布的解析解，這些顯式的函數(shù)關(guān)系式可以更加靈活方便地在計(jì)算程序中調(diào)用并計(jì)算.

實(shí)際工程中地埋管穿透的地層中或多或少都存在著地下水的滲流，尤其是沿海地區(qū)和地下水豐富的地區(qū)，地下水流動(dòng)的影響更為明顯，因此地埋管周圍土壤傳熱過程實(shí)際上是熱滲耦合過程[5].由于中深層（30～150m）地下水運(yùn)動(dòng)以水平方向?yàn)橹?，所以地下水橫向滲流的強(qiáng)弱成為對土壤傳熱的主要影響因素[6].然而，由于這一問題的復(fù)雜性，現(xiàn)有文獻(xiàn)鮮有對這一問題的深入分析.Eskilson[1]利用了Carslaw等[7]給出的移動(dòng)線熱源問題穩(wěn)態(tài)解析解，僅討論了在達(dá)到穩(wěn)態(tài)狀態(tài)以后滲流對地?zé)釗Q熱器的影響.Chiasson等[8]利用有限元法數(shù)值求解了二維滲流問題，對一些實(shí)際問題進(jìn)行了計(jì)算和比較，但未能找到較一般的規(guī)律[9].

文獻(xiàn)[10]根據(jù)多孔介質(zhì)中有滲流時(shí)的能量方程，首次解析求解得到了有均勻滲流時(shí)線熱源引起的二維非穩(wěn)態(tài)溫度響應(yīng)，并用實(shí)驗(yàn)方法進(jìn)行了驗(yàn)證，結(jié)果顯示其理論計(jì)算值與實(shí)測值吻合較好.由于地埋管換熱器涉及的范圍較廣，通常需要數(shù)年甚至更長的時(shí)間才能達(dá)到穩(wěn)定，因此實(shí)際應(yīng)用中必須考慮非穩(wěn)態(tài)問題，這使得文獻(xiàn)[10]給出的考慮水平滲流的地埋管換熱器有限長線熱源解析解模型計(jì)算量非常大，模擬耗時(shí)長.

本文提出了考慮地下水水平滲流的有限長線熱源模型新解析解，與文獻(xiàn)[10]相比，該解析解的理論計(jì)算值在保證合理計(jì)算精度的同時(shí)，顯著提高了模型計(jì)算速度，從而使有限長線熱源模型能更方便地應(yīng)用于實(shí)際工程中.

1 有滲流地埋管傳熱模型

考慮地下水水平滲流的地埋管傳熱過程涉及的幾何條件和物理?xiàng)l件非常復(fù)雜，為了簡化理論分析，做出如下假設(shè)：

1）地下巖土簡化為一個(gè)均勻的多孔介質(zhì)，熱量傳遞是由介質(zhì)（包括固體骨架和其中的流體）的導(dǎo)熱和孔隙中流體（水）的對流傳熱而實(shí)現(xiàn)的[8].

2）土壤為半無限大均勻介質(zhì)，在整個(gè)傳熱過程中土壤的熱物性不變.

3）忽略地表溫度波動(dòng)以及埋管深度對土壤溫度的影響，認(rèn)為土壤溫度均勻一致為T0.

4）忽略鉆孔的幾何尺度而把鉆孔近似為軸心線上有限長度的線熱源，并且每米換熱量均勻恒定.

在無滲流情況下，瞬時(shí)點(diǎn)熱源在無限大介質(zhì)中的溫度響應(yīng)（格林函數(shù)）在直角坐標(biāo)系中可以表示為[7]：

考慮有滲流的情況，比較多孔介質(zhì)中有滲流的能量方程和移動(dòng)熱源問題的熱傳導(dǎo)方程，可以看到二者有相同的形式，并且可以互相轉(zhuǎn)化 .文獻(xiàn)[10]考慮固定于 （x′，y′，z′）的點(diǎn)熱源，假設(shè)地下水沿平行于x軸正方向以速度u′移動(dòng)，由格林函數(shù)法和疊加原理可得到在τ時(shí)刻位于點(diǎn)（x，y，z）在階躍熱流q下的溫度響應(yīng)為：

式中：u為熱源等效移動(dòng)速度；u′為地下水滲流流速；α為地下巖土的熱擴(kuò)散系數(shù)；ρw為水的密度；cw為水的質(zhì)量比熱容；ρ為多孔介質(zhì)（包括水）的平均密度；c為多孔介質(zhì)（包括水）平均質(zhì)量比熱容；θ為過余溫度；T為所求溫度場；T0為地下巖土初始溫度.

g－函數(shù)是由Eskilson[1]提出的，定義為特定的單孔或者多孔地埋管換熱器的鉆孔壁在階躍熱流下的溫度響應(yīng).

式中：ks為巖土導(dǎo)熱系數(shù).

由式（3）可以看出，在單位管長換熱量q和巖土換熱系數(shù)ks確定的情況下，g－函數(shù)值可以確定過余溫度響應(yīng).

根據(jù)Eskilson的定義及公式（3），可得τ時(shí)刻在位于（x，y，z）的點(diǎn)上的g－函數(shù)[10]為：

式中：

引入虛擬熱源法來解決有限長線熱源模型地表面的恒溫邊界條件，即在與線熱源關(guān)于邊界面對稱的位置設(shè)一虛擬線熱匯，其強(qiáng)度與線熱源大小相等，符號相反，且虛擬介質(zhì)中也存在和實(shí)際介質(zhì)相同的均勻滲流，這樣，等溫邊界條件自動(dòng)得到滿足 .當(dāng)z＝H／2時(shí)，可得有限長線熱源模型鉆孔壁中點(diǎn)g－函數(shù)為：

式中：H為埋管深度.

文獻(xiàn)[5]指出，采用鉆孔壁深度方向中點(diǎn)處溫度作為鉆孔壁的代表溫度，雖然比較方便，但在定義上顯得有些隨意性.因此，采用溫度沿整個(gè)鉆孔壁深度方向的積分平均值作為鉆孔壁的代表溫度更為合理.

考慮水平滲流的有限長線熱源模型鉆孔壁積分平均溫度g－函數(shù)為：

2 新型g－函數(shù)

從式（6）中看到鉆孔壁積分平均g－函數(shù)是三重積分形式，這使得逐時(shí)模擬非常耗時(shí)，非常不利于工程應(yīng)用.本文通過調(diào)換參數(shù)積分次序的方法，提出一種計(jì)算速度較快的g－函數(shù)形式.由式（5）得[10]：

式中：

g－函數(shù)對時(shí)間τ的導(dǎo)數(shù)為：

式中：

式（8）對時(shí)間積分可得，有限長線熱源模型新的鉆孔壁中點(diǎn)溫度g－函數(shù)形式為：

分別對A，B，C沿鉆孔深度方向積分可得：

鉆孔壁積分平均溫度g－函數(shù)對時(shí)間τ的導(dǎo)數(shù)為：

式（10）對時(shí)間τ積分可得到鉆孔壁積分平均溫度g－函數(shù)：

3 計(jì)算精度與速度對比

本文提出的新型g－函數(shù)是在文獻(xiàn)[10]提出的g－函數(shù)基礎(chǔ)上利用數(shù)學(xué)方法改進(jìn)的，新的鉆孔中點(diǎn)及積分平均溫度g－函數(shù)值與文獻(xiàn)[10]模型的理論計(jì)算值完全相同.

為評估本文g－函數(shù)與文獻(xiàn)[10]模型的計(jì)算速率，本文分別對單鉆孔中點(diǎn)溫度g－函數(shù)及積分平均溫度g－函數(shù)的新舊2種解析解形式進(jìn)行了分析對比，所有計(jì)算采用Matlab編寫的程序完成，使用電腦性能為雙核1.80GHz.計(jì)算結(jié)果見表1和表2.

由表1和表2可知，新的鉆孔中點(diǎn)和積分平均g－函數(shù)計(jì)算速度較原g－函數(shù)均有了較大提高.模擬單鉆孔地埋管實(shí)際運(yùn)行1 000h，新的中點(diǎn)g－函數(shù)較原g－函數(shù)形式快了131倍；新鉆孔平均溫度g－函數(shù)式（11）較原g－函數(shù)式（6）在20h的計(jì)算中快了9 000倍.原積分平均溫度g－函數(shù)為三重積分的形式，計(jì)算速度過慢，非常不適于工程分析應(yīng)用.

表1 單孔中點(diǎn)溫度g－函數(shù)計(jì)算耗時(shí)比較Tab.1 Comparison of the calculation time for g－function（middle－temperature）of one borehole

表2 單孔平均溫度g－函數(shù)計(jì)算耗時(shí)比較Tab.2 Comparison of the calculation time for g－function（mean－temperature）of one borehole

圖1為滲流對鉆孔中點(diǎn)溫度及平均溫度g－函數(shù)的影響.由圖1可知，有限長線熱源的g－函數(shù)長期運(yùn)行均趨于穩(wěn)定.無滲流的g－函數(shù)均大于有滲流的g－函數(shù)，且無論有無滲流，中點(diǎn)溫度g－函數(shù)均高于積分平均溫度g－函數(shù).這表明在地埋管長時(shí)間模擬中，采用鉆孔壁中點(diǎn)溫度使得鉆孔壁溫度計(jì)算值偏大，這將導(dǎo)致地埋管選型過大，從而使地源熱泵系統(tǒng)的初投資過大.

圖1 滲流對鉆孔中點(diǎn)溫度及平均溫度g－函數(shù)的影響Fig.1 Effect of advection on the g－function of borehole area

4 多鉆孔計(jì)算

對存在多個(gè)鉆孔的鉆井域進(jìn)行分析計(jì)算時(shí)，我們采用Eskilson給出的g－函數(shù)疊加方法進(jìn)行計(jì)算.當(dāng)存在N個(gè)鉆孔時(shí)，任一鉆孔中點(diǎn)和平均溫度g－函數(shù)可分別表達(dá)為[1]：

式中：rb為鉆孔半徑；

下面考察4個(gè)鉆孔在不同滲流速度及滲流方向的情況下鉆孔壁沿鉆孔深度方向積分平均溫度g－函數(shù)及溫度響應(yīng).鉆孔域內(nèi)鉆孔布置及地下水水平滲流方向見圖2，計(jì)算結(jié)果見圖3和圖4.

圖3為不同滲流速度下的g－函數(shù)比較，其中u＝5.4×10－7m·s－1，rb＝0.055m .圖中表明存在水平滲流的g－函數(shù)值較無滲流情況小，且滲流速度越大，g－函數(shù)值越小并能更快地趨于穩(wěn)定.

圖2 鉆孔布置及滲流方向圖Fig.2 Schematic of borehole layout and advection directions

圖3 不同滲流速度下的g－函數(shù)比較Fig.3 Comparison of g－functions under different advection velocity

圖4 不同滲流方向下的g－函數(shù)比較Fig.4 Comparison of g－function under different advection directions

圖4為鉆孔壁1在不同滲流方向下的g－函數(shù)比較，由圖4可知，鉆孔壁沿深度方向積分平均溫度g－函數(shù)，在長時(shí)間運(yùn)行后滲流方向1，2g－函數(shù)值小于滲流方向3，4；當(dāng)趨于穩(wěn)定后滲流方向4的g－函數(shù)比滲流方向2大18.15%.這表明若計(jì)算鉆孔處于地下水滲流方向下游，則g－函數(shù)值較大，相應(yīng)的鉆孔壁過余溫度絕對值較大；若計(jì)算鉆孔處于地下水滲流方向上游，則g－函數(shù)較小，相應(yīng)的鉆孔壁過余溫度絕對值較小.

圖5為在階躍熱流作用下，滲流速度為u＝5.4×10－7m／s，滲流方向1時(shí)鉆孔域的溫度場（即g－函數(shù)值分布）.圖6為在同樣階躍熱流作用下，無滲流時(shí)鉆孔域的溫度場.圖中表明有水平滲流的鉆孔域溫度分布與無滲流情況有較大不同，有滲流鉆孔域溫度分布不對稱，且其溫度場明顯小于無滲流情況.這說明，水平滲流能緩解在階躍熱流長期作用下引起的熱量積累，從而改善鉆孔域的換熱工況.

圖5 有地下水滲流時(shí)鉆孔域溫度場Fig.5 Temperature field of borehole area with groundwater advection

圖6 無地下水滲流時(shí)鉆孔域溫度場Fig.6 Temperature field of borehole area without groundwater advection

5 結(jié) 論

考慮地下水滲流的地埋管換熱模型更接近于實(shí)際的應(yīng)用工況，g－函數(shù)計(jì)算是考慮地下水滲流模型應(yīng)用的關(guān)鍵，在鉆孔單位管長熱流及巖土傳熱系數(shù)ks確定的情況下，g－函數(shù)值唯一決定了鉆孔域的溫度場.針對文獻(xiàn)[10]模型鉆孔壁中點(diǎn)及積分平均溫度g－函數(shù)計(jì)算繁瑣耗時(shí)的問題，本文提出了新的鉆孔壁中點(diǎn)和積分平均溫度g－函數(shù)形式，改進(jìn)后的g－函數(shù)均為一維積分的形式，在保證計(jì)算精度的同時(shí)，計(jì)算速度有了顯著提高.積分平均溫度能準(zhǔn)確地反應(yīng)鉆孔壁溫度，在實(shí)際工程應(yīng)用中可采用新的積分平均溫度g－函數(shù)進(jìn)行鉆孔壁溫度響應(yīng)計(jì)算.

另外，本文還運(yùn)用疊加原理分析了有水平滲流時(shí)4鉆孔方形區(qū)域階躍熱流下鉆孔壁溫度響應(yīng).計(jì)算結(jié)果表明，在階躍熱流作用下，水平滲流能顯著改善鉆孔域的換熱工況，且滲流速度越大，效果越明顯.對于特定的采樣點(diǎn)，其溫度響應(yīng)與地下水滲流方向有關(guān)，處于地下水滲流方向下游的鉆孔壁過余溫度絕對值較大，處于地下水滲流方向上游的鉆孔壁過余溫度絕對值較小 .鉆孔分布形式對多埋管區(qū)域溫度場有較大影響.由于在鉆孔埋管區(qū)域主要存在地下水的水平滲流，所以在實(shí)際工程中采用考慮地下水水平滲流的地埋管換熱模型能更精確地進(jìn)行地埋管設(shè)計(jì)，防止地埋管選型過大.

[1] ESKILSON P.Thermal analysis of heat extraction boreholes[D].Lund，Sweden：Department of Mathematical Physics，U－niversity of Lund，1987：58－69.

[2] HELLSTROM G，SANNER B.Earth energy designer：software for dimensioning of deep boreholes for heat extraction[C]／／Energy Conservaton through Energy Storage Stock Conference CALORSTOCK 94′.Finland Espoo：Energy Conservaton through Energy Storage，1994：185－192.

[3] SPITLER J D.GLHEPRO—a design tool for commercial building ground loop heat exchangers[C]／／Proceedings of the Fourth International Heat Pumps in Cold Climates Conference.Quebec：Oklahoma State University，2000：132－138.

[4] ZENG H，DIAO N，F(xiàn)ANG Z.A finite line－source model for boreholes in geothermal heat exchangers[J].Heat Transfer Asian Research，2002，31（7）：558－567.

[5] 刁乃仁，方肇洪.地埋管地源熱泵技術(shù)[M].北京：高等教育出版社，2006：47－68.DIAO Nai－ren，F(xiàn)ANG Zhao－h(huán)ong.Ground－coupled heat pump technology[M].Beijing：Higher Education Press，2006：47－68.（In Chinese）

[6] 蔡晶晶，陳汝東，王健.地下水滲流對地埋管傳熱影響的理論分析[J].流體機(jī)械，2009，37（12）：62－68.CAI Jing－jing，CHEN Ru－dong，WANG Jian.Analysis of effects of groundwater advection on geothermal heat exchanger[J].Fluid Machinery，2009，37（12）：62－68.（In Chinese）

[7] CARSLAW H S，JAEGER J C.Conduction of heat in solids[M].2nd ed.Oxford：Clarendon Press，1947：258－264.

[8] CHIASSON A D，REES S J，SPITLER J D.A preliminary assessment of the effects of ground water flow on closed－loop ground source heat pump systems[J].ASHRAE Transactions，2000，106（1）：380－393.

[9] 刁乃仁，李琴云，方肇洪.有滲流時(shí)地?zé)釗Q熱器溫度響應(yīng)的解析解[J].山東建筑工程學(xué)院學(xué)報(bào)，2003，18（3）：1－5.DIAO Nai－ren，LI Qin－yun，F(xiàn)ANG Zhao－h(huán)ong.An analytical solution of the temperature response in geothermal heat exchangers with groundwater advection[J].Journal of Shandong University of Architecture and Engineering，2003，18（3）：1－5.（In Chinese）

[10] 刁乃仁.地?zé)釗Q熱器的傳熱問題研究及其工程應(yīng)用[D].北京：清華大學(xué)航空航天學(xué)院，2004：57－64.DIAO Nai－ren.Heat transfer analyses of ground heat exchangers and their engineering applications[D].Beijing：School of Aerospace，Tsinghua University，2004：57－64.（In Chinese）