季冻区渠基冻胀模拟分析

刘 超，郭妍秀

(1.黑龙江省江河流域保护中心，黑龙江 哈尔滨 150069；2.黑龙江省水利科学研究院，黑龙江 哈尔滨 150080)

黑龙江省是我国重要粮食生产基地，2022年产量约占全国总产量的11%，连续七年全国第一。粮食的高产归因于肥沃的黑土地以及充足的水资源，而水渠作为输送水分的主要载体成了灌溉农田的关键，据统计黑龙江省农业用水约占总用水量的4/5，用水效率为25%～40%[1]。同时，黑龙江省也是典型的季节性冻土区，恶劣的气候环境导致渠道发生隆起、滑塌、衬砌开裂等冻胀灾害，据某大型灌区的早年调查报告显示，全省约83%渠系工程受到环境影响发生冻胀破坏[2]。因此，渠道冻胀的研究成了我国学者们的热门选题。

渠道冻胀是渠基土在负温作用下形成温度梯度，导致土内水分向冻结锋面迁移形成冰透镜体，渠基向上隆起的现象。渠道的研究方法包括现场监测、室内试验研究和数值模拟分析，传统的现场监测及室内试验用时长、资源耗量大，数值模拟可有效降低时间和成本，已成为一种普遍的研究手段。1973年，RL Harlan[3]首次提出水热耦合的概念，并给出了相应的数学模型，为冻胀的数值模拟分析奠定了基础。土体冻胀的主要影响因素包含土质、水分、温度[4-6]，只有当三者满足一定要求时，才会发生冻胀，因此数值模拟也围绕着他们展开研究。刘旭东等[7]基于有限元方法探究渠基冻胀敏感性，发现影响冻胀的敏感性因素可依次排序为：渠道断面形式>渠基土冻胀系数>分缝位置>混凝土衬砌厚度>混凝土衬砌板弹性模量。李爽[8]考虑了混凝土衬砌的弹塑性以及冻土、衬砌接触的非线性，对衬砌与渠基土之间滑移和脱落的破坏过程进行数值模拟，研究表明：考虑非线性接触的模拟更加真实可靠，法向冻胀力和切向冻胀力较未考虑非线性接触的模型减小一半。杜民瑞[9]基于流传热公式分别对冬季停水、输水渠道进行温度场、应力场分析，发现输水渠道的边坡法向应力大于停水渠道。陈瑞考[10]基于水热力耦合方程以及有限元方法，建立关于THM的三场耦合渠道模型，针对梯形渠道建立水热耦合模型，分析冬季输水期间渠基内部的水分对冻胀的影响，探讨输水与停水渠道之间的差异。张海晨[11]在考虑高地下水位的基础上，利用有限元软件建立混凝土衬砌模型，分析渠基不同位置处的冻胀变形情况，发现渠坡的冻胀变形较渠底更显著，且最大冻胀量出现在渠坡1/3处。王羿等[12]考虑气温、渠基水分场、地下含水边界和土体渗流特性，构建了关于渠槽形状耦合的渠基土不均匀冻胀的水-热-力耦合数学模型。

目前，有限元法在渠道领域已经得到了普遍应用，但该方法对设备、存储量要求高、无法解决奇异性问题。有限差分法(FDM)妥善地克服了这一缺点，基于能量守恒方程、傅里叶定律等基本理论，利用网格节点上的函数值差商代替控制方程中的导数进行离散，建立以网格节点上的值为因变量的代数方程组的模型[13-15]。本文以黑龙江省松干渠道为原型，通过室内试验获取冻胀模拟的关键参数-热力学参数、冻胀系数、力学参数，采用有限差分方法建立封闭式渠道模型，施加当地冬季温度曲线模拟季冻区渠基冻胀全过程，分析渠基的冻胀特性及一般规律。

1 渠基模型的建立

基于能量平衡方程和傅里叶热传导定律的计算原理，利用有限差分法建立松干渠道渠基模型，模拟冬季停水渠道的冻胀全过程。

1.1 模型尺寸及网格划分

已知该渠基总尺寸为20.00 m×20.00 m×10.15 m，渠顶宽2 m，渠底宽4 m，横断面渠槽深4 m，坡比为1∶1.5，渠底基土厚6 m，纵向总长度为20 m，基土上部为0.15 m厚的衬砌板。将渠基的X轴、Y轴、Z轴分别划分成100、100、51个网格，共510 000个网格，如图1所示。

图1 渠道模型的网格划分

1.2 模型监测点布置

为了实时监测渠基在冻结过程中的温度场、位移场的变化，分别在渠底、渠坡中间的不同深度处设置监测点，共设置24个，其中监测点1#、监测点11#分别位于渠底、渠坡表面，监测点的具体分布情况见表1、表2。

表1 渠底监测点布置情况 m

表2 渠坡监测布置情况 m

1.3 模型的力学初始、边界条件

由于该模型渠道冬季停水，忽略雪荷载等因素的影响，只考虑渠道本身自重，渠基底部和侧向受到来自周边土的挤压，因此模型的侧向及底部分别受到水平约束和竖向约束限制其位移，上部为自由表面。

渠基模型的初始应力场采用弹塑性变形方法，首先将抗拉强度和黏聚力设置为最大值以防止模型发生屈服破坏，待计算完成后将二者恢复成初始数值进行运算，直至达到平衡状态，渠基模型的初始应力场如图2所示。

图2 初始应力场

1.4 模型参数

渠道模拟的关键参数包含：热力学参数、力学参数、冻胀系数，可分别通过导热系数试验、三轴试验、冻胀试验获取。上述试验皆选取含水率为15%的渠基土作为试验土样。

热力学参数依托于导热系数测定仪，试样制备完成后，将导热系数测定仪插入试样中部并设置不同温度。力学参数的获取方法是将含水率为15%的土样制备成Φ100 mm×H50 mm的试样，利用GDS三轴仪进行三轴UU试验(UU表示不固结不排水)。

(1)

式中：η为修正后的冻胀系数，℃-1；η0为初始冻胀系数，℃-1；Ls为试验测得的冻胀增量，mm；L为模拟所得冻胀量，mm。

查阅文献得到衬砌板的力学参数、热力学参数及冻胀系数[16]，结合上述试验所得数据汇总成如表3所示的模型参数。

表3 模型参数r

1.5 渠基模型的温度场施加

温度场演变实质是土质温度分布不均匀导致的热量传输过程。本次模拟采用显式的求解算法，待初始静力达到平衡，对渠基表面施加温度，运用各向同性传导方式模拟季冻区冬季渠基土冻胀全过程。应用FLAC3D进行冻胀模拟计算时，作出以下假设：

(1)土体温度场的热应力分布仅与温度有关。

(2)流体的分布不受结构应力的影响。

(3)可将土体冻结过程中的某一单位时间段内的冻胀率、导热系数视为常量。

1.6 气温变化的选取

由于黑龙江省从11月份开始发生冻结，选取大庆市肇源县2019年11月1日—2020年3月19日的气温作为冻胀模拟时间段。对气温曲线进行离散化处理，以每连续10 d的平均气温作为一组计算气温，共14组。该地区日平均气温和每组的计算气温如图3所示。

图3 模拟时间段的每组计算气温

1.7 温度边界、初始条件

本文所模拟的渠道上边界采用对流热边界条件，即流体温度随时间的变化情况。根据黑龙江水利科学研究院现场实际监测数据测得冻结前的渠基表面温度为3.1 ℃，开始冻结后的温度按图3的每组计算气温施加；渠基底部采用温度边界条件，施加恒定温度12 ℃保持不变；渠基侧面采用热流边界条件，热流密度为零。

2 渠基冬季冻胀模拟结果分析

利用构建的数值模型，模拟季冻区渠基在冬季冻胀的全过程，监测并分析渠底与渠坡的温度场、位移场变化情况。

2.1 冻结全过程中温度场变化规律

从渠基在整个冬季的温度场变化云图可以发现：渠坡和渠底表层土温度变化显著，温度梯度相对较大，渠基深处温度基本不受气温的影响；渠底下部的土层温度梯度近似平行。

根据温度场云图绘制出如图4所示的渠基温度场变化曲线。监测点1#、监测点11#表示渠基底部及坡面的最低温度，分别为-13.84 ℃、-14.80 ℃。渠基温度与深度呈负相关，0～60 d期间，监测点2#～4#、12#～14#的温度随着气温的下降快速降低，其余土层温度降低速率缓慢；61～100 d期间，基土温度随着气温的稳定而缓慢降低；100 d后气温回升，距离表面近的土层受到影响，温度有所升高，但依旧处于负温状态，远距离土层不受影响。

图4 渠基温度变化曲线

2.2 冻结过程中冻胀量变化规律

图5为渠底与渠坡的冻胀量随时间的变化情况，其中监测点1#、监测点11#分别位于渠底和渠坡表面，受温度影响效果最显著，监测结果代表该位置处的最终变形量，分别为2.83 cm、3.07 cm，远离渠基表面的监测点5#～10#、16#～24#，温度未达到初始冻结温度，不发生冻胀。

图5 渠基冻胀量随时间的变化情况

由图5可知，渠底与渠坡的冻胀变形趋势基本一致，前10 d渠基不发生变形，冻胀量为0；11～60 d期间，表面负温下移形成温度梯度，导致冻胀量快速增加，此时位于渠基表面的监测点1#和监测点11#的冻胀量增长最迅速，分别为1.94 cm、2.29 cm；2#～4#、12#～15#监测点由于受到土壤热阻的作用，冻胀量的增长速率减缓。61～100 d期间，随着温度梯度的减小，冻胀增长速率逐渐变得缓慢，直至增长至最大冻胀量，渠底与渠坡的最大冻胀量分别为2.83 cm、3.07 cm；之后气温回升但最高温度仍<0 ℃，未达到渠基融化温度，冻胀量趋于稳定。

3 结 论

本文通过室内试验获取冻胀模拟关键参数，基于有限差分法建立渠基模型，采用温度施加方法模拟松干渠道冬季冻胀全过程，探究渠基温度场、位移场演变规律。研究发现：渠基温度与渠基深度呈负相关，冻胀量变化包含快速冻结阶段、缓慢冻结阶段、稳定阶段，2月份冻胀量达到最大，渠底与渠坡的最大冻胀量分别为2.83 cm、3.07 cm。该方法极大降低了科研成本，为渠基冻胀研究提供了一定参考。