水-熱-力耦合的單管凍結(jié)理論及其數(shù)值模擬

劉紀(jì)峰, 張會芝,余 照, 付佳寧

(1.三明學(xué)院 建筑工程學(xué)院， 福建 三明365004;2.工程材料與結(jié)構(gòu)加固福建省高等學(xué)校重點(diǎn)實(shí)驗(yàn)室(三明學(xué)院)， 福建 三明365004)

0 引言

城市地鐵旁通道施工常處于復(fù)雜的工程地質(zhì)和水文地質(zhì)環(huán)境，控制不好，很容易造成工程事故[1-3].人工凍結(jié)法加固地層具有抗?jié)B性好、強(qiáng)度高、不受深度限制、環(huán)境友好等優(yōu)點(diǎn)，在地下工程施工中得到廣泛應(yīng)用，但其凍脹融沉需加以預(yù)控[4].目前，人工凍結(jié)法的主要研究方法包括試驗(yàn)研究[5,6]、理論分析[7,8]、工程實(shí)測[9,10]和數(shù)值模擬[12-14]等.以上研究方法各有其優(yōu)點(diǎn)和不足.考慮到問題的復(fù)雜性，對人工凍結(jié)的相關(guān)研究要綜合選擇各方法的優(yōu)點(diǎn)，避免單一手段的不足.在以往研究的基礎(chǔ)上，推導(dǎo)考慮水-熱-力耦合的土體凍結(jié)理論公式并基于ANASYS二次開發(fā)進(jìn)行其單管凍結(jié)的數(shù)值模擬，通過對比分析驗(yàn)證該理論的有效性.

1 水-熱-力耦合的軟土凍結(jié)理論

1.1 模型假設(shè)

模型假設(shè)如下[15-17]：(1)水分遷移只以液態(tài)水形式進(jìn)行，忽略氣相和液相之間轉(zhuǎn)變時水分的蒸發(fā)轉(zhuǎn)移；(2)土的凍結(jié)主要是由熱傳導(dǎo)控制；(3)忽略荷載作用下土中冰點(diǎn)的降低；(4)土體固結(jié)已完成，忽略凍結(jié)過程中未凍區(qū)的固結(jié)；(5)凍結(jié)過程中土顆粒體積不變；(6)土體為各相同性體.

1.2 經(jīng)驗(yàn)公式選取和公式推導(dǎo)

忽略單元體凍結(jié)過程中熱對流的作用，傳導(dǎo)產(chǎn)生的凈熱流量，等于水冷卻(熱容量C)時的熱量與水轉(zhuǎn)化為冰時的相變潛熱之和(潛熱L).

熱傳導(dǎo)溫度場方程選用考慮潛熱的Harlan[19]方程：

(1)

熱量與溫度關(guān)系用Fourier定律描述, 其含義為單位時間通過單位面積的熱流量和溫度梯度成正比，即：

(2)

其中：q為熱流密度，J/(m2·s)；λ為導(dǎo)熱系數(shù)，W/(m·℃)；T為溫度，℃；C為體積熱容量，J/(m3·℃);t為時間，s；L為潛熱，J/kg；ρi為冰的密度，kg/m3；θi為冰體積分?jǐn)?shù).負(fù)號表示溫度梯度的方向和熱流方向相反.

將式(2)代入(1)并擴(kuò)展到三維場得到 ：

(3)

在土的凍結(jié)過程中，單元體中由于壓力梯度的建立和冰的形成使得從孔隙中流出的液態(tài)水被保存下來，考慮物質(zhì)遷移平衡過程，可以得到：

(4)

其中：ρl為未凍水(液態(tài))密度，kg/m3；θl為未凍水體積分?jǐn)?shù)；Qx為流量，m3/s.

方程左邊代表通過單元體某方向兩個面的凈物質(zhì)流量，右邊所指的是未凍水及冰形成產(chǎn)生的物質(zhì)流量.

引入描述物質(zhì)擴(kuò)散過程的Fick定律[15]:

(5)

其中：Pl為未凍水壓力，Pa；k為對流傳質(zhì)系數(shù)，m2/(s·Pa).

將式(5)代入式(4)并擴(kuò)展成三維形式，可得穩(wěn)態(tài)和非穩(wěn)態(tài)下的物質(zhì)遷移方程為：

(6)

壓力、密度與溫度關(guān)系引用Clapeyron方程描述：

(7)

其中：Pi為冰壓力，Pa；Pl為液體壓力，Pa；Tk為開爾文溫度(K=℃+273.16)；T0為參考溫度，273.16 K.

假設(shè)在0 ℃冰未產(chǎn)生，所有的孔隙都被未凍水占據(jù)，未凍水量與溫度關(guān)系的半經(jīng)驗(yàn)公式選用線性關(guān)系[16]：

(8)

式中n0為凍土的初始孔隙率.

基于Shen[19]關(guān)于冰壓力分布的建議，假定冰壓力在凍結(jié)前緣是線性分布的，對于二維有：

(9)

式中k為對流傳質(zhì)系數(shù)，m2/(s·Pa).

式(3)和式(6)～(9)有T，θl，θi，Pl，Pi五個未知量，一旦這些量求出，則其中的冰壓力Pi及液體壓力Pl可被作為應(yīng)力分析有限元中的初始應(yīng)力，得出凍脹應(yīng)力為：

σ0=(Pi+Pl)δ,

(10)

式中δ為克羅克內(nèi)爾符號.

凍脹力可表示為：

(11)

式中BT為體積函數(shù).

經(jīng)推導(dǎo)，可得式(12)所示的總體控制微分方程

(12)

假設(shè)一個單元中的溫度變化極小，忽略高階無窮小項(xiàng)，并將式(9)代入式(12)，方程可最終簡化為：

(13)

2 有限元離散化

2.1 整體離散

應(yīng)用Galerkin法對方程(13)進(jìn)行有限元離散化計(jì)算.經(jīng)推導(dǎo)，可得：

(14)

將其寫成矩陣形式：

(15)

式中

對于每一個單元，有

(16)

上述方程為非線性方程，可以采用牛頓迭代法求解.

2.2 流體壓力場

一旦溫度場T(x,y,t)在一個給定的時間t處建立后，可以計(jì)算出流體壓力場Pl和冰壓力場Pi，將式(3)寫成二維形式得到：

(17)

將式(17)代入式(6)并進(jìn)行整理可得：

(18)

對式(18)應(yīng)用Galerkin’s方法得到有限元方程：

(19)

式(19)可寫成式(20)所示的離散形式

ΔtQ/Pl=-(ΔtU/T+R/ΔT) ,

(20)

其中

(21)

可以采用牛頓迭代法求解，流體壓力Pl可以由式(20)計(jì)算出來.

2.3 冰壓力場

研究表明，Claperon方程偏高估計(jì)了冰壓力值，例如，當(dāng)溫度為-5 ℃時計(jì)算得出冰壓力Pi=5.5 MPa，為了得到更合理的值，采用文獻(xiàn)[20]的修正Claperon方程，即Pl=1.091Pi或者Pi=0.9166Pi.

由式(10)及式(11)可以得出由冰/流體壓力導(dǎo)致的荷載，如果將此荷載施加至有限元模型，則可以得到由于熱傳導(dǎo)及水、氣遷移引起的相應(yīng)的位移和應(yīng)力場.

3 數(shù)值實(shí)例驗(yàn)證

研究凍結(jié)壁溫度變化特征和形成規(guī)律的數(shù)值解法的主要優(yōu)點(diǎn)是其能用于比較復(fù)雜的情形，能求解分析解法所不能解決的大量導(dǎo)熱問題.隨著計(jì)算機(jī)技術(shù)的發(fā)展和應(yīng)用，數(shù)值解法解題的范圍、規(guī)模、速度和精確度有了很大的提高.為驗(yàn)證本文模型的實(shí)際效果，在ANSYS中建立單管凍結(jié)的數(shù)值模型.

3.1 參數(shù)選取與建模

模型為一圓形凍結(jié)單管，直徑89 mm，管外部土體半徑為2000 mm，凍結(jié)管在土體圓心位置，土體初始地溫為18 ℃，開始凍結(jié)后，從凍結(jié)管中產(chǎn)生源源不斷的冷量，凍結(jié)管的放熱系數(shù)即熱流為-267 W/m2，凍結(jié)管壁溫度經(jīng)鹽水冷凍至-28 ℃，持續(xù)時間為3 630 000 s(約42 d).材料參數(shù)還包括以下內(nèi)容：

(1) 對于凍土及未凍土可以賦予不同的材料屬性值，如密度、比熱、導(dǎo)熱系數(shù).

(2) 對于土體結(jié)冰釋放的潛熱可以在軟件中通過賦予材料的焓值來實(shí)現(xiàn).

參數(shù)具體取值如表1.其中土體熱交換系數(shù)50 W/(m2·℃)，空氣熱交換系數(shù)10 W/(m2·℃).

表1 各參數(shù)取值表Tab. 1 The parameter values

3.2 ANSYS建模與計(jì)算結(jié)果

ANSYS程序中土體采用PLANE13熱力耦合直接法分析單元.用自帶的APDL語言將所推導(dǎo)的三場耦合理論通過用戶自定義模塊輸入進(jìn)去進(jìn)行計(jì)算.整個凍結(jié)過程中，一些參數(shù)隨溫度、壓力、時間改變，因此要在ANSYS中建立表Table輸入.

3.3 數(shù)值模擬結(jié)果分析

將ANSYS結(jié)果(如圖1-圖6)與實(shí)際工程對比可知：

(1)凍結(jié)溫度場較為吻合，最終凍結(jié)壁厚度約1.5 m，與同條件實(shí)際工程效果較為一致，ANSYS計(jì)算結(jié)果不均勻，溫度梯度跨度較大，計(jì)算精度有待于進(jìn)一步提高.

(2)凍結(jié)壁徑向應(yīng)力受拉最大為0.39 MPa，發(fā)生在凍結(jié)管周邊，受壓最大為0.2 MPa，發(fā)生在凍結(jié)壁周邊.而同條件下實(shí)際工程凍結(jié)土體監(jiān)測結(jié)果為受壓最大為0.4 MPa，基本吻合.

(3)凍結(jié)壁剪應(yīng)力最大為0.36 MPa，最小為-0.27 MPa，分布不規(guī)律，實(shí)際工程中難以監(jiān)測剪應(yīng)力數(shù)值，故無法對比，但其值均在凍土抗剪強(qiáng)度指標(biāo)(1.0 MPa左右)之內(nèi)，未發(fā)生破壞.

(4)脹位移最大為4.15 cm，比實(shí)際工程同條件下聯(lián)絡(luò)通道凍結(jié)土體監(jiān)測結(jié)果大，這是因?yàn)閷?shí)際工程是多管凍結(jié)，且有外在約束，荷載較多，凍脹量應(yīng)比簡化條件下的單管理論計(jì)算小，本模型邊界條件有待進(jìn)一步優(yōu)化.

圖1 網(wǎng)格劃分及邊界條件Fig. 1 Mesh generation and boundary conditions

圖2 溫度場Fig. 2 Temperature field

圖3 Z向(徑向)應(yīng)力Fig. 3 Z-stress

圖4 第三主應(yīng)力Fig. 4 The third principal stress

圖5 剪應(yīng)力Fig. 5 Shear stress

圖6 凍脹位移Fig. 6 Frost heave displacement

4 結(jié)論與展望

推導(dǎo)了考慮水-熱-力耦合的單管凍結(jié)理論計(jì)算公式，用ANSYS數(shù)值軟件對某算例進(jìn)行計(jì)算并與實(shí)際工程進(jìn)行對比分析.結(jié)果表明：計(jì)算結(jié)果與實(shí)際工程監(jiān)測結(jié)果基本吻合，偏差都在解析解及實(shí)際工程變化范圍內(nèi)，考慮到土體凍結(jié)工程本身也是一種不確定問題，難以取得準(zhǔn)確解，因此，推導(dǎo)的理論公式計(jì)算結(jié)果是可行的，但也有以下問題需要改進(jìn)，比如Clapeyron方程作為非飽和正凍土在非穩(wěn)定場作用下的水分遷移的理論依據(jù)造成的計(jì)算誤差問題、土體各種參數(shù)隨環(huán)境條件變化問題、理論公式數(shù)值模擬過程中的優(yōu)化算法問題等，這也是下一步研究努力的方向.