0 引言
中国北方季节性冻土区，水工建筑物经常受冻胀破坏。研究表明，因冻胀作用而破坏的水工建筑物中最多的是渠系构筑物和渠道衬砌[1]。其中，渠道衬砌因厚度小，自重轻而对冻胀作用十分敏感。加之冻胀的不均匀性，使得同一渠道同一断面不同部位的冻胀量不同，导致衬砌开裂，拱起或架空，加剧了渗漏。因此，有必要根据渠基土的冻胀量选择适宜的保温措施和渠道断面形式。王正中等[2]建立了梯形渠道混凝土衬砌抗冻胀破坏的力学模型，提出了衬砌板上各种复杂的冻结力、冻胀量等都可以表示为最大冻结力的函数，该模型简单合理。郭利霞等[3]根据梯形渠道砼衬砌冻胀变形及破坏特征，对其进行了有限元分析，提出了渠道冻胀的二维模型。王正中等[4]提出了弧底梯形渠道混凝土衬砌冻胀破坏的简化力学模型及设计方法，并给出了弧底的法向冻胀量随边坡系数的变化规律。宋清林等[5]提出了混凝土衬砌渠道保温防冻胀最佳结构型式，并基于热传导理论和热阻等效原理建立了相应的保温板厚度计算公式。聂欣岩等[6]在考虑渠道水力最佳断面，侧墙稳定性及土方开挖费用的基础上，设计出能使工程费用最小的梯形渠道边坡系数优选模型。郭靖等[7]通过有限元模拟软件得到了渠道不同位置处合理的保温板铺设厚度，然而，该模型未考虑渠道衬砌与渠基土之间的接触行为，与工程实际不符。为此，参考文献[8]，建立了考虑渠道衬砌与渠基土接触作用，混凝土及冻土塑性屈服的冻胀破坏模型。以黑龙江蛤蟆通灌区总干渠为原型渠道，应用该模型对不同保温板厚度下的渠道进行冻胀数值模拟，并建立了不同边坡系数的梯形渠道模型，分析其冻胀适应性能，进行梯形渠道断面边坡系数优化。
1 材料与方法
1.1 力学模型的建立
1.1.1 热分析
地表温度的周期变化，引起温度梯度作用下热量在地中运输，导致地面温度波向地下传播[9]。为使问题简单化，假设土体是均匀和各向同性的，由于渠基土的冻结是个缓慢的过程，认为此热传导为稳态热传导，并且忽略水分迁移带来的热量和水分相变释放的潜热。将渠道问题视为二维平面应变问题，则稳态热传导控制方程为：
式中：T为温度；λx、λy为冻土沿X、Y向的导热系数；A为渠道计算断面。
1.1.2 材料本构模型
土体冻结过程，实际上是土中温度的变化过程[10]。冬季，随着气温的降低，土层从地表开始向下冻结，土中水分转化为冰，引起土颗粒间的相对位移，使土体体积产生膨胀、地表升高。本文不考虑水分迁移和补给，仅将土体视为“冷胀热缩”的各向同性弹塑性材料，采用Drucker-Prager屈服准则[11]，冻土的本构方程为：
式中：α为线膨胀系数；η为冻胀率；T、T0为当前初始温度（℃）；εx、εy为X、Y方向的正应变；γxy为剪应变；σx、σy为X、Y方向的正应力（MPa）；τxy为剪应力（MPa）；E为弹性模量（MPa）；ν为泊松比；fβ为场变量。
混凝土衬砌的本构关系，是指在外力作用下，混凝土内部应力与应变之间的物理关系[12]。我国北方严寒地区，渠道混凝土衬砌常因冻胀变形发生塑性屈服，因此简单的线性准则并不能准确描述混凝土结构的破坏[13]。有限元模拟软件ABAQUS为混凝土材料提供了一种材料模型-塑性损失，可以有效的模拟混凝土的非弹性行为。
1.1.3 渠基土与衬砌接触面本构关系
本模型通过ABAQUS主菜单中的Interaction菜单为渠基土与衬砌定义接触[14]。渠基土与衬砌接触行为的非线性体现在切向应力和法向应力上。法向考虑渠基土与衬砌之间的法向冻结力[15]，当渠基土与衬砌间的法向拉应力超过二者之间的冻结力，则认为渠基土与衬砌脱开，不存在接触压力。切向考虑渠基土与衬砌之间存在摩擦力，当渠基土与衬砌间的切向应力超过这一摩擦力的值，则认为二者处于脱开状态，否则，认为二者黏结在一起；ABAQUS采用Column定律可计算出法向应力，渠基土与衬砌间的切向应力极限值为二者接触面的稳定冻结强度。以上在ABAQUS中设置接触面间的法向属性、切向属性及热传导率即可。
1.2 有限元模型及参数选取
1.2.1 渠道概况
以蛤蟆通灌区总干渠2+498～2+598处渠道（图1）为原型，蛤蟆通灌区位于黑龙江省宝清县，根据气象资料，最冷时段为11月—次年2月，该时段月平均气温最低为-20℃，多年平均气温2.2℃。渠基土以粉土为主，干体积质量1.44～1.60 g/cm3。地下水埋深大于4 m，最大冻深2 m。宋清林等[5]对其进行了冻胀观测，部分数据见表1、表2。
表1 原型渠道基本情况
表2 原型渠道保温板铺设方案
1.2.2 有限元模型及参数选取
采用ABAQUS软件进行渠道冻胀破坏热力耦合数值模拟，根据不同方案，设置混凝土衬砌板与基土接触非线性，或聚苯板与基土接触非线性，考虑冻土及混凝土塑性性质。有限元模型如图2所示，从渠顶到下边界距离为10 m，认为此处为多年温度不变层，渠顶到左右边界距离皆为2.5 m。渠道的温度上下边界条件采用第一类边界条件。上边界按表1的最低温度进行取值，下边界取多年平均气温2.2℃。土体左右和下边界为固端约束，上边界自由。
图1 蛤蟆通灌区总干渠梯形渠道尺寸（单位：cm）
图2 有限元模型Mesh
在稳定传热条件下，温度场的分布仅与材料的导热系数λ有关。查阅文献[9]，并根据土壤含水率，采用线性差值的方法得，冻土的导热系数为1.56 W/（m·K)，地表以下2 m非冻土，受外界影响小，综合考虑矿物质情况取导热系数为4.7 W/（m·K）。查文献[16]知，混凝土、聚苯板导热系数分别为1.65、0.041 W/（m·K）。
本模型通过线膨胀系数[17]表征渠基冻土冻胀变形能力。线膨胀系数为η/T，η为渠道各部位的冻胀率，T为月平均最低温度，由现场观测实验可得不同保温板厚度渠坡、渠底冻胀率。其他参数见表3—表5。
表3 冻土的弹塑性参数
表4 混凝土塑性损伤参数
表5 材料其他参数
2 结果与分析
2.1 不同方案结果分析
2.1.1 不同方案基土冻深对比
图3为不同方案基土冻深对比图。由图3可知，方案一的渠基土冻深远远大于方案二、三、四的渠基土冻深值，说明采取保温措施后，可以大幅度降低冻深。无论采取保温措施与否，渠道各部位的冻深分布依次为阴坡最大、渠底次之、阳坡最小，与实测结果一致。方案三、方案四最大冻深皆低于0.5 m，满足工程安全性要求。
2.1.2 不同方案法向冻胀量对比
由于渠道阴坡处混凝土衬砌最易发生破坏，故本文选取沿渠道展开的阴坡部分进行冻胀分析。由图4可知，保温试验段相较对比段法向冻胀量大大降低，这说明在我国北方寒冷地区对渠道施以保温防冻胀措施是很有必要的，可以大幅度降低冻胀对渠道衬砌的破坏。保温板厚度越厚，坡板上最大冻胀量出现的位置越靠近坡脚处。方案二、三、四结果表明，在阴坡靠近坡脚约1/3坡板长度处法向冻胀量达到最大，与实测结果分布规律相符，最大冻胀量分别为10.8、9.5、5.0 mm。《渠系工程抗冻胀设计规范》规定[18]，渠道衬砌材料为混凝土的梯形断面允许法向位移值需小于1 cm，方案三、方案四满足要求。
图3 不同方案基土冻深对比
图4 不同方案法向冻胀量对比
2.1.3 衬砌应力分布规律
建立的模型考虑基土与衬砌间接触力学行为，对比段法向冻胀量过大（最大约7 cm）,导致衬砌板与基土发生了脱空与滑移，不能产生法向冻胀力和切向冻结力，并且混凝土衬砌已进入塑性屈服阶段。对于保温试验段，混凝土衬砌板靠近渠顶处塑性应变为0，说明这部分未进入塑性区，在渠坡靠近坡脚处塑性应变接近于0，说明这部分衬砌与基土发生脱空与滑移，已经不存在接触行为，这与工程实际相一致。
2.1.4 不同方案工程经济性分析
实际工程中，在保证渠道工程安全性和可靠性的基础上，还必须考虑经济性的要求。尤其是保温所用聚苯乙烯成本较高，如果保温板厚度不够，则渠道不能满足抗冻胀设计要求；如果保温板厚度过大，那么工程造价也会过高。参考聚苯板市场价格，约为700元/m3。表6为不同保温板铺设方案经济性分析，由表6可知，保温板厚度越厚，相应的工程造价也就越高。在充分考虑工程可靠性、经济性的情况下，方案三是最优方案。
表6 不同保温板铺设方案经济性分析
2.2 考虑不同边坡系数的梯形渠道冻胀模拟
梯形渠道混凝土衬砌破坏的二大主要原因：一是混凝土衬砌结构抗压、抗弯能力较弱，二是由于基土冻胀产生的不均匀变形，加之衬砌板之间的相互约束，使得衬砌板易在法向冻胀力的作用下破坏。
基于上述模型，对对比段和最优方案不同边坡系数的梯形渠道进行数值模拟，进一步研究边坡系数对梯形渠道冻胀适应性能的影响。由图5、图6可知，对未采取任何措施的对比段，改变坡比对渠道最大法向冻胀量的影响不大。在阴坡上部，随着坡角变缓，法向冻胀量明显降低。对采取保温措施的方案三，受坡比变缓的影响，法向冻胀量最大值降低，坡比从1∶2.5变化到1∶3时，法向冻胀量最大值约降低20%。所以，边坡系数为3的梯形渠道的冻胀适应性能更佳。
然而，在实际工程中，梯形渠道边坡系数的改变，不仅影响混凝土衬砌的用量，还影响渠道土方的开挖量、聚苯板的使用量，所以边坡系数的改变对工程造价也有一定的影响[19]。在保证工程可靠性的基础上降低工程造价有实际意义。以本工程为例，对不同边坡系数的渠道进行经济性分析，初步计算每立方米渠道的工程造价。如表7所示，边坡系数为3的渠道工程造价最高，边坡系数为2的渠道工程造价最低。由图5可知，渠道边坡系从2变化到2.5时，渠道的法向冻胀量最大值变化不大，虽然边坡系数变化到3时，渠道的法向冻胀量最大值变化明显，但是相较高昂的造价，显然边坡系数为2更加合理。因此，综合考虑工程的可靠性、经济性，边坡系数为2的渠道更优。
图5 对比段不同坡比渠道法向冻胀量对比
图6 方案3不同坡比渠道法向冻胀量对比
表7 不同边坡系数渠道经济性分析
3 结论与讨论
1）建立考虑混凝土材料非线性特性及基土与衬砌间接触力学行为的梯形衬砌渠道冻胀破坏热力耦合数值模型，将此模型应用到黑龙江省蛤蟆通灌区渠道，结果表明，该模型能反应该地区的渠基土冻胀情况。
2）利用本文模型，对梯形衬砌渠道3种不同的保温措施进行模拟，并与实测结果对比，发现数值模拟结果符合实测数据规律。从冻深、冻胀量、衬砌应力分布的角度对3种保温措施进行分析，当保温板厚度为渠坡10 cm、渠底12 cm和渠坡12 cm、渠底14 cm时，均满足梯形衬砌渠道可靠性要求。在此基础上考虑工程的经济性要求，得出了梯形衬砌渠道最优保温板厚度，即渠坡10 cm，渠底12 cm。
3)对不同边坡系数的梯形衬砌渠道进行数值模拟，定量分析不同边坡系数对梯形衬砌渠道冻胀适应性的影响，结果表明边坡系数为3的渠道冻胀适应性能更佳。当梯形渠道保温板厚度特定，如取渠坡10 cm，渠底12 cm时，从工程可靠性和经济性角度出发，得到梯形渠道的最优边坡系数为2。
所建立的梯形衬砌渠道冻胀破坏热力耦合数值模型综合考虑混凝土材料非线性特性，渠基土与衬砌间接触力学行为，因此结果更加合理、准确。以黑龙江蛤蟆通灌区总干渠为原型渠道，应用该模型对不同保温板厚度下的梯形衬砌渠道进行冻胀数值模拟，并建立了不同边坡系数的梯形渠道模型，分析其冻胀适应性能，对梯形渠道断面边坡系数进行优化。研究结果可为黑龙江及其他相似地区提供参考。
我国北方寒冷地区，聚苯乙烯保温板防冻胀技术是削减渠道衬砌冻胀破坏，减少渗漏，提高渠系水利用率的一项重要措施，借鉴相关有限元方法[20]。利用有限元软件建立数值模型对渠道冻胀破坏进行模拟可为工程设计提供参考。然而衬砌渠道冻胀破坏本身是一个复杂的过程，涉及多方面因素。例如渠基土的土质和水分条件、渠道的走向、渠道衬砌材料及结构形式等。这些因素在一定程度上影响数值模拟过程中的参数选择，使得参数取值比较粗糙，这需要结合实验与理论二方面，进一步准确量化各参数的选择。
从渠道冻胀破坏机理来看，渠道冻胀破坏是温度、应力、水分三场相互作用的结果，是热力学、力学、物理化学的综合问题。而本文建立的考虑混凝土材料非线性特性及基土与衬砌间接触力学行为的梯形衬砌渠道冻胀破坏热力耦合数值模型，只考虑了温度和力二场耦合，忽略了渠基土中水分迁移和补给对渠道冻胀的影响，还需对此作深入的研究。
参考文献：
[1]马巍,王大雁.冻土力学[M].北京：科学出版社,2014：318-322.
[2]王正中.梯形渠道砼衬砌冻胀破坏的力学模型研究[J].农业工程学报,2004,20(3):24-29.
[3]郭利霞,王正中,李甲林,等.梯形与准梯形渠道冻胀有限元分析[J].节水灌溉,2007(4):44-47.
[4]王正中,李甲林,陈涛,等.弧底梯形渠道砼衬砌冻胀破坏的力学模型研究[J].农业工程学报,2008,24(1):18-23.
[5]宋清林,何武全,李根,等.混凝土衬砌渠道保温防冻胀技术研究[J].灌溉排水学报,2015,34(4):43-48.
[6]聂欣岩,侯毅凯.梯形渠道边坡系数优选模型研究与应用[J].中国水利,2012(8):37-38.
[7]郭婧,娄宗科,郭启胜.苯板保温在混凝土衬砌渠道中的应用及数值模拟[J].人民长江,2013,44(5):57-60.
[8]李爽,王正中,高兰兰,等.考虑混凝土衬砌板与冻土接触非线性的渠道冻胀数值模拟[J].水利学报,2014,45(4):497-503.
[9]徐祖,王家澄,张立新.冻土物理学[M].北京：科学出版社,2010：24-26.
[10]叶尔绍夫.冻土学原理:Principles of geocryology.第一册,冻土学的地球物理和化学基础[M].刘经仁，译.兰州：兰州大学出版社,2015：78-85.
[11]王先军,陈明祥,常晓林,等.Drucker-Prager系列屈服准则在稳定分析中的应用研究[J].岩土力学,2009,30(12):3 733-3 738.
[12]郭明.混凝土塑性损伤模型损伤因子研究及其应用[J].土木工程与管理学报,2011,(3):128-132,163.
[13]肖旻,李寿宁,贺兴宏,等.梯形渠道砼衬砌冻胀破坏力学分析[J].灌溉排水学报,2011,30(1):89-93.
[14]张建华,丁磊.ABAQUS基础入门与案例精通[M].北京：电子工业出版社,2012：71-89.
[15]孙杲辰,王正中,王文杰,等.梯形渠道砼衬砌体冻胀破坏断裂力学模型及应用[J].农业工程学报,2013,29(8):108-114.
[16]MARSHALLA L.The thermal properties of concrete[J].Building Science,1972,7(3):167-174.
[17]卢荣胜,费业泰.材料线膨胀系数的科学定义及应用[J].应用科学学报,1996(3):253-258.
[18]中华人民共和国水利部.渠系抗冻胀设计规范:SL23-2006[S].北京：中国水利水电出版社，2006.
[19]王开民,WANG KAIMIN.梯形渠道水力最佳断面边坡系数推求[J].电网与清洁能源,2000,16(2):9-11.
[20]刘倩,石东洋.双相滞热传导方程的一个非协调混合有限元方法[J].河南师范大学学报(自然科学版),2016,44(2):15-21.