冬季冻区混凝土输水渠道冻胀特性的数值模拟

杜民瑞，宋 玲，陈瑞考，魏 鹏

（石河子大学水利建筑工程学院，新疆 石河子 832003）

由于梯形混凝土衬砌渠道断面形式近似圆形，故其在工程实际当中梯形渠道被广泛利用，但是在冬季冻区普遍存在着的渠道衬砌层冻胀破坏现象，又给渠道的建设，乃至整个水利工程设计带来了更大的挑战。人们对于普通的衬砌渠道冻胀模型已经做出了大量的研究，RL Harlan[1]于1973首次提出了水热耦合的概念并得出了耦合模型，Loch[2]提出用克劳修斯–可拉贝龙方程来描述多孔介质冰和水相变过程的热动平衡，Shen Mu和Branko Ladanyi[3]在此基础上于1990年提出了水、热、力三场耦合问题，得出了简化后的耦合模型。Taylor and Luthin[4]基于Harlan模型提出便于求解的有限差分格式，Konrad[5]等通过实验数据以及理论分析，提出分凝势的概念，并把其定义为土体中水分迁移速率与温度梯度的比值。安维东[6]对冻土的温度、水分、应力三场耦合从理论上作了论述，并以某些水利工程的问题为实例进行数值模拟。董江伟[7]考虑了渠基土和所接触复合土工膜结构面之间的相互作用，对大断面刚柔混合衬砌渠道冻胀破坏进行有限元分析，但都是基于冬季停水状态下的研究，还有一类冬季运行的输水渠道，如电站的引水渠、给城市供水的防渗渠道以及跨流域调水工程的输水渠道，此类渠道在运行时不仅存在着水面出现冰盖问题，而且还存在着极大的渠道冻胀破坏的问题。宋玲[8-9]从温度含水率二元复合场的角度，创建了冬季输水渠道冻胀数学模型。对于冬季输水渠道而言，水位线以上是冻结区，图1为新疆玛纳斯河四级电站引水渠上游段横向裂缝及隆起破坏，可见边坡板在渠水位附近产生裂缝。

图1 引水渠道混凝土板局部破坏Fig.1 Local failure of concrete slab in canal

对于冬季输水渠道的冻胀问题可供参考的文献较少。本文利用有限元软件，以新疆玛纳斯河流域红山嘴电站四级电站引水渠道为研究对象，根据传热学以及冻土力学提出了在冻土中采用固液相容且互相传递热量的数学控制方程并对衬砌渠道的温度场、应力场进行数值模拟计算，分析和预测在冻结期渠基土温度场、应力场以及衬砌材料上冻胀力和冻胀变形的变化规律，旨在为冬季冻区输水混凝土渠道的设计、正常输运提出参考。

1 温度含水率二元场

1.1 停水渠道温度含水率二元耦合场

季节性冻土区的运河通道，横截面是一个复合耦合温度场和水分场。

图2 渠基土横截面上温度、含水率二元耦合场中特殊点的温度、含水率示意图Fig.2 Abridged general view of compound temperature and moisture content at special points in the binary variable field

如图2所示，耦合场的临气面是一弯曲的凹形线，而原状渠床的临气面则是一条直线。当其上空温度降至0℃以下时，两者水热耦合场是不同的。倘若A点处以下的土体含水率等于或者大于土的临界冻胀含水率时，该处附近的土层最终会产生冻胀现象和冻胀效应，但与原状渠床的冻胀不同。

1.2 输水衬砌渠道温度含水率二元耦合场

严寒地区的输水衬砌渠道，该横截面仍然是一个“耦合”场——水分场、温度场。存在差异的情况是，在冬季渠道供水状况下，渠堤、渠坡板，底板下的土体含水率梯度走势如图3。

图3 渠基土横截面上温度、含水率二元耦合场中特殊点的温度、含水率示意图Fig.3 Abridged general view of compound temperature and moisture content at special points in the binary variable field

渠坡处C点位置以下的土体含水率远远大于原始状态下的土体含水率，初始渗漏阶段的渗漏强度较大，渗水随着入渗时间的增长而减小，到最终达到稳定状态。这是因为在渗透初期湿润峰的行进距离较短，而水力梯度较大，最后随着水力梯度急剧减小并在后期逐渐趋于稳定造成的。

2 力学模型的建立

2.1 热传导分析

在渠基土冻结的热传导过程中，由于渠基土冻结过程缓慢，且历时较长，故可将其视为稳态传热的过程（忽略相变过程中释放的大量热量）。同时，将渠道冻胀问题简化为平面问题，二维平面应变问题基于稳定状态的接触热传导控制方程为：

式中：T为温度；λx、λy为冻土沿着x、y向的导热系数；A为计算的冻胀区域。

对于冻土，其导热系数随温度变化不大，而随土体含水量变化较大，含水量越小导热系数越小，热阻越大[11]。

2.2 冻胀机理及材料对流传热模型

当土体中的细粒土冻结时，冻结速率很小时，产生的温度梯度作用下促进水分流动。在这种力的作用下，水分不停地涌向冻结锋面并在该处堆积，体积膨胀。土体受土的自身构成、密实度、水分的影响。本文不考虑水分迁移这一复杂的过程，将冻土看作各向同性完全弹性材料，其弹性模量随着温度变化而变化[12]，将此问题看作流体外略平壁时产生对流传热问题，如图4所示。

图4 冬季输水渠道流体对流传热示意图Fig.4 diagram of convection heat transfer in canal with water in winter

考虑物性条件，假设流体速度均匀，且对流只发生在平行壁面方向，在黏性流体的对流传热经典理论中，在水流与壁面不发生滑移的前提下，确定水流与固体壁面之间仅凭借热传导产生热量的互换，因此可以得到对流传热系数表达式为：

式（2）中：λ为流体的导热系数，tw，t∞为来流的温度和壁面的温度。为贴壁处壁面法线方向上的流体温度变化率。

3 原型概况

本文仿真计算原型为新疆红山嘴水电站四级引水渠梯形渠断面，引水渠全长4.153 km，底宽3.0m，渠深4.5m，混凝级土强度为C20的衬砌防渗渠道，渠内水深3.5m，边、底板板厚均为0.2 m，渠道纵坡坡度1/1400，边坡坡度1/1.75。冻土层的最大厚度1.3 m，渠基土与混凝土板的摩擦系数为0.35。文献[14]通过小水壶多次多点取水测平均温度的方法，测得该渠道在冬季运行过程当中，沿流方向水温都在1.5～3.2℃范围内，因此渠内流体在通常情况下不会结冰。

3.1 参数计算

对于稳态的热传导，不同导热系数的材料将影响温度场的分布。混凝土的导热系数λc=1.65W/(m·℃)[15]。依据热传导方程，土壤冻结时导热系数λf=1.9870W/(m·℃)，几何模型需要考虑边界条件等因素建立，并且能简化原渠道模型。

水位面温度设为常温1.5℃，假定渠道模型无限长，转化为平面单元。模型的下边界无位移，x方向为自由边界，上边界两个方向均为自由边界。应力边界条件：左右边界剪应力为0，表1是在不同冰体积含量下，含水量一定时的冻土弹性模量[16]。线膨胀系数按照η/Tmin取值，η为冻胀率，Tmin为相应部位月平均表面温度最小值。

相应材料力学参数[17]的选取见表2。

表1 冻结砂土的等效模量Tab.1 Equivalent modulus of frozen sand

表2 材料力学参数Tab.2 Mechanical parameters of materials

3.2 计算结果

3.2.1 温度场

本文对整个模型采用三角形网格划分，共分成12328个网格单元，6256个节点。边界条件中，下边界温度取10℃，上边界温度取为原型渠道相应部位大气温度最小值，左右边界绝热。以温度场为模拟对象对冬季停水渠道和冬季输水渠道模拟结果进行对比。

图5为冬季停水、输水渠道等温线分布情况，由图5可知：

（1）渠道表层温度梯度大，而随着深度的增加，渠基土所受温度的影响也就越来越小（图5a）。

（2）冬季输水渠道冻结锋面深度不一，从水位面以上的渠坡到渠底的土体温度变化快即温度梯度大，而在远离渠坡处渠堤土的温度分布呈现为几组近似平行的直线，并且由于渠堤土与部分渠坡板下土体与大气负温直接接触，均不同程度产生了不均匀冻胀（图5b）。而与渠道内流体直接接触的渠坡板（水位线下部）和底板下部土体，因为流体的介入阻隔了这部分土体与大气负温的接触，并在热传递和衬砌体渗漏造成土体含水率增加的情况下，该部位土体未产生冻结（图5b）。

图5 冬季停水(a)、输水(b)渠道等温线分布图Fig.5 Isothermal distribution map of temperature canal water cut and with water in winter

3.2.2 应力场计算结果

根据数值计算得出的冬季停水渠道和输水渠道法向冻胀力分布图，绘制出沿衬砌板长的展开图见图6a。从图6a可以看出：

（1）冬季停水渠道底板中部受到的法向冻胀力较小，而在渠底板端部受到的负法向冻胀力较大。输水渠道冻胀力沿衬砌板分布方面产生了较为明显的差异，由于正负温的过渡变化，水位面以下渠坡板、底板的法向冻胀力为0，而水位面以上冻胀力与停水渠道类似。在负温状态下，冬季输水渠道渠坡板的法向冻胀力峰值要比冬季停水渠道渠坡板法向冻胀力峰值高约18%，切向冻胀力大小峰值两者相似。

（2）冬季停水渠道冻胀破坏往往发生在边坡板下1/3范围内，而冬季输水渠道冻胀破坏发生在边坡板中部1/3范围内。这是因为冬季引水渠道水位往往低于渠道设计水位的1/5。此水位对应的渠道水位线在边坡板中部以下1/3处的范围内波动，这个位置以上的渠坡板、渠坡土处于剧烈冻胀区，并且极易发生破坏。

（3）两模型法向冻胀力在水位面处产生了比较大的变化，在负温段，冬季输水混凝土渠道距离水位面1/3处的坡板处两侧受拉力，距渠顶1/3处两侧受压力，由于混凝土衬砌层受力性质抗压强度大抗拉强度低，故距离水位面1/3-1/2的渠坡板处易发生冻胀破坏，这也与实际工程情况相类似。

混凝土刚性衬砌板下表面受到的切向冻胀力分布如图6b所示，其中以指向延渠坡板拔出方向为正切向冻胀力。由图6b可知：

（1）冬季停水渠道阴、阳坡板以距离渠道底部1/3处分界，冻结力由分界点沿着混凝土衬砌板长度方向指向拔出方向，因此该部位极易产生破坏。这一特点被工程实践证明。

（2）冬季输水渠道切向冻胀力以渠水位面上部1/3处分界，更好揭示了经历了一个寒季之后输水渠道高位坡板破坏的现象。

4 结论和讨论

本文不仅首次从场的角度分析输水渠道横断面土体温度-含水率两场的耦合作用，并与冬季停水渠道耦合场进行对比分析，研究了输水渠道渠基土冻胀后的两场变化情况，而且首次考虑固液耦合模型对渠基土温度和衬砌层所受冻胀力的影响，并结合有限元软件对冬季输水渠道进行分析，明确了冬季输水渠道温度分布情况，研究了衬砌板下冻胀力的作用方式与大小，得到了以下结论：

（1）冬季输水渠道中与渠道内流体直接接触的渠坡板（水位线下部）和底板下部土体未产生冻结。

（2）考虑固液相容的模型比不考虑固液相容模型的冻结锋面的位置稍高。

（3）经过一个冬季之后，输水渠道发生冻胀破坏的部位要高于停水渠道。

（4）本文采用有限元软件对冬季输水渠道进行分析，结果基本能反映渠道冻胀变形的受力规律以及冻胀变形规律，故可为冬季输水混凝土衬砌渠道的设计与施工提供参考与建议。

（5）本文存在缺少准确的实测数据、具体的实验验证等不足，因此有待对此类渠道开展更加深入的研究。

[1]Harlan R L.Analysis of Coupled heat-fluid transport in partially frozen soil[J].Water Research，1973，9(5):1314-1323.

[2]Loch J P G.Thermodynamic Equilibrium Between Ice and Water in Porous Media[J].Soil Science，1978，126(2):77-80.

[3]Taylor G S，Luthin J N.A model for coupled heat and moisture transfer during soil freezing[J].Canadian Geotechnical Journal，1978，15(4):548-555.

[4]Shen M，Ladanyi B.Modelling of coupled heat，moisture and stress field in freezing soil[J].Cold Regions Science&Technology，1987，14(3):237-246.

[5]Konrad J M，Morgenstern N R.Prediction of frost heave in the laboratory during transient freezing[J].Canadian Geotechnical Journal，1982，19(3):250-259.

[6]安维东，陈肖柏，吴紫汪.渠道冻结时热质迁移的数值模拟[J].冰川冻土，1987，29(1)：35-46．An W D，Chen X B，Wu Z W.Numerical simulation of heat and mass transfer during channel freezing[J].Journal of Glaciology and Geocryology，1987，29(1):35-36

[7]董江伟，汤骅，姜海波，等.大断面刚柔混合衬砌渠道冻胀有限元分析[J].人民黄河，2016，38(11):149-152.Dong J W，Tang h，Jiang H B.et al.Finite Element Analysis of Rigid Mixed Lining Channels’Frost Heaving of Large Section[J].Yellow River，2016，38(11):149-152

[8]宋玲，欧阳辉，余书超.混凝土防渗渠道冬季输水运行中冻胀与抗冻胀力验算[J].农业工程学报，2015，31（18）:114-120.Song L，Ou Y H，Yu S C.Frozen heaving and capacity of frozen heaving resistance oftrapezoidalconcrete lining canal with water in winter[J].Transactions of the Chinese Society of Agricultural Engineering，2015，31（18）:114-120.

[9]宋玲.季节冻土区刚性衬砌渠道温度、应力及冻胀变形场的计算[J].石河子大学学报(自然科学版)，2008，26(6):763-767.Song L.Studies on computational method about temperature，stress and frost deformation of rigid lined canal in seasonal frost region[J].Journal of Shihezi University(Natural Science)，2008，26(6):763-767.

[10]康绍忠，粟晓玲，沈清林，等.石羊河流域水资源利用与节水农业发展模式的战略思考[J].水资源与水工程学报，2004，15(4):1-8.Kang S Z，Su X L，Shen Q L，et al.Strategy for water resources utilization and water-saving agriculture development models in Shiyanghe river basin[J].Journal of Water Resources&Water Engineering.2004，15(4):1-8.

[11]徐学祖，王家澄，张立新.冻土物理学[M].北京：科学出版社，2010.

[12]陈肖柏，刘建坤，刘鸿绪，等.土的冻结作用与地基[M].北京：科学出版社，2006.

[13]张辉，陈善年，徐益谦.理想流体对流传热问题的理论解[J].热科学与技术，2004(3):200-204.Zhang H，Chen S N，Xu Y Q.Analytical solutions to convective heat transfer of ideal fluid[J].Journal of Thermal Science and Technology，2004(3):200-204.

[14]黄酒林，宗全利，刘贞姬，等.高寒区引水渠道抽水融冰原型试验及分析[J].石河子大学学报(自然科学版)，2014,32（3）:392-396.Huang J L，Zong Q L，Liu Z J，et al.Mechanism and prototype experiment on pumping well water to melt Ice for diversion channel at high altitude and cold regions[J].Jou-rnal of Shihezi University(Natural Science)，2014,32(3):392-396.

[15]刘雄，宁建国，马巍.冻土地区水渠的温度场和应力场数值分析[J].冰川冻土，2005，27（6）:932-938.Liu X，Ning J G，Ma W.Numerical analyses of the temperature and stress Fields of channel in frozen soil regions[J].Journal of Glaciology and Geocryology，2005，27（6）:932-938.

[16]李爽，王正中，高兰兰，等.考虑混凝土衬砌板与冻土接触非线性的渠道冻胀数值模拟[J].水利学报，2014，45(4):497-503.Li S，Wang Z Z，Gao L L，et al.Numerical simulation of canal frost heaving considering nonlinear contact between concrete lining board and soil[J].Journal of Hydraulic Engineering，2014，45(4):497-503.

[17]宁建国，朱志武.含损伤的冻土本构模型及耦合问题数值分析[J].力学学报，2007，39(1):70-76Ning J G，Zhu Z W.Constitutive model of frozen soil with damage and numerical simulation of the coupled problem[J].Chinese Journal of Theoretical and Applied Mechanics，2007，39(1):70-76.