青藏鐵路填土路基熱-力耦合離散元研究

石梁宏，李雙洋，尹 楠

（1.中國科學院西北生態環境資源研究院凍土工程國家重點實驗室，甘肅蘭州730000；2.江蘇蘇邑設計集團有限公司，江蘇南京210041）

0 引言

凍土作為一種特殊土類，與其他巖土體最主要的區別是凍土中有冰的存在，而冰的含量又與外界環境溫度息息相關［1-4］，溫度的變化會使得凍土的物理性質、力學性質、電學性質、滲透性質以及其他性質呈現動態變化的特性，這些變化關系著凍土強度和變形等力學以及熱學等性能的變化，而凍土又在具有海拔高、氣壓低、輻射強、穩定性差和局地差異強等特點的青藏高原地區廣泛分布［5-7］，這勢必會使得青藏高原凍土路基工程的熱、力穩定性面臨更嚴峻的挑戰［8-9］。

為此，可采用不同的路基結構形式或材料，通過調節傳導、對流和輻射的傳熱方式來調控路基溫度場，并且在凍土路基溫度場的長期熱穩定性方面有比較多的研究成果［10-12］，較好地解決了青藏鐵路建設中可能發生的主要路基病害［13］。而在凍土路基力學變形方面，有限元法、有限差分法、有限體積法等數值計算方法被廣泛的應用于凍土路基的熱、力學穩定性研究中。其中，王鐵行等［14］在土力學彈塑性理論和流變理論的基礎上，建立了凍土路基應力及變形的二維數值模型，模擬水、熱及土性的動態變化，重點考慮了土性變化情況下土的流變變形和瞬時變形特征；汪雙杰等［15］針對青藏公路路基下發育多年凍土融化盤的實際情況，應用ABAQUS有限元分析軟件，對凍土路基從修筑到開放交通過程中的路基路面位移及應力進行了分析；許健等［16］以青藏鐵路北麓河試驗段為計算模型，計算凍土路基的溫度場和溫度影響下的應力、變形場，并且與實測溫度和路基變形數據進行了對比分析；毛雪松等［17］建立了凍土路基變形場的二維數值計算模型，并應用有限元的方法求解路基土體凍結時變形場的分布規律，并進一步分析了凍土路基破壞的機理；李寧等［18］基于正凍土三場耦合的理論框架以及所開發的全面考慮正凍土骨架、冰、水三相介質水、熱、力與變形的真正耦合作用的分析系統3G2001，對214國道花石峽試驗進行數值研究，并與路基實測的地溫變化和路基路面變形進行了對比分析驗證；李雙洋［19］、Li等［20］以凍土區青藏鐵路路基某斷面為例，運用傳熱學、凍土物理學和凍土流變學的基本理論，考慮水分場和溫度場的相互作用，以及溫度變化對路基土體力學特性的影響，并引入凍土的蠕變方程，建立了凍土路基的熱、力學（蠕變）穩定性分析模型，對該段路基運營若干年后的熱、力學狀況進行了分析和預測。

但是，以上研究都是基于連續介質力學的分析方法，不能充分考慮顆粒之間的導熱與接觸的相互作用，在分析凍土路基熱-力穩定性的破壞變形上也存在瓶頸。為此，可采用顆粒離散單元法進行分析研究，該方法克服了傳統有限元法難以模擬凍土路基的大變形與破壞等情況。此外，顆粒離散單元法可模擬顆粒間的導熱和點點接觸粘結，并從微細觀水平反映路基的宏觀熱、力變化狀況。

基于目前凍土路基熱-力穩定性研究中存在的不足，本文選取青藏高原五道梁地區某凍土路基斷面為例，首先建立隨凍土路基溫度場變化的熱-力耦合離散元計算模型，在驗證數值模型正確可靠后，對路基的熱、力狀況進行預測和分析，并從微細觀水平分析凍土路基的宏觀變形。本項研究不僅可填補凍土路基離散元研究的空白，而且也拓展了離散單元法在凍土路基工程領域的研究，在理論與實踐上都具有很重要的意義，可更好地為寒區工程提供服務。

1 數學模型及方程

1.1 數學傳熱離散元模型

將物質離散為熱源和熱通道組成的網絡連接體系，從單個熱源出發，建立離散化的熱量傳輸方程。單個顆粒的控制體積為V（通常假定所有顆粒的控制體積等于總物質體積），單位體積中熱流出量可表示為qi的散度，其平均值定義為

運用高斯散度定理，用面積積分代替體積積分，得

式中：ni為面S的單位法向向量（向外）。

假定離散化的網絡連接體系中形成N個熱通道，限定熱量在N個熱通道中流動，則面積積分可用下式表示。

式中：上標（p）為連接通道p的相關變量；Q（p）為在連接通道p中從熱源流出的熱量，qiΔS（p）=Q（p）ni（p）。

將式（3）和式（2）代入式（1），得到

假定溫度對應變變化的影響很小，可以忽略，這個假定應用于涉及到固體和流體的準靜態力學問題，則這種連續介質的熱傳導控制方程為

式中：qv為體積熱源強度（W·m-3）；ρ為質量密度（kg·m-3）；Cv為 質 量 熱 容（J·kg-1·℃-1）；T為 溫度（℃）。

將式（4）代入式（5），得到熱傳導方程。

1.2 顆粒粘結模型

凍土的力學性質主要取決于土中膠結冰的性質［21］，所以在凍土領域運用離散單元法進行仿真研究，需要充分考慮到凍土中冰膠結的重要作用。在離散元法中，可以采用粘結發生在接觸顆粒間有限范圍內的模型來考慮凍土中冰的膠結特性，這種粘結模型能更好地描述冰的膠結性質，如圖1所示的膠結結構［22-23］。取A和B兩個顆粒之間的膠結特性進行研究，如圖2所示。

圖1 凍土中膠結結構示意圖［23］Fig.1 Sketch map of the cementation structure in frozen soils［23］

在三維情況下，顆粒間接觸的相對運動在膠結連接處會產生力和力矩，作用于相互連接的顆粒上，將顆粒A和B相互接觸作用的區域等效為一個圓盤，模型中的接觸力和力矩可分解為法向（Fn、Mn）和切向（Fs、Ms）矢量的形式［23］。在平行粘結模型中，法向抗拉強度和切向抗剪強度分別為-σc、-τc，等效圓盤的半徑為R，則作用在等效圓盤上的最大拉應力和剪應力的計算方法為

圖2 顆粒間的接觸模型示意圖Fig.2 Sketch map of the contact model between particles

式中：A=πR2，為等效圓盤的面積；I為等效圓盤截面的慣性矩；J為等效圓盤橫截面通過接觸點的慣性矩。如果-σmax≥-σ或者-τmax≥-τc，則發生粘結破壞。

2 數值計算模型

根據青藏鐵路沿線鉆探地質資料［24］，建立路基實體模型，圖3所示。模型中填土路面寬度為7.2 m，填土高度為5 m，路基邊坡坡度為1∶1.5。路基模型在坡腳處水平向外各延伸30 m，路基下部土層分為三層：2 m厚的砂土層，6 m厚的粉質黏土層，22 m厚的弱風化巖層。

圖3 路基實體模型Fig.3 Physical model of the embankment

圖4 為凍土路基離散元數值計算模型圖，數值模型的尺寸和各層土質與實體模型一致，外側采用墻體限制。在實際計算中，考慮到路基幾何及邊界條件的對稱性，故取一半路基作為研究對象。

圖4 路基離散元數值計算模型Fig.4 DEM（discrete element method）model of the embankment

考慮全球氣候變暖的影響，取青藏高原未來50年平均氣溫上升2.6℃［25］，由于路基上表面溫度變化不僅與環境空氣有關，而且受到太陽輻射等復雜因素的綜合作用，故可根據附面層理論［26］，對離散元計算模型的熱邊界條件進行如下的設定：

天然地表AB和EF的溫度按下式變化［27］。

式中：th為時間變量。可通過調整α0來改變初始時刻對應的日期，α0=0對應的初始日期為7月15日。

路基斜坡BC和DE的溫度按如下規律變化。

路基頂面CD的溫度變化規律為

邊界ALKJ和FGHI視為絕熱邊界，即

底部邊界JI為熱流邊界，熱流密度為0.06 W·m-2。

離散元的熱學計算需用到5個熱學微觀參數，即密度ρ、線性熱膨脹系數α、質量比熱容Cv、導熱系數λ（或熱阻η）及質量相變潛熱Ls，表1中列出了凍土路基離散元數值計算中的微觀參數。

表1 路基各介質中的離散元熱學微觀參數Table 1 DEM thermal microcosmic parameters in the embankment

在計算中，對于含水介質中相變潛熱問題采用顯熱容法進行處理［28］，假設模型中含水介質相變發生在溫度區間(Tm±ΔT)。當建立等效質量熱容時，應考慮溫度間隔ΔT的影響，同時假設介質在已凍、未凍時的質量比熱容熱容Cf和Cu及導熱系數λf和λu不取決于溫度，此時簡化構造出質量比熱容和導熱系數的表達式［29-30］為

式中：Ls為含水介質單位質量相變潛熱；T為一定范圍內顆粒溫度的平均值。

計算中，顆粒接觸剛度與路基溫度之間的數學關系是基于凍土室內靜三軸試驗，考慮溫度為路基熱-力穩定性的主要影響因素，結合凍土力學、傳熱學、離散元理論以及凍土相關研究成果等理論而建立起來的［20-21，28，31-33］。路基不同土層顆粒接觸剛度與溫度之間的數學關系為：

幾種不同介質土層的顆粒接觸剛度與溫度之間的數學關系，其變化規律具有一致性。即在溫度為正時，顆粒剛度為常數；溫度為負時，顆粒接觸剛度隨溫度的降低而線性增大。

3 結果與分析

3.1 數值模型校核

圖5 路基天然孔地溫變化曲線Fig.5 Ground temperature curves of origin borehole beside the embankment

為了驗證離散元數值模型的準確性，將數值模型的計算結果與已有數據進行對比。從圖5可以看出，離散元計算的路基天然孔地溫曲線與實測結果吻合較好，離散元數值計算的凍土上限（0℃等溫線）與實測的凍土上限分別為-1.96 m和-1.89 m，相差僅為0.07 m，在路基天然孔地溫曲線整體變化趨勢上，兩者無較大的差別。因此，凍土路基離散元數值模型能真實地計算凍土路基的熱變化，模型的準確性較高。

圖6 為路基建成后路肩孔地溫曲線的比較，可以看出，離散元數值計算和測量的地溫曲線總體比較吻合，但是在0～1 m范圍內，可以看出離散元模擬曲線有一個明顯的“突出”，可能的原因是：熱量在粒間通道中傳導，傳導的過程需要時間，而路基剛建成時，路基上部結構與天然地表的溫度存在差異，天然地表顆粒的溫度會影響地表上部一定區域顆粒間的熱量傳輸，導致在天然地表上部1 m范圍內產生明顯的溫度變化，這與真實土顆粒之間在溫度變化下的相互作用很相似，也說明了離散元法在顆粒熱交換相互作用特性方面的模擬更具真實性。總體來看，所建立的凍土路基熱-力耦合離散元數值計算模型能較合理地體現路基中顆粒之間熱-力作用的變化特征，能較真實地反映出上部路基及下部基礎內顆粒體的熱交換作用與溫度分布狀況。

圖6 路基路肩孔地溫變化曲線（圖中天然地表面為縱坐標的原點）Fig.6 Ground temperature curves of shoulder borehole of the embankment（Natural surface is set as the origin of vertical coordinate）

3.2 路基建成第2年的溫度場

圖7 選取了路基建成后第2年四個典型季節的溫度分布圖進行分析。從路基等溫線變化情況來看，四個季節都形成-0.8℃溫度線。從路基融化深度來看，在7月15日，上部路基內融化深度約為2.15 m，而在10月15日，上部路基內融化深度約為3.20 m。

進一步分析來看，在路基頂部、邊坡、天然地表周期變化的溫度邊界以及底部恒定的熱流邊界作用下，多年凍土區路基中顆粒之間熱傳導方式形成的溫度場，在下部基礎內擾動相對較小，相對較穩定；但是在周期溫度變化的邊界下，路基內溫度場擾動性較強，尤其在凍土上限區域，顆粒之間熱傳導形成的溫度場呈現“紊流”狀或“島”狀的分布較明顯，等溫線波動劇烈。這些變化特征也從微觀層面說明，在實際凍土路基中土顆粒之間相互作用和熱交換過程的復雜性，從而會形成較紊亂的溫度場分布，在寒區從事路基工程活動對凍土存在較大的擾動。

3.3 路基建成第20年的溫度場

圖8 中，為凍土路基修筑完成后第20年四個典型季節的溫度場，可見，在這四個典型季節中，-0.8℃溫度線已退化，-0.5℃溫度線退化也很嚴重，而且上部路基內融化深度有一定下移，其中，在7月15日融化深度為2.40 m，下移0.30 m；而在10月15日融化深度為3.80 m，下移0.60 m。從圖8（d）可知，天然地表下凍土上限為-2.70 m，而第2年10月15日的凍土上限為-2.25 m，下降了0.45 m，路基存在凍土退化的問題。這些離散元研究結果與路基已有的研究結果一致，可以說明，運用離散單元法對多年凍土區路基的熱-力耦合數值仿真較合理，也驗證了隨路基溫度場變化的熱-力耦合離散元計算模型的準確性，并且顆粒集合體之間的熱交換作用特性也較真實體現凍土路基土顆粒的實際變化情況。

進一步分析來看，在路基修筑完成后第20年的溫度場仍存在擾動變化特征，其變化規律與第2年的情況類似，但是從路基修筑完成后第20年4月15日溫度分布來看，靠近天然地表區域的顆粒體之間熱交換作用比較劇烈，形成的溫度場變化較復雜，也體現了實際凍土路基長期熱穩定性變化過程的復雜性。

3.4 路基變形分析

圖7 普通路基建成后第2年的地溫分布Fig.7 Distributions of ground temperature of the embankment in four typical seasons in the 2nd service year

圖8 路基建成后第20年的地溫分布Fig.8 Distributions of ground temperature of the embankment in four typical seasons in the 20th service year

圖9 路基建成后第20年10月15日的位移分布Fig.9 Distributions of displacments of the embankment on October 15 in the 20th service year

圖9 為多年凍土區路基建成后第20年位移等值線圖，選取10月15日的位移場進行分析，其余三個季節的變化規律類似。整體來看，在熱-力變化作用下，路基顆粒集合體之間相互接觸、摩擦、熱交換等作用比較充分，體現在路基位移等值線上，不是平滑的曲線，而是呈現“鋸齒狀”波動變化的曲線，這也較真實體現出實際凍土路基復雜變化的特征。從圖9（a）可以看出，路基坡腳處相當于水平位移正負變化特征的分隔點，坡腳以左，路基水平位移為負，表明上部路基對下部基礎向下的壓縮作用較強；坡腳以右，路基水平位移為正，表明下部基礎內被壓縮的顆粒體會向右繼續擠壓產生變形，基礎部分水平向外延伸越遠，則水平壓縮變形作用對顆粒體的影響作用就越弱。從圖9（b）可以看出，一方面，在重力作用下，路基顆粒集合體會發生沉降變形，另一方面，路基修筑完成后，上部路基結構的重量會向下壓縮基礎內的顆粒體，在坡腳以右，會繼續擠壓下部基礎內的顆粒體，因此，從路基中心位置向右，路基沉降變形量會不斷減小，并且，在路基水平向外延伸較遠處與下部基礎較深處，沉降變形作用會減弱。

圖10 為多年凍土區路基建成后第20年10月15日的接觸力鏈分布，其余三個季節變化規律類似。路基中接觸力鏈的粗細可表征路基顆粒間接觸力的大小，也能從微細觀水平反映出路基的宏觀變形特征。可見，在重力作用下，路基形成強力鏈與弱力鏈共同發展的體系。上部路基內的顆粒體會向下與向右運動，擠壓基礎內的顆粒體，在路基坡腳以左，路基中的柱狀力鏈呈現向左的傾斜特征，在路基中心處，柱狀力鏈呈垂直分布，從路基中心至坡腳處，柱狀力鏈向左的傾斜角度逐漸減小，表明了這部分路基中顆粒間擠壓作用向右逐漸減弱；在路基坡腳以右，柱狀力鏈向左的傾斜變化特征逐漸消失，形成柱狀力鏈之間交叉分布的特征，體現了在豎向重力作用下路基顆粒間接觸力的變化特點。而路基中網狀弱力鏈分布在柱狀強力鏈周圍，并與強力鏈處處連接，對強力鏈的穩定起輔助作用，從而使路基發生整體變形。還可知，路基中力鏈的粗細分布不均勻，各向異性很明顯，坡腳以左路基各向異性的程度較坡腳以右路基中的大，表明坡腳以左路基中顆粒間相互作用較強烈。

圖10 路基建成后第20年10月15日的接觸力鏈分布Fig.10 Distributions of contact force chains of the embankment on October 15 in the 20th service year

4 結論

本文選取青藏高原五道梁地區某路基斷面，在考慮凍土路基中復雜的工程狀況基礎上，采用離散單元法，建立了凍土路基溫度場變化的熱-力耦合離散元計算模型，預測和分析了路基的熱、力分布狀況，并從微細觀水平闡釋了凍土路基宏觀變化的復雜性。結論如下：

（1）對比有限單元法和離散單元法可發現，離散單元法突破了連續介質力學中變形條件的限制，可模擬顆粒之間的相互運動，因此能夠較真實反映出路基的熱、力變化特征，可應用于寒區路基工程的長期穩定性研究。

（2）多年凍土區路基離散元仿真研究表明，隨著凍土路基運營時間的增加，凍土上限會發生下移現象，因此路基存在凍土退化的問題。

（3）凍土路基中顆粒之間熱交換過程較復雜，在凍土上限（即0℃等溫線）區域，溫度場呈現“紊流”狀或“島”狀分布，擾動較劇烈，能較真實體現實際凍土路基中顆粒之間的相互作用特性和熱交換情況。

（4）凍土路基在坡腳以左，顆粒體的相互作用較劇烈，路基（沉降）變形主要發生在路基主體結構及下部基礎一定范圍內。

本文的研究具有一定的探索性，所選取斷面的熱邊界條件具有對稱性，并且未考慮陰陽坡效應和水分的作用，因此還需要對該系統進一步的研究。作為初步研究，可為今后凍土工程的穩定性研究和寒區工程的建設提供科學依據。