基于Matlab PDE 工具箱的土體凍結溫度場模擬研究

姚遠宏，張澤雨，劉永成

（中國建筑科學研究院天津分院，天津市 300384）

0 引言

隨著城市的發展，地鐵建設在各地不斷興起，并在城市交通中扮演著日益重要的角色。地鐵中的聯絡通道是連接2 條隧道的逃生通道和排水通道，在施工時主要以凍結法來固結土體后再開挖，因此合理地模擬土體凍結溫度場和凍結壁的發展速度，在提高凍結效率、節約成本、降低能耗方面具有一定的現實意義和經濟意義。

Matlab PDE tool 是一款求解偏微分方程并且能使解圖形化的軟件，可以解決復雜幾何圖形和復雜邊界條件的非穩定熱傳導問題。模型槽土體凍結溫度場是三維的熱傳導問題，而Matlab PDE tool 只能定義在二維空間。徐梓斌等[1]利用Matlab PDE tool解決了一個短圓柱鋼錠置于加熱爐中的三維瞬態非穩定熱傳導問題；王志良等[2]通過坐標轉換將三維導熱問題轉化為二維，并得出4 h 時刻燃燒棒的溫度分布云圖和熱流密度云圖。本文給出運用Matlab PDE tool 進行模型槽土體凍結溫度場模擬的一般步驟并分析驗證模擬值的準確性。

1 Matlab PDE tool 求解方法

Matlab PDE tool 利用有限元的方法求解偏微分方程。應用Matlab PDE tool 可以解決下類拋物線方程：

式中：Ω 為平面有界區域；u 為定義在區域Ω 上的實（或復數）函數；d、c、a 和f 為常數或變量；t 為時間變量。

熱傳導方程的通用偏微分方程為：

式中：u 為預模擬物體內的溫度；k 為導熱系數；qv為熱源；ρ 為密度；cp為比熱容。

在求解模擬區域溫度u 時，對比式（1）與式（2）中的方程系數可以看出，d=cpρ，c=k，f-au=qv。

Matlab PDE tool 定義了2 類邊界條件：

（1）Dirichlet 邊界條件：hu=r。

（2）Neumann 邊界條件：n·(c▽u)+qu=g。

其中：n 為垂直于邊界的單位矢量；h、r、q、c 和g 為常量或與u 有關的變量。熱傳導問題中的第一類邊界條件與Dirichlet 邊界條件對應；第二類和第三類邊界條件與Neumann 邊界條件對應。

運用Matlab PDE tool 求解模型槽土體凍結導熱流程一般步驟可按圖1 進行。

圖1 Ma tla b P DE tool 求解模型槽土體凍結導熱流程圖

2 凍土模型槽熱傳導模擬

2.1 模型槽介紹

本次運用Matlab PDE tool 模擬模型槽中-2 ℃的土，經過120 h 后的溫度場。模型槽尺寸為2.0 m×1.0 m×1.2 m，放置于10 ℃的室內，其側面及底部邊界裝有厚50 mm 的保溫板（硬泡沫板），并在四周鋪設塑料膜防止跑水。試驗土體取自天津某地鐵項目中的原狀土，填充模型槽后補水飽和。凍結管直徑為0.045 m 且水平貫穿寬度方向，在水平方向距離模型箱壁0.8 m，在豎直方向距離模型箱壁0.6 m，兩凍結管水平間距0.4 m。冷凍系統由冷凍機、循環水泵、外接鹽水箱組成，采用鹽水氯化鈣水溶液作為低溫冷媒劑循環吸收周圍土體的熱量。凍結時保溫板溫度為6.6 ℃，凍結管與周圍土體凍結時接觸面溫度為-25 ℃。通過試驗測得土的密度ρ=1860 kg/m3；比熱容cp=1550 J/（kg·℃）；土體導熱系數k=1.26 W/（m·℃）；模型槽土體表面與外界空氣的對流換熱系數為10 W/（m2·℃）。

測溫線采用國標PT1000 鉑熱電阻測溫線，安捷倫采集儀實時采集溫度值。6 個測溫點在同一平面內且測溫點1、2、3、4 所在直線與5、6 所在直線垂直；測溫點2、3 距離凍結管外壁10 cm，測溫點1、4距離凍結管外壁20 cm，測溫點5 距離凍結管外壁20 cm，測溫點6 距離凍結管外壁30 cm。測溫點布置圖見圖2。

圖2 測溫點布置圖

2.2 基本假設

土體凍結溫度場是一個有相變、移動邊界、熱源以及復雜邊界的導熱問題[3-5]，過程比較復雜。為了便于模擬計算，在采用Matlab PDE tool 進行模擬時作如下假設：

（1）土體為均勻連續體，且各向同性。

（2）土體的熱物理學性質在所研究區域內均相同。

（3）凍結過程中相變、水分遷移和凍脹變形均不考慮。

（4）冷量僅沿垂直于凍結管的方向傳導。

2.3 模擬計算

根據測溫點的位置，以凍結管水平貫穿的中點做正截面，以兩凍結管水平間距的中點為原點，以測點1、2、3、4 所在直線為X 軸，測點5、6 所在直線為Y 軸，建立平面直角坐標系。

根據式（1）、式（2）中的方程系數對應關系得d=cpρ=2883 000 J/（m3·℃），c=k=1.26 W/（m·℃），a=0，f=0；模型槽側面及底部有保溫板，不考慮熱損失，在Neumann 邊界條件中設置為q=0，g=0；模型槽土體上部與周圍空氣存在熱交換，因此設為第三類邊界條件，即Neumann 邊界條件，由對應關系得q=10×3600=36000，g=10×36000=360000；凍結管與周圍土體接觸面的邊界條件為Dirichlet 邊界條件，由對應關系得r=-25，h=1。

求得120 h 時模型槽中土體溫度分布圖，見圖3。

圖3 120 h 時模型槽中土體溫度分布圖

2.4 模擬溫度與實測溫度的對比分析

將測溫點2、3 的實測溫度與模擬溫度進行對比，其瞬態溫度變化圖見圖4。

圖4 測溫點2、3 的瞬態溫度變化圖

由圖4 可知：位于2 根凍結管中間部位的2 號測溫點的溫度基本小于與凍結管距離相同的3 號測溫點的溫度；2 個測溫點的溫度在剛開始時均呈快速下降趨勢，這是因為凍結管周圍土體首先吸收冷量，使其自身溫度開始降低；隨著凍結過程的進行，降溫速度逐漸變慢，這是由于土體中水凍結成冰需要釋放大量熱量，導致土體溫度有所升高，所以土體的降溫速度開始變慢。在快速降溫時，測溫點的模擬溫度與實測溫度相差很大，模擬溫度一般都大于實測溫度；當降溫速度逐漸變慢時，測溫點的模擬溫度與實測溫度在接近，當溫度達到一定的低溫時，模擬溫度與實測溫度基本一致。

出現模擬溫度大于實測溫度的原因一般有以下幾方面：（1）測溫儀器誤差；（2）土的導熱系數隨凍結時間逐漸增大，而模擬時假設其不變；（3）模型槽保溫較差，與環境存在熱交換；（4）由于土中水結冰，導致土的比熱容隨溫度降低而減小[6-8]。

將測溫點4 的實測溫度與模擬溫度進行對比，其每隔12 h 的瞬態溫度值見表1。

表1 測溫點4 每隔12 h 的瞬態溫度值

由表1 可知，對于距凍結管較遠的測溫點4 而言，其模擬溫度與實測溫度之間的誤差較大；隨著凍結時間的增加，兩者誤差有所減小并最終趨于穩定。分析原因可能有：（1）模型槽保溫板材料溫度會隨著凍結時間的增加而降低，但在邊界假設時設定的溫度為恒溫6.6 ℃；（2）模擬假設時忽略水分遷移，而水分遷移對溫度變化影響較大；（3）距離凍結管較遠之處的制冷效果沒有達到預設凍結溫度。

3 結語

（1）通過合理的假設和坐標軸建立，可利用Matlab PDE tool 來模擬三維空間的模型槽土體凍結溫度場熱傳導問題，操作簡單并且結果可圖形化。凍結管附近的模擬溫度曲線與實測溫度曲線變化趨勢基本一致，能有效預測實際工程凍結溫度場的趨勢和所需時間。

（2）因模擬時不考慮土內相變潛熱，使得凍結過程中，實測溫度與模擬溫度的誤差過大。隨著凍結時間的增加，土內相變潛熱對溫度的影響降低，使得模擬溫度與實測溫度的差別減小。在模擬時，根據水的狀態進行分階段模擬，可以提高模擬的準確性。