董旭光，董建華，何天虎，代 濤，孫國棟

(1.寧夏大學土木與水利工程學院，寧夏 銀川 750021；2.蘭州理工大學甘肅省土木工程防災減災重點實驗室，甘肅 蘭州 730050；3.蘭州理工大學理學院，甘肅 蘭州 730050)

0 引言

多年凍土占我國國土面積的21.5%，近年來在凍土區(qū)進行了大規(guī)模的土方開挖和填土工程建設，對地質環(huán)境的擾動使很多凍土邊坡穩(wěn)定性變差。如青藏公路K3057段附近發(fā)生融凍泥流滑坡，K3035段滑塌體在5年內向山頂后退約100 m，泥流滑塌體向青藏公路移動，即將影響公路的正常運營[1]。牛富俊等[2]將青藏高原凍土斜坡破壞劃為：崩塌、蠕變、熱融失穩(wěn)和泥流等。邊疆等[3]對多年凍土區(qū)的主要凍害類型和特征進行歸類， 并提出防治措施。馮守中等[4]對寒區(qū)巖質和土質邊坡進行破壞機理分析，并給出了相應穩(wěn)定性計算方法。劉紅軍等[5]對寒區(qū)土質邊坡凍融失穩(wěn)及植被護坡進行了試驗研究，監(jiān)測表明：邊坡凍結過程中，水分向凍結鋒面遷移，木本植物護坡有助于減小邊坡凍融失穩(wěn)。靳德武[6]探討了低角度凍土斜坡的破壞機理，進行了凍融模型試驗，給出了不同滲流條件下凍土邊坡穩(wěn)定性分析方法。LI等[7]對玉樹成都縣支梅村秋季牧場邊坡失穩(wěn)進行了數值模擬，發(fā)現造成失穩(wěn)的主要原因是季節(jié)性霜凍和持續(xù)降水下滲引起的基底抗剪強度降低。趙夢怡等[8]對康定縣新都橋粗顆粒凍土道路邊坡進行變形、氣溫、地溫和地下水等長期觀測表明：在凍結融化和積雪消融條件下，粗顆粒凍土邊坡變形突增，坡面表層1 m內產生滑塌。董旭光[9]針對多年凍土熱融滑塌，基于主動保護凍土理念，提出區(qū)新型框架熱錨管邊坡支護結構，模擬和驗證表明結構降溫效果良好、控制變形能力強。查閱大量文獻可知，目前對凍土邊坡研究還較少，由于凍土是一種對溫度極其敏感的含冰體，凍土邊坡病害與其中水分的凍融和遷移有密切關系。因此，凍融作用下邊坡失穩(wěn)是水熱力共同作用的結果，故基于水熱力耦合研究凍土邊坡的穩(wěn)定具有重要意義。

國內外對水熱力耦合進行了大量研究，如HARLAN[10]提出水-熱耦合模型，隨后又提出分凝勢、水熱力耦合和剛體冰模型等。陳飛熊[11]建立了全面考慮凍土中骨架、冰、水三相介質水、熱、力與變形耦合的數理方程。李洪升等[12]將凍土視為彈性均質體，推導了凍結過程中水分、溫度和應力耦合控制方程。LI等[13]在考慮凍土流變的彈黏塑性本構模型基礎上，建立了描述土體凍融過程的水-熱-力相互作用理論模型，并對某水渠進行了計算。MA等[14]討論了冰-水相變廣義Clasusius-Clapeyron方程的適用條件，并提出了更能真實反映土體凍結過程的動態(tài)模型。這些研究成果為分析凍土邊坡穩(wěn)定提供了堅實基礎，然而水熱力多場耦合計算過程復雜，傳統的軟件還無法直接模擬土體“凍脹融沉”現象，需要通過等效簡化，為工程問題分析帶來不便。因此，開發(fā)多場耦合分析軟件對快速計算凍土邊坡凍融反應特性有重要價值。

本文基于水熱力三場耦合理論，編寫了能反應土體凍脹融沉的多場耦合計算程序，開發(fā)了凍土邊坡水熱力耦合分析軟件，建立了多年凍土邊坡多場耦合計算模型，分析了凍融作用時邊坡的溫度、水分、應力和變形分布規(guī)律，為凍土邊坡工程計算提供方便及參考。

1 凍土邊坡水熱力耦合控制方程

在周期性凍融循環(huán)作用下，凍土骨架和水的熱傳導及冰-水相變作用，則凍土邊坡中溫度滿足帶相變瞬態(tài)溫度場的熱傳導方程為[11]：

(1)

式中:C——體積比熱容；

T——溫度；

t——時間；

x、y——坐標系；

λ——導熱系數；

L——相變潛熱；

ρi——冰的密度；

θi——冰的體積含量。

凍融作用下凍土邊坡中水分遷移以液態(tài)水進行，則水分遷移方程為[12]：

(2)

式中：θu——未凍水體積含量；

D——水分擴散系數；

ρw——水的密度。

凍土未凍水含量與溫度的聯系方程為：

(3)

式中：a、b——與土體性質有關的參數；

θw——總含水量。

將式(1)代入式(2)，并由(3)式的導出關系化簡得到水熱耦合方程為：

(4)

采用顯熱容法來考慮相變面，假設土體相變發(fā)生在Tm附近的一個溫度范圍內(Tm±ΔT)，則簡化構造出體積熱容和導熱系數表達式為[15]：

(5)

(6)

式中：Cf——凍土比熱容；

λf——凍土導熱系數；

Cu——融土比熱容；

λu——融土導熱系數；

Tm——相變區(qū)間中點溫度；

ΔT——相變區(qū)間中點溫度與最高或最低相變溫的絕對差值。

凍土邊坡的靜力平衡微分方程為：

[?]{σ}-{f}=0

(7)

物理方程為[12]：

{σ}=[DT]({ε}-{εv})

(8)

幾何方程為：

{ε}=[?]T{u}

(9)

式中：[?]——矩陣微分算子；

{σ}——應力；

{f}——體力；

[DT]——與溫度有關的彈性矩陣；

{ε}——彈性應變；

{εv}——凍脹或融化引起的應變；

{u}——位移。

其中，

凍結時，

融化時，

式中：E——彈性模量；

μ——泊松比；

θ0——初始含水量；

Δθ——在一定時間內的水分遷移量；

n——土的孔隙率；

A——融沉系數；

a1、b1、a2、b2、m——與土性有關的參數，取值參見文獻[15-16]。

方程(1)～(9)組成了多年凍土邊坡的水熱力耦合方程組。由于土體熱力學參數隨溫度而變化，并在計算模型區(qū)域內有一個隨時間變化的凍融界面，凍融界面產生相變，因此，上述方程描述的是一個大空間大時間尺度的強非線性問題，無法獲得解析解，故采用數值方法進行求解，本文將凍土邊坡控制方程在空間域采用有限元離散，時間域采用有限差分離散，兩者相結合迭代求解，基于Matlab平臺開發(fā)多年凍土邊坡多場耦合計算軟件。

2 程序驗證及軟件開發(fā)

2.1 文獻程序驗證

圖1 PENNER實驗試樣的有限元數值分析模型Fig.1 FEM of PENNER experiment

為了檢驗編制的水熱力耦合程序的正確性，采用PENNER[17]研究凍土冰透鏡體形成過程實驗和MU等[18]對該實驗數值計算，驗證本文水熱力耦合程序的正確性。土體熱學、力學參數和邊界條件均與文獻相同，分析模型見圖1，溫度邊界見圖2。利用開發(fā)的程序計算，選取凍結深度和凍脹量進行對比分析，結果如圖3和圖4所示。

試驗用土的含水量為0.35，未凍水和土體負溫之間的函數關系[18]：

凍土視為線彈性體，凍土和未凍土泊松比均取0.3，未凍時彈性模型11.2 MPa，凍結時彈性模量為：

(11)

土體導水系數用Horiguchi實驗結果：

(12)

溫度初始和邊界條件：土體頂部為暖端；底部為冷端，進行由下向上凍結。邊界條件如圖2所示[17-18]。為了模擬凍結前試樣的初始狀態(tài)，將t=0時刻的上邊界和下邊界溫度作為初始溫度，利用初始溫度條件和土體的熱參數及其它特性參數求解Laplace方程，求得計算區(qū)域溫度場的估計值，修正參數，再繼續(xù)計算直至穩(wěn)定溫度場作為初始溫度場。

圖2 溫度邊界條件Fig.2 Boundary conditions of temperature

水分邊界條件：土體內為飽和土體，凍結時從底部邊界(從暖端)進行水分補給。

位移邊界條件：土體底部約束，左右兩側為法向約束，頂部為自由位移邊界。

圖3 凍結深度隨時間變化的曲線Fig.3 Freezing depths at different time

圖4 凍脹量隨時間的變化曲線Fig.4 Variation of heave deformation with time

由圖3和圖4可知，本文得到不同時刻的凍結深度比實驗值略小，而比MU[18]計算結果更接近試驗值；各時刻凍脹量比實驗值稍小，與MU[18]計算相比前半段偏大，后半段偏小。綜合可知，本文得到的凍結深度和凍脹量的變化規(guī)律與PENNER實驗MU分析結果規(guī)律相符，數值相近，說明編制的程序是可靠和有效的。

2.2 軟件開發(fā)

依托MATLAB[19]平臺，基于水熱力耦合理論，開發(fā)了可獨立運行、面向對象、界面友好和操作簡便的多年凍土邊坡水熱力耦合分析軟件。該軟件包括前處理器(模型尺寸、土體物理參數的輸入和網格劃分)，求解器(迭代求解)和后處理器(結果保存和云圖顯示)，開發(fā)流程圖如圖5所示，軟件主界面、模型生成、參數輸入和計算結果顯示界面見圖6～圖7。該軟件能高效、準確實現外界溫度作用下多年凍土邊坡的溫度、水分、應力和位移計算，對凍土邊坡工程計算有重要的理論意義和應用價值。

圖5 軟件程序設計框圖Fig.5 Design diagram of software

圖6 軟件主界面Fig.6 Main interface of software

3 算例分析

3.1 計算模型

青藏高原大阪山地區(qū)某多年凍土邊坡，坡高5.0 m，坡度45°。該地區(qū)氣溫變化曲線為Ta=-3+12sin(2πt/8 760+π/2)，土層為單一均質飽和黏土，土體物理參數見表2，土體干密度ρd/(kg/m3)，初始含水量θo，凍土未凍水含量與溫度關系參數a和b；彈性模量，泊松比與溫度關系參數a1、b1、a2、b2和m；融沉系數A；熱學參數見表3，融土熱容Cu/(J·kg-1·℃-1)，導熱系數λu/(J·m-1·℃-1)和水分擴散系數Du/(J·m-3·℃-1)；凍土熱容Cf/(J·kg-1·℃-1)，導熱系數λf/(J·m-1·℃-1)和水分擴散系數Df/(J·m-3·℃-1)，有限元計算模型見圖8。

圖7 溫度云圖繪制界面Fig.7 Interface of the temperature nephogram

圖8 有限元計算模型Fig.8 Calculation model of finite element

3.2 邊界條件和初始條件

溫度邊界：根據實測結果，坡面為Ts=0.7+13sin(2πt/8 760+π/2)，坡腳和坡頂為Ts=-1.5+12sin(2πt/8 760+π/2)，兩側絕熱，底邊熱流密度為0.06 W/m2。

水分邊界：不考慮水分與外界蒸騰和入滲交換。

位移邊界：兩側邊為水平支撐，底邊為水平和豎向支撐，其余各邊自由。

求解過程：含水率取凍土邊坡實際值，然后給溫度賦予初始值，結合邊界求解式(4)和(2)，進行凍融循環(huán)計算，直至季節(jié)活動層以下的溫度和水分保持穩(wěn)定，同一時刻各節(jié)點值逐年相同為止，以此時各節(jié)點的溫度、水分作為凍土層穩(wěn)定的初始值，再求解式(7)～(9)得到初始應力，重新計算外界季節(jié)凍融作用下，邊坡土體的響應。

3.3 結果分析

邊坡經歷一個凍融周期，凍結期(寒季)從10月中旬至來年4月中旬(取4月15日)，融化期(暖季)從4月中旬至10月中旬(取10月15日)。凍結期和融化期結束是邊坡在凍融作用下的兩個極端狀態(tài)，選取這兩個時刻進行凍融分析。

表2 邊坡土體力學參數

表3 邊坡土體熱學參數

(1)凍融作用下邊坡溫度場分析

由圖9可知，凍融作用下邊坡的活動層發(fā)生凍融交替變化，不同時刻邊坡溫度場有明顯差別。凍結結束時，邊坡處于完全凍結狀態(tài)，坡面負溫較大的厚度比坡頂和坡腳均薄。融化結束時，坡面融化厚度比坡頂和坡腳都大,凍融對坡面溫度影響較大，熱融可引起邊坡滑移，此與實際觀測相符。

圖9 不同時刻邊坡溫度云圖Fig.9 The temperature of slope at different times

(2)凍融作用下邊坡水分場分析

從圖10可以看出，凍結期結束時，活動層土體的未凍水含量比其他部位低，未凍水含量與溫度同步變化，大小由溫度和土性共同決定。融化結束時，活動層土體的未凍水含量明顯比其他部位高，在凍融交界面處最高，含水量最大值大于初始含水量20%。融化期結束時未凍水含量略大于凍結時，說明凍融作用下土體發(fā)生了水分遷移，總含水量增加。

圖10 不同時刻邊坡的未凍水含量Fig.10 The unfrozen water of slope at different times

(3)凍融作用下邊坡剪應力分析

圖11 不同時刻邊坡剪應力云圖Fig.11 The shear stress of slope at different times

由圖11可知，無論凍結，還是融化結束時，坡體凍融交界面附近剪應力達到最大且呈帶狀分布，條帶與坡面大致平行，但在坡頂和坡腳附近發(fā)生轉折、小范圍延伸而未貫通坡頂和坡腳，其余部分剪應力較小且分布均勻。凍結時坡面的剪應力大于融化時。原因在于：坡面上限以下土體常年處于凍結穩(wěn)定狀態(tài)，以上土體處于凍融交替狀態(tài)，上限附近水分積聚，在自重作用下，坡面上限為凍土邊坡潛在滑動面，故剪應力遠大于其他部位。凍結時土體的彈性模量大于融化時，導致凍結時剪應力較大。

(4)凍融作用下邊坡水平位移分析

由圖12可知，凍結期結束時，坡面水平位移大小基本相同，中下部稍大，水平位移主要由土體凍脹引起；融化期結束時，邊坡水平位移為上小、下大，呈“凸肚”狀。邊坡頂部受雙向凍、融影響，頂部位移略大于頂部之下小范圍土體的位移。凍結時邊坡的水平位移比融化時小很多，原因在于：凍結時(即冬季)，土體強度參數提高，且形成了凍結土骨架，主要為凍脹位移，坡體變形較小，穩(wěn)定性較好；熱融時(即暖季)土體強度參數降低，未凍水含量提高，邊坡在自重和熱融作用下穩(wěn)定性降低，可能發(fā)生滑移，工程設計時應重視暖季熱融失穩(wěn)。

圖12 不同時刻邊坡水平位移Fig.12 Horizontal displacement of slope at different times

4 結論

(1)基于水熱力耦合理論，依托MATLAB編制了能反應凍融作用下土體多場耦合機制的有限元程序，并與經典的PENNER[17]試驗和MU等[18]模擬結果進行了對比分析，驗證了程序的可靠性，開發(fā)了多年凍土邊坡水熱力耦合分析軟件。

(2)凍融作用下邊坡的活動層發(fā)生凍融交替變化，不同時刻凍土邊坡溫度場有明顯差別，與坡面大致平行的凍融交界面出現剪應力最大條帶。凍結期結束時，坡面水平位移大小基本相同，中下部稍大；融化結束時，邊坡水平位移為上小、下大，呈“凸肚”狀。融化結束時邊坡的未凍水含量，水平位移比凍結時大，剪應力最大值比凍結時小；暖季多年凍土穩(wěn)定性較差。

封面照片說明