寒旱区梯形渠道衬砌冻胀破坏数值模拟

张海晨　李振　张茂林　姚志鹏　朱兴林　王铭岩

摘要：针对刚性衬砌渠道冻胀破坏问题，选取宁夏盐环定扬黄灌区衬砌渠道为研究对象，在考虑相变及水分迁移模型基础上，考虑高地下水位对冻胀的影响，建立考虑水分迁移及地下水影响的数学物理模型，并采用Comsol Mutiphysics有限元软件对梯形渠道混凝土衬砌进行冻胀数值模拟，分别从温度场、位移场和应力场对刚性衬砌渠道冻胀破坏等方面进行了研究分析。结果表明：阴坡、阳坡冻胀量较大，渠底较小，渠底中部、渠坡1/3坡板长度处冻胀量分别达到最大值。渠底中部的法向冻胀力较小，两边逐渐增大，阴坡、阳坡的法向冻胀力分布均匀;渠底的切向冻结力成线性分布。阴坡、阳坡靠近坡脚处切向冻结力较大，左、右两端坡顶处切向冻结力较小，渠底衬砌板和两侧衬砌板属于压弯组合变形构件。

关键词：渠道;衬砌;数值模拟;冻胀破坏;寒旱区

中图分类号：TV67

文献标志码：A

doi：10.3969/j .issn. 1000- 1379.2019.03.032

渠系工程在灌区工程中担任着主要的输水配水任务，但由于设计标准低，工程配套差，老化失修严重等问题，每年渠道渗漏导致的灌溉水损失占农业用水量的45%，占全国总用水量的33%.严重制约渠系工程的安全高效运行[1]。衬砌渠道防渗效果明显，其整体可减少渗漏损失的70% - 90%.是渠系工程中主要的防渗措施。我国冻土面积广阔，在北方灌区，每年冻胀融沉导致的渠道衬砌破坏现象屡见不鲜。衬砌渠道冻胀是影响渠系工程安全运行的一大因素，受到国内外相关科研工作者的广泛关注。Everett[2]通过试验探索，提出固、液表面的两相压力差是造成水分凝冰并引起土体冻胀的主要因素，即毛细理论：Miller[3]最早提出冻结缘的概念，解释了冰透镜体的形成过程，即第二冻胀理论：徐学祖等[4]通过研究冻土与融土的热力学特性及冻胀特性，得出了温度梯度与水分迁移的关系。这些学者的研究成果为渠道衬砌的冻胀研究提供了有益参考和理论支撑，鉴于试验统计的差异性及冻土冻胀实际情况的复杂性，使得研究结果存在一定的局限性，因此开展渠道衬砌冻胀的数值模拟分析就显得十分必要。Harlan[5]基于非饱和冻土的水分迁移理论，最早提出了水一热耦合模型：Konrad[6]根据冻土中自由水含水量与温度之间的关系，采用数值方法对冻土中的参数进行预报。然而，这些模型对于高地下水位对冻胀影响的考虑却较少[7]。鉴于此，本文针对渠道衬砌冻胀破坏现象，选取宁夏盐环定扬黄灌区衬砌渠道为研究对象，在考虑相变及水分迁移模型的基础上，考虑高地下水位对冻胀的影响，建立考虑水分迁移以及地下水影响的数学物理模型，并运用Comsol Mu-tiphysics有限元软件进行模拟，得到温度场、位移场、应力场的分布，并将数值模拟结果与试验结果进行对比，验证数值模拟的合理性与正确性，以期为相似环境的渠道衬砌冻胀破坏治理与防护提供理论依据。

1 数学物理模型的建立

1.1 建立模型的基本假设

渠道防渗衬砌工程是否产生冻胀破坏，与当地的土质、土的含水量、负温度及工程结构等因素有关[1]。渠基土冻结时，土体、水和冰之间的相互作用的微观结构及动态过程相当复杂，对于工程设计与优化计算来说，准确模拟整个冻胀过程是非常复杂且没有必要的[8]。为了便于实际应用，需对其进行适当简化，主要假设如下[9]。

（1）根据现场及室内试验研究，假设冻土是均匀连续各向同性体。

（2）水分对冻胀的影响表现在两方面：其一是土壤水分冻结释放相变潜热，影响渠基土温度场及冻深，从而影响冻胀量：其二是水分冻结形成分凝冰，体积增加，特别是地下水不断迁移结冰引起的膨胀。针对第一方面，建立考虑水分迁移热效应的冻土冻胀模型[10];对第二方面，参考《渠系工程抗冻胀设计规范》[11]中考虑地下水迁移影响的冻深一冻胀量关系，对冻土本构关系进行修正。

（3）根据试验研究，假定相变温度在同一种土中和同种外力条件下为常值，即暂取相变温度为0℃[12]。

（4）渠道为线状结构，沿水流方向长度远大于垂直水流方向宽度，可以作为二维平面应变问题进行建模。

（5）不考虑孔隙水压的影响。

1.2 考虑相变潜热的热传导方程

大量的试验结果证明土壤中未冻水在负温度梯度作用下由未冻区向冻结锋面迁移、集聚并冻结释放潜热[13]。本文将这部分相变潜热作为渠基土冻土物理场的热源看待，修正后的热传导公式为

1.3 考虑地下水影响的冻土本构方程

冻土在冻结过程中水冻结成冰，除原位水冻结膨胀外，还有从未冻结区向冻结区迁移的水分冻结成冰。工程中一般将原位水以及迁移水冻结成冰的体积膨胀以冻胀率η表示，其取值不仅与土质、冻深相关，也与地下水埋深有关[14]，《渠系工程抗冻胀设计规范》给出了地下水埋深与冻胀量、冻深的关系，见图1。

渠基土冻胀是由土壤孔隙中水分冻结成冰膨胀引起的，为便于工程应用，将渠基土视为冷胀热缩材料。平面应变问题中考虑地下水位的冻土本构方程为

以二维非稳态热传导方程建立的温度场模型与以本构方程建立的力学模型耦合，最终建立热力耦合分析模型。

2 渠道冻胀数值模拟模型及参数选取

2.1 渠道概况

该工程水源为黄河水，工程位于宁夏回族自治区盐环定揚黄灌区共用工程八干渠上段，为预制混凝土方形板衬砌梯形渠道，渠道断面见图2。边坡系数m=1.5，衬砌板厚度为6 cm，渠基土质为砂壤土，地下水埋深为渠道底板以下20 cm。工程所在区域属中温带干旱区，多年平均降水量290 mm.多年平均蒸发能力1 340 mm;l月（最冷月）平均气温-8.9℃，日照时数长，全年日照时数2 867.9 h。对该渠道进行冻胀数值仿真计算，并与现场实测数据对比。

2.2 有限元模型

为了验证上述数理方程的合理性与实用性，运用数学物理方程并采用Comsol Mutiphysics有限元软件模拟宁夏盐环定扬黄灌区渠道的冻胀过程，有限元模型是在原型渠道基础上的简化，基础从底板向下取250 cm.左、右边界各取75 cm。在边界条件中，上边界温度取原型渠道相应部位表面温度，阴坡、阳坡、渠底分别取-4.92、-4.56、-3.55℃，下边界温度取10℃，左、右边界近似隔热，其位移约束条件是x与z方向位移为零;渠底基土下边界加y向约束。梯形渠道有限元网格见图3。

2.3 参数选取

将混凝土衬砌板看成各向同性材料，弹性模量取2.4x104 MPa，将冻土与非冻土视为各向同性弹性体。未冻土弹性模量取15 MPa，冻土的弹性模量随着温度的改变而改变，取值见表1。其他材料计算参数及冻土的冻胀率与地下水埋深数据见表2、表3。

3 结果分析

3.1 温度场分析

由已构建模型对衬砌渠道的冻土温度场进行模拟，结果见图4。

由图4可知，接近于渠底的温度分布接近于一组平行的直线：0℃等温线以上，未冻水冻结导致的剧烈相变使等温线较密集：0℃等温线以下，等温线较稀疏，温差较大，容易造成渠基土的不均匀冻胀，阴坡冻深为98 cm，渠底冻深为55 cm，阳坡的冻深为91 cm。渠坡及渠底表层的温度变化快，即温度梯度大：而在离渠坡较远的深部，温度分布呈一组几乎平行的直线，基本上不受渠坡边界温度的影响。从温度场的分布和变化来看，三维边界、界面几何形状以及基土中热流对温度分布均有显著影响。从整体计算结果与前人的研究结果比较来看，模拟温度场与实测结果基本一致[15]，即上部冻深较大、底部冻深较小，阴坡冻深较大、阳坡冻深较小。土的最大冻胀率出现在地表下约1/3冻结深度处，而最大冻胀量发生在地表层。一般把冻胀量占总冻胀量的70% - 90%部分土体称为主冻胀区。主冻胀区一般分布于地表。八干渠渠道不同朝向的渠坡接受太阳辐射的能力不同，导致同一地区渠道阴、阳坡平均地温有一定的差异，同一区域内相同下垫面条件的太阳福射差异主要是渠道（地形因子）中的坡度和坡向不同引起的，总的来说阳坡的年平均地温高于阴坡的年平均地温，所以冻土往往具有明显的非对称性。

当冷却强度很大导致土内温度梯度很大时，表层冻结面迅速向下推移，土中水分来不及从下卧层向冻结面转移，就在原地冻结成冰。在这种条件下形成的冰，一般均匀散布于土孔隙中或土粒接触处。此类冻土一般无明显冻胀。当冷却强度小时，冻结面向下推移慢，甚至会因水结冰时放出的潜热的阻挡而较长时间停留在某一深度。此时，下卧层的迁移水流在克服了沿途的阻力后，还来得及到达冻结面形成冰层。有外部水源补给时，冻结面停留的时间越长，则形成的冰层越厚。这种冻土有明显的冻胀。

3.2 位移场分析

由已构建模型对衬砌渠道的冻土位移进行计算，并与试验结果进行比对.结果见图5。

由图5可知：

（1）模拟结果与模型试验结果分布规律一致，阴坡、阳坡冻胀量较大，渠底较小。原因是阴、阳坡冻深较大，发生冻胀的土体较多。因坡脚处衬砌板对基土的约束，故临近坡脚处的冻胀量很小。温度逐渐降低，土体中未冻结的水在负温度梯度的动力驱使下向冻结锋面迁移并最终冻结成冰，加上土体中原位水的冻结，造成土体的体积膨胀以及冻胀量变大，并最终造成了衬砌板的法向位移逐渐增大。

（2）在渠底中部、渠坡1/3长度处冻胀量分别达到最大值，其中：渠底2.16 cm，阴坡4.76 cm.阳坡5.01cm，与李安国模型计算结果基本一致[15]。阳坡冻胀量稍大于阴坡的，原因是在低温条件下，温度越低，冻结向下推进越快，水分越难向冻结边缘附近迁移，致使冻结区由水分迁移形成的冰的含量较少，因此阳坡冻胀量大于阴坡的。

3.3 应力场分析

3.3.1 法向冻胀力分析

渠基土负温条件下冻结膨胀会对渠道衬砌产生法向冻胀力以及切向冻结力。取垂直于衬砌板下表面的应力张量可得到法向冻胀力分布情况，见图6。

渠底两端由于衬砌板对于基土的约束以及阴坡、阳坡的衬砌板对渠底衬砌板的约束作用，渠道底部两端的法向冻胀力较大。渠底中部的法向冻胀力较小，两边的逐渐增大，其原因是基土与衬砌板底部的约束得到释放;阴坡、阳坡的法向冻胀力分布均匀，但在渠坡靠近坡脚处法向冻胀力较大。

3.3.2 切向冻结力分布

根据沿衬砌板下表面的应力张量可得到切向冻结力分布，见图7。

由图7可知，渠底的切向冻结力分布受两端约束，两端切向冻结力较大。渠底左侧受拉，右侧受压，因此渠底的切向冻结力成线性分布;阴坡、阳坡靠近坡脚切向冻结力较大，左、右两端坡顶处切向冻结力较小。

4 结论

（1）在考虑相变及水分迁移模型基础上，考虑高地下水位对冻胀的影响，利用Comsol Mutiphysics软件将冻土与渠道衬砌体作为一个整体进行热力耦合数值仿真计算，计算结果与试验结果基本一致，揭示了刚性衬砌渠道冻胀变形的基本规律及受力规律，可为渠道冻胀破坏研究与治理提供参考。

（2）温度场数值分析表明渠底的温度分布接近于一组平行的直线。0℃等温线以上，等温线较密集;0℃等温线以下，等温线较稀疏：渠坡及渠底表层的温度变化大，温度梯度大：深部温度分布不受渠坡边界温度的影响。上部冻深较大，底部冻深较小;阴坡冻深较大，阳坡冻深较小。

（3）位移场数值分析表明阴坡、阳坡冻胀量较大，渠底较小。渠底中部、渠坡1/3坡板長度处冻胀量分别达到最大值，阳坡冻胀量稍大于阴坡冻胀量，衬砌板冻胀量总体分布不均匀。

（4）应力场数值分析表明渠道底部两端的法向冻胀力较大。渠底中部的法向冻胀力较小，向两边逐渐增大，阴坡、阳坡的法向冻胀力分布均匀，都表现在渠坡靠近坡脚处法向冻胀力较大。渠底的切向冻结力成线性分布，左侧受拉，右侧受压，坡脚两端切向冻结力较大，阴坡、阳坡靠近坡脚处切向冻结力较大，左、右两端坡顶处切向冻结力较小。

参考文献：

[1] 王文杰，冻土水热力三场耦合的衬砌渠道冻胀数值模拟研究[D].杨凌：西北农林科技大学，2013：1-5.

[2] EVERETT D H.The Thermodynamics of Frost Damage toPorous Solids[J].Transactions of Faraday Society， 1961， 57（5）：1541-1551.

[3]MILLER R D.Freezing and Heaving of Saturated and Unsat-urated Soils[C]//Highway Research Record. Washington：Highway Research Board， 1972：1-11.

[4]徐学祖，邱维林，温度梯度诱导薄膜水迁移的冻胀机理[J].科学通报，1997，42（9）：956-959.

[5]HARLAN R L.Analysis of Coupled Heat-Fluid Transport inPartially Frozen Soil[J].Water Resources Research， 1973，9（5）：1314-1323.

[6]KONRAD J M.The Influence of Heat Extraction Rate in Freez-ingsoils[J].Cold Regions Science&Technology，1987， 14（2）：129-137.

[7] 刘月，王正中，王羿，等，考虑水分迁移及相变对温度场影响的渠道冻胀模型[J].农业工程學报，2016，32（17）：83-88.

[8] 朱志武，宁建国，马巍，土体冻融过程中水、热、力三场耦合本构问题及数值分析[J].工程力学，2007，24（5）：138-144.

[9] 陈立杰，王正中，刘旭东，等，高地下水位灌排渠道衬砌结构抗冻胀数值模拟[J].长江科学院院报，2009，26（9）：66-70.

[10] 刘月，考虑迁移水热效应的土体冻胀模型及双膜渠道衬砌防冻胀研究[D].杨凌：西北农林科技大学，2016：21-26.

[11] 中华人民共和国水利部，渠系工程抗冻胀设计规范：SL23-2006 [S].北京：中国水利水电出版社，2007：5-13.

[12]柯洁铭，杨平，冻土冻胀融沉的研究进展[J].南京林业大学学报（自然科学版），2004，28（4）：105-108.

[13]曾桂军，正冻土水分迁移及冻胀模型研究[D].北京：中国科学院大学，2015：28-31.

[14] 郭瑞，考虑分凝冰效应的饱和土冻胀模型及“适应型”渠道衬砌结构研究[D].杨凌：西北农林科技大学，2016：17-23.

[15] 李安国，陈瑞杰，渠道冻胀模拟试验及衬砌结构受力分析[J].水利与建筑工程学报，2000，6（1）：5-16.