冻土路基温度场随机有限元分析

冻土路基温度场随机有限元分析

1冻土温度场随机有限元分析青藏铁路的建设是西藏青藏公路之后的另一条经济和社会路线。青藏铁路穿越550km多年冻土区,由于冻土的特殊性,多年冻土路基稳定性一直是人们所关心的重点。冻土是一种对温度和冰比较敏感的不稳定土体,温度的变化会引起路基土体工程性质产生变化,导致路基各部分产生不均匀变形,当这种变形超过一定限度时,路基便会遭到破坏,这种破坏是冻土路基病害的一个主要原因。如何使温度场的计算更能反映实际情况,对冻土路基的稳定性分析有重要的意义。在传统的温度场计算中,为了计算方便,通常把环境温度和冻土的热学参数作为确定的量来求解。然而在实际工程中,环境温度和路基土的热学参数表现出很大的变异性。在这种情况下传统的温度场求解方法就会使结果出现偏差,为了使温度场的求解更加精确,必须借助概率统计的方法。随机有限元法结合了概率统计方法和有限元法的优点,可以灵活地考虑系统的材料、荷载等的随机性,比较适合用来解决此类问题。冻土路基温度场随机有限元分析的实质就是在有限元法的基础上,考虑土体物热学参数和边界条件的随机性,计算出路基内部结点温度的统计特征,为进一步计算路基稳定的可靠度做准备。在运用随机有限元法进行温度场的计算方面,国内外已有学者进行了研究。Suzuki,Oktayy以及Shigekazu等人率先用随机有限元估算结构的随机温度场;刘宁、刘光廷等提出了大体积温度场的随机有限元算法。但在冻土温度场的计算方面,关于随机有限元的应用国内尚少有人涉及。本文在借鉴前人成果基础上,鉴于一阶Taylor展开随机有限元法具有计算量小和计算速度快的特点,推导了基于一阶Taylor展开的冻土温度场随机有限元法的基本公式,编制了相应的程序,并结合青藏铁路北麓河试验段路基实例进行了冻土随机温度场以及温度方差的计算。2随机变量均值的函数一阶Taylor展开随机有限元法是将数值模型中的控制量在随机变量均值点处进行一阶Taylor级数展开,经过适当的数学处理得出所需要的计算公式。假定X=(X1,X2,…,Xn)为随机变量,f(X)为随机变量的函数,则f(X)在随机变量均值ˉ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的协方差。3冷冻土基温场的随机限制法3.1变相储-显热容法的求解冻土中水的相变潜热多年冻土温度场与一般温度场的区别在于多年冻土中存在一个冻融活动相界面,它的冻结和融化过程是一个具有相变的热传导过程,其温度场数学模型可以用如下的二维热传导方程来表示: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)相变界面温度连续条件:Τf=Τu=Τc(x,y)∈h(t)(6)相变界面热流连续条件:(λf∂Τf∂t-λu∂Τu∂t)=Lρ(W-Wu)dh(t)dt(x,y)∈h(t)(7)式中,f,u表示路基土的冻、融两相状态,Tc为土体冻结临界温度,λX,λy为x,y方向的导热系数,ρ为路基土干容重,c为容积热容量,t为时间,Ω是计算区域,h(t)表示随时间变化的相变界面,W和Wu分别为土体的重量含水量和重量未冻水含量,L代表水的相变潜热。因为在求解域存在一个随时间变化的活动相界面,这个热传导问题的求解比一般的变系数导热问题要复杂的多。为了处理这个问题,本文采用了显热容法,它的基本思路是:由于冻土的相变不是严格地在某一特定的温度条件下发生,而是在剧烈相变区的一个较小的温度范围内发生,因此可以把冻土中水的相变潜热看作是一个小的温度范围内的大热容。这样就可以把式(4)～(7)转化为在整个求解区域上适用的非线性导热方程: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是土中水发生相变的温度,ΔT为土中水发生剧烈相变温度区间的长度。3.2初支地温阵风下的土壤温度波根据附面层原理,可得路基上表面边界条件:Τ(t)=Τ0+R0t+A0sin(2π365t-35π)(11)其中T0为附面层底的初始年平均地温,R0为气候变化引起的附面层底的地温变化率,A0为附面层底地温振幅。计算区域两侧铅垂面边界取为绝热边界:∂Τ∂n=0(12)n表示侧面边界的外法线方向。下边界取恒定热流边界条件:∂Τ∂n=Gs(13)Gs为计算区域下边界地温梯度。3.3热学参数及特性参数求解温度场利用边界条件(11)-(13)和计算区域土体的热学参数及特性参数求解温度场的Laplace方程(14),求得的温度场作为初始温度场。∂∂x(λx∂Τ∂y)+∂∂y(λy∂Τ∂y)=0(14)3.4边界热流密度由于热容和导热系数都与温度有关,方程(8)是非线性的,其定解问题可以转化为以下任意时刻范函的极小值:Ι(Τ)=12∫∫Ω{λx(∂Τ∂x)2+λy(∂Τ∂y)2+cΤ∂Τ∂t}dxdy(15)将求解区域离散为三角形单元,可得在二维情况下温度场T(x,y,t)计算的有限元支配方程为:ΚΤ+Ν∂Τ∂t=Ρ(16)Κ=∑e∫∫Ωe(λx∂Νi∂x∂Νj∂x+λy∂Νi∂y∂Νj∂y)dxdyΝ=∑e∫∫ΩecΝiΝjdxdyΡ=∑e∫ΓeqΝidΓq为边界热流密度。把式(16)在时间域采用Crank-Nicolson法进行离散,可得:(Ν+Δt2Κ)Τn+1=(Ν-Δt2Κ)Τn+Δt2(Ρn+1+Ρn)(17)式(17)两边取均值:(Ν¯+Δt2Κ¯)Τ¯n+1=(Ν¯-Δt2Κ¯)Τ¯n+Δt2(Ρ¯n+1+Ρ¯n)(18)把(17)式两边同时对随机变量X=(X1,X2,…Xn)求导,整理得:∂Τ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)可以求出某一时刻温度场的均值Τ¯n,由方程(19)求出T′n之后,某一时刻任一点温度可以用式(20)求得:Τn=Τ¯n+∑i=1n∂Τn∂Xi(Xi-X¯i)(20)任一点温度的方差为Var[Τkn]=∑i=1n∑j=1n∂Τkn∂Xi|Xi=X¯i⋅∂Τkn∂Xj|Xj=X¯jCov(Xi,Xj)(21)任意两点温度的协方差为:Var[Τkn,Τln]=∑i=1n∑j=1n∂Τkn∂Xi|Xi=X¯i⋅∂Τln∂Xj|Xj=X¯jCov(Xi,Xj)(22)4温度场计算结果根据上述计算方法,作者用FORTRAN语言编制了冻土路基温度场随机有限元计算程序,取青藏铁路北麓河试验段冻土路基剖面进行了随机温度场及标准差计算。计算区域路基为不对称结构,阳坡低阴坡高。取路面以下20m,宽40m,其中路基高5m,路面宽12m,路基以下土层按土性共分为3层,土的干容重和热学参数作为随机变量,分层情况以及各层土的相关参数值见表1。计区域用自动剖分程序共剖分为1695个单元,906个结点,有限元网格见图1。随机场单元采用不规则四边形离散,每个随机场单元包含两个有限元单元。对于初始边界条件,本文作为常量来处理。计算区域左右铅垂面边界取为绝热边界,路基底部地温梯度值为0.03℃·m-1,上边界取为第三类边界条件,初始年平均气温为-3.2℃,阳坡附面层底年平均地温为2.5℃,阴坡为-1℃,附面层底地温增温率为0.02℃·a-1,边界地温振幅为12.75℃。按照以上计算条件,作者对该路基完成以后的随机温度场进行了计算。考虑到短期内的温度场计算值受初始值的影响较大,同时因为实测资料主要局限在第二年,为了对比分析的方便,本文主要考虑第二年的计算结果。图2分别为选取的有代表性的负温期、稳定融化期和正温期末即第二年3月19日、7月19日、10月19日路基的温度场和标准差等值线。对比中国科学院寒旱所的现场实测资料,本文温度场计算结果与实测值基本吻合。其中路基中心温度钻孔数据见表2。从图2中可以看出,在冻土路基完成一段时间以后,温度标准差的分布趋势与温度场的分布趋势大致相同,主要表现在垂直方向上的变化。但两者的分布略有不同,标准差的最大值主要出现在坡肩和靠近路面附近,而在路基下部则比较均一,基本都在0.05以下,表明路面附近和路肩附近温度变化比较剧烈,而路基下部比较平稳。主要原因在于坡肩和路面直接与外部环境接触,温度受外界环境温度随机性的影响较大,而随着深度的增加,外界环境温度的影响逐渐减弱,标准差主要受路基土参数随机性的影响,由于土参数的随机性相对外部环境来说要小的多,所以标准差较小。本文假定路基土参数的标准差均为0.2,所以温度标准差的值在路基下部也基本相同。从本文结果可以看出,对冻土路基温度场影响最大的是外部环境的温度,路基土温度变异性最大的范围在路基上部靠近路面和路肩的地方,因此在实际的工程实践中,应该采取在路基上部铺设保温材料和片石护坡等措施,尽量减小外界环境温度对下部冻土的影响,可望有效控制冻土路基温度的变异,保护冻土,提高路基稳定性。5不同温度场的计算结果冻土路基温度场受到很多随机因素的影响,传统的温度场计算方法往往不能够反映实际情况,本文引进一阶