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

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

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

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

1 地下?lián)Q熱器管土結(jié)構(gòu)凍脹變形過(guò)程數(shù)值模型的建立

1.1 控制方程

土壤凍脹擠壓埋管變形的過(guò)程主要用土壤孔隙增長(zhǎng)率函數(shù)、凍土應(yīng)力-應(yīng)變關(guān)系以及相變傳熱方程進(jìn)行描述。

(1) 土壤孔隙增長(zhǎng)率函數(shù)：凍土體積膨脹速率可用孔隙增長(zhǎng)率函數(shù)進(jìn)行描述，文獻(xiàn)[18]已對(duì)其進(jìn)行驗(yàn)證?？紫对鲩L(zhǎng)率N可表達(dá)為：

(1)

(2) 凍土應(yīng)力-應(yīng)變關(guān)系：可把凍土變形視作彈性變形，彈性應(yīng)變?cè)隽亢蛢鼋Y(jié)過(guò)程中孔隙增長(zhǎng)引起的應(yīng)變?cè)隽繕?gòu)成了總應(yīng)變?cè)隽俊雒浺鸬目倯?yīng)變?cè)隽恐械母鱾€(gè)分量可表達(dá)為：

(2)

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

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

(3)

1.2 計(jì)算流程

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

圖1 基于ABAQUS軟件的地下?lián)Q熱器管土結(jié)構(gòu)凍脹變形過(guò)程的計(jì)算流程Fig.1 Calculation flow of frost heave process of pipe soil structure of underground heat exchanger based on ABAQUS software

1.3 數(shù)值模型與參數(shù)

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

圖2 單U型地下?lián)Q熱器管土結(jié)構(gòu)數(shù)值模型Fig.2 Model of pipe soil structure of the single U-shaped underground heat exchanger

地下?lián)Q熱器管土結(jié)構(gòu)模型中換熱管采用高密度聚乙烯材料，凍結(jié)土壤由水、冰和土顆粒組成，其基本參數(shù)設(shè)定見(jiàn)表1。

表1 地下?lián)Q熱器管土結(jié)構(gòu)數(shù)值模型的基本參數(shù)

已針對(duì)該數(shù)值模型在土體凍結(jié)溫度場(chǎng)和換熱管變形應(yīng)變兩方面進(jìn)行了試驗(yàn)驗(yàn)證，詳見(jiàn)參考文獻(xiàn)[19]。

2 結(jié)果與分析

2.1 不同恒溫工況對(duì)地下?lián)Q熱器管土結(jié)構(gòu)凍脹變形的影響

針對(duì)恒溫工況，本文模擬了不同低溫運(yùn)行模式對(duì)管土結(jié)構(gòu)凍脹變形的影響，模擬設(shè)置100 h的運(yùn)行時(shí)間，對(duì)進(jìn)、出水管內(nèi)壁分別施加三種運(yùn)行模式HW1、HW2和HW3，其運(yùn)行溫度見(jiàn)表2。

表2 地下?lián)Q熱管恒溫運(yùn)行模式

2.1.1 土體凍結(jié)區(qū)變形

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

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

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

圖4為不同恒溫運(yùn)行模式下管間區(qū)域土體特征點(diǎn)O處的X向應(yīng)力隨土體平均凍結(jié)直徑的dxy變化曲線。

圖4 不同恒溫運(yùn)行模式下管間區(qū)域土體特征點(diǎn)O處的 X向應(yīng)力隨土體平均凍結(jié)直徑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可見(jiàn)，隨凍結(jié)直徑發(fā)展，O點(diǎn)的壓應(yīng)力不斷增大，壓應(yīng)力在三種模式間呈現(xiàn)出明顯且逐漸增大的差異，表現(xiàn)為HW3

2.1.2 管體變形

模擬過(guò)程中換熱管的直徑變化量見(jiàn)圖5。

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

圖5 不同恒溫運(yùn)行模式下?lián)Q熱管直徑的變化量隨 運(yùn)行時(shí)間的發(fā)展曲線Fig.5 Curves of pipe diameter variation with operation time under different operation modes of constant temperature

2.1.3 管截面變化

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

(4)

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

當(dāng)Cs>0時(shí)，表明管體變形使管截面的流通面積減?。划?dāng)Cs<0時(shí)，則表明管體變形使管截面的流通面積增大。

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

圖6 不同恒溫運(yùn)行模式下管截面當(dāng)量流通面積變化 率(Cs)隨土體平均凍結(jié)直徑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可見(jiàn)，可以發(fā)現(xiàn)Cs值均大于0，因此Cs值可視為流通截面減小率，隨著凍結(jié)區(qū)發(fā)展，該值基本呈線性增大趨勢(shì)，管體變形使其流通截面不斷減小。盡管進(jìn)水管的管徑變化量略小于出水管，但前者的Cs值卻稍大于后者，因此即使管徑變化量稍大，其流通截面減小率不一定大。

隨著凍結(jié)區(qū)域發(fā)展，Cs值的增大幅度依次為HW1>HW2>HW3。當(dāng)凍結(jié)直徑增至150 mm時(shí)，HW1和HW2運(yùn)行模式的Cs值小于HW3運(yùn)行模式，而當(dāng)凍結(jié)直徑達(dá)到200 mm時(shí)，HW2和HW3運(yùn)行模式的Cs值開始小于HW1運(yùn)行模式，根據(jù)發(fā)展趨勢(shì)可以推斷，HW3運(yùn)行模式的Cs值逐漸會(huì)小于HW2運(yùn)行模式。因此對(duì)于低溫運(yùn)行工況，當(dāng)凍結(jié)范圍較小時(shí)，低溫運(yùn)行模式的管體流通截面減小程度要大于高溫運(yùn)行模式，而當(dāng)凍結(jié)范圍超過(guò)一定值時(shí)，高溫運(yùn)行模式減小程度將超過(guò)低溫運(yùn)行模式。

2.2 不同進(jìn)、出水溫差對(duì)地下?lián)Q熱器管土結(jié)構(gòu)凍脹變形的影響

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

表3 進(jìn)、出水溫差運(yùn)行模式

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

2.2.1 土體凍結(jié)區(qū)變形

不同進(jìn)、出水溫差(DT1、DT2和DT3)運(yùn)行模式下隨運(yùn)行時(shí)間的發(fā)展土體平均凍結(jié)直徑dxy的變化曲線見(jiàn)圖7。

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

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

本文以管間區(qū)域土體特征點(diǎn)O處的應(yīng)力變化來(lái)衡量?jī)雒浟Φ淖兓匦?，得到不同進(jìn)、出水溫差(DT1、DT2和DT3)運(yùn)行模式下，土體特征點(diǎn)O處X向應(yīng)力隨運(yùn)行時(shí)間的變化曲線，見(jiàn)圖8。

圖8 不同進(jìn)、出水溫差(DT1、DT2和DT3)運(yùn)行模式 下管間區(qū)域土體特征點(diǎn)O處X向應(yīng)力隨運(yùn)行 時(shí)間的發(fā)展曲線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可見(jiàn)，三種不同溫差模式下O點(diǎn)位置的壓應(yīng)力有所差異?？梢园l(fā)現(xiàn)， DT1運(yùn)行模式的凍脹力最小，DT2運(yùn)行模式略大，DT3運(yùn)行模式最大，運(yùn)行最終，該位置的壓應(yīng)力分別達(dá)2.31 MPa、2.37 MPa和2.44 MPa，溫差增大2℃，可使壓應(yīng)力增大0.13 MPa。由此可見(jiàn)，較大溫差的運(yùn)行模式可在凍結(jié)區(qū)產(chǎn)生較大的凍脹力。

2.2.2 管體變形

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

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

由圖9可見(jiàn)，相比之下，出水管的變形量更為明顯；對(duì)于進(jìn)水管，管徑變化表現(xiàn)為DT3

2.2.3 管截面變化

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

圖10 不同進(jìn)、出水溫差(DT1、DT2和DT3)運(yùn)行模式下 管截面當(dāng)量流通面積變化率(Cs)隨運(yùn)行時(shí)間的 發(fā)展曲線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可見(jiàn)，Cs值均大于0且隨運(yùn)行時(shí)間不斷增大，變形使得管體流通截面不斷減小。對(duì)于進(jìn)水管[見(jiàn)圖10(a)]，三種運(yùn)行模式的Cs值在運(yùn)行初期有所不同，表現(xiàn)為DT1

3 結(jié) 論

本文通過(guò)建立地下?lián)Q熱器管土結(jié)構(gòu)凍脹變形過(guò)程數(shù)值模型，以0℃以下運(yùn)行的單U型地下?lián)Q熱器管土結(jié)構(gòu)為研究對(duì)象，模擬分析了不同恒溫溫度和不同進(jìn)、出水溫差對(duì)管土結(jié)構(gòu)凍脹變形的影響，得到以下結(jié)論：

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

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