0.   前言

中国北方地区水资源矛盾突出，同时也是引调水工程分布密集的区域。以整体式（现浇混凝土，混凝土板间设置勾缝）或预制装拼式（以六棱块预制衬砌拼接，勾缝充填砂浆）衬砌渠道为主要结构形式的输水明渠，是区域引调水工程、灌溉工程的主要水资源配送结构形式^[1-2]。然而，因其结构单薄、刚度有限，对渠基土的冻胀变形既不能约束也不能很好适应，导致衬砌渠道破坏问题显著^[3]。

冻结状态下衬砌结构与土体的相互作用决定了渠道的破坏形式。近年来，结构与冻土界面的剪切特性正逐步获得关注，一些学者通过单元实验建立了接触面剪切特性的本构模型^[4-5]。然而，考虑冻土与结构物相互作用的数值化工作相对偏缓，目前大多数数值仿真研究的重点放在渠基土水-热-力耦合效应上，而将衬砌板与表层渠基土视为一个整体，或将衬砌混凝土简化为线弹性材料，这种处理方式就计算而言较为便利，但并不能反映低温条件下衬砌与表层渠基的相互作用所呈现的渠道冻胀特征。

实践表明，衬砌渠道因其造价低、输水效率高、施工便利易管理等特点，目前仍是经济可行，建管便利的输水结构形式。鉴于此，本文重点考察渠道衬砌与表面渠基土相互作用接触行为的数值模拟方法，兼顾计算便利，为衬砌渠道冻害处治手段提供技术支持。

1.   接触面黏滞行为的数值化方法

对于渠道衬砌与表面渠基土的接触行为，研究表明，当温度、表面粗糙度一定时，冻土与结构接触面抗剪强度与法向应力大体呈线性关系，符合莫尔-库仑强度准则事实上，由于法向冻结力的存在，接触面间的抗剪强度是动态的，与温度密切相关^[6]。Zhao等^[7]曾建议冻土接触面抗剪强度的表达式为

式中：${\varphi _{\text{s}}}$为静摩擦力对应的摩擦角分量，${\varphi _{\text{f}}}$为冻结力对应的摩擦角分量，c为黏聚力。文献[7]表明，剪切破坏发生时最终状态下的摩擦角随温度变化较小，而黏聚力随温度变化显著。

根据抗剪强度的特点，计算中可以将“摩擦行为”分离，并引入反映冻结作用的“黏滞行为”。沿着这一思路，引入一种线弹性“接触黏滞行为”^[8]，既法线方向产生黏滞力，同时又对切线方向的摩擦作用产生影响，与冻土接触面抗剪强度的作用特点吻合。线弹性的“接触黏滞行为”假设一切破坏都是以线弹性作为起始属性，通过将接触面的名义应力同名义应变联系起来组成弹性矩阵，名义应力t的表达式为

式中：t为启动黏滞行为之前的接触面即时应力；t⁰为单一的法向分离或切向分离行为开始前，界面接触压力的最大值；K为刚度；δ为分离值；角标n，s分别代表法向和切向。

在黏滞行为启动之后通过弱化材料的刚度来达到对破坏发展的定义，具体来说，可以假定一个损伤变量值D，这个值描述了材料整体的破坏程度并且可以捕捉其他与破坏有关的动态的力学影响。随着破坏的发展，D将会随着接触力的变化逐步从0增长到1。破坏启动后的法向应力可以定义如下：

式中：${\bar t_{\text{n}}}$为根据未破坏的即时应变δ_n计算出的法向应力分量；切向应力分项只需将式（4）中的n替换为s即可。根据有效分离值定义估计准则，损伤变量值D表达式为

式中：$\delta _{\text{m}}^{}$为有效分离值，上标f，max，0分别表示完全失效时、加载历史中峰值、以及损伤开始阶段。根据式（1）~（4），使用时需要定义峰值压力，以及接触面有效分离值$\delta _{\text{m}}^{\text{f}}$，$\delta _{\text{m}}^{\text{0}}$。除此以外，计算时还需给定初始应力参数$t_{\text{n}}^{\text{0}}$和$t_{\text{s}}^{\text{0}}$。

近年来，国内学者已完成了一些关于冻土与结构接触面间的剪切损伤模型的研究，可为上述参数的提取提供参考。例如文献[9]记录的数据，在一定温度水平和法向压力初始值条件下，得到了剪应力-法向位移、应力–应变等有意义的数据。由于样本缺乏，这里仅给出设置方法的建议：即按照式（1）的原则，分别对接触面设置切向摩擦接触和法向黏滞接触。摩擦接触即采用常见的莫尔-库仑本构模型；黏滞接触采用式（2）~（5）给出的定义，其中可将试验记录的界面未发生滑动或脱开前的抗剪强度记为$t_{\text{s}}^{\text{0}}$，此时的黏聚力记为初始冻结力t_n⁰。而将位移开始阶段和终止阶段分别记作$\delta _{\text{m}}^{\text{0}}$，$\delta _{\text{m}}^{\text{f}}$。表 1给出了设置示例，当法向压力为200 kPa时，给出不同负温下的接触面参数。其中摩擦行为采用μ标记，并设为定值，目的是突出黏滞行为的影响。将上述参数输入通用有限元软件中，得到了该接触模型下的剪切位移–剪应力关系，如图 1所示。可以看出，通过设置“接触黏滞行为”，可以反映接触面在峰值应力前剪应力随温度的下降而提升的特性。

表  1  不同负温下摩擦及黏滞行为的接触参数取值（法向压力200 kPa）

Table  1.  Contact parameters of friction and viscous behaviors at different negative temperatures (pressure of 200 kPa)

温度	接触面行为
	摩擦行为		黏滞行为
	$\mu$		$t_{\text{n}}^{\text{0}}$/kPa	$t_{\text{s}}^{\text{0}}$/kPa	$\delta _{\text{m}}^{\text{0}}$/mm	$\delta _{\text{m}}^{\text{f}}$/mm
-2℃	0.3		65	240	0.15	1.5
-4℃	0.3		85	270	0.15	2.0
-6℃	0.3		110	300	0.15	2.5

下载: 导出CSV 

| 显示表格

图  1  不同温度、法向应力下剪应力–位移关系的模拟

Figure  1.  Simulation of shear stress-displacement under pressure of 200 kPa and different temperature

下载: 全尺寸图片 幻灯片

2.   模型的建立

以新疆北疆地区实际渠道为原型建立计算模型。该渠道总长近130 km，跨度约19.1 m，断面尺寸远小于渠道长度，因此可将渠道冻胀作用视为平面应变问题。设渠基为均质土，取渠道一半为计算模型。考虑到工程实际冻深大约在1.0 m^[10]，计算时取表面以下3 m为底部边界，整个模型高6 m，宽9.5 m，AE=1 m，DF=5.4 m，HF=4 m，如图 2所示。衬砌材料参数按C20混凝土取值，其中导热系数按低温潮湿条件下取值，有关参数见表 2。

图  2  有限元计算模型示意图

Figure  2.  Schematic diagram of finite element model

下载: 全尺寸图片 幻灯片

表  2  混凝土材料参数

Table  2.  Parameters of concrete materials

密度/ (kg·m-3)	弹性 模量/Pa	泊松比	导热系数/ (W·m·k)	膨胀 系数
2400	2.4×1010	0.2	1.65	1.1×10-5

下载: 导出CSV 

| 显示表格

进一步地，在上图所示建模基础上，将“BD”、“DF”段按上述参数设置黏滞接触，并按μ=tan(0.75φ)设置摩擦接触，而“ED”段（即渠堤、渠底衬砌接触段）仅设置摩擦接触，此时模拟的结构为分段整体式衬砌，记为工况1；同时，将“BD”、“DF”仅设置摩擦工况，记为工况1-1用于对比。在此基础上，追加两种工况，一是整体式衬砌，即将“ED”段作为刚性连接，记为工况2；二是装拼式衬砌，渠堤衬砌按竖直方向每1m设置摩擦接触，记为工况3，用以考察不同衬砌结构形式的破坏特征。表面温度边界设为-4℃。底部温度设为10℃，右边界绝热。上述接触参数的设置均按表 1取值。

对于土体的冻结、膨胀行为，采用较为常见的“冷胀热缩”处理方法，即弱化水分迁移这一复杂过程，将土体冻胀作用简化为热膨胀的常规材料^[11]。有关软件中提供了一种本构模型，既包含弹性行为，也包含热膨胀行为：

式中$\alpha$(θ, $f_\beta ^{}$)为热膨胀系数；无量纲；θ为当前温度；${\theta _{^{\text{I}}}}$为初始温度；$f_\beta ^{}$为当前场变量值；$f_\beta ^{\text{I}}$为初始场变量值；${\theta _{^{\text{0}}}}$为参考温度。使用时软件根据用户输入的在某一温度下的膨胀系数$\alpha$（割线斜率）计算真实膨胀系数$\alpha '$（切线斜率）。根据${\text{d}}{\varepsilon _{{\text{th}}}}$=$\alpha '$(θ)dθ得到该温度下材料的应变值。为此，计算时需要明确材料的热膨胀系数$\alpha$，以及土体的名义弹性模量E。

热膨胀系数按渠基土自由冻胀率试验确定。参考蔡正银等^[12]的试验结果可知，渠基土在0~-5℃间冻胀率增长明显，而在-5℃后冻胀变形基本稳定，因此，计算中认为温度达-5℃后基土冻胀率η保持不变，而0~-5℃区间内的冻胀率按线性内插确定，根据膨胀系数α=η/Δθ，可求得膨胀系数α。取参考温度（相变温度）θ₀=0℃，不同温度下渠基冻土的热膨胀系数如表 3所示。

表  3  渠基土热膨胀系数设定

Table  3.  Setting of thermal expansion coefficient of foundation soil of canal

序号	温度/℃	热膨胀系数α
1	0	0
2	-1	-0.0016
3	-2	-0.0032
4	-3	-0.0049
5	-4	-0.0066
6	-5	-0.0085

下载: 导出CSV 

| 显示表格

另一方面，可以根据不同温度的渠基冻土单轴抗压强度，确定不同温度的冻土的弹性模量，如表 4所示。

表  4  渠基土冻结状态下名义弹性模量设定

Table  4.  Setting of nominal elastic modulus of foundation soil of canal in freezing state

温度/℃	0	-5	-10	-15	-20
弹性模量/MPa	5.5	16.7	43	59	181

下载: 导出CSV 

| 显示表格

计算时采用“先热后力”的计算策略，即先计算模型的温度场，再讲温度场计算结果加载至地应力平衡后的模型中进行应力–应变场计算。

3.   衬砌渠道冻胀破坏特征

3.1   温度场分布

终态温度场计算结果如图 3所示。渠坡及渠底表层的温度梯度大，随着深度的增大温度梯度越来越小，在接近下边界处温度等值线越接近于水平直线。模拟得到冻深分别为，渠底0.7 m，渠坡1.37 m。接触面的热传导热量按q=k(θ_A-θ_B)计算，其中θ_A为从面温度，θ_B为主面温度。从图 3中可以看出，主面混凝土衬砌的温度主要受边界条件影响，而衬砌下方基土受渠顶温度边界和接触导热界面共同影响，因此造成了不同位置处衬砌上下表面温度不一致（表现为渠顶表层基土温度较渠坡、渠底更低）的现象。

图  3  渠基温度场示意

Figure  3.  Temperature fields of model in terminal state

下载: 全尺寸图片 幻灯片

2.2   整体式衬砌结构的破坏特点

工况1、工况1-1的衬砌下表面所受法向冻胀位移及接触力分布如图 4所示。以沿衬砌外法线和方向为正。从看出工况1中渠顶表面、渠坡接近渠底1/3处以及坡脚附近底板一侧存在较大的法向接触力，而底板中部法向接触力较小。另一方面，工况1处渠坡板与渠底板发生了脱开，导致该点处法向接触力为0；渠坡板两端所受接触力为负，说明牵引黏滞力发挥了作用，且由于顶部温度较低，这种牵引作用更大；黏滞冻结力的存在使渠坡板两端受到约束，从而产生挤压作用，导致渠坡中部出现较大的冻胀力，但未出现滑移。渠底板靠近坡脚处呈现为冻胀力，而渠底板中部呈现为冻结力，导致渠底中部衬砌与基土共同上抬而在坡脚位置处相互挤压。设置黏结行为的工况1的模拟结果与实践相符^[13]，反映了渠道衬砌发生冻胀破坏的力学特性。相比较而言，工况1-1中渠坡板均表现为正的冻胀力；渠顶处的法向冻胀力与工况1-1相当，渠坡中部附近冻胀力偏小，而在坡脚处由于挤压作用出现了负的冻胀力。

图  4  衬砌间法向位移及接触力分布

Figure  4.  Distribution of normal displacement and contact force between linings

下载: 全尺寸图片 幻灯片

切向剪应力如图 5所示。以沿衬砌下表面渠底至渠顶为正方向，工况1中上方渠坡剪应力指向渠底，而下方渠坡剪应力指向渠顶，在渠坡板中部靠近1/3处形成分界，导致这一部位衬砌发生较大的法向位移，符合工程原型的破坏特征。相比较而言，工况1-1剪应力方向始终由渠底指向渠顶，且渠坡中部剪应力偏小，没有形成对渠道衬砌的约束作用，无法模拟渠道衬砌冻胀破坏的实际受力特征。

图  5  衬砌切向剪应力分布

Figure  5.  Distribution of tangential shear stress of linings

下载: 全尺寸图片 幻灯片

工况2——整体式衬砌，渠道法向冻胀位移以及衬砌切向剪应力分布如图 6所示。从冻胀位移分布中可以看出，衬砌最大位移出现在渠顶表面，为2.5 cm。渠道衬砌任意位置处均未脱开，渠底板对渠坡板的挤压作用消失，接触界面的黏滞力使渠道衬砌与基土一同变形，这种现象对渠道衬砌抵抗冻胀破坏反而是有利的。然而，根据接触面剪应力曲线，渠道坡脚处存在剪应力突变，因此该处为此类破坏情形的薄弱点。

图  6  整体式衬砌间法向位移及衬砌剪应力分布

Figure  6.  Distribution of normal displacement and shear stress monolithic linings

下载: 全尺寸图片 幻灯片

虽然在整体性较好的渠道中，接触界面的黏滞力有利于防止衬砌破坏，但黏滞力较大的接触面基土含水量往往较大，容易诱发水胀破坏^[14]，渠道衬砌防渗仍是防治衬砌破坏的必要措施。

3.3   装拼式衬砌结构的破坏特点

工况3——装拼式衬砌，衬砌法向冻胀位移如图 7所示。渠底板冻胀位移较小而渠坡板冻胀量较大。渠底板各位置的法向冻胀量基本相同，表现为整体上抬，但上抬位移较小，仅为0.35 cm；渠坡板最大法向位移位于渠顶，为4.0 cm，法向位移沿渠顶至渠底方向逐渐减小，并呈“阶梯”状，说明渠坡衬砌板间均出现了一定程度的相对位移。根据剪应力分布可知，衬砌间连接处两侧的接触压力和剪应力均存在突变，且突变幅值沿渠道自上而下衰减，而位于衬砌间接触点的接触压力和剪应力均为0。对于单个衬砌，靠近渠底一侧的接触力表现为冻结力，而另一侧表现为冻胀力。单个衬砌除受冻结力和冻胀力外，还受衬砌间接触摩擦力约束，最上方渠坡板顶端临空，抵抗冻胀能力最为薄弱，因此冻胀位移最大；下方各衬砌板同时受到两侧衬砌的摩擦约束，由于基土冻胀作用逐渐减弱，因此约束作用逐渐增强，最终呈现出“独立”地、“阶梯”状的位移衰减现象。

图  7  装拼式衬砌法向位移及衬砌剪应力分布

Figure  7.  Distribution of normal displacement and shear stress of fabricated linings

下载: 全尺寸图片 幻灯片

4.   结论

本文采用数值计算方法分析了寒冷地区典型衬砌渠道的冻胀破坏特征，得到以下2点结论。

（1）充分考虑渠道衬砌与渠基土间的冻结力，通过设置“接触黏滞行为”模拟冻结力，并以破坏标量值判定混凝土衬砌板与冻土之间是否会发生脱开现象；有效分离值、初始应力等参数可根据冻结接触面剪切试验确定

（2）对于渠底板和渠坡板间存在接触，而渠坡板整体性较好的衬砌渠道，破坏形式表现为渠坡板隆起破坏，最大隆起位于渠坡板靠近渠底的1/3处。对于整体式衬砌渠道，当冻结强度发挥作用时，渠道衬砌会与基土共同上抬；而当冻结强度未发挥作用时，渠道衬砌表现为整体上抬破坏；对于装拼式渠道，渠坡板上方的衬砌易发生坍塌、脱落破坏，衬砌位移的大小取决于基土的冻胀位移以及衬砌间的摩擦强度。

本文有关计算均在已有研究成果基础上得到。如能有足够的试验结果，得到对应负温条件下渠基土与衬砌接触面的剪切强度，可使研究对象的计算结果更加准确。