冻土地区水渠的温度场和应力场数值分析

文章编号:100020240(2005)0620932207冻土地区水渠的温度场和应力场数值分析 收稿日期:2005210220;修订日期:2005211230 基金项目:国家自然科学基金项目(10202001);中国科学院知识创新工程重大项目(KZCX12SW 204)资助 作者简介:刘雄(1966-),男,江西萍乡人,高级工程师,1991年在南京航空学院获硕士学位,现为北京理工大学在职博士生,主要从事冻土和高速公路的研究工作.E 2mail :xiong0208@刘雄1,2,　宁建国1,　马巍3(1.北京理工大学爆炸科学与技术国家重点实验室,北京　100081;2.长沙理工大学,湖南长沙　410076;3.中国科学院寒区旱区环境与工程研究所冻土工程国家重点实验室,甘肃兰州　730000)摘　要:建立了含损伤的冻土弹性本构模型,对于在不同冰体积含量和不同温度下的冻结砂土,由该损伤本构模型计算的结果与实测的应力2应变曲线相吻合.利用自行开发的程序进行数值模拟,通过对渠道冻结水分场、温度场、应力场三场的数值计算结果分析,得到了比较准确且符合实际的温度场与应力场、位移场、应变场耦合的计算结果,与前人计算结果相符合,规律一致.算例表明,此程序能够计算冻土材料的相关物理量,且能很好地描述它们之间的关系,可为工程设计和施工提供必要的参考.关键词:冻土;损伤;耦合;温度场;应力场中图分类号:P642.14文献标识码:A1　引言 冻土的强度特性是冻土力学最重要的研究内容之一,是寒区工程建设中地基和基础设计的基本依据.虽然人们对冻土的强度等力学特性作了很多工作,但对于冻土中最基本的力学机理还没有认识得很清楚,对冻土力学机理的研究还存在很大缺陷和不足[1,2].在寒冷地区修筑公路铁路、铺设输油管道、修筑房屋等都离不开季节冻土和多年冻土方面的岩土工程技术.Foriero et al.[3]用有限元模拟了冻土边坡的蠕变;Klein [4]根据冻土的蠕变特性提出了一种弹性有限元方法;Rai et al.[5]导出了流体、热流和土体变形完全耦合的控制微分方程,给出了基于双孔隙率多孔热力学公式的解析解.Mas 2ters giers [6]把位移、孔隙压力和湿度值当作初始未知量,研究了土体变形、流体温度和孔隙压力之间的耦合规律.Lai Yuanming et al.[7,8]对寒区隧道的温度场和渗流场的耦合问题进行了非线性分析.何平等[9]、陈飞熊等[10]基于理论分析分别建立了饱和正冻土中的水、热、力场耦合模型.但是关于土体冻融过程中水、热、力三场耦合的力学机理问题的研究还很不深入. 本文针对冻土工程中急待解决的土体冻融过程中水、热、力三场耦合的力学机理问题开展研究,从复合材料细观力学出发,建立含损伤的冻土本构模型.根据传热学,渗流理论和冻土力学,建立了冻土温度场、水分场、应力场耦合问题的数学力学模型,并用自行开发的有限元程序,对冻土地区水渠的温度场和应力场进行数值分析研究.2　冻土的损伤本构模型研究2.1　冻土的弹性模量研究 首先将冻土看作是完全由土和冰组成的复合体单元,而暂时不考虑水的影响,并忽略如孔隙率等因素的影响.假设冻土体中的土颗粒为均匀连续体,土颗粒与冰之间完全粘结.令土和冰的弹性模量分别为E s 、E i ;剪切模量为G s ,G i ;泊松比为v s 、v i .下面各式中的下标s 和i 分别表示冻土中的土和冰的应力分量与应变分量.于是,由细观力学的混合律理论可以得到冻土复合材料的等效弹性模量K 和等效剪切模量G 为:K =<(K s ,K i )(1)第27卷　第6期2005年12月冰　川　冻　土J OU RNAL OF G L ACIOLO GY AND GEOCR YOLO GYVol.27　No.6Dec.2005G =φ(G s ,G i )(2)式(1)、(2)中的K s 、G s 、K i 、G i 分别为土和冰的体积模量和剪切模量. 根据各向同性材料的弹性常数之间的关系,有K =E/3(1-2v ).G =E/2(1+v )经推导后可以得到由土、冰的弹性模量和泊松比表示的冻土的等效弹性模量E 和等效泊松比v 的表达式为:E =[c s E s (1-2v i )+c i E i (1-2v s )]・[c s E s (1+v i )+c i E i (1+v s )]c s E s (1+v i )(1-2v i )+c i E i (1+v s )(1-2v s )(3)v =[c s E s v s (1-v i )(1-2v i )]+c i E i v i (1+v s )(1-2v s )c s E s (1+v i )(1-2v i )+c i E i (1+v s )(1-2v s )(4)式中:c s 和c i 分别为土和冰的体积含量(也称为体积分数),且满足关系式:c s +c i =1.2.2　冻土的损伤本构关系 根据Lemait re 等效应力原理,即应力σ作用在受损材料上引起的应变与等效应力 σ作用在无损材料上引起的应变等价.受损冻土在一维状态下的应力2应变关系为:ε=σ E = σE =σE (1-D )(5)式中: σ为损伤力学中的等效应力,此等效应力相当于土力学中的有效应力;E 为无损材料的弹性模量,即初始弹性模量,由式(3)确定; E 为受损材料的弹性模量,即有效弹性模量;D 为损伤变量. 因此,根据损伤力学理论[11],受损材料(D ≠0)的本构关系可以从无损材料(D =0)的本构方程导出,只要用损伤后的有效应力来取代无损材料本构关系中的名义应力即所谓的Cauchy 应力即可,由此得到冻土材料内部损伤型本构关系为:σ=(1-D )E ε(6)2.3　损伤演化规律 在外载作用过程中,由于材料强度服从概率统计上的Weibull 分布,因而可以认为材料的损伤变量D 也服从该分布.在此我们借助于混凝土损伤的一些研究结果[12],假定冻土材料的损伤变量D 服从双参数的Weibull 分布,即:D =1-exp -εan(7)式中:ε为应变;n 和a 分别为形状参数和尺度参数,均为非负数.经推导,可得:D =1-exp -1n εεf n(8) 上式即为以应变作为损伤演化控制量的损伤演化方程;εf 为应力峰值所对应的应变值;n 为表征材料损伤演化特征的材料参数.将(8)式代入(6)式中,从而得到冻土的损伤本构关系为:σ=E εexp -1n εεfn(9) 从上面的推导可以看出,冻土中任意点的应力σ与冻土的弹性模量、极限强度、应变峰值及该点的应变有关.冻土的损伤形状参数n 是弹性模量E 和割线模量E m 的函数,一般可由实验或细观力学方法加以确定.从理论上讲,n 值的大小与冻土材料的性态及外表破坏模式有一定的关系,即:n 值越大,冻土越趋向于脆性破坏;n 值越小,冻土越趋向于塑性破坏.2.4　计算实例及分析 表1是在不同温度下,选取不同的冰的体积含量c i ,利用式(3)和式(4)计算得到的孔隙比e 为0.4～0.5、天然含水量W 为15%～18%的冻结砂土表1　不同温度下冻结砂土的等效模量Table 1　Modulus of elasticity of f rozen sand soils under different temperatures冰体积含量/%冰温/℃冰弹性模量/MPa砂土弹性模量/MPa 冻结砂土弹性模量/MPa-2.04910.0291.5460.05-5.46342.046.0363.285-8.06820.0387.221-2.04910.0534.9570.1-5.46342.046.0678.237-8.06820.0726.0580-2.04910.01021.2580.2-5.46342.046.01307.694-8.06820.01403.3033396期刘雄等:冻土地区水渠的温度场和应力场数值分析等效弹性模量. 图1是利用表1的弹性模量值,由式(9)计算得到的应力2应变曲线,与在相同温度下试验得出的数据[8]作出的应力2应变曲线的比较.图1　冻结砂土应力2应变曲线(K i =0.05)Fig.1　Curves of stress -strain of f rozen sand soils under different temperatures (K i =0.05) 通过比较可以看出,冻土温度是影响冻土变形的重要因素.随着土温的降低,土中未冻水含量减小,由于冻土中冰含量的增加,固体颗粒间的胶结力增加,整个冻土体的变形减小.另外,随着冻土温度的降低,冻土体的破坏形式也由塑性破坏逐渐转向脆性破坏.在相同的应力水平下,冻结砂土的变形随着土温的降低而减小,弹性变形在总变形中所占的比例却随着土温的降低而增大,所占的残余变形减小.从图中我们还可以知道,弹性变形尽管随着冻土温度的降低而有很大增长,但其所占的比例仍然很小,特别是在冻土温度较高的情况下图2　考虑损伤后计算得到的应力2应变曲线1.K i =0.05;2.K i =0.1;3.K i =0.2Fig.2　Curves of calculated stress -strain underdifferent temperatures(一般高于-5℃范围内),弹性变形只占总变形量的10%～25%.在冻土温度较低时,它的弹性变形所占比例也不会有很大提高. 图2是不同冰体积含量下由式(9)计算的应力2应变曲线.可以看出,在不同的冻土温度下,由于冻土中冰含量的增加,冻土体的损伤强度也相应提高.但在相同的冻土温度下,随着应变的增大,冻土体的压缩应力也逐渐达到最高点(屈服点)并已经开始下降,并且不同冰体积含量的冻土的应力屈服点也非常接近,说明在冻土材料破坏(屈服)时冰的体积含量对冻土体的屈服强度影响不大.3　渠道冻结三场耦合数值分析 用自行开发的有限元程序为平台,加入本文推导出的本构模型,对实际工程应用中冻土的水分场、温度场及应力场之间的相互耦合问题进行数值模拟,与工程测量和前人所做的实验结果进行对比分析.3.1　控制微分方程 正冻结过程中渠道基土的温度分布,可表述为下列二维非线性变系数抛物型偏微分方程:C ・9T 9t =99x λ・9T 9x +99yλ・9T 9y (10)式中:C =C +L ・ρs ・9W u9T(11) λ=λ+L ・ρs ・KT a(12) 未冻水含量与温度的关系由室内试验确定:W u =f (T )C =C+12・(C ++C -)C- T >T pT b ≤T ≤T p T <T b (14)λ=λ+12・(λ++λ-)λ- T >T pT b ≤T ≤T p T <T b (15)K =K 3K 2・exp (A ・T )K 1 T >T pT b ≤T ≤T p T <T b(16)439 冰 川 冻 土 27卷　 综合导水系数K在剧烈相变区内随温度呈指数衰减.其初始条件和边界条件为: T(x,y,0)=φ(x,y)(17) T(x,0,t)=F1(t)(18) T(x,L1,t)=F2(t)(19) 9T(0,y,t)9x=0(20) 9T(L2,y,t)9x=0(21) 上述问题为一个温度、水分影响下以稳态温度场为初始温度条件的瞬态热传导问题,如果再加上外力边界条件和应力场的各种方程,特别是用我们提出的本构模型作为物理方程,则可以进行温度、水分影响下的应力计算. 渠道静力平衡方程为:[L]{σ}=0(22)几何方程为:{ε}=[L]{u}(23)其中:[L]=99x0099y099z99y099x99z00099z099y99x(24)物理方程(本构方程)为:{σ}=[D]{ε}即:σx=E1-v[εx-v(εy+εz)]σy=E1-v[εy-v(εz+εx)]σz=E1-v[εz-v(εx+εy)]τxy=G・γxyτyz=G・γyzτzx=G・γzx(25) 上述方程即组成了描述冻土渠道内含有相变的耦合方程组.3.2　定解条件和各参数的取值 根据文献[13,14],在计算中我们取冰的体积含量c i为0.05,粘土的弹性模量E s为12.0M Pa,冰的弹性模量取为4910.0M Pa,粘土的泊松比和冰的泊松比均取为0.3. 整个平面上的初始温度为0.3℃.上、下边界温度变化取为[14,15]:F1(t)=-11.2sin(π・t/3120.0) (d)F2(t)=10.9e-0.014(t+336.0) (d)(26) 计算模型以不良工程地质条件下渠道基土为例进行数学描述[14].假定渠道无限长,垂直截一断面为平面问题,认为土质均匀、水热蒸发耗热量及其它势场忽略不计.整个渠道模型宽为3.0m,高为2.0m,顶端渠堤宽为1.0m,渠底宽0.5m,渠底与下底边距离为1.0m.模型的位移边界条件为:下边界y方向限制移动、x方向自由,左、右边界x 方向限制移动、y方向自由,上边界的x和y方向均自由. 应力边界条件为:左右边界剪应力均为0.模型上边界受相当于渠道衬切板重量和与衬切板垂直的冻结力的均布压力P=1400N,方向垂直于边界线,如图3所示.图3　本文计算时外载荷施加位置Fig.3　Location of the appliedexternal load,calculated图4　文献[9]的外载荷施加位置Fig.4　Location of the applied externalload,From Reference[9] 各相关参数的取值如下[14,16]:C+=3076.7, C-=2279.0　(K J・m-3・℃-1);λ+=5.53,λ-=5396期刘雄等:冻土地区水渠的温度场和应力场数值分析　6.58　(K J ・m -3h -1・℃-1);K 1=1.8E -0.7　(m・h -1);K 2=0.0036・exp (0.551・(T +0.30))　m・h -1;K 3=3.6E -3　(m ・h -1);L =334.7　(K J・kg -1);T p =-0.30,T b =-0.75　(℃);W u =0.0664・exp (0.0551・T )0.3058+0.596(T +0.30)0.3058T <-0.75-0.75≤T ≤-0.30T >-0.303.3　计算结果分析与讨论 对整个模型采用三角形网格划分,共分成12484个三角形网格单元和6463个节点.为了对整个渠道冻结过程中的温度、位移、应力、应变进行仔细的观察,我们取了较长的时间来进行计算,时间单位为h ,总共计算时间为2000h ,时间步长取为5h.3.3.1　温度场分析 在渠道截面上,不同冻结时间的0℃等温线分布计算结果如图5所示,将其与图6计算的结果比较可以看出,0℃等温线分布相似.图5计算得到的0℃等温线随时间的变化Fig.5　The 0℃isotherm changingwith time ,calculated图6　文献[9]0℃等温线随时间的变化Fig.6The 0℃isotherm changing with time ,f rom Reference [9] 冻结2000h 后,上边界的负温达到最大值-10.12℃(如图7所示),这与文献[9]计算得到的冻结1536h 后上边界负温最大值-10.7℃基本一致.从图7中还可以看出,从渠坡到渠底的温度变化快,也即温度梯度大,而在离渠坡较远的渠堤下,温度分布呈一组几乎平行的直线,基本上不受渠坡边界温度的影响.从温度场的分布和变化来看,二维边界、截面几何形状、底水补给以及地中热流对温度分布均有显著影响.从计算结果与前人的计算结果比较来看,温度场的计算结果是令人满意的.图7　冻结2000h 时的等温线分布(℃)Fig.7　The distribution of isotherms (℃)after the soil frozen for 2000h3.3.2　位移场分析 大量的工程实践表明,在低冻结速率和有充足水分补给的条件下,产生较强的冻胀作用.在渠道底部和渠坡下部,温度较高的地下水起到降低冻结速率和补给冻结锋面水分的双重作用,在渠底和渠坡下部发生较强烈的冻胀.图8是2000h 时的冻胀位移(x 2+y 2)分布,可以看出,最大冻胀位移发生在渠堤顶端和渠坡上部,渠底和渠坡下部位移量也比较大,这与同时段的观测资料(图9)对比,可见位移的渠道截面上的分布规律是基本一致的.图8　冻结2000h 时的等位移线分布Fig.8　The distribution of displacement isolines after the soil f rozen for 2000h639 冰 川 冻 土 27卷　图9　宁夏大新渠试验段实测表面位移分布Fig.9　The observed surface displacements of the studied section of Ningxia Daxin Channel3.3.3　应力场分析图10～13是渠道冻结2000h 后的水平方向应力σx 和水平方向应变εx 、垂直方向应力σy 和垂直方向应变εy 的等值线分布图.从应力图上可以看出,在渠底和渠坡下部的冻胀力都比较大,在渠坡与渠底的交界点处出现了应力集中,最大冻胀力也在此出现:最大σx 约为3.943×105Pa ,最大σy 约为0.85×105Pa ,最大τxy 约为1.376×105Pa.图10　冻结2000h 时等бx 线分布Fig.10　The distribution of бx isolinesafter the soil f rozen for 2000h图11　冻结2000h 时等εx 线分布Fig.11　The distribution of εx isolinesafter the soil f rozen for 2000h 从渠道内的应力、应变分布图看出,应力和应变存在着一一对应关系,等应力线和等应变线的变化趋势相同,说明所用的弹性计算模型是正确的.图12　冻结2000h 时等бy 线分布Fig.12　The distribution of бy isolinesafter the soil frozen for 2000h图13　冻结2000h 时等εy 线分布Fig.13　The distribution of εy isolinesafter the soil frozen for 2000h4　结论 (1)建立了处于某一特定温度情况时土、冰的弹性常数与整个冻土体的弹性常数之间的关系,并在此基础上,将损伤力学的方法应用于冻土力学研究中,建立了考虑土体损伤的冻土本构模型. (2)从温度场的分布和变化来看,二维边界、截面几何形状、底水补给以及地中热流对温度分布均有显著影响.在渠坡和渠底的拐角处出现了应力集中,靠近这个区域点的应力变化比较剧烈. (3)冻胀应力沿渠道基土深度方向的分布状态,在渠道基土上部渠堤处由于外载荷的约束作用,基土中的冻胀应力呈衰减趋势.随着深度的增加,外载引起的变形逐渐减弱而由冻胀引起的变形逐渐增强,在超过某一深度后基土中的应力又会随着深度增加而增大,从而使基土应力随深度的增加呈现衰减—增大—衰减的规律.参考文献(R eferences):[1]　Ma Wei.Review and prospect of t he studies of ground freezingtechnology in China [J ].Journal of G laciology and Geocryolo 2gy (冰川冻土),23(3):218-224.[2]　Zhang Lijie ,Jiang Hao ,Li Lei.Study of calculation of soil7396期刘雄等:冻土地区水渠的温度场和应力场数值分析heat conduction:Progress and Prospect[J].Journal of G laci2 ology and Geocryology,2004,26(5):569-575.[张立杰,江灏,李磊.土壤中热量传输计算的研究进展与展望[J].冰川冻土,2004,26(5):569-575.][3]　Foriero A,Ladaryi B,Dallimore S R,et al Modelling of deepseated hill slope creep in permafrost[J].Can.Geotech.J., 1998,35(4):560-578.[4]　K lein J.Finite element met hod for time-dependent problemsof frozen soils[J].Int.J.for Numerical and Analytical Met h2 ods in Geomechanics,1981,5:263-283.[5]　Rai Mand,Roegiers J C.Fluid flow and heat flow in deform2able fractured porous media[J].International Journal of Engi2 neering Science,1994,32:1615-1633.[6]　Masters I,Pao W K S,Lewis R W.Coupling temperature to adouble-porosity model of deformable porous media[J].In2 ternational Journal for Numerical Met hods in Engineering, 2000,49:421-438.[7]　Lai Yuanming,Wu Ziwang,Zhu Yuanlin.Nonlinear analysisfor t he coupled problem of temperature and seepage fields in cold region tunnels[J].Cold Regions Science and Technolo2 gy,1999,29:89-96.[8]　Lai Yuanming.Numerical simulation for t he coupled problemof temperature and seepage fields in cold region dams[J].Journal of G laciology and Geocryology(冰川冻土),2000.22 (Suppl.):171-176.[9]　He Ping,Cheng Guodong,Yu Qihao,et al.A couple modelof heat,water and stress fields of saturated soil during freezing [J].Journal of G laciology and Geocryology,2000,22(2):135 -138.[何平,程国栋,俞祁浩,等.饱和正冻土中的水、热、力场耦合模型[J].冰川冻土,2000,22(2):135-138.][10]Chen Feixiong,Li Ning,Cheng Guodong.The t heoreticalframe of multi2phase porous medium for t he freezing soil[J].Chinese Journal of Geotechnical Engineering,2002,24(2): 213-217.[陈飞熊,李宁,程国栋.饱和正冻土多孔多相介质的理论构架[J].岩土工程学报,2002,24(2):213-217.] [11]Li Hongsheng,Du Xiaozhen.Constitutive model of damageand fracture of ice[J].Journal of G laciology and Geocryology, 2003,25(Suppl.):304-307.[李洪升,杜小振.冰体材料损伤与断裂破坏的本构理论[J].冰川冻土,2003,25(Suppl.): 304-307.][12]Wu Zheng.Study of concrete’s constitute model based on dam2age mechanics[J].Technology of Water Conservancy and Wa2 ter and Electricity,1995,11:58-63.[吴政.基于损伤的混凝土拉压全过程本构模型研究[J].水利水电技术,1995,11: 58-63.][13]Wu Ziwang,Ma Wei.Creep and Strengt h of Frozen Soil[M].Lanzhou:Lanzhou University Press,1994.48-50.[吴紫汪,马巍.冻土强度与蠕变[M].兰州:兰州大学出版社,1994.48 -50.][14]An Weidong,Wu Ziwang,Ma Wei,et al.Interaction amongTemperature,Moisture and Stress Fields in Frozen Soil[M].Lanzhou:Lanzhou University Press,1989.184-211.[安维东,吴紫汪,马巍,等.冻土的温度、水分应力及其相互作用[M].兰州:兰州大学出版社,1989.184-211.][15]Deng Y ousheng,He Ping,Zhou Chenglin.An experimentalresearch on t he t hermal conductivity coefficient of saline soil [J].Journal of G laciology and Geocryology,2004,26(3): 319-323.[邓友生,何平,周成林.含盐土导热系数的试验研究[J].冰川冻土,2004,26(3):319-323.][16]Chen Ji,Cheng Guodong,Zhang Xifa.Study on t he observa2tion of delamination frost2heave amount[J].Journal of G laci2 ology and Geocryology,2004,26(4):466-473.[陈继,程国栋,张喜发.分层冻胀量的观测方法研究[J].冰川冻土, 2004,26(4):466-473.]Numerical Analyses of the T emperature and StressFields of Channel in Frozen Soil R egionsL IU Xiong1,2,　N IN G Jian2guo1,　MA Wei3(1.S tate Key L aboratory of Ex plosion Science and Technolog y,B U T,B ei j ing100081,China;2.Changsha Uni versit yof Science and Technolog y,Changsha H unan410076,China;3.S tate Key L aboratory of FrozenS oil Engi neeri ng,CA R E ER I,CA S,L anz hou Gansu730000,China)Abstract:In t his paper,an elasticity constit utive model of damaged f rozen soil is established.A cal2 culation is made for f rozen sand soil under different ice content and temperat ure conditions,using dam2 age constit utive model,which agrees wit h t he curves of measured st ress2st rain.Numerical com2 p utations are done by using t he p rocedure devel2 oped in t his st udy,and t he comp utation result s of moist ure,temperat ure and st ress fields of t he channel are analyzed.It is found t hat t he result s are more exact and agree wit h act ual temperat ure field,st ress field,moist ure field and st rain field, and is consistent wit h t he result s of predecessors and t he test result s.The procedure developed in t his st udy could calculate t he related p hysical prop2 erties of frozen soil.In addition,it could describe t he relation between t he related p hysical quantities well,and p rovide a necessary reference for engi2 neering design and const ruction.K ey w ords:f rozen soil;damage;coupling;temperat ure field;stress field839 冰 川 冻 土 27卷　。