基于非傅里叶定律的近场动力学热力耦合模型及花岗岩热损伤破裂模拟

周露明; 朱珍德; 谢兴华; 吕茂淋

doi:10.11779/CJGE20230769

基于非傅里叶定律的近场动力学热力耦合模型及花岗岩热损伤破裂模拟 English Version

周露明^{1, 2,},
朱珍德^{1, 2},
谢兴华³,
吕茂淋^{1, 2}

1.
河海大学岩土力学与堤坝工程教育部重点实验室, 江苏南京 210029
2.
河海大学江苏省岩土工程技术工程研究中心, 江苏南京 210029
3.
南京水利科学研究院水灾害防御全国重点实验室, 江苏南京 210029

基金项目:

国家自然科学基金项目 51878249

国家自然科学基金项目 41831278

详细信息

作者简介:
周露明(1994—)，男，博士研究生，主要从事岩石力学与地下能源工程方面的研究工作。E-mail: 190204030004@hhu.edu.cn

中图分类号: TU45
计量
- 文章访问数: 283
- HTML全文浏览量: 38
- PDF下载量: 57
出版历程
- 收稿日期: 2023-08-13
- 网络出版日期: 2024-03-24
- 刊出日期: 2024-10-31

Non-Fourier law-based peridynamic thermo-mechanical coupling model and simulation of thermal damage and fracture in granite

ZHOU Luming^{1, 2,},
ZHU Zhende^{1, 2},
XIE Xinghua³,
LÜ Maolin^{1, 2}

1.
Key Laboratory of Ministry of Education for Geomechanics and Embankment Engineering, Hohai University, Nanjing 210029, China
2.
Jiangsu Research Center for Geotechnical Engineering Technology, Hohai University, Nanjing 210029, China
3.
The National Key Laboratory of Water Disaster Prevention, Nanjing Hydraulic Research Institute, Nanjing 210029, China

摘要

摘要: 研究岩石的热损伤破裂特征对于地热开采等工程具有重要意义。在传统常规态基近场动力学理论框架内，通过引入双相滞后（DPL）模型，提出了基于非傅里叶热传导定律推导得到的热力耦合模型。通过平板瞬态热传导问题及Lac du Bonnet（LdB）花岗岩的热损伤破裂试验对模型进行了验证。分析了温度梯度弛豫时间和热流弛豫时间对岩石热损伤破裂特征的影响。结果表明：近场动力学热力耦合模型模拟结果很好地反映了LdB花岗岩的热损伤破裂特征及温度分布的不连续性；温度梯度弛豫时间对热传导起促进作用，试件热损伤破裂程度随温度梯度弛豫时间的增加而增大，随热流弛豫时间的增加而减小。研究为深入理解岩石的热损伤破裂行为、优化地热能开采工程提供了有益的探索。
- 花岗岩 /
- 近场动力学 /
- 热损伤破裂 /
- 热力耦合 /
- 热弛豫时间
Abstract: It is of great significance to study the thermal damage and fracture characteristics of rocks for deep rock projects, such as geothermal exploitation. Within the framework of the traditional classical ordinary state-based peridynamic theory, a thermo-mechanical coupling model based on the non-Fourier heat conduction law is proposed by introducing a dual-phase-lag model. The model is validated through the transient heat conduction problems in the plate and thermal damage and fracture tests on LdB granite. It is found that the simulated results accurately reflect the thermal damage and fracture characteristics, as well as the discontinuity in temperature distribution, of LdB granite. Furthermore, through numerical analysis, it is observed that the temperature gradient relaxation time promotes heat conduction, leading to an increase in the degree of the thermal damage and fracture with an increase in the temperature gradient relaxation time, while it decreases with an increase in the heat flux relaxation time. This research provides valuable insights for a better understanding of the thermal damage and fracture behavior of rocks and offers beneficial explorations for optimizing the geothermal energy extraction projects.
- granite /
- peridynamics /
- thermal damage and fracture /
- thermo-mechanical coupling /
- thermal relaxation time

HTML全文

0. 引言

混凝土面板坝是以堆石料或砂砾石料为主要筑坝材料，以钢筋混凝土面板及其分缝止水作为防渗结构的高土石坝坝型^[1]。2000年以来，国内外面板坝进入200 m级特高坝时代，巴西于2006年建成最大坝高202 m的Campos Novos面板坝^[2]；老挝于2005年建成Nam Ngum 3水电站，其面板坝高度达到220 m^[3]；马来西亚于2009年建成Bakun面板坝，最大坝高205 m^[4]；中国则相继建成了天生桥一级、洪家渡、三板溪、水布垭、江坪河、猴子岩等一系列200 m级高混凝土面板坝^[5]。

混凝土面板是面板坝防渗体系的关键防线，其应力变形特性是影响大坝安全的关键。作为弹性薄板结构，面板上作用的法向水压力使其与垫层料之间产生相对位移，从而产生沿接触面的切向摩擦力。显然，面板应力变形主要受控于支撑其工作的堆石料在水压力作用下的应力变形特性及面板与堆石料之间的摩擦接触特性。目前，堆石料的本构模型研究取得大量成果，已能够比较全面地考虑其剪胀剪缩^[6-7]、颗粒破碎^[8]、湿化流变^[9]等复杂特性。

接触面模型研究同样长期受到关注，20世纪60年代，Goodman等^[10]就针对裂隙岩体问题提出无厚度接触面单元，成为解决复杂接触问题的基本方法；70年代初期，Wayne等^[11]提出了模拟土体与混凝土接触面剪应力与剪切位移关系的双曲线模型，并广泛应用于岩土工程^[12]。近年来，接触面模型研究主要集中在模拟单调和循环荷载作用下接触面法向胀缩、峰后软化等复杂行为^[13-17]，这类模型的本质是模拟接触面附近土体强剪切带的应力变形特性^[18]。

值得指出的是，上述接触面模型大多基于平面应变剪切试验建立，即假定剪应力和相对剪切位移是沿固定方向。实际面板坝工程大多修建于V形河谷，三维效应比较显著，面板与垫层之间接触应力和相对位移量值与方向的变化远较平面应变条件复杂。Evgin等^[19]曾开展接触面双向剪切试验，发现两个剪切方向耦合效应显著。本文简要介绍了双曲线接触面模型及其三维形式；提出了三维条件下接触面模型应满足的条件；提出了一个考虑双向剪切特性的接触面模型，并以水布垭面板坝为例，展示了接触面模型在面板坝应力变形计算中的应用。

1. 双曲线接触面剪切模型

Wayne等^[11]假定恒定法向应力σ_n条件下，接触面相对剪切位移u和平均剪应力τ之间符合双曲线关系，如图 1所示，即

。

$\tau = \frac{u}{{(1/{G_0}) + (u/{\tau _\infty })}} {\rm{。}}$

(1)

图 1 堆石料与混凝土接触面试验结果

Figure 1. Test results of rockfill-concrete interface

下载: 全尺寸图片幻灯片

式中： ${G_0}$ ， ${\tau _\infty }$ 分别为初始剪切模量和剪应力渐进值（图 1），

。

$\left. \begin{array}{l}{G}_{0}=\underset{u\to 0}{\mathrm{lim}}\left(\frac{\tau }{u}\right)\text{ }\text{, }\\ {\tau }_{\infty }=\underset{u\to \infty }{\mathrm{lim}}\tau \text{ }。\end{array} \right\}$

(2)

式（1）两侧对剪切位移u求导可得接触面切线剪切模量G_t，即

${G_{\text{t}}} = \frac{{{\text{d}}\tau }}{{{\text{d}}u}} = {\left( {1 - \frac{\tau }{{{\tau _\infty }}}} \right)^2}{G_0} = {\left( {1 - {R_{\text{f}}}\frac{\tau }{{{\tau _{\text{f}}}}}} \right)^2}{G_0} {\rm{。}}$

(3)

式中：R_f为破坏比，R_f = τ_f/τ_∞；初始剪切模量G₀与接触面法向应力σ_n有关，

${G_0} = k \cdot {\gamma _{\text{w}}} \cdot {\left( {\frac{{{\sigma _{\text{n}}}}}{{{p_{\text{a}}}}}} \right)^n} {\rm{。}}$

(4)

式中： ${\gamma _{\text{w}}}$ 为水的重度；p_a为大气压力；k，n为参数。

式（3）中， ${\tau _{\text{f}}}$ 为给定法向应力时接触面的破坏剪应力，一般采用Coulomb破坏准则，即

${\tau _{\text{f}}} = c + {\sigma _{\text{n}}} \cdot \tan \delta {\rm{。}}$

(5)

式中：c，δ分别为接触面的黏聚力和内摩擦角。

式（3）~（5）给出了单一方向剪切时接触面剪切模量表达式，对于双向剪切问题，需对其进行扩展。目前，最常用的扩展表达式为^{[12, 20]}

$\left. \begin{array}{l}{G}_{zx}={\left(1-{R}_{\text{f}}\frac{\left|{\tau }_{zx}\right|}{c+{\sigma }_{z}\mathrm{tan}\delta }\right)}^{2}\cdot k\cdot {\gamma }_{\text{w}}\cdot {\left(\frac{{\sigma }_{z}}{{p}_{\text{a}}}\right)}^{n}\text{, }\\ {G}_{zy}={\left(1-{R}_{\text{f}}\frac{\left|{\tau }_{zy}\right|}{c+{\sigma }_{z}\mathrm{tan}\delta }\right)}^{2}\cdot k\cdot {\gamma }_{\text{w}}\cdot {\left(\frac{{\sigma }_{z}}{{p}_{\text{a}}}\right)}^{n}{\rm{。}}\end{array} \right\}$

(6)

式中：τ_zx，τ_zy为接触面两个正交方向的剪应力；σ_z为接触面法向应力；x，y分别为接触面的两个切向坐标，z为接触面法向坐标。这样，接触面的应力-相对位移关系式可以表示为

$\left( \begin{array}{l} {\text{d}}{\tau _{zx}} \hfill \\ {\text{d}}{\tau _{zy}} \hfill \\ {\text{d}}{\sigma _z} \hfill \end{array} \right) = \left[ {\begin{array}{*{20}{c}} {{G_{zx}}}&0&0 \\ 0&{{G_{zy}}}&0 \\ 0&0&{{K_{zz}}} \end{array}} \right] \cdot \left( \begin{array}{l} {\text{d}}{u_x} \hfill \\ {\text{d}}{u_y} \hfill \\ {\text{d}}{u_z} \hfill \end{array} \right) {\rm{。}}$

(7)

式中：u_x，u_y，u_z分别为接触面的两个切向相对位移和法向相对位移；G_zx，G_zy为接触面的剪切模量；K_zz为法向接触模量，接触面处于压紧状态时取一高值，使接触面两侧材料近似满足不可嵌入条件；接触面处于拉开状态时取一低值，使接触面两侧材料可以相互独立变形^{[12, 20]}。

式（6）是式（3）的实用性推广，当沿x，y方向的剪应力达到破坏标准时，剪切模量趋近于零；当其中任意方向剪切位移为零时，另一方向剪应力-相对位移关系退化至式（3）。从式（7）可以看出，上述三维接触面模型的两个切线方向是相互独立的，不存在相互影响；接触面法向对切向的影响则通过剪切模量中法向应力体现。

2. 三维接触面剪切模型基本要求

接触面单元的应力位移关系一般在正交局部坐标系中建立，现设整体坐标系为XYZ，接触面单元局部坐标系为xyz，如图 2所示。下面分析说明接触面单元模型应满足的基本要求。

图 2 接触面单元分析的整体和局部坐标系

Figure 2. Global and local coordinate system for interface element analysis

下载: 全尺寸图片幻灯片

2.1 破坏准则

接触面双向剪切试验结果表明^[19]，其峰值强度与应力路径无关，且在（τ_zx，τ_zy）平面上，破坏包线近似为圆形，其半径与法向应力σ_z成正比。因此，在三维应力空间中，接触面单元的破坏面是一个圆锥面，锥面顶角为2δ，如图 3所示。这一破坏准则要求，当应力状态点位于破坏面并进一步增大剪切位移时，其各向剪切模量应趋近于零，即

${G_{zx}}\left( {\frac{{\sqrt {{\tau _{zx}}^2 + {\tau _{zy}}^2} }}{{c + {\sigma _{\text{n}}}\tan \delta }} = 1} \right) = {G_{zy}}\left( {\frac{{\sqrt {{\tau _{zx}}^2 + {\tau _{zy}}^2} }}{{c + {\sigma _{\text{n}}}\tan \delta }} = 1} \right) \to {0_{}} {\rm{。}}$

(8)

图 3 应力空间中接触面单元破坏面

Figure 3. Failure surface of interface in stress space

下载: 全尺寸图片幻灯片

显然，式（6）并不满足式（8）的条件。在式（6）中，当且仅当其各向剪应力分别达到破坏标准时，该方向的剪切模量趋近于零。因此，在应力空间中，其破坏面是顶面为正方形的锥面，故式（6）高估了接触面的抗剪强度。

2.2 剪切耦合效应

接触面受双向剪切荷载作用时，两个剪切方向的耦合效应已由试验证实^[19]，这种剪切耦合效应可通过简单的力学分析加以说明，如图 4所示。设想一接触面受恒定法向应力作用，现对该接触面施加x方向剪应力至应力状态A，该过程不产生y向相对位移；然后保持x向剪应力不变，持续增加y向剪应力，直至应力状态达到破坏面上B点。达到B点后，y向剪应力将不再继续增加，但相对位移将持续增加。显然，应力状态B点处产生的塑性位移不可能仅有y向分量，必然还存在x向分量，接触面塑性滑移方向应与剪应力合力方向一致。y向剪应力施加过程是一个连续过程，接触面的位移响应是连续的，因此不可能在B点处突然出现x向位移增量。也就是说，即使保持x向剪应力不变，连续施加y向剪应力的过程中，x向位移增量应贯穿整个y向加载过程（应力路径AB）。如前所述，式（6），（7）给出的传统模型无法反映上述剪切耦合效应。

图 4 接触面双向剪切耦合效应示意图

Figure 4. Coupling effects in two-way shearing

下载: 全尺寸图片幻灯片

2.3 标架无差异原则

连续介质力学中，有所谓的标架无差异原则，即材料的应力应变关系是其固有属性，不应受到所选取坐标系的影响^[21]。接触面单元同样应符合标架无差异原则，即接触面的剪应力-相对位移关系与所选取的两个切向坐标的方向无关。为说明该原则对接触面剪切模量矩阵的要求，考虑两个不同的切向坐标系xy和xʹyʹ，两坐标系具有相同的法向z，xʹyʹ坐标系是由xy坐标系逆时针旋转θ角度后得到的，如图 5所示。

图 5 接触面切向坐标系旋转示意图

Figure 5. Rotation of tangential directions

下载: 全尺寸图片幻灯片

现设xy坐标系中，接触面应力位移关系可表示为

$\left. \begin{array}{l}\text{d}\mathit{\boldsymbol{\tau}} =\mathit{\boldsymbol{G}}\cdot \text{d}\mathit{\boldsymbol{u}}\text{ }\text{, }\\ \left(\begin{array}{l}\text{d}{\tau }_{zx}\\ \text{d}{\tau }_{zy}\end{array}\right)=\left[\begin{array}{cc}{G}_{xx}& {G}_{xy}\\ {G}_{yx}& {G}_{yy}\end{array}\right]\left(\begin{array}{l}\text{d}{u}_{x}\\ \text{d}{u}_{y}\end{array}\right)\text{ }{\rm{。}}\end{array} \right\}$

(9)

类似地，xʹyʹ坐标系中，接触面的应力位移关系可以表示为

$\left. \begin{array}{l}\text{d}{\mathit{\boldsymbol{\tau}} }^{\prime }={\mathit{\boldsymbol{G}}}^{\prime }\cdot \text{d}{\mathit{\boldsymbol{u}}}^{\prime }\text{ }\text{, }\\ \left(\begin{array}{l}\text{d}{{\tau }^{\prime }}_{zx}\\ \text{d}{{\tau }^{\prime }}_{zy}\end{array}\right)=\left[\begin{array}{cc}{{G}^{\prime }}_{xx}& {{G}^{\prime }}_{xy}\\ {{G}^{\prime }}_{yx}& {{G}^{\prime }}_{yy}\end{array}\right]\left(\begin{array}{l}\text{d}{{u}^{\prime }}_{x}\\ \text{d}{{u}^{\prime }}_{y}\end{array}\right)\text{ }{\rm{。}}\end{array} \right\}$

(10)

式（9），（10）中: G，Gʹ分别为两坐标系中的剪切模量矩阵。若记

$\mathit{\boldsymbol{R}} = \left[ {\begin{array}{*{20}{c}} {\cos \theta }&{\sin \theta } \\ { - \sin \theta }&{\cos \theta } \end{array}} \right] \text{, }$

(11)

则

$\left. \begin{array}{l}\text{d}{\mathit{\boldsymbol{\tau}} }^{\prime }=\mathit{\boldsymbol{R}}\cdot \text{d}\mathit{\boldsymbol{\tau}} \text{ }\text{, }\\ \text{d}{\mathit{\boldsymbol{u}}}^{\prime }=\mathit{\boldsymbol{R}}\cdot \text{d}\mathit{\boldsymbol{u}}\text{ }{\rm{。}}\end{array} \right\}$

(12)

将式（12）代入式（10），并将其与式（9）对比可知，两局部坐标系中剪切模量矩阵应满足下面的关系：

$\left. \begin{array}{l}\mathit{\boldsymbol{G}}={\mathit{\boldsymbol{R}}}^{\text{T}}\cdot {\mathit{\boldsymbol{G}}}^{\prime }\cdot \mathit{\boldsymbol{R}}\text{ }\text{, }\\ {\mathit{\boldsymbol{G}}}^{\prime }=\mathit{\boldsymbol{R}}\cdot \mathit{\boldsymbol{G}}\cdot {\mathit{\boldsymbol{R}}}^{\text{T}}\text{ }{\rm{。}}\end{array} \right\}$

(13)

显然式（6），（7）给出的剪切模量矩阵一般情况下并不符合式（13）中的转换关系，即其应力位移关系受到局部坐标系选择的影响。

为更清楚地说明上述问题，设一接触面受到恒定法向应力σ_z = 500 kPa作用，其后沿x向施加剪应力至τ_zx = 400 kPa。现采用式（6）分别在xy坐标系和xʹyʹ坐标系中研究其位移路径，相关参数如表 1所示^[20]。

表 1 堆石料-混凝土接触面计算参数

Table 1. Parameters of rockfill-concrete interface

R_f	k	n	c/ kPa	δ/(°)
0.90	6000	0.85	0.0	41.5

下载: 导出CSV

| 显示表格

图 6中给出了θ为0°（xy坐标系），30°，60°情况下接触面的应力位移关系和应力路径，在xy坐标系中考察该问题时，接触面沿x向位移达到9.37 mm；但在θ为30°，60°的局部坐标系中计算时，x向位移仅5.18 mm，可见式（6）高估了接触面剪切模量。更不合理之处在于，该案例仅沿x向施加荷载，但在θ为30°，60°的局部坐标系中计算时，会产生y向剪切位移。该案例清楚地表明，式（6），（7）预测的应力位移关系和位移路径与局部坐标系选择有关，违背了标架无差异原则。

图 6 不同局部坐标系下接触面应力位移关系与位移路径

Figure 6. Stress-displacement relations and displacement paths under different local coordinate systems

下载: 全尺寸图片幻灯片

3. 接触面的双向剪切弹塑性模型

根据前文分析，面板坝计算常用的接触面模型存在高估剪切强度、不能考虑剪切耦合效应以及违反标架无差异原则等不足，势必影响面板应力计算结果的可靠性。为解决上述问题，本文在弹塑性理论框架下，构建一个新的接触面模型。

3.1 基本假定

接触面剪切试验所揭示的剪胀剪缩特性本质上是其附近强剪切带内土体的力学行为，原则上可通过剪切带土体的本构模型模拟。因此，本文仅限于模拟所谓无厚度接触面单元在压紧状态时的剪切特性，并采用2点基本假定。

（1）不考虑接触面的剪胀与剪缩，接触面的法向和切向耦合效应仅通过法向应力对剪切模量的影响体现。

（2）接触面的剪切位移增量可分解为弹性位移增量和塑性位移增量，即

$\text{d}\mathit{\boldsymbol{u}}=\text{d}{\mathit{\boldsymbol{u}}}^{\text{e}}+\text{d}{\mathit{\boldsymbol{u}}}^{\text{p}}\text{, }\left(\begin{array}{l}\text{d}{u}_{x}\\ \text{d}{u}_{y}\end{array}\right)=\left(\begin{array}{l}\text{d}{u}_{x}^{\text{e}}\\ \text{d}{u}_{y}^{\text{e}}\end{array}\right)+\left(\begin{array}{l}\text{d}{u}_{x}^{\text{p}}\\ \text{d}{u}_{y}^{\text{p}}\end{array}\right) {\rm{。}}$

(14)

式中：上标e，p分别表示弹性和塑性。下面分别说明弹性和塑性剪切位移的计算方法。

3.2 弹性位移计算

接触面的弹性位移增量采用下式计算：

$\left( \begin{array}{l} {\text{d}}u_x^e \hfill \\ {\text{d}}u_y^{\text{e}} \hfill \end{array} \right) = \left[ {\begin{array}{*{20}{c}} {\frac{1}{{{G_{\text{e}}}}}}&0 \\ 0&{\frac{1}{{{G_{\text{e}}}}}} \end{array}} \right]\left( \begin{array}{l} {\text{d}}{\tau _{zx}} \hfill \\ {\text{d}}{\tau _{zy}} \hfill \end{array} \right) {\rm{。}}$

(15)

式中：G_e为接触面弹性剪切模量，

${G_{\text{e}}} = {k_{\text{e}}} \cdot {\gamma _{\text{w}}} \cdot {\left( {\frac{{{\sigma _{\text{n}}}}}{{{p_{\text{a}}}}}} \right)^n} {\rm{。}}$

(16)

式中：k_e为接触面弹性参数，可取k_e = (2~3)k；n为模量参数，与式（4）中的参数n一致。

3.3 塑性位移计算

接触面处于滑动状态时，其滑动方向应与剪应力合力方向一致。因此，本文假定塑性位移增量方向始终与剪应力方向一致，即

$\left( \begin{array}{l} {\text{d}}u_x^{\text{p}} \hfill \\ {\text{d}}u_y^{\text{p}} \hfill \end{array} \right) = \left( \begin{array}{l} \frac{{{\tau _{zx}}}}{{\sqrt {{\tau _{zx}}^2 + {\tau _{zy}}^2} }} \hfill \\ \frac{{{\tau _{{\text{zy}}}}}}{{\sqrt {{\tau _{zx}}^2 + {\tau _{zy}}^2} }} \hfill \end{array} \right){\text{d}}{u^{\text{p}}} = \left( \begin{array}{l} \frac{{{\tau _{zx}}}}{{\sqrt {{\tau _{zx}}^2 + {\tau _{zy}}^2} }} \hfill \\ \frac{{{\tau _{zy}}}}{{\sqrt {{\tau _{zx}}^2 + {\tau _{zy}}^2} }} \hfill \end{array} \right)\frac{1}{{{G_{\text{p}}}}}{\text{d}}\tau {\rm{。}}$

(17)

式中：du^p为塑性位移增量的合成量；G_p可称为接触面的塑性剪切模量；dτ为剪应力增量的合成量，即

$\begin{array}{l}\text{d}\tau =\text{d}\sqrt{{\tau }_{zx}{}^{2}+{\tau }_{zy}{}^{2}}\\ \ \ \ \ =\left(\frac{{\tau }_{zx}}{\sqrt{{\tau }_{zx}{}^{2}+{\tau }_{zy}{}^{2}}}\text{, }\frac{{\tau }_{zy}}{\sqrt{{\tau }_{zx}{}^{2}+{\tau }_{zy}{}^{2}}}\right)\left(\begin{array}{l}\text{d}{\tau }_{zx}\\ \text{d}{\tau }_{zy}\end{array}\right)\end{array} {\rm{。}}$

(18)

因此式（17）可以改写为

$\left( \begin{array}{l} {\text{d}}u_x^{\text{p}} \hfill \\ {\text{d}}u_y^{\text{p}} \hfill \end{array} \right) = \frac{1}{{{G_{\text{p}}}}}\left[ {\begin{array}{*{20}{c}} {\frac{{{\tau _{zx}}^2}}{{{\tau _{zx}}^2 + {\tau _{zy}}^2}}}&{\frac{{{\tau _{zx}}{\tau _{zy}}}}{{{\tau _{zx}}^2 + {\tau _{zy}}^2}}} \\ {\frac{{{\tau _{zy}}{\tau _{zx}}}}{{{\tau _{zx}}^2 + {\tau _{zy}}^2}}}&{\frac{{{\tau _{zy}}^2}}{{{\tau _{zx}}^2 + {\tau _{zy}}^2}}} \end{array}} \right]\left( \begin{array}{l} {\text{d}}{\tau _{zx}} \hfill \\ {\text{d}}{\tau _{zy}} \hfill \end{array} \right) {\rm{。}}$

(19)

3.4 本构方程

将式（15），（19）代入式（14）可得柔度形式的接触面本构方程，即

$\left( \begin{array}{l} {\text{d}}{u_x} \hfill \\ {\text{d}}{u_y} \hfill \end{array} \right) = \left[ {\begin{array}{*{20}{c}} {\frac{1}{{{G_{\text{e}}}}} + \frac{1}{{{G_{\text{p}}}}}\frac{{{\tau _{zx}}^2}}{{{\tau _{zx}}^2 + {\tau _{zy}}^2}}}&{\frac{1}{{{G_{\text{p}}}}}\frac{{{\tau _{zx}}{\tau _{zy}}}}{{{\tau _{zx}}^2 + {\tau _{zy}}^2}}} \\ {\frac{1}{{{G_{\text{p}}}}}\frac{{{\tau _{{\text{zy}}}}{\tau _{{\text{zx}}}}}}{{{\tau _{zx}}^2 + {\tau _{zy}}^2}}}&{\frac{1}{{{G_{\text{e}}}}} + \frac{1}{{{G_{\text{p}}}}}\frac{{{\tau _{zy}}^2}}{{{\tau _{zx}}^2 + {\tau _{zy}}^2}}} \end{array}} \right]\left( \begin{array}{l} {\text{d}}{\tau _{zx}} \hfill \\ {\text{d}}{\tau _{zy}} \hfill \end{array} \right) {\rm{。}}$

(20)

由式（20）可进一步写出劲度形式的接触面本构方程，这里直接给出式（9）中剪切模量矩阵的表达式为

${\mathbf{G}} = \left[ {\begin{array}{*{20}{c}} {{G_{\text{e}}}}&0 \\ 0&{{G_{\text{e}}}} \end{array}} \right] - \frac{{{G_{\text{e}}}^2}}{{{G_{\text{e}}} + {G_{\text{p}}}}}\left[ {\begin{array}{*{20}{c}} {\frac{{\tau _{zx}^2}}{{\tau _{zx}^2 + \tau _{zy}^2}}}&{\frac{{{\tau _{zx}} \cdot {\tau _{zy}}}}{{\tau _{zx}^2 + \tau _{zy}^2}}} \\ {\frac{{{\tau _{zy}} \cdot {\tau _{{\text{zx}}}}}}{{\tau _{zx}^2 + \tau _{zy}^2}}}&{\frac{{\tau _{zy}^2}}{{\tau _{zx}^2 + \tau _{zy}^2}}} \end{array}} \right] {\rm{。}}$

(21)

式中：右侧第一项满足式（13）；第二项方括号内可表示为τ·τ^T/|τ|²，应用式（12）的转换关系可以证实其满足式（13），故上述模型满足标架无差异原则。从推导过程可以看出，本文模型满足该原则的关键是假定弹性剪切模量各向同性，且塑性位移增量方向与剪应力方向一致。

3.5 塑性模量

完成本构模型构建工作，尚需给出塑性剪切模量表达式。为此，考虑沿一固定方向施加剪切荷载的情形，即

$\frac{{{\text{d}}{u_x}}}{{{\text{d}}{u_y}}} = \frac{{{\tau _{zx}}}}{{{\tau _{zy}}}} = \frac{{{\text{d}}{\tau _{zx}}}}{{{\text{d}}{\tau _{zy}}}} = r {\rm{。}}$

(22)

式中：r为比例系数。

将式（22）代入式（20）可得

${\text{d}}u\left( \begin{array}{l} r \hfill \\ 1 \hfill \end{array} \right) = \frac{{{\text{d}}\tau }}{{{G_{\text{e}}}}}\left( \begin{array}{l} r \hfill \\ 1 \hfill \end{array} \right) + \frac{{{\text{d}}\tau }}{{{G_{\text{p}}}}}\left( \begin{array}{l} r \hfill \\ 1 \hfill \end{array} \right) {\rm{。}}$

(23)

即

$\frac{1}{{{G_{\text{t}}}}} = \frac{{{\text{d}}u}}{{{\text{d}}\tau }} = \frac{1}{{{G_{\text{e}}}}} + \frac{1}{{{G_{\text{p}}}}} {\rm{。}}$

(24)

因此

${G_{\text{p}}} = {\left( {\frac{1}{{{G_{\text{t}}}}} - \frac{1}{{{G_{\text{e}}}}}} \right)^{ - 1}} = \frac{{{G_{\text{e}}} \cdot {G_{\text{t}}}}}{{{G_{\text{e}}} - {G_{\text{t}}}}} {\rm{。}}$

(25)

注意式（25）中切线剪切模量仍采用式（3）的基本形式，但将其扩展为二维剪切情形，即

${G_{\text{t}}} = {\left( {1 - {R_{\text{f}}}\frac{{\sqrt {{\tau _{zx}}^2 + {\tau _{zy}}^2} }}{{c + {\sigma _z}\tan \delta }}} \right)^2} \cdot k \cdot {\gamma _{\text{w}}} \cdot {\left( {\frac{{{\sigma _z}}}{{{p_{\text{a}}}}}} \right)^n} {\rm{。}}$

(26)

显然，无论剪应力分量如何，只要当其合成量满足破坏准则时，切线剪切模量趋近于零（R_f = 1时），因此塑性模量也趋近于零，微小的应力增量即可引起很大的塑性位移。

3.6 加卸载准则

弹塑性模型的加卸载准则一般由屈服函数规定，但本文未定义屈服函数，故需补充加卸载准则。记接触面最大历史剪应力为τ_max，规定若接触面应力和位移增量同时满足下述条件

$\left. \begin{array}{l}\sqrt{{\tau }_{zx}^{2}+{\tau }_{zy}^{2}}\ge {\tau }_{\text{max}}\text{ }\text{, }\\ {\tau }_{zx}\text{d}{u}_{x}+{\tau }_{zy}\text{d}{u}_{y}\ge 0\text{ }{\rm{。}}\end{array} \right\}$

(27)

则接触面处于加载状态，剪切模量矩阵由式（21）给出；否则，接触面处于卸载状态，剪切模量矩阵仅含式（21）右端第一项。式（27）实际上规定：在剪切应力（τ_zx，τ_zy）平面上，若当前剪切应力状态点位于历史最大剪应力包线上，且剪切位移增量方向与剪应力合力方向夹锐角，则接触面处于加载状态；否则，按卸载处理。

3.7 模型验证

采用弹塑性理论，将传统平面应变接触面模型扩展至双向剪切情形，使其满足三维破坏准则和标架无差异原则，并能够考虑两个剪切方向的相互影响。除增加弹性剪切模量参数k_e外，未增加其它参数，且所有参数仍可通过恒定法向应力的单向剪切试验确定；k_e可按经验取(2~3)k。这里，首先采用Evgin等^[19]的双向剪切试验对该模型进行验证。

Evgin等针对密砂与粗糙钢板面开展了恒定法向应力的双向剪切试验，试验时首先施加法向应力至 ${\sigma _z}$ = 100 kPa；然后分别施加τ_zy至0，20，40，60 kPa；其后保持 ${\tau _{zy}}$ 恒定，以位移控制方式施加 ${\tau _{zx}}$ 直至接触面发生剪切破坏。图 7中给出了试验的应力路径和x向剪应力和位移关系，可以看出密砂与粗糙钢板接触面具有明显的峰后软化效应，但本文切线模量表达式（26）无法反映软化现象，故此处模拟接触面的峰值强度特性，表 2列出了根据试验结果确定的模型参数，图 7中曲线是应用表中参数模拟试验的结果，可以看出本文模型较好地反映了沿接触面一个切向施加剪应力后，另一切向剪切强度的降低。

图 7 接触面双向剪切试验结果^[19]

Figure 7. Two-way shear results of interface^[19]

下载: 全尺寸图片幻灯片

表 2 砂-钢板接触面计算参数

Table 2. Parameters of sand-steel interface

R_f	k/10³	k_e/10³	n	c/kPa	δ/(°)
1.0	150	300	1.0	0.0	38.7

下载: 导出CSV

| 显示表格

图 8对比了上述试验和模型计算得到的位移路径，其中不同形状点是试验结果；曲线是模拟结果。尽管模型计算的应力位移曲线与试验结果存在差异，但模型计算的位移路径与试验结果是相当相符的。由于接触面施加x向剪切位移后，很快达到峰值强度（图 7），其后进入软化阶段，故无论是x向剪切位移，还是y向剪切位移，均以塑性滑移占主导。因此，图 8中模拟结果与试验结果相符实际上表明，本文采用的塑性流动准则式（17）是可以接受的。

图 8 接触面双向剪切试验的位移路径^[19]

Figure 8. Displacement paths in two-way shear^[19]

下载: 全尺寸图片幻灯片

为进一步对比本文模型与传统模型（即式（6），（7））的差异，这里运用表 1中所列参数研究两个模型模拟等位移比剪切时的应力路径。假定接触面法向应力为1 MPa，沿7个不同的位移路径施加剪切位移至10 mm，如图 9（a）所示；运用两种模型计算的接触面剪应力路径绘制于图 9（b）中。本文模型计算结果显示，x向和y向剪切位移等比例增加时，两个方向剪应力τ_zx和τ_zy也等比例增加，应力路径与位移路径方向一致，最终应力状态点落入一圆弧上，即最终剪应力合力相等；传统模型因其两个方向剪切模量的不同步变化（式6），计算的剪应力路径明显偏离等应力比路径，最终应力状态点甚至落入强度包线之外，这显然是不合理的。

图 9 等位移比剪切应力路径

Figure 9. Stress paths in constant direction shear

下载: 全尺寸图片幻灯片

4. 接触面模型的运用

前文所述接触面模型仅涉及两个剪切方向，实际运用时仍需补充接触面法向本构关系。最简单的处理方式是按照Goodman等^[10]的方法，根据接触面压紧或张开状态，分别对法向模量赋不同量级的值。本节以水布垭面板坝为例，将本文模型运用于面板位移与应力计算。

4.1 工程概况

水布垭面板坝最大坝高233 m，是中国已建最高面板坝^[22]。大坝从上游至下游主要有5个材料分区，即垫层区（2A）、过渡区（3A）、主堆石区（3B）、次堆石区（3C）以及下游堆石区（3D），如图 10所示。实际施工时，大坝分4期填筑；面板分3期浇筑。计算模拟时对上述施工过程进行适当简化：①堆石料按图 10所示4个分期填筑；②从底至顶一次浇筑面板至405 m高程；③水库水位由176 m逐渐提升至400 m高程；④水库在400 m水位条件下运行5 a。

图 10 大坝材料分区与填筑过程

Figure 10. Dam material zones and filling processes

下载: 全尺寸图片幻灯片

4.2 模型与参数

根据坝体施工过程和面板分缝情况建立了大坝三维有限元模型，如图 11所示。面板和垫层料之间设置接触面单元；面板和趾板之间设置周边缝单元；面板和面板之间设置垂直缝单元。接触面单元参数如表 1所列，并取k_e = 2k；缝单元采用文献[12]中的模型与参数。模型的左右岸截断边界设置坝轴向（x向）位移约束；坝体与建基面接触的结点设置为三向约束。

图 11 大坝的三维有限元网格

Figure 11. Finite element mesh of dam

下载: 全尺寸图片幻灯片

堆石料采用“南水”双屈服面弹塑性本构模型^[7]，并考虑流变影响^[23]，其中垫层料、过渡料、主堆石料、下游堆石料采用相同参数；次堆石料采用另一组参数，如表 3，4所列。

表 3 筑坝堆石料“南水”模型参数^[22]

Table 3. 'Nanshui' constitutive model parameters of damming rockfill materials^[22]

材料	ρ/(g·cm^-3)	φ₀/(°)	Δφ/(°)	R_f	k	n	c_d/ %	n_d	R_d
2A, 3A, 3B, 3D	2.20	54.7	10.4	0.81	994	0.33	0.29	0.84	0.72
3C	2.20	51.3	10.4	0.83	602	0.25	0.28	0.98	0.75

下载: 导出CSV

| 显示表格

表 4 筑坝堆石料流变模型参数^[23]

Table 4. Creep model parameters of damming rockfill materials^[23]

材料	t₀/月	c₁/ %	n₁	c₂/ %	n₂	c₃/ %	n₃	n₄
2A, 3A, 3B, 3D	5.5	0.160	0.45	0.006	0.88	0.032	0.66	0.50
3C	5.5	0.180	0.43	0.007	0.90	0.035	0.55	0.50

下载: 导出CSV

| 显示表格

对于混凝土趾板和面板，采用各向同性线弹性模型模拟，弹性模量为E = 28 GPa；弹性泊松比为 $\nu$ = 0.167。

4.3 计算结果分析

图 12，13中分别绘制了用本文接触面模型计算的蓄水期和运行期面板轴向位移、法向位移（挠度）、轴向应力和顺坡向应力等值线分布。水库蓄水至400 m高程时，两岸面板均向河谷中央位移，因此在河床中部面板中产生轴向挤压应力，最大值约12.71 MPa，两岸坡部位面板则出现轴向拉应力，最大值约2.20 MPa；水压力作用下面板挠度从坝顶至坝底先增大后减小，最大挠度发生于1/3坝高处；面板顺坡向几乎处于全断面受压状态，最大压应力为14.84 MPa。

图 12 本文模型计算的蓄水后面板位移和应力分布等值线

Figure 12. Displacements and stresses of concrete slabs after impounding predicted with proposed model

下载: 全尺寸图片幻灯片

图 13 本文模型计算的运行期面板位移和应力分布等值线

Figure 13. Displacements and stresses of concrete slabs during operation predicted by proposed model

下载: 全尺寸图片幻灯片

运行期因坝体堆石料产生流变，两岸坡部位面板进一步向河床中央挤压，中部面板轴向压应力显著增大至35.15 MPa，峰值区上移至320 m高程以上；因堆石体沉降，面板挠度峰值区上移至坝顶部位，最大值达到84.5 cm，与实测结果基本接近^[24]，与此同时，面板顺坡向压应力也大幅增长至28.00 MPa。在该案例中，堆石流变产生的面板坝轴向压应力增幅超过顺坡向压应力增幅。该坝于2008年11月首次蓄水至正常蓄水位附近，运行过程中河床中部面板出现垂直缝两侧面板表层混凝土挤压破损，部位与本文所计算的压应力极值区基本吻合^[24]，表明本文接触面模型较好地反映了堆石坝体与面板之间的荷载传递机制。

5. 结论

分析了现有三维接触面模型的不足，提出了一个模拟混凝土面板与堆石料双向摩擦接触特性的无厚度接触面模型，该模型与传统模型及同类弹塑性模型相比，具有3个特点。

（1）满足接触面剪切强度准则和力学响应的标架无差异原则，且可以反映两个剪切方向的相互影响。满足标架无差异原则，使该模型实际应用时，计算结果与局部坐标系选择无关。

（2）由于未考虑剪切软化等复杂特性，模型参数仍比较简单，与传统接触面模型相比，仅增加一个弹性参数，只需通过常规的恒定法向应力的接触面剪切试验即可确定全部参数。

（3）模型的数值实现方便，只需简单修改传统接触面模型的剪切模量矩阵即可，且模型模拟结果与实际案例出现的情况定性上相符。

图 1 近场动力学相互作用的物质点

Figure 1. Material points of peridynamic interactions

下载: 全尺寸图片幻灯片

图 2 近场动力学变形示意图

Figure 2. Schematic diagram of peridynamic deformation

下载: 全尺寸图片幻灯片

图 3 PD模型离散示意图

Figure 3. Diagrammatic sketch of PD model discretization

下载: 全尺寸图片幻灯片

图 4 均质方板的几何模型和边界条件

Figure 4. Geometric modeling and boundary conditions of homogeneous square plate

下载: 全尺寸图片幻灯片

图 5 不同时刻m收敛结果

Figure 5. Results of m-convergence at different time

下载: 全尺寸图片幻灯片

图 6 不同时刻δ收敛结果

Figure 6. Results of δ-convergence at different time

下载: 全尺寸图片幻灯片

图 7 LdB花岗岩热损伤破裂模型

Figure 7. Thermal damage fracture model for LdB granite

下载: 全尺寸图片幻灯片

图 8 LdB花岗岩破裂声发射结果^[25]

Figure 8. Acoustic emission results of LdB granite

下载: 全尺寸图片幻灯片

图 9 损伤破裂和温度分布演化过程的PD模拟结果

Figure 9. Simulated results of damage fracture and evolution process of temperature distribution by PD

下载: 全尺寸图片幻灯片

图 10 文献[31，32]的花岗岩热破裂声发射结果

Figure 10. Acoustic emission results by Ishida et al

下载: 全尺寸图片幻灯片

图 11 不同温度梯度弛豫时间条件下试件损伤率变化

Figure 11. Variation of damage rate of specimen under different temperature gradient relaxation time conditions

下载: 全尺寸图片幻灯片

图 12 不同温度梯度弛豫时间时x = 0处的温度分布

Figure 12. Distribution of temperature at x = 0 under different temperature gradient relaxation time conditions

下载: 全尺寸图片幻灯片

图 13 不同热流弛豫时间条件下试件损伤率变化

Figure 13. Variation of damage rate of specimen under different heat flux relaxation time conditions

下载: 全尺寸图片幻灯片

图 14 不同热流弛豫时间时x = 0处的温度分布

Figure 14. Distribution of temperature at x = 0 under different heat flux relaxation time conditions

下载: 全尺寸图片幻灯片

表 1 瞬态热传导离散参数值

Table 1 Discrete parameter values for transient heat conduction

条件	物质点分布	点间距/m	近场半径/m	非局部比
m 收敛	250×250	0.040	0.060	1.5
	400×400	0.025	0.060	2.4
	500×500	0.020	0.060	3.0
δ 收敛	250×250	0.040	0.120	3.0
	400×400	0.025	0.075	3.0
	500×500	0.020	0.060	3.0

下载: 导出CSV

参考文献(36)

[1]	ZHANG J Y, SHEN Y J, YANG G S, et al. Inconsistency of changes in uniaxial compressive strength and P-wave velocity of sandstone after temperature treatments[J]. Journal of Rock Mechanics and Geotechnical Engineering, 2021, 13(1): 143-153. doi: 10.1016/j.jrmge.2020.05.008
[2]	郤保平, 吴阳春, 王帅, 等. 热冲击作用下花岗岩力学特性及其随冷却温度演变规律试验研究[J]. 岩土力学, 2020, 41(增刊1): 83-94. XI Baoping, WU Yangchun, WANG Shuai, et al. Experimental study on mechanical properties of granite under thermal shock and its evolution with cooling temperature[J]. Rock and Soil Mechanics, 2020, 41(S1): 83-94. (in Chinese)
[3]	WANG G Y, YANG D, LIU S W, et al. Experimental study on the anisotropic mechanical properties of oil shales under real-time high-temperature conditions[J]. Rock Mechanics and Rock Engineering, 2021, 54(12): 6565-6583. doi: 10.1007/s00603-021-02624-7
[4]	YANG S Q, HUANG Y H, TIAN W L, et al. Effect of high temperature on deformation failure behavior of granite specimen containing a single fissure under uniaxial compression[J]. Rock Mechanics and Rock Engineering, 2019, 52(7): 2087-2107. doi: 10.1007/s00603-018-1725-5
[5]	ZHANG F, ZHAO J J, HU D W, et al. Laboratory investigation on physical and mechanical properties of granite after heating and water-cooling treatment[J]. Rock Mechanics and Rock Engineering, 2018, 51(3): 677-694. doi: 10.1007/s00603-017-1350-8
[6]	唐世斌, 唐春安, 朱万成, 等. 热应力作用下的岩石破裂过程分析[J]. 岩石力学与工程学报, 2006, 25(10): 2071-2078. doi: 10.3321/j.issn:1000-6915.2006.10.019 TANG Shibin, TANG Chun'an, ZHU Wancheng, et al. Numerical investigation on rock failure process induced by thermal stress[J]. Chinese Journal of Rock Mechanics and Engineering, 2006, 25(10): 2071-2078. (in Chinese) doi: 10.3321/j.issn:1000-6915.2006.10.019
[7]	JIAO Y Y, ZHANG X L, ZHANG H Q, et al. A coupled thermo-mechanical discontinuum model for simulating rock cracking induced by temperature stresses[J]. Computers and Geotechnics, 2015, 67: 142-149. doi: 10.1016/j.compgeo.2015.03.009
[8]	严成增. FDEM-TM方法模拟岩石热破裂[J]. 岩土工程学报, 2018, 40(7): 1198-1204. doi: 10.11779/CJGE201807005 YAN Chengzeng. Simulating thermal cracking of rock using FDEM-TM method[J]. Chinese Journal of Geotechnical Engineering, 2018, 40(7): 1198-1204. (in Chinese) doi: 10.11779/CJGE201807005
[9]	YAO J, ZHANG X, SUN Z X, et al. Numerical simulation of the heat extraction in 3D-EGS with thermal-hydraulic- mechanical coupling method based on discrete fractures model[J]. Geothermics, 2018, 74: 19-34. doi: 10.1016/j.geothermics.2017.12.005
[10]	SILLING S A. Reformulation of elasticity theory for discontinuities and long-range forces[J]. Journal of the Mechanics and Physics of Solids, 2000, 48(1): 175-209. doi: 10.1016/S0022-5096(99)00029-0
[11]	SILLING S A, ASKARI E. A meshfree method based on the peridynamic model of solid mechanics[J]. Computers and Structures, 2005, 83(17/18): 1526-1535.
[12]	GAO C L, ZHOU Z Q, LI L P, et al. Strength reduction model for jointed rock masses and peridynamics simulation of uniaxial compression testing[J]. Geomechanics and Geophysics for Geo-Energy and Geo-Resources, 2021, 7(2): 34. doi: 10.1007/s40948-021-00232-x
[13]	马鹏飞, 李树忱, 袁超, 等. 基于SED准则的近场动力学及岩石类材料裂纹扩展模拟[J]. 岩土工程学报, 2021, 43(6): 1109-1117. doi: 10.11779/CJGE202106014 MA Pengfei, LI Shuchen, YUAN Chao, et al. Simulations of crack propagation in rock-like materials by peridynamics based on SED criterion[J]. Chinese Journal of Geotechnical Engineering, 2021, 43(6): 1109-1117. (in Chinese) doi: 10.11779/CJGE202106014
[14]	马鹏飞, 李树忱, 周慧颖, 等. 岩石材料裂纹扩展的改进近场动力学方法模拟[J]. 岩土力学, 2019, 40(10): 4111-4119. MA Pengfei, LI Shuchen, ZHOU Huiying, et al. Simulations of crack propagation in rock-like materials using modified peridynamic method[J]. Rock and Soil Mechanics, 2019, 40(10): 4111-4119. (in Chinese)
[15]	ZHANG H, QIAO P Z. An extended state-based peridynamic model for damage growth prediction of bimaterial structures under thermomechanical loading[J]. Engineering Fracture Mechanics, 2018, 189: 81-97. doi: 10.1016/j.engfracmech.2017.09.023
[16]	WANG Y T, ZHOU X P, KOU M M. An improved coupled thermo-mechanic bond-based peridynamic model for cracking behaviors in brittle solids subjected to thermal shocks[J]. European Journal of Mechanics-A, 2019, 73: 282-305. doi: 10.1016/j.euromechsol.2018.09.007
[17]	WANG Y T, ZHOU X P. Peridynamic simulation of thermal failure behaviors in rocks subjected to heating from boreholes[J]. International Journal of Rock Mechanics and Mining Sciences, 2019, 117: 31-48. doi: 10.1016/j.ijrmms.2019.03.007
[18]	YANG Z, YANG S Q, CHEN M. Peridynamic simulation on fracture mechanical behavior of granite containing a single fissure after thermal cycling treatment[J]. Computers and Geotechnics, 2020, 120: 103414. doi: 10.1016/j.compgeo.2019.103414
[19]	YANG Z, YANG S Q, TIAN W L. Peridynamic simulation of fracture mechanical behaviour of granite specimen under real-time temperature and post-temperature treatments[J]. International Journal of Rock Mechanics and Mining Sciences, 2021, 138: 104573. doi: 10.1016/j.ijrmms.2020.104573
[20]	CATTANEO C. A form of heat conduction equation which eliminates the paradox of instantaneous propagation[J]. Compte Rendus, 1958, 247: 431-433.
[21]	VERNOTTE P. Paradoxes in the continuous theory of the heat conduction[J]. Compte Rendus, 1958, 246: 3154-3155.
[22]	TZOU D Y. A unified field approach for heat conduction from macro- to micro-scales[J]. Journal of Heat Transfer, 1995, 117(1): 8-16. doi: 10.1115/1.2822329
[23]	AL-NIMR M A, NAJI M, ABDALLAH R I. Thermal behavior of a multi-layered thin slab carrying periodic signals under the effect of the dual-phase-lag heat conduction model[J]. International Journal of Thermophysics, 2004, 25(3): 949-966. doi: 10.1023/B:IJOT.0000034247.32646.d4
[24]	WANG L J, XU J F, WANG J X. A peridynamic framework and simulation of non-Fourier and nonlocal heat conduction[J]. International Journal of Heat and Mass Transfer, 2018, 118: 1284-1292. doi: 10.1016/j.ijheatmasstransfer.2017.11.074
[25]	JANSEN D P, CARLSON S R, YOUNG R P, et al. Ultrasonic imaging and acoustic emission monitoring of thermally induced microcracks in Lac du Bonnet granite[J]. Journal of Geophysical Research: Solid Earth, 1993, 98(B12): 22231-22243. doi: 10.1029/93JB01816
[26]	MADENCI E, OTERKUS E. Peridynamic Theory and Its Applications[M]. New York: Springer, 2014.
[27]	HOSSEINI-TEHRANI P, ESLAMI M R. BEM analysis of thermal and mechanical shock in a two-dimensional finite domain considering coupled thermoelasticity[J]. Engineering Analysis with Boundary Elements, 2000, 24(3): 249-257. doi: 10.1016/S0955-7997(99)00063-6
[28]	WANG Y T, ZHOU X P, XU X. Numerical simulation of propagation and coalescence of flaws in rock materials under compressive loads using the extended non-ordinary state-based peridynamics[J]. Engineering Fracture Mechanics, 2016, 163: 248-273. doi: 10.1016/j.engfracmech.2016.06.013
[29]	WANNE T S, YOUNG R P. Bonded-particle modeling of thermally fractured granite[J]. International Journal of Rock Mechanics and Mining Sciences, 2008, 45(5): 789-799. doi: 10.1016/j.ijrmms.2007.09.004
[30]	GOFFREDO DE PORTU. Introduction to Mechanical Behaviour of Ceramics[M]. Faenza: CNR, 1992.
[31]	ISHIDA T, KINOSHITA N, WAKABAYASHI N. Acoustic emission monitoring during thermal cracking of a granite block heated in a center hole[C]//Proceedings of the 3rd Asian Rock Mechanics Symposium. Kyoto, 2004.
[32]	HUANG X, TANG S B, TANG C A, et al. Numerical simulation of cracking behavior in artificially designed rock models subjected to heating from a central borehole[J]. International Journal of Rock Mechanics and Mining Sciences, 2017, 98: 191-202. doi: 10.1016/j.ijrmms.2017.07.016
[33]	胡学功, 刘登瀛, 周定伟. 细陶瓷棒的脉冲热流波动传播行为[J]. 应用科学学报, 2002, 20(3): 305-308. doi: 10.3969/j.issn.0255-8297.2002.03.019 HU Xuegong, LIU Dengying, ZHOU Dingwei. The behavior of pulse heat flux wave propagation in a fine round ceramic bar[J]. Journal of Applied Sciences, 2002, 20(3): 305-308. (in Chinese) doi: 10.3969/j.issn.0255-8297.2002.03.019
[34]	KAMINSKI W. Hyperbolic heat conduction equation for materials with a nonhomogeneous inner structure[J]. Journal of Heat Transfer, 1990, 112(3): 555-560. doi: 10.1115/1.2910422
[35]	GUO S L, WANG B L. Thermal shock fracture of a cylinder with a penny-shaped crack based on hyperbolic heat conduction[J]. International Journal of Heat and Mass Transfer, 2015, 91: 235-245. doi: 10.1016/j.ijheatmasstransfer.2015.07.081
[36]	GUO S L, WANG B L. Thermal shock cracking behavior of a cylinder specimen with an internal penny-shaped crack based on non-fourier heat conduction[J]. International Journal of Thermophysics, 2016, 37(2): 17. doi: 10.1007/s10765-015-2029-6