CN113326551B - 热力耦合条件下的围岩开挖损伤分析方法及其应用 - Google Patents

Abstract

本发明公开了热力耦合条件下的围岩开挖损伤分析方法及其应用，所述分析方法包括：建立隧洞的开挖***‑卸荷‑降温全过程模型，对其包括热应力、动应力及地应力叠加耦合的过程控制方程进行求解，从而获得模型在不同地质条件与开挖方式下的岩石热‑力耦合响应。本发明的分析方法不仅可确定开挖损伤区的分布范围，还能揭示岩体损伤的形成阶段和形成原因，同时具有计算效率高、所需时间短和耗费人力物力小等优点，可为高温高压条件下的硐室适应性开挖设计提供预测与指导。

Description

热力耦合条件下的围岩开挖损伤分析方法及其应用

技术领域

本发明涉及地下硐室开挖的技术领域。

背景技术

随着地表浅层矿物资源的日益枯竭，资源勘探和开采活动不断走向地球深部。但据Lucier等报道，南非某金矿中实测的地应力量级最高可达100MPa，据Belle和Biffi的研究，澳大利亚某煤矿目前揭露的地温梯度则高达7℃/100m。由此可见，资源赋存深度增加引起的高地应力与高地温问题为深部资源开采、深部能源开发和高放废物地质处置等深部岩土工程实践构成了严俊的技术挑战。

在地下硐室开挖时，保留岩体会由于***荷载冲击和动态应力重分布形成不可恢复的动力损伤。此外，对于地热异常区或深部高温矿井，在高温硐室开挖完成后，通常采用人工通风等方式将冷空气导入地下巷道以降低工作环境温度，因此受冷却作用的岩体将产生温度重分布并引起二次应力场的不断调整，从而可能在已形成的动态开挖损伤基础上产生额外的热损伤，严重威胁硐室稳定和施工安全。因此，建立深埋高温硐室开挖***-卸荷-降温全过程的岩石热-力耦合行为模型和分析方法，进而揭示围岩开挖损伤的形成机制、力学机理和影响范围，可为高地应力与高地温环境下的深部硐室开挖设计提供重要的实践帮助。

为解决上述问题，现有技术中提出了一些可行方案，如罗忆等通过数值模拟法，利用有限元研究了高地应力条件下***荷载和瞬态卸荷诱发的围岩损伤区分布情况，但该方法存在以下缺点：

(1)数值计算通常受到网格划分方式的限制，若模型所划分的网格较密，则计算效率低，所需时间长；若模型网格较粗，计算结果的准确性得不到有效保证；

(2)未考虑高地应力与高地温耦合作用下的岩石热-力耦合响应，因此无法描述深部高温高压条件下的围岩开挖损伤演化机制、无法得到该条件下的准确操作方案。

或如中国专利申请CN201911400238.2提出了一种在隧道开挖区布设声波监测点，并根据获取的声波监测数据计算岩体损伤值的方法，但该方法的存在以下缺点：

(1)声波监测点的布置较为繁琐，监测过程所需时间长，人力物力消耗大。

(2)实际监测法仅能获取围岩损伤区的分布范围，而不能反映***-卸荷-降温过程中的岩石应力演化趋势，进而难以揭示开挖损伤的形成机制和力学机理，无法得到精确的操作方案。

发明内容

针对上述现有技术的缺陷，本发明的目的在于提出一种热力耦合条件下的围岩开挖损伤分析方法，其基于热力耦合解析和应力轨迹投影等过程，可准确判断不同地质条件和开挖方式下的围岩损伤区形成机制、力学机理和影响范围，从而可进一步获得不同条件和方式下的最优操作方案，克服了传统数值计算和原位监测方法不仅计算效率低、所需时间长、耗费人力物力大，而且无法区分围岩损伤的形成机制和力学机理等局限。

本发明的目的还在于提供上述方法的一些具体应用。

本发明首先公开了如下的技术方案：

热力耦合条件下的围岩开挖损伤分析方法，其包括：

建立包含变量参数的隧洞开挖全过程热力耦合模型，所述热力耦合模型包括地应力、动态应力及热应力的耦合；

设置所述全过程热力耦合模型的基础设计参数，所述基础设计参数包括***相关参数及岩体热力学参数；

基于所述全过程热力耦合模型和所述基础设计参数，设置围岩岩体在开挖的隧洞边界处的受力模型；

基于所述全过程热力耦合模型和所述基础设计参数，设置围岩岩体在不同开挖过程下的过程控制方程，所述过程控制方程包括状态方程及边界条件，其中，所述状态方程包括动态应力传导方程和热传导方程，所述边界条件包括所述边界处的受力模型；

求解所述过程控制方程，获得不同过程下围岩岩体受到的动态应力和热应力；

根据求解得到的动态应力、热应力及其与初始地应力的耦合，获得围岩热力耦合模型；

通过应力投影，将所得围岩热力耦合模型投影至应力不变量空间中，获得不同所述变量参数下的、围岩开挖全过程中的围岩岩体应力轨迹曲线；

根据所述应力轨迹曲线获得不同变量参数下围岩岩体开挖损伤的分布范围、形成阶段及形成原因。

其中，***相关参数可进一步包括如：***性质，炮孔布设和装药方式等；岩体热力学参数可进一步包括如：岩石密度ρ_r，杨氏弹性模量E、岩石泊松比υ，不同围压条件下的岩石屈服应力，岩石热膨胀系数β，热传导系数κ和比热容C等。

其中，所述形成阶段(或实施方式所述形成机理)是指的隧洞开挖全过程的各阶段，如后文所述的开挖***、地应力卸荷或卸荷后降温这些阶段。

其中，所述形成原因(或实施方式所述力学机制)是指的产生围岩岩体损失的热力学原因，如通过张拉机制形成、剪胀机制形成或剪缩机制形成。

根据本发明的一些优选实施方式，所述耦合的形式为：由地应力叠加动态应力叠加热应力获得耦合力。

根据本发明的一些优选实施方式，所述隧洞开挖全过程包括开挖***、地应力卸荷及卸荷后降温这三个阶段过程；所述变量参数包括隧洞开挖半径a、卸荷时间t_du、地应力P₀和原岩温度T₀。

在该优选实施方式下，所述隧洞边界处的受力模型可采用如均布***与卸载波压力函数F(t)，其由***设计参数、地应力量级和卸荷时间综合确定。

根据本发明的一些优选实施方式，所述全过程热力耦合模型的建立包括以下设置条件：

所述隧洞为通过全断面钻爆法在无限大的弹性均质岩体中开挖的隧洞，其轴线足够长、断面为半径为a的圆形，其轮廓面为在平面应变条件下一次开挖成型；

所述弹性均质岩体的初始温度为原岩温度T₀，且受到远场静水压力即所述地应力P₀的作用；

开挖过程中，当***波压力衰减至与开挖边界的初始地应力相同时，产生地应力的瞬态卸荷，通过在洞壁上沿线性路径施加反向牵引力模拟地应力P₀在卸荷时间t_du下激发的岩体动态响应；

在硐室轮廓即断面半径为a的圆形隧洞的轮廓形成后，其受到冷空气降温，降温至指定时间时，完成全部开挖过程。

本领域技术人员可以理解的是，尽管该优选实施方式中存在表面上无法实际人为操作的过程，如无限大的弹性均质岩体，轴线足够长的隧洞等，但其在本领域的模型设置中具有可操作性，如基于无限大的弹性均质岩体的设置条件，本领域技术人员可以理解到，在模型的设置或求解中，可采用r→∞作为一个边界条件。即，该优选实施方式中的表述在本领域的技术范畴下，是可以具体实施的。

根据本发明的一些优选实施方式，所述基础设计参数包括：***密度、爆轰波速、爆轰产物等熵指、爆轰产物绝热常数、装药直径、炮孔直径、炮孔间距、装药长度、岩石密度、岩石弹性模量、岩石泊松比、岩石热膨胀系数、岩石热传导系数、岩石比热容和岩石换热系数。

根据本发明的一些优选实施方式，所述求解通过拉普拉斯变换法实现。

更具体的，所述拉普拉斯变化法既可进一步应用于后文所述动应力传导方程的求解中，获得开挖边界周围岩体的动态应力演化特征；也可应用于后文所述热传导方程的求解中，获得隧洞通风降温过程中的岩体热应力分布。

根据本发明的一些优选实施方式，所述隧洞边界处的受力模型包括：

其中，

t_r＝L/V_d，t_d＝2πt_r/ln(2)，F_b0＝d_bF_e0/S；

其中，t表示时间，F_b0表示作用在隧洞轮廓上的等效峰值压力，d_b表示炮孔直径，S表示相邻炮孔间距，F_e0表示炮孔壁峰值压力，δ表示***应力波衰减系数，t_r表示***荷载达到峰值的时间，L表示装药长度，V_d表示爆轰波速，t_d表示***波持续时间，t_u表示卸荷结束时间(即地应力卸载结束时间)，t_du表示卸荷时间，P₀表示地应力。

根据本发明的一些优选实施方式，所述炮孔壁峰值压力通过下式获得：

其中，ρ_e表示***密度，d_e表示装药直径，γ和η分别表示爆轰产物的等熵指数和绝热常数。

根据本发明的一些优选实施方式，所述过程控制方程包括：

其中，r表示柱坐标系下围岩岩体质点距硐室中心的距离，

表示柱坐标系下围岩岩体质点的位移势函数，v_p表示动态应力波在围岩岩体介质中进行传播时的纵波波速，

表示柱坐标系下围岩岩体受到的径向动态应力函数。

根据本发明的一些优选实施方式，所述过程控制方程还包括：

其中，T(r,t)表示岩体温度重分布函数，ρ_r表示岩石密度，C表示岩石比热容，κ表示岩石热传导系数，H＝h/κ、其中h表示岩石换热系数，T(a,t)表示根据岩体温度重分布函数T(r,t)获得的洞壁单元温度，t_u表示卸荷结束时间，T_b表示冷却气体温度。

根据本发明的一些优选实施方式，所述围岩热力耦合模型包括：

其中，σ_r、σ_θ和σ_z分别表示柱坐标系下围岩岩体受到的径向总耦合应力、环向总耦合应力和平面外法线方向总耦合应力；τ_rθ表示耦合总剪应力；

和

分别表示柱坐标系下围岩岩体受到的动态径向应力、环向应力和平面外法线方向应力、

表示动态剪应力，其通过求解所述动态应力传导方程获得；

和

分别表示柱坐标系下围岩岩体受到的径向热应力、环向热应力和平面外法线方向热应力、

表示由温度变化(如降温)产生的剪应力，其通过求解所述热传导方程获得。

根据本发明的一些优选实施方式，

及

通过以下过程获得：

根据式(3)所述的传导方程，利用广义胡克定律获得如下的动应力分量：

其中，λ和μ表示拉梅常数，其可由岩石弹性模量E和泊松比υ计算获得；

其中，位移势函数的求解可通过：

通过拉普拉斯变换法对方程(3)进行求解并代入各边值条件得到复变域下的岩体质点动态位移势函数，如下：

其中，k_d＝p/v_p，v_p表示所述纵波波速，p表示拉普拉斯算子，

表示所述隧洞边界处的受力模型F(t)的拉普拉斯变换结果，K₀和K₁分别表示零阶和一阶的第二类修正贝塞尔函数的导数；

通过Durbin提出的数值逆变换方法进一步求解，得到：

其中，α表示0到Re(p)之间的任意实数，Re()表示取实部，k表示计算次数，N表示最大计算次数，i表示虚数单位，Im()表示取虚部，T_i表示求解时间间隔、且0≤t≤T_i/2。

根据本发明的一些优选实施方式，

和

表示通过以下过程获得：

通过拉普拉斯变换法对方程(7)进行求解，得到隧洞通风过程中的岩体温度重分布如下：

其中，u表示积分变量，T_b表示冷却气体温度，T₀表示原岩温度，J₀、J₁和Y₀、Y₁分别表示第一类和第二类的零阶和一阶贝塞尔函数。

则，各热应力分量可由下式获得：

其中，β表示岩石的线性热膨胀系数。

根据本发明的一些优选实施方式，所述应力投影的过程包括：

a)将柱坐标系下的应力分量转换为笛卡尔坐标系，如下：

σ_x＝σ_rcos²θ+σ_θsin²θ-2τ_rθsinθcosθ

σ_y＝σ_rsin²θ+σ_θcos²θ+2τ_rθsinθcosθ

τ_xy＝(σ_r-σ_θ)sinθcosθ+τ_rθ(cos²θ-sin²θ) (11)；

b)通过下式获得模型中岩石单元的静水压力I₁和主应力差J₂：

C)以I₁为横坐标，

为纵坐标，绘制不同变量参数及基础设计参数下的岩石单元在所述全过程中的应力轨迹。

根据本发明的一些优选实施方式，所述分析方法包括：通过分析开挖轮廓周围岩体应力轨迹与其屈服面之间的相对位置关系，确定围岩开挖损伤的分布范围、形成阶段和形成原因。

本发明进一步提供了上述任一分析方法的一种应用，为将其用于开挖深埋高温硐室中。

根据本发明的一些优选实施方式，所述应用包括：根据不同变量参数下确定的围岩开挖损伤的分布范围、形成阶段和形成原因，选择合适的变量参数组成优选的开挖方案。

本发明具备以下有益效果：

(1)本发明的分析方法建立了如深埋高温硐室等情况下的隧洞开挖***-卸荷-降温全过程的分析模型，并在一些具体实施方式中，通过拉普拉斯变换法和应力叠加求解了该模型在不同地质条件与开挖方式下的岩石热-力耦合响应，可在实际开挖前获得准确的分析、评估方案；

(2)本发明的分析方法通过应力投影技术将不同时间尺度下的复杂岩体力学行为转换到应力不变量空间中，进而可在具体实施中根据应力轨迹与岩石屈服面的相对位置关系直观判断开挖损伤的影响范围、形成机理(如***、卸荷、降温)和力学机制(如张拉，剪胀、剪缩)；

(3)相较于数值计算或实际监测方法，本发明的分析方法不仅可确定开挖损伤区的分布范围，还能揭示岩体损伤的形成机理和力学机制，同时具有计算效率高、所需时间短和耗费人力物力小等优点，可为高温高压条件下的硐室适应性开挖设计提供预测与指导。

附图说明

图1为本发明的一种深埋高温硐室开挖过程的热力耦合模型示意图。

图2为实施例1所得不同开挖尺寸下的围岩应力轨迹演化图。

图3为实施例1所得不同开挖尺寸下的围岩损伤区分布特征。

图4为实施例1所得不同卸荷速率下的围岩应力轨迹演化图。

图5为实施例1所得不同卸荷速率下的围岩损伤区分布特征图。

图6为实施例1所得不同地应力条件下的围岩应力轨迹演化图。

图7为实施例1所得不同地应力条件下的围岩损伤区分布特征图。

图8为实施例1所得不同地应力下的开挖损伤随原岩温度的变化规律图。

具体实施方式

以下结合实施例和附图对本发明进行详细描述，但需要理解的是，所述实施例和附图仅用于对本发明进行示例性的描述，而并不能对本发明的保护范围构成任何限制。所有包含在本发明的发明宗旨范围内的合理的变换和组合均落入本发明的保护范围。

根据本发明的技术方案，一些具体的围岩开挖损伤识别方法，包括以下步骤：

(1)建立挖筑深埋高温硐室的全过程模型，该挖筑全过程包括硐室的开挖***、地应力卸荷及卸荷后降温，即获得深埋高温硐室的开挖***-卸荷-降温全过程模型；

更具体的，可建立一种如附图1所示的全过程模型，其包括以下变量参数：硐室开挖半径a、卸荷时间t_du、地应力P₀和原岩温度T₀，且基于上述变量参数的模型设置如下：

所述硐室为通过全断面钻爆法在无限大的弹性均质岩体中开挖一轴线足够长、断面为半径为a的圆形的隧洞获得，该隧洞轮廓面视为在平面应变条件下一次开挖成型；

该弹性均质岩体的初始温度为T₀，且受到远场静水压力即地应力P₀的作用；

该开挖过程中，当***波压力衰减至与开挖边界的初始地应力相同时，产生地应力的瞬态卸荷，通过在洞壁上沿线性路径施加反向牵引力模拟地应力P₀在卸荷时间t_du下激发的岩体动态响应；

在硐室轮廓即断面半径为a的圆形隧洞的轮廓形成后，其受到冷空气降温，降温至指定时间时，完成全部开挖过程。

在上述设置下，附图1所示的全过程模型，可具体表现为：如图1(a)所述的深埋高温硐室围岩在***-卸荷和降温过程中的热-力耦合响应模型为在如图1(b)所述的初始地应力模型的基础上叠加由图1(c)所示的***和卸载波引起的动态力学响应模型以及由图1(d)所示的硐室通风降温造成的围岩二次应力调整模型。

在上述设置下，本发明进一步包括：

(2)获得所述全过程模型所需的基础设计参数；

更具体的，在附图1所示的全过程模型中，其所需的基础设计参数包括：

***密度、爆轰波速、爆轰产物等熵指数、爆轰产物绝热常数、装药直径、炮孔直径、炮孔间距、装药长度、岩石密度、岩石弹性模量、岩石泊松比、岩石热膨胀系数、岩石热传导系数、岩石比热容和岩石换热系数。

(3)设置所述全过程模型的边界受力模型；

更具体的，在附图1所述的全过程模型中，硐室开挖边界由以下均布***-卸载波压力函数F(t)确定：

其中，

t_r＝L/V_d，t_d＝2πt_r/ln(2)，F_b0＝d_bF_e0/S；

其中，t表示时间，F_b0表示作用在隧洞轮廓上的等效峰值压力，d_b表示炮孔直径，S表示相邻炮孔间距，F_e0表示炮孔壁峰值压力，δ表示***应力波衰减系数，t_r表示***荷载达到峰值的时间，L表示装药长度，V_d表示爆轰波速，t_d表示***波持续时间，t_u表示卸荷结束时间(即地应力卸载结束时间)，t_du表示卸荷时间，P₀表示地应力。

其中，炮孔壁峰值压力F_e0可进一步通过如下的C-J爆轰理论模型计算得到：

其中，ρ_e表示***密度，d_e表示装药直径，γ和η分别表示爆轰产物的等熵指数和绝热常数。

(4)设置所述全过程模型在不同时间段的过程控制方程，所述过程控制方程包括：在开挖***-卸荷过程中动态应力波在围岩中的传播方程，在隧洞通风降温过程中描述围岩温度变化的热传导方程；

更具体的，所述动态应力波在围岩中的传播模型如图1(c)所示，其传播方程设置如下：

其中，r表示柱坐标系下围岩岩体质点距硐室中心的距离，

表示柱坐标系下围岩岩体质点的位移势函数，v_p表示动态应力波在围岩岩体介质中进行传播时的纵波波速，

表示柱坐标系下围岩岩体受到的径向动态应力函数。

更具体的，所述隧洞通风降温过程中的热传导模型如图1(d)所示，其传导方程设置如下：

其中，T(r,t)表示岩体温度重分布函数，ρ_r表示岩石密度，C表示岩石比热容，κ表示岩石热传导系数，T(a,t)表示洞壁单元的温度，t_u表示地应力卸载结束时间，T_b表示冷却气体温度，H＝h/κ表示计算系数，其中h表示岩石换热系数。

(5)求解所述过程控制方程，根据动态应力波在围岩中的传播方程获得不同时间下的围岩岩体受到的动态应力，根据隧洞通风降温过程中的热传导方程获得不同时间下的围岩岩体受到的热应力；

更具体的，对所述传播方程的求解可包括：

根据所述传播方程，利用广义胡克定律获得开挖过程中的岩体动应力分量如下：

其中，

和

分别表示柱坐标系下围岩岩体受到的动态径向应力、环向应力和平面外法线方向应力；λ和μ表示拉梅常数，其可由岩体弹性模量E和泊松比υ计算获得；

表示动态剪应力；

通过拉普拉斯变换法对控制方程(3)进行求解并代入各边值条件得到复变域下的岩体质点动态位移势函数，如下：

其中，k_d＝p/v_p，v_p表示所述纵波波速，p表示拉普拉斯算子，

表示所述受力模型F(t)的拉普拉斯变换结果，K₀和K₁分别表示零阶和一阶的第二类修正贝塞尔函数的导数；

通过Durbin提出的数值逆变换方法进一步求解，得到：

其中，α表示0到Re(p)之间的任意实数，Re()表示取实部，k表示计算次数，N表示最大计算次数，i表示虚数单位，Im()表示取虚部，T_i表示求解时间间隔、且0≤t≤T_i/2。

更具体的，对所述热传导方程的求解可包括：

采用拉普拉斯变换法对热传导控制方程(7)进行求解，得到隧洞通风过程中的岩体温度重分布如下：

其中，u表示积分变量，J₀、J₁和Y₀、Y₁分别表示第一类和第二类的0阶和1阶贝塞尔函数。

则，通风降温过程在保留岩体中产生的热应力可由下式获得：

其中，β表示岩石的线性热膨胀系数，

和

分别表示柱坐标系下围岩岩体受到的径向热应力、环向热应力和平面外法线方向热应力；

表示由降温产生的剪应力；

(5)通过所述围岩岩体受到的动态应力、所述围岩岩体受到的热应力及初始地应力获得围岩热力耦合响应模型；

更具体的，所述热力耦合响应模型设置如下：

其中，σ_r、σ_θ和σ_z分别表示柱坐标系下围岩岩体受到的径向总耦合应力、环向总耦合应力和平面外法线方向总耦合应力；τ_rθ表示耦合总剪应力。

(6)通过应力投影，将所得围岩热力耦合响应模型投影至应力不变量空间中，进而获得含有表征不同地质条件和开挖方式的变量参数及基础设计参数的***、卸荷和降温全过程中的岩石应力轨迹曲线；

更具体的，其可包括

a)将柱坐标系下的应力分量转换为笛卡尔坐标系，如下：

σ_x＝σ_rcos²θ+σ_θsin²θ-2τ_rθsinθcosθ

σ_y＝σ_rsin²θ+σ_θcos²θ+2τ_rθsinθcosθ

τ_xy＝(σ_r—σ_θ)sinθcosθ+τ_rθ(cos²θ-sin²θ) (11)；

b)通过下式获得模型中岩石单元的静水压力I₁和主应力差J₂：

C)以I₁为横坐标，

为纵坐标，绘制不同变量参数及基础设计参数下的岩石单元在***-卸荷-降温全过程中的应力轨迹即可。

(7)通过所得岩石应力轨迹曲线分析开挖轮廓周围岩石应力轨迹与其屈服面之间的相对位置关系，获得所述变量参数与岩体开挖损伤的分布范围、形成机理和力学机制之间的关系，获得不同地质条件和开挖方式下的岩石热力耦合响应和围岩损伤区时空演化规律，进而确定实际开挖方案。

其中，所述形成机理包括***、卸荷或降温，所述力学机制包括张拉、剪胀或剪缩。

当***、卸荷或降温应力轨迹与屈服面相交时，认为围岩分别产生***、卸荷和降温损伤。同时，根据应力轨迹与屈服面的交点位置可判断围岩损伤的力学机制，如：

a.当交点位于拉屈服段，围岩产生张拉损伤；

b.当交点位于剪胀屈服段，围岩产生剪胀损伤；

c.当交点位于帽盖段，围岩产生剪缩损伤。

实施例1

选择深埋高温硐室周围的岩石主要为印第安纳灰岩的情况通过以上具体实施方式所述的具体方案进行仿真实验。

其中，开挖设计方案如下表1所示，且冷却降温过程使用温度T_a为15℃的气流，通风时间为15天：

表1模型方案

基础设计参数如下表2所示：

表2开挖***-卸荷-降温模型基础设计参数

数值逆变换求解过程中，取α＝0.2，T_i＝40ms，迭代计算次数N＝5000。

在以上条件下，可得到如附图2-8所示的损伤分析图，其中：

图2为不同开挖半径参数a下的隧洞边界单元围岩岩体在***、卸荷和降温全过程中的应力轨迹演化，其中实心圆表示不同应力轨迹与初始屈服面的交点，即代表岩石开始形成损伤的点，其对应的开挖***-卸荷-降温造成的围岩损伤分布特征如附图3所示。

图4为不同卸荷速率(卸荷时间t_du)下的隧洞边界单元围岩岩体应力轨迹演化，同样的，其中实心圆表示不同应力轨迹与初始屈服面的交点，即代表岩石开始形成损伤的点，其相应的围岩开挖损伤分布特征如附图5所示。

图6为不同地应力下的下的隧洞边界单元围岩岩体应力轨迹演化，同样的，其中实心圆表示不同应力轨迹与初始屈服面的交点，即代表岩石开始形成损伤的点，，其相应的围岩开挖损伤分布特征如附图7所示。

图8为不同地应力下围岩损伤随原岩温度的变化规律统计图，即其热损伤发展规律图，图中所示温度数值为不同的原岩温度，所述时间数值为隧洞边界单元围岩岩体产生热损伤的时间。

从以上各图中可以看出：

在不同硐室开挖半径参数a下：

由图2可知，在***荷载作用下，开挖半径较小时的隧洞边界单元应力轨迹将超出剪切屈服面，并最终在保留岩体内造成***剪胀损伤(图3(a))。当开挖半径增大时，***加载产生的应力轨迹变得更加扁平进而不与岩石屈服面发生相交，因此***荷载无法造成保留岩体的开挖损伤(图3(b)和(c))。而在地应力瞬态卸荷过程中，较大的开挖半径会导致卸荷应力轨迹更进一步地延伸至剪切屈服面外(图2)，从而在保留岩体中造成更强的卸荷剪胀损伤(图3)。在随后的通风降温过程中，围岩的热损伤范围随着硐室半径的增大而逐渐扩大(图3)。

在不同的卸荷速率(卸荷时间t_du)下：

由图4-5可知，随着卸荷速率增大，应力轨迹向拉屈服面运动的趋势更加明显。当卸荷速率足够快时，岩体中将由于地应力的急剧释放而引发动态拉应力，从而可能导致卸荷拉损伤的形成。径向应力卸荷完成后，由于惯性作用的影响，岩体中将产生附加的环向动态应力，且卸荷速率越快动应力幅值越大，进而使得应力轨迹超出剪切屈服面更多，最终在保留岩体中造成更强烈的卸荷剪胀损伤(图5)。

在不同地应力地应力P₀下：

由图6和7可知，随着地应力量级增大，***应力轨迹将随地应力升高造成的初始应力状态变化而逐渐在应力不变量空间中向右上方平移并不断远离剪切屈服面，进而抑制***损伤的形成。但在地应力释放过程中，围压量级越大会引发更高的环向附加动应力，进而产生更强的卸载损伤。此外，由于地应力可抑制降温造成的岩体压力松弛，因此热损伤范围随着地应力增大呈现单调减小的趋势。

在不同的地热条件(即不同地应力和原岩温度T₀)下：

由图8可知，在给定的地应力条件下，具有更高温度的围岩在相同冷却介质作用下会产生更大范围的热损伤，而地应力量级增大则会在相同地热环境下抑制围岩热损伤的发展。此外，在相同的冷却气流作用下，围岩的初始温度越高则热交换速率越快，因此二次应力调整更为剧烈，进而导致围岩热损伤出现的时机更早(图8)。同时，由于围压的上升抵消了隧洞降温造成的岩体压力松弛，因此在相同的地热条件下，热损伤发生时间随地应力量级的增大而显著延缓(图8)。

综上可以看出，本发明的分析方法有效揭示了不同地压和地温条件以及开挖方式下的围岩损伤形成机理、力学机制和演化特征，可为深部高温高压条件下的硐室开挖设计提供重要的实践辅助。

以上实施例仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例。凡属于本发明思路下的技术方案均属于本发明的保护范围。应该指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下的改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (8)

1.热力耦合条件下的围岩开挖损伤分析方法，其特征在于：包括：

建立包含变量参数的隧洞开挖全过程热力耦合模型，所述热力耦合模型包括地应力、动态应力及热应力的耦合；

设置所述全过程热力耦合模型的基础设计参数，所述基础设计参数包括***相关参数及岩体热力学参数；

基于所述全过程热力耦合模型和所述基础设计参数，设置围岩岩体在隧洞边界处的受力模型；

基于所述全过程热力耦合模型和所述基础设计参数，设置围岩岩体在不同开挖过程下的过程控制方程，所述过程控制方程包括状态方程及边界条件，其中，所述状态方程包括动态应力传导方程和热传导方程，所述边界条件包括所述边界处的受力模型；

求解所述过程控制方程，获得不同过程下围岩岩体受到的动态应力和热应力；

根据求解得到的动态应力、热应力及其与初始地应力的耦合，获得围岩热力耦合模型；

通过应力投影，将所得围岩热力耦合模型投影至应力不变量空间中，获得不同所述变量参数下的、围岩开挖全过程中的围岩岩体应力轨迹曲线；

根据所述应力轨迹曲线获得不同变量参数下围岩岩体开挖损伤的分布范围、形成阶段及形成原因；

其中，所述形成原因包括张拉形成、剪胀形成或剪缩形成；

所述隧洞边界处的受力模型包括：

其中，

t_r＝L/V_d，t_d＝2πt_r/ln(2)，F_b0＝d_bF_e0/S；

其中，t表示时间，F_b0表示作用在隧洞轮廓上的等效峰值压力，d_b表示炮孔直径，S表示相邻炮孔间距，F_e0表示炮孔壁峰值压力，δ表示***应力波衰减系数，t_r表示***荷载达到峰值的时间，L表示装药长度，V_d表示爆轰波速，t_d表示***波持续时间，t_u表示卸荷结束时间，t_du表示卸荷时间，P₀表示地应力；

所述围岩热力耦合模型包括：

其中，σ_r、σ_θ和σ_z分别表示柱坐标系下围岩岩体受到的径向总耦合应力、环向总耦合应力和平面外法线方向总耦合应力；τ_rθ表示耦合总剪应力；

和

分别表示柱坐标系下围岩岩体受到的动态径向应力、环向应力和平面外法线方向应力、

表示动态剪应力，其通过求解所述动态应力传导方程获得；

和

分别表示柱坐标系下围岩岩体受到的径向热应力、环向热应力和平面外法线方向热应力、

表示由温度变化产生的剪应力，其通过求解所述热传导方程获得。

2.根据权利要求1所述的分析方法，其特征在于：所述隧洞开挖全过程包括开挖***、地应力卸荷及卸荷后降温这三个阶段过程；所述变量参数包括隧洞开挖半径a、卸荷时间t_du、地应力P₀和原岩温度T₀。

3.根据权利要求2所述的分析方法，其特征在于：所述全过程热力耦合模型的建立包括以下设置条件：

所述隧洞为通过全断面钻爆法在无限大的弹性均质岩体中开挖的隧洞，其轴线足够长、断面为半径为a的圆形，其轮廓面为在平面应变条件下一次开挖成型；

所述弹性均质岩体的初始温度为原岩温度T₀，且受到远场静水压力即所述地应力P₀的作用；

开挖过程中，当***波压力衰减至与开挖边界的初始地应力相同时，产生地应力的瞬态卸荷，通过在洞壁上沿线性路径施加反向牵引力模拟地应力P₀在卸荷时间t_du下激发的岩体动态响应；

在硐室轮廓即断面半径为a的圆形隧洞的轮廓形成后，其受到冷空气降温，降温至指定时间时，完成全部开挖过程。

4.根据权利要求1所述的分析方法，其特征在于：所述基础设计参数包括：***密度、爆轰波速、爆轰产物等熵指、爆轰产物绝热常数、装药直径、炮孔直径、炮孔间距、装药长度、岩石密度、岩石弹性模量、岩石泊松比、岩石热膨胀系数、岩石热传导系数、岩石比热容和岩石换热系数。

5.根据权利要求1所述的分析方法，其特征在于：所述求解通过拉普拉斯变换法实现。

6.根据权利要求1所述的分析方法，其特征在于：所述过程控制方程包括：

其中，r表示柱坐标系下围岩岩体质点距硐室中心的距离，

表示柱坐标系下围岩岩体质点的位移势函数，v_p表示动态应力波在围岩岩体介质中进行传播时的纵波波速，

表示柱坐标系下围岩岩体受到的径向动态应力函数；

和/或：

其中，T(r,t)表示岩体温度重分布函数，ρ_r表示岩石密度，C表示岩石比热容，κ表示岩石热传导系数，H＝h/κ、其中h表示岩石换热系数，T(a,t)表示根据岩体温度重分布函数T(r,t)获得的洞壁单元温度，t_u表示卸荷结束时间，T_b表示冷却气体温度。

7.根据权利要求1所述的分析方法，其特征在于：所述应力投影的过程包括：

a)将柱坐标系下的应力分量转换为笛卡尔坐标系，如下：

σ_x＝σ_rcos²θ+σ_θsin²θ-2τ_rθsinθcosθ

σ_y＝σ_rsin²θ+σ_θcos²θ+2τ_rθsinθcosθ

τ_xy＝(σ_r-σ_θ)sinθcosθ+τ_rθ(cos²θ-sin²θ) (11)；

b)通过下式获得模型中岩石单元的静水压力I₁和主应力差J₂：

C)以I₁为横坐标，

为纵坐标，绘制不同变量参数及基础设计参数下的岩石单元在所述开挖全过程中的应力轨迹。

8.权利要求1-7中任一项所述的分析方法的一种应用方法，其为将所述分析方法应用于开挖深埋高温硐室上。

Images