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

王雪松， 龔曉南

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

0 引言

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

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

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

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

1 模型基本理論

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

1.1 離散元基本理論

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

式中，F為接觸力合力，Fn；Fs分別為法向接觸力和切向接觸力。

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

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

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

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

1.2 傳熱模型基本理論

1.2.1 傳熱學相關理論

文中中主要考慮顆粒間的熱傳導，熱傳導是介質內無宏觀運動時的一種傳熱方式。根據傅里葉定律可推導出傳熱方程：

式中，T為溫度場向量；ρ為質量密度；c為單位體積內的熱量；t為時間；λ為導熱系數。

1.2.2 離散元顆粒傳熱理論

根據有限差分法可以推導出離散元顆粒的傳熱方程［10］。

推導中做如下假設：

兩相鄰熱源間的溫度差與熱阻間的關系。

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

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

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

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

圖1 顆粒間傳熱

根據以上的假設，推導出顆粒的傳熱方程：

式中，Q（m）為顆粒周圍熱源流入該顆粒的熱量；Qv為內熱源熱量；m為質量；cv為比熱容。

2 能源樁幾何模型建立

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

2.1 現場試驗概況

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

2.2 離散元模型

2.2.1 土體模型

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

表1 顆粒半徑分布 m

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

為通過與現場試驗對比驗證模型有效性，參照現場土體條件設置模型的力學參數，并經過多次試算確定模型的微觀力學參數見表2。

表2 顆粒力學性質

2.2.2 能源樁樁體模型

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

圖2 能源樁模型

2.3 模型建模步驟

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

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

圖3 能源樁建模過程

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

2.4 傳熱過程模擬

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

（1） 先賦予顆粒初始溫度，再給U型換熱管設置相應的加熱溫度。

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

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

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

3 模擬結果和分析

3.1 樁體溫度與分布規律

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

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

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

圖4 樁身溫度分布

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

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

圖5 樁周土體溫度分布

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

3.2 樁周土體位移

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

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

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

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

3.3 樁身縱向應力變化規律

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

圖8 樁身應力分布

4 結語

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

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

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

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

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