三維快速多極邊界元法分析地埋管群傳熱問題

宋子欣, 胡宗軍, 胡 斌, 牛忠榮

(合肥工業大學 土木與水利工程學院, 合肥 230009)

0 引 言

地埋管換熱器的埋置深度一般在30～200 m之間,埋管管徑為25～100 mm,壁厚為2～6 mm,U型管道狹長,涉及到薄壁結構與周圍土體耦合問題,導致換熱分析較為困難．為此,國內外學者對其傳熱問題進行了諸多研究．Guan等[1]將U型管通過等效公式簡化為矩形管,略去了U型管管底區域,在此模型下討論了地下水滲流速度對埋管傳熱的影響．Kerme等[2]利用能量平衡方程分析了具有兩個獨立回路的雙U管的瞬態傳熱過程,得到了沿鉆孔深度分布的傳熱曲線．賀澤群等[3]采用自適應負荷法,以鉆孔壁溫度替代管內壁溫度,對不同鉆孔位置的換熱量和熱堆積問題進行了研究．張荻等[4]采用實驗和數值模擬相結合的方法,研究了層流條件下球窩結構內矩形通道內部的流動與換熱特性．朱利媛等[5]考慮熱量沿軸向的傳遞,利用邊界元法對二維穩態熱傳導下的地埋管換熱器進行了模擬研究．Lei等[6]將U型管簡化為三維曲線,利用有限元軟件COSMOL建立了地埋管群的熱滲耦合三維數值模型,研究了系統運行一年后地層溫度場的變化．

目前,多數學者是利用有限元法(FEM)對地埋管換熱器進行模擬研究．在FEM中有無限單元、管道單元和基于連續體的殼單元等,用于計算薄體和厚體耦合的傳熱問題．其中無限單元一般用于解決土壤無限域問題,和常規有限元一起用來解決更復雜的無界問題．同時,對于模擬接觸問題,基于連續體的殼單元由于考慮了厚度的變化,與常規殼單元相比更加精確．相比于FEM,邊界元法(BEM)只需在結構的表面和界面上劃分網格,縮減了網格的自由度數量,因此在分析薄體問題時具有優勢．但傳統BEM生成的定解方程組的系數矩陣是非對稱滿秩矩陣[7],計算量級過大,難以計算大規模問題,因此傳統BEM無力計算地埋管群的傳熱問題．隨著快速多極算法的發展,應運而生的快速多極邊界元法(FMBEM)提高了傳統BEM的計算效率,并降低了內存的占有量,近二十年被廣泛應用于解決各種大規模問題．由于FMBEM在使用線性單元和高階單元時,八叉樹結構中存在同一單元的節點處于不同葉子的現象,導致單元積分難以處理．因此目前FMBEM主要采用常值單元,在處理薄體結構和高梯度場時,需布置稠密網格以彌補計算精度的不足,導致計算效率較低．徐剛等[8]將傳統常值單元中直接布置在流體計算域表面網格中心上的奇點移到了計算區域外部,以此實現了無奇異化．劉靜等[9]基于自適應分塊技術,提出了一種組合變換法,分別消除了徑向和角度方向積分的近奇異性．侯俊劍等[10]利用擴展單元插值法,在不連續的邊界配置虛擬節點,利用虛擬節點插值邊界上連續和不連續的物理場,將非連續單元變為高階連續單元．胡宗軍等[11]針對二維位勢問題,通過在幾乎奇異積分單元上扣除奇異函數部分消除了幾乎奇異性,并將其應用于二維薄體結構溫度場分析中．

本文基于邊界元法基本理論,引入三維線性單元幾乎奇異積分的正則化算法[12],解決了線性單元跨葉子積分難題,對三維位勢問題建立了三角形線性單元FMBEM和計算程序,大大提高了FMBEM分析大規模三維薄壁熱傳導問題的計算效率和精度．文中采用三角形線性單元FMBEM對地埋管換熱器的換熱性能進行分析,討論了管壁厚度對地埋管傳熱性能的影響,并研究了地埋管群的傳熱效率．本文研究成果可為地源熱泵的設計和運行提供參考．

1 三維位勢問題快速多極邊界元分析與數值驗證

1.1 邊界元法分析穩態熱傳導的基本方程

(1)

(2)

(3)

式中r=|y-x|為源點y到場點x的距離．

采用3節點三角形面單元對計算域邊界進行離散．在直角坐標系Ox1x2x3中建立如圖1所示的局部參考坐標系oξη,對單元Γe的幾何坐標、位勢以及熱流進行插值,單元Γe上場點x的坐標xi可以表達為

圖1 參考坐標系Fig. 1 Reference coordinate systems

(4)

式中形函數Nm為

(5)

其中A為三角形單元面積,am,bm,cm為常數,由3節點局部坐標表示．

ξ-ξ0=ρcosθ,η-η0=ρsinθ．

(6)

將式(6)代入式(4),有

(7)

(8)

因此,源點y與場點x之間的距離r可表示為

(9)

式中

(10)

用邊界單元離散后,邊界積分方程表達式(1)變為

(11)

式中f為邊界單元總數．讓源點y遍歷所有邊界節點,則邊界積分方程(11)轉化為代數方程

Az=B,

(12)

其中A為系數矩陣,z為未知量,B為已知量．引入已知邊界條件,即可由方程(12)求得所有未知的邊界節點位勢u(x)和熱流q(x)．上述方程若采用直接法或迭代法求解,如網格數量眾多,則系數矩陣A需要極大的存儲量,計算量也快速增長．

這里的快速多極邊界元法將式(11)中線性單元積分劃分為兩類： 遠場單元積分和近場單元積分．近場單元積分是指源點離單元較近的積分, 包含有奇異積分、 幾乎奇異積分和非奇異積分．本文采用Gauss數值積分計算非奇異積分; 采用常位勢場法計算奇異積分, 而對于幾乎奇異積分采用半解析算法[12], 具體過程如下．

將式(8)和式(9)代入式(2)、(3),離散后的積分形式為

(13)

隨著式(10)中e1的減小,式(13)的積分產生幾乎奇異性,常規Gauss數值積分失效．對式(13)的變量ρ做積分,令

(14)

并對式(14)反復運用分部積分,因Qn(ρ,θ)為ρ的多項式,在三維位勢問題中,于線性單元而言,對ρ求導必有

(15)

因此,反復運用分部積分后,式(13)轉化為

(16)

其中ρ(θ)由單元Γe的3個邊的極坐標表達式決定．式(16)已將式(13)的面積分轉化成一系列關于變量θ的線積分,積分的幾乎奇異性消失,可利用常規Gauss數值積分計算．

對于遠場單元積分,FMBEM通過對基本解(2)和(3)進行雙球諧函數展開,用自適應樹結構代替傳統矩陣,減少了計算量和存儲量．基本解U*(x,y)在球坐標系下的多極展開式[13]為

(17)

(18)

(19)

將基本解Q*(x,y)按同樣方法展開,有

(20)

(21)

由式(18)—(21)可以看出,場點x由中心點yc傳遞到源點y,避免了每個源點都要對每個場點進行單獨積分計算,因此FMBEM在形成式(12)系數矩陣A時大幅度減少了計算量．

目前FMBEM主要采用常值單元[14],因而沒有幾乎奇異積分計算的難題．本文將三角形面單元幾乎奇異積分處理技術[15]與三維邊界元快速算法結合,創立了三維位勢問題三角形線性面單元FMBEM,并研制了計算程序．本文三角形線性面單元FMBEM算法顯著提高了三維位勢邊界元法的計算精度和效率,能夠勝任大規模三維薄體結構的熱傳導分析．

1.2 多介質傳熱模型FMBEM分析

采用三角形線性面單元FMBEM程序,分析圖2所示雙層圓筒壁結構的熱傳遞問題．r1=8 mm,r2= 10 mm,r3=11 mm,圓筒沿Z方向深20 mm,Ω1,Ω2兩個區域的導熱系數分別為9.93 W/(m·℃)和0.5 W/(m·℃),圓筒最內側壁面溫度T=100 ℃,最外側壁面熱流q=21.33 W/m2,其余表面絕熱,熱流q=0 W/m2．

圖2 雙層圓筒壁結構Fig. 2 The double cylinder wall structure

對結構表面和界面按照三角形面單元進行離散,如圖3所示,為做對比分析,分別給出了3種網格模型： ① 具有10 327個節點、20 864個單元的稀疏網格模型; ② 具有16 036個節點、32 324個單元的常規網格模型; ③ 具有36 793個節點、74 004個單元的加密網格模型．

圖3 雙層圓筒壁網格模型Fig. 3 Double-layer cylindrical wall mesh models

圖4給出深度Z=10 mm位置,沿AB路徑的3種網格模型通過FMBEM程序計算所得溫度TFMBEM以及與解析解溫度Ta的相對誤差Δ,其中r為AB路徑上的點到圓心的距離：

圖4 FMBEM在Z=10 mm處沿AB路徑計算所得溫度TFMBEM及相對誤差ΔFig. 4 Temperature TFMBEM and relative error Δ calculated with the FMBEM along the AB path at Z=10 mm

對稀疏網格模型,本文算法計算得到的各點溫度與解析解相比,在r大于9.9 mm后相對誤差Δ增大,但不超過2.1%;對常規網格模型和加密網格模型,本文算法的溫度計算結果基本一致,且與解析解的相對誤差Δ不超過1.05%．計算結果表明本文算法計算三維熱傳導問題時,采用常規網格即可獲得較高的計算精度．

2 U型地埋管換熱器的FMBEM分析

2.1 U型管計算假定及設計參數

地源熱泵系統中,U型地埋管換熱器和土壤之間的實際傳熱關系復雜,現采用三角形面單元FMBEM程序計算其換熱性能．本文對地源熱泵系統做以下假設：

1) 將土壤看成連續均勻的各向同性物質,土壤的熱物性在整個換熱過程中保持不變．

2) 假設土壤地下溫度保持恒定,不受地面氣溫的影響[16]．

3) 水的熱擴散率比土壤低2～4個數量級[17],因此在地埋管換熱器的傳熱過程中,忽略土壤中水分和溶質的轉移,即不考慮地下水滲流的影響．

4) 忽略管道內流體和地埋管管內壁之間的對流換熱熱阻,假設管內壁溫度為流體溫度,且溫度沿管埋深線性變化．在U型管底部區域溫度恒定,為進水與出水溫度的平均值．

計算模型中,土壤區域和回填區域上表面均絕熱,熱流q=0 W/m2,其余設計參數如表1所示．

表1 地埋管換熱器設計參數

在U型地埋管換熱器中,管壁厚度遠小于回填區域和土壤區域的尺寸,一般計算方法難以考慮U型管管壁厚度對換熱的影響,而邊界元法具有降維特性,FMBEM只需在結構表面和界面上劃分網格,特別適用于厚體和薄體區域的耦合問題．在制冷和制熱兩種工況下,本小節研究了不同的管壁厚度對單U型地埋管換熱量的影響．取管外直徑為32 mm,管材為聚乙烯,導熱系數為0.4 W/(m·℃),其他計算參數如表2所示．U型地埋管的FMBEM計算模型見圖5,共有23 113個節點和46 256個三角形面單元．

表2 不同管壁厚度時地埋管單位井深換熱量(單位： W/m)

采用FMBEM計算得到各點的溫度及熱流密度．然后根據U型管內壁的熱流計算地埋管換熱器單位井深換熱量：

(22)

其中Si為U型管內壁上第i個三角形單元Гi的面積,qi為FMBEM計算所得單元Гi平均熱流密度(三角形單元3個節點熱流均值),m為U型管內壁面所劃分的單元總數,L為地埋管埋置深度．

不同管壁厚度時,通過FMBEM程序計算得地埋管單位井深換熱量如表2所示,其中Qi,Qo分別為進、出水管單位管長換熱量,Q=Qi+Qo為單位井深換熱量．由表2可知,隨著管壁厚度增加,地埋管單位井深換熱量逐漸減小．制冷條件下,當管壁厚度為10 mm時,單位井深換熱量為28.944 W/m,與管壁厚度為3 mm時的單位井深換熱量41.381 W/m相比,減小了30.05%;與不考慮管壁厚度時的單位井深換熱量46.696 W/m相比,減小了38.02%．制熱條件下,當管壁厚度為10 mm時,單位井深換熱量為17.739 W/m,與管壁厚度為3 mm時的單位井深換熱量25.361 W/m相比,減小了30.05%;與不考慮管壁厚度時的單位井深換熱量28.638 W/m相比,減小了38.06%．

按我國對地埋管管道外徑及壁厚尺寸的規定[18],當管外徑小于40 mm時,壁厚應為3 mm．由表2可知,制冷條件下,管壁厚度為3 mm時,單位井深換熱量為41.381 W/m,與不考慮管壁厚度時的單位井深換熱量46.696 W/m相比,減小了11.38%;制熱條件下,當管壁厚度為3 mm時,單位井深換熱量為25.361 W/m,與不考慮管壁厚度時的單位井深換熱量28.638 W/m相比,減小了11.44%．

為反映管壁厚度對換熱的影響,圖6分別給出制冷、制熱兩種工況下,不考慮管壁厚度和管壁厚度為10 mm時,地埋管回填土壤區域在底部邊界的溫度分布局部圖．由圖6可知,當管壁導熱系數小于土壤和回填材料的導熱系數時,管壁熱阻是阻礙地埋管換熱的主要因素之一,討論地埋管換熱器的換熱量時,不應忽略管壁厚度對其產生的影響,否則會造成較大誤差．當管壁導熱系數一定時,管壁越厚,對管內流體和土壤之間的換熱影響越大,計算時考慮管壁厚度產生的影響,會得到更符合工程實際的結果．當管壁導熱系數小于土壤導熱系數時,在保證U型管結構強度滿足要求的前提下,管壁應盡量薄,以減小其熱阻對地埋管換熱器產生的負面影響．

(a) 制冷工況時FMBEM計算結果溫度分布圖(a) Calculation results of temperature during cooling

2.2 多孔管群換熱器的FMBEM分析

在工程應用中,地源熱泵地埋管換熱器系統通常都是以多管井即管群的形式出現．由于各埋管之間存在熱相互作用,因此管群的換熱特性與單U型地埋管相比有明顯區別,單U管的結論不能直接應用于大型管群換熱分析．本文采用FMBEM對地埋管換熱器管群系統進行分析,討論管群之間的熱干擾現象．

考慮4(橫向)×4(縱向)等間距排列管群在不同工況下的運行情況, 利用對稱性取1/4結構, 以分析含16個管井的U型管群換熱性能．取土壤半徑為15 m, 土壤區域內分布有4組U型管井, 各管井中心間距為5 m,井深50 m．模型坐標原點位于對稱中心,靠近坐標原點的管井,其中心距兩坐標軸均為2.5 m,為①號管井．②號管井坐標位于(-7.5 m, 2.5 m),③號管井位于(-2.5 m, 7.5 m),④號管井位于(-7.5 m,7.5 m)．

本例中鉆孔區域內部采用原土回填,即土壤和回填土的導熱系數均為2.0 W/(m·℃)．U型管外徑32 mm,管壁厚3 mm,兩管腿中心間距50 mm．利用三角形面單元對結構表面和界面進行離散,共生成58 785個節點,117 670個單元,如圖7所示．土壤上表面與對稱面絕熱,外表面溫度18 ℃,其余計算參數如表2所示．利用FMBEM程序計算得制冷、制熱兩種工況下各管井單位井深換熱量,如表3所示．其中Qi,Qo分別為進、出水管單位管長換熱量,Q=Qi+Qo為單位井深換熱量．

表3 壁厚3 mm時4×4管群單位井深換熱量(單位： W/m)

由于此管群模型有16個管井,呈對稱形式,因此將管井分為3類,中心的4個管稱為中井,離中心最遠位置的4個管稱為角井,除了角井之外的四周的管稱為邊井．

由表3可知,對于4×4地埋管管群,制冷工況下,中井單位井深換熱量為9.431 W/m, 與角井的單位井深換熱量23.855 W/m相比, 減少了60.47%; 制熱工況下, 中井單位井深換熱量為6.434 W/m, 與角井的單位井深換熱量14.627 W/m相比, 減少了55.78%．邊井與角井相比,單位井深換熱量在制冷工況下分別減少了32.61%和27.57%,在制熱工況下分別減少了31.36%和27.49%．中井位于土壤中心, 周圍存在的其他管井同土壤的換熱導致土壤中心區域溫度變化, 減小了中井內流體與其周圍土壤的溫差, 降低了中井的換熱效率．圖8給出了4×4管群在制冷、制熱兩種工況下的溫度云圖,由溫度云圖可以看出,在兩種工況下,角井內流體與周圍土壤的溫差均大于中井．因此,管井與管井之間產生了熱量的互相干擾,周圍管井對中心管井產生的熱干擾阻礙了中心管井沿徑向傳熱．因角井周圍毗鄰無限大土壤區域,熱擴散條件優于其他區域管井,受到的熱干擾影響最小;邊井附近除了無限大土壤區域,還受到角井和中井的影響,熱擴散條件較差,導致其換熱量下降;中井被其他管井所包圍,換熱量最低,在夏季工況下只有角井換熱量的39.53%,冬季工況也僅達到了角井換熱量的44.22%．

(a) 制冷工況時FMBEM計算結果溫度分布圖(a) The FMBEM calculation results of temperature during cooling

為直觀地比較管群間的熱相互影響,圖9給出了在制冷、制熱兩種工況下,單U型地埋管與①至④號管井及4×4管群單位井深換熱量的對比．由圖9可知,制冷工況下,與單U型地埋管單位井深換熱量40.607 W/m相比,4×4管群的單位井深換熱量為266.556 W/m,提高了556.43%;制熱工況下,與單U型地埋管單位井深換熱量24.883 W/m相比,4×4管群的單位井深換熱量為166.828 W/m,提高了570.99%．中井在制冷工況下的單位井深換熱量為9.431 W/m,相較于單U型地埋管換熱量降低了76.77%;在制熱工況下的單位井深換熱量為6.434 W/m,相較于單U型地埋管換熱量降低了74.01%．角井在制冷工況下的單位井深換熱量為23.855 W/m,比單U管換熱量降低了41.25%,制熱工況下的單位井深換熱量為14.627 W/m,比單U管換熱量降低了41.22%．

圖9 4×4地埋管群單位井深換熱量對比Fig. 9 Comparison of heat transfer fluxes of 4×4 buried pipe groups

由此可見,對于地埋管管群而言,周邊區域埋管的熱擴散條件顯著優于中心區域埋管,在管內流體溫度相同時,周邊區域埋管承擔了較大的換熱負荷．中心區域埋管受到周邊區域埋管熱干擾的影響較大,導致土壤中心區域地埋管換熱器換熱效率低下．

為研究管壁厚度對管群換熱量的影響,表4給出了4×4管群在不考慮壁厚時換熱量的計算結果,其中管徑為32 mm．表4結果表明,制冷工況下,不考慮壁厚時4×4管群總換熱量為289.016 W/m,與壁厚3 mm時4×4管群總換熱量266.556 W/m相比,提高了8.43%;制熱工況下,不考慮壁厚時4×4管群總換熱量為193.060 W/m,與壁厚3 mm時4×4管群總換熱量166.828 W/m相比,提高了15.72%．因此,對管群進行模擬計算時,考慮管壁熱阻十分必要．若忽略管壁熱阻產生的影響,會造成較大的計算誤差．表3、4對比表明,換熱量增加主要是由于中井和角井的貢獻．因此,于管群而言,提高換熱量的主要措施在于降低埋管之間的熱相互作用,同時減少中井和角井的管壁厚度．

表4 不考慮壁厚時4×4管群單位井深換熱量(單位： W/m)

3 結 論

本文基于線性單元,克服單元跨葉子積分難題,針對邊界元法在奇點附近區域因存在奇異積分而產生的誤差,將三維位勢問題快速多極邊界元法與幾乎奇異積分的半解析算法相結合,創立了一種精確高效的快速多極邊界元新算法,使得邊界元法能夠有效分析大規模三維薄體結構熱傳導問題．應用本文算法對考慮了管壁厚度的地埋管換熱器三維模型進行計算,分析討論了單U型地埋管換熱器和管群換熱器的傳熱性能．所得結論如下：

1) 管群中,管與管之間產生的熱干擾會阻礙熱量的傳遞,導致中心區域管井換熱效率低下,且管群規模越大,管與管之間熱干擾現象越強烈,提高管群換熱量的主要措施是降低熱干擾的同時減少壁厚．

2) 管壁的厚度是影響地埋管換熱效率的主要因素之一．在討論地埋管換熱器的換熱量時,對管壁的影響因素加以考慮,方可得到更符合工程實際的結果．

3) 本文中對曲面結構表面采用平面單元進行網格劃分,若采用曲面單元劃分網格,則有可能進一步減小問題的自由度數,同時提高計算精度．