开放系统下冻胀引起地表隆起变形预测的解析方法

doi:10.7522/j.issn.1000-0240.2020.0051

开放系统下冻胀引起地表隆起变形预测的解析方法

耿琳^,¹^,²^,³

1.哈尔滨工业大学土木工程学院, 黑龙江哈尔滨 150090

2.北京工业大学建筑工程学院, 北京 100124

3.中国科学院西北生态环境资源研究院冻土工程国家重点实验室, 甘肃兰州 730000

Analytic method for predicting uplift deformation of ground surface induced by frost heave in an open-system

GENG Lin^,¹^,²^,³

1.School of Civil Engineering，Harbin Institute of Technology，Harbin 150090，China

2.College of Architecture and Civil Engineering，Beijing University of Technology，Beijing 100124，China

3.State Key Laboratory of Frozen Soil Engineering，Northwest Institute of Eco-Environment and Resources，Chinese Academy of Sciences，Lanzhou 730000，China

编委: 周成林

收稿日期: 2019-05-14 修回日期: 2020-02-18

基金资助:

国家自然科学基金青年基金项目.  41702382
冻土工程国家重点实验室开放基金课题.  SKLFSE201709
深部岩土力学与地下工程国家重点实验室开放基金课题.  SKLGDUEK1807
国家自然科学基金重点项目.  41430634
国家重大科研仪器研制项目.  41627801

Received: 2019-05-14 Revised: 2020-02-18

作者简介 About authors

耿琳（1983-），女，黑龙江哈尔滨人，助理研究员，2017年在哈尔滨工业大学获得博士学位，从事岩土工程与冻土工程研究.E-mail：gl_hit@163.com. , E-mail：gl_hit@163.com

摘要

开放系统下土体冻胀引起土工结构变形致使其服役特性面临严峻挑战。基于土体冻胀非线性弹性力学模型，推导了计算半空间无限体中土单元在自重下冻胀的数学表达式，建立了水平场地土冻胀模拟方法，提出了考虑上部压载效应与位移边界变化的场地表面隆起位移及土体应力与应变的计算方法。在此基础上，利用开放系统下粉质黏土逐级降温冻胀室内试验，验证了冻胀引起地表隆起变形预测解析方法的可靠性。经验证，该变形预测解析方法可应用于寒区场地冻胀量的评估，进而为开放系统下土工结构物变形预测提供科学参考。

关键词： 冻胀 ; 地表隆起变形 ; 解析方法 ; 开放系统

Abstract

Considering the engineering design demand， a simple analytical method for calculating soil frost heave deformation of horizontal field under steady-state condition is established based on the nonlinear elastic mechanical model of soil frost heave. Mathematical expression to determine the frost heave of soil element in semi - infinite space body under dead weight is deduced. Also， the prediction methods of soil lifting displacement， soil stress and strain field under different overburden and displacement conditions of horizontal site are proposed. Furthermore， the step-freezing test of Harbin silty clay in open system is carried out， and the reliability of the analytical method for predicting the surface lifting deformation caused by frost heave is verified. Besides， the analytical method is further applied to the evaluation of the frozen site. It can be concluded that the frost heave amount can be used as a main quantitative index to evaluate the frost heave of the site. However， it is still necessary to consider the frost heave sensitivity， the thickness of the frozen soil layer， the initial overburden pressure and the surface temperature.

Keywords： frost heave ; uplift deformation of the ground surface ; analytical methods ; open-system

PDF (2825KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

耿琳. 开放系统下冻胀引起地表隆起变形预测的解析方法[J]. 冰川冻土, 2020, 42(2): 540-549 doi:10.7522/j.issn.1000-0240.2020.0051

GENG Lin. Analytic method for predicting uplift deformation of ground surface induced by frost heave in an open-system[J]. Journal of Glaciology and Geocryology, 2020, 42(2): 540-549 doi:10.7522/j.issn.1000-0240.2020.0051

0 引言

我国是世界第三冻土大国，多年冻土和季节冻土面积约占国土面积的2/3。国家“十三五”期间，计划陆续在东北和中西部地区建设一批高速公路和高速铁路等交通基础设施。同时，在“一带一路”的倡议中，拟着力建设一批重要的基础设施工程，而且这些工程多分布于多年冻土和季节冻土地区。这些工程势必遭受低温环境的作用，可能发生不同程度或类型的冻害，并以工程构（建）筑物基础冻胀引起不均匀变形的危害最突出，如青藏公路和青藏铁路路基病害、大坂山公路隧道围岩冻害、哈大高铁路基冻胀等，冻胀与冻害诱发工程灾害的维护费用高达数亿元^［1］。应该说，基础冻胀问题已显著制约寒区重要土木工程的建设与发展。

大量试验结果与工程实践表明，冻胀敏感性土的场地发生冻结时，土颗粒之间的连结作用增强，土体强度显著提高^［2-5］。可见，冻结作用下土体一般不易发生因土体强度丧失的稳定性问题。但是，冻胀敏感性土在环境温度和外界补水条件下产生的冻胀力或不均匀冻胀，势必会引发一系列基础不均匀变形问题^［6-8］。基础与其下部土体之间的冻胀力相互作用，土体冻胀受地表压载的影响对建（构）筑物的基础有一个向上的抬起作用，可能使得基础因发生不均匀变形而出现不同程度的破坏，甚至完全失效^［9-12］。因此，实际工程中准确进行压载作用下地表冻胀的隆起量，以及地表位移达到工程允许抬升量时在地表需要多大压载的预测，将有助于保证建筑物基础安全^［13-15］。这项工作的实施也可以为基础冻胀变形的防治提供必要的依据。

寒区冻胀敏感性土的冻结本质上是多孔介质水力 - 变形 - 温度 - 应力多物理场耦合问题。冻结过程中，土体是由土颗粒、未冻水、冰和空气组成的四相力学介质，各相介质之间发生耦合作用机制，致使冻土多物理场问题的求解过程十分复杂。作为冻土多场耦合研究的一个里程碑， Harlan^［16］率先开展了冻土水动力模型，以考虑相变潜热和水份流动对传热过程的影响。随后，国内外学者针对冻土水热耦合模型和热 - 水 - 力耦合模型，开展了一系列卓有成效的工作。何平等^［17］基于连续介质力学与热力学原理，提出了冻结过程土体三场耦合方程。Mu等^［18］将固体力学引入水 - 热耦合模型，形成能够考虑应力场和应变场的水 - 热 - 力模型。赖远明等^［19］考虑当冻结土体承受外载作用时，不仅会产生瞬时变形，也会产生随时间变化的变形，冻胀变形也是不可忽略的因素。李洪升等^［20］将冻土视为空间弹性体，建立了冻结过程中土体水分场、温度场、应力场三场耦合的数学模型。许强等^［21］基于水 - 热蠕变模型，实现了三场耦合的目的，得到更为完善的温度场方程和体积应变计算式。朱志武等^［22］根据传热学、渗流理论及冻土力学，提出了带相变的温度场、水分场和应力场耦合问题的数学和力学模型。李宁等^［23］和毛雪松等^［24］利用正冻土三场耦合理论框架，对冻土路基的温度场和变形场进行数值模拟。Zhang等^［25］给出冻土弹塑性本构模型，冻结阶段和融化阶段的强度、刚度敏感性，基于状态函数概念描述塑性变形，模型描述了荷载和温度变化作用下变形特性和强度演化。Neaupane等^［26］基于BBM模型框架建立了饱和冻土弹塑性模型。但是，这些复杂而精细的耦合理论模型描述性的方程均具有强非线性的特征，求解需要多次复杂迭代计算，数值实现困难，难于直接用于工程设计与施工实践。

鉴于此，本文直接以服务工程应用与设计的目的，充分考虑土体在周期性冻结作用下的冻胀特性，开展开放系统下有侧限条件固结土样（k₀固结）的单向冻胀试验，基于k₀状态下土冻胀力学模型，建立水平场地土体冻胀的计算方法，方便地计算得到不同上覆荷载和位移边界条件给定的条件下地表冻胀的隆起位移与土体应力、应变的规律。

1 冻胀引起水平场地地表隆起变形简化分析模型

1.1　基本假定

土体视为各向同性的非线性弹性材料，并假定与水平基岩或坚实土层接触。如此，土体只发生竖向冻胀，而且冻胀只与竖向深度z有关，与水平所处的位置无关^［27］。因此，水平场地冻胀的预测分析可以简化成一维问题，只需选取高H的单位面积的土柱开展研究。冻结过程中，一维土柱通常分为四种情况：（1）土柱的上、下边界完全被约束，其中，上边界为地表，土性自上而下均匀且温度相同；（2）土柱自上而下分为冻结区和未冻区两个区域，冻结区内土性均匀、但温度分布不均匀，而且土柱上边界自由、下边界固定；（3）土柱的上、下边界均被固定，土性和温度自上而下均不均匀；（4）土柱自上而下分为冻结区和未冻结区两个区域，冻结区内土性和温度分布均不均匀，而且，土柱的上边界被固定、下边界自由^［28］。考虑第二种情况在实际工程中更常见，据此开展研究。

考虑冻结过程中土体强度逐渐变大，为此仍将土体应力和应变表示为非线性弹性关系。冻结过程中，土体冻胀模量受到温度梯度、冷端温度和上覆荷载的综合影响，据此也可以将其表示为温度梯度、冷端温度和上覆荷载的函数^［27］。对于不同深度的土体由于上覆荷载不同，冻胀模量不同，土体内部之间也存在相互作用力。本文，将根据开放系统下侧限冻胀试验指标^［27］计算冻胀。而且，主要关注不同上覆荷载和位移边界条件下，土柱的应力和应变随着温度的变化规律及土柱表面冻胀力与隆起位移。

在本文提出的分析方法中做了如下假定：（1）未冻区土体为均匀、等向的半无限空间弹性体。冻结过程中土体应力增量与应变增量呈非线性弹性关系；（2）场地土体冻胀时只产生竖向冻胀变形，侧向不发生变形。在水平土层冻胀变形的计算中，应用实验室测定的开放系统下侧限冻胀试验冻胀指标。该方法适用于冻胀变形较大且引起大应力变化的实际工程中。

1.2　一维土柱冻胀模型

冻结过程中，一维土柱可分为非冻结区和冻结区，见图1（a）。冻结区内土柱分为n个等厚的子土层，非冻结区内土柱分为m - n个等厚子土层。冻结区内自上至下的子土层分别为第1个子土层1_L、第2个子土层2_L、 ···、第i个子土层i_L、及第n个子土层n_L。非冻结区内自上至下的各子土层分别为：第（n+1）个子土层（n+1）_L、第（n+2）个子土层（n+2）_L和第m个子土层m_L。自表面向下相邻的两个子土层之间分界面分别为：第1个接触面1_S、第2个接触面2_S、 ···、第i个接触面i_S、第n个接触面n_S和第（m - 1）个接触面（m - 1）_S。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 一维土柱简化力学模型

Fig.1 One - dimensional soil column for predicting frost heave （a） and its mechanical model （b）

冻结过程中，为了使得冻结区内子土层不发生任何竖向冻胀，必须在该子土层上施加足够大的上覆荷载。而且，该上覆荷载应与子土层因冻结产生的冻胀应力的大小相等、方向相反。这里，将冻胀应力看作子土层特有的属性，且只与自身的土性和温度状态相关^［27］。为此，可将子土层的冻胀应力分解为两部分： 1）平衡子土层的上覆荷载； 2）施加在上、下相邻子土层的有效冻胀应力。

子土层既产生有效冻胀应力，又产生了抑制冻胀的抗力。子土层抑制冻胀的作用效应采用等价土弹簧表示，并将其施加在子土层上。按上述，建立的场地土层体系冻胀分析模型如图1（b）所示。

土柱中第i段子土层的弹簧系数k_i 表示为：

k_{i} = \frac{E_{e, i}^{k_{0}}}{Δ h_{i}}

（1）

式中： ∆h_i 为第i个子土层的厚度； $E_{e, i}^{k_{0}}$ 为第i个子土层的冻胀模量，可以表示为第i个子土层的冷端温度、温度梯度和上覆荷载的综合函数，参数取值参考文献［27］。

其中， ∆h_i 为第i个子土层的厚度； $E_{e, i}^{k_{0}}$ 为第i个子土层的冻胀模量， $E_{e, i}^{k_{0}}$ 的计算方法如下：

E_{e, i}^{k_{0}} = 1 / (\frac{E_{e q, i}^{k_{0}} - E_{c, i}^{k_{0}}}{E_{c, i}^{k_{0}} E_{e q, i}^{k_{0}}} \frac{1}{1 - \frac{σ_{v}}{σ_{e, i}^{k_{0}}}} + \frac{1}{E_{e q, i}^{k_{0}}})

（2）

其中， $E_{e q, i}^{k_{0}} = \frac{X_{2} \sqrt[]{g r a d T_{i}}}{X_{3}} σ_{v} + \frac{{X_{2}}^{2} \sqrt[]{g r a d T_{i}}}{X_{3}}$ ； $E_{c, i}^{k_{0}} = \frac{\frac{- X_{3}}{X_{1} \sqrt[]{g r a d T_{i}} + X_{4} g r a d T_{i}} - X_{2}}{X_{1} + \frac{X_{3}}{X_{2} \sqrt[]{g r a d T_{i}}} + X_{4} \sqrt[]{g r a d T_{i}}}$ ；

$σ_{e, i}^{k_{0}} = \frac{- X_{3}}{X_{1} \sqrt[]{g r a d T_{i}} + X_{4} g r a d T_{i}} - X_{2}$ ； X₁=-0.0006T_c，_i - 0.13； X₂=0.4T_c，_i +105.14； X₃= - 0.076T_c，_i +9.5； X₄= - 0.0007T_c，_i +0.05， T_c，_i 、 gradT_i 和σ_v分别为子土层i的冷端温度、温度梯度和上覆荷载^［27］。

2 表面冻胀量给定时一维土柱应力和应变计算方法

2.1　平衡方程

设土柱表面冻胀量为给定值a。第1个接触面1_S、第2个接触面2_S、 ···、第i个接触面i_S、第n个接触面n_S及第（m - 1）个接触面（m - 1）_S的竖向变形（向下为正值）分别为v₁、 ···、 v_m_{- 1}，并将各接触面上的竖向位移表示为向量形式：

v = {\{v_{1}, v_{2}, \dots, v_{m - 1}\}}^{T}

（3）

根据k₀状态下饱和粉质黏土的冻胀力学模型^［27］，据此建立子土层的接触面上力的平衡方程。冻胀过程中，一维土柱分析模型中，第i个子土层i_L对第i个分界面i_S作用力分解为：有效冻胀力 $σ_{e, i}^{'}$ （相邻土层之间发生的冻胀力 $σ_{e, i}^{'} = σ_{e, i}^{k_{0}} - σ_{v, i}$ ，方向向下）和方向向上的弹簧力 $σ_{k_{i}} = k_{i} (v_{i} - v_{i - 1})$ 。第i个接触面i_S下面相邻的第（i+1）个子土层（i+1）_L对第i个接触面i_S作用向上的有效冻胀力 $σ_{e, i + 1}^{'}$ ，其中 $σ_{e, i + 1}^{'} = σ_{e, i + 1}^{k_{0}} - σ_{v, i + 1}$ 和方向向下的弹簧力 $σ_{k_{i + 1}}$ ，其中 $σ_{k_{i + 1}} = k_{i + 1} (v_{i + 1} - v_{i})$ 。同时，各子土层还存在附加的上覆初始应力 $σ_{v, i} = γ h_{i} + σ_{o, r}$ ， $γ$ 为土的重度， $h_{i}$ 为第i子土层中心距地表的深度； $σ_{o, r}$ 为作用正冻土体表面的上覆荷载。

2.2　土柱应变的求解方程

（1）第1个接触面1_S处力的平衡方程。

试验中，土样下边界固定，即位移为0，上边界（表面）冻胀量等于定值 $a$ 。据此，根据第1个接触面1_S处力的平衡条件，得到：

σ_{k_{1}} - σ_{k_{2}} + σ_{e, 1}^{'} - σ_{e, 2}^{'} = 0

（4）

式中： $σ_{k_{1}} = k_{1} (v_{1} - a)$ ， $σ_{k_{2}} = k_{2} (v_{2} - v_{1})$ 。

进一步整理得到：

- k_{1} a + (k_{1} + k_{2}) v_{1} - k_{2} v_{2} = σ_{e, 1}^{'} - σ_{e, 2}^{'}

（5）

（2）第i个接触面i_S处力的平衡方程。

第i_S（i<n）接触面受力如图2所示，根据第i个接触面i_S处上力的平衡方程（i<n）：

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 第i_S接触面受力示意图

Fig.2 Force diagram at the contact surface of i_S

σ_{k_{i}} - σ_{k_{i + 1}} + σ_{e, i}^{'} - σ_{e, i + 1}^{'} = 0

（6）

整理得到：

- k_{i} v_{i - 1} + (k_{i} + k_{i + 1}) v_{i} - k_{i + 1} v_{i + 1} = σ_{e, i}^{'} - σ_{e, i + 1}^{'}

（7）

（3）第n个接触面n_S处力的平衡方程。

第n个接触面n_S以上为冻结区，第n个接触面n_S以下为非冻结区，可知：

σ_{e, n + 1}^{'} = 0

同理，根据第n个接触面n_S处力的平衡条件，得到：

- σ_{k_{n}} + σ_{e, n}^{'} - σ_{k_{n + 1}} = 0

（8）

整理得到：

- k_{n} v_{n - 1} + (k_{n} - k_{n + 1}) v_{n} + k_{n + 1} v_{n + 1} = σ_{e, n}^{'}

（9）

（4）第（n+1）个接触面（n+1）_S至第（n+2）个接触面（n+2）_S的力的平衡方程。

第n个接触面n_S以下由第m - n个非冻结的子土层构成。从第（n+1）个接触面（n+1）_S至第（m - 2）个接触面（m - 2）_S，冻胀力 $σ_{e, i}^{'}$ 均零。第j个接触面处j_S力的平衡方程（n<j<m - 1）：

σ_{k_{j}} - σ_{k_{j + 1}} = 0

（10）

进一步整理得到：

- k_{j} v_{j - 1} + (k_{j} + k_{j + 1}) v_{j} - k_{j + 1} v_{j + 1} = 0

（11）

（5）第（m - 1）个接触面（m - 1）_S处力的平衡方程。

由于 $v_{m} = 0$ 和 $σ_{k_{m}} = - k_{m} v_{m - 1}$ ，得到第（m - 1）个接触面（m - 1）_S处力的平衡方程：

σ_{k_{m - 1}} - σ_{k_{m}} = 0

（12）

进而得到：

- k_{m - 1} v_{m - 2} + (k_{m - 1} + k_{m}) v_{m - 1} = 0

（13）

（6）矩阵方程与求解。

如前所述，建立第1个接触面1_S至第（m - 1）个接触面上（m - 1）_S的m - 1个力的平衡方程，其中v₁， v₂， …， v_m 个未知量。进而， m个接触面处力的平衡方程采用矩阵方式表示：

K υ = R

（14）

式中： $K$ 为系数矩阵； $R$ 为荷载向量。

K = [\begin{matrix} k_{1} & - k_{1} \\ - k_{1} & k_{1} + k_{2} & - k_{2} \\ ⋱ & ⋱ & ⋱ \\ - k_{i} & k_{i} + k_{i + 1} & - k_{i} \\ ⋱ & ⋱ & ⋱ \\ k_{m - 1} & k_{m - 1} + k_{m} \end{matrix}]

（15）

R = \{\begin{matrix} σ_{e, 1}^{'} - σ_{e, 2}^{'} + k_{1} a \\ σ_{e, 2}^{'} - σ_{e, 3}^{'} \\ ⋮ \\ σ_{e, i}^{'} - σ_{e, i + 1}^{'} \\ ⋮ \\ σ_{e, n - 1}^{'} - σ_{e, n}^{'} \\ σ_{e, n}^{'} \\ 0 \\ ⋮ \\ 0 \\ ⋮ \\ 0 \\ 0 \end{matrix}\}

（16）

可见， $K$ 为三对角矩阵，分别令：

\begin{array}{l} K = [\begin{matrix} b_{1} & c_{1} \\ a_{2} & b_{2} & c_{2} \\ ⋱ & ⋱ & ⋱ \\ a_{m - 2} & b_{m - 2} & c_{m - 2} \\ a_{m - 1} & b_{m - 1} \end{matrix}] \\ R = {\{d_{1}, d_{2}, \dots, d_{m - 1}\}}^{T} \end{array}

（17）

并将其分解为：

K = L H

（18）

其中，

L = [\begin{matrix} 1 \\ l_{2} & 1 \\ ⋱ & ⋱ \\ l_{m - 1} & 1 \end{matrix}]

（19）

H = [\begin{matrix} h_{1} & p_{1} \\ h_{2} & ⋱ \\ ⋱ & p_{m - 1} \\ h_{m - 1} \end{matrix}]

（20）

计算 $L$ 和 $H$ 中其余元素的表达式为：

\{\begin{array}{l} p_{i} = c_{i} (i = 1, 2, 3, \dots, m - 1) \\ h_{1} = b_{1} \\ l_{i} = \frac{a {}_{i}}{h_{i - 1}} (i = 2, 3, \dots, m - 1) \\ h_{i} = b_{i} - l_{i} c_{i - 1} (i = 2, 3, \dots, m - 1) \end{array}

（21）

将式（18）代入式（14）得到：

L H υ = R

（22）

令：

y = H υ

（23）

式中， $y = {\{y_{1}, y_{2}, \dots, y_{m - 1}\}}^{T}$ 。

公式（23）改为：

\{\begin{matrix} y_{1} \\ y_{2} \\ ⋮ \\ y_{m - 1} \end{matrix}\} = [\begin{matrix} h_{1} & c_{1} \\ h_{2} & ⋱ \\ ⋱ & c_{m - 1} \\ h_{m - 1} \end{matrix}] \{\begin{array}{l} v_{1} \\ v_{1} \\ ⋮ \\ v_{m - 1} \end{array}\}

（24）

得到：

L y = R

（25）

即：

[\begin{matrix} 1 \\ l_{2} & 1 \\ ⋱ & ⋱ \\ l_{m - 1} & 1 \end{matrix}] \{\begin{matrix} y_{1} \\ y_{2} \\ ⋮ \\ y_{m - 1} \end{matrix}\} = \{\begin{matrix} d_{1} \\ d_{2} \\ ⋮ \\ d_{m - 1} \end{matrix}\}

（26）

计算 $y$ 的元素公式表示：

\{\begin{array}{l} y_{1} = d_{1} \\ y_{i} = d_{i} - l_{i} y_{i - 1} (i = 2, 3, \dots, m - 1) \end{array}

（27）

最后，采用追赶法分别求解，得到h₁， y₁， l₂， h₂， y₂， …， l_m_{- 1}， h_m_{- 1}和y_m_{- 1}。

$R$ 的元素计算公式为：

\{\begin{array}{l} v_{m - 1} = \frac{y_{m - 1}}{h_{m - 1}}, \\ v_{i} = \frac{y_{i} - c_{i} v_{i + 1}}{h_{i}} (i = m - 1, m - 2, \dots, 1) \end{array}

（28）

2.3　土柱应力的求解方程

2.3.1　土柱表面冻胀力 ${\bar{σ}}_{0}$

根据地表处力的平衡条件，得到：

{\bar{σ}}_{0} = σ_{e, 1}^{'} - σ_{k_{1}}

（29）

式中： $σ_{e, 1}^{'}$ 为第1个子土层1_L的冻胀力； $σ_{k_{1}}$ 为第1个子土层1_L的抵抗力。

试验测得地表处冻胀量等于 $a$ ，得到：

σ_{k_{1}} = k_{1} (v_{1} - a)

（30）

地表处 $σ_{v, 1} = σ_{0, r}$ ，得到：

σ_{e, 1}^{'} = σ_{e, 1} - σ_{o, r}

（31）

将式（30）代到式（29）进一步得到：

{\bar{σ}}_{o} = σ_{e, 1} - σ_{o, r} - k_{1} (v_{1} - a)

（32）

2.3.2　子土层间相互作用力

作用在任意i个子土层分界面上i_L的冻胀力 ${\bar{σ}}_{i}$ 计算公式：

{\bar{σ}}_{i} = σ_{e, i + 1}^{'} - σ_{k_{i + 1}}

（33）

式中： $σ_{e, i + 1}^{'}$ 为第（i+1）个子土层（i+1）_L的冻胀力； $σ_{k_{i + 1}}$ 为第（i+1）个子土层（i+1）_L的抵抗力。

3 表面冻胀力给定时一维土柱应变和应力计算方法

3.1　平衡方程

地表冻胀力已知时，地表冻胀量v₀为未知量，拟求的冻胀变形向量与第2节中冻胀变形向量多了一个未知量v₀。同样，求解的方程组中需要补充一个方程式。由冻胀力、冻结抗力和上覆荷载共同作用下地表（第0个接触面0_S）处力的平衡条件，建立额外一个方程式。

土柱下边界固定，此处位移为0，并假设初始上覆荷载为p，根据第0个接触面0_S处力的平衡条件，得到：

σ_{e, 1}^{'} - σ_{k_{1}} = 0

（34）

进一步得到：

k_{1} v_{0} - k_{1} v_{1} = p - σ_{e, 1}

（35）

然后，逐一得到第1个接触面1_S、第i个接触面i_S、第n个接触面n_S、第（n+1）个接触面（n+1）_S至第（m - 2）个接触面（m - 2）_S、第（m - 1）个接触面（m - 1）_S处力的平衡方程。据此，可以将第1个接触面1_S至第（m - 1）个接触面（m - 1）_S对应的m - 1个平衡方程采用矩阵形式见式（14）。

3.2　应力的求解

同时，给出作用在任意第i个子土层对应的接触面i_L上的冻胀力 ${\bar{σ}}_{i}$ 的计算公式：

{\bar{σ}}_{i} = σ_{e, i + 1}^{'} + σ_{k_{i + 1}}

（36）

式中： $σ_{e, i + 1}^{'}$ 为第（i+1）个子土层（i+1）_L的冻胀力； $σ_{k_{i + 1}}$ 为第（i+1）个子土层（i+1）_L的抗力。

4 水平场地冻胀计算方法验证

为了验证水平场地冻胀的计算方法，实施了开放系统下饱和粉质黏土的冻胀试验，并得到土样试验测得的冻胀时程。

4.1　冻胀试验

仪器采用自主研发的土体单向冻胀试验系统，包括制冷系统、温控系统、补水系统和数据采集系统等，见图3。试验用土为哈尔滨粉质黏土，属于冻胀敏感性土。土样为高15.0 cm，直径10 cm的圆柱体。基于土体冻胀基本特征，考虑上覆荷载对土体冻胀影响，开展开放系统条件下（马里奥特瓶无压补水）逐级降温的单向冻胀试验。试验前，土样暖端温度设置为+2 ℃，冷端温度设置为0 ℃并进行预冷。待温度稳定后，将暖端温度调至恒定的+2 ℃，土试样顶部的冷端温度从初始状态0 ℃逐级降至 - 5、 - 10、 - 15 ℃，每一阶段冻结时间均为72 h。对土试样分别施加0、 0.028、 0.054 MPa的上覆荷载，详见文献［27］。

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 试验系统结构图

Fig.3 Test system structure drawing

试样内布置6支MF5E-2.252F型针式温度传感器，阻值互换精度±0.2%，分别埋设在距试件表面为0、 30、 60、 90、 120 mm和150 mm处。土试样各深度处的温度时程见图4，发现土试样中各位置的温度不同，随着冻结持续，各处温度趋势保持一致，但是达到稳定的时间不相同。这更多归因于土样冷端温度出现变化后，传递过程中热量逐渐被吸收耗掉，热量传递的时效性和过程性决定了下层土体温度改变具有滞后性且速度较慢。

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 冻结过程中土体温度时程

Fig.4 Soil temperature changing with time at various depths

各级荷载作用下试样表面冻胀时程，见图5。当冷端温度继续降低并进入下一级冻结阶段之后，冻胀继续增大，但增加幅度变小。不同上覆荷载条件下，冻胀规律与趋势基本相似。冻结后，上覆荷载分别为0、 0.028、 0.054 MPa时，土样相应冻胀分别为43.39、 30.43、 23.46 mm。

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 土试样表面冻胀时程

Fig.5 Frost heave of soil surface changing with time under various loads

4.2　模型验证

选取冻结过程中24、 48、 72、 96、 120、 144、 168、 192 h和216 h共9个时刻的冻胀量作为验证点，将试验测得的冻胀量与冻胀变形预测解析法计算得到的冻胀量进行对比。首先，采用试验过程中温度场结果（图4），计算冻结过程中冻深的变化曲线，并应用第3节的方法，计算地表冻胀力给定时土中应变和应力；然后，与图5中试样冻胀时程进行对比。

图6为试样表面冻胀时程试验值与模拟值的对比图。可见，试验结果和预测结果吻合很好，并随着冻结持续，试验和计算得到的土样表面冻胀呈现相似的规律，只是计算值偏大。这应该由于计算值为假定冻结时间足够长时，冻结达到稳定状态对应的冻胀量。然而，试验中每一级冻结阶段只持续72 h，土柱表面冻胀还未达到稳定状态。综上，新建立的水平场地土体冻胀计算方法［式（36）］可以较好地预测地表冻胀。

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 实测和计算的土试样表面冻胀量时程曲线

Fig.6 Surface frost heave of specimens changing with time under various loads， calculated and experimented

5 冻结过程中水平场地地表抬高量计算及场地评估

5.1　算例

假定场地为水平成层的土层，土体水平均匀且范围较大，即可认为冻胀土层处于k₀状态。冻土层厚度的相关数据，根据所在场地的现场监测获得。假定地表温度为 - 5 ℃时，冻土层厚度为333 mm。地表温度为 - 10 ℃时，冻土层厚度666 mm。地表温度为 - 15 ℃时，冻土层厚度为999 mm。地表处作用上覆荷载范围为0 ~ 0.3 MPa。地表温度从0 ℃逐级降至 - 5、 - 10、 - 15 ℃。本研究假定温度场稳定时，土体内部温度沿着深度呈线性分布。土的重度γ=1.75×10^-5 N·mm^-3。设置的上覆荷载应力分别为0、 0.075、 0.15、 0.225 MPa和0.3 MPa，以考察上覆荷载的影响效应。计算解析解得到不同地表温度、上覆荷载作用下场地表面冻胀抬升量变化规律，如图7所示。

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 地表的冻胀与温度关系

Fig.7 The relationship between frost heave and surface temperature under various overlying loads

预测场地表面冻胀量与表面温度之间的关系，地表温度在 - 10 ℃以上时，场地表面冻胀抬升较为明显，最快抬升速率为40 mm·℃^-1，低于该温度时各个工况下的冻胀抬升速率明显降低，其中，在0.3 MPa上覆荷载应力下的抬升速率接近0 mm·℃^-1。

图8给出场地表面冻胀量与初始上覆荷载之间变化关系，可以看出，场地表面冻胀量随着冷端温度的降低而增大，并随着上覆荷载增大而逐渐减小。

图8

新窗口打开| 下载原图ZIP| 生成PPT

图8 地表冻胀与初始上覆荷载关系

Fig.8 The relationship between frost heave and initial overburden pressure at the soil surface under various surface temperatures

6 结论

本文基于土体冻胀非线性弹性力学模型，推导了开放系统下半空间无限体中土单元在自重下冻胀的数学表达式，提出了考虑上部压载效应与位移边界变化的场地表面隆起位移及土体应力与应变的计算方法，主要结论如下：

（1）采用建立的k₀状态下的冻胀力学模型，推导了自重下半空间无限体的土体冻胀的数学表达式，建立了水平场地土体冻胀的计算方法，给出了考虑上部压载和位移边界变化的地表隆起位移与土体应力、应变的解析方法，基于开放系统下粉质黏土逐级降温冻胀试验数据对解析解计算结果进行了验证，结果表明该方法具有较强的可靠性。

（2）计算解析解与试验所得结果一致性较高。具体而言，随冻结持续，试验和计算得到的土样表面冻胀量呈现相似的规律。由于试验过程冻结时间持续72小时，与计算假设存在差异，导致计算值整体大于试验值。

（3）预测场地表面冻胀量与表面温度之间的关系，地表温度在 - 10 ℃以上时，场地表面冻胀抬升较为明显，最快抬升速率为40 mm·℃^-1，低于该温度时各个工况下的冻胀抬升速率明显降低，其中，在0.3 MPa上覆荷载应力下的抬升速率接近0 mm·℃^-1。场地表面冻胀随着冷端温度的降低而增大，并随着上覆荷载增大而逐渐减小。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Wang

Xiaogang

， Xing

Yichuan

， Zhao

Jianming

， et al.

Geotechnical engineering problems associated with construction of water resources projects in western China

［J］. Chinese Journal of Geotechnical Engineering， 2007， 29（8）： 1129 - 1134.