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

姚遠(yuǎn)宏，張澤雨，劉永成

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

0 引言

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

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

1 Matlab PDE tool 求解方法

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

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

熱傳導(dǎo)方程的通用偏微分方程為：

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

在求解模擬區(qū)域溫度u 時(shí)，對(duì)比式（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)的變量。熱傳導(dǎo)問題中的第一類邊界條件與Dirichlet 邊界條件對(duì)應(yīng)；第二類和第三類邊界條件與Neumann 邊界條件對(duì)應(yīng)。

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

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

2 凍土模型槽熱傳導(dǎo)模擬

2.1 模型槽介紹

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

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

圖2 測(cè)溫點(diǎn)布置圖

2.2 基本假設(shè)

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

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

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

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

（4）冷量?jī)H沿垂直于凍結(jié)管的方向傳導(dǎo)。

2.3 模擬計(jì)算

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

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

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

圖3 120 h 時(shí)模型槽中土體溫度分布圖

2.4 模擬溫度與實(shí)測(cè)溫度的對(duì)比分析

將測(cè)溫點(diǎn)2、3 的實(shí)測(cè)溫度與模擬溫度進(jìn)行對(duì)比，其瞬態(tài)溫度變化圖見圖4。

圖4 測(cè)溫點(diǎn)2、3 的瞬態(tài)溫度變化圖

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

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

將測(cè)溫點(diǎn)4 的實(shí)測(cè)溫度與模擬溫度進(jìn)行對(duì)比，其每隔12 h 的瞬態(tài)溫度值見表1。

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

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

3 結(jié)語

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

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