冻土损伤的本构模型和耦合问题的数值模拟

冻土损伤本构模型和耦合问题的数值模拟作者：朱志武①，宁建国②，马巍③作者单位：①西南交通大学力学与工程学院，成都610031；②北京理工大学爆炸科学与技术国家重点实验室，北京100081；③中国科学院寒区旱区环境与工程研究所冻土工程国家重点实验室，兰州730000收稿日期：2010-01-06；接受日期：2010-02-23摘要从材料的微观力学机理出发，建立了冻土损伤的弹性本构模型。对于不同温度和冰体积含量下的冻结砂土，由该模型计算的结果与实测的应力－应变曲线较吻合，建立的冻土损伤本构模型能够较好地描述实际冻土材料的力学性质。利用自行开发的有限元程序，加入所推出的本构关系和所建立的数学力学模型，通过对渠道冻结和冻土路基的水分场、温度场、应力场进行数值模拟计算，得到了比较准确、详尽的符合实际的温度场与应力场、位移场、应变场耦合的计算结果，与前人的计算及实测结果相吻合，规律一致。算例表明该程序能够计算冻土材料的相关物理量，且能很好地描述它们之间的关系。研究在前人的工作基础上，结合青藏铁路路基工程，通过对冻土水、热、力三场耦合机理这一固体力学学科领域的难点问题研究，为冻土工程的设计、施工和维护提供有效的科学依据、分析模型及必要的参考。关键词冻土，微观力学，路基，耦合，应力场1引言冻土是一种对温度敏感和易变的特殊低温地质体，多年冻土的面积约占全球陆地面积的23%，其中我国多年冻土地区约占国土面积的21。5%。随着社会的发展，冻土地区的道路工程、水利工程、建筑工程、矿山工程、能源工程等应运而生并迅猛发展，特别是随着青藏铁路的建成，关于冻土的研究显得越来越重要而急迫。冻胀融沉最本质的要素是土体冻融过程中水、热、力三场耦合问题，但由于问题的复杂性和研究手段的限制，目前研究仅以水、热两场耦合问题研究的最多，忽略应力场的问题，这样会导致冻融过程研究的不完整，进而影响到冻胀融沉预测的准确性。因此，水，热、力三场耦合问题是冻土研究领域的关键问题所在，也是国际前沿研究课题，研究结果可用于冻土构筑物的设计和稳定、安全评估中，促进冻土材料的大力发展，具有广阔的工程应用前景。目前人们对冻土中冰-土的相互作用已经有了比较深刻的认识和研究。Fremond[1]对正冻土的热力学性质进行了初步描述；Gary[2]对冻融循环进行了系统的分析；Konard[3]提出了一个描述正冻土的冰晶形成与水分迁移的模型；Dennis[4]对冻结过程进行了深入系统的试验研究；我国学安维东等人[5]先后对冻土的水分迁移与热质迁移，水热力耦合及其本构问题，进行了较深入的试与理论研究。然而以往的研究大多从热力学、混合物理论等角度出发建立起冻土的各种力学模型，多重于两场(温度场与水分场)的耦合作用，力场只是在分析冻结温度时，作为一项计算相变温度的指标被引进，专门从力学机理出发考虑冻土的三场耦合及其本构关系的研究尚未见报道。而在寒区的工程建设，，人们不可避免会遇到诸如冻胀融沉等冻土问题并开展相应的研究工作。Foriero等人[6]用有限元模拟了冻土边坡的蠕变；Su等人[7]对多年冻土区通风路基进行了数值分析研究；Dempsey[8]建立了非饱和土中水热耦合运移的数学预报模型；Rai等人[9]导出了流体、热流和土体变形完全耦合的控制微分方程；Masters[10]把位移、孔隙压力和湿度值当作初始未知量，研究了土体变形、流体温度和孔隙压力之间的耦合规律；Zhang等人[11]对寒区碎石路基在冻结过程中边界条件的影响进行了深入研究。对寒区工程中地基基础进行温度、水分、应力等的分析具有十分重要的现实意义，可以为实际施工和道路的正常运营提供有益的参考。因此，本文针对冻土工程中急待解决的土体冻融过程中水、热、力三场耦合的力学机理问题开展研究，结合青藏铁路路基工程，从材料细观力学出发，建立含损伤的冻土本构模型，在对冻土本构模型的研究中走出了一条新路。根据传热学，渗流理论和冻土力，，建立了冻土温度场、水分场、应力场耦合问题的数学力学模型，并对三场之间的耦合作用进行相应的数值模拟研究。2冻土的损伤本构模型2。1冻土的弹性模量从细观力学的角度出发，首先将冻土看作由土和冰组成的复合体单元。其次，把在整个冻土体中所占比例很大的土颗粒作为骨架，而把冰作为填充体看待，根据复合材料理论中的经典混合律思想将其耦合为冻土的本构关系，最后考虑加入损伤的影响。假设土颗粒为均匀连续体，土颗粒与冰之间完全黏结。令土和冰的弹性模量分别为Es和Ei，剪切模量为Gs和Gi，泊松比为νs和νi。下面各式中的下标s和i分别表示冻土中土和冰的各个分量。于是由细观力学的混合律理论可得到冻土复合材料的等效弹性模量K和等效剪切模量G为式中：和为土和冰的体积模量；和为土和冰的剪切模量。根据各向同性材料弹性常数之间的关系：K=E3(1−2ν)，G=E2(1+ν)。经推导后可得到由土、冰的弹性模量和泊松比表示的冻土等效弹性模量E和等效泊松比ν的表达式为[12]式中：和分别为土和冰的体积含量(也称为体积分数)。2。2冻土的损伤本构关系根据Lemaitre等效应力原理[13]，即应力σ作用在受损材料上引起的应变与等效应力σ作用在无损材料上引起的应变等价。受损冻土在一维状态下的应力-应变关系为上式中：σ为损伤力学中的等效应力，此等效应力相当于土力学中的有效应力；E为无损材料的弹性模量，即初始弹性模量，由(3)式确定；E为受损材料的弹性模量，即有效弹性模量；D为损伤量。因此，根据损伤力学理论[14]，得到冻土材料内部损伤型本构关系为σ=(1−D)Eε(6)2。3损伤演化规律在外载作用过程中，由于材料强度服从概率统计上的Weibull分布，因而可以认为材料的损伤量D朱志武等：基于损伤的冻土本构型及水、热、力三场耦合数值模拟研究也服从该分布。在此借助混凝土损伤的研究结果[15]，用双参数的Weibull分布表示冻土的损伤量：式中：ε为应变，n和a分别为形状参数和尺度参数，均为非负数。经推导，可得上式即为以应变作为损伤演化控制量的损伤演化方程，fε为应力峰值所对应的应变值，n为表征材料损伤演化特征的材料参数。将(8)式代入(6)式中，从而得到冻土的损伤本构关系为可以看出，冻土中任意点的应力σ与冻土的弹性模量、极限强度、应变峰值及该点的应变有关。冻土的损伤形状参数n是弹性模量E和割线模量Em的函数，一般可由实验或细观力学方法加以确定。从理论上讲，n值的大小与冻土材料的性态及外表破坏模式有一定的关系，即：n值越大，冻土越趋向于脆性破坏；n值越小，冻土越趋向于塑性破坏。2。4计算实例及分析根据已知土和冰的弹性模量，选取冰不同的体积含量，利用(3)和(4)式得到了冻结砂土在不同温度下的等效弹性模量如表1，冻结砂土的孔隙比e为0。4~0。5、天然含水量(w%)为15~18。利用表1中的弹性模量值，由(9)式计算得到的应力-应变曲线，与在相同温度下单轴压缩试验得到的应力-应变曲线[16]的比较，如图1和2所示。可以看出，温度是影响冻土变形的重要因素。随着土温的降低，土中未冻水含量减少，冻土中冰含量增加，固体颗粒间的胶结力增加，整个冻土体的变形式也由流变和塑性破坏逐渐转向脆性破坏。在相同的应力水平下，冻结砂土的变形随着土温的降低而减小，弹性变形在总变形中所占的比例却随着土温的降低而增大，所占的残余变形减小。另外，从图中还可以知道，弹性变形尽管随着冻土温度的降低而有很大增长，但其所占的比例仍然很小，特别是在冻土温度较高的情况下(一般高于−5℃范围内)，弹性变形只占总变形量10%~25%。在冻土温度较低时，弹性变形所占比例也不会有很大提高。为了更清楚地体现冰的体积百分比ic对冻土损伤强度的影响，我们将(9)式得到的在不同冰体积含量的应力-应变曲线进行比较，如图2所示。可以看出，在不同温度下，由于冻土中冰含量的增加，冻土体的损伤强度也相应提高。但在相同温度下，随着应变的增大，冻土体的压缩应力也逐渐达到最高点(屈服点)并已经开始下降，并且不同冰体积含量的冻土的应力屈服点也非常接近，说明在冻土材料破坏(屈服)时冰的体积含量对冻土体的屈服强度影响不大。3土体冻结水、热、力三场耦合数值模拟土体在冻结或融化过程中，温度、水分、应力三场的相互作用是一个极其复杂的热力学、物理化学和力学的综合问题。由于水-冰之间的相变以及水分迁移聚冰现象的存在，土体在冻融过程中温度场和水份场的变化将伴随着应力场的变化，同时温度场的变化也会引起水分迁移和应力状态的改变；反过来，应力场的变化又引起温度场和水分场的重分布。这三个场相互制约的关系存在于整个土体的冻融过程中，成为冻土材料本身所特有的水、热、力三场耦合问题。以自行开发的有限元程序为平台，加入本文推导出的本构模型，对实际工程中冻土的水分场、温度场及应力场之间的相互耦合问题进行数值模拟，与工程实测和前人所做的实验结果进行对比分析。3。1渠道冻结三场耦合数值分析3。1。1控制微分方程正冻结过程中渠道基土的水热耦合，可表述为二维非线性变系数抛物型偏微分方程：式中，C和λ为路基土的等效容积热容量和等效导热系数，其值分别为：未冻水含量uW与温度T的关系由室内试验确定[5]：式中，C+和C−为融土和冻土的容积热容量；λ+和λ−为融土和冻土的导热系数；Tp和Tb为剧烈相变区的上下温度边界；K为土体的导水系数，在剧烈相变区内随温度呈指数衰减，其初始条件和边界条件为：上述问题为温度、水分影响下以稳态温度场为初始温度条件的瞬态热传导问题，如果再加上外力边界条件和应力场的各种方程，特别是用我们提出的本构模型作为物理方程，则可以进行温度、水分影响下的应力计算，即进行初步的水、热、力三场的耦合计算。渠道静力平衡方程为：上述方程即组成了描述冻土渠道内含有相变的温度场、水分场、应力场的三场耦合方程组。3。1。2定解条件和各参数的取值在计算中取冰的体积含量ci为0。05，土的弹性模量Es为12。0MPa，冰的弹性模量取为4910。0MPa，土和冰的泊松比均取为0。3。整个平面上的初始温度为−0。3℃。上、下边界温度变化为：计算模型以不良工程地质条件下渠道基土为例进行数学描述[5]。渠道垂直截一断面为平面问题，认为土质均匀，水热蒸发耗热量及其它势场忽略不计。整个渠道模型宽为3。0m，高为2。0m，顶端渠堤宽为1。0m，渠底宽0。5m，渠底与下底边距离为1。0m。模型的位移边界条件为：下边界Y方向限制移动、X方向自由，左、右边界X方向限制移动、Y方向自由，上边界的X和Y方向均自由。应力边界条件为：左右边界剪应力均为0，右边界上X方向应力为0。模型上边界受相当于渠道衬切板重量和与衬切板垂直的冻结力的均布压力P=1400N，方向垂直于边界线，如图3所示。为了对比，图4给出了文献[5]计算时的外载荷施加位置图。各相关参数的取值为C+=3076。7，C−=2279。0(kJm−3℃−1)，λ+=5。53，λ−=6。58(kJm−1h−1℃−1)，图3本文计算时外载荷施加位置图4文献[5]的外载荷施加位置K1=1。8×10−0。7，K2=0。0036exp(0。551(T+0。3))，K3=3。6×10−3(mh−1)，L=334。7(kJkg−1)，=−0。30，pT=−0。75bT(℃)，0。0664exp(0。0551)，0。75，0。30580。596(0。30)，0。750。30，0。3058，0。30。uTTWTTT⋅⋅−⎧⎪=++−≤≤−⎨⎪⎩−3。1。3计算结果分析与讨论采用三角形网格划分，共分成12484个三角形网格单元和6463个节点，为了对整个渠道冻结过程中的温度、位移、应力、应变进行仔细的观察，我们取了较长的时间来进行计算，时间单位为小时，总共计算2000。0个时间，时间步长取为5。0。(1)温度场分析。在渠道截面上，不同冻结时间的0℃等温线分布计算结果如下图5所示，将其与图6的文献[5]计算的结果比较可以看出，在几乎相同的冻结时间，0℃等温线分布相似图5本文计算得到的0℃等温线随时间的变化1~8分别为50，100，200，300，500，800，1000，1500h