CN119249552B - 考虑热水力全耦合效应的深部圆形洞室开挖响应预测方法 - Google Patents

Abstract

本发明提供了一种考虑热水力全耦合效应的深部圆形洞室开挖响应预测方法，包括假定弹塑性解界面处的状态；得到弹塑性界面处的径向和切向有效应力；由弹塑性界面处的应力状态向隧道方向逐层差分获得第i层圆环计算单元的应力，并得到该层圆环计算单元的应变；更新得到第i+1层圆环计算单元损伤变量等计算参数；若r不大于R₀，向下一层圆环计算单元继续重复以上计算过程；若r大于等于R₀，检查开挖边界处有效应力、孔隙水压和温度的差分计算值和设定值之间的误差，若误差过大，按梯度下降算法更新弹塑性界面处状态假定值，持续反演直至误差符合预期，获得最终解。本发明解决了目前深部热水力耦合环境洞室开挖设计缺乏有效计算模型的问题。

Description

考虑热水力全耦合效应的深部圆形洞室开挖响应预测方法

技术领域

本发明属于地下洞室施工技术领域，尤其涉及一种考虑热水力全耦合效应的深部圆形洞室开挖响应预测方法。

背景技术

地下空间的使用日渐进入深部空间，面临热水力耦合环境。热水力耦合环境下，洞室开挖引起的卸荷效应与非耦合环境不同，总的来说热水力耦合环境会加剧围岩的卸荷效应，对工程的直接影响就是收敛变形大幅上升，围岩破裂深度更大。原本来讲，洞室开挖后的收敛位移和破裂区半径等的预测对于地下工程的建设就尤为重要，是指导洞室支护的关键参数，而目前对于热水力耦合环境下洞室开挖响应的评估仍主要依赖数值模拟，缺乏有效的理论分析模型。不可否认，数值模拟是地下工程建设的有力工具，然而，解决大规模工程问题，特别是在热水力耦合环境下，由于问题的非线性，往往需要很长的计算时间。对于这类复杂的问题，通常采用建立简化模型并解析求解的方法。事实上，初步设计和大部分设计都是基于相对简单的分析模型。因此，建立一个能够有效反映热水力耦合环境下地下洞室开挖分析模型于工程而言是十分必要和迫切的。

在深部地下工程建设中，经常遇到热水力耦合问题。在应力、孔隙水压力和岩石温度的共同作用下，不仅岩体的物理性质发生了相对于常规地质的变化，而且地下开孔周围围岩的应力调节过程也与仅在应力条件下有明显的不同。准确预测洞室开挖后围岩的稳定性对地下工程施工至关重要，是保证施工人员安全及后续支护设计的前提。尽管众多学者对隧道力学解析解进行了延伸，但仍多在非耦合场的背景下。然而，随着地下工程向深部发展，围岩多处于耦合场环境，此类工程建设经验表明，耦合环境下围岩卸荷效应将加剧。为此，一些学者通过各场域之间的简单叠加进行近似求解，然而，水压、温度对围岩卸荷效应的影响不仅仅是简单的应力、渗透力、热应力的叠加，而是三场之间的实时动态耦合，即应力、水压、温度之间的相互作用。一般而言，可以通过建立岩石损伤与各物理量之间的动态变化关系考虑这种动态耦合，研究表明热水力耦合场作用下的围岩行为预测考虑围岩损伤与各场之间的耦合作用是十分必要的。但由于各场之间的相互作用，此列问题的求解十分困难，因此，目前关于热水力耦合环境下的围岩稳定性分析多利用数值模拟进行。

显然，数值模拟是解决地下工程设计和施工问题的有力工具，然而数值软件求解类似于热水力耦合这类强烈非线性问题的工程大尺度问题通常需要较长的运行时间。因此，对于此类复杂工程问题，通常可以建立简化模型并进行解析解求解，进而为之后的支护初步设计提供指导。

发明内容

针对现有技术中的上述不足，本发明提供的一种考虑热水力全耦合效应的深部圆形洞室开挖响应预测方法，解决了目前深部热水力耦合环境洞室开挖设计缺乏有效计算模型的问题。

为了达到以上目的，本发明采用的技术方案为：一种考虑热水力全耦合效应的深部圆形洞室开挖响应预测方法，包括以下步骤：

S1、在深部圆形洞室的预设范围内假定弹塑性解界面处的状态；

S2、计算获得弹性计算区域的水压分布函数系数项C1和温度分布函数系数项C3；

S3、根据假定的弹塑性解界面处的状态，计算得到弹塑性界面处的径向有效应力与切向有效应力的和；

S4、根据径向有效应力与切向有效应力的和以及屈服准则，分别计算得到弹塑性界面处的径向有效应力和切向有效应力；

S5、将塑性区域离散化，获得塑性区若干个等间距的圆环计算单元，共n+1个圆环计算节点，将S4的计算结果作为第1个圆环计算节点状态值；

S6、根据第i层圆环计算单元的状态值，计算第i+1层圆环计算单元的状态值；

S7、判断当前计算圆环单元距离洞室中心的距离是否大于开挖半径，若是，则将第i+1层圆环计算单元视为新的第i层圆环计算单元，并返回步骤S6，否则，进入步骤S8；

S8、判断开挖边界处通过差分计算得到的孔隙水压、温度以及径向有效应力与设定条件之间的相对误差是否均满足要求，若是，则完成对深部圆形洞室开挖的响应预测，否则，返回步骤S1，并使用梯度下降法计算获得新的弹塑性界面状态假定值。

进一步地，所述径向有效应力以及切向有效应力的和的表达式如下：

其中，G(R_p)表示弹塑性界面处的径向有效应力以及切向有效应力的和，σ'₀表示初始地应力，R_e表示边界距离洞室中心的距离，R_p表示弹塑性界面距离洞室中心的距离，E表示围岩的弹性模量，α表示围岩单元的线弹性密度，C₁和C₃分别表示弹性区域围岩中温度分布系数项和孔隙水压分布系数项，v表示泊松比，T_p和T₀分别表示弹塑性界面处的地温和初始地温，P_w,p和P_w,0分别表示弹塑性界面处的孔隙水压和初始孔隙水压。

再进一步地，所述屈服准则的表达式如下：

其中，F表示屈服面函数，F(σ_θ,σ_r)＝0表示围岩单元达到塑性，σ_θ和σ_r分别表示切向有效应力和径向有效应力，表示瞬时摩擦角，c表示围岩瞬时粘聚力，c_ini表示初始粘聚力，γ_p表示等效塑性剪切应变，n'表示粘聚力演变曲线常数，表示切向塑性应变，表示径向塑性应变，D表示损伤变量，c_res表示残余粘聚力，ω表示瞬时孔隙率，ω_ini表示初始孔隙率，ω_res表示残余孔隙率，k表示瞬时渗透率，k_ini表示初始渗透率，λ表示围岩的瞬时导热系数，λ_s表示岩石骨架的导热系数，λ_w表示水的导热系数。

再进一步地，所述第i层圆环计算单元的径向有效应力和切向有效应力的表达式分别如下：

其中，表示第i+1层圆环计算单元的径向有效应力，表示第i+1层圆环计算单元的切向有效应力，表示第i层圆环计算单元的径向有效应力，rⁱ⁺¹表示第i+1层圆环计算单元距离隧道中心的距离，E表示围岩的弹性模量，β表示围岩的比奥系数，α表示围岩热膨胀系数，rⁱ表示第i层圆环计算单元距离隧道中心的距离，Tⁱ⁺¹表示第i+1层圆环计算单元的温度，Tⁱ表示第i层圆环计算单元的温度，βⁱ⁺¹表示第i+1层圆环计算单元的比奥系数，表示第i+1层圆环计算单元孔隙水压，表示第i层圆环计算单元的径向有效应力，表示i层圆环计算单元的切向有效应力，表示第i层圆环计算单元的孔隙水压，Tⁱ表示第i层圆环计算单元的温度，表示第i层圆环计算单元的摩擦角，c表示围岩瞬时粘聚力；

第i+1层圆环计算单元的切向塑性应变和径向塑性应变的表达式分别如下：

其中，表示第i+1层圆环计算单元的切向塑性应变，表示第i层圆环计算单元的切向塑性应变，表示第i层圆环计算单元的径向塑性应变，表示第i层圆环计算单元的径向弹性应变，K_ψ表示塑性乘子，rⁱ⁺¹表示第i+1层圆环计算单元距离隧道中心的距离，rⁱ表示第i层圆环计算单元距离隧道中心的距离，表示第i+1层圆环计算单元的切向弹性应变，表示第i层圆环计算单元的切向弹性应变，表示第i+1层圆环计算单元的径向塑性应变，表示第i层圆环计算单元的径向塑性应变，ψ表示剪胀角。

再进一步地，所述根据第i层圆环计算单元的状态值，计算第i+1层圆环计算单元的状态值，其具体为：

基于第i层圆环计算单元孔隙水压和温度，计算第i+1层圆环计算单元的孔隙水压和温度；

结合第i层圆环计算单元应力状态、计算参数以及第i+1层圆环计算单元的孔隙水压和温度，得到第i+1层圆环计算单元的有效应力和塑性应变；

利用第i+1层圆环计算单元的塑性应变，更新第i+1层圆环计算单元的计算参数；

将原第i+1层圆环计算单元看作新的第i层圆环计算单元，继续进行新的第i+1层圆环计算单元的状态值，直至圆环计算单元距隧道中心的距离不大于隧道开挖半径，完成对第i+1层圆环计算单元的状态值的计算，其中，状态值包括有效应力、塑性应变、孔隙水压以及温度。

再进一步地，所述第i+1层圆环计算单元的孔隙水压的表达式如下：

其中，表示第i+1层圆环计算单元的孔隙水压，表示第i层圆环计算单元的孔隙水压，Q表示渗流通量，rⁱ表示第i层圆环计算单元距离隧道中心的距离，rⁱ⁺¹表示第i+1层圆环计算单元距离隧道中心的距离，kⁱ表示第i层圆环计算单元的渗透率；

所述第i+1层圆环计算单元的温度的表达式如下：

其中，Tⁱ⁺¹表示第i+1层圆环计算单元的温度，Tⁱ表示第i层圆环计算单元的温度，Φ表示热通量，λⁱ表示第i层圆环计算单元的导热系数。

再进一步地，所述新的弹塑性界面状态假定值的表达式如下：

LF(R_p,P_w,p,T_p)＝||(P_s,P_w,s,T_s)_calu-(P_s,P_w,s,T_s)_actu||

其中，(R_p,P_w,p,T)_j表示上一轮弹塑性界面状态假定解，(R_p,P_w,p,T)_j+1表示新的弹塑性界面状态假定解，η表示学习速率，表示损失函数的梯度，LF(R_p,P_w,p,T_p)表示损失函数，(P_s,P_w,s,T_s)_calu表示利用假定的(R_p,P_w,p,T)计算获得的开挖边界处的状态，(R_p,P_w,p,T)表示弹塑性界面状态假定解，(P_s,P_w,s,T_s)_actu表示根据实际工况设定的开挖边界状态参数，R_p表示塑性区半径，R_e表示弹性区半径，T表示温度，T_p表示弹塑性边界的温度，P_s表示开挖边界处的径向支护力，P_w,s表示开挖边界处的孔隙水压，T_s表示开挖边界处的温度。

本发明的有益效果：

本发明拟为热水力耦合环境地下工程建设提供一种理论分析手段。针对目前热水力耦合环境地下工程建设分析模型未能考虑热-水-力三个物理场之间的全耦合效应，本发明提供了一种计算算法，通过采用梯度下降算法反演参数中关键参数，即涌水量、热通量以及塑性区半径从而获得考虑热-水-力全耦合效应的计算结果。

附图说明

图1为本发明的方法流程图。

图2为本实施例中深部耦合环境洞室开挖响应预测弹塑性分析模型示意图。

图3为本实施例中热水力全耦合效应示意图。

图4为本实施例中塑性区离散化示意图。

具体实施方式

下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

实施例

本发明拟为热水力耦合环境地下工程建设提供一种理论分析手段。针对目前热水力耦合环境地下工程建设分析模型未能考虑热-水-力三个物理场之间的全耦合效应，本发明提供了一种计算算法，通过采用梯度下降算法反演参数中关键参数，即涌水量、热通量以及塑性区半径从而获得考虑热-水-力全耦合效应的计算结果。

如图1所示，本发明提供了一种考虑热水力全耦合效应的深部圆形洞室开挖响应预测方法，其实现方法如下：

S1、在深部圆形洞室的预设范围内假定弹塑性解界面处的状态；

S2、计算获得弹性计算区域的水压分布函数系数项C1和温度分布函数系数项C3；

S3、根据假定的弹塑性解界面处的状态，计算得到弹塑性界面处的径向有效应力与切向有效应力的和；

S4、根据径向有效应力与切向有效应力的和以及屈服准则，分别计算得到弹塑性界面处的径向有效应力和切向有效应力；

S5、将塑性区域离散化，获得塑性区若干个等间距的圆环计算单元，共n+1个圆环计算节点，将S4的计算结果作为第1个圆环计算节点状态值；

S6、根据第i层圆环计算单元的状态值，计算第i+1层圆环计算单元的状态值，其具体为：

基于第i层圆环计算单元孔隙水压和温度，计算第i+1层圆环计算单元的孔隙水压和温度；

结合第i层圆环计算单元应力状态、计算参数以及第i+1层圆环计算单元的孔隙水压和温度，得到第i+1层圆环计算单元的有效应力和塑性应变；

利用第i+1层圆环计算单元的塑性应变，更新第i+1层圆环计算单元的计算参数；

将原第i+1层圆环计算单元看作新的第i层圆环计算单元，继续进行新的第i+1层圆环计算单元的状态值，直至圆环计算单元距隧道中心的距离不大于隧道开挖半径，完成对第i+1层圆环计算单元的状态值的计算，其中，状态值包括有效应力、塑性应变、孔隙水压以及温度；

S7、判断当前计算圆环单元距离洞室中心的距离是否大于开挖半径，若是，则将第i+1层圆环计算单元视为新的第i层圆环计算单元，并返回步骤S6，否则，进入步骤S8；

S8、判断开挖边界处通过差分计算得到的孔隙水压、温度以及径向有效应力与设定条件之间的相对误差是否均满足要求，若是，则完成对深部圆形洞室开挖的响应预测，否则，返回步骤S1，并使用梯度下降法计算获得新的弹塑性界面状态假定值。

本实施例中，本发明首先建立热水力耦合环境下洞室开挖响应分析模型，如图2所示，图中，T表示温度，r表示围岩单元距离洞室中心的距离，P_w表示孔隙水压，σ'₀表示初始地应力，P_w,0表示初始孔隙水压，T₀表示初始温度，R₀表示洞室开挖半径，R_p表示塑性区半径，R_e表示模型计算边界，即弹性区半径，认为在热水力耦合环境下进行洞室开挖，围岩达到稳定状态后，从远端到洞室可划分为弹塑性区域和塑型区域。认为水压场、温度场、应力场的影响在距离洞室中心足够远处停止，此足够远距离在图中标记为R_e。围岩的强度参数(粘聚力，摩擦角)、耦合参数(渗透系数、导热系数)仅在塑性区内发生动态演变，忽略弹塑性区域内改变。

由此，建立了以上深部耦合环境下的弹塑性分析模型以预测洞室开挖造成的围岩响应。由图2可见，热水力耦合环境下，围岩中除了存在机械应力场，还存在热应力场和渗流应力场，这也正是热水力耦合环境问题的复杂性。在洞室开挖之后，非耦合环境下围岩状态的重新调整仅需要满足机械应力场的平衡。通常该问题的求解可以通过联立弹性力学中的Lame解析解和岩石的屈服准则从而获得弹塑性界面上的应力状态。而热水力耦合环境下，由于热应力和渗流应力的参与，以上的求解思路无法达成。

实际上，机械应力场与岩石的强度参数(粘聚力、摩擦角)密切相关，而热应力场和渗流应力场也分别与围岩的渗透系数和导热系数相关。而在塑性区内，粘聚力、摩擦角、渗透系数、导热系数持续演化，因而机械应力场、热应力场和渗透应力场之间动态演变，直至最终达到机械应力场、热应力场和渗透应力场各自独立平衡，且三个应力场之间满足一定的物理定律，此为热水力全耦合效应。通常可以通过损伤变量D将三个物理场之间的状态相互联系，其思想为如图3所示。

本实施例中，损伤变量D的定义因人而异，本发明依据粘聚力的衰退程度建立方程，其依据是，岩石破坏后强度的衰退主要是粘聚力的消散引起的，摩擦角的减小可以忽略不计。于是，整体上可以获得：

其中，F表示屈服面函数，F(σ_θ,σ_r)＝0表示围岩单元达到塑性，σ_θ和σ_r分别表示切向有效应力和径向有效应力，表示瞬时摩擦角，c表示围岩瞬时粘聚力，c_ini表示初始粘聚力，γ_p表示等效塑性剪切应变，n'表示粘聚力演变曲线常数，表示切向塑性应变，表示径向塑性应变，D表示损伤变量，c_res表示残余粘聚力，ω表示瞬时孔隙率，ω_ini表示初始孔隙率，ω_res表示残余孔隙率，k表示瞬时渗透率，k_ini表示初始渗透率，λ表示围岩的瞬时导热系数，λ_s表示岩石骨架的导热系数，λ_w表示水的导热系数。

通过在图2计算模型中嵌入以上方程以考虑热水力全耦合效应。此外，因在塑性区域内围岩发生力学参数、耦合参数的动态演变，需要将计算模型塑性区域进行离散化，如图4所示。于是，开始借助弹塑性理论进行模型求解。具体求解过程分为弹塑性区和塑性区两个部分：

1)弹塑性区

围岩的应力平衡方程为：

其中，σ'_r表示径向有效应力，σ'_θ表示径向有效应力，r表示围岩单元距离洞室中心的距离，E表示围岩的弹性模量，ΔT表示温度的该变量，ΔT＝T-T₀,T表示瞬时温度，T₀表示初始温度，β表示比奥模量，P_w表示孔隙水压，α表示岩石的热膨胀系数。

而温度场和孔隙压力场的分布形式为：

其中，C₁和C₂分别表示弹性区温度场分布函数的系数项和常数项，C₃和C₄分别表示孔隙水压场分布函数的系数项和常数项。

根据边界条件和可以解得：

以及：

其中，T_p和T₀分别表示弹塑性边界的温度和初始温度，P_w,p和P_w,0分别表示弹塑性边界和模型计算边界处的孔隙水压，R_p表示塑性区半径，即弹塑性边界到洞室中心的距离，R_e表示模型计算边界半径，即模型计算边界到洞室中心距离。

将式(3)到(5)代入式(2)，从而获得：

其中，σ'_r表示径向有效应力，σ'_θ表示径向有效应力，E表示弹性模量，α表示热膨胀系数，β表示比奥系数，r表示围岩单元距离洞室中心的距离。

另一方面，平面应变问题中应变相容方程为：

其中，ε_θ表示切向塑性应变，ε_r表示径向塑性应变，r表示围岩单元距离洞室中心的距离。

应变与应力之间满足虎克定律：

其中，v为泊松比。将式8代入式7中，可以获得：

将式(9)代入式(2)，可以获得：

求解以上方程，可以获得：

其中，σ'_r表示径向有效应力，σ'_θ表示径向有效应力，C₅和C₆均表示应力分布系数项。

将式(11)代入式(6)中获得：

认为r＝R_e处，开挖扰动的影响近似为0。于是代入弹塑性区边界条件于是：

为了确定弹塑性区围岩的位移，需要使用单元几何方程、胡克方程和应力平衡方程。岩石单元的几何方程为：

其中，ε_r表示径向应变，ε_θ表示切向应变，u表示径向位移。r表示围岩单元距离洞室中心的距离。

将式(14)和式(8)代入式(2)中，可以获得：

求解得到：

其中，u表示径向位移，C₇和C₈均表示弹性区位移分布函数系数项。

2)塑性区

进入塑性区之后，由于在距离隧道不同距离的空间的岩体的力学参数、耦合参数都发生持续不断地演变，因此，无法给出封闭的解。因此，将塑性区范围离散成为n个圆环单元，当n足够大的时候，认为每一环中的力学参数和耦合参数不发生改变，如图4所示。由此，塑性区范围内，第i+1层圆环计算单元的孔隙水压和岩温可通过i层圆环计算单元的空隙水压和岩温分别通过达西定律以及傅里叶定律积分得到：

其中，表示第i+1层圆环计算单元的孔隙水压，表示第i层圆环计算单元的孔隙水压，Q表示渗流通量，rⁱ表示第i层圆环计算单元距离隧道中心的距离，rⁱ⁺¹表示第i+1层圆环计算单元距离隧道中心的距离，kⁱ表示第i层圆环计算单元的渗透率，Tⁱ⁺¹表示第i+1层圆环计算单元的温度，Tⁱ表示第i层圆环计算单元的温度，Φ表示热通量，λⁱ表示第i层圆环计算单元的导热系数。

而有效应力的计算则相对复杂。在塑性区内，即R₀<r<R_p区间内，由于划分了n个厚度很小的环单元，可近似认为在圆环内耦合参数不发生改变。此时，围岩单元的应力平衡方程式(2)可改写为差分形式：

其中，表示第i+1层圆环计算单元的径向有效应力，表示第i层圆环计算单元的径向有效应力，rⁱ⁺¹表示第i+1层圆环计算单元距隧道中心的距离，rⁱ表示第i层圆环计算单元距隧道中心的距离，表示第i层圆环计算单元切向有效应力，E表示弹性模量，α表示热膨胀系数，β表示比奥系数，ΔTⁱ⁺¹表示第i+1层圆环计算单元的温度改变量，ΔTⁱ表示第i层圆环计算单元的温度改变量。

联立屈服方程式(1)，得到：

其中，表示第i+1层圆环计算单元的径向有效应力，表示第i+1层圆环计算单元的切向有效应力，表示第i层圆环计算单元的径向有效应力，rⁱ⁺¹表示第i+1层圆环计算单元距离隧道中心的距离，E表示围岩的弹性模量，β表示围岩的比奥系数，α表示围岩热膨胀系数，rⁱ表示第i层圆环计算单元距离隧道中心的距离，Tⁱ⁺¹表示第i+1层圆环计算单元的温度，Tⁱ表示第i层圆环计算单元的温度，βⁱ⁺¹表示第i+1层圆环计算单元的比奥系数，表示第i+1层圆环计算单元孔隙水压，表示第i层圆环计算单元的径向有效应力，表示i层圆环计算单元的切向有效应力，表示第i层圆环计算单元的孔隙水压，Tⁱ表示第i层圆环计算单元的温度，表示第i层圆环计算单元的摩擦角，c表示围岩瞬时粘聚力。

同理，相容方程式(7)的差分形式为：

其中，表示第i+1层圆环计算单元的切向塑性应变，表示第i+1层圆环计算单元的径向塑性应变，表示第i层圆环计算单元的径向塑性应变，表示第i+1层圆环计算单元的径向塑性应变，表示第i层圆环计算单元的径向塑性应变，表示第i+1层圆环计算单元的径向塑性应变，表示第i层圆环计算单元的切向塑性应变，表示第i层圆环计算单元的切向弹性应变，表示第i层圆环计算单元的径向弹性应变，K_ψ表示塑性乘子。

于是，当第i层圆环计算单元应力已知时，进一步可以得到i+1层圆环计算单元的塑性应变计算方程：

其中，表示第i+1层圆环计算单元的切向塑性应变，表示第i层圆环计算单元的切向塑性应变，表示第i层圆环计算单元的径向塑性应变，表示第i层圆环计算单元的径向弹性应变，K_ψ表示塑性乘子，rⁱ⁺¹表示第i+1层圆环计算单元距离隧道中心的距离，rⁱ表示第i层圆环计算单元距离隧道中心的距离，表示第i+1层圆环计算单元的切向弹性应变，表示第i层圆环计算单元的切向弹性应变，表示第i+1层圆环计算单元的径向塑性应变，表示第i层圆环计算单元的径向塑性应变，ψ表示剪胀角。

以上得到了塑性区内通过i层圆环计算单元温度、水压、应力计算i+1层环温度、水压、应力的计算公式。然而，由于弹塑性界面处的温度、水压、应力无法直接获取，所以需要借助某些算法获得模型最终解。

为了获得模型的解，将弹塑性界面处的温度T_p、孔隙水压P_w,p以及塑性圈半径R_p看作模型中的未知参数，所需要做的就是利用模型中边界条件解出此三个参数。联立式(11)和式(12)，可知在弹塑性区内径向应力与切向应力的和是关于此处到硐室中心距离的函数，即：

其中，C₆表示应力分布函数系数项，利用式(13)计算。

于是可以获得在r＝R_p处，

联立屈服方程式(1)，可以获得此时弹塑性界面处的有效应力此时获得了弹塑性界面上的温度、应力、水压。将式17和式18中边界条件改换为模型弹塑性区域边界条件，可以计算获得模型中的渗流量Q和热通量Φ。

可以看出，通过上述计算过程，一组初始条件对应一组开挖边界状态(T，σ_r,P_w)_r＝R0，需要确定正确的初始条件使计算得到的开挖边界径向应力、温度和孔隙压力与实际边界条件相匹配。采用梯度下降法可以得到正确的初始条件。梯度下降法是一种迭代优化算法，用于寻找函数的最小值。它利用损失函数负梯度的方向来确定每次迭代的搜索方向，保证每次迭代使目标函数的值逐渐减小。为此，损失函数可定义为：

LF(R_p,P_w,p,T_p)＝||(P_s,P_w,s,T_s)_calu-(P_s,P_w,s,T_s)_actu|| (23)

因此，在每一次迭代计算后，新的弹塑交界面处的计算参量为：

其中，(R_p,P_w,p,T)_j表示上一轮弹塑性界面状态假定解，(R_p,P_w,p,T)_j+1表示新的弹塑性界面状态假定解，η表示学习速率，表示损失函数的梯度，LF(R_p,P_w,p,T_p)表示损失函数，(P_s,P_w,s,T_s)_calu表示利用假定的(R_p,P_w,p,T)计算获得的开挖边界处的状态，(R_p,P_w,p,T)表示弹塑性界面状态假定解，(P_s,P_w,s,T_s)_actu表示根据实际工况设定的开挖边界状态参数，R_p表示塑性区半径，R_e表示弹性区半径，T表示温度，T_p表示弹塑性边界的温度，P_s表示开挖边界处的径向有效应力，即径向支护应力，P_w,s表示开挖边界处的孔隙水压，T_s表示开挖边界处的温度。

最终，获得了考虑热水力全耦合效应的弹塑性求解模型的解。

综上，本发明为热水力耦合地下工程建设提供了理论分析手段，可以通过模型计算获得塑性区半径、洞壁收敛位移等工程关键参数，为地下工程建设的支护设计提供参数支撑。尤为重要的是，研究表明热水力耦合工程中若不考虑全耦合效应则对于工程而言是危险的，原因是若忽略全耦合效应则围岩的卸荷效应将减弱，造成支护关键参数的减小。如塑性区半径往往作为锚杆支护中锚杆长度的关键设计参数，而忽略全耦合效应将造成塑性区半径计算结果偏小，导致洞室开挖过程中打设的锚杆偏短，从而造成支护失效，延误施工工期，甚至是造成人员伤亡。而本发明提出的计算模型考虑了全耦合效应，因此计算结果对工程而言是偏安全的。

Claims (7)

1.一种考虑热水力全耦合效应的深部圆形洞室开挖响应预测方法，其特征在于，包括以下步骤：

S1、在深部圆形洞室的预设范围内假定弹塑性解界面处的状态；

S2、计算获得弹性计算区域的水压分布函数系数项C1和温度分布函数系数项C3；

S3、根据假定的弹塑性解界面处的状态，计算得到弹塑性界面处的径向有效应力与切向有效应力的和；

S4、根据径向有效应力与切向有效应力的和以及屈服准则，分别计算得到弹塑性界面处的径向有效应力和切向有效应力；

S5、将塑性区域离散化，获得塑性区若干个等间距的圆环计算单元，共n+1个圆环计算节点，将S4的计算结果作为第1个圆环计算节点状态值；

S6、根据第i层圆环计算单元的状态值，计算第i+1层圆环计算单元的状态值；

S7、判断当前计算圆环单元距离洞室中心的距离是否大于开挖半径，若是，则将第i+1层圆环计算单元视为新的第i层圆环计算单元，并返回步骤S6，否则，进入步骤S8；

S8、判断开挖边界处通过差分计算得到的孔隙水压、温度以及径向有效应力与设定条件之间的相对误差是否均满足要求，若是，则完成对深部圆形洞室开挖的响应预测，否则，返回步骤S1，并使用梯度下降法计算获得新的弹塑性界面状态假定值。

2.根据权利要求1所述的考虑热水力全耦合效应的深部圆形洞室开挖响应预测方法，其特征在于，所述径向有效应力以及切向有效应力的和的表达式如下：

其中，G(R_p)表示弹塑性界面处的径向有效应力以及切向有效应力的和，σ'₀表示初始地应力，R_e表示边界距离洞室中心的距离，R_p表示弹塑性界面距离洞室中心的距离，E表示围岩的弹性模量，α表示围岩单元的线弹性密度，C₁和C₃分别表示弹性区域围岩中温度分布系数项和孔隙水压分布系数项，v表示泊松比，T_p和T₀分别表示弹塑性界面处的地温和初始地温，P_w,p和P_w,0分别表示弹塑性界面处的孔隙水压和初始孔隙水压。

3.根据权利要求1所述的考虑热水力全耦合效应的深部圆形洞室开挖响应预测方法，其特征在于，所述屈服准则的表达式如下：

其中，F表示屈服面函数，F(σ_θ,σ_r)＝0表示围岩单元达到塑性，σ_θ和σ_r分别表示切向有效应力和径向有效应力，表示瞬时摩擦角，c表示围岩瞬时粘聚力，c_ini表示初始粘聚力，γ_p表示等效塑性剪切应变，n'表示粘聚力演变曲线常数，表示切向塑性应变，表示径向塑性应变，D表示损伤变量，c_res表示残余粘聚力，ω表示瞬时孔隙率，ω_ini表示初始孔隙率，ω_res表示残余孔隙率，k表示瞬时渗透率，k_ini表示初始渗透率，λ表示围岩的瞬时导热系数，λ_s表示岩石骨架的导热系数，λ_w表示水的导热系数。

4.根据权利要求1所述的考虑热水力全耦合效应的深部圆形洞室开挖响应预测方法，其特征在于，所述第i层圆环计算单元的径向有效应力和切向有效应力的表达式分别如下：

其中，表示第i+1层圆环计算单元的径向有效应力，表示第i+1层圆环计算单元的切向有效应力，表示第i层圆环计算单元的径向有效应力，rⁱ⁺¹表示第i+1层圆环计算单元距离隧道中心的距离，E表示围岩的弹性模量，β表示围岩的比奥系数，α表示围岩热膨胀系数，rⁱ表示第i层圆环计算单元距离隧道中心的距离，Tⁱ⁺¹表示第i+1层圆环计算单元的温度，Tⁱ表示第i层圆环计算单元的温度，βⁱ⁺¹表示第i+1层圆环计算单元的比奥系数，表示第i+1层圆环计算单元孔隙水压，表示第i层圆环计算单元的径向有效应力，表示i层圆环计算单元的切向有效应力，表示第i层圆环计算单元的孔隙水压，Tⁱ表示第i层圆环计算单元的温度，表示第i层圆环计算单元的摩擦角，c表示围岩瞬时粘聚力；

第i+1层圆环计算单元的切向塑性应变和径向塑性应变的表达式分别如下：

其中，表示第i+1层圆环计算单元的切向塑性应变，表示第i层圆环计算单元的切向塑性应变，表示第i层圆环计算单元的径向塑性应变，表示第i层圆环计算单元的径向弹性应变，K_ψ表示塑性乘子，rⁱ⁺¹表示第i+1层圆环计算单元距离隧道中心的距离，rⁱ表示第i层圆环计算单元距离隧道中心的距离，表示第i+1层圆环计算单元的切向弹性应变，表示第i层圆环计算单元的切向弹性应变，表示第i+1层圆环计算单元的径向塑性应变，表示第i层圆环计算单元的径向塑性应变，ψ表示剪胀角。

5.根据权利要求1所述的考虑热水力全耦合效应的深部圆形洞室开挖响应预测方法，其特征在于，所述根据第i层圆环计算单元的状态值，计算第i+1层圆环计算单元的状态值，其具体为：

基于第i层圆环计算单元孔隙水压和温度，计算第i+1层圆环计算单元的孔隙水压和温度；

结合第i层圆环计算单元应力状态、计算参数以及第i+1层圆环计算单元的孔隙水压和温度，得到第i+1层圆环计算单元的有效应力和塑性应变；

利用第i+1层圆环计算单元的塑性应变，更新第i+1层圆环计算单元的计算参数；

将原第i+1层圆环计算单元看作新的第i层圆环计算单元，继续进行新的第i+1层圆环计算单元的状态值，直至圆环计算单元距隧道中心的距离不大于隧道开挖半径，完成对第i+1层圆环计算单元的状态值的计算，其中，状态值包括有效应力、塑性应变、孔隙水压以及温度。

6.根据权利要求5所述的考虑热水力全耦合效应的深部圆形洞室开挖响应预测方法，其特征在于，所述第i+1层圆环计算单元的孔隙水压的表达式如下：

其中，表示第i+1层圆环计算单元的孔隙水压，表示第i层圆环计算单元的孔隙水压，Q表示渗流通量，rⁱ表示第i层圆环计算单元距离隧道中心的距离，rⁱ⁺¹表示第i+1层圆环计算单元距离隧道中心的距离，kⁱ表示第i层圆环计算单元的渗透率；

所述第i+1层圆环计算单元的温度的表达式如下：

其中，Tⁱ⁺¹表示第i+1层圆环计算单元的温度，Tⁱ表示第i层圆环计算单元的温度，Φ表示热通量，λⁱ表示第i层圆环计算单元的导热系数。

7.根据权利要求1所述的考虑热水力全耦合效应的深部圆形洞室开挖响应预测方法，其特征在于，所述新的弹塑性界面状态假定值的表达式如下：

LF(R_p,P_w,p,T_p)＝||(P_s,P_w,s,T_s)_calu-(P_s,P_w,s,T_s)_actu||

其中，(R_p,P_w,p,T)_j表示上一轮弹塑性界面状态假定解，(R_p,P_w,p,T)_j+1表示新的弹塑性界面状态假定解，η表示学习速率，表示损失函数的梯度，LF(R_p,P_w,p,T_p)表示损失函数，(P_s,P_w,s,T_s)_calu表示利用假定的(R_p,P_w,p,T)计算获得的开挖边界处的状态，(R_p,P_w,p,T)表示弹塑性界面状态假定解，(P_s,P_w,s,T_s)_actu表示根据实际工况设定的开挖边界状态参数，R_p表示塑性区半径，R_e表示弹性区半径，T表示温度，T_p表示弹塑性边界的温度，P_s表示开挖边界处的径向支护力，P_w,s表示开挖边界处的孔隙水压，T_s表示开挖边界处的温度。