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

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

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
计量
- 文章访问数: 328
- HTML全文浏览量: 46
- PDF下载量: 65
出版历程
- 收稿日期: 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]。例如，在深部地热开采过程中，岩石的热损伤和破裂形成了连通的裂隙网络，可作为工作流体的流动通道。通过理解和利用热损伤破裂特征，可以增强裂隙的扩展，提高采热效率。近年来，众多学者进行了关于岩石热损伤破裂的试验研究。郤保平等^[2]通过对花岗岩进行热冲击破裂处理，总结了其宏观力学性质的演变规律；Wang等^[3]利用SEM和CT技术分析了油页岩的各向异性力学性能；Yang等^[4]研究了含预制裂隙高温花岗岩在自然冷却后的力学性能；Zhang等^[5]研究了加热和水冷过程对花岗岩物理力学性能以及微裂纹演化的影响。然而，由于设备的限制，同时获得岩石损伤破裂与温度的演化过程仍较为困难。同时，对不同数值的变量分析要进行大量平行试验，成本较高。

数值模拟是解决以上问题的有效途径。唐世斌等^[6]利用基于有限元法（FEM）的RFPA系统模拟了含单个内嵌颗粒岩石的热破裂行为。Jiao等^[7]利用非连续变形分析法（DDA）、严成增^[8]基于有限元-离散元方法模拟了圆盘型岩石的热破裂过程。然而，传统的连续介质力学数值方法，例如FEM，在模拟岩石破裂问题时存在应力奇异性、网格依赖性等局限。常见的非连续介质力学数值方法，如离散元法（DEM）^[9]，需要进行繁琐的参数标定，这对计算效率和精度均会产生影响。DDA和数值流形法（NMM）也各自存在网格依赖性和裂纹不能扩展到单元内部等局限。

Silling等^[10-11]提出了采用积分形式控制方程的近场动力学（PD）理论。与传统连续介质力学理论的“局部”思想不同，PD理论将固体离散为空间域内一系列包含所有物性信息的物质点，并假定每一点的对点力或热流密度受到有限半径区域内其它所有物质点的影响，即这是一种“非局部”理论。积分方程的特性使其有效避免了应力奇异性、网格依赖性等局限。PD理论从最初只能模拟固定泊松比材料的键基PD模型发展至后来不受泊松比限制的常规态基模型和非常规态基模型。Gao等^[12]、马鹏飞等^[13-14]的研究证明了PD模型在模拟荷载作用下岩石断裂问题上的有效性。

PD模型同样适用于模拟热力耦合作用下的材料损伤破裂问题。Zhang等^[15]提出了一种考虑热效应的常规态基PD模型，成功预测了陶瓷等双材料结构的损伤演化路径。Wang等^[16-17]模拟了单中心钻孔和两个对称偏心钻孔岩石在加热条件下的热破裂行为以及陶瓷淬火试验中热致裂纹的周期性和层次性特征。Yang等^[18-19]构建了基于常规态基PD理论的热力耦合模型，研究了实时高温和高温冷却处理后花岗岩在压缩荷载作用下的断裂特征。迄今为止，PD热力耦合模型基本都是基于傅里叶定律推导得到的。

傅里叶定律将导热视为一种扩散行为，假定基于热扰动的传播速度是无限大的，适用于模拟稳态传热过程。但对于某些极端环境下的强瞬态传热问题，例如在干热岩地热能开采过程中，低温冷水与高温岩体的温差很大，二者接触时岩体表面温度急剧变化，新的热平衡建立需要一定的时间，即温度场重新达到平衡的时间晚于热扰动施加的时间（这一时间差称为热驰豫时间），此时傅里叶定律严格上是不成立的。在岩石矿物颗粒非均匀膨胀以及温度梯度的作用下，脆性断裂在很短的时间内发生，因此非傅里叶定律或许更适合描述此时的热传导及岩石的热损伤破裂过程。

为克服傅里叶定律的局限，Cattaneo^[20]和Vernotte^[21]提出了热流和温度梯度之间存在热弛豫时间的双曲型热传导方程，但是该模型只反映了热流滞后于温度梯度，且在强瞬态热作用下单项延迟的解可能不够准确。Tzou^[22]引入热流密度滞后项和温度梯度滞后项，提出了反映非傅里叶导热现象的双相滞后（DPL）模型，这一模型从宏观角度描述了材料微观结构在很小的时间尺度内对热传导的影响。Al-Nimr等^[23]模拟了在周期性热载荷作用下，基于DPL模型的多层板热响应特征。Wang等^[24]将非傅里叶定律引入PD传热理论中，建立了相应的控制方程并验证了温度解的唯一性。但据笔者所知，目前尚未见基于非傅里叶定律条件下岩石热损伤破裂问题的PD研究。

基于此，首先通过引入双相滞后模型，提出了一个基于非傅里叶热传导定律的热力耦合模型，该模型考虑了破裂对热传导的影响。以Jansen等^[25]的试验结果为参照，采用热力耦合PD模型模拟了Lac du Bonnet（LdB）花岗岩的热损伤破裂过程以及温度场演化过程。然后，研究了温度梯度弛豫时间、热流弛豫时间对岩石热损伤破裂特征及温度分布的影响，以期为深部干热岩地热开采等工程提供理论参考。

1. 近场动力学热力耦合基本理论

1.1 PD运动方程

如图 1所示，近场动力学理论假设在任一时刻t，物体区域R内任一物质点x与其周围一定区域H_x内的其他任意物质点x'通过键相互作用，区域H_x称为近场域，其半径δ为两个物质点间的最大相互作用距离，定义为 ${H_x} = \left\{ {x \in R:\left| {x - x' \leqslant \delta } \right|} \right\}$ 。

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

Figure 1. Material points of peridynamic interactions

下载: 全尺寸图片幻灯片

F_f (x - x', t)和F_f (x' -x, t)分别为物质点x对物质点x'作用和物质点x'对物质点x作用的不包含温度作用的对点力密度矢量。而F_T (x -x', t)和F_T (x' -x, t)则是由于温差产生的物质点x和物质点x'之间相互作用的对点力密度矢量。

如图 2所示，物质点x和物质点x'的初始相对位置矢量和相对位移矢量定义为ξ = x - x'和η = u_x -u_x'. u_x和u_x'为物质点x和x'的位移矢量。物质点y和y'为物质点x和x'变形后的位置，有ξ + η = y - y'。

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

Figure 2. Schematic diagram of peridynamic deformation

下载: 全尺寸图片幻灯片

物质点x在时刻t的运动方程为

$\rho \frac{{{\partial ^2} \boldsymbol{u}(x,t)}}{{\partial {t^2}}} = \int_{{H_x}} {\left\{ {\left[ {{ \boldsymbol{F}_{\text{f}}}(x' - x,t) - {\boldsymbol{F}_{\text{T}}}(x' - x,t)} \right] - } \right.} \\ \text{ } \left. {\left[ {{\boldsymbol{F}_{\text{f}}}(x - x',t) - {\boldsymbol{F}_{\text{T}}}(x - x',t)} \right]} \right\}{\text{d}}{V_{x'}} + \boldsymbol{b}(x,t) 。$

(1)

式中：ρ为质量密度；b为体力密度矢量；V_x'为物质点x'的体积。

对于二维各向同性弹性材料，在时刻t的常规态基PD模型可表示为^[26]

$\rho \frac{{{\partial ^2}\boldsymbol{u}(x,t)}}{{\partial {t^2}}} = \int_{{H_x}} {2\delta \left( {\frac{{y' - y}}{{\left| {y' - y} \right|}}} \right)\left\{ {\left( {\frac{{\varPsi d}}{{\left| \xi \right|}}} \right)\left[ {{a_1}({\theta _x} + {\theta _{x'}}) - } \right.} \right.} {\text{ }} \\ \left. {\left. {\frac{1}{2}{a_2}({T_x} + {T_{x'}})} \right] + b\left[ {2s - \alpha ({T_x} + {T_{x'}})} \right]} \right\}{\text{d}}{V_{x'}} + \boldsymbol{b}(x,t) 。$

(2)

式中：T_x，T_x'分别为物质点x和x'的温度变化量；α为线性热膨胀系数；θ_x和θ_x'分别为物质点x和x'的体积应变^[26]，

$\left.\begin{array}{l} \theta_x=\int_{H_x} \delta \mathrm{~d}\left(s-\alpha T_x\right) \varPsi \mathrm{d} V_{x^{\prime}}+2 \alpha T_x, \\ \theta_{x^{\prime}}=\int_{H_{x^{\prime}}} \delta \mathrm{d}\left(s-\alpha T_{x^{\prime}}\right) \varPsi \mathrm{d} V_x+2 \alpha T_{x^{\prime}}, \\ \varPsi=\left(\frac{y^{\prime}-y}{\left|y^{\prime}-y\right|}\right) \cdot\left(\frac{x^{\prime}-x}{\left|x^{\prime}-x\right|}\right)。 \end{array}\right\}$

(3)

s为在物质点x和x'间的键伸长率，

$s = \frac{{\left| {y' - y} \right| - \left| {x' - x} \right|}}{{\left| {x' - x} \right|}} 。$

(4)

式（2）中的a₁，a₂，b，d为近场动力学计算参数，通过通过等效连续介质力学理论框架与近场动力学理论框架内的应变能密度，可将它们用连续介质力学中的弹性模量E、泊松比 $\nu$ 以及剪切模量S表示为^[26]

$\left.\begin{array}{l} a_1=\frac{1}{2}\left(\frac{E}{2(1-v)}-2 S\right), a_2=4 \alpha a_1, \\ b=\frac{6 S}{\pi h \delta^4}, d=\frac{2}{\pi h \delta^3} 。 \end{array}\right\}$

(5)

1.2 基于傅里叶定律的PD热传导方程

傅里叶导热定律是指在单位时间内通过单位截面的热量与温度梯度成正比，其表达式为

$\boldsymbol{q}(x,t) = - {k_{\text{T}}}\nabla \varTheta (x,t) 。$

(6)

式中：q为热通量矢量；k_T为宏观导热系数； $\nabla \varTheta$ 为温度梯度。

在任意时刻t，PD热传导方程如下^[26]：

$\rho c\frac{{\partial \varTheta (x,t)}}{{\partial t}} = Q(x{\text{,}}t) + {b_{\text{h}}}(x{\text{,}}t) \text{，}$

(7)

$Q(x{\text{,}}t) = \int_{{H_x}} {\left[ {{F_h}(x' - x,t) - {F_{\text{h}}}(x - x',t)} \right]{\text{d}}{V_{x'}}} 。$

(8)

式中：c为比热容；Q为热流量；F_h (x'- x, t), F_h (x - x', t)分别为物质点x和x'的热流密度；b_h为热源密度，即单位时间单位体积的热源产生的热量。

假定物质点x和x'之间的热通量仅是这两个点之间温度差的函数，则有^[26]

${F_{\text{h}}}(x' - x{\text{,}}t) - {F_{\text{h}}}(x' - x{\text{,}}t) = {\kappa _{\text{T}}}\frac{{\varTheta (x'{\text{,}}t) - \varTheta (x{\text{,}}t)}}{{\left| {x' - x} \right|}} 。$

(9)

式中：Θ(x, t)，Θ(x', t)分别为物质点x和x'的温度； ${\kappa _{\text{T}}}$ 为PD微导热系数，其与宏观导热系数 ${k_{\text{T}}}$ 的关系为^[26]

${\kappa _{\text{T}}} = \frac{{6{k_{\text{T}}}}}{{{\text{π }}h{\delta ^3}}} 。$

(10)

2. 基于非傅里叶定律的近场动力学热力耦合模型

2.1 基于DPL模型的PD热传导方程

Tzou^[22]扩展了经典傅里叶导热定律，提出了同时考虑温度梯度和热流驰豫时间的双相滞后模型，这是非傅里叶导热定律中的一种模型，其表达式为

$\boldsymbol{q}(x\text{,}t+{\tau }_\text{q})=-{k}_\text{T}\nabla \varTheta (x\text{,}t+{\tau }_\text{T})\text{ (}{\tau }_\text{q} > 0,\text{ }{\tau }_\text{T} > 0\text{) }。$

(11)

式中： ${\tau _\text{q}}$ ， ${\tau _\text{T}}$ 分别为热流弛豫时间和温度梯度弛豫时间。

将DPL模型引入PD热传导理论框架中，则式（7），（8）可改写为

$\rho c\frac{{\partial \varTheta (x{\text{,}}t + {\tau _\text{q}})}}{{\partial t}} = Q(x{\text{,}}t + {\tau _\text{q}}) + {b_\text{h}}(x{\text{,}}t + {\tau _\text{q}}) \text{，}$

(12)

$Q(x\text{,}t+{\tau }_\text{q})={\displaystyle {\int }_{{H}_{x}}[{F}_\text{h}({x}^{\prime }-x\text{,}t+{\tau }_\text{T})-{F}_\text{h}(x-{x}^{\prime }\text{,}t+{\tau }_\text{T})]\text{d}{V}_{{x}^{\prime }}}。$

(13)

将式（13）代入式（12）中，并对其左右两端关于t的函数项进行一阶泰勒展开，经过整理可得

$\begin{array}{l} \rho c\frac{{\partial \varTheta (x,t)}}{{\partial t}} + {\tau _{\text{q}}}\rho c\frac{{{\partial ^2}\varTheta (x,t)}}{{\partial {t^2}}} \hfill \\ = \int_{{H_x}} {\left[ {{F_\text{h}}(x' - x,t) - {F_\text{h}}(x - x',t)} \right]\text{d}{V_{x'}}} + \hfill \\ {\tau _{\text{T}}}\int_{{H_x}} {\frac{{\partial \left[ {{F_\text{h}}(x' - x,t) - {F_\text{h}}(x - x',t)} \right]}}{{\partial t}}\text{d}{V_{x'}} + } \hfill \\ {b}_\text{h}(x,t)+{\tau }_\text{q}\frac{\partial {b}_\text{h}(x,t)}{\partial t}\text{ }。 \end{array}$

(14)

式（14）即为基于DPL模型推导得到的PD热传导方程，其与式（2）PD运动方程共同组成了基于非傅里叶定律的PD热力耦合模型。

2.2 断裂准则

PD理论不需要引入外部准则来确定材料是否破坏和损伤路径。通过键的断裂来定义损伤，物质点x处的局部损伤值φ定义为^[26]

$\varphi (x,t) = 1 - \frac{{\int_{{H_x}} {\mu (x,x',t) \text{d}{V_{x'}}} }}{{\int_{{H_x}} {\text{d}{V_{_{x'}}}} }} 。$

(15)

式中：μ为表征键是否断裂的标量函数，

$\mu (x,x',t) = \left\{ \begin{gathered} 1{\text{ (}}s < {s_0}) \hfill \\ 0{\text{ (}}s \geqslant {s_0}) \hfill \\ \end{gathered} \right. 。$

(16)

式中：s₀为键的临界伸长率，与经典断裂力学中的能量释放率有关。当s小于临界伸长率s₀时，两物质点间存在正常相互作用，键保持完整，而s超过临时伸长率s₀时，键将被破坏且无法恢复。根据各物质点的局部损伤值便可以得到岩石宏观损伤破裂路径。

2.3 数值实现

如所示，将整个模型均匀地离散成多个子域，每一个子域的物理性质则通过子域的几何中心点来传递，任意两个中心点的距离为 $\Delta x$ 。在某一时刻t，物质点x与其近场域范围内的点x'相互作用，因此式（2），（14）可以用黎曼和的形式改写为

$\begin{array}{c} \rho \frac{{{\partial ^2}\boldsymbol{u}_i^n}}{{\partial {t^2}}} = \sum\limits_{j = 1}^{{N_i}} {\left[ {2\delta \left( {\frac{{{\boldsymbol{y}_j} - {\boldsymbol{y}_i}}}{{\left| {{\boldsymbol{y}_j} - {\boldsymbol{y}_i}} \right|}}} \right)\left\{ {\left( {\frac{{\Psi d}}{{\left| \xi \right|}}} \right)} \right.} \right.\left[ {{a_1}({\theta _i} + {\theta _j}) - \frac{1}{2}{a_2}({T_i} + {T_j})} \right] + } \\ \left. {b\left[ {2s - \alpha ({T_i} + {T_j})} \right]} \right\}\left. {{V_j}} \right] + \boldsymbol{b}_i^n \text{，} \end{array}$

(17)

$\begin{array}{c} \rho c\frac{{\partial \varTheta _i^n}}{{\partial t}} + {\tau _\text{q}}\rho c\frac{{{\partial ^2}\varTheta _i^n}}{{\partial {t^2}}} = \sum\limits_{j = 1}^{{N_i}} {\left[ {{\kappa _\text{T}}\frac{{\varTheta _j^n - \varTheta _i^n}}{{\left| {{x_j} - {x_i}} \right|}} + {\tau _\text{T}}\frac{{\partial {F}_{hij}}}{{\partial t}}} \right]{V_j}} + \\ b_{{\text{h}}i}^n + {\tau _\text{q}}\frac{{\partial b_{{\text{h}}i}^n}}{{\partial t}} 。 \end{array}$

(18)

图 3 PD模型离散示意图

Figure 3. Diagrammatic sketch of PD model discretization

下载: 全尺寸图片幻灯片

式中：i为材料的第i个物质点；j为物质点x_i近场域范围内的第j个物质点；N_i为物质点x_i近场域范围内子域的总数；n为时间步数；V_j为点x_j的体积。对于图 3中靠近近场域中心的其它物质点来说，其所代表的体积完全处在近场域内，然而对于处于近场域边界上的物质点，只有部分子域位于近场域范围内，需要根据边界上的物质点与近场域半径之间的关系按比例进行折减，节点体积V_x'可以表示为

$V_{x^{\prime}}= \begin{cases}(\Delta x)^2 & (|\xi| \leqslant(\delta-r)) \\ \left(\frac{\delta+r-|\xi|}{2 r}\right)(\Delta x)^2 & ((\delta-r) <|\xi| \leqslant \delta)。 \\ 0 & (|\xi|>\delta)\end{cases}$

(19)

式中：r为网格间距∆x的一半。

通过向前差分技术求解式（17），（18）：

$\left.\begin{array}{l} \frac{\partial \boldsymbol{u}_i^{n+1}}{\partial t}=\frac{\partial^2 \boldsymbol{u}_i^n}{\partial t^2} \Delta t_M+\frac{\partial \boldsymbol{u}_i^n}{\partial t}, \\ \boldsymbol{u}_i^{n+1}=\frac{\partial \boldsymbol{u}_i^{n+1}}{\partial t} \Delta t_M+\boldsymbol{u}_i^n, \end{array}\right\}$

(20)

$\left.\begin{array}{l} \frac{\partial \varTheta_i^{n+1}}{\partial t}=\frac{\partial^2 \varTheta_i^n}{\partial t^2} \Delta t_T+\frac{\partial \varTheta_i^n}{\partial t}, \\ \varTheta_i^{n+1}=\frac{\partial \varTheta_i^{n+1}}{\partial t} \Delta t_T+\varTheta_i^n 。 \end{array}\right\}$

(21)

式中：∆t_M，∆t_T分别为位移场和温度场时间步长。

3. 模型可行性验证

3.1 瞬态热传导模拟

为验证PD模型在模拟瞬态热传导问题上的有效性，对一个尺寸为10 m×10 m的二维均质方板受到热冲击作用的问题进行了模拟，得到了点A（-3, 0）处的温度随时间变化的曲线，并与Hosseini-Tehrani等^[27]使用边界元法（BEM）得到的结果对比。图 4为均质方板的几何模型。材料密度ρ = 1 kg/m³，比热容c= 1 J/kg/K，导热系数为k_T = W/m/K。温度时间步长∆t_T = 5×10^-4 s。模型初始温度为0℃，左侧边界受到随时间变化的热冲击作用，其余边界热流密度为0。初始和边界条件如下：

$\begin{array}{l}\varTheta (x,y,t=0)=0\text{ (}\varTheta (x=-5,t)=5t\cdot {e}^{-2t}\text{)}\text{，}\\ {\varTheta }_{,x}(x=5,y)=0\text{ (}{\varTheta }_{,y}(x,y=\pm 5)=0,\text{ }t > 0\text{)}。\end{array}\}$

(22)

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

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

下载: 全尺寸图片幻灯片

在PD模拟中，为了测试数值收敛性与模拟结果的准确性，需要进行两种数值收敛分析，分别称为m收敛和δ收敛，它们之间的关系式为

$\delta = m \cdot \Delta 。$

(23)

式中：m为近场域半径δ与物质点间距 $\Delta$ 的比值。

分别模拟了δ固定为0.06 m及m固定为3时，点A的温度变化曲线，具体方案表 1所示。

表 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

| 显示表格

图 5显示了不同m值情况下点A处的温度值随时间的变化曲线。可以看到当近场域尺寸δ固定为0.06 m时，随着非局部比m的增大，PD结果与BEM结果的偏离程度逐渐减小。图 6显示了不同δ值情况下点A处的温度值随时间的变化曲线。可以看到当m固定为3时，随着δ的减小，PD结果与BEM结果逐渐趋于一致。

图 5 不同时刻m收敛结果

Figure 5. Results of m-convergence at different time

下载: 全尺寸图片幻灯片

图 6 不同时刻δ收敛结果

Figure 6. Results of δ-convergence at different time

下载: 全尺寸图片幻灯片

3.2 LdB花岗岩热损伤破裂模拟

如图 7（a）所示，Jansen等^[25]选择LdB花岗岩试件，在其中心凿一个圆孔，试件直径为30 cm，中央圆孔直径为3 cm。通过中央圆孔处的加热器对其进行加热，并通过声发射技术和超声波断层扫描技术监测其在无侧限条件下的热损伤破裂位置。PD数值模拟采用的几何模型如图 7（b）所示：将LdB花岗岩试件简化为一个带有中央圆孔的二维圆形模型，以中央圆孔中心为坐标系原点，内外径尺寸与试验相同。材料性质参数如下：弹性模量E = 100 GPa，密度ρ = 2650 kg/m³，泊松比 $\nu$ = 0.33，导热系数k_T = 3.5 W/m/K，比热容c= 880 J/kg/K，热膨胀系数α= 7.5×10^-6 1/K。将试样离散为30000个物质点，节点间距∆x = 1 mm，近场域半径δ = 3.015 mm。参考以往典型的岩石破裂PD研究^{[13, 28]}，采用Weibull随机分布函数反映试件的非均质性，假定临界伸长率服从尺度参数为0.002，形状参数为10的Weibull分布。温度时间步长∆t_T = 5×10^-2 s，位移时间步长∆t_M = 5×10^-8 s. 根据试验设置，试件整体初始温度Θ₀= 20℃，中央圆孔边界温度Θ_r= 200℃，试验外边界为自由边界条件。

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

Figure 7. Thermal damage fracture model for LdB granite

下载: 全尺寸图片幻灯片

图 8为试验得到的声发射结果。图 8中的圆点为声发射活动发生的位置，表示变形破裂过程的发生。可以看到，LdB花岗岩外边界处的宏观热损伤与中央圆孔相连接。Goffredo^[30]从连接介质力学的角度解释了类似的现象，认为在温度梯度的驱动下，材料外部温度较低的部分会受到拉伸应力的作用，并且内部区域处于压缩状态。而岩石等脆性材料的抗拉强度远小于其抗压强度，因此热破裂首先从外边界表面处产生。

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

Figure 8. Acoustic emission results of LdB granite

下载: 全尺寸图片幻灯片

图 9为PD模拟得到的热损伤破裂演化过程和温度场分布结果。可以看到热损伤破裂从试件中部正上方产生，然后逐渐向中央圆孔扩展。除此之外，模拟结果还显示靠近试件外边界的损伤面积明显大于内部损伤面积。这是因为热损伤首先出现在岩石外边界处，在损伤向内扩展的过程中，损伤萌生部位仍发生变形和破裂。因此，试件外边界的损伤程度最高，并且沿着扩展路径逐渐减小。另外，尽管模拟的几何条件与边界条件是对称的，但由于本文采用了Weibull分布函数来描述岩石的非均质性，在受到温度变化的影响时，非均质性导致应力分布的局部差异，因此结果呈现出非对称性。在Ishida等^[31]、Huang等^[32]进行的类似的花岗岩中心孔加热破裂试验中，也观察到了相同的特征（图 10），表明这种现象并非偶然的。

图 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

下载: 全尺寸图片幻灯片

另外，从图 9中的温度场演化结果可以看到，在第350秒和第375秒的灰色椭圆标记部位，温度梯度出现了明显的不连续。这是因为在模拟中材料物质点间键的断裂导致热流密度中断，体现在宏观现象中则是空气的导热系数相比于固体材料（岩石）可以忽略不计，因此热破裂导致岩石裂纹面之间无法传热。而其它的数值方法例如FEM、DEM等则难以捕捉到因破裂产生的温度不连续特征，从这个意义上讲，基于非局部思想的PD模型可以得到更真实的温度场。

总之，尽管由于非均质性、试验误差等因素没有进行应力等变量的定量对比，但PD模拟结果仍然从多个定性角度重现了热力耦合作用下花岗岩热损伤破裂演化过程的主要特征：①热损伤破裂从岩石的外边界表面处开始萌生，然后向岩石内部中心圆孔扩展。②由于岩石材料的非均质性，率先破裂的位置导致应力重分布，因此岩石的热损伤破裂自始至终只沿一条路径扩展，与试验结果一致。③在热损伤向内扩展的过程中，起裂处仍然会发生变形和破裂，因此起点处的损伤程度最高，并沿着扩展路径逐渐减小，与试验结果特征相符。④模拟结果中观察到温度分布出现断续现象，反映了裂纹面无法传热的特征，也凸显了PD方法在模拟热力耦合损伤破裂问题上的优越性。

4. 热弛豫时间影响分析

热驰豫时间有助于预测和评估岩石对瞬态温度变化的响应速度，从而评估岩石的损伤断裂特征。以3.2节中的LdB花岗岩为例，进行不同热驰豫时间条件下的岩石热损伤破裂PD模拟。DPL模型中包括热流弛豫时间τ_q和温度梯度弛豫时间τ_T。一般金属的热驰豫时间小于10^-11量级，但岩石作为一种非均质材料，其热弛豫时间并不局限于分子或晶格水平^[33]，因此热弛豫时间相对较大。目前鲜有关于岩石热弛豫时间测定的研究，Kaminski^[34]测定得到同为多孔介质的砂的热弛豫时间最大可达26.4 s。为得到区分度相对明显的结果，本次模拟中τ_q和τ_T的值分别取0，1.0，1.5，2.0 s。

引入损伤率（产生局部损伤的物质点数量占总数量的百分比）作为衡量试件热损伤破裂程度的指标。图 11为τ_q固定时，不同τ_T情况下试件损伤率的变化情况。可以看到损伤率随着τ_T的增加而增加。图 12为不同τ_T情况下第375秒时试件中轴（x = 0）上各点的温度值。可以看到由于试件热损伤破裂的发生，沿轴两侧的温度分布并不是对称的，在y > 0一侧的温度分布曲线呈现出现明显的断续特征。

图 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

下载: 全尺寸图片幻灯片

有研究表明温度梯度弛豫时间对热传导在一定空间内其促进作用^[35-36]。从图 12可以看出，当τ_q固定时，τ_T的增加对热传导的促进作用越来越明显，试件内部温度的响应变快，相同位置处的温度越来越高。这与Tzou^[22]得到的热驰豫时间对热传导影响的结论一致。结合图 11看出：各位置升温速率加快使得相应的温度梯度增大，导致热应力变大，因此试件的起裂时间更早且损伤程度加剧。另外，当温度梯度弛豫时间不存在时时，热流弛豫时间超过1.5 s后会导致试件损伤的时间明显延迟，相应的温度曲线变得平滑，没有明显的跳跃。

图 13为τ_T固定时，不同τ_q情况下试件损伤率的变化情况。图 14为不同τ_q情况下第375秒时试件中轴（x = 0）上各点的温度值。当τ_T固定时，热流弛豫时间的增加使得热量传输路径上的温度响应出现更多的延迟，温度梯度弛豫时间被热流弛豫时间的影响抵消，因此图 13中显示试件的起裂时间更晚并且损伤破裂程度减弱。然而从图 14中可以看出，当τ_T固定时，τ_q的增加仍然使得x = 0处各点位置的温度升高，推测这是因为相比于温度梯度弛豫时间，相同差值的热流弛豫时间对热传导的影响很小，虽然温度值增加但相应的温度梯度减小所导致的。

图 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

下载: 全尺寸图片幻灯片

5. 结论

（1）通过引入双相滞后热传导模型，在近场动力学理论框架内推导并建立了基于非傅里叶定律的热力耦合模型。通过模拟平板瞬态热传导问题和LdB花岗岩热损伤破裂问题验证了模型的可行性。

（2）LdB花岗岩首先从外边界处开始破裂，然后沿一个方向向内部中心圆孔扩展，损伤面积沿扩展路径从外向内逐渐减小。PD热力耦合模型模拟结果很好地反映了上述特征及温度分布的不连续现象。

（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