凍融作用下多年凍土邊坡局部穩(wěn)定性研究

褚志成，雷勝友，原喜忠，郭長(zhǎng)志，王漢文，陳鵬輝，陳虹宇

(長(zhǎng)安大學(xué) 公路學(xué)院, 陜西 西安 710064)

隨著我國(guó)經(jīng)濟(jì)的蓬勃發(fā)展，為滿足人民日益增長(zhǎng)對(duì)美好生活的需求，我國(guó)在寒區(qū)基建的投入也逐年增加。在全球氣候變暖的大環(huán)境下，我國(guó)多年凍土分布整體呈現(xiàn)出南界北移的趨勢(shì)，寒區(qū)邊坡的失穩(wěn)滑塌問(wèn)題也逐漸增多，給人們的生產(chǎn)生活帶來(lái)了極大的隱患。究其緣由，周期性氣溫冷暖交替引起的反復(fù)凍融是造成多年凍土邊坡失穩(wěn)破壞的重要原因。在凍融過(guò)程中，冷季隨著氣溫的下降，季節(jié)活動(dòng)層土體水分在凍結(jié)鋒面積聚，而到了暖季，春融期積雪融化以及原凍結(jié)土層融化使得土體水分處于過(guò)飽和狀態(tài)[1-3]。就多年凍土邊坡而言，暖季氣溫的升高直接影響到坡體溫度場(chǎng)的分布，邊坡處于融化期，活動(dòng)層水分的增加降低了潛在滑動(dòng)面的抗剪強(qiáng)度，造成邊坡土體承載力急劇下降；當(dāng)顆?？障锻耆錆M水時(shí)，應(yīng)力的傳遞方式發(fā)生了變化，由起初的顆粒傳遞變?yōu)榭紫端畟鬟f，導(dǎo)致土體的粘結(jié)作用降低，從而降低了土體的強(qiáng)度。

目前，凍土邊坡凍融滑塌穩(wěn)定性分析已有大量的研究成果，武鶴等[4]考慮了滑坡面的端摩阻力作用，并利用基線平衡原理推導(dǎo)出了寒區(qū)土質(zhì)邊坡安全系數(shù)的計(jì)算公式。劉義高[5]利用有限差分法從土體強(qiáng)度的角度分析了西藏地區(qū)邊坡的穩(wěn)定性，并揭示了該地區(qū)邊坡的破壞機(jī)理。葛琪等[6]考慮了水分遷移對(duì)季凍區(qū)邊坡穩(wěn)定性的影響，結(jié)合有限差分法，分析了季凍區(qū)正融土坡的失穩(wěn)狀態(tài)。高檣等[7]利用有效應(yīng)力極限平衡法在不考慮凍土融化過(guò)程中的孔隙水壓力增大的影響的情況下，推導(dǎo)出了斜坡穩(wěn)定性系數(shù)的計(jì)算公式。曾韜睿等[8]考慮了融化土體的滲透力，基于修正的傳遞系數(shù)法推導(dǎo)了凍融邊坡穩(wěn)定性的計(jì)算公式，并表明融化深度是凍土邊坡穩(wěn)定性的主要影響因素。孫國(guó)棟等[9]考慮凍融循環(huán)以及滲流對(duì)寒區(qū)高陡邊坡的作用，提出了基于凍融折線型滑移面的概念，并結(jié)合框架錨桿支護(hù)邊坡得到了寒區(qū)高陡邊坡穩(wěn)定性評(píng)價(jià)方法?？煽闯?，對(duì)于寒區(qū)土質(zhì)邊坡凍融穩(wěn)定性分析中，大多采用靜力平衡法、傳遞系數(shù)法、強(qiáng)度折減等整體分析方法，而且大多未考慮凍融循環(huán)對(duì)土體的損傷作用以及凍土邊坡的局部?jī)鋈诘奶厥庑?，所以上述方法在多年凍土邊坡穩(wěn)定性分析中不能完全適用?；诖?，本文以東北大興安嶺地區(qū)某一級(jí)公路路塹土質(zhì)邊坡為研究對(duì)象，考慮凍融循環(huán)對(duì)邊坡土體的損傷作用以及考慮多年凍土邊坡的局部?jī)鋈?，建立有限元邊坡穩(wěn)定性分析方法，為多年凍土地區(qū)邊坡的穩(wěn)定性評(píng)價(jià)以及邊坡防護(hù)治理提供參考。

1 計(jì)算方法

1.1 帶有相變的非穩(wěn)態(tài)平衡方程

根據(jù)熱傳導(dǎo)和質(zhì)量遷移原理，僅考慮對(duì)流換熱和土體在凍融過(guò)程中的相變問(wèn)題，則二維非穩(wěn)態(tài)熱傳導(dǎo)方程表示為:

(1)

式中：ρ為土體密度；T表示土體溫度；q為內(nèi)熱源強(qiáng)度；在相變區(qū)為[TbTp]時(shí)，土體的導(dǎo)熱系數(shù)λ及容積熱容量C的取值如式(2)、式(3)所示：

(2)

(3)

式中：Cf、Cu分別表示凍土和融土的比熱；λf、λu分別表示凍土和融土的導(dǎo)熱系數(shù)；L為水在結(jié)晶或融化時(shí)釋放的潛熱，故相變熱Q=Lρd(w-wu)，其中wu為未凍水含量，ρd為土體的干密度。在土體融化后，未凍水含量認(rèn)為是定值，如式(4)所示[10]：

(4)

1.2 土體破壞準(zhǔn)則的選取

土體的本構(gòu)模型是邊坡凍融穩(wěn)定性分析的重要前提，本文選取D-P(Drucker-Prager)準(zhǔn)則匹配Mohr-Coulomb條件。如式(5)和式(6)所示，在平面應(yīng)變情況下，當(dāng)采用關(guān)聯(lián)流動(dòng)法則時(shí)，此準(zhǔn)則成為DP4準(zhǔn)則，較好的克服了Mohr-Coulomb屈服面存在棱角不好收斂的問(wèn)題。

(5)

(6)

式中：α、k是土體內(nèi)摩擦角與黏聚力有關(guān)的常數(shù)；c、φ是土體黏聚力和內(nèi)摩擦角。

1.3 基于有限元的熱-力耦合多年凍土邊坡穩(wěn)定性計(jì)算方法

凍融作用下多年凍土邊坡穩(wěn)定性是一個(gè)復(fù)雜的動(dòng)態(tài)平衡問(wèn)題，極端的氣候條件使得邊坡土體性質(zhì)具有明顯的時(shí)效性，并且該地區(qū)邊坡有著整體穩(wěn)定性較好，坡體淺層易受極端氣候反復(fù)凍融的影響，而永凍層相對(duì)較穩(wěn)定的特點(diǎn)?；诖耍P者考慮土體凍融過(guò)程中，環(huán)境溫度及凍融循環(huán)對(duì)邊坡土體的綜合作用，利用多物理場(chǎng)有限元軟件強(qiáng)大的函數(shù)編輯功能，編寫(xiě)了熱-力耦合動(dòng)態(tài)if函數(shù)，即利用溫度來(lái)界定凍融交界面，以此解決在數(shù)值計(jì)算中難以考慮局部?jī)鋈诘膯?wèn)題。

強(qiáng)度折減法有著不需假設(shè)滑動(dòng)面以及能夠適應(yīng)各種復(fù)雜情況的優(yōu)點(diǎn)。就邊坡而言，其本質(zhì)就是在有限元計(jì)算中，使得巖土體的抗剪強(qiáng)度參數(shù)黏聚力以及內(nèi)摩擦角不斷折減同樣的倍數(shù)，直至最終巖土體達(dá)到臨界的破壞狀態(tài)，此時(shí)的折減系數(shù)即為邊坡的穩(wěn)定性系數(shù)[11-12]，即：

cF=c/F
φF=tan-1(tanφ/F)

(7)

式中：cF、φF分別為折減后的黏聚力和內(nèi)摩擦角；F為折減系數(shù)即邊坡的穩(wěn)定性系數(shù)。

值得注意的是，上述強(qiáng)度折減法是針對(duì)邊坡整體穩(wěn)定性而言的，對(duì)僅受到局部?jī)鋈诘亩嗄陜鐾吝吰露圆荒芡耆m用。本文是在熱-力耦合動(dòng)態(tài)if函數(shù)的基礎(chǔ)上，利用強(qiáng)度折減原理僅對(duì)發(fā)生凍融作用的淺層土體不斷折減其強(qiáng)度，以坡頂位置處位移發(fā)生突變作為邊坡失穩(wěn)破壞的判斷標(biāo)準(zhǔn)，以此分析邊坡穩(wěn)定性以及計(jì)算其安全系數(shù)。

2 模型建立過(guò)程

2.1 工程概況

某一級(jí)公路位于我國(guó)東北大興安嶺內(nèi)蒙古境內(nèi)，整體呈東西走線。根據(jù)氣象數(shù)據(jù)得到該區(qū)年平均氣溫為-3℃，最大凍結(jié)深度為3.1 m，根據(jù)中國(guó)季節(jié)性凍土標(biāo)準(zhǔn)凍深線，該地區(qū)屬于多年凍土區(qū)。根據(jù)課題組在該段鉆孔實(shí)測(cè)數(shù)據(jù)顯示地表大多覆蓋0.5 m～2.0 m的草炭土，以下多為粉質(zhì)黏土、碎石土以及砂礫。人工邊坡坡高度在8 m～15 m之間，坡率大多集中在1∶1.1～1∶3之間，土質(zhì)邊坡為黏性土，這里選取坡率為1∶1.2、坡高為15 m的土質(zhì)邊坡進(jìn)行分析。

2.2 凍融循環(huán)對(duì)土體強(qiáng)度的影響

(1) 對(duì)抗剪強(qiáng)度的影響?？辜魪?qiáng)度指標(biāo)在邊坡穩(wěn)定性分析中的基本參數(shù)。目前大量的研究成果表明，土體在長(zhǎng)期的凍融作用下，其強(qiáng)度都有不同程度的損傷[13-14]。文獻(xiàn)[14]給出了土體在凍融作用下，土體抗剪強(qiáng)度指標(biāo)損傷系數(shù)與凍融循環(huán)次數(shù)的關(guān)系：

(8)

式中：kc、kφ分別是黏聚力和內(nèi)摩擦角的凍融損傷系數(shù)；n為凍融循環(huán)次數(shù)，這里指時(shí)間周期，單位為a。則得到經(jīng)過(guò)凍融循環(huán)后的抗剪強(qiáng)度參數(shù)為：

c′=ckc,φ′=φkφ

(9)

式中：c、φ分別是土體初始黏聚力和內(nèi)摩擦角；c′、φ′分別為土體損傷后的黏聚力和內(nèi)摩擦角。

(2) 對(duì)彈性模量的影響。彈性模量是巖土工程中變形和穩(wěn)定性分析中的一個(gè)重要參數(shù)。文獻(xiàn)[15]以軸向應(yīng)變?yōu)?.0%時(shí)所對(duì)應(yīng)的應(yīng)力作為土體彈性模量受凍融作用影響的標(biāo)準(zhǔn)，得到土體的彈性模量隨凍融循環(huán)次數(shù)先降低而后增大規(guī)律。文獻(xiàn)[16]針對(duì)東北地區(qū)典型黏性土通過(guò)不固結(jié)不排水三軸試驗(yàn)，得到了黏性土應(yīng)力-應(yīng)變關(guān)系。本文選取軸向應(yīng)變?yōu)?.0%時(shí)所對(duì)應(yīng)的的偏差應(yīng)力增量與軸向應(yīng)變?cè)隽康谋戎底鳛橥馏w的彈性模量，即：

(10)

式中：E表示土體的彈性模量；σ、ε表示在三軸試驗(yàn)時(shí)的應(yīng)力與應(yīng)變，腳標(biāo)表示該應(yīng)變時(shí)的應(yīng)力值。

根據(jù)上式計(jì)算得到黏性土在經(jīng)歷凍融循環(huán)后的彈性模量的變化趨勢(shì)，計(jì)算結(jié)果如圖1所示?？煽闯鲳ば酝猎诮?jīng)歷凍融循環(huán)13次后彈性模量都發(fā)生了不同程度的下降，前兩次的凍融對(duì)土體彈性模量損傷幅度最大，而后隨著凍融循環(huán)次數(shù)的增加，彈性模量表現(xiàn)出先增大后減小的趨勢(shì)?？傮w而言，彈性模量的平均降低幅值為未凍融土彈性模量的38%。

圖1 黏性土彈性模量隨凍融循環(huán)次數(shù)變化曲線

2.3 模型建立

(1) 幾何模型。通過(guò)對(duì)該段某處邊坡實(shí)際形態(tài)的簡(jiǎn)化處理建立邊坡分析模型，邊坡高度為15 m，坡率1∶1.1，重度17.8 kN/m3，初始黏聚力39.3 kPa,初始內(nèi)摩擦角20.9°，初始彈性模量20 MPa，泊松比為0.3，幾何模型如圖2所示。

圖2 邊坡幾何模型(單位:m)

(2) 參數(shù)選取。凍融作用對(duì)土體強(qiáng)度的影響如2.2節(jié)所述，土體凍結(jié)時(shí)力學(xué)參數(shù)抗剪強(qiáng)度、泊松比以及彈性模量與土溫的關(guān)系如式(11)所示[17],其熱力學(xué)參數(shù)參見(jiàn)表1和表2。

E=a1+b1|T|m
v=a2+b2|T|
c=a3+b3|T|
φ=a4+b4|T|

(11)

式中:ai、bi是試驗(yàn)參數(shù)，各個(gè)參數(shù)的取值如表1所示：在土體溫度T>0℃時(shí)，bi均為0；m取0.6。

表1 土體力學(xué)試驗(yàn)參數(shù)[18]

表2 土體熱物理參數(shù)

(3) 邊界條件。氣候環(huán)境對(duì)邊坡土體的作用主要取決于邊坡表面與大氣的對(duì)流換熱過(guò)程。根據(jù)近60年?yáng)|北地區(qū)氣候變化趨勢(shì)特征，預(yù)測(cè)東北地區(qū)年平均氣溫未來(lái)升溫速率為0.05℃/a，在考慮大氣逐漸升溫的情況下，氣溫和地表溫度隨時(shí)間變化如公式(12)所示[19]：

(12)

式中:t為時(shí)間，單位s；Ta為氣溫；Ts為地表溫度；Tm為年平均氣溫；A為年平均氣溫較差；g(t)為考慮全球變暖升溫項(xiàng)。

(4) 初始條件。初始溫度場(chǎng)對(duì)整個(gè)坡體狀態(tài)的影響較為重要，將年平均地溫作為整個(gè)坡體初始溫度場(chǎng)計(jì)算的初始值，把g(t)=0時(shí)的天然地表溫度作為整個(gè)坡體表面的溫度邊界條件，最后得到熱-力耦合階段的初始溫度場(chǎng)。

為消除自重對(duì)邊坡位移的影響，應(yīng)首先平衡地應(yīng)力，排除由土體自重對(duì)坡體應(yīng)力及位移的影響，把此時(shí)的應(yīng)力狀態(tài)作為熱-力耦合階段的初始應(yīng)力場(chǎng)。

2.4 天然溫度場(chǎng)驗(yàn)證

根據(jù)實(shí)測(cè)鉆孔數(shù)據(jù)，該段天然地層從上依次為0.3 m草炭土、2.2 m粉質(zhì)黏土、2.2 m碎石土以下為砂礫，按照實(shí)際地層情況建立80 m×40 m的矩形幾何模型，以此模擬天然地表地溫變化情況。限于篇幅，幾何模型不再展示。

圖3所示為天然地表鉆孔實(shí)測(cè)溫度與數(shù)值模擬值地溫對(duì)比曲線，由圖3可看出，模擬值與實(shí)測(cè)值隨深度的變化曲線吻合度較高，在7月30日和11月18日兩個(gè)時(shí)間點(diǎn)除了在接近地表位置有所差異(可能是接近地表位置測(cè)溫鋼管內(nèi)的溫度受到外界空氣的影響較大所致)，在凍土上限以下，基本一致，證明了本模型的適用性及合理性。也得到該段年平均地溫為-1℃，年變化深度在地表以下9 m處。

3 多年凍土邊坡穩(wěn)定性分析

3.1 邊坡凍融演變過(guò)程

(1) 溫度場(chǎng)演變過(guò)程。隨著氣溫對(duì)土體周期性的作用，坡面季節(jié)活動(dòng)層隨之經(jīng)歷著反復(fù)的融化和凍結(jié)，坡體溫度場(chǎng)演變過(guò)程中，融化深度呈現(xiàn)出先增加后減少的趨勢(shì)。在5月初出現(xiàn)淺層融化現(xiàn)象直至9月份達(dá)到最大融化深度，而后進(jìn)入冷季，氣溫的下降，冷量向坡體傳遞，原來(lái)的融化層逐漸開(kāi)始縮小，直至次年1月末融化層完全消失，之后坡體一直處于凍結(jié)狀態(tài)。如果將坡體融化深度增長(zhǎng)的過(guò)程認(rèn)為邊坡處于融化期，坡體融化深度減小的過(guò)程認(rèn)為邊坡處于凍結(jié)期，則邊坡融化期為5月—9月份，凍結(jié)期為10月—次年4月份。通過(guò)“0.00℃”等值線位置可以明顯看出坡體淺層融化夾層的存在，在10月份出現(xiàn)到次年2月份坡體中部會(huì)形成保持長(zhǎng)達(dá)5個(gè)月的融化夾層，融化夾層在10月份出現(xiàn)并逐漸發(fā)展直至次年2月份消失。

圖3 溫度實(shí)測(cè)值與模擬值對(duì)比曲線

以坡頂位置作為特征點(diǎn)，提取坡頂特征點(diǎn)以下的融化層厚度，得到邊坡在一個(gè)溫度周期的融化層厚度隨時(shí)間的變化曲線見(jiàn)圖4。由圖4可看出，邊坡融化層厚度在溫度周期內(nèi)和環(huán)境溫度變化呈現(xiàn)出很強(qiáng)的相關(guān)性，并隨著環(huán)境溫度的變化而變化，但與環(huán)境溫度在時(shí)間上存在一定的滯后性。環(huán)境溫度在4月份到7月份是升溫階段，坡體表層從5月份開(kāi)始逐漸融化直到10月份達(dá)到最大融化深度2.9 m；而后環(huán)境溫度下降，坡體出現(xiàn)融化夾層，坡體內(nèi)部融化層厚度慢慢減小，隨之進(jìn)入冷季，到次年2月份融化夾層完全消失；在此期間坡體融化層厚度經(jīng)歷先增大后減小最終完全消失的變化趨勢(shì)?？傮w上，坡體融化層厚度與環(huán)境溫度在時(shí)間上存在約1個(gè)月的滯后期。

圖4 邊坡融化深度與環(huán)境溫度隨時(shí)間變化曲線

(2) 應(yīng)力及位移場(chǎng)分析。凍融作用下，坡體應(yīng)力場(chǎng)分布表現(xiàn)出明顯的時(shí)效性。在整個(gè)溫度周期內(nèi)最大壓應(yīng)力都出現(xiàn)在邊坡底部；在邊坡融化期，拉應(yīng)力區(qū)集中在邊坡淺層，并且隨著溫度的升高，拉應(yīng)力區(qū)不斷增加；在凍結(jié)期，邊坡拉應(yīng)力區(qū)隨著溫度的降低不斷減少。而凍融循環(huán)對(duì)邊坡淺層土體強(qiáng)度的損傷作用，是邊坡在融化期發(fā)生失穩(wěn)破壞不可忽視的原因。

圖5為坡頂特征點(diǎn)豎向位移和水平位移隨時(shí)間發(fā)展的變化曲線。從位移變化情況看，在整個(gè)循環(huán)周期內(nèi)，特征點(diǎn)豎向位移和水平位移隨時(shí)間變化趨勢(shì)基本一致，整條曲線呈“臺(tái)階”狀，可分為快速發(fā)展階段和穩(wěn)定發(fā)展階段。具體的，在溫度周期循環(huán)下，隨著氣溫的逐漸升高，季節(jié)活動(dòng)層開(kāi)始融化，土體融化使得含水量逐漸增加，導(dǎo)致土體強(qiáng)度的急劇下降，此時(shí)位移出現(xiàn)第一個(gè)快速發(fā)展階段；當(dāng)邊坡處于凍結(jié)期，活動(dòng)層逐漸凍結(jié)，土體在凍結(jié)狀態(tài)時(shí)的強(qiáng)度遠(yuǎn)遠(yuǎn)大于融化狀態(tài)時(shí)的強(qiáng)度，并且隨著溫度的降低在不斷增大，在此期間特征點(diǎn)位移進(jìn)入穩(wěn)定發(fā)展階段并且豎向位移出現(xiàn)輕微的回升。如此反復(fù)特征點(diǎn)位移在不斷增加，以位移發(fā)生突變作為判定邊坡發(fā)生破壞的標(biāo)準(zhǔn)，從圖中可看出，在經(jīng)歷8個(gè)周期循環(huán)后位移發(fā)生突變，此時(shí)數(shù)值計(jì)算也出現(xiàn)了不收斂的情況，可以判定此時(shí)邊坡達(dá)到了破壞的臨界狀態(tài)。

圖5 坡頂位移隨時(shí)間變化過(guò)程

3.2 安全系數(shù)的討論

依據(jù)第1節(jié)所述多年凍土邊坡穩(wěn)定性分析方法以及3.1(2)所述，可認(rèn)為在第9年即經(jīng)歷8次凍融循環(huán)后邊坡達(dá)到失穩(wěn)的臨界狀態(tài)。即可根據(jù)強(qiáng)度折減原理得到該邊坡安全系數(shù)為1.494。

(1) 滑動(dòng)面的確定。如圖6所示為邊坡安全系數(shù)F=1.494時(shí)的位移等值線分布圖，可以明顯看出有位移突變點(diǎn)的存在，圖中粗實(shí)線為位移等值線沿著突變點(diǎn)位移增大的方向形成的實(shí)線，該方向上的特征節(jié)點(diǎn)全部都發(fā)生了突變，此時(shí)坡體已經(jīng)是失穩(wěn)狀態(tài)，故此可以以該方向作為邊坡的最危險(xiǎn)滑動(dòng)面。

圖6 F=1.494時(shí)的位移分布圖

(2) 對(duì)比極限平衡法。極限平衡法根據(jù)不同的假設(shè)條件有如表3所示的幾種方法。本節(jié)利用GEO-SLOPE極限平衡法方法來(lái)驗(yàn)證有限元方法的合理性，計(jì)算模型如圖7所示。因凍土邊坡的失穩(wěn)破壞常常發(fā)生在暖季，所以可僅考慮邊坡處于融化期的狀態(tài)。考慮凍融影響，將有限元計(jì)算得到的邊坡處于臨界狀態(tài)時(shí)的土體強(qiáng)度參數(shù)作為邊坡融化層的參數(shù)，融化層厚度按照邊坡最大融化深度設(shè)置；永凍層的力學(xué)參數(shù)按照凍結(jié)狀態(tài)設(shè)置。計(jì)算結(jié)果如表3所示，可看出極限平衡法的5種方法與本文所得安全系數(shù)相差不大，證明了本文方法的準(zhǔn)確性。

圖7 GEO-SLOPE計(jì)算模型

表3 極限平衡法計(jì)算結(jié)果

4 結(jié) 論

(1) 在氣溫周期作用下，邊坡季節(jié)活動(dòng)層經(jīng)歷著反復(fù)凍結(jié)和融化，在9月份達(dá)到最大融化深度2.9 m，坡體中部會(huì)形成保持長(zhǎng)達(dá)5個(gè)月的融化夾層，并隨環(huán)境溫度的下降而逐漸消失；融化層厚度與環(huán)境溫度有較強(qiáng)的相關(guān)性，但存在約1個(gè)月的滯后期。

(2) 邊坡坡頂特征點(diǎn)豎向位移和水平位移隨時(shí)間變化趨勢(shì)基本一致，整個(gè)曲線呈“臺(tái)階”狀，可分為快速發(fā)展階段和穩(wěn)定發(fā)展階段，直至邊坡出現(xiàn)失穩(wěn)破壞。