搜索
中国科学D辑地球科学2006,36(6):587~592587冻土路基的随机温度场*刘志强①②赖远明①**张明义①张学富③陆昊②(①中国科学院寒区旱区环境与工程研究所冻土工程国家重点实验室,兰州730000;②兰州交通大学土木工程学院,兰州730070;③重庆交通学院,重庆400074)摘要在导出温度场的变分原理基础上,利用摄动技术将冻土热力学参数及边界条件的随机性引入到温度场泛函的变分中,得到了随机温度场的变分原理和有限元公式.以此编制了计算程序,并计算了寒区路基在随机边界条件和材料参数影响下的随机温度场.关键词冻土路基随机温度场有限元青藏铁路收稿日期:2005-06-10;接受日期:2005-08-30*国家杰出青年科学基金(批准号:40225001)、中国科学院知识创新工程重要方向项目(批准号:KZCX3-SW-351)、中国科学院“百人计划”资助项目、中国科学院知识创新工程重大项目(批准号:KZCX1-SW-04)和中国科学院全国百篇优秀博士论文专项资金项目共同资助**联系人,E-mail:ymlai@ns.lzb.ac.cn寒区工程建设,尤其是青藏公路、铁路、隧道的建设中,提出许多积极保护冻土的工程措施,如抛石、通风管、热棒、旱桥、保温板、遮阳板等等.对每种形式的工程措施进行热状况及热稳定性的理论计算、分析、预测是决策、施工及养护的重要依据.从现有文献看,目前对寒区工程温度场的计算都进行的是确定性计算、分析[1,2],即模型的各种物性参数、尤其是边界条件都是确定的,计算结果也是确定的.然而,冻土的物性参数和边界条件有很大的随机性、变异性,其结果也应是具有一定的变异性.真实的反映这些随机因素,采用概率的方法,计算结构的随机温度场,估算可能出现的摆幅,对进行工程可靠度计算是十分必要的.Sluzalec采用对温度场有限元公式直接进行小参数摄动的方法,得到随机有限元公式,这种方法简单明了,但不能将各种随机因素(如环境温度)都考虑到有限元公式中[3].刘宁在考虑环境温度、混凝土绝热升温及材料热学参数随机性的影响下,对混凝土重力坝的随机温度场进行了计算[4];Hien和Kleiber则在建立温度场泛函的变分基础上,应用摄动法得到随机温度场的变分原理,从而导出随机温度场的有限元公式[5];Kaminski和Hien用同样的方法对复合材料的随机温度场进行了计算[6].本文以冻土路基为对象,将边界条件和材料参数的随机性模拟为随机场或随机变量,采用摄动随机有限元法对其随机温度场进行了计算.1温度场的变分原理为了将材料的热力学参数、边界条件等随机因素直接引入到有限元公式中,首先建立温度场的变分原理,在此基础上建立随机温度场的变分原理.平面瞬态温度场导热微分方程为(不计内热源项):SCIENCEINCHINASer.DEarthSciences588中国科学D辑地球科学第36卷TTCtxxyyλλ⎛⎞∂∂∂∂∂⎛⎞=+⎜⎜⎟∂∂∂∂∂⎝⎠⎝⎠T⎟,(1)根据所研究的问题,边界条件为(只考虑温度边界条件和热流边界条件):在Γ1上,(2)wTT=在Γ2上Tqnλ∂−=∂,(3)在Γ3上(aTTTnλα∂−=−∂),(4)式中:C为容积热容量,Tw为边界温度,Ta为外界温度,α为对流换热系数,λ为导热系数,q为边界热流密度.初始条件:00tT==T.(5)考虑温度场的任意满足强制边界条件(2)的变化δT,由(1)式可得dTTTCtxxyyλλδΩ⎡⎤⎛⎞∂∂∂∂∂⎛⎞−−Ω⎢⎜⎟⎜⎟∂∂∂∂∂⎝⎠⎝⎠⎣⎦∫0T=⎥,(6)对上式的后两项应用分部积分和格林公式,引入边界条件(3),(4)式可得TTTTCTtxxxyyδδλδλΩ⎡⎤⎛⎞⎛⎞∂∂∂∂∂∂∂⎛⎞⎛⎞++⎢⎥⎜⎟⎜⎟⎜⎟⎜⎟∂∂∂∂∂∂∂⎝⎠⎝⎠⎝⎠⎝⎠⎣⎦∫Ty0,(7)23dd()daTqTTTδαδΓΓ⋅Ω+Γ+−Γ=∫∫上式即为温度场泛函的一阶变分方程δJ=0,温度场泛函J的表达式为22d2TTTJCTtxyλΩ⎧⎫⎡⎤⎛⎞∂∂∂⎪⎪⎛⎞⎢⎥=++⎨⎬⎜⎟⎜⎟∂∂∂⎝⎠⎢⎥⎝⎠⎪⎪⎣⎦⎩⎭∫Ω232d2aTTqTTαΓΓ⎛⎞+Γ+−⎜⎝⎠∫∫dΓ⎟(8)2随机温度场的变分原理设br,r=1,…,R,为R个随机变量.将热容C,导热系数λ,对流换热系数α,外界温度Ta和温度T等作为随机变量的函数,在br的均值rb处进行一阶泰勒展开,并写成摄动形式有:1ΔrRrrrbTΤTbbTTγγ=⎛⎞∂′=+=+⎜⎟∂⎝⎠∑(9)1ΔrRrrrbCCCbCCbγγ=⎛⎞∂′=+=+⎜⎟∂⎝⎠∑,(10)1rRrrrb__bbλλλγλγλ=⎛⎞∂′=+Δ=+⎜⎟∂⎝⎠∑,(11)1ΔrRrrrb_qqqqbqqbγγ=⎛⎞∂′=+=+⎜⎟∂⎝⎠∑,(12)1ΔrRaaaraarrb__TTTbTTbγγ=⎛⎞∂′=+=+⎜⎟∂⎝⎠∑,(13)式中,(irbbbγγΔ=−)r为随机变量br在均值处的一阶变量,(′)表示对随机函数的一阶变量.将(9)~(13)带入(7)式,比较同此项的系数,得随机温度场的变分原理,它们是:零次变分原理dTTTTTCTtxxyyδλδδΩ⎧⎫⎡⎤⎛⎞∂∂∂∂∂⎪⎪⎛⎞++Ω⎨⎬⎢⎥⎜⎟⎜⎟∂∂∂∂∂⎝⎠⎝⎠⎪⎪⎣⎦⎩⎭∫()23daqTTTTδαδΓΓd=−Γ−−∫∫Γ;(14)一次变分原理33()d()ddaaTTTTTCTtxxyyTTTTTTTTTTTCTtxxyyδλδδαδαδδλδδΩΓΓΩ⎧⎫⎡⎤′′′⎛⎞∂∂∂∂∂⎪⎪⎛⎞d++Ω⎨⎬⎢⎥⎜⎟⎜⎟∂∂∂∂∂⎝⎠⎝⎠⎪⎪⎣⎦⎩⎭′′′=−Γ−−Γ⎧⎫⎡⎤⎛⎞∂∂∂∂∂⎪⎪⎛⎞′′−++⎨⎬⎢⎥⎜⎟⎜⎟∂∂∂∂∂⎝⎠⎝⎠⎪⎪⎣⎦⎩⎭∫∫∫∫Ω2dqTδΓ′−Γ∫.(15)3随机温度场的有限元列式将温度场离散为Ne各单元、Np个结点,有1pNiiiTN==∑TT,(16)1pNiiiTN=′′=∑,(17)式中,Ni为整体形函数,对于指定的单元e,i是该单元SCIENCEINCHINASer.DEarthSciences第6期刘志强等:冻土路基的随机温度场589结点时,有;否则,,为单元形函数.将(16),(17)带入(14),(15)式整理可得温度场随机有限元列式eiNN=i0iN=eiN零次方程+=CTKTQ;(18)一次方程()rrrrr′′′′′+=−+CTKTQCTKT,,(19)1,2,,rR=(18),(19)式中各项的具体表达式dijijCCNNΩ=∫Ω,(20)3ddjjiiijijNNNNKNxxyyλαΩ∂∂⎛⎞∂∂=+Ω+⎜⎟∂∂∂∂⎝⎠∫∫NΓΓ,(21)32diaiiQTNqNαΓΓ=Γ−∫∫dΓ,(22)3()dd,jjiirijrrjjNNNNKxxyyλΩ∂∂⎛⎞∂∂′′′=+Ω+⎜⎟∂∂∂∂⎝⎠∫∫NNαΓΓΩ(23),(24)()drijrijCCNNΩ′′=∫3()()ddriraaririQTTNqNααΓΓ′′′′=−Γ−∫∫2Γ,(25)温度场的协方差矩阵为11(,)(,)ijRRijklklklTTCovTTCovbbbb==∂∂=∂∂∑∑.(26)(,)klCovbb是任意两个随机变量的协方差,根据给定的相关函数以及随机场单元的大小、间距,可以计算任意随机场单元之间的协方差.(26)式计算的是结点温度的协方差矩阵,其主对角元素即为结点的方差.需要注意的是,由于冻土可能发生相变,当考虑材料热学参数的随机性时,随机变量的协方差是温度的函数,需要根据不同时刻计算相应的协方差矩阵,使得计算量显著增加,这和无相变材料只需计算一次协方差矩阵是不同的.本文对相变的处理,用焓来构造热容,焓定义为热容对温度的积分,具体采用DelGiudice等提出的方法,详见参考文献[7].热容的表达式为22HTHTxxyyCTTxy∂∂∂∂+∂∂∂∂=⎛⎞∂∂⎛⎞+⎜⎟⎜⎟∂∂⎝⎠⎝⎠,(27)上式中,H代表焓.4随机场的离散早期有限元计算中考虑随机因素时,将随机参数视为单一的随机变量而不是模拟成随机场,不能称为严格意义上的随机有限元.土性参数具有明显的空间变异性,宜将土性参数分布模拟为随机场而不是随机变量.在随机有限元的求解中,除了结构要进行有限元离散外,随机场也要进行离散.上面各式中的随机变量既是离散后的随机变量.目前随机场的离散方法主要有:(1)中心点离散法,(2)随机场的局部平均,(3)随机场的插值,(4)随机场的加权积分,(5)随机场的正交展开.局部平均法是将随机场单元的属性用随机场在单元上的局部平均来表征.局部平均法具有对随机场相关结构不敏感的特点,可以降低对数据的要求,目前采用较多.目前,一维和二维随机场的局部平均应用较多,三维也主要是将其转化为一维或二维的组合来处理.本文采用随机场的局部平均理论对随机场进行离散.5冻土路基随机温度场的计算及结果分析根据以上理论分析,笔者采用VC++6.0,应用面向对象的编程技术编制了“冻土路基随机温度场计算程序”,该程序可同时考虑冻土的材料参数、边界条件等随机因素,对冻土工程的随机温度场进行预测、计算.本文对青藏铁路某一海拔4500m的路基随机温度场进行了计算,路基示意图如图1所示(为简化计算,未考虑道碴层).图中区域为1填土-碎石砂砾,区域2为亚黏土,区域3为弱风化岩,它们的热力学参数见表1.计算两种随机因素作用下路基的随机温度场:随机的边界条件和随机材料参数.随机的边界条件主要考虑路基上边界(即图1中的A~F)气温的随机性,材料参数随机性视路基砂砾填土部分1和亚黏土2的比热容和导热系数为随机场.(W·m−1·℃−1)Cf(J·m−3·℃−1)λu(W·m−1·℃−1)Cu(J·m−3·℃−1)L(J·m−3)砂砾1.9801.913×1061.9192.227×10620.4×106亚黏土1.3511.879×1061.1252.357×10660.3×106弱风化岩1.8241.846×1061.4742.099×10637.7×106图1简化的路基结构示意图本文所讨论的路基为南北走向,不考虑其阴阳坡效应,故可认为问题具有对称性,取图1左半部为计算对象.随机场单元的离散和有限元单元相同.采用四结点等参元,单元划分都如图2所示.计算模型划分为876个单元,942个结点.图2单元划分为考虑边界条件的随机性,根据气温资料,通过统计回归分析,将气温模拟为正弦函数[1].计算考虑气候变暖,在未来50年气温上升2.6℃[8],气温表示为2ππ2.6sin876025036524hhTABtt⎛⎞=+++⎜⎟××⎝⎠,(28)式中A为年平均气温,B为平均年较差温度,第三项是升温项.A,B也作为相互独立的随机场来反映边界条件的随机性,根据曲线上的局部平均理论将其离散为随机变量,在本文中,随机场单元采用和有限元单元相同的划分.具体的,在本文考虑地区,年平均气温为A=−4.3℃,气温平均年较差24℃,即B=12℃,方差分别取为21Aσ=,22Bσ=.由附面层理论假定计算模型的边界条件如下[9]:天然地表AB和EF边的温度按下式变化nh2ππ2.6(3.3)sin876025036524TABtt⎛⎞=++++⎜⎟××⎝⎠h,(29)路堤斜坡BC和DE的温度变化规律为sh2ππ2.6(4.3)(1)sin876025036524TABtt⎛⎞=+++++⎜⎟××⎝⎠h,(30)路基中面CD的温度ph2ππ(5.5)(2.0)sin87602TABt⎛⎞=++++⎜⎟⎝⎠h2.65036524t+××(31)在ML边界上,地温梯度为0.03℃·m-1,AM,FL为绝热边界.对于土性参数的随机性,我们取路基填土和亚黏土的比热容和导热系数为相互独立的随机场.表1为各参数的均值,变异系数均取为0.1.并且只考虑随机场在竖直方向的相关性,标准相关函数取指数型,()exp(||/)ρξξθ=−,假定θ=10m.假设路基7月15日完工,图3是路基