凍融-滲流耦合作用下土質邊坡穩定研究進展

汪恩良,靳婉瑩,商舒婷,趙 曦

(東北農業大學 水利與土木工程學院，黑龍江 哈爾濱 150030)

研究表明，飽和-非飽和滲流是誘發邊坡失穩的一個重要原因。經過多次凍融循環后，土體物理力學性質以及滲透特性發生了很大的變化。特別是在春融期，已融土和凍土層界面存在的一個高含水量、冰水混合的區域，水分在此處不斷富集，降低了土體抗剪強度、加大土體融化深度、降低斜坡的整體穩定性，引起滑坡發生。為了對春融期滑坡機理進行研究，文章分別從凍融循環對土體滲透系數的影響，飽和-非飽和土滲流理論，降雨及融雪入滲過程，以及凍融過程中水-熱-力多場耦合等多個方面進行了大量的研究和探討，在前人研究的基礎上進行總結與分析。

1 凍融循環對土體滲透系數的影響

1.1 飽和及非飽和土滲透系數

根據土體的飽和程度，一般將土體分為飽和土與非飽和土。在進行凍融土的滲透試驗研究時，飽和土滲透系數一般采用“變水頭滲透試驗儀器”進行測量，滲透系數一般采用式(1)進行計算[1]。對于非飽和土的滲透特性，一般采用實測方法或用大量的試驗數據確定經驗參數，然后制定出土水特征曲線再得出非飽和土的滲透系數[2-4]。目前，對非飽和土的物理力學性能實驗主要采用GDS三軸儀，它可以同時控制孔隙氣壓、孔隙水壓、軸向和徑向參數，很好的完成大多數線性變化的非飽和土試驗[5-6]。當前常用的測量非飽和滲透系數的模型有Mualem模型、Genuchten模型以及Fredlund&Xing模型等[7-10]。而最常用的預測非飽和土滲透系數的解析函數是Van Genuchten-Mualem[11]，如公式(2)所示。

(1)

(2)

式中：Ks為飽和土壤導水率,cm/s；K(h)為非飽和土壤導水率,cm/s；a,A分別為變水頭試管與試樣的截面積,cm2；h1,h2分別為起始水頭和終止時水頭,cm；L為試樣高度,cm；t為時間,s；Se為非飽和土的有效飽和度；m為土壤特征曲線的經驗參數。

1.2 凍融對土體物理力學性質影響

作為一種強風化作用，凍融循環使土體顆粒重新排列嚴重影響土體內部結構[12-14]。對重塑土試驗表明，凍融循環顯著改變土體各層含水率。即凍融循環后土層上部含水量增加，下部含水量減少，而在試樣的中間部位會出現含水量很大的土層[15-16]。對重塑土試驗研究表明，凍融作用對土的干密度具有雙向作用，即松散土密度增大，而密實土密度降低，多次凍融循環后，土體干密度會趨向于一個特定值，而這一特定值與土的種類有關[17-18]。研究表明無論凍融循環后土體密度如何變化，其滲透系數都會呈現增大趨勢。多數研究都發現土體滲透系數在經過凍融循環后，會增大1～2個數量級，土的滲透性和密度經3～5個凍融循環趨于穩定[19-24]。

2 飽和-非飽和滲流理論研究

2.1 飽和-非飽和滲流方程

長期以來，眾多學者針對飽和滲流研究相對較多，由于試驗條件限制以及數值模擬技術的不足，非飽和滲流研究較淺。飽和-非飽和數值分析的基本理論是建立在廣義達西定律式和質量守恒方程上的Richards方程式，如式(3)和式(4)所示[25]。

v=kJ

(3)

(4)

式中：v為滲透流速,cm/s；k為滲透系數,cm/s;J為滲透坡降；θ為土體含水量；t為時間,s；ψ為土水勢,kPa。

對于非飽和滲流計算，其參數的確定極其重要。1980年Van Genuchten等[26]提出了VG模型，隨后嚴飛等人[27]在2004年提出“雙參數VG模型”用來確定非飽和滲流參數。在無法開展試驗的條件下，可以選用函數曲線來替代非飽和滲流試驗所測得的數據線性對應關系，用VG模型參數的代表值進行計算，計算式如公式(5)和公式(6)。

(5)

(6)

式中：α,n,m為土壤特征曲線的經驗參數；h為壓力水頭,cm；θ為土壤含水率；θs為土壤飽和含水率；θr為土壤殘余含水率。

2.2 飽和-非飽和滲流研究方法

對飽和-非飽和數值分析研究方法中，常用的方法包括有限差分法和有限單元法。隨著計算機技術的進步，基于有限差分法的FLAC軟件也在飽和-非飽和滲流分析中得到了越來越廣泛的應用[28]。廖紅建[29]以Morgenstern-Price條分法和極限平衡原理為基礎對壩坡產生滑移的最小安全系數變化規律進行了分析，從而得出滲流-滑坡的內在聯系。吳良驥等[30]提出的有限差分法數值模型以及相應的應用程序，這種方法不僅具有Narasimhan提出的積分有限差分法的優點，并且應用了辛甫生數值積分提高了質量平衡精度。朱學愚等[31]將SUPG有限元方法應用于非飽和流動的Richards方程，并導出了數值表達式，對一維入滲問題進行了數值計算。陳虹等[32]提出一個簡易的曲線擬合方法，從有限的觀測試驗中獲取所需參數，這種程序可以應用到具有復雜地質邊界條件的實際滲流工程。隨著試驗方案以及試驗儀器的不斷改進，人們對于飽和-非飽和區中滲流問題的數值模型的研究愈加趨向于完善。

3 降雨及融雪入滲模型研究

3.1 降雨入滲理論模型研究

一般情況下春融期降雨量較小，但春融期坡面土壤尚未完全融化，由于土壤的單項凍結雙向融化的性質，在邊坡內部會存在一個解凍不完全、高含水率、滲透性差的薄弱帶，雨水入滲很容易破壞邊坡的應力平衡體系，此時的降雨很容易導致土壤流失或者滑坡的發生。對于凍土邊坡，土體內大部分水是以冰的形式存在的，所以在降雨初期入滲量較少，溫度較高的雨水入滲使得固態冰融化，降雨入滲率隨之增大，直到邊坡土體達到飽和產生地表徑流，降雨入滲量隨之減小，然后趨于穩定。

降雨型邊坡失穩數值模擬研究的關鍵是從理論上揭示雨水入滲后邊坡應力的變化過程、雨水在邊坡中的滲透特性和滲透過程，以及暫態附加水荷載問題等[33]。王佳佳[34]等利用GIS平臺得到滑坡土體飽和因子在時間和空間上的分布情況，并進行了水文模型與無限斜坡模型的耦合。Chue Y S[35]利用GIS平臺描述了滑坡潛在滑坡圖，并研究了極端降雨事件對斜坡滑坡和滑坡發育特征的影響。唐揚等[36]以Mein-Larson降雨入滲模型為基礎，通過假設初始土體含水率隨深度(坡面垂直的方向)呈線性分布，推導出一種適用于初始含水率不均勻分布的滑坡降雨入滲函數。鐘佩文等[37]研究得出降雨入滲使土體的抗剪強度顯著降低，邊坡穩定安全系數隨降雨入滲的進行而逐漸減小。王建新[38]通過對降雨非飽和入滲過程的水勢描述及理論模型研究，提出一種新的邊界條件處理方法，并推導出降雨自由入滲和壓力入滲模型。常用的降雨入滲模型有Green-Ampt積水入滲模型、Horton入滲模型以及著名的Philip入滲公式[39-40]。其中Green-Ampt入滲模型可表示為：

i=Ks(1+Hm/Lf)

(7)

(8)

式中：i為土壤入滲率,cm/s；Hm為濕潤峰處的平均吸力,m；Lf為概化的濕潤峰深度,m；hd為進氣吸力。

3.2 融雪入滲理論模型研究

與降雨不同，積雪融水引起邊坡失穩是一個長時間的過程，而融化后的雪水絕大部分入滲地表松散沉積物，徑流較少。影響融雪入滲的因素有很多，除溫度外，地表下墊面情況也是影響融雪入滲的一個重要因素。融雪入滲是一個復雜的過程，寒區土壤包氣帶融雪入滲過程可分為四個階段，即積雪層入滲、包氣帶融層入滲、包氣帶凍層入滲、凍層至飽水帶的入滲[41-42]。逄淑然[43]對寒區春季融雪入滲規律進行監測試驗與分析，研究了不同階段融雪水入滲深度與入滲特征。魏召才[44]通過能量平衡方程和水量平衡方程數值模擬，建立了一個具有物理基礎的基于單點雙層融雪模型；Lilbak G等[45]提出一種關于積雪離子入滲的理論模型，該模型表明離子濃度與浸潤體積之間的關系是非線性的，具有正協方差。俞鑫穎等[46]、房世峰等[47]基于遙感(RS)和地理信息系統(GIS)技術基于能量平衡和水量平衡原理對融雪產匯流進行模擬計算，提高了模擬的精度；余文君等[48]、孟現勇等[49]以SWAT模型(土壤和水評估工具)改善了融雪徑流的計算方法，提高該物理模型整體模擬精度。孟現勇等指出能量平衡融雪模型參數較多、模型復雜但精度較高，溫度因子法融雪模型結構簡單、合理確定參數后效果更好；付強等[50]、李天霄等[51]對融雪過程進行分析，得出積雪阻礙了環境因素對于土壤水熱遷移過程的影響，在一定時間段內，融雪入滲會有延后現象。錢曉慧等[52]對寒區凍融土壤的各項參數進行監測，指出邊坡滲透系數以及邊坡安全系數與積雪深度、邊坡坡度和氣溫溫度相關。譚娟[53]針對多種不同水保措施坡面觀測其融雪、侵蝕特征，指出對于融雪侵蝕，生物措施>工程措施>耕作措施。

4 多場耦合模型研究

凍土體在凍融過程中的熱量傳輸、水分遷移與相變過程并不是由單獨某個因素所造成的。凍融循環過程中土體的應力與變形計算與溫度和含水量有著密切的關系。應力、變形與水、熱是相互影響的，在進行凍融對渠道邊坡穩定性分析計算中應考慮水、熱、力的耦合計算[54]。

武建軍等[55]做了飽水狀態下凍結土體的滲流、溫度、應力耦合機理分析。金棟[56]基于水、熱、力耦合的FLAC 3D溫度場分析基本原理，模擬寒區邊坡的凍融循環過程，得到相應的安全系數。徐軼等人[57]借助COMSOL對三維滲流場進行數值分析，綜合考慮滲流、應力、溫度，實現了三場(THM)耦合。許孝臣等[58]則針對單裂隙滲流實現了滲流-應力-溫度-化學四場耦合過程，模擬了地應力、溫度、滲透性較高，水化學環境復雜條件下的裂隙滲透特性的變化規律。王世梅等[59]利用GeoStudio seep/w與slope/w兩個模塊分別進行滲流分析與穩定性分析，得到了水位升降速率和邊坡穩定性的關系。何敏等[60]研發出一種可以用于凍土中水、熱、力三場耦合分析的平臺3GEXFEM，該模型能全面考慮凍土中土骨架、冰、水三相介質水、熱、力與變形真正的耦合作用的數理方程，分析得到的溫度場、水分場與變形場與試驗結果較一致，具有較好的應用前景。

在實際工程中，多場耦合數值模擬難度仍然較大，尤其是在春融過程中，上層已融土和下層未融土交界面處的物理力學性能尤其復雜。這要求我們不斷地改進實驗方案并在試驗設計初期考慮更多的因素，以得到更加完善準確的模擬參數進行多場耦合模型的研究。

5 基于強度折減法土質邊坡穩定研究

5.1 強度折減法原理

由于土質邊坡發生失穩破壞過程，往往在邊坡內部形成一塑性區，隨著塑性區的不斷擴大，貫穿整個邊坡，邊坡隨之失穩。因此目前土質邊坡穩定性分析方法主要基于巖土的彈塑性理論，考慮邊坡土質的應力應變狀況進行分析。有限元強度折減法在邊坡穩定性分析中有著非常廣泛的的應用，它可以利用多種概化模型邊坡求解邊坡的安全系數[61]。強度折減法的基本原理如公式(9)所示：

c′=c/Ftrial,

φ′=arctan(tanφ/Ftrial)

(9)

式中：c,φ分別為黏聚力和內摩擦角;Ftrial為強度折減系數;c,φ分別除以Ftrial得到一組新的c′和φ′然后作為新的參數進行運算，直到運算不收斂，此時對應的Ftrial被稱為坡體的最小穩定安全系數。

5.2 強度折減法在邊坡穩定性研究的應用

寇天等[62]基于ABAQUS軟件，通過動態非等比例的雙強度折減法理論，研究出一種收斂快速、邊坡穩定性系數和位移均小于傳統強度折減法的動態非等比例的雙強度折減法。王俊杰等[63]基于強度折減有限元的非均質土坡失穩判據分析了非均質土坡失穩破壞的過程，提出了以特征點位移矢量突變結合塑性區貫通的失穩判據。江勝華等[64]提出針對特征點位移突變的變步長折減法，探討了強度折減系數的折減幅度，并對基于位移變化率的強度折減有限元法給出了相關的建議。張永明等[65]基于強度折減有限元法的土質邊坡穩定性分析，將強度折減有限元法與ABAQUS軟件相結合，借助ABAQUS建立了平面破壞型邊坡有限元分析模型。趙煉恒等[66]基于雙強度折減策略的邊坡穩定性分析方法探討，研究結果表明，邊坡穩定性影響因素眾多，抗剪強度參數c′和φ′的折減系數之間不存在唯一確定的函數關系。王創業等[67]通過強度折減法在邊坡穩定性中的對比分析，研究表明雙強度折減得到的臨界滑動面范圍要比等比例折減的范圍大，在不同安全系數定義中采用c-tanφ曲線分析得到的結果，物理意義更明確、更符合邊坡實際情況。用強度折減法分析邊坡穩定性時，不需要事先假定滑動面的位置和形狀，因此具有較強的適用性。

6 結論與展望

(1) 凍融作用一定程度上破壞土體的結構，使土體變疏松，滲透性增強，在凍融循環過程中，周期性的凍結與融化持續改變著土體內部結構，使土的孔隙比與滲透性發生相應的變化。凍融循環后土樣孔隙比隨凍融次數增加呈指數增加趨勢，但經過多次凍融循環土體結構強度逐漸趨向于穩定的殘余強度，這時土體的孔隙比與滲透系數逐漸趨于穩定。

(2) 飽和-非飽和滲流參數大多數是通過實驗室試驗測得的，這與實際工程中測得的參數不完全一致，針對這方面還需進行更多的現場測量以得到更多的實測資料，使計算結果更加的具體可靠。對于室內模型試驗時，需要結合現場勘測的工況進行分析。目前對滲流場的研究大多數是對飽和滲流場進行研究，而沒有進一步研究凍融期非飽和滲流場，該項工作有待后續進一步加強。

(3)春融期融雪入滲引起邊坡失穩涉及到水動力學、土工程學、計算機技術等多個學科。融雪過程復雜且影響因素極多，不同時間段融雪入滲過程不盡相同。在進行模擬計算的同時，更應該結合現場實測資料，針對不同實際工況進行分析。

(4)凍融循環過程中應力、變形與水、熱之間相互影響。目前寒區渠道多場耦合數值模擬難度仍然較大，尤其是在春融期，凍土的滲透系數很小，特別是飽和狀態下凍土滲透系數可忽略不計。春融期坡面土壤內部會存在一個解凍不完全、高含水率、滲透性差的薄弱帶，已融層水大部分受重力作用積蓄在不透水界面上，上層已融土和下層未融土交界面處的物理力學性能尤其復雜，這也是未來需要解決的又一難題。

(5)研究表明渠基土體在凍融循環后結構損傷，或者組成顆粒發生重新排列，導致土體抗剪強度發生變化。但c,φ值變化不是同步的，未來研究可以著重于抗剪強度凍融損傷因子的有限元不等比例折減系數法計算c′和φ′值，根據強度折減法分析邊坡穩定性，分析渠基土體在凍融循環后結構損傷量以及土體抗剪強度發生的變化。