自由約束條件下能源樁的離散元研究

王雪松， 龔曉南

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

0 引言

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

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

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

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

1 模型基本理論

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

1.1 離散元基本理論

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

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

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

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

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

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

1.2 傳熱模型基本理論

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

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

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

1.2.2 離散元顆粒傳熱理論

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

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

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

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

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

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

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

圖1 顆粒間傳熱

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

式中，Q（m）為顆粒周圍熱源流入該顆粒的熱量；Qv為內(nèi)熱源熱量；m為質(zhì)量；cv為比熱容。

2 能源樁幾何模型建立

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

2.1 現(xiàn)場試驗概況

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

2.2 離散元模型

2.2.1 土體模型

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

表1 顆粒半徑分布 m

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

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

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

2.2.2 能源樁樁體模型

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

圖2 能源樁模型

2.3 模型建模步驟

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

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

圖3 能源樁建模過程

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

2.4 傳熱過程模擬

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

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

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

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

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

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

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

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

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

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

圖4 樁身溫度分布

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

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

圖5 樁周土體溫度分布

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

3.2 樁周土體位移

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

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

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

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

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

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

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

4 結(jié)語

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

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

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

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

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