冻胀反力系数在渠道衬砌冻胀弹性地基梁模型中的应用

李宗利，姚希望，李云波，吴正桥，肖帅鹏，刘士达

（1. 西北农林科技大学旱区农业水土工程教育部重点实验室，杨凌 712100；2. 西北农林科技大学水利与建筑工程学院，杨凌 712100；3. 中水北方勘测设计研究有限责任公司，天津 300222）

0 引 言

弹性地基梁理论因其反映基础和结构协调变形的优越性，已被诸多学者用于分析寒冷地区土体冻胀和上部结构之间的力学响应[1-5]。Rajani等[1-4]基于弹性地基梁理论分析管道在基础冻土冻胀时的力学响应。董建华等[5]基于Winkler弹性地基梁理论，研究了寒区锚杆与格构梁复合结构在冻胀作用下力学响应。相比于渠道衬砌冻胀材料力学模型[6-9]，弹性地基梁力学模型能反映出衬砌与冻土之间的协调变形[10]。肖旻等[11]基于冻胀力和冻胀强度成线性关系并结合弹性地基梁理论，建立了考虑冻土与结构相互作用的梯形渠道冻胀破坏弹性地基梁模型。李宗利等[12-13]基于自由冻胀量和弹性地基梁理论，结合梯形渠道衬砌冻胀变形特点，给出合理的边界条件，分别建立了基土均匀自由冻胀和不均匀自由冻胀时的梯形渠道衬砌冻胀弹性地基梁力学模型，计算结果与已有试验和数值模拟吻合较好。

将弹性地基梁理论引入到冻土冻胀和上部结构相互作用分析中，虽解决了考虑冻土与上部结构之间相互作用的方法问题，但如何合理确定基床系数仍是核心问题。本文用弹性地基梁理论分析冻土冻胀和衬砌板之间的相互作用时，将弹性地基梁理论中的基床系数称为冻胀反力系数。冻胀反力系数与基床系数存在一定相似性，但却有所不同：结构对地基施加一定的荷载，地基产生相应的压缩位移，从而地基给结构一定的反力，地基反力与压缩位移的比值就是基床系数；基土产生冻胀变形，其上部结构对冻胀变形存在约束，从而形成一定的冻胀反力，冻胀反力和被约束冻胀量的比值即冻胀反力系数。若直接用基床系数分析冻土冻胀和上部结构之间的力学响应有所不妥。同时，考虑到衬砌板下不同点的冻胀反力系数随被约束冻胀量的不同而发生变化，且每一点冻胀反力系数具体值又是被约束冻胀量的非线性函数，若采用常冻胀反力系数进行计算，计算结果和实际结果会有所偏差。

冻胀反力系数的合理取值影响到基于弹性地基梁理论分析冻胀问题的合理性。目前对基床系数的研究相对较多[14-15]，但对冻胀反力系数研究相对较少。同时，基于弹性地基梁理论在渠道衬砌冻胀力学响应上的分析[11-13]尚未考虑到冻土冻胀变形过程中的非线性。本文基于冻土三轴试验结果，建立考虑围压和温度的邓肯-张本构模型；参考室内三轴试验测定基床系数方法，基于数值模拟建立变冻胀反力系数计算式；基于有限差分法离散弹性地基梁平衡微分方程，由此建立变冻胀反力系数梯形渠道衬砌冻胀弹性地基梁力学模型。探究冻胀反力系数分别为变量与常量时在梯形渠道衬砌冻胀力学响应计算结果上的差异，以期为大型梯形渠道衬砌抗冻胀设计提供参考。

1 基于数值模拟的冻胀反力系数确定

1.1 基床系数

1.1.1 基本概念

基床系数是地基土在外力作用下产生单位变形时所需的应力，一般可表示为

式中k为基床系数，MPa/m；P为地基土所受应力，MPa；s为地基变形，m。

1.1.2 三轴试验确定法

现阶段直接通过现场试验确定地基土的基床系数方法尚不成熟，而室内三轴试验操作简单，可控性强。操作过程模拟现场K30原位平板荷载试验[16]。根据相关规范[16-17]，基床系数室内三轴试验测定法是将土样经饱和处理后，在有侧向静止土压力状态下进行排水固结，侧限围压应按下式计算：

式中σ1为轴向压力，MPa；σ3为侧向围压，MPa；γg为土体容重，kN/m3；hg为土壤所处的深度，m；K0为土体静止侧压力系数，可按下式计算：

式中ν为泊松比。

固结稳定后，控制围压增量Δσ3与主应力增量Δσ1的比值n为某一固定数值，得到Δσ1～Δh0曲线，求得初始切线斜率或某一割线斜率定义为基床系数k。

1.2 冻胀反力系数

冻土顶面存在衬砌等结构约束，当冻结达到稳定状态时，可能出现如图1所示几种情况。在上部结构约束作用下自由冻胀量不能被释放，部分冻胀量被约束，甚至冻土出现压缩情况，被约束冻胀量的大小直接影响着冻胀反力大小。本文根据被约束冻胀量和自由冻胀量之间的关系将冻土冻胀类型分为自由冻胀、部分约束冻胀、完全约束冻胀和超约束冻胀。

因冻土受力变形具有较强的非线性[18]，被约束冻胀量和冻胀反力并不是成正比的关系，如图2所示。为了建立冻胀反力Pf和被约束冻胀量y之间关系，本文定义冻胀反力Pf和被约束冻胀量y之比为冻胀反力系数kf，即图2中割线的斜率，如式（4）所示。

式中kf为冻胀反力系数，MPa/m；Pf为冻胀反力，MPa；y为被约束冻胀量，m。

本文定义初始冻胀反力kf0为图2曲线的初始切线斜率；完全约束冻胀反力系数kf2则定义为图1中完全约束冻胀情况下冻胀反力Pf2与被约束冻胀量y2（y2=Δh）的比值。

本文参考基床系数室内三轴试验确定法[16-17]，基于数值模拟对冻土的冻胀反力系数进行研究。

1.3 冻土数值模型建立

冻土冻胀是一个较为复杂的物理现象，涉及到温度传递、水分运移、冰水相变、水结冰后的体积膨胀引起的冻胀。通过数值模拟方法模拟冻土冻胀的全部过程较为困难。因此本文研究的重点是冻土冻结和冻胀达到某种稳定状态时冻土非线性变形特性对基床系数的影响，暂不考虑冻土冻结和冻胀过程中冻土和衬砌板之间的相互作用。本文以封闭系统饱和冻土为研究对象，忽略水分迁移，认为水分在原位发生完全冻结，温度对冻土力学性能的影响占主导作用。

1.3.1 温度传导方程

冻土冻结和冻胀达到稳定状态时可按稳态温度传导方程对温度场进行求解。稳态热传导的偏微分方程为

式中T为温度，℃；λ为冻土等效导热系数，W/(m·℃)；∇为拉普拉斯算子。

对于冻土冻结和冻胀的稳定状态，可认为冻土中的水全部冻结为冰。因此对于本文研究的封闭系统饱和冻土可认为是由冰和土颗粒两相组成。导热系数等效值可按两相含量来计算[19]，表达式为

式中λi和λs分别为冰和土颗粒的导热系数，W/(m·℃)；θi和θs分别是体积含冰量和体积土颗粒含量，m3/m3。

1.3.2 冻土本构模型

邓肯-张模型[18]最早是基于非冻土的变形特性提出的，后来一些学者[20-22]将其应用到冻土的变形研究中。本文采用邓肯-张模型来描述冻土在发生冻胀时的应力应变关系。

式中a和b为试验常数，εa为轴向应变。

冻土的应力张量S可以分解为球应力张量Tr(εs)和偏应力张量Dev(εs)

式中Tr(εs)为球应变张量；Dev(εs)为偏应变张量；δ为单位张量；K为体积模量，Pa；Gs为割线剪切模量，Pa；εs为应变张量。

邓肯-张模型的割线剪切模量Gs和体积模量K可按下式计算：

式中G为剪切模量，Pa；γ为剪切应变；qult为极限偏差应力，Pa；E0为初始弹性模量，Pa。

初始弹性模量E0和极限偏差应力qult与邓肯-张模型中a、b参数有关。当a、b参数与温度T和围压σ3c有关时，初始弹性模量E0和极限偏差应力qult可按下式计算：

式中T为温度，℃；σ3c为围压，Pa。

在考虑冻土冻胀时，应变张量εs可按下式计算：

式中ε为总应变张量；εf为冻胀应变张量。

本文研究的封闭系统的饱和冻土冻胀应变[23]εf为

式中θ为含冰量，m3/m3；ρw为水的密度，ρw=1 000 kg/m3；ρi为冰的密度，ρi=917 kg/m3。

1.3.3 数值模型参数取值及验证

温度场较为简单，仅涉及到冰和土颗粒的导热系数。冰的导热系数λi=2.22 W/(m·℃)，土颗粒的导热系数

λs=1.5 W/(m·℃)[19]。

温度和围压对冻土的邓肯-张模型参数有显著的影响，对已有试验数据[22]重新整理，建立考虑温度和围压的冻土邓肯-张本构模型，再将其带入COMSOL软件计算，与试验结果[22]进行对比，验证其合理性。

1）冻土三轴试验

文献[22]对冻结粉质黏土进行温度分别为-5、-10和-15 ℃，围压分别为0.6、1.0和1.4 MPa的三轴剪切试验。试样为Φ50 mm×100 mm的土样。先将土样置于真空饱和缸中10 h，使其充分饱和；再对其进行给定围压下排水固结；固结完成后进行给定温度下不少于24 h冻结；冻结结束进行三轴剪切试验，试验结果如图3所示。

2）邓肯-张模型参数确定及验证

对该试验结果进行重新处理，得到不同围压和不同温度下的邓肯-张模型参数，如表1所示。

对表1中的数据进行回归分析，可得

式中A1、A2、A3、A4、A5、A6、B1、B2、B3、B4、B5和B6为多项式系数，如表2所示。

表2 邓肯-张模型a、b参数拟合多项式系数Table 2 a, b parameter fitting polynomial coefficients of Duncan-Chang model

将式（16）～式（17）代入式（9）～式（13），再代入式（8），便可得到考虑围压σ3c和温度T的冻土邓肯-张本构模型。将其代入COMSOL软件对该试验[22]进行数值模拟，结果如图3所示，可以看出数值模拟结果和试验结果吻合较好，决定系数R2平均为0.97，说明本文建立的邓肯-张模型能够很好地反映文献[22]试验冻土样随温度和围压变化的力学行为。

1.4 基于数值模拟对冻土冻胀反力系数的测定

1.4.1 三轴试验数值模型

对于常规三轴试验，试样尺寸相对较小，试验结果只能代表冻土层某一个点的力学特性。本文认为渠道衬砌板所受到的冻胀力全部由冻土层冻胀变形后被约束冻胀量产生。为使测得的冻胀反力系数能综合反映出整个冻土层的力学行为，数值模拟模型以整个冻土层为研究对象，荷载边界依然采用规范中室内三轴试验测基床系数方法中的类似边界。

如图4所示，取高度1.5 m、直径0.5 m的冻土圆柱。土柱为封闭系统饱和土样，孔隙率约0.45。土柱上表面环境温度分别为-5、-10和-15 ℃，底部为冻深线所在处，可认为该处为0 ℃。如图4a所示，数值模拟模型中考虑土柱自重G对应力的影响，在土柱的侧面施加K0γgh的围压，模拟土柱冻胀前的受力状态。如图4b所示，参考基床系数室内三轴试验测定方法，在土柱冻结冻胀的同时，向土柱顶部施加σ1荷载，在侧面施加σ3荷载，逐级增加荷载σ1和σ3，保持增幅Δσ3=mΔσ1，m=0.1[24-25]；记录被约束冻胀量和σ1，此时σ1就是冻胀反力，如图5所示；根据式（4）和图2，图5中曲线的割线斜率即为冻胀反力系数，如图6所示。根据数值模拟计算结果，土样在自由冻胀时，垂直方向将会出现61.4 mm的自由冻胀量。

由图5和图6可以看出，相同冻结温度下，被约束冻胀量每增加1 cm，冻胀反力平均增加0.15 MPa，冻胀反力系数平均减小1.12 MPa/m。相同约束冻胀量下，冻结温度每降低5 ℃，冻胀反力平均增大0.21倍。采用双曲函数反映这种变化规律，方便工程应用。根据基床系数k取值方法，一般以初始切线斜率或某一割线斜率定义基床系数k[16-17]，若以此方法对冻胀反力系数kf取值，明显会使冻胀反力计算结果偏大。因此，本文认为冻胀反力系数随被约束冻胀量的改变而发生变化。

1.4.2 冻胀反力系数计算模型

王洪新等[26]认为基坑被动区土体的非线性弹簧地基模型满足双曲函数关系，得到不同变形的基床系数计算表达式。该成果对本文整理冻胀反力系数计算结果具有一定的借鉴意义。随着被约束冻胀量的增加，冻胀反力的增幅逐渐减小，这种规律可通过双曲线函数反映。因此，本文基于双曲函数模型对所得冻胀反力和冻胀反力系数结果进行拟合，其拟合公式如式（18）和式（19）所示，拟合参数如表3所示。

式中a′和b′为拟合参数。表3 为本研究土样所得到的结果，对于其他特定冻土试样可参考本文的方法进行测定。若条件允许可，通过物理试验更准确地对冻胀反力系数进行测定。

表3 冻胀反力系数双曲函数模型拟合参数Table 3 Fitting parameters of hyperbolic function model of frost heave reaction force coefficient

2 基于变冻胀反力系数的梯形渠道衬砌冻胀力学模型

弹性地基梁的基床系数为非常数时，通常很难给出解析解，只能转而寻找近似的数值解。蔡四维[27]首次将有限差分法用于求解弹性地基梁，用一组差分方程替代微分方程，求解线性代数方程组，使求解过程大大简化。本节采用有限差分法和牛顿法，通过MATLAB编程，建立变冻胀反力系数的梯形渠道衬砌冻胀弹性地基梁力学模型。

2.1 弹性地基梁基本微分方程

当弹性地基梁上分布荷载q(x) = 0时，弹性地基梁基本微分方程为

式中y为弹性地基梁的挠度，对于本文研究的渠道衬砌冻胀问题即为前述被约束冻胀量，m；E为渠道衬砌板弹性模量，Pa；I为渠道衬砌板截面惯性矩，m4。

如图7所示，根据有限差分法原理[28]，地基梁和其边界可以被等分成n+4段，共计n+5个节点。设梁挠曲线方程为y=f(x)（被约束冻胀量方程），在挠曲线上取等距离的5个点，间距为Δx，其编号分别为(i-2)、(i-1)、i、(i+1)和(i+2)。

根据有限差分法[28]和材料力学[29]相关原理，节点四阶微分和内力的差分方程分别为

式中M为弯矩，N·m；V为剪力，N；Δx为单段长度，m。

将式（18）代入式（20）再代入式（21），可得

式中α=(Δx)4/(EI)。

则图7中梁上的每一个节点均满足式（24）。

2.2 梯形渠道衬砌冻胀边界条件离散

2.2.1 边坡衬砌板

根据文献[12]边坡衬砌板在冻胀时边界条件，有

式中Ms为边坡衬砌板弯矩，N·m；Vs为边坡衬砌板剪力，N；Fs为边坡衬砌板在坡脚处受到渠底衬砌板的法向约束力，N；Ls为边坡衬砌板长度，m。

将式（25）～式（28）代入式（22）～式（23），则有

当x=0（i=0）时

当x=Ls（i=n）时

联立式（29）和式（30）可得

式中βs=2Fs(Δx)3/(EI)。

联立式（31）和式（32）可得

2.2.2 渠底衬砌板

根据文献[12]，渠底衬砌板在冻胀时边界条件为

式中Mb为渠底衬砌板弯矩，N·m；Vb为渠底衬砌板剪力，N；Fb为渠底衬砌板在坡脚处受到边坡衬砌板的法向约束力，N；Lb为渠底衬砌板长度，m。

和边坡衬砌板同理，可得

式中βb=2Fb(Δx)3/(EI)。

2.3 矩阵方程建立

2.3.1 边坡衬砌板

将式（33）～（36）带入式（24）中，图7中的每一个节点有

式（45）为非线性方程组，可写成下面矩阵形式

式中M s为边坡衬砌板系数矩阵；ys为边坡衬砌板挠度矩阵；βs为边坡衬砌板边界矩阵。各矩阵表达式为

2.3.2 渠底衬砌板

和边坡衬砌板同理，将式（41）～式（44）代入式（24），整理成矩阵形式，则有

式中M b为渠底衬砌板系数矩阵；yb为渠底衬砌板挠度矩阵；βb为渠底衬砌板边界矩阵。各矩阵表达式为

2.4 矩阵方程求解

式（46）和式（50）为非线性方程组，可以通过牛顿法进行计算[30]。牛顿法实际上是一种线性化方法，基本思想是将非线性方程组F(y)=0逐步归结为线性方程组的求解。

设方程组F(y)=0有近似根ky，将方程组F(y)在ky处展开，有

于是方程组F(y)=0可近似地表示为

这是一个线性方程组，记根为k+1y，则k+1y的计算公式为

将本文式（46）和式（50）代入式（56），则有

式（57）便可采用迭代进行计算，ky和k+1y分别为旧的和新的迭代计算结果；kJ为雅克比矩阵，按下式计算：

对于首次迭代计算

将式（59）代入式（57）进行迭代计算。通过判断k+1y和ky之差的一次范数是否达到给定精度10-6进行停止迭代，通过MATLAB编程实现这一计算过程。本文数值模型收敛速度快，精度高。如图8所示，在梁长度1 m、分段数100情况下进行计算，只需经过5次迭代计算，便可完成计算，且具有非常高的计算精度。

2.5 数值模型验证

文献[12]根据弹性地基梁短梁解析解，建立梯形渠道衬砌冻胀力学模型。基于该解析解对本文数值解进行验证。由于弯矩和冻胀反力可根据挠度（冻胀量）进行求解，因此只需对挠度（冻胀量）进行验证即可。对文献[12]中的工程算例采用解析解和数值解进行分别计算，结果如图9所示。可以看出本文基于数值模型计算出的结果和解析解基本吻合。

3 冻胀反力系数对梯形渠道衬砌冻胀弹性地基梁力学模型计算结果的影响

衬砌板长度与厚度比值体现了衬砌的抗弯刚度，反映了衬砌板对渠基冻土冻胀约束能力大小，影响冻胀量分布，而渠道衬砌板厚度取值范围不大，因此结合衬砌长度来分析冻土非线性变形特性影响规律，即衬砌刚度对渠道衬砌冻胀的影响。

本节以某寒区梯形衬砌渠道为例，探究冻胀反力系数对渠道衬砌冻胀的影响。渠道基土自由冻胀量5 cm，环境温度-15 ℃，衬砌板厚度10 cm，衬砌混凝土弹性模量24 GPa。冻胀反力系数为变量时，选取本文-15 ℃计算结果；冻胀反力系数为常量时，选取-15 ℃计算结果初始值。边坡衬砌板和渠底衬砌板计算结果如图10和图11所示。

3.1 边坡衬砌板

3.1.1 冻胀量

如图10a所示，常冻胀反力系数冻胀量计算结果比变冻胀反力系数冻胀量计算结果略大，但分布规律基本一致。当边坡衬砌板长度为1 m时，冻胀量分布接近直线，随着长度增加，冻胀量曲率逐渐增加。这主要是因为在自由冻胀量一定，边坡衬砌板坡脚约束不变的情况下，随着边坡衬砌板长度的增加其整体刚度逐渐变小。

3.1.2 冻胀反力

如图10b所示，冻胀反力系数分别为常量和变量时，计算结果在法向冻结力的分布和大小上基本相同，但在法向冻胀力上却有较大差异。变冻胀反力系数计算出的法向冻胀力明显比常冻胀反力结果小，同时距离坡脚越近，被约束冻胀量越大，计算结果差异越明显。本文对梯形渠道衬砌冻胀过程简化，认为在渠底衬砌板横向约束作用下，边坡衬砌板坡脚处的自由冻胀量完全被约束。自由冻胀量相同时，不同长度边坡衬砌板在坡脚处被约束冻胀量相同。由Winkler地基理论[28]可知，若地基梁上某点变形相同，则产生的反力相同。但由于冻胀反力系数的差异，常量冻胀反力系数计算出的最大冻胀反力是变量的1.43倍。

3.1.3 弯矩

由图10c可以看出变冻胀反力系数弯矩计算结果明显小于常冻胀反力系数计算结果；变冻胀反力系数计算弯矩最大值位置相对于常冻胀反力系数整体略远离坡脚。随着长度增加，弯矩最大值逐渐变大后趋于稳定。这主要是因为在渠基冻土在自由冻胀量一定时，随着边坡衬砌板长度增加，整体刚度逐渐变小，远离坡脚区域的衬砌板对渠基冻土冻胀约束较小，随基土冻胀一起发生上抬位移；在坡脚处由于渠底衬砌板对边坡衬砌板坡脚处的约束，使得该处产生较大的冻胀反力，但边坡衬砌板只有坡脚点处的冻胀受到来自渠底衬砌板的强迫约束，其影响范围相对有限，且当边坡衬砌板长度增加到一定程度时，影响范围基本不变；在被约束冻胀量一定的情况下，影响范围内产生的法向冻胀力也基本相同，因此当边坡衬砌板长度增加到一定值时，其弯矩最大值基本保持不变。常量冻胀反力系数计算出的弯矩最大值平均是变量的1.12倍。

3.2 渠底衬砌板

3.2.1 冻胀量

如图11a所示，在自由冻胀量一定的情况下，随着衬砌板长度的增大，中间部位被约束冻胀量逐渐减小，变冻胀反力系数和常冻胀反力系数的差异逐渐减小，导致随着渠底衬砌板长度的增加，变冻胀反力系数和常冻胀反力系数的冻胀量计算结果差异逐渐变小；中间部位冻胀量逐渐增大，但增幅逐渐变小。

3.2.2 冻胀反力

如图11b所示，当渠底衬砌板长度较小时，变冻胀反力系数的冻胀反力计算结果和常冻胀反力系数的冻胀反力计算结果相差较大，反之相差较小。由于在坡脚被约束冻胀量都是5 cm，因此不同长度的渠底衬砌板计算的冻胀反力在坡脚处相同，和边坡衬砌板原因相同。结合图11a在渠底衬砌板两边自由冻胀量被约束，中间渠道基土冻胀发生了一定的变形，释放了一定的冻胀力，因此渠底衬砌板冻胀反力分布为两边大，中间小。当渠底衬砌板长度为3和4 m时，在渠底衬砌板中部局部产生了不大于0.15 MPa的法向冻结力，这主要是由于衬砌板长度相对较长时，中间区域产生的冻胀量相对较大，使得衬砌板中间有远离基土的趋势，但渠道基土和衬砌板冻结在一起，因此，在中间区域产生了一定的法向冻结力。

3.2.3 弯矩

由图11c可以看出随着渠底衬砌板长度的增加，变冻胀反力系数的弯矩计算结果和常冻胀反力系数的弯矩计算结果差异逐渐变小，和冻胀反力出现这种现象的原因基本一致。当渠底衬砌板长度增加到一定程度时，在渠底衬砌板中部出现法向冻结力，如图11b所示，对中部弯矩有一定削弱作用。因此，当渠底衬砌板长度超过一定值时，随着长度增加，中部弯矩开始减小，弯矩最大值出现在两边。

4 结 论

1）渠道基土在冻胀过程中具有非线性变形特性，基于室内三轴试验确定基床系数方法是可行的。提出了双曲函数反映基床系数随约束冻胀位移的变化规律，方便工程应用。

2）探究了冻胀反力系数分别为变量与常量时，在梯形渠道衬砌冻胀力学响应计算结果上的差异。结果表明，对于边坡和渠底衬砌板，常量冻胀反力系数计算出的最大冻胀反力是变量的1.43倍，计算出的弯矩最大值平均是变量的1.12倍。因此在采用弹性地基梁理论分析渠道衬砌冻胀问题时，若冻胀反力系数采用为常量，不考虑冻土的非线性变形，会使得计算结果偏大。

3）本文所建立的冻胀反力系数计算数值模型未能很好考虑水分迁移，得到冻胀反力系数比实际值偏小，但不影响所建立的梯形渠道衬砌冻胀弹性地基梁力学模型的应用。同时，尚未考虑衬砌板与地基分离和衬砌板的冻缩对弯曲影响，计算结果和实际值存在偏差，未来模型仍需进一步完善。