低溫模式下地下換熱器管土結構凍脹變形特性研究

王有鏜，鄭 斌，孔喜磊，李成宇，王春光

(山東理工大學交通與車輛工程學院，山東 淄博 255049)

近些年，地下蓄能以及地源熱泵技術已在許多領域與地區得到了廣泛應用[1-5]，作為換熱系統中的重要組成部分，地下換熱器的低溫凍脹變形問題一直制約著其在嚴寒地區的適用性[6-8]。在諸如油氣儲運[9-10]、供水[11]、渠道[12]和樁基[13]等行業同樣會出現土壤凍脹管路變形的問題，針對該問題，許多研究工作已經開展。在研究過程中，雖然原位工程的測試試驗研究結果真實可靠[14]，但限于工程規模、周期和測試手段，對于更為系統深入的研究，利用實驗方法顯然難以實現。因此，從20世紀70年代開始，逐步出現了很多描述巖土凍脹過程的數值模型，用以進行更為細致的研究，較為著名的有水力模型[15]、剛性冰體模型[16]、分凝勢模型[17]和熱力模型[18]等。其中，熱力模型主要研究凍土整體冰晶生長與分布，本文所采用模型正是基于熱力模型的改進。

對于換熱結構凍脹變形的研究目前仍處于初級階段，因此本文通過建立地下換熱器管土結構凍脹變形過程的數值模型，在驗證模型可行性的基礎上，對地下換熱器管土結構凍脹變形的基本特性展開了初步研究[19]。對于實際應用，換熱管內流體溫度的變化取決于地上供熱負荷或蓄能需求的變化，其溫度運行工況必然呈現多樣性，然而不同的溫度運行工況將通過影響凍結溫度場，從而對地下換熱器管土結構的凍脹變形帶來不同的影響。為此，本文以單U型地下換熱器管土結構為研究對象，模擬其不同恒溫溫度和不同進、出水溫差工況下的凍脹變形過程，分析凍脹土體應力、管體變形和管截面的變化特征，為合理選擇地下換熱器的運行方式以減輕凍脹危害提供參考依據。

1 地下換熱器管土結構凍脹變形過程數值模型的建立

1.1 控制方程

土壤凍脹擠壓埋管變形的過程主要用土壤孔隙增長率函數、凍土應力-應變關系以及相變傳熱方程進行描述。

(1) 土壤孔隙增長率函數：凍土體積膨脹速率可用孔隙增長率函數進行描述，文獻[18]已對其進行驗證。孔隙增長率N可表達為：

(1)

(2) 凍土應力-應變關系：可把凍土變形視作彈性變形，彈性應變增量和凍結過程中孔隙增長引起的應變增量構成了總應變增量。凍脹引起的總應變增量中的各個分量可表達為：

(2)

式中:1、2、3為由熱流方向所確定直角坐標系的坐標方向;σ為土體的應力(Pa)；E為土體的彈性模量(Pa)；μ為泊松比；ξ為土體的凍脹各向異性系數；γ為土體的剪切應變；τ為土體的剪切應力(Pa)；G為土體的剪切模量(Pa)，G=E/[2×(1+μ)];t為時間(s)。

(3) 凍土相變傳熱：基于能量守恒定律和導熱微分方程，水在固液相變過程中釋放出潛熱的導熱方程表達式如下：

(3)

1.2 計算流程

本文基于ABAQUS仿真平臺對地下換熱器管土結構的凍脹變形過程進行數值模擬。先利用ABAQUS軟件中CAE模塊進行幾何建模及參數條件設置；再利用Standard模塊中的溫度-位移分析步計算傳熱和應力場。對于凍土變形過程，ABAQUS軟件本身并沒有對其進行描述的本構模型，因此需將該本構模型編輯為子程序并關聯ABAQUS軟件的主程序。首先采用Fortran語言設計編寫子程度UMATHT計算凍土相變傳熱方程[式(3)]，獲取溫度場分布，再編寫UEXPAN子程序求解土體孔隙增長率函數[式(1)]，確定土體孔隙率增量；然后將子程序編輯為for格式文件，在ABAQUS/Command程序中直接輸入指令使求解器調用子程序，以此完成凍變本構方程計算，同時實現主程序和子程序計算數據的耦合，見圖1。對土壤凍脹過程中的孔隙變化、相變傳熱以及土顆粒、水、冰各組分體積變化依次進行迭代循環，最終完成凍土溫度場、冰含量、換熱結構位移和管體變形等參數的計算。

圖1 基于ABAQUS軟件的地下換熱器管土結構凍脹變形過程的計算流程Fig.1 Calculation flow of frost heave process of pipe soil structure of underground heat exchanger based on ABAQUS software

1.3 數值模型與參數

U型彎管以外的大部分換熱器管段沿軸向沒有結構差異，因此本文應用二維結構來模擬單U型地下換熱器管土結構的凍脹變形，同時為減少計算時間，結合管土結構和熱量傳遞的對稱特征，建立半圓形幾何結構，如圖2(a)所示。在該半圓形幾何結構中，將其半徑R1設置為5 m，兩換熱管半徑R2設置為16 mm，管壁厚度設置為3.5 mm，管間距L設置為60 mm。忽略管土結構的初始應力，將其初始溫度場設置為4℃，將模型的外邊界進行限位設置，并給予恒溫4℃的邊界條件，相應地，X軸上的邊界設置為對稱絕熱邊界。模擬采取分區法劃分模型單元網格，對于溫度梯度較大的換熱管及其周圍土壤區域，增加該重點區域單元網格數量以提高計算精度，而對于剩余土壤區域則可適當降低網格數量，網格單元采用雙線性平面四節點位移-溫度耦合結構(CPE4T)，如圖2(b)所示。為衡量凍結區域范圍，定義X向凍結直徑dx為沿X軸方向上兩凍結鋒面之間的距離，定義Y向凍結直徑dy為沿Y軸方向上兩凍結鋒面之間的距離，如圖2(c)所示，同時定義凍結區平均直徑dxy為dx和dy的平均值，以此衡量凍結區發展。

圖2 單U型地下換熱器管土結構數值模型Fig.2 Model of pipe soil structure of the single U-shaped underground heat exchanger

地下換熱器管土結構模型中換熱管采用高密度聚乙烯材料，凍結土壤由水、冰和土顆粒組成，其基本參數設定見表1。

表1 地下換熱器管土結構數值模型的基本參數

已針對該數值模型在土體凍結溫度場和換熱管變形應變兩方面進行了試驗驗證，詳見參考文獻[19]。

2 結果與分析

2.1 不同恒溫工況對地下換熱器管土結構凍脹變形的影響

針對恒溫工況，本文模擬了不同低溫運行模式對管土結構凍脹變形的影響，模擬設置100 h的運行時間，對進、出水管內壁分別施加三種運行模式HW1、HW2和HW3，其運行溫度見表2。

表2 地下換熱管恒溫運行模式

2.1.1 土體凍結區變形

為研究運行溫度對管土結構凍脹變形的影響，本文基于相同的管土換熱量，即在相同凍結范圍條件下對比分析各模式凍脹變形的情況。圖3為不同恒溫運行模式下巖土平均凍結直徑dxy的發展曲線。

圖3 不同恒溫運行模式下土體平均凍結直徑dxy隨 運行時間的發展曲線Fig.3 Curves of average freezing diameter (dxy) of soil mass with operation time under different operation modes of constant temperature

由圖3可見，三種恒溫運行模式下土體平均凍結直徑dxy值均發展至337 mm時，HW1、HW2和HW3三種恒溫運行模式分別耗時100 h、36 h和23 h，可發現運行溫度越低，傳遞相同換熱量的用時也越短。

圖4為不同恒溫運行模式下管間區域土體特征點O處的X向應力隨土體平均凍結直徑的dxy變化曲線。

圖4 不同恒溫運行模式下管間區域土體特征點O處的 X向應力隨土體平均凍結直徑dxy的變化曲線Fig.4 Curves of soil feature point stress in X direction at O point between the underground heat exchange pipes with average freezing diameter (dxy) under different operation modes of constant temperature

由圖4可見，隨凍結直徑發展，O點的壓應力不斷增大，壓應力在三種模式間呈現出明顯且逐漸增大的差異，表現為HW3

2.1.2 管體變形

模擬過程中換熱管的直徑變化量見圖5。

由圖5可見，換熱管發生了X向直徑減小、Y向直徑增大的橢圓化變形。隨凍結直徑發展，換熱管截面變形不斷增大，換熱管X向直徑變化量要大于Y向，出水管直徑變化量略大于進水管。三種模式下的管直徑變化量表現為HW3

圖5 不同恒溫運行模式下換熱管直徑的變化量隨 運行時間的發展曲線Fig.5 Curves of pipe diameter variation with operation time under different operation modes of constant temperature

2.1.3 管截面變化

土體凍脹擠壓換熱管發生橢圓化變形，其管截面流通面積隨之改變，由此可對循環產生阻力。為分析管體變形對其流通面積的作用，定義管體流通截面變化率Cs，可表示為：

(4)

式中:Si和Sd分別為管體變形前、后管截面的流通面積(m2),其中Sd值可由橢圓形面積公式Sd=abπ計算得到[a、b分別為管體變形后橢圓形截面的短軸半徑(m)和長軸半徑(m)]。

當Cs>0時，表明管體變形使管截面的流通面積減小；當Cs<0時，則表明管體變形使管截面的流通面積增大。

根據公式(4)，可計算得出不同恒溫運行模式下的Cs值，并繪制三種恒溫運行模式下Cs值隨土體平均凍結直徑的變化曲線，見圖6。

圖6 不同恒溫運行模式下管截面當量流通面積變化 率(Cs)隨土體平均凍結直徑dxy的變化曲線Fig.6 Curves of equivalent flow area variation ratio of pipe section with average freezing diameter (dxy) under different operation modes of constant temperature

由圖6可見，可以發現Cs值均大于0，因此Cs值可視為流通截面減小率，隨著凍結區發展，該值基本呈線性增大趨勢，管體變形使其流通截面不斷減小。盡管進水管的管徑變化量略小于出水管，但前者的Cs值卻稍大于后者，因此即使管徑變化量稍大，其流通截面減小率不一定大。

隨著凍結區域發展，Cs值的增大幅度依次為HW1>HW2>HW3。當凍結直徑增至150 mm時，HW1和HW2運行模式的Cs值小于HW3運行模式，而當凍結直徑達到200 mm時，HW2和HW3運行模式的Cs值開始小于HW1運行模式，根據發展趨勢可以推斷，HW3運行模式的Cs值逐漸會小于HW2運行模式。因此對于低溫運行工況，當凍結范圍較小時，低溫運行模式的管體流通截面減小程度要大于高溫運行模式，而當凍結范圍超過一定值時，高溫運行模式減小程度將超過低溫運行模式。

2.2 不同進、出水溫差對地下換熱器管土結構凍脹變形的影響

隨著埋設深度增加，單U型地下換熱器的進、出水管溫差會逐漸減小。為考察換熱器不同埋設深度的凍脹變形差異，本文通過設置不同的進、出水管溫差，進行模擬分析。設置三種溫差運行模式DT1、DT2和DT3，各運行模式均以10℃的線性降溫運行100 h，其溫差分別為1℃、2℃和3℃，如表3所示。

表3 進、出水溫差運行模式

由表3可知，三種運行模式下進、出水兩管的平均運行溫度相同，均為從0℃線性降至-10℃。

2.2.1 土體凍結區變形

不同進、出水溫差(DT1、DT2和DT3)運行模式下隨運行時間的發展土體平均凍結直徑dxy的變化曲線見圖7。

圖7 不同進、出水溫差(DT1、DT2和DT3)運行模式下 土體平均凍結直徑dxy隨運行時間的發展曲線Fig.7 Curves of average freezing diameter of rock with operation time under different operation modes of inlet and outlet temperature difference

由圖7可見，在低溫運行100 h后，不同溫差模式的dxy值均達到380 mm，各曲線基本重合，這是由相同的管體平均運行溫度所致。

本文以管間區域土體特征點O處的應力變化來衡量凍脹力的變化特性，得到不同進、出水溫差(DT1、DT2和DT3)運行模式下，土體特征點O處X向應力隨運行時間的變化曲線，見圖8。

圖8 不同進、出水溫差(DT1、DT2和DT3)運行模式 下管間區域土體特征點O處X向應力隨運行 時間的發展曲線Fig.8 Curves of soil feature point stress in X direction at O point between the underground heat exchange pipes with operation time under different operation modes of inlet and outlet temperature difference

由圖8可見，三種不同溫差模式下O點位置的壓應力有所差異。可以發現， DT1運行模式的凍脹力最小，DT2運行模式略大，DT3運行模式最大，運行最終，該位置的壓應力分別達2.31 MPa、2.37 MPa和2.44 MPa，溫差增大2℃，可使壓應力增大0.13 MPa。由此可見，較大溫差的運行模式可在凍結區產生較大的凍脹力。

2.2.2 管體變形

降溫過程中，不同進、出水溫差(DT1、DT2和DT3)運行模式相比，進水管與出水管的直徑表現出不同的變形特性。圖9為不同進、出水溫差(DT1、DT2和DT3)運行模式下進、出水管在X向的管體直徑減小量和Y向的管體直徑增大量變化曲線。

圖9 不同進、出水溫差(DT1、DT2和DT3)運行模式下 換熱管管體直徑的變化量隨運行時間的發展曲線Fig.9 Curves of pipe diameter variation with operation time under different operation modes of inlet and outlet temperature difference

由圖9可見，相比之下，出水管的變形量更為明顯；對于進水管，管徑變化表現為DT3

2.2.3 管截面變化

根據公式(4)，可計算得出不同進、出水溫差運行模式下Cs值，并繪制三種運行模式下管截面當量流通面積變化率隨運行時間的發展曲線，見圖10。

圖10 不同進、出水溫差(DT1、DT2和DT3)運行模式下 管截面當量流通面積變化率(Cs)隨運行時間的 發展曲線Fig.10 Curves of equivalent flow area variation ratio of pipe section with operation time under different operation modes of inlet and outlet temperature difference

由圖10可見，Cs值均大于0且隨運行時間不斷增大，變形使得管體流通截面不斷減小。對于進水管[見圖10(a)]，三種運行模式的Cs值在運行初期有所不同，表現為DT1

3 結 論

本文通過建立地下換熱器管土結構凍脹變形過程數值模型，以0℃以下運行的單U型地下換熱器管土結構為研究對象，模擬分析了不同恒溫溫度和不同進、出水溫差對管土結構凍脹變形的影響，得到以下結論：

(1) 對于恒溫運行模式，在土體凍結范圍基本一致的情況下，高溫將比低溫運行模式可產生更大的凍脹力，且其差異隨土體凍結范圍的增大而增大，當土體凍結直徑發展至150 mm和337 mm時，運行溫度增加10℃，可使土體凍脹應力分別增大0.31 MPa和0.81 MPa；與之對應，高溫運行可產生更大的管體變形，運行溫度增加10℃，可有0.55 mm的管體變形增量；當土體凍結范圍較小時，低溫比高溫運行模式有更大的管截面減小程度，而當土體凍結范圍超過200 mm時，高溫運行模式下管截面減小程度將大于低溫運行模式。

(2) 對于平均運行溫度一致的降溫模式，不同的進、出水溫差可對換熱器管土結構凍脹變形產生一定的影響：進、出水溫差增大2℃，可使土體凍脹應力增大0.13 MPa；進、出水溫差增大2℃，可使進水管管體變形量減少0.57 mm，出水管管體變形量增加1 mm，且出水管的管體變形量更為明顯；進、出水溫差增大2℃，出水管截面流通面積減小1.3%，而對進水管的影響不大。