人工凍結(jié)砂土熱力耦合模型分析

高亞飛，曹虎生，麻世壟

(1.陜西省一八五煤田地質(zhì)有限公司，陜西 榆林 719000;

2.安徽理工大學(xué)土木建筑學(xué)院，安徽 淮南 232001)

0 引 言

在工程施工中，土體的凍脹、融沉?xí)o道路和建筑物造成很大的危害，有必要對(duì)凍土施工中土體的凍脹融沉進(jìn)行研究。何敏等[1]以飽和凍土多孔介質(zhì)傳熱理論為基礎(chǔ)，運(yùn)用克拉伯龍方程揭示溫度梯度對(duì)凍結(jié)區(qū)水分遷移的影響規(guī)律，簡(jiǎn)化數(shù)學(xué)計(jì)算模型，建立了實(shí)用性的水熱力耦合數(shù)學(xué)模型；白青波等[2]基于熱傳導(dǎo)理論和達(dá)西定律理論，通過(guò)Comsol軟件二次開發(fā)，運(yùn)用偏微分方程對(duì)水熱耦合求解，實(shí)現(xiàn)凍土水分場(chǎng)和溫度場(chǎng)全耦合數(shù)值計(jì)算；劉建軍等[3]建立多年凍土區(qū)的水熱力耦合模型，考慮了水分遷移和冰水相變的影響，運(yùn)用Comsol數(shù)值計(jì)算軟件對(duì)多年凍土進(jìn)行水熱力耦合計(jì)算，得出了多年凍土的位移曲線；韓建澤[4]基于海相軟土降溫凍結(jié)過(guò)程中出現(xiàn)的二次相變現(xiàn)象，對(duì)溫度場(chǎng)、水分場(chǎng)、應(yīng)力場(chǎng)和鹽分場(chǎng)控制方程進(jìn)行修正，建立了水-鹽-熱-力四場(chǎng)耦合的人工凍結(jié)海軟土凍脹預(yù)報(bào)模型；王永濤等[5]基于飽和青藏粉土在不同溫度條件下的單向凍結(jié)試驗(yàn)，對(duì)凍結(jié)過(guò)程中凍脹速率和凍結(jié)厚度的影響因素進(jìn)行研究分析，得出影響土體凍脹率變化的內(nèi)在機(jī)理為土樣凍結(jié)過(guò)程中由冰透鏡體分凝冰所導(dǎo)致的未凍區(qū)固結(jié)和已凍區(qū)凍脹共同作用的結(jié)果；李銳陽(yáng)[6]通過(guò)凍土理論和水熱遷移的試驗(yàn)?zāi)M及雙圈管凍結(jié)模型分析，對(duì)深厚地層土壤水熱遷移進(jìn)行系統(tǒng)分析；高建強(qiáng)等[7]在室內(nèi)單向凍結(jié)試驗(yàn)的基礎(chǔ)上，分析了在考慮覆蓋層和補(bǔ)水條件下，溫度、土質(zhì)和溫度邊界條件對(duì)非飽和粗顆粒填料凍脹特性的影響。

大多數(shù)學(xué)者認(rèn)為，引起土體的凍脹融沉主要來(lái)源于凍結(jié)土層內(nèi)部水分遷移以及因水分遷移而導(dǎo)致的水分重分布。然而，水分重分布主要為冰水相變，導(dǎo)致冰水相變的主要因素為土體內(nèi)部的凍結(jié)溫度。本文基于前人研究的基礎(chǔ)，研究人工凍結(jié)砂土的熱力耦合機(jī)理，通過(guò)Comsol Multiphysics數(shù)值計(jì)算軟件中多孔介質(zhì)傳熱模塊、固體力學(xué)模塊，建立多物理場(chǎng)耦合模型，模擬不同凍結(jié)溫度下的飽和砂土內(nèi)部?jī)鼋Y(jié)溫度場(chǎng)的分布及凍脹量大小，對(duì)采用人工凍結(jié)法施工的工程具有一定的參考意義。

1 凍結(jié)砂土熱力耦合分析

1.1 熱力耦合模型基本假定

假設(shè)砂土在未凍結(jié)以及凍結(jié)狀態(tài)下均為連續(xù)、各向同性的彈性材料；凍結(jié)砂土內(nèi)部被液態(tài)水和固態(tài)冰填滿。

冰水相變的熱量的改變僅以熱傳導(dǎo)的方式進(jìn)行，忽略熱輻射和熱對(duì)流；試驗(yàn)為不補(bǔ)水試驗(yàn)，冰水相變只考慮原含水率；溫度載荷以實(shí)際試驗(yàn)循環(huán)的凍結(jié)溫度直接加在模型頂部邊界；數(shù)值計(jì)算模型兩側(cè)面邊界為對(duì)流換熱邊界條件，底部為絕熱邊界條件；在力學(xué)場(chǎng)中，模型頂面為自由邊界，兩側(cè)面為輥支撐約束，底部為固定約束。

1.2 溫度場(chǎng)控制方程

在未凍結(jié)區(qū)，二維溫度場(chǎng)的控制微分方程為

(1)

在凍結(jié)區(qū)，二維溫度場(chǎng)的控制微分方程為

(2)

凍土中未凍水相變釋放的相變潛熱為

(3)

式中，Tu和Tf分別為未凍結(jié)區(qū)與凍結(jié)區(qū)的溫度；ku和kf分別為未凍結(jié)區(qū)與凍結(jié)區(qū)的導(dǎo)熱系數(shù)；Cu和Cf分別表為未凍土與凍結(jié)土的比熱容；q0為源項(xiàng)；Q為相變潛熱；θi為體積含冰率；t為時(shí)間；pi為冰的密度；L為冰水相變潛熱，為335 kJ/kg。

初始條件為

T|t=0=T0

(4)

式中，T0為土體的初始溫度。

模型中，直接以試驗(yàn)循環(huán)的冷媒劑的溫度加在試樣的頂部邊界，因此頂部邊界條件為

T|y=yc=Tc(t)

(5)

式中，yc為模型上邊界的縱坐標(biāo)；Tc(t)為循環(huán)的冷媒劑的溫度。

模型中，兩側(cè)面考慮對(duì)流傳熱的影響，有對(duì)流換熱邊界條件，即

q0=h(t)×(Text-T)

(6)

式中，h(t)為對(duì)流換熱系數(shù)；Text為外部溫度。

以上式(1)～(6)構(gòu)成二維凍結(jié)溫度場(chǎng)的定解問(wèn)題。

1.3 熱力耦合控制方程

熱力耦合情況下的力學(xué)場(chǎng)為平面應(yīng)變問(wèn)題，其微分控制方程為

(7)

(8)

(9)

式中，σx為平面x方向正應(yīng)力；εx為平面x方向正應(yīng)變；σy為平面y方向正應(yīng)力；εy為平面y方向正應(yīng)變；τxy為平面剪應(yīng)力；εxy為平面剪應(yīng)變；T為溫度差；E(T)為考慮溫度變化的土體彈性模量；u為泊松比；α為熱膨脹系數(shù)。

2 數(shù)值模擬

試驗(yàn)所用砂為榆林某路基砂土，試驗(yàn)試樣為直徑50 mm、高度100 mm的圓柱體。為簡(jiǎn)化計(jì)算，數(shù)值計(jì)算模型采用二維平面模型，所建試樣模型尺寸為寬50 mm、高100 mm的矩形，網(wǎng)格劃分采用自由剖分三角形網(wǎng)格，上下邊界為試驗(yàn)邊界條件，左右邊界考慮相應(yīng)的熱散失。幾何模型由1個(gè)域、4個(gè)邊界、4個(gè)頂點(diǎn)組成。以試樣低端為0 ℃面，在20、40、60、80 mm處埋設(shè)溫度探針。網(wǎng)格劃分見圖1。

圖1 網(wǎng)格劃分(單位：mm)

數(shù)值模擬計(jì)算時(shí)，試樣上端為冷端，試樣下端為暖端。通過(guò)數(shù)值模擬計(jì)算，得到試樣在不同冷端溫度下試樣溫度達(dá)到穩(wěn)定時(shí)的溫度云圖，見圖2。從圖2可知，初始含水率一定時(shí)，在上部不同冷端溫度下，隨著冷端溫度的降低，試樣沿高度方向溫度區(qū)域劃分更明顯。因凍結(jié)端面向下傳遞的冷量受土體深度影響，當(dāng)土體深度一定時(shí)，隨著冷端溫度的降低，熱阻造成的冷量損失減小，故土層溫度穩(wěn)定時(shí)整體溫度偏低。

圖2 不同冷端溫度下試樣溫度云圖(單位：℃)

在數(shù)值模擬計(jì)算過(guò)程中，通過(guò)時(shí)間計(jì)算步長(zhǎng)可知，當(dāng)冷端溫度為-5 ℃時(shí)，試樣低端溫度達(dá)到0 ℃所需時(shí)間為128 h；當(dāng)冷端溫度為-10 ℃時(shí)，試樣低端溫度達(dá)到0 ℃所需時(shí)間為50 h；當(dāng)冷端溫度為-15 ℃時(shí)，試樣低端溫度達(dá)到0 ℃所需時(shí)間為38 h，說(shuō)明隨著試樣冷端溫度的降低，試樣低端溫度降至0 ℃所需時(shí)間變短。這是因?yàn)殡S著冷端溫度的降低，凍結(jié)端面施加的負(fù)溫增量逐漸增大，單位時(shí)間內(nèi)土體凍結(jié)深度變大，故達(dá)到0 ℃所需時(shí)間變短。

基于熱力耦合模型，本文所建立的模型試驗(yàn)不考慮外部荷載作用，此處主要以凍脹量來(lái)描述熱力耦合中熱交換對(duì)力的影響。通過(guò)數(shù)值模擬計(jì)算，得到在不同冷端溫度下試樣溫度達(dá)到穩(wěn)定時(shí)的位移場(chǎng)云圖，見圖3。從圖3可知，在無(wú)外界水源補(bǔ)給的試驗(yàn)條件下，當(dāng)冷端溫度為-5 ℃時(shí)，試樣位移場(chǎng)達(dá)到穩(wěn)定時(shí)土體的凍脹量為1.49 mm；當(dāng)冷端溫度為-10 ℃時(shí)，試樣位移場(chǎng)達(dá)到穩(wěn)定時(shí)土體的凍脹量為2.34 mm；當(dāng)冷端溫度為-15 ℃時(shí)，試樣位移場(chǎng)達(dá)到穩(wěn)定時(shí)土體的凍脹量為1.91 mm。這是因?yàn)殡S著冷端溫度的降低，試樣內(nèi)部未凍水含量降低，導(dǎo)致

圖3 不同冷端溫度下試樣凍脹量云圖(單位：mm)

凍結(jié)時(shí)產(chǎn)生的冰晶體數(shù)量相應(yīng)增長(zhǎng)，故法向凍脹力增長(zhǎng)，從而使土體的凍脹量增大；但當(dāng)凍脹力達(dá)到一定值時(shí)，土體內(nèi)部的冰晶體破壞，從而土體的凍脹量減小。實(shí)際工程中，由于存在地下水的補(bǔ)給條件，凍結(jié)土體的凍脹性更明顯。

3 數(shù)據(jù)分析

根據(jù)溫度和凍脹量云圖，選取-10 ℃的各測(cè)點(diǎn)模擬值和試驗(yàn)值進(jìn)行對(duì)比。試樣測(cè)點(diǎn)溫度實(shí)測(cè)值與模擬值對(duì)比見圖4。從圖4可知，試樣測(cè)點(diǎn)溫度實(shí)測(cè)值和數(shù)值模擬值變化規(guī)律一致，距離冷端的距離越近，土體降溫速率越大；距離冷端越遠(yuǎn)，土體降溫速率越小。在凍結(jié)試驗(yàn)前10 h，測(cè)點(diǎn)實(shí)測(cè)溫度和模型模擬溫度下降速率最大，當(dāng)凍結(jié)時(shí)間達(dá)到48h后，各測(cè)點(diǎn)實(shí)測(cè)溫度和模型模擬溫度都逐漸趨于穩(wěn)定。對(duì)比試驗(yàn)值和模擬值發(fā)現(xiàn)，模擬得到的穩(wěn)定溫度值略低于試驗(yàn)值，但整體效果較好。因試驗(yàn)過(guò)程中，試樣受周圍環(huán)境的影響比較大，試樣一部分能量用于熱交換，故在測(cè)點(diǎn)處模型模擬值較試驗(yàn)值溫度低。

圖4 試樣測(cè)點(diǎn)溫度實(shí)測(cè)值與模擬值對(duì)比

試樣實(shí)測(cè)凍脹量與模擬凍脹量對(duì)比見圖5。從圖5可知，熱力耦合模擬值和實(shí)際測(cè)量值變化趨勢(shì)一致，熱力耦合模型較好地模擬了凍結(jié)過(guò)程中凍脹量的變化趨勢(shì)。在凍結(jié)時(shí)間的前16 h，凍結(jié)砂土土體凍脹量隨時(shí)間的增大幾乎呈線性增大；在16 h以后，凍脹量變化幾乎呈波浪變化并逐漸趨于定值。表明所建立的模型適用于榆林某路基砂土，所建模型合理。

圖5 試樣實(shí)測(cè)凍脹量與模擬凍脹量對(duì)比

4 結(jié) 語(yǔ)

本文通過(guò)能量守恒定律和多孔介質(zhì)傳熱理論，分析砂土在單向凍結(jié)時(shí)試樣內(nèi)部熱量分布情況；運(yùn)用固體力學(xué)理論，分析土體在自重情況下凍脹量的大小；運(yùn)用廣義胡克定律和克拉伯龍方程，考慮砂土在單向凍結(jié)時(shí)土體內(nèi)部的熱力耦合。榆林某路基砂土熱力耦合模型模擬結(jié)果和試驗(yàn)實(shí)測(cè)結(jié)果對(duì)比表明，此模型能很好地模擬土體在凍結(jié)過(guò)程中土體內(nèi)部的熱量變化和凍脹量情況，可應(yīng)用于工程實(shí)際。