整体式小U形渠道混凝土衬砌结构冻胀破坏数值模拟研究

温 彤，王红雨

0 引言

混凝土衬砌是输水渠道普遍采用的防渗措施，也是提高农田灌溉用水效率的有效途径，然而，在我国北方高寒地区，冬季混凝土衬砌板因渠道基土冻胀而产生裂缝、断裂，甚至隆起破坏的现象较为普遍，水资源浪费严重，持续制约着灌区农业发展。文献资料显示，宁夏引黄灌区受干湿冻融交替等因素影响，渠系建筑物破损率达40%以上,农渠输水损失在50%左右[1]。渠道冻害防治已成为国内外学者长期探索的热点课题之一[2]。

针对渠道冻胀破坏问题，吴颖等[3]研究表明柔性渠道结构的抗冻胀效果优于混凝土渠道结构。张茹等[4]研究表明与各向同性模拟结果相比，横观各向同性模拟结果与实测成果更加吻合，从而验证了数值模型是合理的。高靖等[5]研究表明混凝土衬砌板的厚度或者倾角越大，冻胀不均匀性越小。

本文结合小型渠道衬砌结构原型试验，考虑渠道同一横剖面上阴坡和阳坡冻胀变形特征，提出了整体式小U形混凝土衬砌渠道冻胀数值计算模型，通过对试验段混凝土衬砌结构的温度、应变和应力模拟分析，结果与观测结果基本吻合。以期为高寒地区混凝土衬砌渠道抗冻胀设计提供参考依据。

1 基本假定及理论

1.1 基本假定

参考相关文献，结合宁夏引黄灌区小型渠道衬砌结构原型观测试验情况，对整体式小U形渠道混凝土衬砌结构在温度作用下的受力与变形模拟分析做以下基本假定：

1）将混凝土衬砌板看作各向同性线弹性材料，渠道基础冻土看作横观各向同性线弹性材料；

2）原型观测试验发现，小型渠道的冻结深度小于其地下水埋深，可以近似认为冻结期无地下水补给、无水分迁移；

3）鉴于原型观测试验区渠道基土为沙壤土，冻结温度接近0℃，所以渠道基土冻结基准温度取0℃，暂不考虑相变；

4）由于渠道属于细长结构，可不考虑沿跨度方向的温度差异，仅考虑渠道横断面的热传导，所以分析过程中可将衬砌渠道冻胀看作平面应变问题进行分析研究，即二维应变模拟分析。

1.2 基本理论

1.2.1 热传导方程

不考虑相变，假定整个冻结期、融沉期过程相对比较缓慢，放热、吸热为零，属于稳态热传导过程，因此渠道基土冻结稳态热传导方程为[6]：

式中:λx、λy、λz分别为渠道基础冻土沿 X、Y 和 Z 方向的导热系数，T 为温度，x、y、z分别为 X、Y、Z 方向的坐标。

边界条件：T(L，t)=TL，其中 t为时间、L 为冻结问题的边界[7]。

1.2.2 本构方程

将冻土看作横观各向同性材料，其本构方程如下所示[8]:

式中：T为温度，｛ε ｝为应变矩阵，｛σ ｝为应力矩阵，｛G ｝为刚度矩阵。具体表示如下:

式中：εx、εy和 εz为正应变，γyz、γzx和 γxy为切应变，σx、σy和 σz为正应力，τz、τzx和 τxy为切应力，其中 β1=（G11+G12）α1+G13α3），β3=2G13α1+G33α3，β2=β1，α1和 α3分别是 x（或 y）和 z 方向的线膨胀系数。

其中［S ］为柔度矩阵，［S］和［G ］互为逆矩阵，求解后得：

式中：S

而横观各向同性材料是正交各向异性材料的特例，还应该满足:

式中:x、y、z轴为弹性主应力方向，xy面为冻结面，z轴为温度梯度方向，且垂直于冻结面，Ex、Ey、Ez分别为渠道基础冻土沿x、y、z轴主方向的弹性模量，μij为渠道基础冻土的泊松比。

2 有限元模型及相应参数选取

2.1 原型渠道基本概况

原型观测试验地点为宁夏青铜峡市沙湖村灌区，选用整体式U形农渠横断面如图1；经过2016年、2017年两个冬季的冻融循环期原型观测，得到冻胀基本资料见表1，表中最大冻深是以冻结基准温度0℃为界实地测量取得。

图1 沙湖村整体式小U形农渠横断面图(单位:mm)

2.2 有限元模型及单元网格划分

结合前面的基本假定，应用ABAQUS软件进行混凝土衬砌结构热应力耦合分析时，将冻土与混凝土衬砌板看作一个整体考虑，冻胀分析作为平面应变问题进行研究，所以本文所提出的有限元模型是二维模型。

有限元模型计算时，冻结深度取值为原型观测结果，其中：阴坡最大冻深50 cm，渠底最大冻深40 cm，阳坡最大冻深30 cm，详见图2（衬砌板尺寸同图1）。

单元网格划分时，采用四边形作为有限元模型的单元形状进行网格划分，划分后网格数为13889个，节点数为14428个。

2.3 材料的参数选取

取渠道基土冻结后1.8373 W/m·℃[9]；在进行热应力耦合分析计算时，依据冻土热缩冷胀特性，冻土的热膨胀系数α取负值（α=-η/T）；鉴于渠道基土为沙壤土（冻结温度接近0℃），渠道基土冻结基准温度取0℃；冻胀率取表1中的最大值，即η=1.06%；冻土材料的弹性模量、泊松比、剪切模量等，参考相关文献[10-11]进行选取（见表 2）。

表2 冻土材料参数

混凝土材料参数参考文献[11]进行选取，其中：导热系数λ=2.3 W/m·℃，弹性模量 E=26 GPa，热膨胀系数 α=1.1×10-5，泊松比 μ=0.167。

3 数值计算

3.1 温度边界条件

温度:取原型试验段渠道衬砌结构表面各部位的平均最低温度作为有限元模型的上边界温度，其中阴坡取-8.5℃，渠底取-7.0℃，阳坡取-5.0℃，将其分别施加在混凝土衬砌板表面及其两侧渠顶地表面，下边界施加温度取冻结基准温度0℃。边界范围见图2。

图2 有限元模型边界范围与取值(单位:mm)

3.2 位移边界条件

位移:有限元模型下边界受两个方向的约束（μx=0，μy=0），上边界自由冻胀不受约束。

3.3 计算结果分析

3.3.1 温度场

经计算，阴坡板上面的土层冻结深度10 cm、20 cm、30 cm部位的温度值分别是-7.1℃、-5.7℃、-4.3℃，阴坡板下面的土层冻结深度10 cm、20 cm、30 cm部位的温度值分别是-6.4℃、-5.0℃、-3.5℃，渠底板下面的土层冻结深度10 cm、20 cm、30 cm部位的温度值分别是-5.0℃、-3.5℃、-2.1℃，阳坡板下面的土层冻结深度10 cm、20 cm、30 cm部位的温度值分别是-3.5℃、-2.1℃、-0.7℃，阳坡板上面的土层冻结深度10 cm、20 cm、30 cm部位的温度值分别是-3.5℃、-2.1℃、-0.7℃。计算结果与实测结果基本吻合。

3.3.2 变形场

表3 两种情况下的冻胀量对比

由表3可知，阴坡、渠底、阳坡的冻胀量计算值分别是5.27 mm、1.97 mm、2.63 mm，而实测值见表1，两者误差分别为0.5%、6.2%、6.1%，计算结果基本满足精度要求，所以两者基本吻合。

3.3.3 应力场

经计算，U形衬砌板弧底底部附近应力值最大，这个位置容易出现裂缝；而整体式U形渠道衬砌板的裂缝位置通常出现在渠底中部附近[12]，模拟的位置与实际的位置比较接近。

4 结语

（1）本文分析结果表明，基于有限元软件ABAQUS,根据当地渠道原位实测资料合理选用参数以及边界条件，对整体式小U形混凝土衬砌渠道进行数值模拟计算的方法是可行的。

（2）考虑渠道同一横剖面上阴坡和阳坡冻胀变形特征，提出了整体式小U形混凝土衬砌渠道冻胀数值模型，模拟结果与观测结果基本吻合。

（3）为了使模拟计算结果更符合实际情况，建议在后续的数值模拟研究中，通过大量试验数据验证，提出模拟计算参数选取方法及其计算结果修正系数，以便更好的为季冻区衬砌渠道抗冻胀设计提供参考依据。

[1]王茜,杨建全.宁夏灌区节水影响因素及节水途径研究[J].安徽农业科学,2013,41(24):10101-10103.

[2]建功.我国渠道防渗工程的冻害及其防治措施[J].防渗技术,2002,8(2):1-5.

[3]吴颖,周成,何宁,等.废旧轮胎柔性护岸减缓渠道冻胀数值模拟[J].地下空间与工程学报,2014,10(增2):1777-1781+1793.

[4]张茹,王正中,牟声远,等.基于横观各向同性冻土的U形渠道冻胀数值模拟[J].应用基础与工程科学学报，2010,18(5):773-783.

[5]高靖,田军仓,王斌.不同U型结构混凝土渠道冻胀数值模拟分析[J].灌溉排水学报，2015，34(4):38-42.

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

[7]闫长城,王正中,刘旭东,等.季节性冻土区玻璃钢防渗渠道抗冻胀性能初探[J].人民黄河,2011,33(3):140-142.

[8]王正中,袁驷,陈涛.冻土横观各向同性非线性本构模型的实验研究[J].岩土工程学报，2007(8):1215-1218.

[9]张洁,徐伟,钟建驰.润扬大桥南锚碇超深基坑围护冻结法施工温度场全过程的数值模拟[J].建筑施工,2004,26(1):15-17.

[10]李宁,徐彬,陈飞雄.冻土路基温度场、变形场和应力场的耦合分析[J].中国公路学报，2006,19(3):1-7.

[11]王正中,沙际德,蒋允静,等.正交各向异性冻土与建筑物相互作用的非线性有限元分析[J].土木工程学报，1999(3):55-60.

[12]李翠玲,唐少容,王红雨.整体式U形渠道混凝土衬砌结构冻胀破坏力学模型[J].水利水电技术，2015,46(4):148-152.