一种高温岩石细观节理本构模型构建方法

1.本发明涉及岩石力学数值模拟技术领域，尤其涉及一种高温岩石细观节理本构模型构建方法。


背景技术：

2.随着地球浅部资源的日益枯竭，对深部资源开采和深部地下空间开发利用的需求不断增加。而深部岩体在高温高压条件下的有着异于浅部岩体的物理性质和力学行为，这给增强型地热系统、高放射性废弃物处置等工程应用带来了许多挑战。而岩石作为工程应用中最为常见的工程材料之一，理解它在热力耦合条件下的力学本构行为是我们继续向地球深部进军过程中的一个关键问题。
3.然而到目前为止，对高温岩石的热力耦合行为的研究仍不充分，主要原因是在实际高温下对岩石进行力学试验的实验技术手段十分缺乏。因此，将数值模拟作为替代技术以复现岩石的热力耦合行为的实时过程成了一种可行之道。作为数值模拟手段之一，离散元法等不连续方法以离散方式完成微裂纹的成核和扩展。而这种微裂纹同时受到滑移效应和温度效应的影响，需要一种同时考虑裂纹滑移效应和温度效应的岩石细观节理本构模型来描述其热力学行为，以在离散元软件中更精确地复现整个岩石模型的热力耦合行为。


技术实现要素：

4.本发明要解决的技术问题是针对上述现有技术的不足，提供一种高温岩石细观节理本构模型构建方法，实现高温岩石细观节理本构模型的构建。
5.为解决上述技术问题，本发明所采取的技术方案是：一种高温岩石细观节理本构模型构建方法，包括以下步骤：
6.步骤1、建立考虑裂纹滑移效应与温度效应的岩石细观节理本构关系的初步数学表达式；
7.步骤1.1、考虑热力耦合作用下热应力裂纹的滑移效应，将库仑摩擦滑动准则的残余内聚力和残余内摩擦角修改为滑移位移依赖性；
8.热应力裂纹的剪切参数受滑移变形的影响表达式为：
[0009][0010][0011]
其中，和是受裂纹滑移距离影响的残余内聚力和残余内摩擦角，和是对应的残余内聚力和残余内摩擦角的滑移弱化或强化系数方程，c0和是岩石细观节理未受温度影响的初始内聚力和内摩擦角；
[0012]
步骤1.2、考虑热力耦合作用下的温度效应，将库仑摩擦滑动准则的内聚力、内摩擦角、抗拉强度修改为温度依赖性；
[0013]
考虑温度效应，将库仑摩擦滑动准则节理的内聚力c0、内摩擦角抗拉强度σ
t0
修
改为温度依赖性，表达式为:
[0014]ctemp
＝f
temp-c
c0ꢀꢀꢀ
(3)
[0015][0016]
σ
temp-t
＝f
temp-t
σ
t0
ꢀꢀꢀ
(5)
[0017]
其中，c
temp
、σ
temp-t
分别为受温度效应影响的内聚力、内摩擦角、抗拉强度；f
temp-c
、f
temp-t
分别为内聚力、内摩擦角、抗拉强度的温度依赖性系数方程；
[0018]
将高温作用下，库仑摩擦滑动准则的岩石细观节理最大抗拉强度t
max
和最大抗剪切强度s
max
修改为温度依赖性，如下公式所示：
[0019]
t
temp-max
＝-σ
temp-t ac＝f
temp-t
σ
t0 acꢀꢀꢀ
(6)
[0020][0021]
其中，t
temp-max
为节理的温度依赖性最大抗拉强度，s
temp-max
为节理的温度依赖性最大抗剪切强度，ac为节理面积，fn为节理法向力；
[0022]
步骤1.3、建立同时考虑裂纹滑移效应与温度效应的岩石细观节理本构关系的初步数学表达式；
[0023]
对于张拉裂纹，设定残余抗拉强度为0；对于剪切裂纹，其力学行为同时受到温度和裂纹滑移变形的影响，因此，剪切裂纹相应的残余剪切参数的力学关系式变为如下公式所示：
[0024][0025][0026][0027]
其中，和是同时考虑节理滑移效应和温度效应的残余内聚力和残余内摩擦角,是同时考虑裂纹滑移效应和温度效应的最大残余剪切强度；
[0028]
步骤2、制备多个岩石试样，分别对岩石试样进行测量分析和高温岩石力学试验，获取岩石的宏观物理力学性质；
[0029]
步骤2.1、对目标岩石试样进行矿物成分测量；
[0030]
对目标岩石试样进行矿物成分测量，确定岩石试样的矿物晶体组成成分和比例；
[0031]
步骤2.2、对岩石试样的所有矿物晶体的热力学性质进行收集整理；
[0032]
对岩石试样的所有矿物晶体的热力学性质进行收集整理，其中，热学参数包括线性热膨胀系数α
t
，热导系数k，比热容c
p
；力学参数包括：密度ρm，泊松比νm，杨氏模量em，内聚力cm，内摩擦角抗拉强度σ
tm
；矿物晶体间接触面的力学性质取相邻矿物晶体力学参数的平均值；
[0033]
步骤2.3、对岩样试样进行多种高温岩石力学试验，得到岩石试样的应力应变曲线和抗拉强度；
[0034]
根据未经热处理岩石试样的应力应变关系和抗拉强度，确定岩石试样数值几何模
型的初始力学参数，包括密度ρ，泊松比ν，杨氏模量e，内聚力c，内摩擦角和抗拉强度σ
t
；
[0035]
根据应力应变曲线反算出基于实验室试验结果的内聚力温度依赖性系数方程f
temp-c
、内摩擦角温度依赖性系数方程抗拉强度温度依赖性系数方程f
temp-t
，残余内聚力的裂纹滑移弱化或强化系数方程残余内摩擦角的裂纹滑移弱化/强化系数方程
[0036]
对岩样试样进行多种高温岩石力学试验时需要对岩石试样先进行加热-冷却预处理；
[0037]
步骤3、在离散元数值模拟软件中建立岩石试样的数值几何模型，将目标岩石的热力学参数和初步本构关系输入到数值几何模型的单元体和单元体接触面；
[0038]
步骤3.1、根据要建模的实验室岩石试样的尺寸，在离散元数值模拟软件中生成1:1的数值几何模型；
[0039]
步骤3.2、根据不同矿物成分在多晶体岩石试样中所占比例，对数值几何模型中的数值单元体进行遍历，将单元按照步骤2.1测得的真实矿物比例随机分组，并定义成相应的矿物晶体；
[0040]
步骤3.3、对数值几何模型中的矿物晶体和晶内、晶体间接触面进行识别，将热力学参数和本构关系赋予晶体单元和晶体间接触面；
[0041]
其中，晶体单元设置为弹性本构模型，晶体间接触面的本构关系设置为步骤1中获得的考虑裂纹滑移效应与温度效应的岩石细观节理本构关系的初步数学表达式；
[0042]
步骤4、复制实验室试验条件，在离散元数值模拟软件中进行平行数值模拟试验；
[0043]
复制实验室试验条件，对离散元数值模拟软件中所建模型进行平行加热加载数值模拟试验；
[0044]
步骤5、对比数值模拟结果和实验室试验结果，利用反分析法，获得对应岩石试样的细观节理本构关系所需系数方程，完成高温岩石细观节理本构模型的构建；
[0045]
步骤5.1、根据步骤4进行的加热加载数值模拟试验，对比数值几何模型和试件在同样加热-冷却预处理后受载荷作用下的应力应变关系和抗拉强度，基于模拟结果和试验结果的差别，修改数值几何模型的初始力学参数，再次进行加热加载平行数值模拟试验，并再次对比试验结果，由此反复，直至数值模拟结果与试验结果相一致或在设定误差范围内，由此确定不同温度下的力学参数与初始力学参数的定量关系；
[0046]
步骤5.2、基于获得的不同温度下的定量关系，建立不同温度下岩石试样的内聚力温度依赖性系数方程f
temp-c
、内摩擦角温度依赖性系数方程抗拉强度温度依赖性系数方程f
temp-t
，残余内聚力的裂纹滑移弱化/强化系数方程残余内摩擦角的裂纹滑移弱化/强化系数方程
[0047]
步骤5.3、将所得到的内聚力温度依赖性系数方程f
temp-c
、内摩擦角温度依赖性系数方程抗拉强度温度依赖性系数方程f
temp-t
、残余内聚力的裂纹滑移弱化/强化系数方程数方程残余内摩擦角的裂纹滑移弱化/强化系数方程代入剪切裂纹相应的残余剪切参数的力学关系式中，完成高温岩石细观节理本构模型的构建。
[0048]
采用上述技术方案所产生的有益效果在于：本发明提供的一种高温岩石细观节理本构模型构建方法，结合高温岩石力学试验，使用离散元数值模拟软件，利用反分析法构建高温岩石细观节理本构模型，能够同时考虑裂纹滑移效应与温度效应对岩石细观节理的影响，可在离散元软件中进一步还原岩石受到高温作用后细观节理的成核和扩展机制，经该方法构建的本构模型通过数值模拟试验验证，所得到的结果与实验室试验结果吻合度较高，可应用于岩石力学数值模拟领域。
附图说明
[0049]
图1为本发明实施例提供的一种高温岩石细观节理本构模型构建方法的流程图；
[0050]
图2为本发明实施例提供的岩石试样在不同温度处理后进行单轴压缩试验得到的应力应变曲线图；
[0051]
图3为本发明实施例提供的通过实验室试验得到的内聚力、内摩擦角和抗拉强度的温度依赖性关系图，其中，(a)为内聚力和内摩擦角与温度的依赖关系，(b)为抗拉强度与温度的依赖关系；
[0052]
图4为本发明实施例提供的通过实验室试验得到的残余内聚力和残余内摩擦角的滑移弱化或强化系数方程；
[0053]
图5为本发明实施例提供的通过udec建立的岩石试样的数值几何模型，其中，(a)为圆柱试样的数值几何模型，(b)为圆盘试样的数值几何模型；
[0054]
图6为本发明实施例提供的对数值几何模型进行加热-冷却处理的示意图；
[0055]
图7为本发明的一个实施例中的对数值几何模型进行加载模拟试验的示意图，其中，(a)为进行单轴压缩试验的圆柱模型，(b)为进行巴西劈裂的圆盘模型；
[0056]
图8为本发明的一个实施例中的对数值几何模型进行单轴压缩试验得到的应力应变曲线图。
具体实施方式
[0057]
下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。
[0058]
本实施例中，一种高温岩石细观节理本构模型构建方法，如图1所示，包括以下步骤：
[0059]
步骤1、建立考虑裂纹滑移效应与温度效应的岩石细观节理本构关系的初步数学表达式；
[0060]
步骤1.1、考虑热力耦合作用下热应力裂纹的滑移效应，将库仑摩擦滑动准则(coulomb friction criteria)的残余内聚力和残余内摩擦角修改为滑移位移依赖性；
[0061]
热应力裂纹的剪切参数受滑移变形的影响表达式为：
[0062][0063][0064]
其中，和是受裂纹滑移距离影响的残余内聚力和残余内摩擦角，和是对应的残余内聚力和残余内摩擦角的滑移弱化或强化系数方程，c0和是岩石细
观节理未受温度影响的初始内聚力和内摩擦角；
[0065]
步骤1.2、考虑热力耦合作用下的温度效应，将库仑摩擦滑动准则的内聚力、内摩擦角、抗拉强度修改为温度依赖性；
[0066]
考虑温度效应，将库仑摩擦滑动准则节理的内聚力c0、内摩擦角抗拉强度σ
t0
修改为温度依赖性，表达式为:
[0067]ctemp
＝f
temp-c
c0ꢀꢀꢀ
(3)
[0068][0069]
σ
temp-t
＝f
temp-t
σ
t0
ꢀꢀꢀ
(5)
[0070]
其中，c
temp
、σ
temp-t
分别为受温度效应影响的内聚力、内摩擦角、抗拉强度；f
temp-c
、f
temp-t
分别为内聚力、内摩擦角、抗拉强度的温度依赖性系数方程；
[0071]
将高温作用下，库仑摩擦滑动准则的岩石细观节理最大抗拉强度t
max
和最大抗剪切强度s
max
修改为温度依赖性，如下公式所示：
[0072]
t
temp-max
＝-σ
temp-t ac＝f
temp-t
σ
t0 acꢀꢀꢀ
(6)
[0073][0074]
其中，t
temp-max
为节理的温度依赖性最大抗拉强度，s
temp-max
为节理的温度依赖性最大抗剪切强度，ac为节理面积，fn为节理法向力；
[0075]
步骤1.3、建立同时考虑裂纹滑移效应与温度效应的岩石细观节理本构关系的初步数学表达式；
[0076]
对于张拉裂纹，可认为裂纹在垂直于裂纹面方向上没有抵抗力，因此设定残余抗拉强度为0；对于剪切裂纹，其力学行为同时受到温度和裂纹滑移变形的影响，因此，剪切裂纹相应的残余剪切参数的力学关系式变为如下公式所示：
[0077][0078][0079][0080]
其中，和是同时考虑节理滑移效应和温度效应的残余内聚力和残余内摩擦角,是同时考虑裂纹滑移效应和温度效应的最大残余剪切强度；
[0081]
剪切裂纹相应的残余剪切参数的力学关系式即为岩石细观节理本构关系的初步数学表达，该表达式中的系数方程均是待定的。
[0082]
步骤2、制备多个岩石试样，分别对岩石试样进行测量分析和多种高温岩石力学试验，获取岩石的宏观物理力学性质；
[0083]
步骤2.1、对目标岩石试样行矿物成分测量；
[0084]
对目标岩石试样进行矿物成分测量，确定岩石试样的矿物晶体组成成分和比例；
[0085]
步骤2.2、对岩石试样的所有矿物晶体的热力学性质进行收集整理；
[0086]
对岩石试样的所有矿物晶体的热力学性质进行收集整理，其中，热学参数包括线性热膨胀系数α
t
，热导系数k，比热容c
p
；力学参数包括：密度ρm，泊松比νm，杨氏模量em，内聚
力cm，内摩擦角抗拉强度σ
tm
；矿物晶体间接触面的力学性质取相邻矿物晶体力学参数的平均值；
[0087]
步骤2.3、对岩样试样进行多种高温岩石力学试验，得到岩石试样的应力应变曲线和抗拉强度；
[0088]
同时，根据未经热处理岩石试样的应力应变关系和抗拉强度，确定岩石试样数值几何模型的初始力学参数，包括密度ρ，泊松比ν，杨氏模量e，内聚力c，内摩擦角和抗拉强度σ
t
；
[0089]
根据应力应变曲线反算出基于实验室试验结果的内聚力温度依赖性系数方程f
temp-c
、内摩擦角温度依赖性系数方程抗拉强度温度依赖性系数方程f
temp-t
，残余内聚力的裂纹滑移弱化/强化系数方程残余内摩擦角的裂纹滑移弱化/强化系数方程
[0090]
对岩样试样进行多种高温岩石力学试验前需要对岩石试样进行加热-冷却预处理，处理方法如下：
[0091]
利用电热炉，对岩石试样先进行缓慢加热，加热速度不超过10℃/min，达到目标温度后，继续保温1-3小时，然后在加热炉内进行缓慢降温至室温；
[0092]
本实施例中，所用岩石试样均按国际岩石力学学会建议试验标准制备，分别是3件直径50mm，长110mm的单轴压缩试验用圆柱试件和3件直径50mm，厚25mm的巴西劈裂试验用圆盘，两种试件均为艾本斯托克花岗岩。
[0093]
本实施例首先对目标岩石试样进行矿物成分测量，测得其中长石成分占比50％，石英石成分占比44％，云母成分占比6％。
[0094]
然后对岩石试样的所有矿物晶体的热力学性质参数进行收集整理，其中力学参数包括：密度ρm，泊松比νm，杨氏模量em，内聚力cm，内摩擦角抗拉强度σ
tm
，如表1所示。
[0095]
表1岩石试样的所有矿物晶体力学参数值
[0096][0097][0098]
矿物晶体间接触面的力学参数取相邻晶体力学性质的平均值，如表2所示。
[0099]
表2矿物晶体间接触面的力学参数值
[0100][0101]
热学参数包括：线性热膨胀系数α
t
，热导系数k和比热容c
p
，如表3所示。
[0102][0103]
收集好各项参数后，需要对岩石试样进行加热-冷却预处理，本实施例以单轴压缩试验用圆柱试件为例：
[0104]
首先利用电热炉，对岩石试样先进行缓慢加热，加热速度不超过10℃/min，然后在加热炉内进行缓慢降温至室温，将3个试样分别加热至25℃、400℃、600℃，其中加热至400℃和600℃的试样需在目标温度保温36小时。巴西劈裂试验用圆盘试样的加热处理方式相同。
[0105]
待试样冷却至室温后，对圆柱体试样进行单轴压缩试验，对圆盘试样进行巴西劈裂试验，加载速率分别是0.1mm/min和0.05mm/min，记录不同温度处理后的试件的巴西劈裂强度和单轴压缩应力应变曲线，如图2所示。
[0106]
最后反算出基于实验室试验结果的内聚力温度依赖性系数方程f
temp-c
、内摩擦角温度依赖性系数方程抗拉强度温度依赖性系数方程f
temp-t
，如图3所示；残余内聚力的裂纹滑移弱化或强化系数方程残余内摩擦角的裂纹滑移弱化/强化系数方程如图4所示。
[0107]
步骤3、在离散元数值模拟软件中建立岩石试样试样的数值几何模型，将目标岩石的热力学参数和初步本构关系输入到数值几何模型的单元体和单元体接触面；
[0108]
步骤3.1、根据要建模的实验室岩石试样的尺寸，在离散元数值模拟软件中生成1:1的岩石试样数值几何模型；
[0109]
步骤3.2、根据不同矿物成分在多晶体岩石试样中所占比例，对数值几何模型中的数值单元体进行遍历，将单元按照步骤2.1测得的真实矿物比例随机分组，并定义成相应的矿物晶体；
[0110]
步骤3.3、对数值几何模型中的矿物晶体和晶内、晶体间接触面进行识别，将热力学参数和本构关系赋予晶体单元和晶体间接触面；
[0111]
其中，晶体单元设置为弹性本构模型，晶体间接触面的本构关系设置为步骤1中获得的考虑裂纹滑移效应与温度效应的岩石细观节理本构关系的初步数学表达式，初步数学表达式中系数方程如图3、图4所示；
[0112]
本实施例中，在离散元数值模拟软件udec中建立1：1的岩石试样数值几何模型，即直径0.05，长0.11的圆柱体模型，和直径0.05，厚0.025的圆盘模型；
[0113]
遍历几何模型中的数值单元体，将单元体按照石英44％、长石50％、云母6％的比例随机分组，并将单元体定义为相应的矿物晶体，如图5所示。
[0114]
步骤4、复制实验室试验条件，在离散元数值模拟软件中进行平行数值模拟试验；
[0115]
包括加热-冷却预处理模拟，对预处理后的模型进行单轴/三轴加载与巴西劈裂试验模拟；
[0116]
本实施例，复制实验室试验条件，对离散元数值模拟软件中所建模型进行平行加热加载数值模拟试验。对数值几何模型进行加热-冷却预处理，在模型外表面施加初始温度25℃，并以10℃/min速率增加的温度边界条件，在边界温度条件达到400℃、600℃时分别保存模型；重建保存的模型，并将温度保持在400℃、600℃，直至模型内温度分布统一；然后以-10℃/min速率降低边界温度条件直至25℃，保持25℃的温度边界条件至模型温度完全降低至25℃，整个过程如图6所示。
[0117]
对数值几何模型进行单轴压缩和巴西劈裂模拟试验，如图7所示，在圆柱模型顶部和底部边界施加2.1
×
10-2m/s的加载速率，圆盘模型则被施加了2.2
×
10-3m/s的加载速率。
[0118]
步骤5、对比数值模拟结果和实验室试验结果，利用反分析法，获得对应岩石试样的细观节理本构关系所需系数方程，完成高温岩石细观节理本构模型的构建；
[0119]
步骤5.1、根据步骤4进行的加热加载数值模拟试验，对比数值几何模型和试件在同样加热-冷却预处理后受载荷作用下的应力应变关系和抗拉强度，基于模拟结果和试验结果的差别，修改数值几何模型的初始力学参数，再次进行加热加载平行数值模拟试验，并再次对比试验结果，由此反复，直至数值模拟结果与试验结果相一致或在设定误差范围内(小于10％)，由此确定不同温度下的力学参数(温度下内聚力，内摩擦角，抗拉强度，残余内聚力，残余内摩擦角等)与初始力学参数的定量关系；本实施例，对数值几何模型进行单轴压缩试验得到的应力应变曲线如图8所示，将图8与图2对比可知，数值模拟结果与试验结果吻合良好。
[0120]
步骤5.2、基于获得的不同温度下的定量关系，建立不同温度下岩石试样的内聚力
温度依赖性系数方程f
temp-c
、内摩擦角温度依赖性系数方程抗拉强度温度依赖性系数方程f
temp-t
，残余内聚力的裂纹滑移弱化/强化系数方程残余内摩擦角的裂纹滑移弱化/强化系数方程
[0121]
步骤5.3、将所得到的内聚力温度依赖性系数方程f
temp-c
、内摩擦角温度依赖性系数方程抗拉强度温度依赖性系数方程f
temp-t
、残余内聚力的裂纹滑移弱化/强化系数方程数方程残余内摩擦角的裂纹滑移弱化/强化系数方程代入剪切裂纹相应的残余剪切参数的力学关系式(8)、(9)、(10)中，完成高温岩石细观节理本构模型的构建。
[0122]
最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明权利要求所限定的范围。