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

姚遠宏，張澤雨，劉永成

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

0 引言

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

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

1 Matlab PDE tool 求解方法

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

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

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

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

在求解模擬區(qū)域溫度u 時，對比式（1）與式（2）中的方程系數(shù)可以看出，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 有關(guān)的變量。熱傳導問題中的第一類邊界條件與Dirichlet 邊界條件對應；第二類和第三類邊界條件與Neumann 邊界條件對應。

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

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

2 凍土模型槽熱傳導模擬

2.1 模型槽介紹

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

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

圖2 測溫點布置圖

2.2 基本假設

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

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

（2）土體的熱物理學性質(zhì)在所研究區(qū)域內(nèi)均相同。

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

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

2.3 模擬計算

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

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

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

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

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

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

圖4 測溫點2、3 的瞬態(tài)溫度變化圖

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

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

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

表1 測溫點4 每隔12 h 的瞬態(tài)溫度值

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

3 結(jié)語

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

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