凍土路基溫度場(chǎng)隨機(jī)有限元分析

凍土路基溫度場(chǎng)隨機(jī)有限元分析

1凍土溫度場(chǎng)隨機(jī)有限元分析青藏鐵路的建設(shè)是西藏青藏公路之后的另一條經(jīng)濟(jì)和社會(huì)路線。青藏鐵路穿越550km多年凍土區(qū),由于凍土的特殊性,多年凍土路基穩(wěn)定性一直是人們所關(guān)心的重點(diǎn)。凍土是一種對(duì)溫度和冰比較敏感的不穩(wěn)定土體,溫度的變化會(huì)引起路基土體工程性質(zhì)產(chǎn)生變化,導(dǎo)致路基各部分產(chǎn)生不均勻變形,當(dāng)這種變形超過一定限度時(shí),路基便會(huì)遭到破壞,這種破壞是凍土路基病害的一個(gè)主要原因。如何使溫度場(chǎng)的計(jì)算更能反映實(shí)際情況,對(duì)凍土路基的穩(wěn)定性分析有重要的意義。在傳統(tǒng)的溫度場(chǎng)計(jì)算中,為了計(jì)算方便,通常把環(huán)境溫度和凍土的熱學(xué)參數(shù)作為確定的量來求解。然而在實(shí)際工程中,環(huán)境溫度和路基土的熱學(xué)參數(shù)表現(xiàn)出很大的變異性。在這種情況下傳統(tǒng)的溫度場(chǎng)求解方法就會(huì)使結(jié)果出現(xiàn)偏差,為了使溫度場(chǎng)的求解更加精確,必須借助概率統(tǒng)計(jì)的方法。隨機(jī)有限元法結(jié)合了概率統(tǒng)計(jì)方法和有限元法的優(yōu)點(diǎn),可以靈活地考慮系統(tǒng)的材料、荷載等的隨機(jī)性,比較適合用來解決此類問題。凍土路基溫度場(chǎng)隨機(jī)有限元分析的實(shí)質(zhì)就是在有限元法的基礎(chǔ)上,考慮土體物熱學(xué)參數(shù)和邊界條件的隨機(jī)性,計(jì)算出路基內(nèi)部結(jié)點(diǎn)溫度的統(tǒng)計(jì)特征,為進(jìn)一步計(jì)算路基穩(wěn)定的可靠度做準(zhǔn)備。在運(yùn)用隨機(jī)有限元法進(jìn)行溫度場(chǎng)的計(jì)算方面,國內(nèi)外已有學(xué)者進(jìn)行了研究。Suzuki,Oktayy以及Shigekazu等人率先用隨機(jī)有限元估算結(jié)構(gòu)的隨機(jī)溫度場(chǎng);劉寧、劉光廷等提出了大體積溫度場(chǎng)的隨機(jī)有限元算法。但在凍土溫度場(chǎng)的計(jì)算方面,關(guān)于隨機(jī)有限元的應(yīng)用國內(nèi)尚少有人涉及。本文在借鑒前人成果基礎(chǔ)上,鑒于一階Taylor展開隨機(jī)有限元法具有計(jì)算量小和計(jì)算速度快的特點(diǎn),推導(dǎo)了基于一階Taylor展開的凍土溫度場(chǎng)隨機(jī)有限元法的基本公式,編制了相應(yīng)的程序,并結(jié)合青藏鐵路北麓河試驗(yàn)段路基實(shí)例進(jìn)行了凍土隨機(jī)溫度場(chǎng)以及溫度方差的計(jì)算。2隨機(jī)變量均值的函數(shù)一階Taylor展開隨機(jī)有限元法是將數(shù)值模型中的控制量在隨機(jī)變量均值點(diǎn)處進(jìn)行一階Taylor級(jí)數(shù)展開,經(jīng)過適當(dāng)?shù)臄?shù)學(xué)處理得出所需要的計(jì)算公式。假定X=(X1,X2,…,Xn)為隨機(jī)變量,f(X)為隨機(jī)變量的函數(shù),則f(X)在隨機(jī)變量均值ˉX處的一階Taylor展開式為:f(X)=f(ˉX)+n∑i=1?f?Xi|Xi=ˉXi(Xi-ˉXi)(1)取均值得:E[f(X)]=f(ˉX)(2)f(X)的方差為:Var[f(X)]=n∑i=1n∑j=1?f?Xi|Xi=ˉXi??f?Xj|Xj=ˉXjCov(Xi,Xj)(3)式中,Cov(Xi,Xj)為Xi,Xj的協(xié)方差。3冷凍土基溫場(chǎng)的隨機(jī)限制法3.1變相儲(chǔ)-顯熱容法的求解凍土中水的相變潛熱多年凍土溫度場(chǎng)與一般溫度場(chǎng)的區(qū)別在于多年凍土中存在一個(gè)凍融活動(dòng)相界面,它的凍結(jié)和融化過程是一個(gè)具有相變的熱傳導(dǎo)過程,其溫度場(chǎng)數(shù)學(xué)模型可以用如下的二維熱傳導(dǎo)方程來表示:cf?Τf?t=??x(λxf?Τf?x)+??y(λyf?Τf?y)t>0,(x,y)∈Ωf(4)cu?Τu?t=??x(λxu?Τu?x)+??y(λyu?Τu?y)t>0,(x,y)∈Ωu(5)相變界面溫度連續(xù)條件:Τf=Τu=Τc(x,y)∈h(t)(6)相變界面熱流連續(xù)條件:(λf?Τf?t-λu?Τu?t)=Lρ(W-Wu)dh(t)dt(x,y)∈h(t)(7)式中,f,u表示路基土的凍、融兩相狀態(tài),Tc為土體凍結(jié)臨界溫度,λX,λy為x,y方向的導(dǎo)熱系數(shù),ρ為路基土干容重,c為容積熱容量,t為時(shí)間,Ω是計(jì)算區(qū)域,h(t)表示隨時(shí)間變化的相變界面,W和Wu分別為土體的重量含水量和重量未凍水含量,L代表水的相變潛熱。因?yàn)樵谇蠼庥虼嬖谝粋€(gè)隨時(shí)間變化的活動(dòng)相界面,這個(gè)熱傳導(dǎo)問題的求解比一般的變系數(shù)導(dǎo)熱問題要復(fù)雜的多。為了處理這個(gè)問題,本文采用了顯熱容法,它的基本思路是:由于凍土的相變不是嚴(yán)格地在某一特定的溫度條件下發(fā)生,而是在劇烈相變區(qū)的一個(gè)較小的溫度范圍內(nèi)發(fā)生,因此可以把凍土中水的相變潛熱看作是一個(gè)小的溫度范圍內(nèi)的大熱容。這樣就可以把式(4)～(7)轉(zhuǎn)化為在整個(gè)求解區(qū)域上適用的非線性導(dǎo)熱方程:c?Τ?t=??x(λx?Τ?x)+??y(λy?Τ?y)(8)c={cfΤ<Τm-ΔΤ/2Lρ(W-Wu)ΔΤ+cu+cf2Τm-ΔΤ/2<Τ<Τm+ΔΤ/2cuΤ>Τm-ΔΤ/2(9)λ={λfΤ<Τm-ΔΤ/2λf+(λf+λu)ΔΤ(Τ-(Τm-ΔΤ/2))Τm-ΔΤ/2<Τ<Τm+ΔΤ/2λuΤ>Τm-ΔΤ/2(10)Tm是土中水發(fā)生相變的溫度,ΔT為土中水發(fā)生劇烈相變溫度區(qū)間的長度。3.2初支地溫陣風(fēng)下的土壤溫度波根據(jù)附面層原理,可得路基上表面邊界條件:Τ(t)=Τ0+R0t+A0sin(2π365t-35π)(11)其中T0為附面層底的初始年平均地溫,R0為氣候變化引起的附面層底的地溫變化率,A0為附面層底地溫振幅。計(jì)算區(qū)域兩側(cè)鉛垂面邊界取為絕熱邊界:?Τ?n=0(12)n表示側(cè)面邊界的外法線方向。下邊界取恒定熱流邊界條件:?Τ?n=Gs(13)Gs為計(jì)算區(qū)域下邊界地溫梯度。3.3熱學(xué)參數(shù)及特性參數(shù)求解溫度場(chǎng)利用邊界條件(11)-(13)和計(jì)算區(qū)域土體的熱學(xué)參數(shù)及特性參數(shù)求解溫度場(chǎng)的Laplace方程(14),求得的溫度場(chǎng)作為初始溫度場(chǎng)。??x(λx?Τ?y)+??y(λy?Τ?y)=0(14)3.4邊界熱流密度由于熱容和導(dǎo)熱系數(shù)都與溫度有關(guān),方程(8)是非線性的,其定解問題可以轉(zhuǎn)化為以下任意時(shí)刻范函的極小值:Ι(Τ)=12∫∫Ω{λx(?Τ?x)2+λy(?Τ?y)2+cΤ?Τ?t}dxdy(15)將求解區(qū)域離散為三角形單元,可得在二維情況下溫度場(chǎng)T(x,y,t)計(jì)算的有限元支配方程為:ΚΤ+Ν?Τ?t=Ρ(16)Κ=∑e∫∫Ωe(λx?Νi?x?Νj?x+λy?Νi?y?Νj?y)dxdyΝ=∑e∫∫ΩecΝiΝjdxdyΡ=∑e∫ΓeqΝidΓq為邊界熱流密度。把式(16)在時(shí)間域采用Crank-Nicolson法進(jìn)行離散,可得:(Ν+Δt2Κ)Τn+1=(Ν-Δt2Κ)Τn+Δt2(Ρn+1+Ρn)(17)式(17)兩邊取均值:(Νˉ+Δt2Κˉ)Τˉn+1=(Νˉ-Δt2Κˉ)Τˉn+Δt2(Ρˉn+1+Ρˉn)(18)把(17)式兩邊同時(shí)對(duì)隨機(jī)變量X=(X1,X2,…Xn)求導(dǎo),整理得:?Τn+1?Xi=(Ν+Δt2Κ)-1[(?Ν?Xi-Δt2?Κ?Xi)Τn+(Ν-Δt2Κ)?Τn?Xi-(?Ν?Xi+Δt2?Κ?Xi)Τn+1+Δt2(?Ρn+1?Xi+?Ρn?Xi)](19)由方程(18)可以求出某一時(shí)刻溫度場(chǎng)的均值Τˉn,由方程(19)求出T′n之后,某一時(shí)刻任一點(diǎn)溫度可以用式(20)求得:Τn=Τˉn+∑i=1n?Τn?Xi(Xi-Xˉi)(20)任一點(diǎn)溫度的方差為Var[Τkn]=∑i=1n∑j=1n?Τkn?Xi|Xi=Xˉi??Τkn?Xj|Xj=XˉjCov(Xi,Xj)(21)任意兩點(diǎn)溫度的協(xié)方差為:Var[Τkn,Τln]=∑i=1n∑j=1n?Τkn?Xi|Xi=Xˉi??Τln?Xj|Xj=XˉjCov(Xi,Xj)(22)4溫度場(chǎng)計(jì)算結(jié)果根據(jù)上述計(jì)算方法,作者用FORTRAN語言編制了凍土路基溫度場(chǎng)隨機(jī)有限元計(jì)算程序,取青藏鐵路北麓河試驗(yàn)段凍土路基剖面進(jìn)行了隨機(jī)溫度場(chǎng)及標(biāo)準(zhǔn)差計(jì)算。計(jì)算區(qū)域路基為不對(duì)稱結(jié)構(gòu),陽坡低陰坡高。取路面以下20m,寬40m,其中路基高5m,路面寬12m,路基以下土層按土性共分為3層,土的干容重和熱學(xué)參數(shù)作為隨機(jī)變量,分層情況以及各層土的相關(guān)參數(shù)值見表1。計(jì)區(qū)域用自動(dòng)剖分程序共剖分為1695個(gè)單元,906個(gè)結(jié)點(diǎn),有限元網(wǎng)格見圖1。隨機(jī)場(chǎng)單元采用不規(guī)則四邊形離散,每個(gè)隨機(jī)場(chǎng)單元包含兩個(gè)有限元單元。對(duì)于初始邊界條件,本文作為常量來處理。計(jì)算區(qū)域左右鉛垂面邊界取為絕熱邊界,路基底部地溫梯度值為0.03℃·m-1,上邊界取為第三類邊界條件,初始年平均氣溫為-3.2℃,陽坡附面層底年平均地溫為2.5℃,陰坡為-1℃,附面層底地溫增溫率為0.02℃·a-1,邊界地溫振幅為12.75℃。按照以上計(jì)算條件,作者對(duì)該路基完成以后的隨機(jī)溫度場(chǎng)進(jìn)行了計(jì)算?？紤]到短期內(nèi)的溫度場(chǎng)計(jì)算值受初始值的影響較大,同時(shí)因?yàn)閷?shí)測(cè)資料主要局限在第二年,為了對(duì)比分析的方便,本文主要考慮第二年的計(jì)算結(jié)果。圖2分別為選取的有代表性的負(fù)溫期、穩(wěn)定融化期和正溫期末即第二年3月19日、7月19日、10月19日路基的溫度場(chǎng)和標(biāo)準(zhǔn)差等值線。對(duì)比中國科學(xué)院寒旱所的現(xiàn)場(chǎng)實(shí)測(cè)資料,本文溫度場(chǎng)計(jì)算結(jié)果與實(shí)測(cè)值基本吻合。其中路基中心溫度鉆孔數(shù)據(jù)見表2。從圖2中可以看出,在凍土路基完成一段時(shí)間以后,溫度標(biāo)準(zhǔn)差的分布趨勢(shì)與溫度場(chǎng)的分布趨勢(shì)大致相同,主要表現(xiàn)在垂直方向上的變化。但兩者的分布略有不同,標(biāo)準(zhǔn)差的最大值主要出現(xiàn)在坡肩和靠近路面附近,而在路基下部則比較均一,基本都在0.05以下,表明路面附近和路肩附近溫度變化比較劇烈,而路基下部比較平穩(wěn)。主要原因在于坡肩和路面直接與外部環(huán)境接觸,溫度受外界環(huán)境溫度隨機(jī)性的影響較大,而隨著深度的增加,外界環(huán)境溫度的影響逐漸減弱,標(biāo)準(zhǔn)差主要受路基土參數(shù)隨機(jī)性的影響,由于土參數(shù)的隨機(jī)性相對(duì)外部環(huán)境來說要小的多,所以標(biāo)準(zhǔn)差較小。本文假定路基土參數(shù)的標(biāo)準(zhǔn)差均為0.2,所以溫度標(biāo)準(zhǔn)差的值在路基下部也基本相同。從本文結(jié)果可以看出,對(duì)凍土路基溫度場(chǎng)影響最大的是外部環(huán)境的溫度,路基土溫度變異性最大的范圍在路基上部靠近路面和路肩的地方,因此在實(shí)際的工程實(shí)踐中,應(yīng)該采取在路基上部鋪設(shè)保溫材料和片石護(hù)坡等措施,盡量減小外界環(huán)境溫度對(duì)下部凍土的影響,可望有效控制凍土路基溫度的變異,保護(hù)凍土,提高路基穩(wěn)定性。5不同溫度場(chǎng)的計(jì)算結(jié)果凍土路基溫度場(chǎng)受到很多隨機(jī)因素的影響,傳統(tǒng)的溫度場(chǎng)計(jì)算方法往往不能夠反映實(shí)際情況,本文引進(jìn)一階