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

石梁宏，李雙洋，尹 楠

（1.中國科學(xué)院西北生態(tài)環(huán)境資源研究院凍土工程國家重點實驗室，甘肅蘭州730000；2.江蘇蘇邑設(shè)計集團(tuán)有限公司，江蘇南京210041）

0 引言

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

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

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

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

1 數(shù)學(xué)模型及方程

1.1 數(shù)學(xué)傳熱離散元模型

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

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

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

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

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

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

假定溫度對應(yīng)變變化的影響很小，可以忽略，這個假定應(yīng)用于涉及到固體和流體的準(zhǔn)靜態(tài)力學(xué)問題，則這種連續(xù)介質(zhì)的熱傳導(dǎo)控制方程為

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

將式（4）代入式（5），得到熱傳導(dǎo)方程。

1.2 顆粒粘結(jié)模型

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

圖1 凍土中膠結(jié)結(jié)構(gòu)示意圖［23］Fig.1 Sketch map of the cementation structure in frozen soils［23］

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

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

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

2 數(shù)值計算模型

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

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

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

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

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

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

式中：th為時間變量?？赏ㄟ^調(diào)整α0來改變初始時刻對應(yīng)的日期，α0=0對應(yīng)的初始日期為7月15日。

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

路基頂面CD的溫度變化規(guī)律為

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

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

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

表1 路基各介質(zhì)中的離散元熱學(xué)微觀參數(shù)Table 1 DEM thermal microcosmic parameters in the embankment

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

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

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

幾種不同介質(zhì)土層的顆粒接觸剛度與溫度之間的數(shù)學(xué)關(guān)系，其變化規(guī)律具有一致性。即在溫度為正時，顆粒剛度為常數(shù)；溫度為負(fù)時，顆粒接觸剛度隨溫度的降低而線性增大。

3 結(jié)果與分析

3.1 數(shù)值模型校核

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

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

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

圖6 路基路肩孔地溫變化曲線（圖中天然地表面為縱坐標(biāo)的原點）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年四個典型季節(jié)的溫度分布圖進(jìn)行分析。從路基等溫線變化情況來看，四個季節(jié)都形成-0.8℃溫度線。從路基融化深度來看，在7月15日，上部路基內(nèi)融化深度約為2.15 m，而在10月15日，上部路基內(nèi)融化深度約為3.20 m。

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

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

圖8 中，為凍土路基修筑完成后第20年四個典型季節(jié)的溫度場，可見，在這四個典型季節(jié)中，-0.8℃溫度線已退化，-0.5℃溫度線退化也很嚴(yán)重，而且上部路基內(nèi)融化深度有一定下移，其中，在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，路基存在凍土退化的問題。這些離散元研究結(jié)果與路基已有的研究結(jié)果一致，可以說明，運用離散單元法對多年凍土區(qū)路基的熱-力耦合數(shù)值仿真較合理，也驗證了隨路基溫度場變化的熱-力耦合離散元計算模型的準(zhǔn)確性，并且顆粒集合體之間的熱交換作用特性也較真實體現(xiàn)凍土路基土顆粒的實際變化情況。

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

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

圖10 為多年凍土區(qū)路基建成后第20年10月15日的接觸力鏈分布，其余三個季節(jié)變化規(guī)律類似。路基中接觸力鏈的粗細(xì)可表征路基顆粒間接觸力的大小，也能從微細(xì)觀水平反映出路基的宏觀變形特征?？梢姡谥亓ψ饔孟拢坊纬蓮娏︽溑c弱力鏈共同發(fā)展的體系。上部路基內(nèi)的顆粒體會向下與向右運動，擠壓基礎(chǔ)內(nèi)的顆粒體，在路基坡腳以左，路基中的柱狀力鏈呈現(xiàn)向左的傾斜特征，在路基中心處，柱狀力鏈呈垂直分布，從路基中心至坡腳處，柱狀力鏈向左的傾斜角度逐漸減小，表明了這部分路基中顆粒間擠壓作用向右逐漸減弱；在路基坡腳以右，柱狀力鏈向左的傾斜變化特征逐漸消失，形成柱狀力鏈之間交叉分布的特征，體現(xiàn)了在豎向重力作用下路基顆粒間接觸力的變化特點。而路基中網(wǎng)狀弱力鏈分布在柱狀強力鏈周圍，并與強力鏈處處連接，對強力鏈的穩(wěn)定起輔助作用，從而使路基發(fā)生整體變形。還可知，路基中力鏈的粗細(xì)分布不均勻，各向異性很明顯，坡腳以左路基各向異性的程度較坡腳以右路基中的大，表明坡腳以左路基中顆粒間相互作用較強烈。

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

4 結(jié)論

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

（1）對比有限單元法和離散單元法可發(fā)現(xiàn)，離散單元法突破了連續(xù)介質(zhì)力學(xué)中變形條件的限制，可模擬顆粒之間的相互運動，因此能夠較真實反映出路基的熱、力變化特征，可應(yīng)用于寒區(qū)路基工程的長期穩(wěn)定性研究。

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

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

（4）凍土路基在坡腳以左，顆粒體的相互作用較劇烈，路基（沉降）變形主要發(fā)生在路基主體結(jié)構(gòu)及下部基礎(chǔ)一定范圍內(nèi)。

本文的研究具有一定的探索性，所選取斷面的熱邊界條件具有對稱性，并且未考慮陰陽坡效應(yīng)和水分的作用，因此還需要對該系統(tǒng)進(jìn)一步的研究。作為初步研究，可為今后凍土工程的穩(wěn)定性研究和寒區(qū)工程的建設(shè)提供科學(xué)依據(jù)。