双梯度驱动的非饱和冻土水热耦合模型_吕清硕.pdf

青藏高原冻土地区的地质灾害与冻土因冻融作用引起的水分重分布密切相关，其实质是复杂水热耦合作用的结果。已有的冻土水热耦合模型较多关注未冻水含量梯度驱动引起的水分迁移过程，而缺乏对温度梯度驱动效应的探讨。基于经典热传导方程和非饱和土体渗流理论，考虑未冻水含量梯度和温度梯度的共同作用，建立了双梯度联合驱动作用下的非饱和冻土水热耦合模型。在采用已有试验数据对模型有效性进行验证的基础上，分别对、和 种环境气温条件下土体的水热响应开展了数值模拟。结果表明：温度梯度在冻结过程中对于水分迁移的驱动作用不可忽略；冻结过程中由于冰水相变释放潜热造成冻结速率逐渐减慢；土中水分聚集的位置处于冻结初期形成的冻结锋面处，且外界气温越低，内外温差越大，则土体水分发生聚集的位置越深，水分迁移量相对也越大，冻结过程中土体水分随深度呈 型分布；、和 环境温度下含水率极值分别位于、和 深度处，含水率增量分别为 、和 。关键词：冻土；水热耦合模型；冻融效应；温度梯度；未冻水含量梯度中图分类号：；文献标识码：文章编号：（），（，）：，：；研究背景青藏高原冻土地区特有的融冻泥流和热融滑塌等地质灾害问题长期受到广泛关注，此类灾害的诱因与冻土对环境条件极其敏感的特性密切相关，其实质是在气温、降雨等气候环境作用下，土体内部的复杂水热耦合作用造成水热重分布而打破了原有的平衡。为了准确揭示高寒地区融冻泥流和热融滑塌等地质灾害发生的机理，获得复杂气候环境作用下土体内部温度场和水分场的演化规律，亟需建立有效的非饱和冻土水热耦合数学模型，从而为地质灾害评价和解决寒区工程问题提供理论依据。在冻土水热耦合过程方面，科研工作者们已开展了大量的研究。李述训等经试验研究发现，土中水分在冻结过程中发生迁移聚集而出现明显的聚冰现象是冻土区别于常规土体的显著特性，因此引起土中液态水迁移的驱动力成为冻土在冻融过程中水热耦合模型研究的重点。对于耦合模型的研究，最早可以追溯到 年 率先提出的水热耦合概念，他基于热传导方程和质量守恒方程建立了由土水势驱动水分迁移的水热耦合模型，但未给出土水势的具体影响因素。随后 年，等对于 模型中驱动水分迁移的土水势开展了进一步讨论，认为土水势是由于土中未冻水含量存在梯度而引起，并认为土水扩散系数是与未冻水含量有关的函数，由此建立了由未冻水含量梯度驱动水分迁移的冻土水热耦合模型。李洪升等、等、等认为水分迁移是重力势和冰压力共同作用的结果，并通过引入 方程建立了考虑重力势和冰压力作用的冻土水热耦合模型。白青波等、等采用固液比表示冻土中孔隙冰体积与未冻水体积之比，引入冰阻抗因子考虑冻土中孔隙冰对于未冻水迁移的阻滞作用，基于非饱和土体渗流特性进一步发展了冻土水热耦合模型。目前，对于冻土水热耦合模型的研究已取得了丰硕的成果，但土体冻融过程中水分迁移的驱动力十分复杂，仍未有统一定论，其中对未冻水含量梯度的驱动作用已形成共识，但温度梯度的影响在现有水热耦合模型中鲜有体现。土体冻结过程中，上下边界存在的温度梯度将诱发水分迁移。等 基于不可逆热力学深入分析了温度梯度引起的孔隙水运动机制。温度梯度导水作用强弱常用温度导水系数表示，取值范围为 （）。等通过试验发现只有当温度导水系数大于 （）时，温度梯度作用才会显著影响土体孔隙水流动。本文在已有冻土水热耦合模型的基础上，引入温度导水系数 描述温度梯度的导水作用，进一步扩展液态水迁移的驱动力，由此建立了考虑未冻水含量梯度和温度梯度共同驱动的非饱和冻土水热耦合模型，采用多物理场耦合数值模拟软件 进行求解、验证，并对不同外界环境温度条件下的冻土内部水热耦合过程开展了数值模拟，获得了温度及水分随时间的响应规律，揭示了冻结过程中的液态水迁移规律及温度传导演化过程。双梯度驱动的非饱和冻土水热耦合模型 基本假设（）土体为均质、连续且各向同性；（）土骨架无收缩与膨胀变形；（）冻融过程中水分以液态形式运移，忽略水分运移过程中冻融过程对土中未冻水和热传导的影响；（）冻土中未冻水含量与负温保持动态平衡。温度场控制方程在冻融过程中，冻土内部冰水相变作用对温度场分布的影响较为显著。在经典热传导方程的基础上，考虑冰水相变诱发的温度效应，可得冻土在冻融过程中的二维热传导方程：（）（）式中：为土体比热，（）；为土体密度，；为土体温度，；为时间，；为土体导热系数，（）；为冰水相变潜热，取值为 ；为冰的密度，；为含冰率。水分场控制方程冻土在冻融过程中水分迁移的驱动力除未冻水含量梯度外，温度梯度对液态水流动的影响亦不可第 期 吕清硕，等：双梯度驱动的非饱和冻土水热耦合模型忽略，。为更加完整和准确地揭示水分迁移过程和行为，基于非饱和土渗流理论，在 方程的基础上，引入 来考虑温度梯度引起的孔隙水流动，进而描述土体在冻结过程中未冻水含量梯度和温度梯度双梯度作用下的水分迁移现象，则水分场控制方程可以表示为：（）（）（）（）式中：为土中未冻水含水率；为土水扩散系数，；为温度导水系数，（）；（）为重力方向土体渗透率，；和分别为水和冰的密度，；为哈密顿算子。土水扩散系数定义为非饱和土渗透率（）与土壤比水容量（）的比值，然而在冻土中因孔隙冰的存在使得液态水的迁移受到阻碍，引入冰阻抗因子 表示冻土内孔隙冰的阻滞作用，于是非饱和冻土中土水扩散系数可以表示为：（）（）（）（）公式（）中非饱和土体渗流参数可由 模型（）确定：（）（）（）（）（）（）（）（）式中：为饱和含水率，；为残余含水率，；为饱和度；为饱和导水系数，；和 均为渗流参数。耦合变量温度场方程与水分场方程之间的耦合关系通过土体中的未冻水含量与温度之间的动态平衡关系建立。徐斅祖等基于大量试验数据，得到了冻土中未冻水含量与土体温度的动态平衡关系：（）（）（）式中：为初始含水率，；为土体冻结温度，；为土性参数，经试验获得，若缺乏试验数据，可按砂土 ，粉土 ，黏土 取值。通过将公式（）（）联立则构成了双梯度驱动的非饱和冻土二维水热耦合模型。该模型通过引入未冻水含水率与温度间的定量关系，实现了水分场和温度场的全耦合。求解公式（）得到冻结过程中土体内部温度分布后，结合公式（）得到土体内部未冻水含量，通过公式（）计算得到土体内部水分场分布变化。冰水转化过程中潜热的释放再次影响温度场分布，此时通过公式（）中相变潜热项计算实现温度场和水分场的全耦合过程。非饱和冻土水热耦合模型验证 软件具备强大的多物理场耦合求解功能，故本文选用该软件中的偏微分方程模块（）对所建模型开展数值求解并进行模型验证。采用文献开展的冻结试验对所建模型的有效性和适用性进行对比验证。冻结试验采用封闭系统，试验土样为内蒙古黏土，高度为 ，冻结方向自下向上，冻结时间为 ，冻结试验装置如图 所示，数值模拟所采用的计算参数见表、。需注意的是该试验未测定温度导水系数，因此本节参考已有研究成果，将其取值为 （）。图 冻结试验装置示意图表 试验条件初始含水率 顶部温度 底部温度 冻结温度 表 试验土性参数饱和渗透系数 （）土性参数 残余含水率 饱和含水率 非饱和土渗流参数 非饱和土渗流参数 土体冻结温度 水 资 源 与 水 工 程 学 报 年 图 给出了仅考虑未冻水含水率梯度与考虑双梯度驱动下的含水率随土层深度变化规律的数值模拟结果及冻结试验结果。因土样自下向上冻结，顶部为恒温 ，底部为恒温 ，致使土样在 后仍然无法完全冻结。图 表明，冻融交界面在约 深度处形成，最终的冻结区域深度约为 ，未冻结区域深度约为 。冻结过程中，液态水冻结相变转化成冰，未冻水含量减少，形成未冻水含量梯度，并且由于试样顶部与底部温差始终存在，温度梯度同样起到导水作用，造成未冻结区域内水分不断向冻结区域迁移，冻结区域内含水率明显增加。通过与冻结试验结果对比可以明显看出，考虑未冻水含水率梯度与温度梯度联合作用的耦合模型数值模拟结果与冻结试验结果更为接近，两者吻合度较高。因此，温度梯度在冻结过程中对于水分的迁移具有驱动作用，所建立的双梯度驱动非饱和冻土水热耦合模型能够准确揭示冻结过程中的水分迁移过程和冰水相变转化过程。复杂环境温度下冻土水热耦合过程数值模拟 模拟工况高寒冻土地区海拔高，气压低，空气稀薄，大气逆辐射弱造成保温能力较差，故高海拔地区气温的特点主要表现为零上温度低，正负温差大，降温速率快，因此高寒地区冻土经历着复杂环境气温作用下的冻结过程。为揭示环境温度引起的土体内部水热响应，选取、和 个特征环境温度，采用所建耦合模型开展数值模拟，各模拟工况均为自上而下冻结，模型参数取值见表。图 土样含水率随土层深度变化规律的数值模拟及冻结试验结果（冻结 ）模拟结果与分析在 、和 的环境温度条件下冻结、后，土体温度及含水率随土体深度变化的模拟计算结果如图 所示。表 数值模型参数取值土体密度 （）冰密度（）水密度（）冰水相变潜热（）土体热容 （）土体导热系数（）温度导水系数（）饱和渗透系数（）土性参数 残余含水率 饱和含水率 非饱和土参数 非饱和土参数 土体冻结温度 注：表中各参数取值分别参考文献 、。图 不同环境温度和冻结时间条件下土体温度随土体深度的变化曲线第 期 吕清硕，等：双梯度驱动的非饱和冻土水热耦合模型图 不同环境温度和冻结时间条件下土体含水率随土体深度的变化曲线 由图 中土体温度分布可知，深度范围的土体在 种不同外界环境温度条件下，经过 后温度均降至 以下。外界环境温度越低则土体冻结速度越快，当外界环境为 时，范围内土体温度降低至 所需时间从 缩短为约 。在外界负温条件下，土体内部温度逐渐降低，土中以液态形式存在的孔隙水亦随之出现相变转化成冰，并释放一定的潜热，进而阻延土体冻结进程。从图（）中可以发现，由于冻结过程中冰水相变潜热的释放，深度处土体降温速度逐渐放缓，第 个 内土体温速降低了；第 个 内土体温速仅降低了，速度放缓了；第 个 内土体温速降低幅度进一步减小为，可见冻结过程中水分场与温度场之间存在强耦合作用。在土体自上而下的冻结过程中，土中水分相变成冰，造成了仍以液态形式存在的未冻水在土体内部形成未冻水含量梯度，驱动未冻水由高含量处向低含量处迁移。同时土体在冻结过程中时刻存在温度的不均匀分布，由于土体内部存在温差而形成的温度梯度同样会使土中水分向冻结锋面处迁移，因此土中水分在温度梯度和未冻水含量梯度两种驱动力的共同作用下出现了重分布现象。对比图 中土体含水率分布可知，外界环境温度越低，土体水分明显聚集的位置越深：的环境温度下，土体水分在 深度处聚集；的环境温度下，在 深度处聚集；的环境温度下，在 深度处聚集。冰水相变的发生是动态转化过程，土体温度降低至 以下后土中冰开始形成，液态水含量降低，液态水的迁移开始发生。土体冻结后，渗透系数大幅降低，发生渗流的区域明显减小，结合温度场分布来看，水分明显迁移并聚集的位置均位于冻结初期土体温度开始进入 以下后形成的冻结锋面附近。在自上而下冻结的初始阶段，土体快速冻结，土中液态水迅速相变转化为冰态，未冻水含水率迅速减少，水分迁移现象并不明显。外界温度越低，土体冻结速度越快，深度越深，水分分布未发生扰动的范围相应越大。随着冰水相变潜热的释放，冻结速度开始放缓，未冻水开始在温度梯度和未冻水含量梯度的共同驱动作用下在冻结锋面处聚集，因此造成了液态水在一定深度区域聚集，而聚集处的上、下区域含水率变化均较小。在 的冻结过程中，水分出现明显的聚集现象，体积含水率随深度呈 型分布，外界环境温度越低，则作用在土体内部形成的温度梯度越大，从而引起水分迁移量的增加，如：在 的环境温度下，深度处含水率从 增加至（图（）；环境温度下，深度处含水率从 增加至 （图（）；环境温度下，深度处含水率从增加至 （图（）。的进一步冻结过程中，水分聚集处体积含水率有一定程度的减小，此时随着迁移而来的液态水继续冻结，水相变成冰后体积增大约，土体中水分的状态应是以冰的形式填充土中孔隙，而体积含水率定义为液态水和等效为液态水的冰共同组成，因此冰水相变后的体积变化是造成此时含水率降低的主要原因。结 论本文在现有冻土水热耦合模型研究的基础上，综合考虑未冻水含量梯度和温度梯度对于土体在冻结过程中水分迁移的共同作用，建立了双梯度驱动的非饱和土体水热耦合模型，在完成模型有效性验证的基础上，开展了土体水热响应的数值模拟，得到以下主要结论：（）建立的考虑未冻水含量梯度和温度梯度共同作用的水热耦合模型能够有效反映土体在冻结过程中发生的水分重分布现象，且温度梯度在冻结过程中对于水分迁移的驱动作用不可忽略。水 资 源 与 水 工 程 学 报 年（）土体冻结过程中冰水相变释放的潜热一定程度上会减缓土体的冻结速度，当外部环境为 时，深度处冻结速度放缓幅度达到。（）土体与环境温差越小，越有利于水分迁移，液态水分布发生扰动的深度越浅；相反土体与环境温差越大，则土体会迅速向下冻结，未扰动深度越深。当土体与环境温差增加至 时，土体水分分布发生扰动的位置加深 。（）所建立的未冻水含水率梯度和温度梯度联合作用的耦合模型能够真实反映冻结过程中土体的水热耦合过程，对研究青藏高原地区冻土退化及其诱发的次生灾害问题具有一定的参考意义。参考文献：邱国庆，程国栋 中国的多年冻土 过去与现在第四纪研究，（）：靳德武 青藏高原多年冻土区斜坡稳定性研究 西安：长安大学，牛富俊，马立峰，靳德武 多年冻土地区斜坡稳定性评价问题 工程勘察，（）：戴福初，邓建辉 青藏高原东南三江流域滑坡灾害发育特征 工程科学与技术，（）：李述训，程国栋，刘继民，等 兰州黄土在冻融过程中水热输运实验研究 冰川冻土，（）：，（）：，（）：李洪升，刘增利，梁承姬 冻土水热力耦合作用的数学模型及数值模拟 力学学报，（）：，：，：，白青波，李 旭，田亚护，等 冻土水热耦合方程及数值模拟研究 岩土工程学报，（）：，：，（）：毛雪松，李 宁，王秉纲，等 多年冻土路基水热力耦合理论模型及数值模拟 长安大学学报（自然科学版），（）：，：毛雪松，马 骉 基于水热耦合效应的冻土路基稳定性研究 北京：人民交通出版社，：张明礼，温 智，董建华，等 考虑降雨作用的气温升高对多年冻土活动层水热影响机制 岩土力学，（）：，：，（），（）：，：，（）：，（）：，（）：徐斅祖，王家澄，张立新 冻土物理学 北京：科学出版社，张志红，韩 林，田改垒 饱和土体热水力化全耦合一维溶质运移模型 东南大学学报（自然科学版），（）：，：，：，：第 期 吕清硕，等：双梯度驱动的非饱和冻土水热耦合模型