中深層U型埋管換熱器數值傳熱模型精準度影響因素研究

蔣向明，朱余良，郭曉凱，鮑玲玲，李振興，郭海明

（1.中國煤炭地質總局水文地質局，河北 邯鄲 056004；2.河北工程大學 能源與環境工程學院，河北 邯鄲056038）

0 引言

地熱資源具有儲能豐富、清潔環保、穩定可靠、受季節和環境變化影響小等特點[1]，充分地開發和利用地熱資源，對實現我國“雙碳”目標具有重大意義。

地熱能資源根據埋藏深度可分為3類：地表以下200 m以內為淺層地熱能；200～3 000 m為中深層地熱能；3 000 m以下為深層地熱能[2]。我國陸地和海域的大地熱流分布呈現為西南部（平均熱流值為90 mW/m2）高于全球平局值（平均熱流值為65 mW/m2）[3]，東部持平，中部（平均熱流值為55 mW/m2）和西北部（平均熱流值為50 mW/m2）偏低的格局。受到能源分布以及現有地下勘探、人工造儲和深層鉆井技術的限制，我國深層地熱能開采難度較大[4]。同時，淺層地埋管地源熱泵系統因其占地面積大、全年冷熱負荷平衡要求高、長期運行后系統性能顯著下降等原因發展受限[5]，因此中深層地熱資源的開發利用逐步成為研究的焦點。中深層地源熱泵系統取熱的核心部件是井下地埋管換熱器，換熱器主要分為同軸套管式地埋管換熱器及U型對接井地埋管換熱器。國內外學者對同軸套管式地埋管換熱器的傳熱機理、取熱特性、系統運行模式、影響因素等進行了分析[6]～[9]，并已在工程上得到了驗證[10]。

目前，少有對中深層U型埋管換熱器數值傳熱模型模擬計算條件以及簡化幾何模型精準度影響的相關研究，本文通過建立的中深層U型埋管換熱器數值傳熱模型，定量分析了控制方程的坐標變化系數等模擬條件變化對模擬結果的影響，同時對數值計算所用的不同簡化幾何模型的精準性進行了計算、對比，為中深層U型對接井換熱器數值計算的模型研究和工程計算提供了參考。

1 傳熱模型

1.1 幾何模型

U型對接井地埋管換熱器因其優秀的取熱能力而受到關注[11]，其結構如圖1所示。

由于受到鉆井工藝的限制，它由一個垂直井和一個L型對接井構成，L型井圓弧連接段的彎曲半徑由鉆井工藝和鉆井直徑共同決定。受限于鉆井工藝，中深層U型對接井的設計施工幾何尺寸并不能直接用于數值計算，需要對其進行合理而必要的簡化，本文沿用國內學者[12]針對幾何模型的簡化方式進行模型簡化。

1.2 數學模型

對模型進行合理假設是模型建立與求解的關鍵，對整體換熱系統作出如下假設。

①下行管、水平管與上行管3部分及其換熱邊界內巖土體之間不存在熱干擾。

②下降管、水平管和上升管內流體的流動和傳熱均視為一維導熱問題，不考慮流體在管道橫截面上的速度和溫度分布，只強調平均溫度和速度。

③循環流體、固井材料的熱物性不變；下行管、水平管與上行管周圍換熱區域內同層巖土體物性均勻、一致、熱物性不變；不考慮換熱區域內巖土體中地下水滲流對換熱的影響。

④換熱器套管外壁面、回填材料以及井孔壁面三者緊密接觸，忽略其接觸熱阻。

⑤數值模擬計算中設置的換熱邊界足夠遠，邊界外巖土體溫度恒定。

⑥計算區域每一位置的熱流量恒定。

⑦井孔內固井材料及循環流體的溫度一致，軸向導熱相對于徑向導熱可忽略。

⑧換熱器中循環水沒有泄漏或其它質量交換，下降管、水平管、上升管的流速相同。

根據假設條件，將換熱區域分為井孔內和井孔外兩部分，井孔內傳熱可視為一個線熱源或者是柱熱源在無限大均勻介質中的非穩態傳熱過程，其控制方程為

式中：C為井眼通道的熱容量，包括管內的流體、管壁和回填材料，J/（m3·K）；Tf為循環流體溫度，℃；Tb為地熱井壁面巖土溫度，℃；t為時間，s；M為單位時間流經截面的流體質量，kg/s；z為軸向坐標；R為循環流體和井孔壁之間的傳熱熱阻，K/W；ri，ro，ris，rb分別為套管的內半徑、外半徑、保溫材料外半徑和井壁外半徑，m；λp，λis，λb分別為套管、保溫材料及固井材料的導熱系數，W/（m·K）；cf，cp，cb分別為循環流體、管道壁及固井材料的體積比熱容，J/（kg·℃）；ρf，ρp，ρb分別為循環流體、管道壁及固井材料的密度，kg/m3；h為循環流體與管壁的對流換熱系數，W/（m2·K）。

巖土體之間的傳熱可以簡化為圓柱坐標系下的二維非穩態導熱問題，控制方程為[13]

式中：α為土壤的熱擴散系數，m2/s；T為巖土溫度，℃；r為徑向坐標。

考慮到換熱過程中巖土體徑向溫度梯度與其到井壁的距離成反比，在邊界處趨近于零，因此，沿徑向考慮采用不等距網格劃分的方法，在近地熱井孔壁面處網格尺寸盡量要小，而足夠遠處的邊界網格可以大一些，從而減小計算的數量，提高計算速度。具體方法如下：

式中：x為新坐標系下x向坐標；rc為坐標變換基數，m。

則新坐標系下控制方程轉化為

在新坐標系下控制方程中x向取步長Δx劃分網格，同時令

式中：r1，r2···rn為從管壁開始的第n個網格徑向長度，m；β為坐標變換比例系數。

原圓柱坐標系中r向為不等距網格劃分，即：

1.3 邊界條件及參數選擇

在地表邊界溫度設定以及巖土體初始溫度計算過程中，一些學者以地表空氣年平均溫度[13]或地表年平均溫度[14]為基礎進行初始溫度場計算，本文認為應以巖土體恒溫層溫度為基礎進行初始溫度場計算，主要原因如下。

①巖土體恒溫層溫度基本不受外部環境影響，可認為常年恒定，是理想的邊界層。

②冬季采暖項目以全年地表平均溫度作為計算標準，在時間區間上不符合實際情況。

③地表年平均溫度和恒溫層之間的土壤溫度梯度與恒溫層以下的土壤溫度梯度方向是否一致，大小是否在可接受的范圍內，需要進一步考證。

實驗項目位于河北省邯鄲市，圖2為該地區淺層地下巖土體全年溫度分布。

圖2 邯鄲地區淺層地下巖土體全年溫度分布Fig.2 Annual temperature distribution of shallow underground rock and soil mass in Handan

由圖2可以看出，淺層巖土體恒溫層出現在地下15 m處[15]，恒溫層溫度為15.7℃。因此，本文中數學模型以地表下15 m作為模型上邊界，以地埋管換熱器下100 m作為模型下邊界，徑向距地埋管換熱器中心100 m為右邊界，均為第一類邊界條件。

2 數值求解

按實驗項目的實際情況，整理得到實驗區域基本參數（表1），并以此作為討論基準。

表1 系統參數Table 1 System parameter

表1中，第1層巖土體為0～420 m，第2層巖土體為420～1 040 m，第3層巖土體為1 040～1 540 m，第4層巖土體為1 540～2 290 m，第5層巖土體為2 290～2 600 m。

2.1 模型驗證

利用交替方向隱式（Alternating Dorectopm Implicit，ADI）格式建立離散方程，采用追趕法直接對節點方程進行求解。通過與中科院汪集旸院士團隊領銜開發并完成的“地熱計算器”軟件模擬結果對比發現，本文數值模型計算結果與地熱計算器計算的結果在換熱器溫度分布及換熱器出口循環水溫度變化趨勢上是完全一致的，換熱器出口循環水溫度差為0.47℃，驗證了數值計算模型、計算方法以及Matlab程序的可靠性。

2.2 計算條件分析

首先以坐標變換比例β=1.2，rc=rb，徑向節點數N=40、軸向步長dz=20 m、時間步長dt=180 s為基準進行計算，以模型網格的第15，35，65，95和125層進行分析，結果如圖3所示。

圖3 運行不同時長下不同土層徑向節點位置、編號與溫度對應關系Fig.3 The corresponding relationship between the position and number of radial nodes in different soil layers and temperature under different operating conditions

圖4為連續運行不同時長下,徑向遠邊界節點號與巖土體深度、溫度對應關系。

圖4 連續運行24個月，徑向遠邊界節點編號與巖土體深度、溫度對應關系Fig.4 The corresponding relationship between the node number of radial far boundary and the depth and temperature of rock and soil mass when the system operates continuously for 24 months

由圖4可知，不同物性巖土體徑向遠邊界不同，徑向遠邊界在420，1 040，1 540，2 280 m附近會發生躍遷，徑向邊界最遠處的節點為N=40。這是由于上述4處均位于巖層物性變化的位置，其導熱系數及體積熱容等參數突變導致的。同時，巖土體的徑向溫度遠邊界隨連續運行時間的增加而增加。

2.3 模擬條件對計算精度的影響分析

為進一步確定換熱器出口溫度與巖土體徑向遠邊界之間的關系，針對巖土體徑向節點的數量（N=15，20，25，26，27，28，29，30，35，40，45）變化進行了1個采暖季（2 880 h）的模擬，結果如表2所示。

表2 徑向巖土體節點數及遠邊界對溫度模擬計算精度的影響Table 2 Influence of radial rock and soil node number and far boundary on temperature simulation calculation accuracy

由表2可知，換熱器出口循環水溫度隨徑向遠邊界的擴大而減小，最終趨于穩定。當N＜25時，換熱器出口溫度隨其變化波動較大；當N＞25時，換熱器出口溫度隨其變化波動微小；當N≥29時，換熱器出口溫度穩定于19.820℃，取熱量為1 030 kW。巖土體分層界面的突變以及徑向遠邊界附近（26≤N≤32）巖土體傳熱對換熱器出口溫度幾乎沒有影響。針對計算模型中時間步長dt變化進行了模擬計算，結果見表3。由表3可知，當dt≤1 500 s時，換熱器循環水出口溫度的模擬結果隨著時間步長的增加基本不變；當dt≥1 800 s時，下行管中循環水出口溫度即出現振蕩。得出在模型的模擬過程中，出口溫度隨時間步長的增加而發生不收斂的狀況可能是突變的，這一點將在后續的研究中進一步驗證。

表3 時間步長對模擬計算精度的影響Table 3 Influence of time step on simulation accuracy

針對計算模型中垂直管、水平管軸向步長變 化進行了模擬計算，結果見表4，5。

表4 平管軸向步長對模擬計算精度的影響Table 4 Influence of horizontal tube axial step size on simulation calculation accuracy

表5 垂直管軸向步長對模擬計算精度的影響Table 5 Influence of vertical tube axial step size on simulation calculation accuracy

由于管道軸向步長與每層巖土體的厚度相關，為方便比較，垂直管部分管道軸向步長僅取值10，20 m，同時將水平管的管道計算總長設置為600 m。可見對于垂直管，軸向步長dz從10 m增加至20 m，上行管的出口溫度降低約0.4%，但是模擬計算耗時減少70%。對于水平管，隨著dz的增加，水平管的出口溫度同樣逐漸降低，降低幅度小于0.02%/10 m。

針對計算模型中坐標變化比例系數β的變化進行了模擬計算，結果見表6。由表6可知，坐標變化比例系數β在一定范圍內波動時，換熱器出口溫度隨β取值增加而緩慢降低。但當坐標變化比例系數β取值較小時，時間步長的取值相應的也需要減小，否則計算結果離散。

表6 坐標變化比例系數對模擬計算精度的影響Table 6 Influence of coordinate change proportion coefficient on simulation calculation accuracy

綜上所述，為提高計算速度，同時保證模擬計算結果的準確性和可靠性，推薦相近項目的模型研究計算參數：坐標變換比例β=1.20，軸向步長dz=20 m，時間步長dt=1 200 s，徑向遠邊界為55～65 m；工程設計的計算參數：坐標變換比例β=1.40，軸向步長dz=20 m，時間步長dt=1 200 s，徑向遠邊界為15～20 m。

活動是德育課程實施的重要載體，也是學生主體性生成和發展的源泉。《新課程標準》指出：“課程設計與實施注重聯系學生的生活實際，引導學生在實踐中發現和提出問題，在親身參與豐富多樣的社會活動中，逐步形成探究意識和創新精神。”這從課程理念的高度，充分肯定了活動是德育課教學之必需。隨著課標的深入落實，活動成了德育課教與學的中介，在課堂上教師通過活動來創設教學情境，使抽象的概念具體化，使傳授的方式趣味化，使認知與情感融合化，使思維與形象統一化。活動符合小學生的認知特點，有利于學生主體人格的形成，從而促進學生知、情、意、行的全面發展。

3 簡化模型精準度分析

3.1 現有模型分析

針對前文數值計算依據的幾何模型，本文認為地埋管換熱器計算管長、水平管與下行管、上行管連接處及其周邊的熱干擾是影響模擬計算結果的重要因素。

3.1.1地埋管換熱器計算管長

受限于鉆井工藝以及地下巖土狀況，L型對接井中必然存在的圓弧連接段，而數值模型中將L型井部分直接簡化為下行管+水平管進行計算，導致數值模型計算得出的循環水出口溫度、地熱井提取熱量等結果高于真實值。

3.1.2水平管與下行管、上行管連接處及其周邊的熱干擾

通過前文對傳熱區域進行邊界分析，模擬運行1個采暖季時，水平管與下行管、上行管連接處及其周邊巖土體換熱區域內存在熱干擾，如圖5所示。

圖5 數值計算模型水平管與垂直管連接處熱干擾Fig.5 Thermal interference at the connection between horizontal and vertical pipes in numerical calculation model

因此，數值模型增加的下行管與水平管連接處以及原垂直井與L型對接井連接位置的無熱干擾加熱，均增加了理論計算的循環水出口溫度。

3.2 簡化模型

圖6為依據圖1幾何模型的簡化過程。

圖6 幾何模型簡化過程及尺寸Fig.6 Geometric model simplification process and dimensions

3.2.1簡化幾何模型A

已知本模型中距離管中心15 m以內的巖土中的溫度變換占整體溫度變化的99%以上。因此針對L型對接井的圓弧連接段，以沿管道軸向20 m長度段為單元，以15 m為半徑繪制了換熱單元斷面圖，見圖7。由于圓環內外側斷面弧長與換熱器管道單元長度相比，變化均小于4%，所以將此圓環體換熱單元視為圓柱體換熱單元。模擬過程中以1個軸向步長即20 m的長度差分別對內外弧長的圓柱取熱區域進行了計算，換熱器出口循環水溫度與以管中心長度計算的出口溫度差約為±0.065℃，取熱量差約為±6.8 kW，占比0.7%，偏差可接受。

圖7 圓環體換熱單元斷面Fig.7 Cross-section view of heat exchange unit of torus

3.2.2簡化幾何模型B

在簡化模型A的基礎上，忽略拐點處沿換熱器軸線且遠離換熱器方向巖土體的傳熱，將整個換熱器視為一根直管，得到簡化模型B。

3.2.3簡化幾何模型C，D

L型對接井本身并沒有連接拐點，圖1中的數值計算尺寸不僅增加了計算管長，同時增加了計算拐點，但是由于兩處拐點循環流體溫度不同，這將對其與周圍巖土體的換熱產生影響，因此分為兩步對其簡化，得到簡化模型C，D。

3.3 條件分析

為便于描述本章節中涉及模型的初始條件和邊界條件，以換熱器入口為起點，以模型A，C中換熱器拐點、模型B，D換熱器出口為終點的管道長度作為縱坐標進行描述：

其中：

對于模型A，C：

對于模型B，D：

式中：Tg為巖土體恒溫層溫度，℃；L為管道計算深度，m；L1為物理模型中下行管豎直管段長度，m；L2為物理模型中下行管豎直管段與弧形管道長度之和，m；L3為物理模型中下行管豎直管段、弧形管段與水平管段長度之和，m；TL1，TL2，TL3分別為從入口起至L1，L2，L3長度后的流體溫度，℃；rL為物理模型中水平管與下行管之間圓弧連接的半徑，m；Lh為圓弧連接段之后與上行管之間水平管的長度，m；Lm為管道總長，m。

當τ≥0，rb≤r≤rbd，L=0時，T（r，L，τ）為巖土體地表溫度邊界條件，當τ≥0，rb≤r≤rbd，L=Lm時，T（r，L，τ）為巖土體地底溫度邊界條件，當τ≥0，rb=rbd，T（r，L，τ）為巖土體徑向換熱遠邊界溫度邊界條件。

3.4 模擬結果比較

為了便于比較，A～D 4組模型依然以坐標變換比例β=1.2、rc=rb、徑向節點數N=35、軸向步長dz=20 m、時間步長dt=180 s為基準進行計算，其它參數按表1執行，模擬結果見表7。

表7 不同幾何模型的數值計算結果Table 7 Numerical results of different geometric models

由表7可知，模型A，B循環水溫度數值計算結果在水平管出口及上行管出口相差小于0.01℃，說明拐點附近熱擾動在現有條件下對模擬結果的影響確實可以忽略不計。但是在水平管出口處，模型A的循環水溫度略低于模型B，而在上行管出口處，模型A的循環水溫度略高于模型B。這是因為當存在拐點時，拐點既參加了水平管的計算也參加了垂直管的計算，水平管出口溫度處拐點僅參與了水平管部分計算，因此模型A水平管出口溫度略低于模型B，但是經過上行管進口拐點計算后，至上行管出口處，循環水溫度反而會高于模型B。模型A，B循環水溫度數值計算結果在水平管出口及上行管出口處均低于圖1模型，相差約0.52℃，結合表1參數，取熱量差約為54 kW，約占總體取熱量的5%。模型C，D與圖1中數值計算模型的模擬結果同樣驗證了拐點處的熱擾動對換熱器出口處循環水溫度影響可以忽略這一假設，同時與模型A，B的影響趨勢一致。

L型井連接處彎曲半徑的大小既與鉆井工藝及井下地質條件相關，同時影響整個系統的取熱性能，因此針對其變化進行了模擬，得出L型井連接段計算彎曲半徑依次為150，200，300，400 m時，循環水提取熱量依次為1 007.2，1 001.1，983.5，970.5 kW，與圖1模型計算所得的循環水提取熱量比較依次減少17.5，23.6，41.2，54.2 kW，兩種模型取熱量差值隨著彎曲半徑的增加而增加。

4 結論

①中深層地埋管換熱器與其周圍巖土體的傳熱數值計算中，當巖土體上邊界（表面邊界）采用第一類邊界條件時，邊界溫度應選用巖土體恒溫層溫度作為邊界條件；當巖土體上邊界（表面邊界）采用第三類邊界條件時，以中深層恒定熱流與表面對流換熱系數之比納入土壤溫度場計算是不合適的。巖土體的主要溫度變化發生在近管區域。即使隨著運行時間的增加，巖土體的換熱遠邊界不斷擴大，但是在實驗參數下同深度巖土體溫度變化的99%始終集中在距離管中心15 m范圍內。

②通過一系列模擬計算分析，給出了適當鉆井直徑下數值計算的基本參數。推薦相近項目的模型研究參數：坐標變換比例β=1.2，軸向步長dz=20 m，時間步長dt=1 200 s，徑向遠邊界rbd=55～65 m；工程設計參數：坐標變換比例β=1.4，軸向步長dz=20 m，時間步長dt=1 200 s，徑向遠邊界rbd=15～20 m。應用數值模型對換熱系統關于熱干擾假設進行了定量計算，垂直管與水平管連接拐點處的熱干擾對換熱器出口循環水溫度的影響小于0.01℃。

③將中深層U型對接井換熱器中L型井部分的圓弧連接段直接投影成垂直井+水平井進行數值計算是不合適的，這直接導致計算管長增加了0.429rL（實驗項目增加了172 m）。因此，在實驗參數下，針對L型井圓弧連接段以圓柱換熱單元等效圓弧換熱單元進行數值計算，其計算結果與將其投影成垂直井+水平井的數值計算結果取熱量相差54 kW，約占總體取熱量的5%，偏差較大。因此當L型井連接處彎曲半徑≥200 m時，推薦以幾何模型A，B進行中深層地埋管換熱器數值計算。