岩石疲劳应力等效化及非线性疲劳变形本构模型

唐欣; 俞缙; 林立华; 高海东; 李刚; 林植超

doi:10.11779/CJGE202101012

岩石疲劳应力等效化及非线性疲劳变形本构模型 English Version

唐欣^1,,
俞缙^1, ,,
林立华²,
高海东³,
李刚³,
林植超²

1.
华侨大学福建省隧道与城市地下空间工程技术研究中心,福建　厦门　361021
2.
厦门路桥工程投资发展有限公司,福建　厦门　361026
3.
中铁十八局集团有限公司,天津　300222

基金项目:

国家自然科学基金项目 51874144

国家自然科学基金项目 51679093

国家自然科学基金项目 42077254

华侨大学研究生科研创新能力培育计划项目 18014086024

详细信息

作者简介:
唐欣(1993— ),男,江西抚州人,硕士研究生,主要从事岩石力学研究工作。E-mail: Tx14046341@outlook.com

通讯作者:
俞缙, E-mail: bugyu0717@163.com

中图分类号: TU43
计量
- 文章访问数: 314
- HTML全文浏览量: 23
- PDF下载量: 138
出版历程
- 收稿日期: 2020-05-27
- 网络出版日期: 2022-12-04
- 刊出日期: 2020-12-31

Equivalent fatigue stress and non-linear constitutive model for fatigue of rock

1.
Fujian Research Center for Tunneling and Urban Underground Space Engineering, Huaqiao University, Xiamen361021, China
2.
Xiamen Road & Bridge Engineering Investment and Development Co., Ltd., Xiamen361026, China
3.
China Railway 18th Bureau Group Co., Ltd., Tianjin 300222, China

摘要

摘要: 传统的岩石疲劳本构模型对循环加卸载次数敏感性方面存在明显不足,无法准确描述残余变形与疲劳损伤进程的关系。在蠕变与疲劳试验基础上,通过对岩石蠕变和疲劳的时间（循环次数）与变形归一化,分析验证了两者间存在某种“同源性”。运用流变理论描述疲劳荷载下岩石变形全过程,基于等寿命Gerber曲线方程与Goodman直线方程,提出了通过疲劳寿命与应力关系的方式,定义了一种新的疲劳应力等效方法。用村山黏弹塑性体替换Kelvin黏弹性体,并加入开关函数和岩石完整度参数对Bingham黏塑性体进行改造,建立了一种新的7参数岩石非线性疲劳变形本构模型。该模型能较好地表征疲劳荷载下岩石减速、等速和加速疲劳变形的全过程,其加速疲劳变形段的完整度参数随着循环次数N增加而减小,能够体现疲劳损伤对加卸载次数因素的敏感性。
- 疲劳变形 /
- 应力等效 /
- 西原模型 /
- 村山体 /
- 非线性本构模型
Abstract: The traditional damage model is insensitive to cycles of loading and unloading, and it is obviously inadequate for the characterization of fatigue deformation and failure of rock. Based on the creep and fatigue tests, the "homology" hypothesis of the two kinds of deformation is verified by normalizing the time (cycle number) and deformation of rock creep and fatigue. The rheological theory is used to describe the whole process of fatigue deformation of rock under fatigue loads, and susequently the improved fatigue stress equivalence is put forward based on the equal life curve equation, Gerber and Goodman linear equation. The Kelvin visco-elastomer is replaced by the Murayama visco-elastomer, and the switch function and rock integrity parameters are added to non-linear treatment of Bingham viscoplastic body to establish a 7-parameter non-linear fatigue model for rock. The results show that the model can characterize the whole process of rock deceleration, constant velocity and accelerated fatigue deformation under the fatigue loads. The integrity parameter decreases as the number of cycles, N, increases, and this trend reflects the sensitivity of the fatigue damage to factors of loading and unloading times. The research put a foundation for the long-term deformation prediction of engineering rock mass.
- fatigue deformation /
- stress equivalence /
- Nishihara model /
- Murayama /
- non-linear constitutive model

HTML全文

0. 引言

饱和土土颗粒和孔隙流体常具有热胀冷缩的特性,由于孔隙流体的热膨胀系数通常大于土颗粒的热膨胀系数,温度升高会导致土体内产生超静孔压,超静孔压的消散与土体的热胀冷缩都将导致土体的变形,这一现象在岩土工程中常常被忽略。但随着中国不断推进地热资源开发和地下能源存储等工程的建设,温度场影响下的岩土体热–水–力耦合响应分析已成为研究的热点。

许多学者针对温度变化引起的饱和弹性土体热固结响应展开了解析研究。Booker等^[1-2]最早系统进行了相关工作,如通过Laplace变换、Fourier变换先后提出了点热源和球热源周围土体的热固结解析解。Mctigue^[3]将由热传导方程得到的温度场解代入渗流场方程,得到了半无限空间饱和土体表面透水或不透水时,土体内部超静孔压在表面恒定温度荷载作用下的热响应解。Giraud等^[4]建立了内部存在指数衰减热源的半无限空间饱和土体简化计算模型,并通过Laplace变换求解得到了温度、超静孔压和位移随深度的分布。Smith等^[5-6]提出了力场–渗流场–温度场全耦合的饱和土固结理论,并提出了在Laplace变换域上解决热弹性固结问题的边界元法。Bai等^[7]通过有限Fourier变换给出了单侧排水的单层饱和土体热固结解,分析了控制方程非耦合、部分耦合及全耦合情况下温度、超静孔压与位移的变化区别。Zhou等^[8-9]建立了考虑热渗效应与等温热流效应的热固结全耦合控制方程,并用Laplace变换给出了线热源、柱热源与球热源作用下的饱和土体热固结解析解。Blond等^[10]考虑了土颗粒和孔隙流体的可压缩性,给出了半无限空间土体在表面简谐温度作用下超静孔压变化最大值的解,并确定了孔压最大值所在位置。白冰^[11-13]利用Laplace变换和数值逆变换,给出了半无限饱和介质在周期变化温度作用时,土体内部一维热响应的半解析解。Liu等^[14-15]建立了热–水–力耦合的动力形式控制方程,并采用Laplace变换得到了轴对称和球对称情况的土体热动力响应。He等^[16]与Yang等^[17]分别假设土骨架与孔隙流体处于局部非热平衡态,用Laplace变换的方法求解土体内渗流场与温度场分布,并与局部热平衡态的热固结解进行对比。Ai等^[18-20]、Wang等^[21-22]在Laplace变换域引入解析层元法,对层状多孔介质的热–水–力耦合响应进行了大量研究。

从前人的研究了解到,饱和土热固结的解析解可用于揭示饱和弹性土体受温度荷载作用的固结响应规律,也可用于验证相关数值解的正确性,具有十分重要的意义。但学者们广泛应用了Laplace变换对热固结问题进行求解,复杂边界条件下得到的Laplace变换域解无法直接进行逆变换得到其对应的时域表达式,数值逆变换的参数又常难以确定。此外,学者们对半无限空间饱和土体在季节性温度荷载作用下的一维热固结问题研究较少,表面透水或不透水对热固结的影响未进行对比分析。因此,本文借鉴钮家军等^[23]提出的多场耦合方程求解方法,采用正余弦变换与Hamilton–Cayley定理,给出两类非齐次边界条件下,半无限饱和土体一维热固结精确时域显式积分解的求解方法,避免了Laplace数值逆变换可能带来的误差。最后通过算例验证本文解的正确性,并对比研究表面透水和不透水的半无限空间饱和土体在季节性温度荷载作用下的一维热固结响应。

1. 饱和土热–水–力耦合控制方程

Smith等^[5-6]曾提出饱和土热–水–力耦合的平衡方程、质量守恒方程与能量守恒方程；Zhou等^[8]进一步提出考虑了热渗效应、等温热流效应及物性参数变化的饱和多孔介质全耦合理论模型；Coussy^[24]从热力学第一定律和第二定律角度,推导得到了对称形式的饱和多孔介质耦合控制方程；Bai等^[25-26]考虑了渗流过程对热传导的影响,给出了饱和土多场全耦合的控制方程。本文首先给出土体的本构关系及修正的Darcy定律和Fourier定律：

$σ_{i j} = 2 G ε_{i j} + λ ε_{v} δ_{i j} - 3 K_{0} α_{s} T δ_{i j} - α p δ_{i j}$ ,

(1)

$q_{w} = - \frac{k_{w}}{γ_{w}} \nabla p - S_{w} \nabla T$ ,

(2)

$q_{T} = - S_{T} \nabla p - K_{T} \nabla T$ 。

(3)

式中 $ε_{v} = ε_{x} + ε_{y} + ε_{z} ； K_{0} = λ + \frac{2}{3} G ； α = 1 - \frac{K_{0}}{K_{s}} ； S_{T} \approx$ $T_{0} S_{w}$ ； $σ_{i j}$ 为总应力；λ,G为Lame常数； $ε_{v}$ 为体应变；K₀为土骨架体积模量； $α_{s}$ 为土颗粒的线热膨胀系数；p为土体内超静孔压；T为区别于初始温度T₀的土体温度改变量；α为Biot系数； $γ_{w}$ 为孔隙流体的重度；k_w和K_T分别为土体的渗透系数与热传导系数；S_w和S_T分别为土体的热渗效应系数与等温热流效应系数。

根据修正的Darcy定律可知,当式（2）右侧两项的差值不能忽略不计时,即土渗透系数较小或温度梯度较大时,应当考虑热渗效应的影响,这与Carnahan^[27]、Horseman等^[28]的研究结论一致。等温热流效应同理。将考虑热渗效应、等温热流效应的Darcy定律和Fourier定律代入Coussy^[24]提出的饱和多孔介质热–水–力耦合控制方程,可以得到

$G \nabla^{2} u + (λ + G) \nabla ε_{v} - α \nabla p - β \nabla T = 0$ ,

(4)

$\frac{k_{w}}{γ_{_{w}}} \nabla^{2} p + S_{w} \nabla^{2} T = α \frac{\partial ε_{v}}{\partial t} + α_{p} \frac{\partial p}{\partial t} - Y \frac{\partial T}{\partial t}$ ,

(5)

$\frac{S_{T}}{T_{0}} \nabla^{2} p + \frac{K_{T}}{T_{0}} \nabla^{2} T = β \frac{\partial ε_{v}}{\partial t} - Y \frac{\partial p}{\partial t} + \frac{c_{v}}{T_{0}} \frac{\partial T}{\partial t}$ 。

(6)

其中,

$β = 3 α_{s} K_{0}, Y = 3 (α - n) α_{s} + 3 n α_{w}$ ,

$α_{p} = n / K_{w} + (α - n) / K_{s}, c_{v} = (1 - n) ρ_{s} c_{s} + n ρ_{w} c_{w}$ 。

式中 u为土骨架的位移；t为时间；n为土体的初始孔隙率；K_s和K_w分别为土颗粒和孔隙流体的体积模量；ρ_s和ρ_w分别为土颗粒和孔隙流体的密度；α_s和α_w分别为土颗粒和孔隙流体的线热膨胀系数；c_s和c_w分别为土颗粒和孔隙流体的质量比热容。

考虑一维情况,将式（4）左右两侧分别积分,可以得到

$ε_{v} = \nabla \cdot u = \frac{α p + β T + q (t)}{λ + 2 G} = \frac{α p + β T + q (t)}{E_{c}}$ 。

(7)

式中 q(t)为关于时间的函数,由位移边界确定,常与外荷载有关；ε_v为体应变,以拉为正。

将式（7）代入式（5）和（6）,可以得到以超静孔压p与温度改变量T为未知量的饱和土体热–水–力耦合控制方程,其一维形式为

$\begin{array}{l} \frac{k_{w}}{γ_{w}} \frac{\partial^{2} p}{\partial x^{2}} + S_{w} \frac{\partial^{2} T}{\partial x^{2}} \\ = (\frac{α^{2}}{E_{c}} + α_{p}) \frac{\partial p}{\partial t} + (\frac{α β}{E_{c}} - Y) \frac{\partial T}{\partial t} + \frac{α}{E_{c}} \frac{\partial q (t)}{\partial t} \end{array}$ ,

(8)

$\begin{array}{l} \frac{S_{T}}{T_{0}} \frac{\partial^{2} p}{\partial x^{2}} + \frac{K_{T}}{T_{0}} \frac{\partial^{2} T}{\partial x^{2}} \\ = (\frac{α β}{E_{c}} - Y) \frac{\partial p}{\partial t} + (\frac{β^{2}}{E_{c}} + \frac{c_{v}}{T_{0}}) \frac{\partial T}{\partial t} + \frac{β}{E_{c}} \frac{\partial q (t)}{\partial t} \end{array}$ 。

(9)

为便于表示,将式（8）,（9）简写为

$g_{11} \frac{\partial^{2} p}{\partial x^{2}} + g_{12} \frac{\partial^{2} T}{\partial x^{2}} + c_{11} \frac{\partial p}{\partial t} + c_{12} \frac{\partial T}{\partial t} = q_{1} (t)$ ,

(10)

$g_{21} \frac{\partial^{2} p}{\partial x^{2}} + g_{22} \frac{\partial^{2} T}{\partial x^{2}} + c_{21} \frac{\partial p}{\partial t} + c_{22} \frac{\partial T}{\partial t} = q_{2} (t)$ 。

(11)

式中, $g_{11} = \frac{k_{w}}{γ_{w}}, g_{12} = S_{w}, c_{11} = - \frac{α^{2}}{E_{c}} - α_{p}, c_{12} = Y - \frac{α β}{E_{c}}$ , $g_{21} = \frac{S_{T}}{T_{0}}, g_{22} = \frac{K_{T}}{T_{0}}, c_{21} = Y - \frac{α β}{E_{c}}, c_{22} = - \frac{β^{2}}{E_{c}} - \frac{c_{v}}{T_{0}}$ 。

2. 半无限空间饱和土一维热固结

2.1 半无限空间饱和土一维热固结定解问题

为研究季节性温度变化对半无限空间饱和土体固结的影响,假设土体表面作用有大范围的均布温度荷载,饱和土的温度传导、渗流与变形均发生在沿深度方向,简化为一维热固结问题。设土体表面为坐标原点,坐标轴x方向垂直向下,建立如图1所示的半无限空间饱和土一维热固结计算模型。

图 1 半无限空间饱和土一维热固结计算模型

Figure 1. Model for one-dimensional thermal consolidation of semi-infinite saturated soils

下载: 全尺寸图片幻灯片

设土体位移、超静孔压与温度变化量的初始条件与无穷远处的自然边界条件分别为

$u (x, 0) = 0, p (x, 0) = 0, T (x, 0) = 0$ ,

(12)

$u (\infty, t) = 0, p (\infty, t) = 0, T (\infty, t) = 0$ 。

(13)

考虑上表面为自由边界,超静孔压和温度变化量为以下两类边界条件：

$p (0, t) = f_{1} (t), T (0, t) = f_{2} (t) (边界条件 (A))$ ,

(14)

$p_{, x} (0, t) = f_{1} (t), T_{, x} (0, t) = f_{2} (t) (边界条件 (B))$ 。

(15)

式中,f₁(t)和f₂(t)为任意关于时间t的函数。边界条件（A）代表已知半无限空间土体表面超静孔压p和温度改变量T随时间变化的工况；根据修正的Darcy定律和Fourier定律,边界条件（B）代表已知半无限空间土体表面渗流流量q_w和热流流量q_T随时间变化的工况。

式（10）～（15）分别给出了饱和土热固结的控制方程,边界条件和初始条件,构成了半无限空间饱和土一维热固结的定解问题。

2.2 求解过程

采用正弦变换和余弦变换方法对上述定解问题进行求解,正弦变换、余弦变换及其相应逆变换的定义分别为

$\begin{array}{l} s (ω) = \int_{0}^{\infty} f (x) \sin (ω x) d x, \\ f (x) = \frac{2}{π} \int_{0}^{\infty} s (ω) \sin (ω x) d ω 。 \end{array}}$

(16)

$\begin{array}{l} c (ω) = \int_{0}^{\infty} f (x) \cos (ω x) d x, \\ f (x) = \frac{2}{π} \int_{0}^{\infty} c (ω) \cos (ω x) d ω 。 \end{array}}$

(17)

（1）边界条件(A)的求解

对超静孔压p和温度改变量T实行正弦变换,定义为

${\bar{p}}_{s} (ω, t) = \int_{0}^{\infty} p (x, t) \sin (ω x) d x$ ,

(18)

${\bar{T}}_{s} (ω, t) = \int_{0}^{\infty} T (x, t) \sin (ω x) d x$ 。

(19)

根据正弦变换的性质,结合边界条件,可将式（10）,（11）转化为如下形式：

$\begin{matrix} g_{11} [- ω^{2} {\bar{p}}_{s} (ω, t) + ω p (0, t)] + g_{12} [- ω^{2} {\bar{T}}_{s} (ω, t) + \\ ω T (0, t)] + c_{11} \frac{\partial {\bar{p}}_{s} (ω, t)}{\partial t} + c_{12} \frac{\partial {\bar{T}}_{s} (ω, t)}{\partial t} = {\bar{q}}_{s1} \end{matrix}$ ,

(20)

$\begin{array}{l} g_{21} [- ω^{2} {\bar{p}}_{s} (ω, t) + ω p (0, t)] + g_{22} [- ω^{2} {\bar{T}}_{s} (ω, t) + \\ ω T (0, t)] + c_{21} \frac{\partial {\bar{p}}_{s} (ω, t)}{\partial t} + c_{22} \frac{\partial {\bar{T}}_{s} (ω, t)}{\partial t} = {\bar{q}}_{s2} \end{array}$ 。

(21)

式（20）,（21）记为矩阵微分方程形式,即

$C \frac{d X (t)}{d t} = ω^{2} G X (t) + H_{s} (t)$ ,

(22)

式中,

$X (t) = {\begin{matrix} {\bar{p}}_{s} (ω, t) \\ {\bar{T}}_{s} (ω, t) \end{matrix}}, C = [\begin{matrix} c_{11} & c_{12} \\ c_{21} & c_{22} \end{matrix}], G = [\begin{matrix} g_{11} & g_{12} \\ g_{21} & g_{22} \end{matrix}]$

$H_{s} = {\begin{matrix} {\bar{q}}_{s1} - g_{11} ω f_{1} (t) - g_{12} ω f_{2} (t) \\ {\bar{q}}_{s 2} - g_{21} ω f_{1} (t) - g_{22} ω f_{2} (t) \end{matrix}}$ 。

对应的初始条件的正弦变换为

${\bar{p}}_{s} (ω, 0) = 0, {\bar{T}}_{s} (ω, 0) = 0$ 。

(23)

（2）边界条件(B)的求解

与边界条件(A)的求解类似,对超静孔压p和温度改变量T进行余弦变换,记为

${\bar{p}}_{c} (ω, t) = \int_{0}^{\infty} p (x, t) \cos (ω x) d x$ ,

(24)

${\bar{T}}_{c} (ω, t) = \int_{0}^{\infty} T (x, t) \cos (ω x) d x$ 。

(25)

根据余弦变换的性质,结合边界条件,可将式（10）,（11）转化为如下形式：

$\begin{array}{l} g_{11} [- ω^{2} {\bar{p}}_{c} (ω, t) - p_{, x} (0, t)] + g_{12} [- ω^{2} {\bar{T}}_{c} (ω, t) - \\ T_{, x} (0, t)] + c_{11} \frac{\partial {\bar{p}}_{s} (ω, t)}{\partial t} + c_{12} \frac{\partial {\bar{T}}_{s} (ω, t)}{\partial t} = {\bar{q}}_{c1} \end{array}$ ,

(26)

$\begin{array}{l} g_{21} [- ω^{2} {\bar{p}}_{c} (ω, t) - p_{, x} (0, t)] + g_{22} [- ω^{2} {\bar{T}}_{c} (ω, t) - \\ T_{, x} (0, t)] + c_{21} \frac{\partial {\bar{p}}_{s} (ω, t)}{\partial t} + c_{22} \frac{\partial {\bar{T}}_{s} (ω, t)}{\partial t} = {\bar{q}}_{c2} \end{array}$ 。

(27)

将式（26）,（27）记为矩阵微分方程形式,即

$C \frac{d X (t)}{d t} = ω^{2} G X (t) + H_{c} (t)$ ,

(28)

式中,

$X (t) = {\begin{matrix} {\bar{p}}_{c} (ω, t) \\ {\bar{T}}_{c} (ω, t) \end{matrix}}, C = [\begin{matrix} c_{11} & c_{12} \\ c_{21} & c_{22} \end{matrix}], G = [\begin{matrix} g_{11} & g_{12} \\ g_{21} & g_{22} \end{matrix}]$ ,

$H_{c} = {\begin{matrix} {\bar{q}}_{c1} + g_{11} f_{1} (t) + g_{12} f_{2} (t) \\ {\bar{q}}_{c 2} + g_{21} f_{1} (t) + g_{22} f_{2} (t) \end{matrix}}$ 。

对应的初始条件的余弦变换为

${\bar{p}}_{c} (ω, 0) = 0, {\bar{T}}_{c} (ω, 0) = 0$ 。

(29)

（3）常微分方程的解

式（22）,（28）均为关于时间t的常微分方程,其可统一表示为标准形式：

$\frac{d}{d t} X (t) = ω^{2} N X (t) + H (t)$ ,

(30)

式中,N=C^-1G,对边界条件(A),H=C^-1H_s；对边界条件(B),H=C^-1H_c。

方程（30）有以下形式的解：

$X (t) = e^{ω^{2} N t} X (0) + e^{ω^{2} N t} \int_{0}^{t} e^{- ω^{2} N ξ} H (ξ)d ξ$ 。

(31)

根据Hamilton-Cayley定理, $e^{N t}$ 写为如下形式：

$e^{N t} = α_{0} (t) I + α_{1} (t) N$ ,

(32)

其中,

$α_{0} (t) = \frac{r_{1} e^{r_{2} t} - r_{2} e^{r_{1} t}}{r_{1} - r_{2}}, α_{1} (t) = \frac{e^{r_{1} t} - e^{r_{2} t}}{r_{1} - r_{2}}$ ,

(33)

式中,r₁和r₂为矩阵N的特征值。当r₁=r₂时,其表达式应为

$α_{0} (t) = (1 - r_{1} t) e^{r_{1} t}, α_{1} (t) = t e^{r_{1} t}$ 。

(34)

因此,将初始条件（23）,（29）, $e^{N t}$ 与H(t)代入式（31）,逐步计算即可求得超静孔压p和温度改变量T在变换域上的解,再根据相对应的正弦逆变换（16）或余弦逆变换（17）即可得到其对应的时域精确显式积分解。将超静孔压p和温度改变量T的解代入式（7）,即可得到土体体应变 $ε_{v}$ 的精确显式积分解。将体应变沿深度方向积分,即可得到饱和土体各深度处的位移解u。

3. 算例分析

3.1 饱和土体热–水–力耦合物性参数

假定土体的初始温度T₀=300 K,表面外荷载q(t) =0。参考Bai^[11]给出的土体物性参数,本文土体的热–水–力基本参数取值见表1。

表 1 饱和土体的材料参数

Table 1. Material parameters of saturated soils

孔隙率n	泊松比μ	土体弹性模量E/MPa	土颗粒体积模量K_s/GPa	土颗粒密度ρ_s/(kg·m^-3)	土颗粒比热容c_s/(J·kg^-1·K^-1)	土颗粒热膨胀系数α_s/K^-1	孔隙流体体积模量K_w/GPa	孔隙流体密度ρ_w/(kg·m^-3)	孔隙流体比热容c_w/(J·kg^-1·K^-1)	孔隙流体热膨胀系数α_w/K^-1	土体热传导系数K_T/(J·s^-1·m^-1·K^-1)	热渗系数S_w/(m²·s^-1·K^-1)
0.4	0.3	5	20	2600	800	1.5×10^-5	5	1000	4200	2×10^-5	1	5×10^-11

下载: 导出CSV

| 显示表格

因为土体某深度处温度升高越快,渗流越慢,超静孔压响应越大；相反,温度升高越慢,渗流越快,超静孔压响应越小,所以饱和土中的超静孔压响应是渗流和温升快慢的耦合作用结果。Bai^[12]在研究中引入了无量纲固结系数与热扩散系数的比值c/κ来表示两者的相对速度。本文类似引入简化的无量纲相对扩散数 $v_{1} / v_{2}$ ,其中v₁为太沙基一维固结理论定义的固结系数, $v_{1}$ =k_wE_s/ $γ_{w}$ ； $v_{2}$ 为根据热传导方程定义的温度变化系数, $v_{2}$ =K_T/c_v。

3.2 解析解正确性验证

Mctigue^[3]提出直接求解非耦合的热传导方程,得到表面温度突增时,半无限饱和土体内的温度响应,再将温度场的解代入渗流场方程,得到土体内超静孔压的解析解。为验证本文解析解的正确性,将本文解析方法计算得到的土体热固结响应结果与Mctigue^[3]方法的结果进行对比。所以,将本文温度场控制方程简化为热传导方程,土体的初始条件和无穷远处的边界条件不变,表面边界条件设为透水边界,施加突增温度荷载：

$p (0, t) = 0, T (0, t) = 25$ 。

(35)

图2,3给出了分别采用本文的积分解方法和Mctigue^[3]方法得到的不同时刻土体温度改变量T和超静孔压p随深度的分布。由图3可以看出,两种方法计算得到的结果几乎相同,所以本文解析解的正确性可以得到验证。图2,3还反映了突增温度荷载作用下的半无限空间饱和土体的一维热固结响应特性,即温度从表面向深处传播,较深处温度增长比较浅处慢；各深度处温度上升引起的超静孔压峰值几乎相同,并呈波动特征向土体深处传递。

图 2 表面温度突增时不同时刻温度分布图

Figure 2. Distribution of temperature change along depth at different moments subjected to step temperature increase

下载: 全尺寸图片幻灯片

图 3 表面温度突增时不同时刻超静孔压分布图

Figure 3. Distribution of excess pore pressure along depth at different moments subjected to step temperature increase

下载: 全尺寸图片幻灯片

3.3 半无限空间饱和土表面季节性温度荷载

根据半无限空间饱和土体表面是否有覆盖层分为两类边界条件。当土体表面无覆盖层时,土体表面为透水边界,由于土体表面与大气直接接触,受阳光直射,故假定土体表面温度变化与气温相似,为分析季节性周期变化的温度荷载对饱和土体的热固结影响,取土体表面温度周期变化,幅值为25℃,周期为T_p=π×10⁷ s,约为1 a,即边界条件(A)。当土体表面有覆盖层时,如机场跑道、高等级路面等表面有道面板等不透水结构,设定土体表面为不透水边界,又因为覆盖层的存在,外界温度变化无法直接改变土体表面温度,而是通过热量向下传递的形式产生影响,因此根据Fourier定律,假定土体上表面温度荷载以周期变化的热流量表示,即边界条件(B),热流量变化周期与边界条件(A)相同,为T_p=π×10⁷ s,约为1 a。为分析季节性温度荷载作用下,表面渗透条件对饱和土热固结响应的影响,通过计算使热流边界条件(B)情况下得到的土体表面温度变化响应与边界条件(A)表面周期变化温度近似相同：

$p (0, t) = 0, T (0, t) = 25 \sin (2 t / 10^{7}) (边界条件 (A))$ ,

(36)

$p_{x} (0, t) = 0, K_{T} T_{x} (0, t) = - 20 \cos (2 t / 10^{7}) (边界条件 (B)) 。$

(37)

3.4 半无限空间饱和土体一维热固结特性分析

（1）边界条件(A)

图4,5分别给出了半无限饱和土表面透水时,不同深度处温度改变量T和超静孔压p的变化过程。由于表面温度周期性变化,土体不同深度处的温度改变量T和超静孔压p均呈周期性变化,从温度荷载施加后的瞬态响应渐渐过渡到稳定状态,其稳态的变化周期与表面温度周期相同。随着深度增加,温度改变量T变化的幅值减小,且有滞后的趋势。深度0 m处温度变化幅值为25 K,深度5 m处的温度变化不超过±1.98 K。

图 4 边界条件(A)不同深度处温度变化过程

Figure 4. Evolution of temperature change with time at different depths under boundary condition A

下载: 全尺寸图片幻灯片

图 5 边界条件(A)不同深度处超静孔压变化过程

Figure 5. Evolution of excess pore pressure with time at different depths under boundary condition A

下载: 全尺寸图片幻灯片

由于表面为透水边界,温度变化产生的超静孔压可以向两个方向消散,超静孔压p的变化幅值沿深度先增大后减小。深度2 m处,超静孔压稳态变化幅值最大约为6.76 kPa。图6还给出了表面处(x=0)位移随时间的变化,显然位移亦随时间周期性波动,表面最大位移可达4.4 mm。

图 6 边界条件(A)表面位移变化过程

Figure 6. Evolution of surface displacement of soils with time under boundary condition A

下载: 全尺寸图片幻灯片

（2）边界条件(B)

图7,8分别给出了半无限空间饱和土表面为边界条件（B）时,不同深度处温度改变量T和超静孔压p的变化过程。由图7,8可知,土体不同深度处的温度改变量T和超静孔压p受表面周期变化的热流影响而呈现周期变化的特性,稳态变化的周期与表面热荷载周期相同。不同于边界条件（A）的结果,土体温度改变量T和超静孔压p周期变化的幅值都随着深度增大而减小。深度0 m处温度稳态周期变化的幅值最大,约为26.1 K,较边界条件（A）深度0 m处幅值25 K高约4.4%。而深度0 m处超静孔压稳态周期变化的幅值最大,约为16.46 kPa,为边界条件（A）深度2 m处稳态幅值6.76 kPa的2.43倍。原因是深度0 m处温度波动幅值最大,且表面为不透水边界,温度产生的超静孔压只能向无穷远处消散,排水路经长导致了表面的超静孔压相对较大。图9反映了表面处（x=0）位移随时间变化的趋势,表面位移随时间周期性变化,变化最大值约为8.18 mm,比边界条件（A）表面的位移最大值高出85.9%。

图 7 边界条件(B)不同深度处温度变化过程

Figure 7. Evolution of temperature change with time at different depths under boundary condition B

下载: 全尺寸图片幻灯片

图 8 边界条件(B)不同深度处超静孔压变化过程

Figure 8. Evolution of excess pore pressure with time at different depths under boundary condition B

下载: 全尺寸图片幻灯片

图 9 边界条件(B)表面位移变化过程

Figure 9. Evolution of surface displacement of soils with time under boundary condition B

下载: 全尺寸图片幻灯片

综合边界条件（A）和边界条件（B）的结果可知,虽然两种边界条件下的土体温度变化特性类似,但表面不透水边界较透水边界,土体浅层的超静孔压更难以消散,孔压的累积更有可能导致浅层土体发生强度破坏,在机场道路,高等级路面等地基设计与施工过程中应引起重视。

（3）扩散数 $v_{1} / v_{2}$ 的影响

由于在表面不透水情况下,浅层土体更易产生较大的超静孔压和变形,因此本部分保持 $v_{2}$ 不变,采用不同的固结系数 $v_{1}$ ,对边界条件（B）情况下,土体的热固结响应进行进一步研究。图10～12分别给出了不同扩散数 $v_{1} / v_{2}$ 条件下,深度0 m处土体温度改变量T、超静孔压p和位移u的响应。可以发现不同的 $v_{1}$ 对表面温度变化和位移几乎没有影响,但扩散数增大5倍后,超静孔压p的峰值16.46 kPa减小为9.55 kPa,降低了约41.2%；扩散数减小10倍后,超静孔压p的峰值增大至27.48 kPa,提高了约66.9%。不同 $v_{1}$ 对表面位移没有影响的原因在于同样的温度变化产生的总超静孔压是相同的, $v_{1}$ 的不同仅导致了超静孔压沿深度方向更加均匀的分布,如图13。由于假设了土体沿深度方向的压缩模量E_s不变,上表面位移是体应变 $ε_{v}$ 沿深度方向的积分,所以扩散数对土体表面位移几乎无影响。

图 10 边界条件(B)不同扩散数表面温度随时间变化

Figure 10. Evolution of temperature of surface with time with different coefficients under boundary condition B

下载: 全尺寸图片幻灯片

图 11 边界条件(B)不同扩散数表面超静孔压随时间变化

Figure 11. Evolution of excess pore pressure of surface with time with different coefficients under boundary condition B

下载: 全尺寸图片幻灯片

图 12 边界条件(B)不同扩散数表面位移随时间变化

Figure 12. Evolution of displacement of surface with time with different coefficients under boundary condition B

下载: 全尺寸图片幻灯片

图 13 边界条件(B)不同扩散数超静孔压随深度变化

Figure 13. Distribution of excess pore pressure along depth with different coefficients under boundary condition B

下载: 全尺寸图片幻灯片

4. 结论

本文先给出了考虑热渗效应和等温热流效应的饱和土热–水–力耦合固结控制方程,再利用正弦变换和余弦变换,提出了两类任意非齐次边界条件下,半无限空间饱和土体一维热固结响应的精确积分形式解的求解方法,避免了复杂边界Laplace数值逆变换带来的误差。通过算例验证了本文解的正确性,分析了表面透水和不透水的半无限空间饱和土体在季节性周期变化温度荷载作用下的热响应,得到以下4点结论。

（1）土体内部温度从瞬态响应逐渐过渡到稳态周期波动,波动的周期与表面温度荷载相同。温度波动的幅值随深度增加而减小,且深处温度变化滞后于浅处温度变化。

（2）土体内部超静孔压响应呈周期性波动,变化的周期与表面温度荷载相同。表面透水时,孔压波动的幅值随深度增大先增大后减小；表面不透水时,孔压波动的幅值随深度增大而减小。表面不透水情况的最大孔压比表面透水情况高出1倍多。

（3）土体上表面位移随时间周期性变化。表面不透水情况的位移最大值明显高于表面透水情况。

（4）表面不透水且外界温度荷载相同时,土体渗透性越好,深度0 m处的超静孔压波动幅值越低,且超静孔压沿深度方向的分布更趋于均匀化,使得上表面的位移几乎不受影响。

图 1 不同幅值下疲劳荷载的应力-应变曲线

Figure 1. Stress-strain curves of fatigue loads under different amplitudes

下载: 全尺寸图片幻灯片

图 2 不同应力比下蠕变变形-时间、疲劳残余变形-荷载循环次数曲线

Figure 2. Curves of creep deformation-time and residual deformation- cycle times of fatigue loads under different stress ratios

下载: 全尺寸图片幻灯片

图 3 归一化后对数关系图

Figure 3. Graph of normalized logarithmic relation

下载: 全尺寸图片幻灯片

图 4 应力-寿命曲线图

Figure 4. Stress-life curves

下载: 全尺寸图片幻灯片

图 5 等寿命曲线图

Figure 5. Haigh diagram

下载: 全尺寸图片幻灯片

图 6 不同应力幅下的等寿命曲线图

Figure 6. Constant-life curves under different stress amplitudes

下载: 全尺寸图片幻灯片

图 7 初步等效应力曲线图

Figure 7. Curves of preliminary equivalent stress

下载: 全尺寸图片幻灯片

图 8 疲劳元件力学模型

Figure 8. Mechanical model for fatigue element

下载: 全尺寸图片幻灯片

图 9 非线性黏弹塑性疲劳模型

Figure 9. Nonlinear viscoelastic-plastic fatigue model

下载: 全尺寸图片幻灯片

图 10 频率f特殊值下的应力-时间曲线

Figure 10. Time-stress curves under special value of frequency

下载: 全尺寸图片幻灯片

图 11 不同应力条件下疲劳变形模拟

Figure 11. Fitting of fatigue deformation of stone under different stresses

下载: 全尺寸图片幻灯片

图 12 不同条件下岩石进入加速段非线性变形时完整度 $η_{0}$

Figure 12. Integrity $η_{0}$ of rock entering nonlinear deformation in acceleration section under different conditions

下载: 全尺寸图片幻灯片

表 1 各类岩样的荷载条件

Table 1 Load conditions of various types of rock samples

文献	岩石种类	单调加载强度 $σ_{c}$ /MPa	破坏门槛值 $σ_{s}$ /MPa	变形门槛值 $σ_{ve}$ /MPa	循环荷载值/MPa		频率 $f$ /Hz
文献	岩石种类	单调加载强度 $σ_{c}$ /MPa	破坏门槛值 $σ_{s}$ /MPa	变形门槛值 $σ_{ve}$ /MPa	$σ_{\min}$	$σ_{\max}$	频率 $f$ /Hz
文献[11]	花岗岩A56	112.35	89.88	3.00	28.09	106.73	0.2
文献[11]	花岗岩A47	112.35	89.88	3.00	28.09	106.73	1.0
文献[11]	花岗岩A22	112.35	89.88	3.00	22.47	101.11	1.0
文献[12]	砂岩S2	87.55	70.04	2.00	10.00	80.00	0.1
文献[12]	砂岩S2	87.55	70.04	2.00	10.00	80.00	1.0
文献[12]	砂岩S10	134.08	107.26	2.00	10.00	117.30	0.1
文献[12]	砂岩S10	134.08	107.26	2.00	10.00	117.30	1.0
文献[14]	白砂岩^#5	119.60	95.68	3.00	23.92	101.66	0.2
文献[14]	白砂岩^#3	119.60	95.68	3.00	35.88	107.64	0.2
文献[14]	白砂岩^#1	119.60	95.68	3.00	35.88	95.68	0.2

下载: 导出CSV

表 2 各参数的理论模型计算值

Table 2 Calculated values of theoretical model for each parameter

文献	岩石种类	$η$ /(MPa·s^-1)	$Ε_{Μ}$ /GPa	$Ε_{K}$ /GPa	$η_{Κ}$ /(GPa·s^-1)	$η_{M}$ /(GPa·s^-1)	$n$
本文	大理岩	6000	85	200	30000	45000	2.1
文献[11]	花岗岩	10000	50	120	5000	15000	3.2
文献[12]	砂岩S2	15000	55	120	8000	25000	3.5
文献[12]	砂岩S10	16000	68	110	2000	15000	3.5
文献[14]	白砂岩	12000	23.5	60	15000	18000	3.0

下载: 导出CSV

参考文献(23)

[1]	罗嗣海, 钱七虎, 周文斌, 等. 高放废物深地质处置及其研究概况[J]. 岩石力学与工程学报, 2004, 23(5): 831-838. doi: 10.3321/j.issn:1000-6915.2004.05.023 LUO Si-hai, QIAN Qi-hu, ZHOU Wen-bin, et al. Deep geological disposal of high-level nuclear waste and its research outlines[J]. Rock and Soil Mechanics, 2004, 23(5): 831-838. (in Chinese) doi: 10.3321/j.issn:1000-6915.2004.05.023
[2]	SHEN Y, GAO B, YANG X, et al. Seismic damage mechanism and dynamic deformation characteristic analysis of mountain tunnel after Wenchuan earthquake[J]. Engineering Geology, 2014, 180: 85-98. doi: 10.1016/j.enggeo.2014.07.017
[3]	SHI T, CHEN Z, LUO Z, et al. Mechanism of groundwater bursting in a deep rock salt mine region: a case study of the Anpeng trona and glauber mines, China[J]. Environmental Earth Sciences, 2013, 68(1): 229-239. doi: 10.1007/s12665-012-1734-8
[4]	CHEN Y, ZUO J P, GUO B H, et al. Effect of cyclic loading on mechanical and ultrasonic properties of granite from Maluanshan Tunnel[J]. Bulletin of Engineering Geology and the Environment, 2020, 79(1): 299-311. doi: 10.1007/s10064-019-01563-6
[5]	XIAO J Q, FENG X T, DING D X, et al. Investigation and modeling on fatigue damage evolution of rock as a function of logarithmic cycle[J]. International Journal of Rock Mechanics & Mining Sciences, 2011, 35: 1127-1140.
[6]	SUN B, ZHU Z, SHI C, et al. Dynamic mechanical behavior and fatigue damage evolution of sandstone under cyclic loading[J]. International Journal of Rock Mechanics and Mining Sciences, 2017, 100(94): 82-89.
[7]	WANG Y, MA L, FAN P, et al. A fatigue damage model for rock salt considering the effects of loading frequency and amplitude[J]. International Journal of Mining Science and Technology, 2016, 26(5): 955-958. doi: 10.1016/j.ijmst.2016.05.054
[8]	蔡燕燕, 唐欣, 林立华, 等. 疲劳荷载下大理岩累积损伤过程的应变速率响应[J]. 岩土工程学报, 2020, 42(5): 827-835. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC202005007.htm CAI Yan-yan, TANG Xin, LIN Li-hua, et al. Strain rate response of damage accumulation under fatigue loading on marble[J]. Chinese Journal of Geotechnical Engineering, 2020, 42(5): 827-835. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC202005007.htm
[9]	LIU X S, NING J G, TAN Y L, et al. Damage constitutive model based on energy dissipation for intact rock subjected to cyclic loading[J]. International Journal of Rock Mechanics & Mining Sciences, 2016, 85: 27-32.
[10]	LI N, CHEN W, ZHANG P, et al. The mechanical properties and a fatigue-damage model for jointed rock masses subjected to dynamic cyclical loading[J]. International Journal of Rock Mechanics & Mining Sciences, 2001, 38: 1071-1079.
[11]	MOMENI A, KARAKUS M, KHANLARI G R, et al. Effects of cyclic loading on the mechanical properties of a granite[J]. International Journal of Rock Mechanics & Mining Sciences, 2015, 77: 89-96.
[12]	LIU E L, HUANG R Q, HE S M. Effects of frequency on the dynamic properties of intact rock samples subjected to cyclic loading under conﬁning pressure conditions[J]. Rock Mechanics and Rock Engineering, 2012, 45(1): 89-102. doi: 10.1007/s00603-011-0185-y
[13]	葛修润, 蒋宇, 卢允德, 等. 周期荷载作用下岩石疲劳变形特性试验研究[J]. 岩石力学与工程学报, 2003, 22(10): 1581-1585. doi: 10.3321/j.issn:1000-6915.2003.10.001 GE Xiu-run, JIANG Yu, LU Yun-de, et al. Testing study on fatigue deformation law of rock under cyclic loading[J]. Chinese Journal of Rock Mechanics and Engineering, 2003, 22(10): 1581-1585. (in Chinese) doi: 10.3321/j.issn:1000-6915.2003.10.001
[14]	冯春林, 吴献强, 丁德馨. 周期荷载作用下白砂岩的疲劳特性研究[J]. 岩石力学与工程学报, 2009, 28(增刊1): 2749-2754. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX2009S1026.htm FENG Chun-lin, WU Xian-qiang, DING De-xin. Investigation on fatigue characteristics of white sandstone under cyclic loading[J]. Chinese Journal of Rock Mechanics and Engineering, 2009, 28(S1): 2749-2754. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX2009S1026.htm
[15]	葛修润. 岩石疲劳破坏的变形控制律、岩土力学试验的实时X射线CT扫描和边坡坝基抗滑稳定分析的新方法[J]. 岩石力学与工程学报, 2008, 30(1): 1-20. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC200801003.htm GE Xiu-run. Deformation control law of rock fatigue failure, real-time X-ray CT scan of geotechnical testing, and new method of stability analysis of slopes and dam foundations[J]. Chinese Journal of Rock Mechanics and Engineering, 2008, 30(1): 1-20. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC200801003.htm
[16]	FUENKAJORN K, PHUEAKPHUM D. Physical model simulation of shallow openings in jointed rock mass under static and cyclic loadings[J]. Engineering Geology, 2010, 113: 81-89.
[17]	KHALEDI K, MAHMOUDI E, DATCHEVA M, et al. Stability and serviceability of underground energy storage caverns in rock salt subjected to mechanical cyclic loading[J]. International Journal of Rock Mechanics & Mining Sciences, 2016, 86: 115-131.
[18]	王军保, 刘新荣, 黄明, 等. 低频循环荷载下盐岩轴向蠕变的Burgers 模型分析[J]. 岩土力学, 2014, 35(4): 934-941. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201404004.htm WANG Jun-bao, LIU Xin-rong, HUANG Ming, et al. Analysis of axial creep properties of salt rock under low frequency cyclic loading using Burgers model[J]. Rock and Soil Mechanics, 2014, 35(4): 934-941. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201404004.htm
[19]	郭建强, 黄质宏. 循环荷载作用下岩石疲劳本构模型初探[J]. 岩土工程学报, 2015, 37(9): 1698-1703. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201509025.htm GUO Jian-qiang, HUANG Zhi-hong. Constitutive model for fatigue of rock under cyclic loading[J]. Chinese Journal of Geotechnical Engineering, 2015, 37(9): 1698-1703. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201509025.htm
[20]	杨新华, 陈传尧. 疲劳与断裂[M]. 武汉: 华中科技大学出版社, 2018. YANG Xin-hua, CHEN Chuan-yao. Fatigue and Fracture[M]. Wuhan: Huazhong University of Science and Technology Press, 2018. (in Chinese)
[21]	BAGDE M N, PETROSˇA V. Fatigue properties of intact sandstone samples subjected to dynamic uniaxial cyclical loading[J]. International Journal of Rock Mechanics & Mining Sciences, 2005(42): 237-250.
[22]	薛东杰, 周宏伟, 王子辉, 等. 不同加载速率下煤岩采动力学响应及破坏机制[J]. 煤炭学报, 2016, 41(3): 595-602. https://www.cnki.com.cn/Article/CJFDTOTAL-MTXB201603011.htm XUE Dong-jie, ZHOU Hong-wei, WANG Zi-hui, et al. Failure mechanism and mining-induced mechanical properties of coal under different load-ing rates[J]. Journal of China Coal Society, 2016, 41(3): 595-602. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-MTXB201603011.htm
[23]	孙钧. 岩石流变力学及其工程应用研究的若干进展[J]. 岩石力学与工程学报, 2007, 26(6): 1081-1106. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX200706001.htm SUN Jun. Rock rheological mechanics and its advance in engineering applications[J]. Chinese Journal of Rock Mechanics and Engineering, 2007, 26(6): 1081-1106. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX200706001.htm