凍融循環作用下露天煤礦內排土場土體力學特征及強度劣化機理

張合勇，王雪冬，朱永東，王海鵬，漆利輝

(遼寧工程技術大學 礦業學院，遼寧 阜新 123000)

我國露天煤礦的采礦規模持續擴大，開采過程中產生的棄土礦渣形成了大量平臺緊實、邊坡松散高陡的排土場[1]。由于我國大部分露天煤礦位于高寒生態脆弱地區，受凍融作用影響，土體的力學性質會隨著溫度和含水率的變化表現出復雜的動態特性[2]，常導致排土場凍脹、融凍滑塌、剝蝕等變形破壞。因此，研究凍融循環作用下土體力學特征及強度劣化機理對礦山安全生產和生態修復等工作具有十分重要的價值。

近年來，凍融循環作用影響土體物理力學性質的研究主要集中在不同凍融循環次數中應力-變形、抗剪強度、黏聚力和內摩擦角等方面的變化規律[3]。隨著凍融次數的增加，土體的抗剪強度和黏聚力多數表現出先降低后保持穩定或先降低后增加而后保持穩定的規律[4-5]，內摩擦角的變化則沒有明顯的規律性；應力-應變曲線表現出不受影響或由軟化型向硬化型過渡的變化[6-7]。為了研究凍融循環對土體內部的影響，采用核磁共振測試、掃描電鏡觀測和壓汞試驗等[8-9]，從微觀角度觀察到土體的變化不只有土顆粒的凍縮和水的相變[10]，還存在水分運移。在水分場、溫度場的相互作用及變化規律的控制下，土體原本的穩定結構被破壞，造成強度劣化。

國內外學者通過對凍土內部水分場、溫度場和應力場的研究，分析引起土體強度劣化的機理，建立了多種條件下的理論模型[11-13]，但理論模型難以說明水、熱、力三場相互作用下土體的劣化過程，而數值模擬則可以彌補理論模型的不足[14]。考慮到離散元法可以通過顆粒流理論和有限差分法等思想較好地從細觀角度模擬土的力學行為和內部結構的相關變化，桑宏偉等[15]采用離散元法實現熱量在顆粒間的傳遞，模擬出能源樁與周圍土層換熱的過程。尹楠等[16]基于離散元法模擬不同圍壓、不同溫度下凍土的靜三軸試驗。Le Tiancheng 等[17]基于離散元軟件，模擬出土樣蒸發、收縮和開裂的過程，均取得了良好的研究成果。因此，筆者選擇離散元法模擬水分場、溫度場和應力場的相互作用來分析凍融循環作用下土體強度的劣化機理。

筆者以內蒙古元寶山露天煤礦內排土場典型黏土為研究對象，通過進行天然含水率下不同凍融循環次數的室內試驗，獲取凍融作用對土料宏觀力學性質的影響規律。采用離散元數值模擬，研究凍融循環中溫度場、水分場的變化導致土體強度劣化的影響機理，以期為凍區排土場邊坡穩定性分析和其他工程建設提供一定的參考。

1 試驗材料與試驗方案

1.1 試驗材料

試樣所用土料取自內蒙古元寶山露天煤礦東南幫內排土場，取土地理位置及顆粒級配曲線如圖1 所示。按照試驗規程[18]進行相關室內土性試驗得到土料的天然干密度為1.63 g/cm3，天然含水率為17.6%，液限和塑限分別為33.8%、19.3%，塑性指數為14.5。通過塑性指數判定土料類型為黏土，根據顆粒級配、細顆粒含量和液塑限認為該土料屬于凍融敏感性土[19]。

圖1 取土位置和顆粒級配曲線Fig.1 Sampling location and particle gradation curve

1.2 試樣制備

由于排土場是人工堆積而成的松散體，具有質地疏松、強度較低等特點。土料的離散性較大，為獲取具有代表性的土樣，在多個平臺上收集一定量的原狀土料，通過室內試驗測得土料的天然含水率和天然干密度，依據上述條件進行重塑，制備直徑39.1 mm、高度80 mm 的三軸試樣。為了探究無外界水分交換條件下凍融作用引起的土體內部水分運移及相變對力學性質的影響，試驗采用保鮮膜將試樣包裹后再放入密封袋，模擬無外界水分交換的封閉系統。具體重塑和制樣步驟及過程如圖2 所示。

圖2 土料重塑、制樣步驟及圖片Fig.2 Soil remolded and sample preparation steps and pictures

1.3 試驗方案

試驗采用的凍融設備包括冰箱和恒溫恒濕箱。冰箱最低凍結溫度為-40℃；恒溫恒濕箱溫控范圍0～50℃，精度±0.1℃。考慮到研究區2011-2022 年中最低的月平均氣溫為-18℃，本試驗條件設置為：凍結溫度-20℃、融化溫度+20℃，凍結和融化時間均為12 h，即24 h 為一個周期。試驗共設置7 組，每組設置3 個試樣及1 個備用，共需28 個，循環次數分別為0、1、3、5、10、15、25 次。試樣達到凍融次數后采用TSZ10-1 型三軸儀進行不固結不排水試驗，試驗設置圍壓分別為100、200 和300 kPa，軸向應變速率設為0.8 mm/min，當應變量達到15%時停止剪切，取峰值為最大破壞強度。

2 凍融循環試驗結果與分析

2.1 對試樣質量和宏觀變形的影響

試樣在凍融過程中質量的變化從圖3 可以看出，試樣質量的損失率很低，雖然呈緩慢增長，但最終損失率僅為0.124%，即損失0.23 g 水分，說明在試驗中采用的密封處理方法有效避免了水分流失對試驗結果的影響。通過高度和直徑變化可以看出，首次凍融后兩者的凍縮均相對明顯，分別為0.095%和0.347%；高度的凍縮率在前期增長較快，至第5 次時增加到0.276%，隨后變化減緩并基本穩定在0.3%；而直徑凍縮率經過首次顯著變化后，開始緩慢增長并逐漸趨于穩定。經分析，在凍融循環試驗中，凍結時土骨架收縮，因土體含水率較低，此過程以孔隙體積被壓縮為主；融化后土骨架的收縮無法完全恢復，所以宏觀上表現為凍縮，此外，土在凍結過程中形成的冷生構造會產生微裂隙[20]，且試樣在試驗時為平放，受到重力影響，所以直徑收縮率大于高度收縮率；而高度方向的變化受重力影響較小，主要是凍融作用的結果，說明該土體受凍融作用的影響主要體現在前期。首次體積凍縮率最大，達到0.786%，隨著凍融次數增加，增長逐漸減緩并趨于穩定；可以看出，由于土體的凍融敏感性，首次凍融作用下土體結構、顆粒間的膠結和孔隙變化幅度最大[21]，造成的體積縮小最為明顯，隨著凍融次數增加，凍融作用的影響減小，內部結構逐漸調整完成，體積趨于穩定。

圖3 試樣的質量和宏觀變形與凍融次數的關系Fig.3 The relationship between the quality and macroscopic deformation of the sample and the number of freeze-thaw cycles

2.2 對土體力學性質的影響

2.2.1 應力-應變關系

不同凍融次數、不同圍壓下的應力-應變曲線如圖4 所示，其中F-T 表示凍融次數。從圖4 可以發現，圍壓100 kPa 時應力-應變曲線呈現出應變軟化型，經過多次凍融循環，有向應變硬化型轉變的趨勢，剪切破壞時試樣產生破裂面；圍壓200、300 kPa 時表現出應變硬化型，破壞形式也由脆性轉變為塑性，試樣發生剪脹破壞。通過不同凍融次數的應力-應變曲線可以看出，前期凍融循環作用對土體強度的影響十分敏感，尤其是第1 次和第3 次凍融結束后，不同圍壓下的應力-應變曲線均明顯低于其他凍融次數，隨著凍融次數增加，強度有所增大并趨于穩定。

圖4 不同凍融次數、不同圍壓的應力-應變曲線Fig.4 Stress-strain curves of different freeze-thaw times and different confining pressures

2.2.2 破壞強度

凍融次數與破壞強度的關系如圖5 所示，可以看出，經過3 次的凍融循環后土體結構產生損傷，破壞強度急劇降低。根據圖3 試樣的凍縮變化率發現，凍融1～3 次的高度凍縮率曲線的斜率是3～5 次的1.97 倍，凍融作用的影響明顯減弱。凍縮速率的快速減小代表內部結構將不會再進行劇烈調整，隨著土體的凍縮，固結度提高，結構強度逐漸恢復；另外，第5 次凍融后的試樣抗剪強度發生大幅回升，說明第4 次凍融后土體抗剪強度應增加而非降低。因此，凍融循環的劣化作用發生在前3 次。雖然多次凍融后強度趨于穩定，但仍低于未凍融時的強度，說明凍融作用對土體造成的損傷具有不可恢復性。

圖5 不同凍融次數下破壞強度變化Fig.5 Change of failure strength under different freeze-thaw times

2.2.3 抗剪強度指標

內摩擦角和黏聚力在不同凍融次數下的變化規律如圖6 所示。內摩擦角經歷1 次凍融之后便達到平衡；黏聚力與強度變化相似，先減小后增加并達到穩定。經分析，初次凍融后，由于水分的相變和遷移，體積膨脹擠壓鄰近土顆粒，破壞土顆粒間的膠結，降低咬合摩擦力，引起內摩擦角和黏聚力均減小[22]。試樣的凍縮作用增大固結度，顆粒團聚體解體，小顆粒增多，咬合程度提高，所以內摩擦角受后續的凍融影響較小[3]。土體結構凍融過程中的水分運移、溶蝕、沖刷黏性礦物和產生的冷生結構等削弱了顆粒間的膠結力，導致黏聚力繼續降低，隨著凍融次數增加，水分運移通道趨于穩定，新的土體結構逐漸形成，因此黏聚力也不再發生較大變化。

圖6 不同凍融次數下抗剪強度參數變化Fig.6 The change of shear strength parameters under different freeze-thaw times

3 理論分析與離散元數值模擬

依據室內試驗結果，本文將通過離散元數值模擬研究前3 次凍融循環作用引起強度損傷的細觀影響機理和作用過程。

3.1 離散元單元間接觸模型

顆粒離散元法中，巖土體是由一系列遵循牛頓運動定律的離散單元堆積而成，單元間依靠可斷裂的彈簧表示法向力(Fn)、切向力(Fs)和摩擦力(Ff)等效真實世界中巖石或土壤中礦物顆粒之間的膠結作用。接觸模型參數中Kn為法向剛度，Ks為切向剛度，Xb為斷裂位移，Xn為法向相對位移，Xs為切向相對位移，Fs0和μp分別為切向初始抗剪力和摩擦因數，如圖7 所示。

圖7 離散元顆粒接觸關系Fig.7 Contact of discrete element

單元間法向力與法向剛度的關系可用下式表示：

相鄰單元連接時相互之間存在彈簧的拉壓作用，若位移Xn超過斷裂位移Xb時彈簧斷裂，單元間將不會產生拉力，但是發生擠壓時仍然可以產生壓力作用。

單元間切向力與切向剛度的關系為：

根據摩爾庫倫準則，單元間連接完整時最大剪切力Fsmax表示為：

當剪切力超過Fsmax時，彈簧發生斷裂，初始抗剪力Fs0消失，發生滑動，產生的摩擦力Ff為-μpFn，與真實顆粒相比，膠結斷裂意味著法向和切向之間的連接均會斷開。

由于在物理世界中巖土體的顆粒并不是完全彈性體，內部存在水分、節理或微裂隙等因素阻礙顆粒的運動，動能會逐漸轉化為熱能而衰減。在數值模擬中，單元上的黏滯阻力(Fv)表征了真實世界中各種因素對動能的消耗，MatDEM 通過引入該力避免了動能在系統中的積累[23]。黏滯阻力由下式給定：

式中：ψ為阻尼系數，與模型的體積、顆粒直徑、顆粒質量和顆粒剛度系數有關；μ為顆粒的速度。

3.2 離散元多場耦合理論

在離散元數值模擬中，模型中的顆粒代表由一定體積黏土顆粒和水組成的集合體，其中水是虛設值，并不真實存在。顆粒除了具有一定的力學性質外，還可以通過虛設一個溫度值，賦給顆粒后讓其具有溫度屬性。在本研究中，真實土體的條件是含水率為17.6%，溫度為20℃，將這2 種條件賦到顆粒上便可以通過建立半徑與溫度、含水率和應力之間的關系，實現凍融循環中溫度場、水分場和應力場的耦合模擬。假定顆粒的含水率表示土顆粒及其周圍區域的平均含水率，溫度和水分通過顆粒間的接觸進行傳遞。

由于水凍結時會釋放潛熱出現過冷現象[24]，模擬時假設將0℃至過冷溫度的區域設為凍結緣，且水分僅在該區域內運移[25]。

3.2.1 溫度場和水分場控制方程

在數值模擬中，首先進行溫度傳遞，單元顆粒經過一個時間步的溫度變化量?T，通過利用溫度場控制方程[26]轉化為下式進行計算：

式中：C為土的體積熱容，kJ/(m3·K)；λ為土與水的平均導熱系數，W/(m·K)；η為顆粒間溫度傳遞影響系數，與時間步數和循環次數有關；L為冰-水相變潛熱，一般取值為334.56 kJ/kg；T為單元顆粒初始溫度，℃；Tj為周圍接觸到的單元顆粒溫度，j=1,2,···,N；ρi為冰密度，取值0.917×103kg/m3；Vi為溫度低于0℃的顆粒單元體積；V為組成模型的單元顆粒總體積。

隨著溫度的降低，顆粒間逐漸出現水分運移現象。每個時間步的水分運移量?ω參考水分場控制方程[27]，轉換為如下公式進行計算，升溫時也按照此公式進行水分回遷模擬。

式中：D為水分擴散系數；ρw為水密度，取 值1×103kg/m3；ζ為顆粒溫度差與水分運移量的系數，通過時間步數和循環次數等模擬訓練獲得；ωk為與周圍接觸到的單元顆粒的含水率，k=1,2,···,b；Tf為過冷溫度，取-1.5℃[28]。

水分擴散系數D計算公式為：

式中：a、m、l為水土特征曲線擬合參數，引用基于室內試驗數據擬合的V-G 模型和凍土物理學[29]；θi為冰的體積含量；S為相對飽和度；θw為液態水的體積含量；ws為飽和含水率；wc為殘余含水率；Kb為飽和滲透系數。

3.2.2 單元顆粒半徑變化及對強度的影響

在凍融過程中，由于溫度的變化，土骨架凍縮、回彈和水分的運移及相變均會對土體體積產生影響，為了在數值模擬中體現這一過程，本文將上述變化細化到單元顆粒的變形上，把單元顆粒的變形原因視為土顆粒和水分變化兩部分。假設顆粒的半徑按照土顆粒和水分之間的質量比值進行劃分，即水分所占半徑的大小為含水率的值ω，土顆粒占半徑的(1-ω)。對于土顆粒，計算時采用土顆粒所占半徑的比例乘以土顆粒收縮系數，模擬土顆粒受溫度影響引起的顆粒半徑變化。對于水分，考慮了水分運移引起的含水率增加或減少以及水分凍結和融化導致的顆粒半徑脹縮變化。為了展示變化過程，給出了顆粒在含水率增加時的凍融過程圖如圖8 所示。

圖8 含水率增加時顆粒變化Fig.8 Particle changes when the moisture content increases

因此，當顆粒溫度低于過冷溫度或高于0℃時，顆粒半徑R的變化采用如下公式計算：

式中：R0為初始單元顆粒半徑；R1為融化時單元顆粒半徑；ω0為初始含水率；ωd、ωr分別為凍結完成和融化完成后的含水率；γ為土顆粒收縮相關系數，根據試樣凍縮量取值。溫度低于過冷溫度視為凍結完成，高于0℃時視為融化完成，顆粒半徑不變。

依據應力場基本方程中應力σ與彈性模量E、泊松比v和體積變形有關[30]，模型中顆粒間的應力場計算公式如下：

式中：RB為凍結前顆粒半徑；RL為凍結或融化后的顆粒半徑。

由于顆粒半徑R和泊松比v與顆粒間的法向剛度和切向剛度之間存在關系，從Liu Chun 等[31]給出的參數轉換公式可知：

通過式(1)、式(2)可知法向力和切向力與法向剛度和切向剛度相關，所以模型中顆粒的正應力 σn和切應力σs的計算公式為：

式中：FnB、FsB與FnL、FsL分別為顆粒半徑為RB和RL時的法向力和切向力；xnB、xnL為顆粒的法向相對位移；xsB、xsL為切向相對位移。

根據室內試驗結果，試樣強度劣化發生在前3 次的凍融循環，將彈性模量、泊松比與凍融次數結合，擬合公式如下：

式中：B為凍融次數。

在數值模擬單次迭代中，首先利用式(5)實現顆粒間的熱量傳遞；利用式(6)實現水分在溫度的影響下發生運移；然后篩選溫度低于凍結溫度或高于0℃的顆粒，利用式(9)或式(10)分別實現凍結和融化時溫度、水分運移及相變對顆粒半徑產生的影響；顆粒半徑的增大和縮小引起顆粒間的擠壓和移動，改變了顆粒間的力學性質，通過式(14)或式(15)計算當前狀態下的應力；最后平衡模型。由于顆粒位置發生重排列，顆粒間的接觸狀態被改變，進而影響下一次迭代中溫度傳遞、水分運移和應力變化。經過以上步驟的循環迭代實現凍融條件下土體內部的水熱力三場耦合。

3.3 離散元模型的建立

3.3.1 堆積模型的建立

本模擬試驗基于南京大學自主研發的三維離散元模擬軟件MatDEM[32]。通過該軟件首先建立一個內徑3.91 cm 的圓管，在圓管內生成平均半徑0.1 cm 的隨機顆粒，然后對單元顆粒施加隨機速度進行壓密沉積，建立直徑0.039 1 m、高0.08 m 的模型如圖9 所示。由圓管、頂板及底板構成的恒溫結構如圖10 所示，溫度設定為20℃或-20℃，不參與水分運移和不對模型產生約束，僅作為溫度邊界。模型共由24 466 個活動單元組成，其中試樣模型由17 930 個單元顆粒組成。

圖9 試樣模型Fig.9 Sample model

圖10 恒溫結構Fig.10 Thermostatic structure

3.3.2 模型的材料設置

采用離散元法建模時，為了獲得與真實土體宏觀力學性質相似的堆積模型，在MatDEM 軟件中的材料自動訓練模塊，利用室內試驗獲取的土體力學參數如：彈性模量(E)、泊松比(v)、抗壓強度(Cu)、拉伸強度(Tu)和摩擦因數(μp)，進行連續5 次的訓練優化，獲得建立模型所需的宏觀力學參數，然后根據Liu Chun 等[31]給出的宏微觀參數轉換公式，計算反映模型顆粒間性質的微觀力學參數，如：法向剛度(Kn)、剪切剛度(Ks)、初始抗剪力(Fs0)、斷裂位移(Xb)和摩擦因數(μp)。完成材料訓練后，將獲得的力學參數賦給堆積的試樣模型，進行單軸無側限壓縮模擬試驗，并與室內無側限抗壓試驗結果對比，如圖11 所示。室內試驗所獲的土體力學參數和通過對模型材料自動訓練獲得的宏微觀力學參數見表1。

表1 土體室內試驗力學參數和模型材料的宏微觀力學參數Table 1 Macroscopic parameters of soil test maro-micro parameters of model meterials

圖11 室內與數值模擬的無側限抗壓強度Fig.11 Indoor and numerical simulation of unconfined compressive strength

數值模擬中溫度控制方程和水分傳導方程的相關參數見表2，其中η和 ζ經過多次訓練分別取0.1 和0.05，水和冰的導熱系數λw和λi均取自規范值[29]，式(7)、式(8)中ws、wc、Kb及式(5)中C、λ根據已有的室內土工試驗成果進行取值。

表2 分析計算的相關參數Table 2 Macroscopic and microscopic mechanical parameters of soil

4 數值模擬結果分析

4.1 溫度場演變規律

試樣的凍融過程如圖12 所示。根據模型中心區域的3 個溫度節點：0℃、過冷溫度和-20℃，將凍結過程劃分為溫度快速下降、緩慢相變過程、繼續降溫、溫度穩定4 個階段。第一階段，如圖12a-圖12e 所示，在50 min 時由于試樣表層與恒溫結構之間溫差大，顆粒快速凍結，形成明顯凍結區；凍結區內側分布著溫度位于0℃附近的一層較薄區域，形成凍結緣區[33]；距離冷端較遠的顆粒溫度未明顯下降，形成未凍土區。至 170 min 時，因凍結面與未凍區溫差變小，所以凍結區擴大較慢，但是內部溫度降低較快，模型中心區域溫度基本全部降至0℃左右，結合圖12c 和圖12d 可以發現，凍結鋒面推進較多，凍結緣變得越來越厚。凍結進入第二階段，如圖12f 和圖12g 所示，此時凍結緣內冰水共存，水相變時釋放潛熱，減緩降溫速度。直到250 min 后模型核心溫度全部低于過冷溫度，可視為模型已完全凍結。第三階段開始，此時模型中溫度的傳遞速度有所提高，顆粒溫度繼續降低，到400 min后(圖12h)，降溫過程結束，進入第四階段，顆粒與恒溫結構沒有溫差，達到溫度穩定狀態。

圖12 凍融過程中的溫度傳遞Fig.12 Temperature transfer during freeze-thaw

融化過程與凍結過程類似，同樣經過4 個階段。將圖12i-圖12k 所處時間點的溫度變化與凍結過程相比，區別在于融化過程中0℃以下的變溫速率大于0℃以上的變溫速率，直至全部融化。第2、第3 次凍融循環中的溫度傳遞過程與第1 次基本一致。

4.2 水分場演變規律

依據4.1 節中對凍融過程的階段劃分，可以看出水分運移范圍跟隨凍結區向模型的內部逐漸轉移。在凍結過程的第一階段中，如圖13a-圖13e 所示，由于冷端溫度與未凍區之間溫差大，土中水分成冰較快，水分運移量很少。第二階段中，如圖13f-圖13g 所示，凍結區發展緩慢，但凍結緣增厚，土中水分有較長的時間進行遷移，所以圖13e-圖13g 中含水率升高的顆粒明顯增多。第三階段中，模型已完全凍結，含水率基本不變，如圖13h 所示。以上模擬水分運移過程說明，溫度梯度的存在導致引起水分運移的土水勢出現，水分從高溫向低溫處遷移[34]；溫差越大，土體凍結越快，水分以原位凍結為主；溫度梯度較小時，凍結較慢，水分運移增加。凍結區邊緣顆粒的體積會通過周圍顆粒的水分補給逐漸增大，擠壓周圍顆粒，說明水分運移會引起水分重分布，局部含水率升高，水分凍結成冰后導致土體發生不均勻變形，改變顆粒排列方式，破壞內部結構。

圖13 凍融過程中水分運移Fig.13 Moisture migration during freeze-thaw

融化過程中，在溫度梯度和土水勢的作用下，土體水分發生回遷現象，高含水率顆粒中的水分向周圍低含水率的顆粒轉移，如圖13i-圖13l 所示。升溫結束后，發現水分并未完全回遷，導致含水率分布不均，說明顆粒凍結后位置被動調整，與周圍顆粒的連接狀態發生改變，所以不能按照水分運移的路徑回遷。在第2、第3 次凍融循環中，水分的遷移會繼續引起土體結構的變化，含水率升高或降低的顆粒逐漸增加，進而對溫度傳遞、水分運移的路徑和結構強度產生影響。

4.3 單元顆粒半徑變化及宏觀變形

圖14 為17 930 個模型顆粒半徑在首次凍結后和融化后的大小分布圖，可以發現顆粒半徑主要集中在0.000 75～0.001 10 m，凍融后的平均半徑低于0.001 m，發生細觀角度上土顆粒和水分集合體的凍縮現象。受水分運移的影響，凍融前后相同半徑區間的單元顆粒數均有差別，融化后的分布曲線偏低，顆粒大小更為分散。隨著凍融次數增加，含水率升高和降低的顆粒增多，最大半徑和最小半徑的顆粒也隨之增長。顆粒半徑的變化引起周圍顆粒發生位移，原本的接觸狀態隨之改變，形成位移場，如圖15 中顆粒在X方向的位移。通過圖15a、圖15b 對比發現，前后位移場存在一定的差異，說明模型在凍融作用下內部結構受顆粒的脹縮影響已經發生變化。

圖14 凍結、融化后顆粒半徑分布Fig.14 Particle radius distribution after freezing and thawing

圖15 凍結、融化后X 方向位移場Fig.15 X-directional displacement field after freezing and thawing

通過對模型凍結和融化后的直徑和高度變化數據的提取，凍結后含水率升高的顆粒半徑明顯增加，擠壓周圍顆粒，引起模型凍脹；融化后，顆粒半徑縮小，在顆粒間膠結力、黏聚力等的作用下向內側移動，變形逐漸恢復(圖16)。整體上可以看出模型在凍融過程中體積先凍脹而后融沉的現象。由于顆粒代表一定體積的土水集合體，顆粒直徑占高度的1/80，直徑的1/39，所以宏觀變形上較難達到一致，但是模擬的最高誤差未超過0.172%，可以有效說明凍融作用對土體變形的影響。

圖16 凍結、融化后高度和直徑變化率Fig.16 Rate of change in height and diameter after freezing and thawing

4.4 溫度和水分綜合作用對應力和強度的影響

凍融作用引起的內部應力變化如圖17 所示。凍結前顆粒間的應力受重力影響由上至下逐漸增加，但是變化并不明顯。經過凍融作用，在溫度、水分運移及相變的影響下引起部分顆粒的脹縮讓模型內產生局部的變形，引起顆粒重新排列，導致顆粒斷開連接或與其他顆粒建立新連接，改變了顆粒間的受力狀態，破壞了土體的內部結構，最終造成應力明顯降低。

圖17 凍結前、融化后X 方向應力場Fig.17 X-direction stress field before freezing and after melting

將經過0、1 和3 次凍融過程的模型進行圍壓為100 kPa 的單軸壓縮模擬試驗，與室內相同條件的試驗結果進行對比，結果如圖18 所示，模型破壞強度大小與室內試驗結果基本相同。在凍融循環模擬中，顆粒的脹縮、冰-水相變、水分聚集與回遷現象的發生，改變了土體原有的穩定結構。隨著凍融次數的增加，顆粒之間的孔隙、位置、連接狀態和應力大小不斷變化。其次，凍結時產生的水分運移令部分顆粒持續獲得水分而明顯增大，與土體中局部孔隙得到水的補充后，凍結成冰發生膨脹現象一致，改變了顆粒間的相對位置及接觸狀態，最終造成難以恢復的結構損傷，導致土體強度降低。

圖18 不同凍融次數的強度對比Fig.18 Strength comparison of different freeze-thaw crcles

利用MatDEM 軟件中的材料訓練模塊和能對單元顆粒添加虛設值的特點，獲得與真實土體溫度和含水率相同、力學性質相似的模型，并通過溫度與含水率變化影響顆粒半徑，引起的模型強度和變形與室內試驗結果基本一致，說明數值模擬的結果是可靠的，且具有實際參考價值。

5 結論

a.土體在凍融循環作用下發生凍縮，體積減小；抗剪強度及其參數的劣化主要受前3 次凍融作用的影響，第3 次劣化程度最大，隨著凍融次數的增加，抗剪強度及其參數回升并趨于穩定。

b.建立的水、熱、力和半徑變化的計算方程，適用于非飽和黏性土；可較好地模擬溫度、水分和應力的動態變化規律，獲得的試樣變形和破壞強度與室內試驗結果吻合度也較高；可用于描述凍結鋒面的移動軌跡和含水率的變化過程。

c.溫度傳遞過程中可劃分為溫度快速下降、緩慢相變過程、繼續降溫、溫度穩定4 個階段；水分運移主要發生在溫度傳遞的前2 個階段，緩慢相變過程中的水分運移量居多。

d.凍融作用導致顆粒整體縮小且分散性增大；土體主要由于土顆粒的脹縮、冰-水相變、水分運移和冷生結構等共同作用引起位移場和應力場改變，產生不可逆的結構損傷而造成抗剪強度劣化。