王雪松， 龔曉南

（浙江大學(xué)建筑工程學(xué)院濱海中心，杭州 310012）

0 引言

氣候變化和能源危機(jī)是如今面臨的兩大巨大挑戰(zhàn)，當(dāng)前形勢(shì)下新能源以其具有可再生、無(wú)污染、形式種類(lèi)豐富等優(yōu)點(diǎn)深受各國(guó)青睞。

最早由日本的Morino［1］提出樁埋管換熱器，而能源樁是一種新型樁基礎(chǔ)形式，該樁基礎(chǔ)結(jié)合地埋管換熱器的塑料換熱管與建筑物的混凝土樁基以實(shí)現(xiàn)結(jié)構(gòu)物承載及換熱功能。在國(guó)外方面，1948年Kelvin［2，3］提出線(xiàn)熱源理論，該理論假定豎直埋管換熱器為一恒定熱流密度的無(wú)限長(zhǎng)線(xiàn)熱源。Lee等［4］根據(jù)半離散化網(wǎng)絡(luò)的研究方法，提出三維有限差分的數(shù)值模型。Yoon等［5］基于有限元模擬分析，對(duì)比研究了不同埋管形式的能量樁在3個(gè)月內(nèi)通過(guò)間歇運(yùn)行和連續(xù)運(yùn)行下的換熱性能。Bezyaii等［6］通過(guò)三維數(shù)值模擬研究了不同換熱管的能量樁的換熱性能，發(fā)現(xiàn)螺旋型埋管換熱性能優(yōu)于單U型和W型埋管。

在國(guó)內(nèi)方面，石磊［7］推導(dǎo)出了熱源模型的溫度分布解析解，并進(jìn)一步提出了線(xiàn)圈模型。桑宏偉［8］嘗試使用二維離散元模型對(duì)能源樁的傳熱特性進(jìn)行研究，采用MatDEM離散元軟件對(duì)能源樁建立模型并通過(guò)與現(xiàn)場(chǎng)工程監(jiān)測(cè)結(jié)果對(duì)比，驗(yàn)證了離散元法在能源樁換熱性能研究方面的有效性。

文中提出能源樁的離散元模型，并通過(guò)MatDEM軟件對(duì)能源樁的傳熱性能和承載能力特性進(jìn)行研究。

1 模型基本理論

離散元法由Cundall［9］首先提出，并將其主要用于計(jì)算分析非連續(xù)介質(zhì)的各種力學(xué)行為。離散元受限于實(shí)際工程顆粒數(shù)量需求巨大導(dǎo)致的計(jì)算量較大等問(wèn)題，在巖土工程實(shí)際應(yīng)用中相對(duì)較少。而MatDEM軟件采用全新的算法較好地解決了計(jì)算量問(wèn)題，目前已經(jīng)可實(shí)現(xiàn)應(yīng)力場(chǎng)-位移場(chǎng)-溫度場(chǎng)的多場(chǎng)耦合，給離散元的工程實(shí)際應(yīng)用帶來(lái)了可能。

1.1 離散元基本理論

在文中離散元模型中采用線(xiàn)彈性赫茲接觸模型如下：

式中，F(xiàn)為接觸力合力，F(xiàn)n；Fs分別為法向接觸力和切向接觸力。

顆粒間的法向彈簧可以表示法向力Fn和法向變形Xn的關(guān)系，見(jiàn)式（2）。

式中，Kn為法向剛度；Xn是法向相對(duì)位移；Xb為斷裂位移。

同理切向彈簧可以表示顆粒間的剪切力Fs和剪切變形Xs的關(guān)系，見(jiàn)式（3）。

式中，Ks為切向剛度；Xs為切向位移。

1.2 傳熱模型基本理論

1.2.1 傳熱學(xué)相關(guān)理論

文中中主要考慮顆粒間的熱傳導(dǎo)，熱傳導(dǎo)是介質(zhì)內(nèi)無(wú)宏觀運(yùn)動(dòng)時(shí)的一種傳熱方式。根據(jù)傅里葉定律可推導(dǎo)出傳熱方程：

式中，T為溫度場(chǎng)向量；ρ為質(zhì)量密度；c為單位體積內(nèi)的熱量；t為時(shí)間；λ為導(dǎo)熱系數(shù)。

1.2.2 離散元顆粒傳熱理論

根據(jù)有限差分法可以推導(dǎo)出離散元顆粒的傳熱方程［10］。

推導(dǎo)中做如下假設(shè)：

兩相鄰熱源間的溫度差與熱阻間的關(guān)系。

式中，ΔT為溫度改變；η為單位長(zhǎng)度的熱阻；L為熱管長(zhǎng)度(此處簡(jiǎn)化為兩顆粒中心距離)；n→i為熱管方向向量；q→i為兩顆粒間的熱通量。

在模型中，假設(shè)溫度的變化只影響顆粒體積，且顆粒根據(jù)溫度變化線(xiàn)性膨脹，如式（6）所示。

式中，ΔR為顆粒因溫度改變的半徑；α為熱膨脹系數(shù)，m/℃；R為顆粒初始半徑；ΔT為溫度改變量。

在離散元體系中將顆粒視為獨(dú)立的熱源和鄰近顆粒的熱管構(gòu)成，熱管作為熱量的傳導(dǎo)媒介，以顆粒是否接觸來(lái)判斷熱量是否傳遞，若顆粒間無(wú)接觸則不傳到熱量，若兩顆粒重疊則視為顆粒間的熱管進(jìn)行傳熱。如圖1所示，T0顆粒與T1、T2、T3顆粒分別接觸，顆粒形心之間形成3個(gè)熱管，熱管長(zhǎng)度為形心間的距離Li=|R0Ri|，熱流流經(jīng)的截面Ai視為接觸面積，ni為接觸面指向低溫物體的法向方向，T0顆粒每次通過(guò)3個(gè)熱管與周?chē)w粒實(shí)現(xiàn)熱傳導(dǎo)，如圖1所示。

圖1 顆粒間傳熱

根據(jù)以上的假設(shè)，推導(dǎo)出顆粒的傳熱方程：

式中，Q（m）為顆粒周?chē)鸁嵩戳魅朐擃w粒的熱量；Qv為內(nèi)熱源熱量；m為質(zhì)量；cv為比熱容。

2 能源樁幾何模型建立

文中根據(jù)以上離散元傳熱模型，并根據(jù)郭易木等［11］論文中PHC能源樁原型試驗(yàn)建立能源樁幾何模型并驗(yàn)證模型有效性。

2.1 現(xiàn)場(chǎng)試驗(yàn)概況

現(xiàn)場(chǎng)勘查顯示地下水位于1.0m深左右，于0～3.6m為硬殼層，土層分界線(xiàn)位于地下11.1m左右，分界線(xiàn)以上土層薄弱，其下土層強(qiáng)度較高，土層的熱導(dǎo)率在1.15～1.17W/（m·K）范圍內(nèi)?，F(xiàn)場(chǎng)試驗(yàn)中PHC能源樁樁長(zhǎng)為24m，樁外徑為500mm，內(nèi)徑為280mm，樁體材料為C80高強(qiáng)混凝土，換熱管為單U型埋管。樁體彈性模量4.57×104MPa，熱膨脹系數(shù)αc=11.57×10-6/℃。

2.2 離散元模型

2.2.1 土體模型

文中建立能源樁的三維離散元模型，模型盒長(zhǎng)、寬均為4m，高度為28m，為保證計(jì)算的精確度并降低模型計(jì)算量，文中采用從中心到模型盒邊緣粒徑擴(kuò)大的方式建立土體模型，顆粒半徑分布如表1所示。

表1 顆粒半徑分布 m

模型總顆粒數(shù)為571940個(gè)，其中活動(dòng)顆粒為499739個(gè)，邊界單元72201個(gè)，在模型中將邊界單元定義為隔溫光滑邊界，不計(jì)算熱傳導(dǎo)。

為通過(guò)與現(xiàn)場(chǎng)試驗(yàn)對(duì)比驗(yàn)證模型有效性，參照現(xiàn)場(chǎng)土體條件設(shè)置模型的力學(xué)參數(shù)，并經(jīng)過(guò)多次試算確定模型的微觀力學(xué)參數(shù)見(jiàn)表2。

表2 顆粒力學(xué)性質(zhì)

2.2.2 能源樁樁體模型

能源樁直徑D1=500mm，樁長(zhǎng)L=24m，為保證樁體光滑，樁體采用顆粒半徑為0.03m的顆粒以一定系數(shù)重疊堆積而成，并設(shè)置為clump團(tuán)簇，能源樁的U型傳熱管在離散元中以顆粒豎向堆疊為兩根直管并緊貼樁體設(shè)置為clump團(tuán)簇，以此進(jìn)行傳熱模擬，能源樁模型如圖2所示。

圖2 能源樁模型

2.3 模型建模步驟

（1） 初始模型堆積，先建立1.1m×1.1m×28m的模型盒B1，隨機(jī)生成顆粒半徑為0.03m的顆粒經(jīng)過(guò)重力沉積和壓縮。再建立2.1m×2.1m×28m的模型盒B2，隨機(jī)生成顆粒半徑為0.04m的顆粒經(jīng)過(guò)重力沉積和壓縮，再切割B1模型盒生成半徑0.5m的圓柱體結(jié)構(gòu)體插入B2模型盒并平衡模型。同理依次生成不同顆粒半徑、不同尺寸的模型盒，并不斷切割結(jié)構(gòu)體插入新的模型盒。

（2） 根據(jù)現(xiàn)場(chǎng)試驗(yàn)條件對(duì)地層進(jìn)行簡(jiǎn)單劃分，以便模型建立完成后賦予不同地層不同顆粒力學(xué)性質(zhì)，地層劃分如圖3（a）所示。

圖3 能源樁建模過(guò)程

（3） 建立能源樁及U型換熱管并插入堆積好的土體模型中，給樁體及土體賦予不用的力學(xué)性質(zhì)，先斷開(kāi)土體膠結(jié)，讓土體在重力作用下填充間隙，模型充分平衡后再膠結(jié)地層與能源樁并最終實(shí)現(xiàn)模型平衡，如圖3（b）所示。

2.4 傳熱過(guò)程模擬

離散元的傳熱過(guò)程是通過(guò)給不同顆粒設(shè)置不同的溫度，以形成顆粒間的溫度差來(lái)進(jìn)行熱傳導(dǎo)模擬。能源樁的傳熱過(guò)程通過(guò)以下幾步進(jìn)行模擬：

（1） 先賦予顆粒初始溫度，再給U型換熱管設(shè)置相應(yīng)的加熱溫度。

（2） 在每個(gè)時(shí)間步內(nèi)根據(jù)有限差分法思想計(jì)算每個(gè)顆粒不同時(shí)間的溫度值。

（3） 根據(jù)顆粒溫度的改變，通過(guò)考慮顆粒線(xiàn)性膨脹更新顆粒半徑。

（4） 根據(jù)每一個(gè)時(shí)間步的顆粒半徑的改變，對(duì)模型進(jìn)行堆積和平衡。

3 模擬結(jié)果和分析

3.1 樁體溫度與分布規(guī)律

模型試驗(yàn)為對(duì)比現(xiàn)場(chǎng)試驗(yàn)，通過(guò)設(shè)置U型換熱管與樁土顆粒的溫度差進(jìn)行168h的傳熱模擬，將U型換熱管顆粒設(shè)為恒溫40℃，依據(jù)現(xiàn)場(chǎng)試驗(yàn)中土體及樁體溫度對(duì)模型中樁土顆粒初始溫度進(jìn)行設(shè)置，溫度設(shè)置見(jiàn)表3。

表3 土體及樁體顆粒初始溫度

圖4給出了數(shù)值模擬試驗(yàn)和郭易木等現(xiàn)場(chǎng)試驗(yàn)0、24、72、168h的樁身受熱后的溫度對(duì)比圖。根據(jù)模擬結(jié)果可得，能源樁加熱初始階段能源樁內(nèi)土體溫度及能源樁樁體溫度上升較快，但樁內(nèi)土體溫度分布具有一定差異，與桑宏偉論文中觀測(cè)現(xiàn)象一致，其原因一方面是靠近U型換熱管顆粒的土體顆粒溫度上升較快，另一方面是傳熱樁內(nèi)的土體顆粒分布不均勻，無(wú)法像連續(xù)介質(zhì)一樣完全均勻傳熱。當(dāng)樁體溫度上升超過(guò)30℃后傳熱趨于平穩(wěn)，樁身溫度增長(zhǎng)減緩，樁內(nèi)土體溫度分布趨于穩(wěn)定，且分布較傳熱初始階段更為均勻。

圖4 樁身溫度分布

168h時(shí)樁身平均溫度增長(zhǎng)18.8℃，樁體最高溫度達(dá)到38℃，樁身溫度總體呈現(xiàn)兩端低中部高的形式，原因分析為一方面由于樁體初始溫度分布呈現(xiàn)兩端低中間高現(xiàn)象；另一方面由于樁兩端散熱條件良好。但隨著傳熱，頂部樁體溫度升溫較快，原因分析為傳熱模型未考慮頂部土體及樁體向上方空氣中散熱的影響，在后續(xù)的建模中可對(duì)此進(jìn)行修正，模擬結(jié)果的樁體溫度與分布規(guī)律與現(xiàn)場(chǎng)試驗(yàn)結(jié)果基本一致。

樁外土體每隔5m依次取靠近樁體0.1m的土體平均溫度，不同時(shí)刻溫度分布如圖5所示，可見(jiàn)土體顆粒溫度在0～48h樁外土體溫度上升較小，說(shuō)明初始階段能源樁傳熱對(duì)樁外土體影響較小。

圖5 樁周土體溫度分布

樁外土體在48～168h持續(xù)穩(wěn)定上升，在168h時(shí)樁外土體最高溫度達(dá)到26.6℃，分布規(guī)律整體呈現(xiàn)中間高兩端低，是由于樁端以下土體散熱空間更大，沒(méi)有周?chē)鸁崃炕ハ嘤绊?，因此樁底土體散熱較快，而模型頂部因該離散元模型未考慮土體及樁體與空氣交換熱量而導(dǎo)致頂部土體溫度上升較快。

3.2 樁周土體位移

能源樁模型中通過(guò)U型換熱管傳熱，使樁、土體溫度發(fā)生改變并根據(jù)式（7）產(chǎn)生相應(yīng)的顆粒熱膨脹，取樁體及土體溫度變化處即埋深16.3m水平面的土體顆粒分析樁周土體位移變化規(guī)律。

如圖6、圖7所示，樁體內(nèi)部回填土體因受樁體限制顆粒位移較小，而隨著熱量傳遞可見(jiàn)樁外土體產(chǎn)生相應(yīng)的位移。樁外土體沿X、Y方向均產(chǎn)生對(duì)稱(chēng)、相反的位移分布，越靠近樁體的土體相對(duì)位移越大。而模型中U型管緊貼樁體分布，靠近U型管的樁體、土體溫度上升較快，鄰近土體相對(duì)其他位置土體沿著樁體徑向產(chǎn)生更大的位移，即樁外土體沿X方向位移相對(duì)沿Y方向位移絕對(duì)值較大。

圖6 樁周土體沿X方向位移

圖7 樁周土體沿Y方向位移

3.3 樁身縱向應(yīng)力變化規(guī)律

樁身受熱產(chǎn)生由中間向兩端的線(xiàn)性熱膨脹，而樁周土體及樁底土體對(duì)樁身產(chǎn)生約束，能源樁以鉆孔放入土體顆粒中進(jìn)行平衡的形式，故可視為樁體除自身重力作用外無(wú)初始應(yīng)力，能源樁的縱向附加熱應(yīng)力如圖8所示。模擬結(jié)果為與郭易木等現(xiàn)場(chǎng)試驗(yàn)對(duì)照，取現(xiàn)場(chǎng)試驗(yàn)應(yīng)變片位置的顆粒進(jìn)行附加應(yīng)力計(jì)算，可見(jiàn)能源樁隨著加熱產(chǎn)生軸向壓應(yīng)力，傳熱結(jié)束后樁體應(yīng)力總體呈現(xiàn)兩端小、中間大的趨勢(shì)，但相對(duì)現(xiàn)場(chǎng)試驗(yàn)測(cè)得結(jié)果較小，原因分析為與現(xiàn)場(chǎng)試驗(yàn)相比模型顆粒孔隙度較大、土體給樁體約束較小。在幾個(gè)分析點(diǎn)中樁體的最大附加應(yīng)力達(dá)4.12MPa，位于16.3m處，最小附加應(yīng)力為樁頂部，與現(xiàn)場(chǎng)試驗(yàn)現(xiàn)象基本一致。

圖8 樁身應(yīng)力分布

4 結(jié)語(yǔ)

文中使用南京大學(xué)自主研發(fā)的MATDEM軟件建立能源樁三維離散元傳熱模型模擬PHC能源樁（U型換熱管）傳熱過(guò)程，并與郭易木等論文中PHC能源樁原型試驗(yàn)的模擬結(jié)果進(jìn)行驗(yàn)證對(duì)比，得出以下結(jié)論：

（1） 在保證計(jì)算精度的前提下，可利用粒徑優(yōu)化方法，從模型中心向邊緣逐層堆積不同粒徑離散元顆粒來(lái)降低離散元模擬的計(jì)算量，文中驗(yàn)證了該方法在能源樁的傳熱模型中的可行性。

（2） 模型中樁體溫度上升趨勢(shì)與郭易木等現(xiàn)場(chǎng)試驗(yàn)結(jié)果基本一致，樁體溫度呈現(xiàn)兩端低中間高的變化，但由于離散元模型中無(wú)法考慮土體及樁體向空氣中散熱，因此模擬試驗(yàn)中樁體溫度相對(duì)現(xiàn)場(chǎng)試驗(yàn)結(jié)果較高，在后續(xù)的離散元模型中可以考慮對(duì)模型改進(jìn)以使結(jié)果更準(zhǔn)確。

（3） 在能源樁傳熱過(guò)程中，能源樁與樁周土體隨著傳熱的進(jìn)行溫度發(fā)生改變，且樁周土體隨著能源樁樁體膨脹而沿著遠(yuǎn)離能源樁方向產(chǎn)生對(duì)稱(chēng)、相反的位移分布。

（4） 將離散元傳熱模型應(yīng)用到實(shí)際工程的研究中具有可行性，文中已驗(yàn)證U型換熱管能源樁的傳熱模擬與現(xiàn)場(chǎng)試驗(yàn)結(jié)果具有較高的重合度，并且根據(jù)文中的建模思想調(diào)整換熱管顆粒模型可以實(shí)現(xiàn)不同形狀換熱管的能源樁的傳熱模擬，如單U型、雙U型、W型、螺旋型換熱管的模擬均可實(shí)現(xiàn)。