吐鲁番市典型灌区输水渠道冬季冰盖多场耦合下力学响应特性分析研究

李 鸣

(吐鲁番市高昌区水管总站，新疆 吐鲁番 838000)

1 概 述

为提升农业灌溉效率，我国许多地区都建有长距离输水渠道，调控水资源至农业生产用水[1-3]。但在西北等冬季高寒地区，冬季输水势必会形成河冰及较厚的冰盖，对渠道安全输水带来较大影响，因而开展输水渠道冰盖体研究具有重要作用[4-5]。已有许多农业水利工程师针对河冰等流冰体开展过渠道结冰预警等研究，还有一些学者通过现场安装监测设备，研究冰盖在升温过程中内部温度场以及应力变形等变化，探讨各参数变化规律，为渠道稳定性设计提供重要参考[6-9]。针对多物理场条件，以有限元数值软件作为计算手段，研究不同工况下冰盖体力学以及温度场特征，从数值分析角度解释冰盖体对渠道安全运营的威胁[10-12]，为渠道安全建设提供理论依据。

2 工程概况

吐鲁番市典型代表灌区内有5条河流，最大年流量达3.31 m3/s，区域内春冬季温差较大，年风速超过1.5 m/s，工程结构需耐冷热冻胀效应影响。灌区内水资源丰枯水期层次性显著，可供应水量之比为2∶1，其中地表水主要来源于5条河流，为流域内农业生产提供灌溉用水。区域内已建有长度超过240 km的输水渠道，均采用格宾石笼针对性衬砌防渗结构，其中干渠有9条，分别支援不同乡镇农业生产用水；支渠超过35条，渠道断面形式为梯形，渠底宽度为40 m，水深设计为4 m，上覆盖层冰盖体平均50 cm，年灌溉水资源超过24 000×104m3，满足2.933 3×104hm2农业土地用水需求。局部地区水资源利用效率较低，灌溉效率每年以2%速率下降，这主要受地下水资源近几年有节制调度使用，而灌溉用水中地下水资源又占比较大。另一方面，灌区输水渠道衬砌结构在低温冬季易形成冰盖体，当温度变化时，冰盖体与渠道之间互相作用，降低水资源在渠道内输送效率，而实质上输水渠道与冰盖体为多物理场耦合作用下(温度场与固体场)，其力学特征关乎水资源输送效率。

为准确评判输水渠道在冬季低温环境下运营效率，针对渠道开展多场耦合数值试验分析。根据现场地质踏勘发现，灌区内地势并不较高，处于盆地低高程区域，高差100～150 m，地质构造主要在东南部丘陵地带可见背斜等地质现象。灌区内钻孔取样资料表明，覆盖土层以人工杂填土、粉质壤土及砂砾土为主，其中人工杂填土在表面分布范围较广，主要为一些种植填土与碎石填土，密实性较差，析水性较强，分布厚度约为2.5 m；粉质壤土主要存在于东南部丘陵进入灌区内的支渠周围，是部分支渠的渠地基，中等承载力，防渗性较好，含水量为18%；砂砾土现场取样得知，最大粒径超过6 mm，级配较差，与下卧砂岩层界面错动接触，松散型较大，覆盖土层地质年代基本均属于第四系。岩层以二叠统砂岩以及角砾石为主，强风化作用，完整性较差。地下水位最深处达16 m，基岩内孔隙水分布较少。另一方面调查发现，冬季输水渠道常出现冰盖体，部分渠道区段内衬砌结构在冬季易受损坏，冰盖体在春季温度升高时转换成流体运动，冰盖体颗粒尺寸不超过20 mm，厚度不等，浅层冰盖体内部具有较多的气泡，常见的冰盖体微观示意图见图1，冰盖体产生的冰压力分布是渠道衬砌结构安全考虑的重要方面。

图1 冰盖体微观示意图

3 多场耦合力学响应理论

冰盖体是一种特殊状态下固体结构，本文在分析过程中认定冰盖体为各向同性天然材料，在冬季至春季温度变化过程中处于流变变形，而描述流变变形的重要手段即是蠕变本构模型。本文以非线性Burger模型作为其蠕变本构模型，其具体表述可用下式：

(1)

式中:cice指常数参数；dref、dice为冰盖体尺寸；Asinha为黏滞系数。

式(1)中，表述了冰盖体所在流变过程中发生的弹性变形与塑性不可逆变形，直至失稳破坏状态。常用的冰盖体蠕变耦合模型还有以冰体流变速率表述的流变本构模型，其表达式为：

(2)

式中：Anor为流变实系数；n为常数参数，取3；Q为冰盖体自身能量；T、R分别为物理场参数。

在输水渠道中，冰盖温度变化过程中产生的冰压力会对衬砌结构产生一定破坏影响，考虑此，应结合温度场-应力场开展渠道冰盖体力学特征分析，其中热力耦合下控制方程为：

(3)

其中应力场变形服从以下方程式：

(4)

冰盖体在耦合场中会出现长期蠕变变形与热变形，其总变形可用下列3项式子表述：

(5)

式(5)中，各向同性横向变形取值参数以矩阵形式表述为：

(6)

上述理论为输水渠道中冰盖体热力耦合下应力分析基本理论，其中还包括有流变应力，要求解上述材料的本构模型方程，在有限元数值软件中需要建立耦合场，本文引入COMSOL Multiphyisics有限元数值软件，通过多物理场的微分方程进行模拟研究环境，并求解应力变形等特征参数。而在水利工程中，常见的多物理场微分方程一般呈非线性，其边界条件及方程式可用下式表述：

(7)

n·(cu+αu-γ)=g-qu+hTμ

(8)

式中：u为物理变化中的因变量；ea、da、c、α、β、γ、f均为耦合场中特定物理参数。

本文计算时忽略渠道冻胀变形影响，且认为冰盖体受热效应影响仅在融点内，忽略流变变形对体积应变张量不产生影响。

4 力学特征分析

4.1 研究工况与建模

模拟工况上下冰盖体温差为10℃，表面温度为-10℃，流变时间为6 d，所选取渠道特征断面见图2。分析计算温度从低温至融点过程中，即升温过程中冰盖体应力响应特征以及渠道结构参数改变对冰盖体应力特征影响。

图2 渠道特征断面

4.2 多场耦合下冰盖体应力特征

根据COMSOL数值有限元软件获得图3所示冰盖体在不同时间段时温度场分布云图。从图3中可看出，冰盖体在流变过程中表面层温度逐渐增大，初始状态时处于-10℃，在流变时间为3 d时，温度增大1.5倍，达-4℃；在流变时间为第5 d时，表面层温度已达到0℃，即温度由最底层完成传递至冰盖体顶层。从整体温度场变化分布来看，在未完成温度传递之前，即第5 d之前，冰盖体整体温度分布呈从底部至顶部，逐渐减小，但在冰盖体完成自上而下的温度传递后，顶部与底部温度为最高，均已达到0℃。但在冰盖体中间较厚区域仍然存在负温，即冰盖体在长期流变过程中温度场分布并不是均匀状态，即使已完成温度热传递，这也表明温度场与应力场耦合效应下，冰盖体应力特征变化过程中，内部温度场的温差作用持续发生。

图3 不同时间段时温度场分布云图

图4为计算获得流变升温过程中冰盖体压应力分布云图。从图4中可看出，初始状态下拉应力为最大，达150 kPa，最大压应力为241 kPa，位于距离冰盖体顶部5 cm区域，而拉应力区域集中在该冰盖体顶部，即初始状态下冰盖体顶部实质上是出于收缩变形才致拉应力出现；随着升温与流变时间推移，顶部逐渐由受拉转变至受压，在第1 d时压应力约为50～100 kPa，最大压应力集中在底部区域，最大拉应力相比初始状态增大46.7%，达314 kPa；升温至第3 d时，压应力逐渐降低，拉应力值亦处于较低水平，仅为7.13 kPa，冰盖体整体上均处于膨胀状态，即压应力分布为主要占比；在第5 d时，冰盖体顶部亦均为压应力分布，但在底部还残留有拉应力，即在温度传递完成之后，整体上拉应力分布亦从顶部至底部转移。

图4 冰盖体压应力分布云图

图5为COMSOL在迭代求解升温过程中冰盖体应力特征分布时获得的最大静冰压力曲线。从图5中可看出，静冰压力在升温过程中约为20～140 kPa，其最大静冰压力为142.3 kPa，与文献[13]计算获得静冰压力范围具有一致性。另外，亦可看出静冰压力变化是分阶段，呈先增后减变化。

图5 静冰压力变化曲线

图6为计算出的蠕变有效位移特征分布云图。从图6中可看出，蠕变有效位移从初始状态至温度完全传递过程中，应变量级增大12个量级，即冰盖体在升温过程中会促进蠕变进展，此亦在文献[14]中得到印证，外界温度升高，岩石等固体材料流变速率会增大，位移值也会相对提高。另从位移在冰盖体上分布来看，主要集中在右侧靠近渠道边坡侧以及顶面靠近，即与渠道相接触截面上蠕变位移持续稳定分布， 在冰盖体顶部右侧区域第5 d相比第1 d位移升高1个量级。

图6 蠕变有效位移特征分布云图

4.3 渠道结构参数对冰盖体应力响应影响

由于渠道结构涉及参数众多，本文以其中渠道侧边坡系数与渠道底宽作为影响冰盖体应力分析，基于COMSOL数值软件设定不同的边坡系数研究工况，对比各工况之间冰盖体结果差异。

4.3.1 边坡系数影响

边坡系数设定为1.5、2.25、3.25、3.5共4个工况，渠道底宽为40 m，其他参数与前述一致，研究升温过程中冰盖体应力变化特征，结果见图7。

图7 极限冰压力变化曲线(边坡系数影响)

从图7中可看出，各边坡系数下极限冰压力变化曲线形态基本一致，均呈先增大后减小，从各边坡系数下极限冰压力水平来看，当边坡系数增大时，极限冰压力峰值均有所增大。边坡系数1.5时，极限冰压力峰值为508.3 kPa；而边坡系数为2.25时，增大9%，达553.53 kPa；但在边坡系数3.25后，极限冰压力峰值逐渐降低，其中边坡系数3.5时极限冰压力峰值相比2.25时降低5.3%，为524.3 kPa。由此表明，边坡系数影响冰盖体极限冰压力峰值亦是先增后减，即边坡系数对冰盖体极限冰压力促进作用具有拐点系数，当超过边坡拐点系数后，冰盖体极限冰压力会逐渐降低。从极限冰压力随升温时间变化曲线可看出，各边坡系数下极限冰压力在后期流变过程中处于一致性，即极限冰压力在二次减小阶段减小速率为一致，均为每小时减小冰压力2 kPa。

4.3.2 渠道底宽影响

渠道底宽是灌区输水渠道设计中重要尺寸，本文以渠道底宽10、20、40和50 m作为研究工况，其他参数均为一致，获得冰盖体极限冰压力变化曲线，见图8。

图8 极限冰压力变化曲线(渠道底宽影响)

从图8中可看出，随着渠道底宽增大，冰盖体极限冰压力水平逐渐减小，渠道底宽10 m时冰盖体极限冰压力峰值为1 220.4 kPa，而底宽40和50 m时分别降低53.7%和62%，且各渠道底宽下冰盖体极限冰压力峰值基本均处于同一升温时刻，即第26 h。在冰盖体极限冰压力二次下降阶段中，各渠道底宽下冰盖体的极限冰压力曲线均为一致性，一次下降阶段中，渠道底宽愈小，则下降速率愈大，底宽10 m时冰盖体在一次下降阶段中每小时减小冰压力15.3 kPa，而底宽40 m时该参数为6.6 kPa。分析表明，渠道底宽尺寸愈小时，温度热效应在小断面内会具有较大变形约束，因而在冰盖体内部产生较大的冰压力，即呈现图8中现象，当完成温度热效应传递后，即第120 h后，此时温差传递对热效应减弱，故而冰盖体极限冰压力逐渐为一致性变化。

5 结 论

针对吐鲁番市典型灌区输水渠道冬季产生的冰盖体，利用COMSOL Multiphyisics有限元数值软件，开展热-力耦合下冰盖体力学响应分析，得到以下几点结论：

1) 研究了多场耦合下冰盖体升温过程中即使已完成温度热传递，温度场分布并不是均匀状态，表面层温度在5 d内完成热传递；冰盖体初始拉应力分布在底部，压应力约为50～100 kPa，热传递后拉应力至顶部，压应力占据冰盖体大部分区域，但升温过程中压应力为先增后减。

2) 获得了多场耦合下冰盖体静冰压力为20～140 kPa，最大静冰压力为142.3 kPa，静冰压力变化呈先增后减；流变有效位移在完成热传递后，增大12个量级，且有效位移位于冰盖体与渠道接触截面上。

3) 分析了渠道边坡系数对冰盖体应力响应影响特征，渠道边坡系数对冰盖体极限冰压力促进作用具有拐点系数2.25，当边坡系数超过2.25后，极限冰压力降低，边坡系数3.5时极限冰压力峰值相比2.25时降低5.3%；各边坡系数下极限冰压力在二次减小阶段减小速率为一致，均为每小时减小冰压力2 kPa。

4) 研究了渠道底宽影响冰盖体应力特征，各渠道底宽下冰盖体的极限冰压力曲线变化均为一致性，随底宽增大，冰盖体极限冰压力水平减小，极限冰压力峰值均处于第26 h；各渠道底宽的冰压力在二次下降阶段变化为相同状态。