寒区土体一维水-热-力耦合模型与数值分析
时间：2023-06-23 14:05:04 来源：雅意学习网 本文已影响 人 
魏道凯，荆 皓，陈 琦，寇海磊
（1.山东交通职业学院公路与建筑系，山东潍坊 261206；
2.中国海洋大学工程学院，山东青岛 266100）
冻土在我国国土面积中占据着极其广阔的一部分，其中多年冻土约占1/4。寒区冻土的冻胀融沉特性使结构物产生不均匀沉降、结构断裂等现象，工程灾害频频发生，对上覆结构物的稳定带来了极大的挑战。冻土冻融作用的机理是当前冻土研究中最受关注的问题之一［1-5］。
在外界温度环境作用下，土体内液态水会因为土体内部温度低于土体的冻结温度引起物相改变，水变成冰，冰的体积较水的体积大，颗粒间的孔隙减小，进而实现土的冻胀。同时，冰的胶结作用使土颗粒间应力增大。故冻土冻融为温度场、水分场和应力场三场耦合的结果。对于三场耦合的研究一直是国内外研究学者关注的重点。Harlan［6］首次将达西定律应用于冻土水分场分析中，建立水热耦合方程。Mu等［7］首次提出冻土中水热两场的变化会引起土体内部应力场的变化。近些年关于冻土冻融机理的研究都是在此基础上进行的。Zhou等［8］基于人工冻结试验将土体冻胀引入水热力耦合模型中，探究了水热作用下土体变形特性。基于水热力三场耦合模型，Qi等［9］对冻融环境下膨胀土边坡的稳定性进行了分析。周志［10］通过水热耦合模型对冻融环境下粉质黏土内温度场以及水分场进行了研究。冉洪伍等［11］对水动力学模型、刚性冰模型以及热力学模型3种耦合模型的应用以及准确性展开了研究。但是已有研究多针对于饱和土体，对于非饱和冻土在水热两场作用下应力场以及冻胀变形的研究较少。
文中基于质量守恒定律以及能量守恒定律推导水分场控制方程及温度场控制方程，由水分场和温度场作用产生土体内部应力变化，实现土体的冻胀，建立与应力场的耦合。以此建立一维土柱数值模型，探究三场耦合的可行性。对冻融条件下土柱温度变化、水分迁移以及应力分布进行了研究。
1.1 基本假设
在保证模型准确性的基础上忽略部分情况，充分考虑冻土中水热力因素，简化模型建立过程，做出以下基本假设：
（1）在土体中，水分迁移主要以液态形式进行迁移，在模型中忽略气态水的迁移；
（2）土体中水分的迁移服从广义达西定律；
（3）土体各向同性，均匀连续；
（4）计算过程无温度损失，冻土内含水量以及温度处于平衡状态。
1.2 水分场控制方程
根据质量守恒定律建立水分场控制方程。二维土体微单元水分迁移模型如图1所示。
图1 土体微单元水分迁移示意图Fig.1 Schematic diagram of water migration in soil micro-units
式中：Δmx为单位时间土单元中x方向水分质量变化；
vx为x方向水流通量；
Δmz为单位时间土单元中z方向水分质量变化；
vz为z方向水流通量；
ρl为水的密度。
则单位时间内土体微单元的水分增量为：
单位时间内土体微单元水分质量变化有可从未冻水体积含量θl以及冰体积含量θi的变化进行分析。
联立方程可得：
式中ρi为冰的密度。
水分迁移服从广义达西定律，则水流通量v：
式中：k（θl）为土体的渗透系数；
ψm代表基质势；
z代表重力势。
最终得到水分场控制方程：
1.3 温度场控制方程
基于能量守恒定律建立土柱内温度场控制方程。
图2表示土体微单元内的热量迁移，与水分场类似，在单位时间内，土体微单元的能量变化ΔQ为：
热通量q又为导热系数λ与温度T的函数：
当温度降至土体冻结温度时，土颗粒间液态水会发生物相改变以释放能量［12-13］，因此土体微单元内能量变化又是由土体热容C以及相变潜热L决定的。
联立方程建立温度场控制方程：
水分场、温度场控制方程中共有θl、θi以及T这3个未知函数因此需要还需一个联系方程实现水热两场的耦合。白青波等［14］提出固液比B（T）的概念，将其作为联系方程：
式中：Tf为土体的冻结温度；
b为固液比系数，砂土通常取值0.61，粉土取值0.47，黏土取值0.56。
文中水分场与温度场之间属强耦合，彼此相互影响。方程（15）以土体冻结温度Tf为分界点可分为2个阶段：当温度超过土体的冻结温度时，土体内不会发生冰水相变过程，此时温度变化与土体的热容以及导热系数有关；
当温度低于土体的冻结温度时，土体内水开始相变成冰。一方面，相应位置含水量减少，基质势发挥作用，土体内水分开始迁移，由液态水含量高的位置迁移至液态水含量低的位置。另一方面，冰水相变所产生的潜热也开始影响温度的变化速率，水分场与温度场之间的耦合开始发挥作用。
1.4 应力场控制方程
土体在冻结过程中，水分场重分布以及水相变成冰产生的冻胀造成土体内部应力重分布［15］。由水热两场耦合产生的体积应变作为土体内部的应力场，从而建立应力场控制方程。
图2 土体微单元热通量示意图Fig.2 Schematic diagram of heat flux in soil micro-unit
静力平衡方程：
几何方程：
物理方程：
在冻结过程中，冻结区内的水分产生了相变，水相变成冰体积相应增大9%。与此同时，在基质势的作用下，未冻结区的液态水向冻结区移动，则土体的应变可以用迁移的液态水含量来表示，即温度作用下由于水分迁移产生的应变εvf为：
式中：θ0为初始含水量；
Δθ为水分迁移量；
θu为未冻水含量；
n为初始孔隙比。
需要说明的是，土体冻融是一个非常复杂的物理过程，温度场以及水分场变化影响应力的重分布。反之，应力重分布又影响着土体温度以及体积含水量的变化［16］。在温度场作用下，冻结区的水分发生相变导致冻结区液态水含量降低，进而在基质势作用下未冻结区液态水迁移至冻结区，其填充土体内孔隙并产生进一步相变，水相变成冰体积增大，土体内部结构产生应力。在应力作用下土体发生应变，进而造成土体的冻胀。
本次数值模型参数取自涂志斌［17］一维土柱室内模型试验。参考试验用土为膨胀土，土体参数如表1所示。试验采用如图3所示装置，试验开始前将土体分层填充至装置中，最终形成直径14 cm、高40 cm的土柱，采用直线位移传感器测量试验过程中土体顶部的冻胀量，采用土壤温湿度一体传感器测量距离土柱顶部5 cm位置处的温度及含水量变化。对土柱顶部进行密封处理，将土柱置于多功能环境箱内，设置环境温度为10℃并维持24 h，以确保土样中的温度均匀分布，然后在试验装置外侧及底部包裹保温棉以达到单向冻结的目的。最后，将试验装置放于多功能环境箱内，控制环境温度为-10℃并维持144 h。
表1 试验用土基本物理参数表［17］Table 1 Table of basic physical parameters of test soil
由于在土柱底板以及侧边施加保温措施，本次数值模拟假定土柱四周及底部为绝热状态且完全侧限，土柱与外界之间的换热仅由顶端发生。由于温度，水分只沿轴向传递，且只发生轴向变形，因此，本次数值模拟可简化为一维土柱水-热-力耦合模型。
数值模拟基于有限元数值软件Comsol。为简便计算，定义冻土相对饱和度S为：
式中：θu为未冻水含率；
θs为饱和含水率，选择S为变量代替θ进行求解。
将上述水分场、温度场控制方程编写成图4（a）所示格式进行求解。其中，水分场控制方程为：
温度场控制方程为：
对于应力场采用Comsol内置固体力学模块进行求解。水分场边界条件设为零通量，温度场边界条件选用狄式边界条件。狄式边界条件可以是随时间变化的变量，也可以是恒定值，方程形式如图4（b）所示，初始温度为10℃，上边界温度为-10℃，下边界设置为绝热。模型计算中所需参数如表2所示。其中，土柱密度、初始体积含水量、弹性模量均根据参考试验取值，土体的导热系数以及热容参数根据《冻土地区建筑地基基础设计规范》（JGJ 118-2011）［18］建议值选取。为简化计算，对土柱进行对称建模，沿土样高度划分为450个域单元以及118个边界单元，有限元数值模型如图5所示。
图3 监测点布置示意图［17］（单位：cm）Fig.3 Layout diagram of monitoring points（Unit：cm）
图4 Comsol内置假设方程形式Fig.4 Equation form in Comsol
图5 有限元数值模型（单位：m）Fig.5 Finite element numerical model（Unit：m）
表2 数值模型参数表Table 2 Table of numerical model parameters
3.1 温度场
图6为距离土柱顶部5 cm位置处温度随冻结时间的变化。分析可知，温度变化曲线可分为2个阶段：在冻结时间12 h内，该位置处温度下降速度较快，由10℃迅速下降至土体的冻结温度（-1.33℃），这是由于土柱顶部位置温度梯度较大，在此阶段无冰晶出现，未发生冰水相变。待温度达到冻结温度后，冰晶成核迅速生长，土体内自由水发生相变并迅速冻结；
在冻结时间超过12 h后，由于水相变吸热，温度下降速度变缓，在递降阶段中，土体中的结合水也开始发生冻结。数据显示，温度场数值模拟结果与实测结果具有较好的吻合度，验证了温度场控制方程以及耦合的准确性。
图7显示了在冷端温度为-10℃时，不同制冷时间下一维土柱的温度场有限元数值模拟云图。随着制冷时间的增加，冻结面逐渐向下传递，在冻结12 h时，冻结面距土柱底部0.35 m处，与图6实测结果有较好的吻合度，在冻结144 h时，冻结面距土柱底部0.24 m。图8表示不同冻结时间下温度沿土柱高度分布。由图可知，冻结区与未冻结区之间形成了冻结锋面，且下移速度在开始的1～48 h内较快，这是由于在冻结初始阶段土
图6 距离土柱顶部5 cm处温度随冻结时间变化Fig.6 Temperature variation with freezing time at 5 cm from the top of soil column
图7 温度场云图Fig.7 The cloud map of temperature field
内温度变化幅度较大。由图可以看出土柱内的温度分布曲线中存在明显的拐点，这是在温度降至冻土的冻结温度时，冻土内液态水在该温度处发生物相改变，冻结区与未冻结区的温度分布出现了明显的差异。在冻结后期，冻结深度变化速度变缓。在冻结144 h后，温度沿土柱高度近似呈线性分布。
3.2 水分场
图9为距离土柱顶部5 cm位置处液态水含量随冻结时间的变化。由图可知，该位置处液态水含量随冻结时间不断降低。在冻结时间12 h内，液态水含量几乎无变化，这是由于该位置处的温度高于土体的冻结温度，液态水未发生物相改变。当冻结时间超过12 h后，温度降低至土体的冻结温度，液态水发生物相改变，液态水含量降低，且前期降低速率较快。与此同时，在基质势作用下液态水开始迁移，由土柱下部迁移至土柱上部。该位置处液态水在相变成冰的同时，又受到来自土柱下部液态水的补充，故后期液态水含量降低速率逐渐变缓。水分场数值模拟结果与实测结果具有较好的吻合度，并证明了水分场控制方程以及水-热耦合的准确性。
图10表示在不同冻结时间下土柱内含水量分布的数值模拟结果。由图分析，在土柱冻结过程中，冻结区含水量明显增加，在冻结区与未冻结区处含水量差异较大，表明未冻结区处的液态水不断向冻结锋面迁移。这是由于冻结区液态水在低温环境下发生相变成冰，导致液态水含量减少，在水分场控制方程中，液态水会在基质势的作用下由含量高的地方迁移至含量低的地方，这就造成在冻结锋面处形成鲜明的S形曲线［19-20］。冻结锋面在一维冻结的条件下随冻结时间逐渐下移，这与温度场的变化规律保持一致。
图8 不同冻结时间土柱温度分布图Fig.8 Temperature distribution map of soil column with different freezing time
图9 距离土柱顶部5 cm处液态水含量随冻结时间变化Fig.9 Water content variation with freezing time at 5 cm from the top of soil column
图10 不同冻结时间含水量分布图Fig.10 Distribution map of water content at different freezing time
3.3 冻胀位移
图11表示土体冻胀量随冻结时间的变化关系。为便于比较，将考虑多场耦合以及不考虑多场耦合冻胀量计算结果同样列于图11中。分析可知，在冻结30 h前，土体冻胀变形发展较快。这是因为冷端温度传递速度快，上部土体孔隙被冰迅速填充，进而发生冻胀。在冻结后期，冻胀变形速度变缓，最终冻胀量为5 mm。进一步分析可得，考虑多场耦合计算出数值模拟结果（即以式（19）计算的结果）与模型试验结果所表现的冻胀量变化趋势相似，且最终冻胀量均为5 mm，如图12所示。然而，未考虑多场耦合计算出的数值模拟结果虽然变化趋势与实测值相似，但计算出的冻胀量明显小于实测值。当不考虑多场耦合时，即温度场以及水分场对应力场无影响，此时的冻胀应变εvf为：
此时，冻胀应变仅仅为初始含水量条件下土体内水相变成冰所产生的应变，不包含温度场以及水分场作用下的水分迁移所导致的应变增加，故未考虑多场耦合所计算出的数值模拟结果远远小于实测值。该对比结果进一步验证了文中水热力三场耦合的准确性以及可靠性。
图11 冻胀量随冻结时间变化Fig.11 Relationship between the frost heave changes and freezing time
图12 冻结144 h后土柱冻胀位移云图Fig.12 Frost heave displacement cloud map of soil column after freezing for 144 h
文中基于质量守恒定律以及能量守恒定律对冻土水分场以及温度场控制方程进行系统的推导，并结合水热条件下冻土内部应力变化建立水热力三场耦合数学模型，通过建立一维非饱和土柱数值模型可得结论：
（1）在冻结过程中，冻结面不断下移，初始冻结时间下移速度较快，后期趋于平稳，最终土柱内温度随高度线性分布；
（2）土柱内水分在冻结过程中发生迁移，在冻结面处形成明显的S型曲线，未冻结区水分不断向冻结区迁移；
（3）土柱冻胀量随冻结时间不断增大，最终冻胀量约为5 mm，数值模拟结果与实测结果趋势变化一致，验证了该模型的有效性。
猜你喜欢土柱液态水冻土降雨条件下植物修复分层尾矿土壤重金属迁移的模拟分析有色金属科学与工程(2022年4期)2022-09-01基于微波辐射计的张掖地区水汽、液态水变化特征分析成都信息工程大学学报(2021年2期)2021-07-22Ka/Ku双波段毫米波雷达功率谱数据反演液态水含量方法研究成都信息工程大学学报(2021年2期)2021-07-22北极冻土在求救小哥白尼(趣味科学)(2020年7期)2020-05-22分层土壤的持水性能研究西南农业学报(2019年9期)2019-11-11冻土下的猛犸坟场小哥白尼(趣味科学)(2018年9期)2018-12-18零下温度的液态水军事文摘·科学少年(2018年4期)2018-09-29火星上发现第一个液态水湖奥秘(创新大赛)(2018年6期)2018-08-24不同化学浸取剂对土壤镉淋溶过程影响化工管理(2017年1期)2017-03-05不同灌水量对2种盐碱土的洗盐效果比较江苏农业科学(2016年7期)2016-10-20
