考虑土中水分迁移路堑边坡冻融失稳数值分析

ClassifiedIndex:TU445U.D.C.:625DissertationfortheMasterDegreeinEngineeringNUMERICALANALYSISOFCUTTINGSLOPESTABILITYWITHMOISTUREMIGRATIONDURINGFREEZE-THAWPROCESSCandidate:XuHaoSupervisor:ProfessorTaoXiaxinAcademicDegreeAppliedfor:MasterofEngineeringSpecialty:StructrualEngineeringAffiliation:Dept.ofCivilEngineeringDateofDefence:June,2009Degree-Conferring-Institution:HarbinInstituteofTechnology 哈尔滨工业大学工学硕士学位论文摘要在季节冻土区，经常在开春的融化季节发生路堑土质边坡的浅层失稳滑塌，轻者造成修复的工程量、经费投入增加，重者造成对正常通行的影响。这种冻融失稳与冻融过程中的水分迁移有密切的关系。为了深入分析路堑土质边坡冻融失稳过程中水分迁移的影响，本文结合作者导师课题组最近发展的冻融过程中水分迁移模型，估计在秋天、冬天、春天的含水量分布、浅部土层抗剪强度参数的变化，借助大型软件ABAQUS分析了一个10米高，初始含水量30%的土坡内力、变形和塑形应变的分布。在几种边坡稳定性计算方法、岩土材料中土的各种屈服准则的归纳基础上，系统阐述了水分迁移的研究进展和本文引用的水分迁移模型，有限单元的基本原理和计算基本步骤。据当地的地表气温统计资料，以10月10日为冻结开始的日期，3月10日为融化开始的日期，即冻结后150天。根据地区经验数据，最大冻结深度为1.6米，对照冻融试验中获取的不同深度的温度随时间的变化曲线，本文建立了冻结深度随时间的变化关系、温度梯度与距冻结锋面距离的关系、融化土层底面随时间的关系、温度梯度与据冻融界面距离的关系。用修正的Fick定律来描述冻结过程中的水分向上迁移，增加了温度梯度的影响项，按每10天为一个步距，计算冻结后含水量的变化。融化中水分向下渗流用修正的非饱和土达西定律描述，也增加了温度梯度的影响项，亦按每10天为一个步距，计算融化后含水量的变化。结果显示，在冻融界面上的土层中含水量可以增加20%。显然，含水量的变化会导致相应土层抗剪强度的变化。本文根据收集到的有限数据，建立了含水量与抗剪强度参数内聚力和内摩擦角的关系。据此，按含水量的变化一步步修改了土层抗剪强度的参数。借助ABAQUS软件，采用M-C模型考虑土的弹塑性，计算了秋天（开始冻结时）、冬天（冻结达到1.6米）和春天（融化达到0.8米）三个时刻的边坡内力、变形。结果显示了由于含水量的变化引起的最大塑性应变的变化。关键词：边坡稳定；有限元；水分迁移；冻融–I– 哈尔滨工业大学工学硕士学位论文AbstractInseasonfrozenregions,thereareshallowslidesonsoilcuttingslopesalonghigh-waysquiteoftenoccurredinspring.Theslightofthemusuallycausealotrepaircost,labourtimeandmaterials,theseveremayinfluencethetraffic,eveninterruptit.Theseshallowslidearecloselyrelatedtothewatermigrationinsoilduringthefreeze-thawpro-cess.Inordertodeeplydealwiththeeffectofwatermigrationduringthefreeze-thawprocessontheslideofsoilcuttingslope,watercontentdistributionsinautumn,winterandspringinsoilslope,variationofshearstrengthparametersinshllowsoilareestimatedbymeansofthewatermigrationmodelinsoilduringfreeze-thawprocessthatdevelopedrecentlybytheteamoftheauthor’ssupervisorinthispaper.Thestresses,displacementsandplasticstrainsinasoilcuttingslopewith10mhightandinitialwatercontent30%arecalculatedbyABAQUSsoftware.Basedonasummaryonstabilityanalysisofsoilslopeandsomeyieldingcriterionsofgeotechnicalmaterials,thewatermigrationmodeladoptedinthispaperandtherelativeresearchachievements,basicprinciplesandmainstepsofthefiniteelementmethodthenareintroducedsystematically.Accordingtothestatisticresultoflocalclimaterecords,October10thisseletedasthedaytostartsoilfreezing,March10thisthedaytothaw,whichis150daysafterfreezestarted.Thelocalmaximumfrozendepthisabout1.6maccordtolocalempiricaldatas.Comparedwiththetemperaturecurvesatsomedifferentshallowdepthswithtimefromthefreeze-thawtest,relationshipsoffrozendepthwithtime,oftemperaturegradientwithdistancetothefreezingfront,ofthaweddepthwithtime,oftemperaturegradientwithdistancetothethawed-frozensurface,areestablished.ThemoistureinsoilmoveupwardduringfreezeprocessisdescribedbyanimprovedFicklawinwhichtheeffectoftemperaturegradientisadded.Thevariationofwatercon-tentsiscalculatedbyevery10days.ThewaterpermeatedownwardduringthawprocessisdescribedbyanimprovedDarcylawinwhichatemperaturegradienttermisalsoadded.Thevariationofwatercontentsiscalculated.Theresultshowsthatthewatercontentinthesoillayerjustabovethefrozen-thawedinterfaceincreased20percentsobviously.–II– 哈尔滨工业大学工学硕士学位论文Ofcource,thechangeofwatercontentofsoilmustcauseamodificationofsoilshearstrength.Arelationofwatercontentwithshearstrengthparameters,isestablishedfromsomecollecteddata.Thesoilshearstrengthismodifiedstepbystepfollowingthewatercontent.Thestresses,displacementsandplasticstrainsinthesoilcuttingslopeinautumn(freezingstarts),winter(frozento1.6m)andspring(thawedto0.8m)arecalculatedbymeansofABAQUSsiftwareandM-Cmodelforsoilelasto-plasiccharacristics.Theresultshowsachangeofmaximumplasicstrainintheslopefromthewatermigrationduringthefreeze-thawprocess.Keywords:stabilityofslope;FEM;moisturemigration;freeze-thaw–III– 哈尔滨工业大学工学硕士学位论文目录摘要..........................................................................IAbstract........................................................................II第1章绪论...................................................................11.1课题背景及研究的目的和意义..........................................11.2路堑边坡冻融失稳研究的发展概况.....................................21.2.1边坡稳定性理论的发展.............................................21.2.2冻土研究的发展....................................................31.2.3冻土边坡稳定性理论研究...........................................41.3本文的主要研究内容及章节安排........................................6第2章粘性土质边坡失稳机理...............................................82.1引言.....................................................................82.2边坡稳定性计算方法....................................................82.2.1土坡稳定的一些传统分析方法......................................82.2.2有限元法............................................................102.3土的屈服准则............................................................102.3.1两个古典屈服准则..................................................102.3.2广义Mises准则（Drucker-Prager准则）............................152.3.3Mohr-Coulomb准则.................................................162.3.4Mohr-Coulomb等面积圆屈服准则..................................192.4本章小结.................................................................20第3章冻融土中水分冻融迁移模型..........................................213.1引言.....................................................................213.2土冻结过程中的水分迁移...............................................213.2.1土冻结过程中的水分重分布........................................213.2.2土冻结中水分迁移微观机理的一种解释............................223.3土冻结的水分迁移模型..................................................233.3.1模型的物理过程描述...............................................233.3.2模型的数学表达....................................................243.4土融化过程中的水分迁移...............................................263.4.1土融化过程中的水分重分布及物理过程描述.......................263.4.2模型的数学表达....................................................273.5本章小结.................................................................30–IV– 哈尔滨工业大学工学硕士学位论文第4章边坡冻融数值模拟和分析.............................................314.1引言.....................................................................314.2有限单元法基本原理....................................................314.2.1有限单元法求解的基本步骤........................................314.2.2计算区域的选择....................................................324.2.3边界条件的设置....................................................324.2.4单元的选择.........................................................334.3软件介绍.................................................................334.4模型建立.................................................................344.4.1假设条件............................................................344.4.2计算思路............................................................344.4.3模型的建立和力学参数的选取......................................354.5本章小结.................................................................36第5章冻融过程中边坡稳定性的计算分析...................................385.1引言.....................................................................385.2冻结深度和冻结过程中边坡土中含水量的变化.........................385.3融化深度和融化过程中边坡土中含水量的变化.........................415.4冻融过程中边坡表层土的抗剪强度参数的变化.........................445.5冻融过程中三个时刻边坡的内力和变形的分布.........................455.6本章小结.................................................................45结论与展望.....................................................................47参考文献.......................................................................49攻读硕士学位期间发表的学术论文及其它成果................................53જᇨⒼ工Ϯ໻学硕士学位论文原߯性ໄᯢ.....................................54જᇨⒼ工Ϯ໻学硕士学位论文Փ用授ᴗк.....................................54致谢..........................................................................55个人简历.......................................................................56–V– 哈尔滨工业大学工学硕士学位论文第1章绪论1.1课题背景及研究的目的和意义在我国幅员广大的季节性冻土地区，进行着许许多多水利工程、工业民用建筑及交通运输工程的建设。在这些工程建设和使用中，势必遇到大量的季节性冻土地区特有的工程技术问题。工业民用建筑方面，地基土的冻胀及融化引起房屋裂缝、倾斜及变形的现象在冻土地区早已是屡见不鲜了；农田水利建设方面，我国季节冻结深度较大的东北、华北及西北地区，许多工程设施由于冻害问题不能发挥良好的应用；交通运输工程方面，路基的冻胀、融沉、翻浆及路堤与路堑边坡的冻融滑坡，影响了公路交通的正常营运与行车安全[1]。目前，我国正处于经济建设的高速发展时期，公路交通正在人民的生产生活中占据着越来越重要的地位。而一直以来，公路工程的边坡稳定性问题直接关系到建设的资金投入和人民的生命财产安全。黑龙江省季节性冻土分布地域广，路堑边坡破坏原因复杂，冻融失稳引起的问题比较突出。一些公路的路段，在边坡开挖、砌筑完工后是稳定的，秋季雨水较多就会在冬季引发土体冻胀使浆砌块石等覆面胀裂，春季随着冻土的融化，路堑就可能出现失稳下滑。路堑边坡失稳影响了公路的正常营运，增加维护的费用，严重时甚至导致交通中断。例如，1997年哈同高速公路（哈尔滨至同江）通车次年春融期出现大量滑坡，治理费6000余万元；2002年黑北公路（黑河至大连公路的黑河市至北安路段）路基施工后次年春融期发生5km整体滑坡，被迫改线，直接经济损失达7000余万元，严重影响了工程的效益，直接延缓了当地的经济发展[2–4]。路堑边坡冻融失稳已经引起了国内外寒冷地区公路工程、岩土工程界的关注，相关单位相继投入人力、财力展开科学研究，也得到了一些共识。路基边坡在冻融过程中的失稳与土中水分在冻结过程中的迁移有密切关系。冬季温度降低，冻结过程中土中的水分向上迁移，使上层土体中水分富集，形成冰晶或冰夹层。冬季路基边坡又被雪所覆盖，进入春季，坡面土层由表及里逐渐融化，冻结面逐步下降，土中的水分（主要来源为土体中冰的融化和坡面及其上方融化的积雪）逐渐向下渗透。由于冻结土层不透水，上层土中渗透下来的水分在冻结面上方积聚，在一个不厚的土层中达到饱和或接近饱和，导致水分聚集层土体的抗剪强度及摩擦力下降。未融化的冻结土层表面形成倾斜滑动–1– 哈尔滨工业大学工学硕士学位论文面，表面摩擦系数小，在重力作用下会发生浅层滑塌破坏。边坡的冻融失稳原因除了与边坡土质条件有关外，还要受地下水、气候、边坡坡度、边坡朝向、春融期冻结层等各方面因素的影响。作者导师领导的课题组持续开展了一系列的研究，最近通过试验研究发展了冻融过程中水分迁移的模型。本文拟在此基础上，借助大型商业软件ABAQUS建立边坡有限元分析模型，计算从秋至冬到春路堑边坡土层冻结–融化过程中含水量的变化，根据含水量调整相应土层的抗剪强度，计算边坡的内力和变形，评价冻融稳定性，为季节性冻土区土质边坡的设计及浅层冻融失稳防治提供参考。1.2路堑边坡冻融失稳研究的发展概况1.2.1边坡稳定性理论的发展边坡稳定性研究由来已久，早期的边坡研究是仅以土体为主要研究对象的，是采用材料力学和简单的均质弹性、弹塑性理论为基础的半经验半理论的研究。由于力学机理的粗浅或假设的不合理，计算结果往往与实际情况差别较大。1959年的法国Malpasset坝左岸坝肩岩体的崩塌及1963年的意大利Vajont坝上游左岸的库岸边坡滑坡等，使人们清醒地认识到了边坡破坏的力学机理研究需要大改进。以弹塑性理论为基础和改进的极限平衡法应用为主的多种方法应运而生。1967年，人们第一次尝试用有限元研究边坡稳定性问题，给定量评价边坡的稳定性创造了条件，发挥数值方法的作用，使边坡稳定性研究进入模式机制和作用过程研究成为可能，随着大量规模巨大工程的开展，我国边坡工程的不断实践与发展，边坡稳定性分析的原理、方法也获得了不断丰富与发展。进入20世纪80年代后，各种复杂的数值计算方法广泛地应用于边坡研究，学科之间的相互渗透，使研究方法各式各样，成果叠出。现今使用的数学模型大致可分为两类：一类是基于极限平衡理论的条分法；另一类是数值分析方法。条分法人为假定的条间作用力并不代表斜坡真实的应力状态，但计算简单，发展历史较长，就均质土坡稳定性而言，结果已可满足实际需要。进入20世纪90年代，边坡问题的研究将传统的边坡工程地质学、现代岩土力学和现代数学相结合，形成了现代边坡工程学；各种现代科学的新技术不断用于边坡问题研究中、从而给边坡的稳定性研究提供了新理论、新方法[5]。–2– 哈尔滨工业大学工学硕士学位论文1.2.2冻土研究的发展冻土，一般是指温度在0◦C或0◦C以下、含有冰的各种岩土。按岩土冻结状态保持时间的长短，冻土又可以分为多年冻土与季节冻土，冻结数年至数万年以上的称为多年冻土。地球上冻土区面积约占陆地总面积70%，多年冻土分布面积35106km2，占陆地面积的25%，主要分布在北半球，包括欧亚大陆的西伯利亚和北美大陆的阿拉斯加及加拿大广阔地区的多年冻土，约占全球多年冻土总面积的63%。我国多年冻土分布面积约21.5105km2，位居世界第三，包括位于欧亚大陆高纬度多年冻土区南缘的东北大、小兴安岭（分布南界达北纬45◦C46◦C），以及西部高山和青藏高原等地。季节性冻土分布面积5.14106km2，约占我国国土面积的53.5%。人类对冻土的认识伴随着工农业生产实践，经历了不断发展的进程。国际上对冻土工程问题最早开展研究进而建立起独立学科的当数俄罗斯[6]。罗蒙诺索夫在1757年发表《冻土地》的科学综述。19世纪后半叶，西伯利亚大铁路的建设极大地推动了冻土研究。19041914年间修建的阿穆尔铁路常因冻融影响而破坏，每年要耗费大量的资金维护。科学家对冻土工程地质做了进一步研究，并在后来的贝加尔–阿穆尔铁路干线（19791984年建设）修建中完成了大量的冻土研究工作。1917年后，前苏联在冻土力学、冻土温度场计算等理论分析方面独到的见解，对全球冻土研究有很大的借鉴作用。由于寒区工程的需要，20世纪30年代初，前苏联对冻胀融沉机理进行了研究，20世纪50年代初，北美开始进行冻土强度理论的研究，20世纪60年代，美国以Hass为代表开始了现场冻胀的观测工作。在冻胀方面的研究主要集中于冻胀与冰透镜体的形成和位置关系、冻胀与冻渗及冻胀过程的数值模拟方面。Miller（1972）提出正冻粘土的冰透镜体底面与冻结锋面间存在着冻结缘，Penner和Ueda（1977）的试验结果也支持该观点。Williams和Wood（1985）通过试验测试了正冻土中的内应力变化。由于缺乏适合的仪器设备和方法，目前通过试验观测冻土中应力的变化仍相当困难。在冻土的强度方面，Butkovich（1954）、Mellor和Smith（1966）指出冻土的力学性质主要取决于胶结冰的存在，而胶结多晶冰的强度受多种因素的影响，主要有温度、压力、应变速率及冰晶的大小、结构和方向。冰的强度随温度的降低而增加，并随冰晶大小、方向、密度和结构的变化而变化。Sayles和Epanchin（1966）、Zhu（1988）等人发现随着应变速率的增加，冰的强度也逐–3– 哈尔滨工业大学工学硕士学位论文渐增大，并且表现在破坏类型上由塑性到脆性的变化；当冻土发生脆性变化时，强度与应变速率或加荷速率关系较小[2]。在冻土的变形的结构效应上，Miller（1960）、Sayles（1973）、Vrachev（1984）指出固体矿物客体间的摩擦力和粘聚力是影响冻土长期强度的主要因素。Ladanyi（1985）、马巍（1993）发现当冻土体受到小于破坏强度的常荷载作用时，冻土体逐渐发生变形，由冰所承担的一部分应力逐渐转移到土骨架上，直至土骨架达到它的长期强度极限，一旦超过此极限，冰在冻土中起到一种粘聚力的作用，土体将继续变形直至破坏。Chamberlain等（1972）、Simonsen等（1974）、Jones和Parameswarm（1983）、Fish（1991）等还指出，当冻土受围压作用时，首先在颗粒接触处的冰发生融化，然后向低应力处迁移；当围压超过某一界限，土体结构发生破坏，孔隙水产生整体压融，从而导致冻土强度降低。1973年，Harlan提出了土体冻结过程中水–热迁移耦合数学模型，从此进入了多场耦合问题的研究阶段，工作核心是基于冻土中的热质迁移来进行数值模拟。随着对冻胀现象的研究，在伴有相变过程的冻土三场耦合问题方面，提出的冻胀预报模型有经验型、半经验型、流体动力学模型、刚性模型、热力学模型。苗天德等（1999）在连续统力学混合物理论框架下研究了冻土力学–热学性质，建立起固、液两相介质伴有相变的水、热二场耦合模型。该模型属于非线性的Burgers方程，相对于经典的Stefan线性热传导方程，可以描述冻结过程中的水热耦合效应。李宁、陈飞熊等（2000）在冻土中骨架、冰、水、气体多孔多相介质的基础上，建立了冻土多孔多相微元体的平衡方程、多孔固液介质的质量守恒方程及多孔多相介质的热、能守恒方程，进一步分析冻土的冻胀与融沉现象[2]。1.2.3冻土边坡稳定性理论研究冻土边坡冻融稳定性分析是介于边坡工程和冻土力学的一个交叉研究领域。1972年前苏联学者叶米里扬诺娃在《滑坡作用的基本规律》一书中较详尽地刻画了边坡融冻泥流现象[7]。国内研究开始较晚，1980年周幼吾提出了冰缘地貌中最值得注意的热融滑塌问题[8]，1989年王文宝研究了冻土地区路基边坡滑塌与整治问题[9]，1990年王绍令研究了冻土边坡热融滑塌问题[10]，1993年张长庆、朱林楠等研究了冻土地区边坡开挖与热融稳定性的问题[11]，郭东信研–4– 哈尔滨工业大学工学硕士学位论文究了融冻泥流的阶地问题[12]。冻土边坡稳定性评价方法发展至今，较为典型的有三种：Weeks和Chandler两人推导的基于“冰阻渗流”导致孔隙水压力增加机理的有效应力分析法[13]；由Hutchinson提出总应力分析法，认为由于冻结锋面冰集聚导致融化时土体含水量增加，土的不排水抗剪强度降低，从而引发斜坡失稳[14]；由McRoberts和Morgenstern提出，基于有效应力和融化–固结理论，认为在斜坡冻土融化–固结过程中，易滑面上超孔隙水压力增加并可引发斜坡失稳[15]。我国冻土滑坡研究起步较晚，《中国冻土》一书，曾对冻土斜坡失稳现象进行了划分和描述。近年来，随着我国西部大开发、青藏铁路的建设，中国科学院寒区旱区环境与工程研究所、兰州大学及长安大学的一些科技工作者对高原多年冻土斜坡的稳定性进行了详细的现场观测、模型试验研究及理论分析，取得了可喜的成果。刘红军、王丕祥（2006）根据实际现场调查情况，论述了土质边坡冻融失稳的主要原因及影响边坡冻融失稳的主要因素。考虑了地下水的渗流作用,建立了冻土边坡冻融失稳稳定性分析几何模型。采用有效应力法，分别把土和水的混合体当作研究对象和把土骨架作为研究对象，得出两种方法在求季冻区土质边坡冻融失稳的稳定性系数时，会得到一致的结果[16]。李东庆、魏春玲、吴紫汪（2000）利用数值模拟方法,分析了在多年冻土地区修筑路堤后,夏季路基边坡积水通过路堤及在其下土壤的渗流过程中,由于对流换热所引起的路堤基底的不稳定现象。通过计算机数值模拟进一步证实,夏季路基边坡积水可导致路基基底之下的多年冻土上限下降,并使路基下形成凹型融化核或扩大已有融化核的深度和范围,从而加大路堤的融化下沉量和冬季的冻胀量,造成路基稳定性程度明显下降[17]。刘兵（2008）建立了土中水分在冻融过程中的迁移机理研究及模型[18]。孟燕燕等人（2007）建立了路堑边坡的有限元分析模型，用大型有限元软件ABAQUS研究路堑边坡冻融稳定性的判别方法，综合考虑在冻融过程中不同深度土的力学性质的变化，易滑面上超孔隙水压力增加致使其上部土层抗剪强度明显降低的影响，据此分析边坡的稳定性，进一步评价了不同工程防治措施的效果。然而，已发表的研究成果均没有能够模拟水分迁移过程对土质边坡失稳的影响，这无疑将成为我国季节性冻土区道路工程中急待研究并解决的突出问题。–5– 哈尔滨工业大学工学硕士学位论文1.3本文的主要研究内容及章节安排本文拟在边坡冻融失稳的有限元数值模拟中，把作者导师和刘兵（2008）新发展的水分迁移模型结合进来，形成一套新的分析程序。具体工作如下：(1)完善冻结和融化过程中土中水分的迁移模型。在分析过程中区分冻结和融化两个过程，应用不同的水分迁移模型。在冻结过程中，温度降低使地表冻结层逐步向深处发展，汽态水分接近冻结层就会凝结成液态水、冻结成冰晶，；邻近冻结层的气态水冻结后，土中孔隙蒸汽压力梯度变化，使下层的汽态水向上扩散，补充上层水分，导致上层水分再冻结，下层水分再扩散补充，造成了冻结过程的水分迁移。用扩散模型来描述土中的水分在冻结过程中的向上迁移。融化过程中，冰融化成水，在重力的作用下，非饱和土中液态水向下渗流，应用渗流模型来描述土中的水分的向下迁移。(2)借助ABAQUS的二次开发，完善边坡稳定性计算的有限元程序，根据含水量的变化替换土层的抗剪强度参数，计算边坡土体内力和变形，判别路堑边坡冻融稳定性。建立地表层单元的抗剪强度参数与含水量的统计关系，根据水分迁移模型顺序计算冻结、融化过程中每一步各个深度的含水量。根据得到的含水量，修正土体的力学性质，主要是抗剪强度，一步步分析边坡的稳定性。模拟多种因素对稳定性的影响。本文的章节安排如下：第一章，阐述选题的背景及研究的目的和意义，总结了路堑边坡冻融失稳、冻土以及冻土边坡稳定性理论研究的发展，分析了原有理论的不足，并针对该问题归纳研究思路第二章，阐述了几种边坡稳定性计算方法，着重归纳了岩土材料中土的各种屈服准则，并对各准则进行了论述比较。第三章，主要叙述了水分迁移的研究进展和本文引用的水分迁移模型，并进行了完善。区分冻结和融化这两个过程，分别采用了修正的Fick定律和修正的达西定律作为水分迁移模型。第四章，主要介绍了有限单元的基本原理，以及其在边坡冻融数值模拟和分析中的应用，和其分析步骤，并在ABAQUS/CAE中建立边坡稳定性计算的有限元模型。在分析过程中区分冻结和融化两个过程，应用不同的水分迁移模型计算含水量变化，依据含水量的和粉质粘土的抗剪强度参数的数学关系得到抗–6– 哈尔滨工业大学工学硕士学位论文剪强度指标，在计算过程中不断更改，分析边坡稳定性。如此反复计算分析边坡的塑性发展，判断边坡是否失稳。第五章，主要根据哈尔滨市气温的统计资料，模拟出冻结和融化过程中融冻深度和过程中边坡土的含水量变化。将秋至冬到春划分成10天一个步距，计算冻结深度、温度梯度和相应的含水量分布，根据每一步距边坡剖面上含水量的变化，依含水量与抗剪强度参数的关系调整参数值，带入相应的边坡有限元模型计算，评价边坡的稳定性。–7– 哈尔滨工业大学工学硕士学位论文第2章粘性土质边坡失稳机理2.1引言边坡是指具有倾斜坡面的土体，通常可分为天然土坡和人工土坡。边坡上的部分土体在自然或人为因素的影响下沿某一滑动界面发生整体向坡下运动的现象称为滑坡。影响土坡滑动的因素复杂多变，从力学上讲，根本原因在于土体内部某个面上的剪应力达到、超过了抗剪强度，使稳定平衡遭到破坏。导致边坡滑动失稳的原因可有以下两种：(1)外界荷载作用或突破环境变化等导致土体内部剪应力加大，例如路堑或基坑的开挖、堤坝施工中上部填土荷重的增加、降雨导致土体饱和增加重度、土体内地下水的渗流力、坡顶荷载过量或由于地震、打桩等引起的动力荷载等。(2)外界各种因素影响导致土体抗剪强度降低，促使土坡失稳破坏，例如超孔隙水压力的产生，气候变化产生的干裂、冻融，粘土夹层因雨水等侵入而软化，以及粘性土蠕变导致的土体强度降低等。边坡稳定性是高速公路、铁路、机场、高层建筑深基坑开挖等土木工程建设中十分重要的问题。土边坡冻融稳定性问题可通过边坡稳定性来分析解决。2.2边坡稳定性计算方法2.2.1土坡稳定的一些传统分析方法土坡稳定分析的方法都是基于岩土塑性理论，土边坡的变形与发展使滑动面附近的土层处于塑性阶段，直至坡体失稳破坏。边坡失稳在力学上主要是一个强度问题，计算上可简化为一个静力平衡问题和一个岩土屈服条件，其解可满足工程的要求。传统的边坡稳定分析方法都循着这种思路，常用的方法有极限平衡法、滑移线场法、极限分析法等。(1)极限平衡法极限平衡法是当今国内外应用最广的边坡稳定分析方法，是传统边坡稳定分析方法的代表。在已知滑移面上对边坡进行静力平衡计算，从而求出边坡–8– 哈尔滨工业大学工学硕士学位论文稳定安全系数。必须事先知道滑移面的位置和形状，对于均质土体可以通过经验或者优化的方法获得滑移面，因而十分适用于土质边坡。当滑移面为一简单平面时，静力平衡计算可采用解析法计算，获得解析解。著名的库仑公式就是一例，一直沿用至今。当滑移面为一圆弧、对数螺线、折线或任意曲线时，无法获得解析解，通常要采用条分法求解。此时坡体为一静不定问题，通过对某些未知量作假定，使方程式的数目与未知数数目相等从而使问题成为静定，这种方法十分简便，而且计算结果能满足工程要求而被广泛应用。由于假设条件与应用的方程不同，条分法分为非严格条分法与严格条分法。在非严格条分法中，通常只满足一个平衡条件，而不管另一个平衡条件，计算结果会有一定误差。非严格条分法有两个未知数（安全系数和条件力的作用方向），但只有一个方程，还需作一个假定。非严格条分法通常是假定条间力的方向，由于假定不同而形成各种方法，有瑞典法、简化Bishop法、简化Janbu法、Sarmar(I)法等。严格条分法满足所有的力平衡条件，它有三个未知数（安全系数、条间作用力方向和作用点）和两个方程，也要做一个假定。随所采用的假设不同，有Morgenstern-Price法、Spencer法、Janbu法、Sarma(II)法和Correia法等。(2)滑移线场法滑移线场法严格满足塑性理论，假定土体为理想塑性体，并将土体分为塑性区与刚性区。塑性区满足静力平衡条件和莫尔–库仑准则。二者结合得一组偏微分方程，采用特征线法求解。然而，严格滑移线场的需要数值解，且也只能解稍微复杂的问题，对于复杂的问题滑移线场法常常无效，而且也只能用于均质土体。(3)极限分析法极限分析法是运用塑性力学中的上、下限定理求解边坡稳定问题。上限法也称能量法，通常需要假设一个滑裂面，并将土体分成若干块，土体视作刚塑性体，然后构筑一个协调位移场。为此，需要假设滑裂面为对数螺线或直线，然后根据虚功原理求解滑体处于极限状态时的极限荷载或稳定安全系数。极限分析下限法的理论基础是下限定理，它在计算过程中需要构造一个合适的静力许可的应力分布，在通常情况下可用应力柱法或者应力不连续法等来求得问题的下限解，偏于安全，可以使用。但只有极少数情况下可以获得下限解。目前，已将其扩展为上限有限元法与下限有限元法，不需假定滑面，从而扩大了应用范围。显然这种方法也适用于土体[19]。–9– 哈尔滨工业大学工学硕士学位论文2.2.2有限元法随着计算机技术、数值分析方法、多媒体技术和网络技术的发展，数值模拟得到了越来越广泛的应用，很多数值计算方法都用到边坡稳定分析中。其中，有限元法是一种十分成熟的方法，它几乎可适用于所有的计算领域。其最大优点是可采用确定的本构关系，分析任何形状的几何体，不但能进行线性分析还可以进行非线性分析，是边坡稳定分析中用得较多的一种数值方法。有限元用于边坡中应力和位移的分析已经比较完善，有许多成功应用的例子。将该方法得到的结果用于边坡稳定性评价却没有引起足够的重视。下面归纳采用有限元法确定边坡的应力和位移，评价边坡稳定性的几种方法。(1)刚体有限元法模型把结构离散成任意凸多边形刚体单元，这些刚性单元由分布在它们之间的弹塑性夹层连接，凸多边形单元中任一点的位移完全由其形心得刚体位移描述，它们可以模拟各种形状的岩块，而用弹塑性夹层模拟软弱结构面。以形心点处的物理量来建立刚体弹塑性夹层模型的有限元方程式，用有限元法进行边坡稳定性分析评价。(2)有限元强度折减法在有限元法中，通过强度折减，直至计算到不收敛为止，其折减的倍数即为稳定系数。这种方法国外在20世纪80年代就采用，其缺点是力学概念不十分明确，而且要受到计算程度及计算精度的限制。利用有限元法，通过强度折减来求边坡稳定安全系数。通过强度折减，使系统达到不稳定状态时，有限元计算将不收敛，此时的折减系数就是安全系数，安全系数的大小与所采用的屈服准则有关。综上所述，这种方法具有广阔的应用前景，计算准确、简便、适应性强，有可能引发岩土工程设计方法的重大改革[5]。2.3土的屈服准则有限元法中岩土材料本构模型采用理想弹塑性模型。选用合理岩土屈服准则十分重要。通常采用的岩土屈服准则是广义Mises准则和Mohr-Coulomb准则。2.3.1两个古典屈服准则(1)Tresca准则–10– 哈尔滨工业大学工学硕士学位论文当最大剪应力达到一定数值时，材料开始进入塑性状态，这是Tresca（1864）根据金属实验提出的屈服准则，称最大剪应力条件，在材料力学中称为第三强度理论，可写成kτmax=(2-1)2式中k———实验常数；τmax———三个主剪应力绝对值的最大者。三个主剪应力与三个主应力有关，如规定σ1σ2σ3，Tresca准则可表示为：σ1σ3kτmax==(2-2)22或σ1σ3=k(2-3)但在一般情况下σ1，σ2，σ3不按大小次序排列，则上式(2-3)可表示为：max(jσ1σ2j,jσ1σ2j,jσ1σ2j)=k(2-4)上面的等式如果能成立，则材料开始发生塑性变形。也可将上式中三个式合并写成222222F=[(σ1σ2)k][(σ2σ3)k][(σ3σ1)k]=0(2-5)在π平面上，有σ1σ3kX=p=p=常数(2-6)66这在30◦Cθ30◦C范围内是一条平行Y轴（或σ′）的直线，将其对称就2成为正六边形，如图2-1a)所示。在主应力空间内，这个屈服面按式(2-5)为三对相互平行的平面组成，并垂直于π平面的一个正六面柱体，如图2-1c)所示。在平面应力状态，并规定σ2=0情况下，则式(2-5)变为σ1σ3=k,σ1=k,σ3=k(2-7)在σ1σ3应力平面上，相当于六条直线，构成如图2-1b)所示的六边形。(2)Mises准则Tresca屈服准则不考虑中主应力的影响；另外当应力处于两个屈服面的棱线处时，在数学上处理时会遇到困难；在主应力方向不知时屈服准则又很复杂。因此Mises（1913）在研究了实验结果之后，提出了另一种屈服准则，即当J2=k1(2-8)–11– 哈尔滨工业大学工学硕士学位论文σ1σ1kO−kσ3Okσ2σ3−ka)Tresca准则在π平面上投影b)Tresca准则在主应力平面上投影a)TrescaCriterioninπplaneb)TrescaCriterioninprincipalstressplaneσ1σ3σ2c)Tresca准则在应力空间中c)TrescaCriterioninstressspace图2-1Tresca准则Fig.2-1TrescaCriterion时，材料就进入屈服。式中k———实验常数；J2———应力张量的第二不变量。早在1904年Huber就提出过同样的准则，故又称Huber-Mises准则，在材料力学中称为第四强度理论。因为–12– 哈尔滨工业大学工学硕士学位论文√2τoct=J2(2-9)3故式(2-8)可写成√√22τoct=J2=k1=k(2-10)331√或F=(σ1σ2)2+(σ2σ3)2+(σ3σ1)2k=0(2-11)3所以，Mises准则可以解释为：当材料中八面体剪应力达到极限时，材料开始屈服。由罗德应力参数对Mises屈服条件转化得：γ=τoct=k=常数(2-12)因此，在π平面上Mises准则则必定为一圆，如图2-2a)所示。在主应力空间是一个以λ线为中心线的正圆柱体，如图2-2c)。在σ2=0的平面应力情况下，Mises准则几何条件可表示成2292σ1σ1σ3+σ3=k(2-13)2其形状如图2-2b)所示，这是一个椭圆。如果在简单拉伸中确定实验常数，当屈服时σ1=σy,σ2=σ3=0，于是k=σy对Tresca准则p2(2-14)k=σy对Mises准则3此时，π平面上Tresca六边形角点到原点的距离为p1k22d=X=pp=σy(2-15)sin60◦C633而Mises圆半径为p2γ=τoctk=σy=d(2-16)3两种屈服条件这时重合，则Tresca六边形将内接于Mises圆，如图2-3a)所示。并在主应力空间中Tresca为内接于Mises正圆柱体的正六角形柱体，如图2-3b)。同理，如果采用纯剪切实验常数，这时两种屈服条件重合，则得到Tresca正六边形（或正六角柱体）将外切于Mises圆（或正圆柱体）。–13– 哈尔滨工业大学工学硕士学位论文σ1σ2σ3OOσ3σ1a)Mises准则在π平面上投影b)Mises准则在主应力平面上投影a)MisesCriterioninπplaneb)MisesCriterioninprincipalstressplaneσ1σ3σ2c)Mises准则在应力空间中c)TrescaCriterioninstressspace图2-2Mises准则Fig.2-2MisesCriterionTresca和Mises准则有一个共同点，即只有粘聚强度C，而没有摩擦角ϕ，故与静水压力无关。一般金属材料和饱和粘土具有这样的特点，两个屈服准则都为实验所证实，而且Mises接近实验成果。–14– 哈尔滨工业大学工学硕士学位论文σ1σ1⎫⎬Tresca⎭MisesOσ3σ2σ3σ2a)古典屈服准则在π平面投影b)古典屈服准则在应力空间中a)ClassicCriterioninπplaneb)ClassicCriterioninstressspace图2-3古典屈服准则Fig.2-3ClassicCriterion2.3.2广义Mises准则（Drucker-Prager准则）广义Mises准则是在Mises准则的基础上，考虑平均压应力而将Mises准则推广成为如下形式：σ2ÆDP1Mohr-CoulombDP3Mohr-CoulombDP5DP4r4r130◦30◦r23rσ3σ1ÆDP2图2-4各屈服准则在π平面上的曲线Fig.2-4Curvesofyieldcriterionsinπplane√F=αI1+J2=k(2-17)–15– 哈尔滨工业大学工学硕士学位论文I1=σ1+σ2+σ3(2-18)1222J2=[(σ1σ2)+(σ1σ3)+(σ2+σ3)](2-19)6式中I1,J2———分别为应力张量的第一不变量和应力张量的第二不变量；α,k———与岩土材料内摩擦角ϕ和粘聚力C有关的常数，不同的α,k在π平面上代表不同的圆（图2-4），各准则的α，k见表2-1。表2-1各准则α,k参数值Table2-1Parametervaluesofα,kofdifferentcriterions编号准则种类αk2sinϕ6ccosϕDP1外角点外接圆，=pp63(3sinϕ)3(3sinϕ)时，受压破坏2sinϕ6ccosϕDP2内角点外接圆，=pp63(3+sinϕ)3(3+sinϕ)时，受拉破坏ppDP3与Mohr-Coulomb截面等面积pp23sinϕpp63ccosϕ圆屈服准则23π(9sin2ϕ)23π(9sin2ϕ)psinϕp3ccosϕDP4平面应变关联法则下Mohr-3(3+sin2ϕ)3(3+sin2ϕ)Coulomb匹配准则sinϕDP5平面应变非关联法则下Mohr-ccosϕ3Coulomb匹配准则此式是1952年Drucker-Prager提出的，因此广义Mises准则也称为Drucker-Prager（D-P）准则。广义Mises准则在主应力空间的屈服面为一圆锥面，在π平面上也圆形，不存在尖顶产生的数值计算问题，目前国际上流行的大型有限元软件一般用它来近似替代Mohr-Coulomb准则。同理，在模型中中取不同θ值，即有不同的α,k值，由此得到大小不同的圆锥屈服面，通常将它统称为广义的Mises屈服条件，把α,k具有式(2-17)的形式称为Drucker-Prager屈服条件。当取θ=时受压破坏，可得DP1准则；当取θ=时受拉破坏，可66得DP2准则；当要求任何π平面上Mohr-Coulomb破坏面的下限，即π平面上的Drucker-Prager的截面圆内切于Mohr-Coulomb的六边形截面，于是得DP4准则，也称之为平面应变关联法则下Mohr-Coulomb匹配准则。由于计算α,k表达式还有很多，不同的公式计算结果也不一样。据Zienkiewicz等人研究，按不同的α,k所求极限荷载甚至相差达4–5倍之多。2.3.3Mohr-Coulomb准则早在1773年Coulomb就提出土的强度理论，其表达式为–16– 哈尔滨工业大学工学硕士学位论文τf=C+σntanϕ(2-20)式中C———土的粘聚强度；ϕ———土的内摩擦角；σn———受剪面上的法向应力，其符号取压应力为正。如果在物体内某一点的一个平面上有τn=τf(2-21)则该点土就破坏（即屈服）。这就是Coulomb屈服准则。该准则也是以剪应力的数值作为判断土的破坏标准，但是这个剪应力并不是最大剪应力，而是相应于土强度破坏时所达到的剪应力极限值τf（抗剪强度）。实质上土强度破坏是决定于滑动面上的剪应力与法向应力之比，或滑动面上作用应力的倾斜角达到其最大值，故Coulomb准则又称为最大倾角理论。Mohr（1900）根据强度破坏时极限应力状态的多个应力圈，求得与这些应力圆相切的包络曲线。虽然Leon（1934）及其他人提出该包络线为抛物线，双曲线型比较接近实际。但在静水压力不大的情况仍用Coulomb的直线表达该包络曲线以简化计算，Mohr-Coulomb准则在二维应力空间中可表示如图2-5。ττxyRMCϕσCyσxσnρ图2-5二维应力空间的Mohr-Coulomb屈服条件Fig.2-5Mohr-CoulombYieldConditionsintwo-dimensionalstressspace该直线方程可写成以下表达式：σ1σ3σ1+σ3=sinϕ+Ccosϕ(2-22)22σ1σ3σ1+σ3或F=sinϕCcosϕ=0(2-23)22这称为Mohr-Coulomb屈服准则。上式为表达平面应力状态的Mohr-Coulomb屈服函数。在应力空间中，所表达的破坏面为以λ线为中心并不规则的六角–17– 哈尔滨工业大学工学硕士学位论文锥，如图2-6c)所示。在π平面上沿坐标轴正负方向分别得出不同的截距，则得到的破坏轨迹为不规则六边形，如图2-6a)所示。这是因为法向受压时摩擦力使抗剪强度增大，而法向受拉时抗剪强度降低，因而三轴压缩与拉伸两者强度不同，屈服曲线沿坐标轴不对称，导致屈服曲线成不规则形状。σ1σ1σ3Oσ2σ3a)M-C准则在π平面投影b)M-C准则在主应力平面投影a)M-CCriterioninπplaneb)M-CCriterioninprincipalstressplaneσ1σ3σ2c)M-C准则在应力空间中c)M-CCriterioninstressspace图2-6Mohr-Couclomb准则Fig.2-6Mohr-CouclombCriterion按图2-5中表示，因为有：–18– 哈尔滨工业大学工学硕士学位论文1MC1MCσn=(σx+σy)Rsinϕ=(σ1+σ3)Rsinϕ(2-24)22√RMC=ccosϕ+ρsinϕ=(σσ)2+τ2=1(σσ)(2-25)xyxy122σxσyρ=(2-26)2MCτ=Rcosϕ(2-27)代入公式(2-20)有：σ1(1+sinϕ)σ3(1sinϕ)=2ccosϕ(2-28)若将主应力换成应力张量的第一不变量I1和应力偏张量的第二不变量J2及罗德角θ，则可得：()11√F=I1sinϕ+cosθpsinθJ2ccosϕ=0(2-29)33式中：θ。66Mohr-Coulomb准则在π平面上的图形为不等角六边形，存在尖顶和菱角，给数值计算带来困难。所以计算程序中采用Mohr-Coulomb准则时常要作一些近似处理，或采用与Mohr-Coulomb准则相应的广义Mises准则。本文采用的是Mohr-Coulomb准则。2.3.4Mohr-Coulomb等面积圆屈服准则在有限元数值计算中，直接采用Mohr-Coulomb屈服准则常引起不便。因而在假定土体或岩体为理想弹塑性体的情况下，与之近似的Drucker-Prager屈服准则被现有许多大型有限元程序采用。实际计算比较表明，按照该准则计算与Mohr-Coulomb屈服准则计算结果存在较大误差。为此，Zienkiewicz以及Pande等人提出了二次型屈服准则[20]去逼近Mohr-Coulomb屈服准则，但仍不够理想。徐干成、郑颖人曾提出等面积圆屈服准则[21,22]，且用工程实例验证了其良好的逼近效果[23–26]。等面积圆准则，徐干成、郑颖人（1990）提出，并成功的应用于边坡的稳定性分析中[21,27,28]。实际上是将Mohr-Coulomb模型转化成等效的Drucker-Prager准则，满足偏平面上的Mohr-Coulomb不等边六角形与Drucker-Prager圆面积相等的条件，其中α,k详见表2-1。–19– 哈尔滨工业大学工学硕士学位论文2.4本章小结本章首先阐述了几种边坡稳定性计算方法，着重归纳了岩土材料中土的各种屈服准则，并对各准则进行了论述比较。–20– 哈尔滨工业大学工学硕士学位论文第3章冻融土中水分冻融迁移模型3.1引言实践证明，含有粘性颗粒的所有湿土，在冻结前后，土体内水分会产生重分布现象。在土体上部，冻结后比冻结前的水分有不同程度增加，在土体下部则有不同程度地减少（无地下水补给时）或增加（有地下水补给时）。融化后，水分还会发生二次重分布。水分的迁移，需要分冻结、融化两个过程描述。3.2土冻结过程中的水分迁移3.2.1土冻结过程中的水分重分布自从认识土体冻结过程中水分重分布的规律来，人们花费了很大精力集中于水分迁移动力的研究。由于水分迁移过程受到许多复杂因素的影响，并随水分的相转换及各种类型水的变化而变化，使这个问题的研究至今仍没有一个统一的、十分完善的、并有定量依据的水分迁移理论。土体冻结过程中水分迁移理论仍是一个热烈争论的问题。世界上许多有关冻土方面的文献都论述到这个问题。许多学者都倾向于以下结论：土体冻结过程中以各种不同机理进行着水分迁移是粘性土冻胀形成的主要原因，薄膜水运动是其中最为主要的水分迁移机理。崔托维奇及奥尔洛夫把许多假说归纳为四种观点：(1)热量、流体动力观点这种观点认为，在热力学梯度作用下，毛细孔隙胶体介质中会发生热量与物质的传递迁移，致使水分在土孔隙中移动，称之为薄膜水迁移。(2)抽吸力、结晶力观点抽吸力观点认为，土矿物骨架的吸力作用造成水薄膜外层压力亏损而引起水分移动。水分迁移流量取决于土的吸力与外荷，以及冰晶体增长所产生的压力差值。而结晶力观点认为，在冰–水系统中，冰晶作用所产生的压力梯度会引起水膜移动。也有人认为，水–冰边界上的吸力与水向冰晶生长方向迁移是取决于土体的孔隙尺寸，也可归结到抽吸力观点之中。事实上，抽吸力是一种–21– 哈尔滨工业大学工学硕士学位论文综合力，它包括吸附力，毛细管作用力，结晶力等等。因此，抽吸力随土体埋藏条件、岩土和化学成分、孔隙率而变化。(3)物理化学观点这种观点认为矿物颗粒表面、同一相内离子交换、不同相分界表面离子吸附的薄膜水，或土颗粒表面自由能等势差作用均是水分迁移的基本动力。这种势差取决于土颗粒大小、土的吸水作用力大小以及交换盐基构成的情况。(4)构造形成观点这种观点认为土粒的凝聚作用、分散作用压缩沉陷等等形成水分迁移和薄膜水结构特征是水分迁移的动力。总之，关于水分迁移的学说很多，无论哪种假说都未能全面地解释土体冻结过程中水分迁移现象。3.2.2土冻结中水分迁移微观机理的一种解释由水分重分布现象可以看出，土冻结过程中，冻结土与非冻结土接触地带及已冻土区内的水分迁移特点各不相同。前者可能以毛细水及薄膜水迁移为主，后者以冻土中未冻结的薄膜水迁移为主。不论哪个区，水分迁移现象是随着亲水性物质的冻结而连续不断地进行着。根据上述解释水分迁移的理论，土体冻结过程中促使水分迁移的各种动力都是与冻土体内电分子力有关，都存在着土粒吸附作用、薄膜水结晶分凝作用，以及土粒（和冰）之间粒面的分子引力相互作用。按照表面现象的化学概念，含水的细分散土体性质的多样性都是由土的固体成分、表面积、孔隙溶液成分、重力及孔隙空间等因素综合作用所造成，其中主要的是土粒表面吸附层的成分与构造。同水分迁移动力有关的现象，如毛细现象、渗透现象、吸附现象、热力转移现象等等，都是与土的矿物成分有关，且以吸附层为媒介同孔隙水发生相互作用而表现出来的。用能量观点研究土与水的相互作用规律，一定量的质量必然是与一定量的任何形式能量之间相互联系着。张君常曾论述了土孔隙空间既是一个引力场又是一个电磁场，而且又是一个能量连续的势场[29]。在这个孔隙空间区域里，任何一个水分子都是在这势场的作用下保持其状态平衡。当水分子为质点运动时，引力所做的功只与质点的起点和终点位置有关，而与路径的形状无关。同样，正如胶体化学家所证实，土胶体粒子是带有电荷的，不论土粒是什么形状，都可以形成一个到几个等势面。在等势面上所有各点都是具有相同的–22– 哈尔滨工业大学工学硕士学位论文电势，这就构成土胶体双电层，即内层———密集部分；外层———松散部分。土粒表面吸附层上是多水分子层，每层水分子层，吸附力较强的一端，向着内部，暴露在外的部分具有较弱的吸引力。这就使得土粒表面附近的水分子紧密地、整齐地排列起来，越接近土颗粒表面就被吸附得越紧，越远则越弱。土粒吸附力是随距离的增加而骤减。超出吸附力范围的水分子就成为自由水。在土粒四周的水膜，称为“水化膜”，紧靠土粒表面的胶结水层，称为强结合水层，余外的为外吸附水层，也称薄膜水层。这种“水化膜”厚度，主要取决于土颗粒表面力场的强弱。土体冻结过程中，自由水结晶也构成相应的质点。随着土温下降，吸附层水膜结晶，厚度减薄，导致土粒表面电力势的变化。每当土的相态平衡被破坏时，土体中冻与未冻界面上的势差就成为未冻结土中水分向冻结锋面迁移的动力。所以，土体的相态平衡破坏及外部温度、压力等的因素发生变化时，水分迁移现象就会连续不断地进行着。土温降低是土体冻结的基本因素，也必然成为水分迁移动力的因素。土温降低引起水结晶，冰分凝，土粒的自由能量增长，使得冷源方向存在着各种分子力，成为土体内部液态水向冷锋面不断迁移的动力。在已冻土中，当冰面及土粒表面存在着未冻水时，在温度梯度的影响下，未冻水迁移也遵循正温土体薄膜水迁移的规律。但是，这种能克服重力的薄膜水向上迁移相当缓慢，使土粒间的公共水化膜是这种薄膜水迁移的基本动力的解释受到质疑。3.3土冻结的水分迁移模型3.3.1模型的物理过程描述作者导师没有遵循上述在冻结过程中土中水分的迁移理论中比较流行的薄膜水迁移理论，考虑非饱和土体中水分迁移量与饱和度密切相关，饱和度降低，水分迁移机制逐步由毛细水、薄膜水向气态水的扩散过渡，和刘兵一起直接采用土中水汽的扩散理论来建立水分迁移模型，应用改进的Fick第一定律来定量地研究冻结过程中水分的迁移。刘兵（2008）通过一系列试验建立的解释中，假定冻结前土中的水分是均匀分布的，土中的各状态水之间处于动态平衡，土表面从正温向设定负温变化需要一段时间，冻结锋面向下移动的比较慢，指出土层中的水分有足够时间发生迁移。迁移机理解释为土空隙中的气态水分接近冻结锋面就被吸附、冻–23– 哈尔滨工业大学工学硕士学位论文结，使得附近水汽的浓度降低，下面土层中水汽的浓度大，在浓度梯度的作用下，下面土层中的气态水向上扩散，到冻结锋面处的水汽又遇冷凝结成水、冻结成冰晶。水汽不断向上移动，使下部土层的气态水浓度降低，部分液态水向气态水转变，补充水汽浓度损失。从冻结锋面往下的某一深度处，浓度梯度引起的水汽迁移力小于土对水汽的阻碍作用时，气态水就不会发生迁移了，对应的深度称为临界面，或冻结锋面的影响深度。随冻结锋面的向下发展，作用的时间加长，影响深度逐渐增大[18]。3.3.2模型的数学表达受土样冻融试验的限制，只能建立一维模型，从最简单竖向的一维问题入手。3.3.2.1模型的假设(1)非饱和土中水分在冻结过程中的迁移是以水汽扩散的形式移动的，液态水通过汽化转化成水汽参加迁移，不考虑具体转化过程。(2)冻结前后土的干密度不变。(3)土中水汽移动的通道是土中的孔隙，孔隙是连通的。3.3.2.2Fick第一定律在模型中的应用Fick第一定律表述为扩散物质通过单位面积的流量与其浓度成比例dMdC=Dv(3-1)dtdz式中dM———单位时间通过单位面积的水汽质量；dtDv———水汽扩散系数；dC———在z方向上水汽的浓度梯度；dzM———在冻结过程中通过单位面积的水汽质量；C———单位土体中水汽质量。试验表明，水分迁移，即水汽扩散与土中温度梯度有明显的关系。为此，在Fick定律的公式中补充一项，式(3-1)变为：dMdTdC=Dv(3-2)dtdzdz式中T———土中一点温度；dT———试样中任意两个截面的温度梯度，取有水汽通过的平均值。dz–24– 哈尔滨工业大学工学硕士学位论文以冻结后的含水量分布曲线为目标，借助遗传算法参数反演了的数值[30]，两组结果列于表3-1。表3-1冻结过程反演的参数值Table3-1Parametervaluesinfreezingprocess反演参数aD(104)hb(103)v试样120.0346079929.927671130.00033570试样160.0112918529.990844450.000046453.3.2.3计算过程(1)上述冻结过程模型，如图3-1所示，坐标沿深度方向为正，用Z表示，以土层顶面为坐标原点。用H表示冻结锋面在坐标系Z中的位置，取z为至冻结锋面的深度。HChLzGdZ图3-1冻结模型示意图Fig.3-1Thesketchoffreezingmodel(2)冻结锋面附近的水汽被吸附、凝结成冰而导致水汽浓度减小，假定锋面处浓度为零，浓度梯度最大，截面越远离锋面，水汽浓度越高、梯度越小，至临界截面，浓度梯度减小至零，并不再变化。随着时间发展，冻结锋面的影响深度逐渐增大。浓度随深度的变化方程在zc坐标系中表示，图形如图3-1所示。设浓度随深度的函数为：1−3bt2−2btC=aez+ahez(3-3)2–25– 哈尔滨工业大学工学硕士学位论文式中h———冻结锋面的初始影响范围；a、b———浓度方程的两个参数。(3)计算试样中任意一截面G在冻结过程中的水汽通过的质量。截面G到土柱顶端的距离为Z，冻结过程中通过截面G的水汽总质量为：∫t(Z)btM(t)dt当zhe时M=M(Z)=∫0(3-4)t(Z)btM(t)dt当z>he时t(Z−hebt)在冻结过程中厚度d的土层中水分的变化量∆m为：∆m=M(Z+d)M(Z)(3-5)(4)冻结结束后，厚度d的土层中水的质量mw为初始水的质量加上冻结过程中该层内水分的变化量。mw=w0ρdAd+∆mA(3-6)则可得冻结后厚度d的土层的含水量mww0ρdAd+∆mA∆mM(Z+d)M(Z)wf===w0+=w0+(3-7)msρdAdρddρ0d式中wf———厚度d土层的含水量；mw———厚度d土层里水的质量；ms———厚度d土层里干水的质量；w0———土层的初始含水量；ρd———土层的干密度；A———土柱的横截面积。根据上述步骤求得的曲线二例，如图3-2所示，从中可见，和试验的结果大体上相符。3.4土融化过程中的水分迁移3.4.1土融化过程中的水分重分布及物理过程描述冻土自地表向深部融化，上部温度高，土中的冰晶迅速融化,液态水在重力作用下，向下渗流。随着融化面向下发展，水分逐渐向下迁移。冻结土层不透水，渗流到达融化-冻结界面会受阻。融化面向下发展的速度小于水分向下渗流–26– 哈尔滨工业大学工学硕士学位论文ω%ω%20.420.721.021.321.621.922.222.517.818.018.218.418.618.819.019.20045810试验含水量曲线试验含水量曲线12计算含水量曲线计算含水量曲线15hcm16hcm202024252830a)试样12的ωh曲线b)试样16的ωh曲线a)Curveofωhfor12thsampleb)Curveofωhfor16thsample图3-2ωh曲线Fig.3-2Curveofωh的速度时，下渗的水分就会在融冻界面之上积聚，形成水分富集层[31,32]，当某一层土的含水量达到持水量时，水分不再向下渗流[33–35]。3.4.2模型的数学表达如前所述，仅考虑竖向渗流，采用下述增加温度梯度影响项的非饱和土的达西定律来建立模型。3.4.2.1考虑温度梯度影响非饱和土的达西定律在融化过程中各土层的截面温度是随时间变化的，在应用非饱和土渗流的达西定律时，考虑温度对水分运移通量的影响。模型中用温度梯度乘以温度梯度引起的水分运动通量来反应温度的影响[36–39]，将一般形式的达西定律修正为：dθdTq=D(θ)K(θ)DT(3-8)dzdz式中q———单位时间单位面积上水分通过的质量；D(θ)———水分扩散系数；d———含水量梯度；dzK(θ)———土的导水系数；DT———温度梯度作用下的水分扩散率；dT———温度梯度；dzθ———分界面出的体积含水量。–27– 哈尔滨工业大学工学硕士学位论文3.4.2.2计算思路(1)据已有的经验公式可知土的水分扩散系数和导水系数都与含水量成指数关系[36,40]，土的水分扩散系数和导水系数可写成如下形式：beD(θ)=a(θ),K(θ)=c(θ)(3-9)式中a，b，c，e———模型参数。D(θ)———水分扩散系数；(2)融化过程的模型如图3-3所示，沿深度方向以土层顶面为坐标原点建立Z坐标轴，用以表示截面的位置，用H表示融冻界面在坐标系Z中的位置。把融化的土层等距离分层若干小层。HZ图3-3融化模型示意图Fig.3-3Thesketchofthawingmodel(3)通过预测的温度场，可以得到融化界面到达各测点的时间，然后拟合成曲线。通过曲线计算融化界面到达各个分层界面的时间，进而计算出融化界面在各个分层内的发展时间∆t。(4)融化过程中，当融化界面到某一深度时，该处截面开始有水分通过，直到试验停止。当分层足够小时，∆t就足够小，认为在∆t内的含水量不变，各分界面处的温度梯度不变，含水量梯度不变。分界面的含水量认为是分界面上层土的含水量，含水量梯度按分界面上下土层的含水量差与间距的比值计算，温度梯度按上下土层的温度差与间距的比值计算。(5)通过式(3-8)和式(3-9)可以计算出融冻界面从第i1分界面到达第i分界–28– 哈尔滨工业大学工学硕士学位论文面处时，第j个分界面处的在时间段内水分通过的量(j