球模仪测试冻土松弛模量的非线性Kelvin解答及其试验研究

doi:10.7522/j.issn.1000-0240.2020.0018

球模仪测试冻土松弛模量的非线性Kelvin解答及其试验研究

黄灿杰^,¹^,², 张泽^,²^,³, 金会军¹^,²^,³, 冯文杰², 金豆豆²^,⁴, MELNIKOV Andrey⁵

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

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

3.东北林业大学土木工程学院/寒区工程与科学技术研究院，黑龙江哈尔滨 150040

4.中国科学院大学，北京 100049

5.俄罗斯科学院西伯利亚分院梅尔尼科夫冻土研究所，俄罗斯雅库茨克 117997

Nonlinear Kelvin model solution of frozen soil relaxation modulus and experimental study based on spherical template indenter

HUANG Canjie^,¹^,², ZHANG Ze^,²^,³, JIN Huijun¹^,²^,³, FENG Wenjie², JIN Doudou²^,⁴, MELNIKOV Andrey⁵

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

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

3.School of Civil Engineering / Institute of Cold Regions Engineering，Science and Technology，Northeast Forestry University，Harbin 150040，China

4.University of Chinese Academy of Sciences，Beijing 100049，China

5.Melnikov Permafrost Institute，Siberian Branch，Russian Academy of Sciences，Yakutsk 117997，Russia

通讯作者: 张泽，副研究员，从事冻土工程与环境研究. E-mail： zhangze@lzb.ac.cn.

编委: 武俊杰

收稿日期: 2019-12-02 修回日期: 2020-04-22

基金资助:

国家自然科学基金项目.  41771078.  41971085
中国科学院国际合作局对外合作重点项目.  131B62KYSB20170012
国家自然科学基金中俄交流项目.  4191101321
Russian基金项目:Foundation for Basic Research.  20-55-53006

Received: 2019-12-02 Revised: 2020-04-22

作者简介 About authors

黄灿杰（1992-），男，广东汕尾人，2017年在广东工业大学获学士学位，现为哈尔滨工业大学在读硕士研究生，从事冻土物理力学研究.E-mail：15014252819@163.com , E-mail：15014252819@163.com

摘要

冻土的松弛模量是进行蠕变力学分析的关键参数。但在理论方面，尚未基于球模仪试验建立较完善的松弛模量表征公式。基于半空间黏弹性理论和分数阶微积分理论，获得了松弛模量公式的非线性Kelvin模型解答并进行了材料参数影响分析。在给出恒载条件下，可以用单轴蠕变试验的分析结果评估基于球形压痕试验的冻土松弛模量来预测准确性。随后分别以冻结细砂的球模仪试验曲线、冻结砂和冻结砂质黏土的单轴蠕变曲线为例，预测了各条试验曲线和冻土的松弛模量，并通过其他试验实测数据进行了检验。结果表明：该松弛模量表征公式的材料参数最少，随各材料参数变化具有较佳单调性；基于球形模板试验，非线性Kelvin解答能对冻结细砂的松弛模量和位移曲线做出合理可靠的预测。非线性Kelvin模型能很好地吻合单轴蠕变试验曲线，基于单轴蠕变试验的冻结砂松弛模量随时间变化规律和数值范围与球模仪试验分析结果较一致，冻结砂质黏土的松弛模量随应力和温度的变化规律再现了已有研究结论。

关键词： 球模仪 ; 分数阶理论 ; 非线性Kelvin解答 ; 参数分析 ; 冻土松弛模量

Abstract

Relaxation modulus of frozen soil is a most important factor in creep analysis. However， there is still no a more perfected formulation of relaxation modulus according to the spherical template test theoretically. Based on half spacious visco-elastic theory and fractional order thesis， this paper had found formulations of relaxation modulus and analyzed the impact of parameters included in formulations according to the nonlinear Kelvin model. We propose a universal method for acquiring relaxation modulus based on dead load spherical template indenter tests and uniaxial creep test. This means relaxation modulus analysis results based on uniaxial creep tests can be applied for evaluating the results based on spherical template indenter tests. Then， the model is used for predicting the spherical template indenter tests curves of frozen fine sand and uniaxial creep curves of frozen sand and frozen sandy clay， as well as relaxation modulus of them， respectively. Besides， other experimental data are provided to evaluate the calculation results. Analytical results show that relaxation modulus characterization function has the least material parameters and good monotony versus material parameters. Nonlinear Kelvin model can predict the reasonable and reliable value of relaxation modulus of frozen fine sand on the basis of spherical template indenter test. The nonlinear Kelvin model has a good description of uniaxial creep curves of frozen sand and frozen sandy clay. The changing regularity of relaxation modulus of frozen sand based on uniaxial creep tests agrees with that based on spherical template indenter test. The changing regularity of frozen sandy clay relaxation modulus with temperature and stress is consistent with existed conclusions.

Keywords： spherical template indenter ; fractional order thesis ; solution of nonlinear Kelvin model ; parameters analysis ; relaxation modulus of frozen soil

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

本文引用格式

黄灿杰, 张泽, 金会军, 冯文杰, 金豆豆, MELNIKOV Andrey. 球模仪测试冻土松弛模量的非线性Kelvin解答及其试验研究[J]. 冰川冻土, 2020, 42(2): 550-561 doi:10.7522/j.issn.1000-0240.2020.0018

HUANG Canjie, ZHANG Ze, JIN Huijun, FENG Wenjie, JIN Doudou, MELNIKOV Andrey. Nonlinear Kelvin model solution of frozen soil relaxation modulus and experimental study based on spherical template indenter[J]. Journal of Glaciology and Geocryology, 2020, 42(2): 550-561 doi:10.7522/j.issn.1000-0240.2020.0018

0 引言

冻土的长期强度和松弛模量为寒区基础工程的建设提供了基本的冻土力学依据，可以通过多种试验方式获取。相对于其他试验仪器，球模仪（球形模板压入仪）具有试验周期短、试验操作简便、试验效果达到工程精度要求、场地适应性强等多个方面的优势^［1］，可确定冻土的长期强度，已经被俄罗斯国家及其建筑行业标准列为最基本的土工试验方法。

将球形压痕试验作为冻土力学性质测试的方法最早由崔托维奇提出。崔托维奇和维亚洛夫随后基于理想黏塑性非强化体假设进行严密推导，建立了通过测试获得冻土长期强度的理论基础^［1］，所测定的长期强度在大量工程实践中被证明是较为可靠的。在崔托维奇和维亚洛夫工作的基础上，冻土力学领域的研究者较为关注冻土的长期强度随冻融循环次数、温度、含水率、土质、土的化学组成等影响因素的变化规律，以此作为冻土长期强度预测的依据^［2-5］，对冻土松弛模量和蠕变柔量的研究还较少。主要原因是在球形模板和冻土样的接触变形问题上，被证明得不到完全闭合的弹塑性解析解^［6］，导致从测试数据当中反分析出试样准确的变形特性参数（包括松弛模量、蠕变柔量、泊松比）较难。因此，学者们基于各种假设提出简化分析方法，以确定冻土的变形特性参数^［7］。

早在1882年，海因里希·鲁道夫·赫兹（Heinrich Rudolf Hertz，电磁波发现者之一）采用半逆解法获得弹性球和半空间弹性体接触问题的解答，在此基础上建立了Hertz接触理论，它在冶金、机械和化工中已得到广泛的运用^［8］。然而Hertz接触理论比较适用于金属等表面比较光滑，且在中小等级荷载下塑性区扩展较小的接触问题，在这些问题中材料主要处于弹性状态。为了将Hertz接触理论的运用范围推广至黏弹塑性材料，分析时采用半经验半理论的方法，基于实际问题引入了简化假设。如：将球形模板和土样均视为弹性体，基于Hertz理论的严格解答， Roman等^［7］认为土样始终与球形模板呈小面积区域紧密接触状态，修正得到Roman解答，其结果和三轴压缩试验的较接近。但是，所提供的松弛模量公式不是时间的显式，在预测松弛模量随时间的变化规律时需要借助经验公式进行拟合。

模拟研究表明^［9-10］，在单级恒载条件下，随着球形压痕压入深度的扩大和渐趋稳定，应力的大小在减少，应力分布范围在扩大并渐趋稳定，这是球形压痕试验不同于其他强度试验的特点之一。它决定从冻土的总位移-时间曲线难以区分出弹性阶段、黏塑性阶段的确切分界点，更不能像其他强度试验一样直接获得应力-应变关系曲线，导致反向分析方法是从球形压痕试验数据提取松弛模量等力学参数的最佳分析方法。Lee等^［11］最早采用基于线性黏弹性模型和球形压痕试验的反向分析法对聚合物的蠕变柔量进行研究，他们的工作对后来的许多研究者有启发作用。Huang等^［12］基于球形压痕试验和线性黏弹性模型，提出蠕变柔量和松弛模量关于时间显式的计算方法，该方法实际上基于Hertz解答和Boltzmann线性叠加原理，他们的研究思想为Martynova^［13］所继承。后者在文献［12］基础上，借助多项式展开提出更简化的蠕变柔量和松弛模量计算方法，分析结果和单轴蠕变试验的较接近。不足之处是所提出的多项式是为分析简便而经验性地引入，待定参数较多，需要进行多次参数拟合工作。实际上，分析过程的繁复和运用线性黏弹性模型有关。这类模型可以统一用Prony级数表示，在模型简单时，其提供的描述能力实在有限。为能和试验数据吻合，不得不增加模型的复杂性，导致推导出的Prony级数的各项系数极为复杂，确定各项系数的工作较繁琐。相反，分数阶黏弹性模型的模型简洁，物理概念清晰，在基于单轴和三轴蠕变试验建立冻土蠕变模型方面已经有较成功的研究结果^［14-15］。

基于以上分析，本文运用分数阶Kelvin黏弹性模型，使用最少材料参数建立蠕变柔量和松弛模量以时间为变量的显式多项式。所提出多项式的各项系数仅需由一条单轴蠕变或球形压痕蠕变试验曲线拟合便能全部确定，减少分析计算过程，便于工程运用。为验证分析结果的准确性，还与单轴蠕变试验的分析结果进行对比，以保证结果可靠。

1 球模仪试验原理

球模仪的结构和试验步骤，在文献［1］已经有详细介绍，本文不再重复。其测定长期强度的原理建立在理想黏性非强化体解答的基础之上。

C = K \frac{F}{π D S}

（1）

式中： C为随时间变化的等效黏聚力（MPa）； F为作用在球形模板的轴向荷载（N）； K为比例系数，取0.18； D为球形模板直径（mm），一般在8 ~ 40 mm，俄罗斯国家标准ΓOCT 12248-96^［16］建议取球模直径22 mm，在后续计算中采用此值； S为球形模板压入位移（mm）。球模压入土样简化模型见图1。

图1

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

图1 球模压入土样简化模型示意图^［10］

Fig.1 Skeleton figure showing spherical template indenter pressing ball into the soil^［10］

利用球形模板试验测试非黏性土时，所测的强度值包含了内摩擦角和黏聚力的影响，称为等效黏聚力C_eq（图2）^［1］。

图2

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

图2 等效黏聚力强度随时间变化示意图

Fig.2 Equivalent cohesive strength C_eqchanging with time

一般把在0 ~ 25 s内测定的等效黏聚力作为瞬时强度C₀（图2），后续在测定松弛模量瞬时值时采用与等效黏聚力瞬时值一致的时间范围。当时间趋于较大值时（如50年），所对应的冻土不发生蠕变破坏的强度为长期强度极限C_∞。

在式（1）中，等效黏聚力与材料的硬度概念类似，是土体力学性能的综合性表征物理量，且不是时间的显式公式。若假设土体服从Mohr-Coulomb准则，需要设法由等效黏聚力分离出黏聚力和内摩擦角，以便于工程设计和理论研究。这是从强度角度分析的研究思路。本文从接触变形出发，根据Hertz弹性解答及遗传蠕变理论导出土体的蠕变柔量、松弛模量和时间的显式关系。据此，实现球形压痕试验和单轴压缩试验两种试验分析结果的相互对比和验证。

2 黏弹性解答

Roman等^［7］根据弹性半空间的Hertz解答进行修正，获得了基于球形模板的松弛模量公式。

E = \frac{3 (1 - ν^{2}) F}{4 S^{1.5} \sqrt[]{D - S}} \approx \frac{3 (1 - ν^{2}) F}{4 S^{1.5} \sqrt[]{D}}

（2）

式中： S为球形模板轴向位移（mm）； D为球形模板直径（mm）； F为外荷载（N）； ν为冻土试样的泊松比； E为冻土试样的松弛模量（MPa）。

Roman等^［7］结合球模试验数据对比该式对含盐冻结黏土变形模量的计算结果和常规三轴试验的结果，证实两者较为接近。但是，式（2）是基于球形模板和土样之间无摩擦接触假设的弹性理论解答，理论上与时间过程无关，因而没法说明蠕变过程中冻土松弛模量随时间衰减的原因。球模仪在加载过程中除瞬时点状接触作用容易产生塑性变形外，随着接触面积扩大，应力值逐渐减小，不再继续保持屈服状态^［17］，在加载后期主要发生黏弹性变形。为便于阐释本文的分析思路，以下用最简单的非线性Kelvin模型为例进行推导（图3）。该模型参数最少，适合描述蠕变过程。

图3

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

图3 非线性Kelvin模型示意图［图中γ为工程剪应变，其余符号意义见式（3）、式（4）注解］

Fig.3 Skeleton figure of nonlinear Kelvin model： where γ stands for engineering shear strain， the meanings of other characters within the figures are illustrated in Eqs. （3） and （4）

采用三维非线性Kelvin模型［图3（b）］，该模型中假定偏应变与偏应力之间为非线性黏弹性，体积应变与体积应力之间为线性弹性，对应的本构关系为

S_{i j} (t) = 2 (G_{0} + η \frac{\partial r}{\partial t^{r}}) e_{i j} (t) = 2 G_{0} e_{i j} (t) + 2 \frac{d}{d t} [\int_{0}^{t} (\frac{1}{Γ (1 - r)} \frac{η}{{(t - τ)}^{r}}) e_{i j} (τ) d τ]

（3）

σ_{i i} (t) = 3 K ε_{i i} (t)

（4）

式中： S_ij 、 e_ij 分别为偏应力张量及偏应变张量； σ_ii 为第一应力不变量； ε_ii 为体积应变； G₀为瞬时剪切模量，不随时间改变； K为体积模量， 3K=2G₀（1+ν₀）/（1-2ν₀），其中ν₀为泊松比初始值； η为Abel黏滞系数； r为分数阶， 0≤r≤1。

式（3）的第二个等式采用Riemann-Liouvile给出的形式^［18］。基于式（2），由黏弹性模型对应定律有

\tilde{F} = \frac{8 \tilde{G}}{3 (1 - \tilde{ν})} D^{0.5} {\tilde{S}}^{1.5}

（5）

式中： $\tilde{F}$ 、 $\tilde{S}$ 、 $\tilde{G}$ 和 $\tilde{ν}$ 分别为经Laplace变换后的轴向荷载、轴向压入位移、剪切模量和泊松比。

将式（3）和式（4）进行Laplace变换得到

2 \tilde{G} = {\tilde{S}}_{i j} / {\tilde{e}}_{i j} = 2 (G_{0} + η δ^{r})

（6）

\tilde{ν} = \frac{3 K - 2 \tilde{G}}{6 K + 2 \tilde{G}} = \frac{3 K - 2 (G_{0} + η δ^{r})}{6 K + 2 (G_{0} + η δ^{r})}

（7）

式中： δ为Laplace变换的时间因子。

由于轴向荷载保持不变，可用单位阶跃函数表示轴向荷载。

F = F_{c r} \cdot H (t) = {\begin{matrix} F_{c r}, t \geq 0 \\ 0, t < 0 \end{matrix}

（8）

式中： F_cr为竖向恒定荷载； H（t）为关于时间t的Hessian函数。

式（8）经Laplace变换后得到 $\tilde{F} = F_{c r} / δ$ ，并与式（6）和式（7）代入式（5），有

{\tilde{S}}^{1.5} = \frac{3 F_{c r}}{16 D^{0.5} δ} [\frac{1}{G_{0} + η δ^{r}} + \frac{3 (1 - 2 ν_{0})}{3 G_{0} + (1 - 2 ν_{0}) η δ^{r}}]

（9）

对式（9）进行Laplace逆变换得到

S^{1.5} (t) = \frac{3 F_{c r}}{16 D^{0.5} G_{0}} {\sum_{i = 1}^{n} {(- 1)}^{i - 1} \cdot {{(\frac{G_{0}}{η})}^{i} + (1 - 2 ν_{0}) \cdot {[\frac{3 G_{0}}{(1 - 2 ν_{0}) η}]}^{i}} \cdot \frac{t^{r}}{Γ (r \cdot i + 1)}}

（10）

式中： Г（x）为Gamma函数， $Γ (x)$ = $\int_{0}^{+ \infty} a^{x - 1}$ ⋅exp（-a）da。

作为特例，当r→0， n→+∞时，对式（10）取极限便为线性Kelvin模型的分析结果。

\begin{array}{l} S^{1.5} (t) = \frac{3 F_{c r}}{16 D^{0.5} G_{0}} {1 - e x p (- \frac{G_{0} t}{η}) + (1 - 2 ν_{0}) \cdot \\ [1 - e x p (- \frac{3 G_{0} t}{(1 - 2 ν_{0}) η})]} \end{array}

（11）

根据遗传蠕变理论，若压入位移和荷载的初始值为0，压入位移和荷载关系满足

S^{1.5} (t) = \frac{3 (1 - ν_{0}^{2})}{4 D^{0.5}} \int_{0}^{t} J (t - τ) d F (τ)

（12）

F (t) = \frac{4 D^{0.5}}{3 (1 - ν_{0}^{2})} \int_{0}^{t} Y (t - τ) d S^{1.5} (τ)

（13）

式中： J（t）为蠕变柔量； Y（t）为松弛模量。蠕变柔量的物理意义为：对施加恒定应力值的一维变形体，在t时刻，单位应力作用下变形体产生的应变值。松弛模量的物理意义为：对施加恒定应变值的一维变形体，在t时刻，变形体发生单位应变所需要的应力值。

可见，蠕变柔量和松弛模量分别是弹性力学中柔度和刚度概念在时间域上的推广。

将式（8）代入式（12），得到

S^{1.5} (t) = \frac{3 F_{c r} (1 - ν_{0}^{2})}{4 D^{0.5}} \int_{0}^{t} J (t - τ) δ (τ) d τ = \frac{3 F_{c r} (1 - ν_{0}^{2}) J (t)}{4 D^{0.5}}

（14）

式中： δ（t）为Dirac函数，是Hessian函数的广义导数， δ（t）=dH（t）/dt。

对比式（10）与式（14）得到

J (t) = \frac{1}{4 (1 - ν_{0}^{2}) G_{0}} {\sum_{i = 1}^{n} {(- 1)}^{i - 1} \cdot {{(\frac{G_{0}}{η})}^{i} + (1 - 2 ν_{0}) \cdot {[\frac{3 G_{0}}{(1 - 2 ν_{0}) η}]}^{i}} \cdot \frac{t^{r}}{Γ (r \cdot i + 1)}}

（15）

将式（13）进行Laplace变换得到

\tilde{F} = \frac{4 D^{0.5}}{3 (1 - ν_{0}^{2})} \tilde{Y} \cdot {\tilde{S}}^{1.5} \cdot δ

（16）

将式（16）与式（5）对比，有

\tilde{Y} = \frac{2 (1 - ν_{0}^{2})}{δ} \frac{\tilde{G}}{1 - \tilde{ν}}

（17）

将式（6）和式（7）代入式（17），整理得到

\tilde{Y} = \frac{2 (1 - ν_{0}^{2}) (G_{0} + η δ^{r})}{δ} \frac{3 G_{0} + (1 - 2 ν_{0}) η δ^{r}}{3 (1 - ν_{0}) G_{0} + 2 (1 - 2 ν_{0}) η δ^{r}}

\begin{array}{l} = \frac{(1 - ν_{0}^{2})}{δ} [\frac{G_{0} (5 - ν_{0})}{2 (1 - 2 ν_{0})} + η δ^{r} - \\ \frac{3 (1 + ν_{0})^{2} G_{0}^{2}}{2 (1 - 2 ν_{0})} \frac{1}{3 (1 - ν_{0}) G_{0} + 2 (1 - 2 ν_{0}) η δ^{r}}] \end{array}

\begin{array}{l} \approx \frac{(1 - ν_{0}^{2})}{δ} {\frac{G_{0} (5 - ν_{0})}{2 (1 - 2 ν_{0})} + η δ^{r} + \frac{3 (1 + ν_{0})^{2} G_{0}^{}}{2 (1 - 2 ν_{0}) (1 - ν_{0})} \cdot \\ {\sum_{i = 1}^{n} {(- 1)}^{i} \cdot {[\frac{3 (1 - ν_{0}) G_{0}}{2 (1 - 2 ν_{0}) η}]}^{i} δ^{- r \cdot i}}} \end{array}

（18）

对式（18）进行Laplace逆变换，得到

Y (t) = (1 - ν_{0}^{2}) \cdot {\frac{G_{0} (5 - ν_{0})}{2 (1 - 2 ν_{0})} + \frac{η}{t^{r} \cdot Γ (1 - r)} + \frac{3 (1 + ν_{0})^{2} G_{0}^{}}{2 (1 - 2 ν_{0}) (1 - ν_{0})} \cdot {\sum_{i = 1}^{n} {(- 1)}^{i} \cdot {[\frac{3 (1 - ν_{0}) G_{0}}{2 (1 - 2 ν_{0}) η}]}^{i} \frac{t^{r \cdot i}}{Γ (r \cdot i + 1)}}}

（19）

在单轴蠕变试验和单轴松弛试验中，根据遗传蠕变理论可知如下关系^［19］成立。

σ (t) = ε_{0} \cdot Y (t) + \int_{0}^{t} Y (t - τ) d ε (τ)

（20）

ε (t) = σ_{0} \cdot J (t) + \int_{0}^{t} J (t - τ) d σ (τ)

（21）

式中： ε₀为应变初值； σ₀为应力初值； σ（t）为应力； ε（t）为轴向应变。

利用式（10）、式（15）、式（19）、式（20）和式（21）可以实现恒定荷载的球形压痕试验和恒定应力单轴蠕变试验的数据对比，以验证球形压痕试验分析结果的准确性。

3 参数影响

由于瞬时剪切模量G₀、 Abel黏滞系数η、泊松比初始值ν₀、分数阶r与时间t或δ无关，分析参数取值变化对松弛模量Y的影响时可以基于式（18）依次对G₀、 η、 r和ν₀进行求偏导。整理得到

\begin{array}{l} \frac{\partial \tilde{Y} (δ)}{\partial G_{0}} = \frac{(1 - ν_{0}^{2})}{δ} {\frac{(5 - ν_{0})}{2 (1 - 2 ν_{0})} + \frac{(1 + ν_{0})^{2}}{4 (1 - 2 ν_{0}) (1 - ν_{0})} \\ {{[\frac{3 G_{0} (1 - ν_{0})}{3 G_{0} (1 - ν_{0}) + 2 (1 - 2 ν_{0}) η δ^{r}} - 1]}^{2} - 1}} \end{array}

\geq \frac{2 (1 - ν_{0}^{2})}{δ} \frac{(5 - ν_{0}) (1 - ν_{0}) - (1 + ν_{0})^{2}}{4 (1 - 2 ν_{0}) (1 - ν_{0})} = \frac{2 (1 + ν_{0})}{δ} > 0

（22）

\begin{array}{l} \frac{\partial \tilde{Y} (δ)}{\partial η} = \frac{(1 - ν_{0}^{2})}{δ^{1 - r}} {1 + \frac{3 (1 + ν_{0})^{2}}{(1 - 2 ν_{0})^{2}} \\ {[\frac{(1 - 2 ν_{0})}{3 G_{0} (1 - ν_{0}) + 2 (1 - 2 ν_{0}) η δ^{r}}]}^{2}} > 0 \end{array}

（23）

\begin{array}{l} \frac{\partial \tilde{Y} (δ)}{\partial r} = \frac{η \cdot (1 - ν_{0}^{2}) \cdot l n δ}{δ^{1 - r}} {1 + \frac{3 (1 + ν_{0})^{2}}{(1 - 2 ν_{0})^{2}} \\ {[\frac{(1 - 2 ν_{0})}{3 G_{0} (1 - ν_{0}) + 2 (1 - 2 ν_{0}) η δ^{r}}]}^{2}} < 0 \end{array}

（24）

\frac{\partial \tilde{Y} (δ)}{\partial ν_{0}} = \frac{18 (1 - ν_{0}) (1 - 3 ν_{0}) G_{0}^{2} (G_{0} + η δ^{r})}{δ [3 G_{0} (1 - ν_{0}) + 2 (1 - 2 ν_{0}) η δ^{r}]} [1 + \frac{3 (1 - 2 ν_{0})^{2} (- 4 ν_{0}^{3} + 17 ν_{0}^{2} - 12 ν_{0} + 3)}{9 (1 - ν_{0}) (1 - 3 ν_{0})} (\frac{η δ^{r}}{G_{0}}) - \frac{4 ν_{0} (1 - 2 ν_{0})^{2}}{9 (1 - ν_{0}) (1 - 3 ν_{0})} {(\frac{η δ^{r}}{G_{0}})}^{2}]

（25）

式（22）的正负号通过放缩法判定，式（23）可以直接判断正负号，得出式（24）正负号的原因是δ远小于1。这三个不等式表明，松弛模量随瞬时剪切模量和Abel黏滞系数单调递增，随分数阶单调递减。这种规律性从力学概念上是极好理解的，瞬时剪切模量G₀和Abel黏滞系数η越大，材料越刚化，抵抗变形能力越强。而当r=0时， Abel黏壶代表理想弹簧体， r为1时对应理想牛顿流体，分数阶r从0至1变化，材料的黏塑性增强，导致松弛模量或强度减小，从这个角度可以论证上述关系是正确的。殷德顺等^［20］基于Holloman强度准则和分数阶微积分理论，分析得到单轴恒定应变率加载试验所得的岩石应变硬化指数n和分数阶r具有反比例关系，硬化指数n越大，材料的抗变形能力越强，再次说明式（24）推断正确。式（25）的正负号待定。取表1中的参数组合代入式（19），绘制松弛模量随r和ν₀的变化如图4所示。

表1 参数变化表

Table 1 The four parameters changing with r and ν₀

参数	仅r变化	仅ν₀变化
G₀/MPa	0.6	0.6
η/（MPa⋅s ^r ）	55	55
$r$	0.1 ~ 0.2	0.15
ν₀	0.3	0.3 ~ 0.45

新窗口打开| 下载CSV

图4

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

图4 松弛模量随参数变化曲线

Fig.4 The curves of relaxation modulus Y changing with fractional order and time （a） and initial Poisson’s ratio ν₀ and time （b）

图4从数值角度表明上述结论正确，且分数阶的微小变化使松弛模量数值波动较大，而泊松比初始值的变化对松弛模量数值影响较小，最大值和最小值相差在10%以内。以上关于参数影响的研究对试验数据分析过程具有理论指导作用。

4 试验过程与结果

4.1　样品制备

球模仪测试适用于无块状碎屑物质，且具有斑状、层状和细网状结构的冻结黏土、冻结粉土和冻结砂土^［1］。以标准砂制备试样，放入圆柱状容器中进行击实，然后利用真空机进行24 h抽真空至饱和。采用规格为ϕ61.8 mm×20 mm的环刀在圆柱体上切取土样，制备平行样，并用保鲜膜将试样严密包裹以防止水分散失。在球形模板试验前，先进行物理性质测试，测得土样密度为1.84 g⋅cm^-3，含水率为14.5%，土颗粒密度为2.65 g⋅cm^-3。

4.2　试验步骤

试样和球模仪均被放入恒温箱中，考虑到起始冻结温度的要求，恒温箱环境温度设定为-8.4 ℃，并保温24 h以上，以便土样冻结。考虑表层冻结冰层的影响和冻胀效应，之后将土样从恒温箱中取出，快速削平土样表面并用保鲜膜严密包裹。在土样要压的一面上边标记好接触点号。调节球形模板压头高度至稍微和土样表面接触，并用固定螺栓将压头锁定。为消除扰动，继续在恒温箱中冻结2 h以上。球形模板试验开始时，快速打开恒温箱并松动固定螺栓和关闭恒温箱舱门。用位移计以1 s为采样步长，记录前1 h的压入位移，在压入位移增长速度缓慢以后，改用步长为10 min记录数据。压入试验持续24 h左右，直到压入位移渐趋于稳定值。根据文献［1］，对冻土，球形模板的第15 min压入位移需要满足下式的要求。据此多次试验后，确定荷载F为50 N。

0.005 D \leq S_{15 m i n} \leq 0.05 D

（26）

式中： S_15min为球形模板第15 min的竖向压入位移（mm）； D为球形模板直径，试验时为22 mm。

对冻结黏土、冻结亚黏土和冻结砂的试验统计结果表明，可以通过动弹性模量换算获得静弹性模量（下式），而动泊松比与静泊松比之间无固定关系^［21］。

E_{s} = 0.06 + 0.16 E_{d} + 0.01 E_{d}^{2}

（27）

式中： E_s为冻土静弹性模量（GPa）。

4.3　试验结果

等效黏聚力按照式（1）计算，球形模板的试验结果经整理如图5所示，该图表明冻结砂的球形模板压入位移曲线具有衰减蠕变的特点。在开始试验的前3 h内，轴向压入位移变化较为剧烈。在这之后，压入位移的增长速率极为缓慢，在24 h后基本趋近于稳定值。

图5

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

图5 球形模板试验结果

Fig.5 Tests result of the spherical template indenter

5 分析与讨论

5.1　球形压痕试验数据分析与验证

球形模板压入试验接近点荷载加载试验，不可避免地存在离散性问题。为解决离散性问题，一方面需要至少在6个不同接触位置测试获得平均压入位移^［1］；另一方面，需要避开瞬时接触的不稳定阶段，可以把5 ~ 25 s内所测松弛模量当作瞬时松弛模量。另外，球模仪所施加外荷载是恒定值，测试时采用F=50 N，对应的球模直径D=22 mm（由俄罗斯国家标准ГΟCT 12248-96^［16］建议）。取式（10）展开式的前5项和式（11），对冻结砂轴向位移曲线进行拟合，所得结果如表2及图6（a）所示。可以看出，对试验前3 h压入位移曲线的拟合，线性Kelvin模型（r=1）不如非线性Kelvin模型（r=0.142）。对压入稳定阶段位移曲线的描述两者较为接近，但是线性Kelvin模型的结果比实测值略偏低，而非线性Kelvin模型的结果略偏高。

表2 冻结砂轴向位移曲线拟合结果

Table 2 The fitted results of the axial displacement curves of frozen sand

r	R²	G₀/MPa	η/（MPa⋅s ^r ）	ν₀
1	0.718	5.143	1 544.00	0.349
0.142	0.948	0.561	54.99	0.300

新窗口打开| 下载CSV

图6

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

图6 轴向压入位移与松弛模量随时间的变化

Fig.6 Axial pressing displacement （a） and relaxation modulus （b） changing with time

造成位移随时间非线性变化的原因为：开始段接触应力集中在小区域接触面范围内，引起接触面附近冰的压融，冰和薄膜水被挤出，起支撑作用的土颗粒相应发生位移。随着接触面积扩大，接触应力的集中程度逐渐减弱，应力在冻土内部发生扩散过程。因接触面处的冻土薄膜水被挤出，土颗粒被压实而导致局部抵抗变形能力提升，这是压入后期位移曲线更趋近稳定的原因。

压入过程中冻土复杂的黏弹塑性变形表现出非线性特点。线性Kelvin仅有有限个时间尺度，定量描述非线性变形过程的精度较差。相反，非线性Kelvin模型因为Abel黏滞系数的存在，相当于引入多个不同时间尺度的线性元件，能更合理地刻画冻土衰减蠕变的特点。

通过比较后，线性Kelvin模型比较差，因此图6（b）仅按照式（2）和式（19）给出冻结细砂松弛模量的预测值。冻结砂的泊松比一般取值在0.3 ~ 0.35^［19］。根据式（2）和式（19）计算松弛模量时，因泊松比取值不同产生的相对误差在10%以内。由此可取非线性Kelvin模型拟合得到的泊松比。

在运用Roman弹性球解答时，考虑到位移速率较小，因此每一步都视为准静态加载，这样以压入位移回代到式（2），便得到冻结细砂松弛模量的一系列离散点［图6（b）］。为了便于局部对比，在图6（b）中已经把Roman弹性球解答给出的松弛模量初值184.8 MPa省略。由此可知，除了初值不同外，弹性球Roman解答和非线性Kelvin模型预测的松弛模量随时间的变化规律较为一致。

对不同温度的冻结黏土，弹性球Roman解答已经被证明和常规三轴试验的结果较为接近^［7］。表3的数据可证实这点。单轴试验和三轴试验所得的结果比球模仪试验所得的结果略大。测试存在差异的原因不仅与土样、测试条件和加载方式有关，还与应变值密切相关。

表3 不同试验结果对比

Table 3 Comparison of different experimental results

试验土样	温度/℃	总含水率/%	试验方式	松弛模量/MPa	文献来源	备注
冻结中砂	-7.0	12.6	单轴压缩	245.0	［22］
冻结中砂	-7.0	18.4	单轴压缩	203.7	［22］
冻结细砂	-7.0	15.0	动三轴	752.0	［23］	由式（27）换算
冻结砂	-8.0	15.0 ~ 18.0		726.1	［24］	理论分析
冻结砂	-8.0	15.0 ~ 18.0		754.0	［25］	理论分析
冻结细砂	-8.4	14.5	球模仪	184.8	本文	由式（2）计算

新窗口打开| 下载CSV

常规三轴试验或单轴试验的冻土应变值在10^-3 ~ 10^-2［26］，以最大压入位移和土样高度之比作为冻土接触面中心点处轴向压应变的表征值，本文得到的冻土最大轴向压应变为0.03，与单轴和三轴的应变值较为接近。球形模板试验与单轴蠕变试验中，前者的冻土在接触区域产生局部的黏弹塑性变形，后者产生整体的黏弹塑性变形，两者具有相似之处。

5.2　单轴蠕变试验验证

5.2.1　算例一

该土为内蒙古某矿区的砂土，含水率为9.62%。重塑砂土分别在-5 ℃、 -10 ℃及-15 ℃温度条件下施加0.3倍、 0.5倍的单轴抗压强度进行恒定应力蠕变试验^［27］。取式（15）和式（21）的前五项作为计算基础，冻结砂的泊松比取0.3，分数阶取值在0 ~ 1，对它们取值区间进行限定并通过MATLAB自编程序计算，分析结果见图7。

图7

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

图7 -5 ℃、 -10 ℃、 -15 ℃条件下冻结砂单轴蠕变与松弛模量曲线

Fig.7 The uniaxial creep curves under -5 ℃ （a）， -10 ℃ （b）， -15 ℃ （c） and the relaxation modulus curves under -5 ℃ （d）， -10 ℃ （e）， -15 ℃ （f） of frozen sand

图7（a） ~ （c）表明，非线性Kelvin模型对不同温度、不同应力条件下冻结砂单轴蠕变实测曲线的预测能力较佳。因此，进一步将拟合获得的各参数（表4）回代式（19），得到松弛模量曲线。

表4 冻结砂的力学参数

Table 4 Mechanics parameters of the frozen sands

力学参数	-5 ℃		-10 ℃		-15 ℃
力学参数	σ=0.62 MPa	σ=1.05 MPa	σ=1.09 MPa	σ=1.82 MPa	σ=1.33 MPa	σ=2.22 MPa
G₀/MPa	2.395	3.251	4.383	7.460	1.170	7.486
η/（MPa⋅s ^r ）	71.68	91.18	151.40	137.80	187.90	98.74
r	0.074	0.062	0.079	0.069	0.205	0.066
ν₀	0.3	0.3	0.3	0.3	0.3	0.3

新窗口打开| 下载CSV

图7（d） ~ （f）表明在相同温度条件下，随着应力值增大，冻结砂的松弛模量增大。在各个温度和应力条件下，非线性Kelvin模型预测的冻结砂松弛模量随时间的变化规律和数值范围与球模仪试验的分析结果较一致。这表明该模型能为球形模板试验和单轴蠕变试验的数据分析过程提供统一、科学的计算方法。温度变化与应力变化对冻结砂松弛模量的影响相互耦合，规律不明显。

5.2.2　算例二

所取试验数据为含水率12%的冻结砂质黏土分别在-5 ℃、 -10 ℃及-15 ℃的恒定应力单轴蠕变试验数据^［28］。加载应力值分别是单轴抗压强度的0.3、 0.5及0.7倍。冻结黏性土的泊松比在0.35 ~ 0.4^［19］，按照算例一的相同步骤计算，结果如图8所示。力学参数见表5。

图8

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

图8 -5 ℃、 -10 ℃、 -15 ℃条件下冻结砂质黏土的单轴蠕变与松弛模量曲线

Fig.8 The uniaxial creep curves under -5 ℃ （a）， -10 ℃ （b）， -15 ℃ （c） and the relaxation modulus curves under -5 ℃ （d）， -10 ℃ （e）， -15 ℃ （f） of frozen sandy clay

表5 冻结砂质黏土的力学参数

Table 5 Mechanics parameters of the frozen sandy clay

温度/℃	轴向应力σ/MPa	G₀/MPa	η/（MPa⋅s ^r ）	r	ν₀
-5	1.4	14.210	1 098.0	0.155	0.35
	2.3	10.280	884.8	0.171	0.35
	3.4	1.495	655.5	0.135	0.35
-10	1.7	6.563	2 338.0	0.206	0.35
	2.9	14.140	2 126.0	0.253	0.35
	4.1	11.580	1 360.0	0.217	0.35
-15	2.6	14.620	2 588.0	0.234	0.35
	4.4	13.440	2 409.0	0.264	0.35
	6.2	12.740	1 367.0	0.211	0.35

新窗口打开| 下载CSV

将各力学参数回代式（19），得到冻结砂质黏土的松弛模量曲线如图8（d） ~ （f）所示。

图8（a） ~ （c）表明，非线性Kelvin模型对恒定应力冻土单轴蠕变的试验具有较强描述能力。由图8（d） ~ （f）可知，随着温度降低，冻结砂质黏土松弛模量的瞬时值增加，符合已有试验规律。在相同温度条件下，低应力下的冻结砂质黏土松弛模量值高于处在高应力下的值，该变化趋势和冻结砂的正好相反。这表明对不同的土质而言，其松弛模量是应力和温度的函数，不同土质的松弛模量随温度和应力的变化各不相同，很难找到共性的规律。

6 结论

本文基于半空间黏弹性理论和分数阶理论，获得了冻土压入位移和松弛模量的非线性Kelvin解，其特例为线性Kelvin解，并进行参数分析。提出将恒载单轴蠕变试验和恒载球形压痕蠕变试验统一分析的方法，可以实现两种试验所得松弛模量函数分析结果的相互验证。随后以冻结细砂、冻结砂及冻结砂质黏土为例，对比了球形压痕试验及单轴蠕变等试验的分析结果。得到以下结论：

（1）非线性Kelvin模型具有众多优势，可用于基于球模仪试验和单轴蠕变试验的冻土松弛模量提取。该模型中的参数最少且物理意义明确，其解答所得的松弛模量函数随参数变化具有较佳单调特性，随瞬时剪切模量和Abel黏滞系数单调递增，随分数阶单调递减，而受泊松比变化的影响较小，泊松比在经验值范围内变化时，最大值和最小值相差在10%以内。这便于进一步建立各参数随冻土的温度、含水率等物理参数的关系，由此直接用冻土的各物理参数估算松弛模量。

（2）对于冻结细砂，非线性Kelvin模型比线性Kelvin模型更吻合试验所得的位移-时间曲线，所给出的冻结细砂松弛模量变化规律和弹性球解答的基本相同。特别是对工程上较为关注的冻土长期松弛模量，能获得较好的预测值。

（3）非线性Kelvin模型预测值十分吻合冻结砂及冻结砂质黏土的单轴蠕变曲线，不同土质的松弛模量随温度和应力的变化规律各不相同，基于单轴蠕变试验的冻结砂松弛模量随时间的变化规律和数值范围与球形模板试验的较为一致，冻结砂质黏土的松弛模量随应力和温度的变化规律再现了已有研究结论。这表明非线性Kelvin模型解答能为单轴蠕变试验和球模仪试验提供科学、统一的冻土松弛模量估算方法，便于实现基于这两种试验的分析数据之间的相互对比和校正。

本文研究的不足之处在于分数阶黏弹性力学模型是基于试验现象，且不考虑球形模板与土样之间复杂的接触摩擦行为，这需要在后续研究中进一步完善。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Zhang

， Ma

Wei

， Zhang

Zhongqiong

， et al.

Application of spherical template indenter to long-term strength tests for frozen soil

［J］. Rock and Soil Mechanics， 2012， 33（11）： 3516 -3520.

[本文引用: 7]

张泽，马巍，张中琼，等.

球形模板压入仪在冻土长期强度测试中的应用

［J］. 岩土力学， 2012， 33（11）： 3516 -3520.

[本文引用: 7]

[2]

Roman

L T

， Zhang

Effect of cycles of freezing and thawing on the physical and mechanical properties of moraine loam

［J］. Soil Mechanics and Foundation Engineering， 2010， 47（3）： 96 - 101.