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

1、冻土损伤的本构模型与耦合问题的数值模拟冻土损伤本构模型和耦合问题的数值模拟作者：朱志武， 宁建国， 马巍作者单位： 西南交通大学力学与工程学院， 成都 610031； 北京理工大学爆炸科学与技术国家重点实验室， 北京 100081； 中国科学院寒区旱区环境与工程研究所冻土工程国家重点实验室， 兰州 730000收稿日期： 2010-01-06； 接受日期： 2010-02-23摘要 从材料的微观力学机理出发，建立了冻土损伤的弹性本构模型。对于不同温度和冰体积含量下的冻结砂土，由该模型计算的结果与实测的应力应变曲线较吻合，建立的冻土损伤本构模型能够较好地描述实际冻土材料的力学性质。 利用自行开发

2、的有限元程序，加入所推出的本构关系和所建立的数学力学模型，通过对渠道冻结和冻土路基的水分场、温度场、应力场进行数值模拟计算，得到了比较准确、详尽的符合实际的温度场与应力场、位移场、应变场耦合的计算结果，与前人的计算及实测结果相吻合， 规律一致。算例表明该程序能够计算冻土材料的相关物理量，且能很好地描述它们之间的关系。 研究在前人的工作基础上，结合青藏铁路路基工程， 通过对冻土水、热、力三场耦合机理这一固体力学学科领域的难点问题研究， 为冻土工程的设计、施工和维护提供有效的科学依据、分析模型及必要的参考。关键词 冻土， 微观力学， 路基， 耦合， 应力场1 引言冻土是一种对温度敏感和易变的特殊低

3、温地质体， 多年冻土的面积约占全球陆地面积的23%， 其中我国多年冻土地区约占国土面积的21。5%。 随着社会的发展，冻土地区的道路工程、水利工程、建筑工程、矿山工程、能源工程等应运而生并迅猛发展，特别是随着青藏铁路的建成， 关于冻土的研究显得越来越重要而急迫。冻胀融沉最本质的要素是土体冻融过程中水、热、力三场耦合问题， 但由于问题的复杂性和研究手段的限制， 目前研究仅以水、热两场耦合问题研究的最多， 忽略应力场的问题，这样会导致冻融过程研究的不完整， 进而影响到冻胀融沉预测的准确性。 因此，水，热、力三场耦合问题是冻土研究领域的关键问题所在，也是国际前沿研究课题，研究结果可用于冻土构筑物的设

4、计和稳定、安全评估中，促进冻土材料的大力发展，具有广阔的工程应用前景。目前人们对冻土中冰-土的相互作用已经有了比较深刻的认识和研究。 Fremond1对正冻土的热力学性质进行了初步描述； Gary2对冻融循环进行了系统的分析； Konard3提出了一个描述正冻土的冰晶形成与水分迁移的模型； Dennis4对冻结过程进行了深入系统的试验研究；我国学安维东等人5先后对冻土的水分迁移与热质迁移， 水热力耦合及其本构问题，进行了较深入的试与理论研究。 然而以往的研究大多从热力学、混合物理论等角度出发建立起冻土的各种力学模型， 多重于两场(温度场与水分场)的耦合作用， 力场只是在分析冻结温度时，作为一项

5、计算相变温度的指标被引进， 专门从力学机理出发考虑冻土的三场耦合及其本构关系的研究尚未见报道。而在寒区的工程建设， 人们不可避免会遇到诸如冻胀融沉等冻土问题并开展相应的研究工作。Foriero 等人6用有限元模拟了冻土边坡的蠕变；Su 等人7对多年冻土区通风路基进行了数值分析研究；Dempsey8建立了非饱和土中水热耦合运移的数学预报模型； Rai 等人9导出了流体、热流和土体变形完全耦合的控制微分方程；Masters10把位移、孔隙压力和湿度值当作初始未知量， 研究了土体变形、流体温度和孔隙压力之间的耦合规律；Zhang 等人11对寒区碎石路基在冻结过程中边界条件的影响进行了深入研究。对寒区

6、工程中地基基础进行温度、水分、应力等的分析具有十分重要的现实意义， 可以为实际施工和道路的正常运营提供有益的参考。因此， 本文针对冻土工程中急待解决的土体冻融过程中水、热、力三场耦合的力学机理问题开展研究， 结合青藏铁路路基工程，从材料细观力学出发，建立含损伤的冻土本构模型， 在对冻土本构模型的研究中走出了一条新路。根据传热学，渗流理论和冻土力， 建立了冻土温度场、水分场、应力场耦合问题的数学力学模型，并对三场之间的耦合作用进行相应的数值模拟研究。2 冻土的损伤本构模型2。1 冻土的弹性模量从细观力学的角度出发， 首先将冻土看作由土和冰组成的复合体单元。 其次，把在整个冻土体中所占比例很大的土

7、颗粒作为骨架，而把冰作为填充体看待， 根据复合材料理论中的经典混合律思想将其耦合为冻土的本构关系， 最后考虑加入损伤的影响。假设土颗粒为均匀连续体，土颗粒与冰之间完全黏结。 令土和冰的弹性模量分别为Es 和Ei， 剪切模量为Gs 和Gi， 泊松比为s 和i。 下面各式中的下标s 和i分别表示冻土中土和冰的各个分量。于是由细观力学的混合律理论可得到冻土复合材料的等效弹性模量K 和等效剪切模量G 为式中：和为土和冰的体积模量；和为土和冰的剪切模量。根据各向同性材料弹性常数之间的关系： K =E 3(1 2 )， G = E 2(1+ )。 经推导后可得到由土、冰的弹性模量和泊松比表示的冻土等效弹性

8、模量E和等效泊松比 的表达式为12 式中： 和分别为土和冰的体积含量(也称为体积分数)。2。2 冻土的损伤本构关系根据Lemaitre 等效应力原理13，即应力 作用在受损材料上引起的应变与等效应力 作用在无损材料上引起的应变等价。受损冻土在一维状态下的应力-应变关系为上式中： 为损伤力学中的等效应力，此等效应力相当于土力学中的有效应力； E 为无损材料的弹性模量，即初始弹性模量，由(3)式确定； E为受损材料的弹性模量， 即有效弹性模量； D 为损伤量。因此， 根据损伤力学理论14， 得到冻土材料内部损伤型本构关系为 = (1 D)E(6)2。3 损伤演化规律在外载作用过程中， 由于材料强度

9、服从概率统计上的Weibull 分布， 因而可以认为材料的损伤量D朱志武等： 基于损伤的冻土本构型及水、热、力三场耦合数值模拟研究也服从该分布。 在此借助混凝土损伤的研究结果15，用双参数的Weibull 分布表示冻土的损伤量：式中： 为应变， n 和a 分别为形状参数和尺度参数，均为非负数。 经推导，可得上式即为以应变作为损伤演化控制量的损伤演化方程， f 为应力峰值所对应的应变值， n 为表征材料损伤演化特征的材料参数。 将(8)式代入(6)式中， 从而得到冻土的损伤本构关系为可以看出， 冻土中任意点的应力 与冻土的弹性模量、极限强度、应变峰值及该点的应变有关。 冻土的损伤形状参数n 是弹

10、性模量E 和割线模量Em 的函数， 一般可由实验或细观力学方法加以确定。 从理论上讲， n 值的大小与冻土材料的性态及外表破坏模式有一定的关系，即： n 值越大，冻土越趋向于脆性破坏； n 值越小，冻土越趋向于塑性破坏。2。4 计算实例及分析根据已知土和冰的弹性模量， 选取冰不同的体积含量， 利用(3)和(4)式得到了冻结砂土在不同温度下的等效弹性模量如表1， 冻结砂土的孔隙比e 为0。40。5、天然含水量(w%)为1518。利用表1 中的弹性模量值， 由(9)式计算得到的应力-应变曲线， 与在相同温度下单轴压缩试验得到的应力-应变曲线16的比较， 如图1 和2 所示。可以看出， 温度是影响冻

11、土变形的重要因素。 随着土温的降低， 土中未冻水含量减少， 冻土中冰含量增加， 固体颗粒间的胶结力增加， 整个冻土体的变形式也由流变和塑性破坏逐渐转向脆性破坏。 在相同的应力水平下， 冻结砂土的变形随着土温的降低而减小， 弹性变形在总变形中所占的比例却随着土温的降低而增大， 所占的残余变形减小。 另外， 从图中还可以知道， 弹性变形尽管随着冻土温度的降低而有很大增长， 但其所占的比例仍然很小， 特别是在冻土温度较高的情况下(一般高于5范围内)， 弹性变形只占总变形量10%25%。 在冻土温度较低时，弹性变形所占比例也不会有很大提高。为了更清楚地体现冰的体积百分比i c 对冻土损伤强度的影响，

12、我们将(9)式得到的在不同冰体积含量的应力-应变曲线进行比较， 如图2 所示。 可以看出，在不同温度下， 由于冻土中冰含量的增加， 冻土体的损伤强度也相应提高。 但在相同温度下， 随着应变的增大， 冻土体的压缩应力也逐渐达到最高点(屈服点)并已经开始下降，并且不同冰体积含量的冻土的应力屈服点也非常接近， 说明在冻土材料破坏(屈服)时冰的体积含量对冻土体的屈服强度影响不大。3 土体冻结水、热、力三场耦合数值模拟土体在冻结或融化过程中， 温度、水分、应力三场的相互作用是一个极其复杂的热力学、物理化学和力学的综合问题。 由于水-冰之间的相变以及水分迁移聚冰现象的存在， 土体在冻融过程中温度场和水份场

13、的变化将伴随着应力场的变化， 同时温度场的变化也会引起水分迁移和应力状态的改变； 反过来，应力场的变化又引起温度场和水分场的重分布。 这三个场相互制约的关系存在于整个土体的冻融过程中， 成为冻土材料本身所特有的水、热、力三场耦合问题。以自行开发的有限元程序为平台， 加入本文推导出的本构模型， 对实际工程中冻土的水分场、温度场及应力场之间的相互耦合问题进行数值模拟， 与工程实测和前人所做的实验结果进行对比分析。3。1 渠道冻结三场耦合数值分析3。1。1 控制微分方程正冻结过程中渠道基土的水热耦合， 可表述为二维非线性变系数抛物型偏微分方程：式中， C 和 为路基土的等效容积热容量和等效导热系数，

14、 其值分别为：未冻水含量u W 与温度T 的关系由室内试验确定5：式中， C +和C 为融土和冻土的容积热容量； +和为融土和冻土的导热系数； Tp和Tb为剧烈相变区的上下温度边界； K 为土体的导水系数， 在剧烈相变区内随温度呈指数衰减， 其初始条件和边界条件为：上述问题为温度、水分影响下以稳态温度场为初始温度条件的瞬态热传导问题， 如果再加上外力边界条件和应力场的各种方程， 特别是用我们提出的本构模型作为物理方程， 则可以进行温度、水分影响下的应力计算， 即进行初步的水、热、力三场的耦合计算。渠道静力平衡方程为：上述方程即组成了描述冻土渠道内含有相变的温度场、水分场、应力场的三场耦合方程组

15、。3。1。2 定解条件和各参数的取值在计算中取冰的体积含量ci为0。05， 土的弹性模量Es为12。0 MPa， 冰的弹性模量取为4910。0 MPa， 土和冰的泊松比均取为0。3。整个平面上的初始温度为0。3。 上、下边界温度变化为：计算模型以不良工程地质条件下渠道基土为例进行数学描述5。 渠道垂直截一断面为平面问题， 认为土质均匀， 水热蒸发耗热量及其它势场忽略不计。整个渠道模型宽为3。0 m， 高为2。0 m， 顶端渠堤宽为1。0 m， 渠底宽0。5 m， 渠底与下底边距离为1。0 m。 模型的位移边界条件为： 下边界Y 方向限制移动、X 方向自由， 左、右边界X 方向限制移动、Y 方向

16、自由， 上边界的X 和Y 方向均自由。应力边界条件为： 左右边界剪应力均为0， 右边界上X 方向应力为0。 模型上边界受相当于渠道衬切板重量和与衬切板垂直的冻结力的均布压力P=1400N， 方向垂直于边界线， 如图3 所示。 为了对比， 图4给出了文献5计算时的外载荷施加位置图。各相关参数的取值为C+=3076。7， C=2279。0 ( kJ m3 1)， + = 5。53， = 6。58 ( kJ m1 h1 1)，图3 本文计算时外载荷施加位置图4 文献5的外载荷施加位置K1=1。8100。7， K2=0。0036exp(0。551(T+0。3)，K3=3。6103 (m h1)， L=

17、334。7 (kJ kg1)，= 0。30， p T = 0。75 b T ()，0。0664 exp(0。0551 )， 0。75，0。3058 0。596( 0。30)， 0。75 0。30，0。3058， 0。30。uT TW T TT 3。1。3 计算结果分析与讨论采用三角形网格划分， 共分成12484 个三角形网格单元和6463 个节点， 为了对整个渠道冻结过程中的温度、位移、应力、应变进行仔细的观察， 我们取了较长的时间来进行计算， 时间单位为小时， 总共计算2000。0 个时间， 时间步长取为5。0。(1) 温度场分析。 在渠道截面上， 不同冻结时间的0等温线分布计算结果如下图5

18、 所示， 将其与图6 的文献5计算的结果比较可以看出， 在几乎相同的冻结时间， 0等温线分布相似图5 本文计算得到的0等温线随时间的变化18 分别为50， 100， 200， 300， 500， 800， 1000， 1500 h 时的数值图6 文献5 0等温线随时间的变化17 分别为80， 320， 560， 800， 1040， 1280 和1600 h 时的数值冻结2000 h 后， 上边界的负温达到最大值10。12(如图7 所示)， 这与文献5计算得到的冻结1536 h 后上边界负温最大值10。7基本一致。 从图7中还可以看出， 从渠坡到渠底的温度变化快， 也即温度梯度大， 而在离渠坡

19、较远的渠堤下， 温度分布呈一组几乎平行的直线， 基本上不受渠坡边界温度的影响。 从温度场的分布和变化来看， 二维边界、截面几何形状、底水补给以及地中热流对温度分布均有显著影响。 从以上比较来看， 温度场的计算结果是令人满意的。(2) 位移场分析。 大量的工程实践表明， 在低冻结速率和有充足水分补给的条件下， 产生较强的冻胀作用。 在渠道底部和渠坡下部， 温度较高的地下水起到降低冻结速率和补给冻结锋面水分的双重作用，在渠底和渠坡下部发生较强烈的冻胀。 图8 是2000 h时的冻胀位移( x2 + y2 )分布， 可以看出， 最大冻胀位移发生在渠堤顶端和渠坡上部， 渠底和渠坡下部位移量也比较大，

20、与同时段的观测资料(图9)对比，可见渠道截面上的位移分布规律是基本一致的。图7 冻结2000 h 时的等温线分布图8 冻结2000 h 时的等位移线分布图9 宁夏大新渠试验段(南北走向)实测表面位移分布(3) 应力场分析。 图1015 是渠道冻结2000 h 后，各方向的应力和应变等值线分布图。 在渠底和渠坡下部的冻胀力都比较大， 在渠坡与渠底的交界点处出现了应力集中， 最大冻胀力也在此出现： 最大x 为3。943105 N/m2， 最大y 为0。85105 N/m2， 最大xy 为1。376105 N/m2。从图中可以看出， 冻胀应力沿渠道基土深度方向的分布状态。 在渠道基土上部渠堤处， 由

21、于外载荷图10 冻结2000 h 的等x 线分布图11 冻结2000 h 的等x 线分布图12 冻结2000 h 的等y 线分布的约束作用， 基土中的冻胀应力呈衰减趋势。 随着深度的增加， 外载引起的变形逐渐减弱而由冻胀引起的变形逐渐增强， 在超过某一深度后基土中的应力又会随着深度的增加而增大， 而使基土应力随深图13 冻结2000 h 的等y 线分布图14 冻结2000 h 的等xy 线分布图15 冻结2000 h 的等xy 线分布度的增加呈现衰减-增大-衰减的规律， 这与文献5中计算得到的结论相同。从以上各图还可以看出， 应力和应变存在着一一对应关系， 等应力线和等应变线的变化趋势相同说明

22、我们所用的弹性计算模型是正确的。为将计算结果看得更清楚一些， 我们选定了平面内的九个点作为研究对象， 分别来看各点的温度、位移、应力、应变随时间的推移而产生的变化。图16 所选点的位置示意图图17 各点的温度随时间变化曲线图18 各点y 方向的位移随时间变化曲线从图1620 可以看出， 在离渠道中心线(右边线)越远处， 温度变化得越剧烈， 在靠近渠坡处各点的位移变化比较明显。 另外， 由于左右两边界上位移条件的限制， 各个点水平方向的位移不是很大， 而在垂直图19 各点y 方向的应力随时间变化曲线图20 各点y 方向的应变随时间变化曲线方向的位移变化则比较大。 靠近渠坡与渠底的交界处， 应力变

23、化比较明显， 量值也比距渠堤较近处点的应力值大， 表明在渠坡和渠底的拐角处出现了应力集中， 靠近这个区域点的应力变化比较剧烈。 而同一点处在同方向的应力和应变值， 在变化趋势上是相同的， 它们之间存在一一对应关系。3。2 路基冻结三场耦合数值算例3。2。1 控制微分方程控制微分方程同上面算例一样， 连同方程(22)(25)即组成了描述冻土路基内含有相变的温度场、水分场、应力场的三场耦合方程组。3。2。2 定解条件取路基横断面以二维问题来处理， 并取路基的一半作为计算区域。 设计算区域由S1， S2， 􀀢， S6 六层均匀介质组成(图21)， 其中水泥路面的半宽为4。5m(包括

24、路肩1。0 m)， 路堤高3。0 m。(1) 边界条件。 根据青藏高原几十年的观测资料和附面层原理可得到固定边界上的边界条件如下17图21 路基温度场计算模型S1S6 为路基计算区域； B1B6 为路基区域边界( ) = + + 0 0 01 2 3( ) sin 2 3 ，8760 5， ，f t T R t A tx y B B B(27)= s T G ，n( ) 5 x， y B ， (28)=T 0，n( ) 4 6 x， y B B ， (29)式中T0 为下附面层底的初始年平均地温(地表以下0。5 m 处的实测温度值)， 其取值与文献17的取值一致(见表2)， R0 为气候变暖引

25、起的下附面层底地温增温率， A0 为下附面层底地温振幅， n 为边界界面的外法线方向， Gs 为计算区域下边界地温梯度， Gs=0。03/m。路基模型的位移边界条为： 下边界B5 限制y方向移动， x 方向自由； 左边界B4 和右边界B6 均限制x 方向移动， y 方向自由； 上边界B1， B2 和B3 的x方向和y 方向均自由。应力边界条件为： 左右边界B4 和B6 上的剪应力均为0， 右边界B5 上x 方向应力为0， 上边界B1、B2 和B3 上无载荷作用， 因而在温度变化过程中产生的应力全部为热应力， 土的热膨胀系数取为1。0105。表2 典型表面附面层参数值类型 天然地表 砂土边坡 水

26、泥路面附面层总温度增量() 2。5 3。7 4。7(2) 初始条件。 根据边界条件(27)(29)和计算区域土体的热参数及特性参数求解Laplace 方程：e T e T 0，x x y y + = (30)得到季节融化开始时计算区域温度场的估计值。3。2。3 参数的选取根据工程实际， 本算例的计算区域由六种不同的土质组成(见图21 所示)， 其热参数及土壤特性参数的取值如下表3 所示。 表中C 为路基土的容积热容表3 计算用土体热参数及特性参数土类参数 S1砂砾与碎石土 S2草炭亚黏土 S3含泥炭亚黏土S4 含碎石黏土S5 含碎石亚黏土 S6 弱风化基岩u(J m1 1) 6875。1 33

27、76。7 5110。8 2052。6 3124。2 7164。8f(J m1 1) 9396。0 4684。1 6249。1 3162。1 4376。2 9689。2Cu (J m3 1) 2。183106 2。258106 2。342106 2。166106 2。446106 2。784106Cf (J m3 1) 1。694106 1。781106 1。932106 1。539106 1。932106 2。192106d(kg m3) 1800 1200 1400 700 1300 2250W(wt%) 10 24 22 56 24 5。6Wu(wt%) 3 7 7 20 7 0。6D 9

28、。35106 1。65106 4。66106 1。21106 3。73106 3。44106量， 为路基土的导热系数， 下标f 和u 分别表示路基土体的冻、融两相状态， d 为土体的干容重(干密度)，W 和Wu 分别为路基土的重量含水量和重量未冻水含量， D 为路基土的水分扩散系数， 一般认为土的x 方向和y 方向的水分扩散系数相等， 即Dx=Dy。未冻水含量与温度的关系根据下式确定5：( )0。0644e0。0551 ， 0。75，0。3058 0。596 0。30 ， 0。75 0。30，0。3058， 0。30。TwTT TT (31)3。2。4 计算结果与分析计算采用三角形网格单元对整

29、个路基区域进行网格划分， 共分为11800 个单元， 6061 个节点。(1) 温度场分析。 根据边界条件(27)(29)和计算区域土体的热参数及特性参数求解拉普拉斯方程(30)， 得到季节融化开始时计算区域温度场的估计值，如图22 所示。在上述拉普拉斯方程解出的初始温度分布基础上， 利用温度条件表达式(27)中R0=0。0 时作为上边界的温度条件逐时段求解方程(10)， 直到年变化层以下温度场基本保持稳定， 而且年变化层以上相同位置的温度值在同一时间逐年相同时为止， 此时， 各个节点的温度值作为三场耦合计算时的初始条件， 温度分布如图23 所示为了对整个路堤冻结过程中的温度、应力、应变和位移

30、进行详尽的观察， 我们取了较长的时间来进行计算， 时间单位为小时(h)， 总共计算2000。0 个时间， 时间步长取为1。0。 图2427 是冻结不同时间后路基的温度等值线变化图。可以看到， 由于路基上边界有随时间变化的温度边界条件， 下边界有热流交换， 所以在冻结初始段路基内部的温度变化比较剧烈， 最低温度在冻结50 h后就由初始时刻的9。229下降到9。734， 在500h 达到了10。24； 最高温度也由初始时的0下降到了500 h 的3。988。 由于上边界是周期变化的温度条件， 所以经过一段时间后路基平面内的温度又有所回升， 最低温度由500 h 的10。24上升到了1000 h 的

31、9。821和2000 h 时的9。145。 此外， 从整个温度场的变化可以看出， 温度变化比较剧烈的区域都是在路基上部距离天然地表约5。0 m6。0 m 的深度范围内， 在此深度以下的大部分区域里温度变化都不明显。 在地表下约1。0 m3。0 m 的深度范围内存在一个温度变化更为频繁、剧烈的区域， 这就是通常所说的剧烈相变区域， 在此区域内水、冰之间由于冻融而产生相变转化。 这些现象与实际工程应用中的观察结果一致， 一般受气温变化影响的地表深度也是在5。0 m10。0 m 以内， 比这个深度更深的区域基本不受地表气温变化的影响。(2) 位移场分析。 由于本节对路基进行耦合计算时模型不受外载荷的作用