一种考虑疲劳和时间效应的岩体稳定性计算方法与流程

1.本发明涉及岩体工程计算领域，尤其涉及一种考虑疲劳和时间效应的岩体稳定性计算方法。


背景技术：

2.近年来，我国的交通、水利、采矿等诸多领域的岩体工程建设方兴未艾，例如边坡、水电厂房、隧道、地基等。岩体是完整岩石和结构面组成的地质结构体，存在着大量的断层、节理以及软弱夹层等。由于地质体所处环境复杂性，与一般的介质相比，岩体具有非均质性和时间效应，且受到环境疲劳作用，岩体工程稳定性问题日益突出，因此，研究考虑流变和疲劳作用的岩体稳定性计算方法具有重要意义，对于保证复杂环境下岩体工程的建设运营安全具有重要的指导性作用。
3.目前对于岩体工程的稳定性计算一般采用理想弹塑性的宏观计算模型，传统岩体稳定性计算方法中，采用线性强度准则的本构模型，无法反映岩体强度的非线性特征；其将各个地层看成均质的材料赋计算参数，采用宏观计算无法反映岩体介质的非均质性；采用元件组合的流变模型，无法反映荷载、时效和疲劳综合作用下岩体劣化的本质，以及对应的稳定性变化。对于综合考虑流变和疲劳作用的岩体稳定性存在较大的误差。


技术实现要素：
:
4.本发明提供一种考虑疲劳和时间效应的岩体稳定性计算方法，以克服上述技术问题。
5.为了实现上述目的，本发明的技术方案是：
6.一种考虑疲劳和时间效应的岩体稳定性计算方法，包括如下步骤：
7.s1：将岩石划分为若干基元v，以获取所述基元的力学参数；
8.所述基元的力学参数包括反映岩体特征的经验参数mb、岩体破碎及完整程度参数s、岩体质量参数a和损伤前的弹性模量e0；
9.s2：根据所述基元的力学参数，获取所述基元的hoek-brown弹塑性损伤模型；
10.s3：获取基元的时效损伤变量d(t)；
11.s4：获取基元的疲劳损伤变量dr；
12.s5：根据荷载损伤变量dm、所述时效损伤变量d(t)和所述疲劳损伤变量dr，获取耦合损伤变量dc：
13.s6：根据所述基元的hoek-brown弹塑性损伤模型和所述耦合损伤变量，获取t
m+1
时间步下的名义应力张量σ'
m+1
；
14.s7：根据第t
m+1
时间步下的名义应力张量σ'
m+1
，获取第t
m+1
时间步下的耦合损伤变量值。
15.有益效果：本发明的一种考虑疲劳和时间效应的岩体稳定性计算方法，针对岩体强度的非线性和细观的非均匀性特点，将岩石划分为若干基元，获取能够反映岩体特征的
经验参数、岩体破碎及完整程度参数、岩体质量参数和损伤前的弹性模量；并根据这些参数建立hoek-brown弹塑性损伤模型，然后基于岩体的时间效应特点，获取基元的时效损伤变量，并结合疲劳损伤变量和载荷损伤变量得到基元的耦合损伤变量，最终得到随时间变化的耦合损伤变量的值，实现流变效应和疲劳效应作用下非均质性岩体随着时间的稳定性计算。
附图说明
16.为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图做一简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
17.图1为本发明的考虑疲劳和时间效应的岩体稳定性计算方法流程图；
18.图2a为本发明的实施例中x0为10时的rve力学性质的概率密度分布图；
19.图2b为本发明的实施例中ζ为0.5时的rve力学性质的概率密度分布图；
20.图3为本发明的实施例中的poytin-thomson模型示意图；
21.图4为本发明的实施例中的岩体损伤变量随疲劳次数变化规律图；
22.图5为本发明的实施例中的细观分布下岩石试样计算模型图；
23.图6a为本发明的实施例中ζ为0.5时的mb值分布云图；
24.图6b为本发明的实施例中ζ为2时的mb值分布云图；
25.图6c为本发明的实施例中ζ为4时的mb值分布云图；
26.图6d为本发明的实施例中ζ为6时的mb值分布云图；
27.图6e为本发明的实施例中ζ为8时的mb值分布云图；
28.图6f为本发明的实施例中ζ为10时的mb值分布云图；
29.图7a为本发明的实施例中ζ为0.5时的mb值区域rve统计示意图；
30.图7b为本发明的实施例中ζ为2时的mb值区域rve统计示意图；
31.图7c为本发明的实施例中ζ为4时的mb值区域rve统计示意图；
32.图7d为本发明的实施例中ζ为6时的mb值区域rve统计示意图；
33.图7e为本发明的实施例中ζ为8时的mb值区域rve统计示意图；
34.图7f为本发明的实施例中ζ为10时的mb值区域rve统计示意图；
35.图8a为本发明的实施例中ζ为0.5时的mb和损伤值演化规律图；
36.图8b为本发明的实施例中ζ为4时的mb和损伤值演化规律图；
37.图9为本发明的实施例中不同应力水平下岩石的蠕变曲线图；
38.图10a为本发明的实施例中ζ为6时的岩样损伤带随时间变化曲线；
39.图10b为本发明的实施例中ζ为6时的岩样损伤带随时间变化曲线；
40.图11为本发明的实施例中边坡计算有限元模型图；
41.图12a为本发明的实施例中ζ为2时边坡弹性模量及损伤区变化规律图；
42.图12b为本发明的实施例中ζ为6时边坡弹性模量及损伤区变化规律图；
43.图12c为本发明的实施例中ζ为10时边坡弹性模量及损伤区变化规律图；
44.图13a为本发明的实施例中时间月为0(月)时的边坡损伤带示意图；
45.图13b为本发明的实施例中时间月为4(月)时的边坡损伤带示意图；
46.图13c为本发明的实施例中时间月为8(月)时的边坡损伤带示意图；
47.图13d为本发明的实施例中时间月为12(月)时的边坡损伤带示意图；
48.图13e为本发明的实施例中时间月为16(月)时的边坡损伤带示意图；
49.图13f为本发明的实施例中时间月为20(月)时的边坡损伤带示意图；
50.图14为本发明的实施例中的坡面安全系数随时间变化规律图；
51.图15为本发明的实施例中不同冻融次数下边坡最大耦合损伤随时间变化规律图(t/月)；
52.图16a为本发明的实施例中冻融80次后时间为0(月)的边坡位移示意图；
53.图16b为本发明的实施例中冻融80次后时间为4(月)的边坡位移示意图；
54.图16c为本发明的实施例中冻融80次后时间为8(月)的边坡位移示意图；
55.图16d为本发明的实施例中冻融80次后时间为12(月)的边坡位移示意图；
56.图16e为本发明的实施例中冻融80次后时间为16(月)的边坡位移示意图；
57.图16f为本发明的实施例中冻融80次后时间为20(月)的边坡位移示意图；
58.图17为本发明的实施例中坡面监测点位移随时间变化规律图；
59.图18a为本发明的实施例中锚索预应力为200kn时的边坡损伤区示意图；
60.图18b为本发明的实施例中锚索预应力为300kn时的边坡损伤区示意图；
61.图18c为本发明的实施例中锚索预应力为400kn时的边坡损伤区示意图；
62.图18d为本发明的实施例中锚索预应力为500kn时的边坡损伤区示意图。
具体实施方式
63.为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
64.本发明实施例提供一种考虑疲劳和时间效应的岩体稳定性计算方法，如图1所示；包括如下步骤：
65.s1：为从连续介质损伤力学的角度来研究如空洞、晶体的影响导致的微观不连续性，需要对岩体材料微观力学性质进行均匀化处理。首先将岩石划分为若干基元v，其次利用统计方法对基元v的物理力学性质进行描述，以获取所述基元的力学参数；
66.所述基元的力学参数包括但不限于反映岩体特征的经验参数mb、岩体破碎及完整程度参数s、岩体质量参数a和损伤前的弹性模量e0等。
67.具体的为，使用具有门槛值的幂函数律来描述岩体材料强度极值分布律。假设基元中的力学参数pz可以用weibull统计分布函数进行描述，
68.所述基元的力学参数获取方法如下：
69.pz＝p0·
xzꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
70.式中，p0为岩石试样的宏观性质；xz为weibull随机数；pz为基元的力学参数；
71.weibull随机数xz获取方法如下：
72.首先根据weibull分布的累积分布函数和概率密度函数如下：
[0073][0074][0075]
式中：w(x)为weibull分布的累积分布函数(cumulative distribution function,cdf)；w(x)为weibull分布的概率密度函数(probability density function,pdf)；x0》0是比例参数；ζ为形状参数；x为随机变量；
[0076]
式(3)中，随着ζ的增大，岩石中
[0077]
xz＝x0(-ln(um))
1/ζ
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0078]
式中：xz为weibull随机数，um为均匀分布的随机数；
[0079]
s2：根据所述基元的力学参数，获取所述基元的hoek-brown弹塑性损伤模型(霍克-布朗弹塑性损伤模型)，以反映岩体强度非线性、塑性软化特性；
[0080]
具体的，通常在已有的岩石力学屈服函数基础上引入荷载损伤变量dm，用于表征损伤对强度的弱化作用。本实施例假定受荷过程中损伤对参数mb和s产生弱化作用，提出了考虑塑性损伤的所述基元的hoek-brown弹塑性损伤模型，其表达式如下：
[0081][0082][0083]
式中，f为屈服函数；j2为第二偏应力不变量；θ
σ
为罗德角；dm为荷载损伤变量；σ
ci
为完整岩石的单轴抗压强度；p为静水压力；m
bg
为mb对应的塑性势函数参数；sg为s对应的塑性势函数参数；ag为a对应的塑性势函数参数；g为塑性势函数；当m
bg
＝mb、sg＝s、ag＝a时，为关联流动法则，否则为非关联流动法则。
[0084]
其中，荷载损伤变量dm获取如下：由于工程开挖产生的应力重分布效应会使岩体内部产生塑性损伤，损伤随着等效塑性应变地累积而加剧，荷载损伤变量dm可以表示为的幂指数形式：
[0085][0086]
式中，α为决定损伤后岩石材料软化曲线的初始斜率，取值范围为(0,+∞)；β为决定岩石最大损伤值的参数，取值范围为[0,1)；为等效塑性应变；
[0087][0088]
式中，为空间坐标系(即主应力空间)x轴方向上的主塑性应变；为空间坐标系(即主应力空间)y轴方向上的主塑性应变；为空间坐标系(即主应力空间)z轴方向上的主塑性应变；
[0089]
s3：获取时效损伤变量d(t)，方法如下；
[0090]
当轴压小于长期强度值时，所述基元的力学模型可以用图3所示的poytin-thomson模型进行表达，建立基元的蠕变方程为：
[0091][0092]
式中，ε为应变；σ0为轴向的应力；k1、k2均为弹性元件的力学参数；η为黏壶的粘滞系数；
t
为流变的时间；在本实施例中，poytin-thomson模型为现有技术，其中包括两个弹性元件，因此k1、k2分别为模型中的两个弹性元件的力学参数。
[0093]
化简式(9)可得
[0094]
ε＝σ0[p+qexp(-rt)]
ꢀꢀꢀꢀꢀꢀꢀꢀ
(10)
[0095]
式中，p、q、r均为中间参数，其中，
[0096]
从公式(10)可以得到，弹性模量随时间t发生变化的表达式为：
[0097]
e(t)＝p+qexp(-rt)
ꢀꢀꢀꢀꢀꢀꢀꢀ
(11)
[0098]
e(t)为弹性模量随着时间t变化的函数；
[0099]
根据连续损伤力学理论，材料模量的损伤程度可以由弹性模量法进行测定，即：
[0100][0101]
式中，d(t)为时效损伤变量，即由时间效应引起的损伤变量；e0为损伤前的弹性模量；
[0102]
s4：获取疲劳损伤变量dr，如下：
[0103]
假设第n
t
次疲劳扰动作用下断面的损伤值为d
t
可以测得，则：
[0104][0105]
式中：d
t
为第n
t
次疲劳扰动作用下岩石断面的损伤值；dr为疲劳损伤变量dr，即第n
t
次疲劳扰动之前的n次疲劳损伤变量；n
t
为疲劳扰动次数的编号；n为第n
t
次疲劳扰动之前的疲劳扰动总次数；w和u均为材料参数；
[0106]
经过第n
t
次疲劳扰动作用下岩石断面的损伤值d
t
可以通过岩石孔隙率和纵波波速获得，计算公式如式(14)所示：
[0107][0108]
式中，n0为岩石疲劳作用前的孔隙率；n为岩石疲劳作用后的孔隙率；v
p
为岩石疲劳作用后的纵波波速；根据通过式(13)拟合疲劳作用次数与疲劳损伤之间的关系，结果如图4所示。
[0109]
步骤5：根据荷载损伤变量dm、所述时效损伤变量d(t)和所述疲劳损伤变量dr，获取耦合损伤变量dc：
[0110]
对于服役过程中的北方寒区工程，岩体在荷载的长期作用下产生时效损伤和疲劳
损伤，也会产生冻融疲劳损伤。对上述多种损伤进行耦合，来分析多种因素影响下的岩体稳定性问题。综合考虑疲劳损伤变量dr、荷载损伤变量dm和时效损伤变量d(t)，则此时获取耦合损伤变量如下为：
[0111]
dc＝d(t)+dm+d
r-d(t)d
m-d
rdm-d
rdt
+d(t)d
mdr
ꢀꢀꢀꢀꢀ
(15)
[0112]
式中，dc为耦合损伤变量；dm为荷载损伤变量；
[0113]
s6：根据所述基元的hoek-brown弹塑性损伤模型和所述耦合损伤变量，获取t
m+1
时间步下的名义应力张量σ'
m+1
；
[0114]
s61：获取考虑时间效应-疲劳-荷载耦合损伤的应力-应变本构关系为：
[0115][0116]
式中，σ为考虑时间效应-疲劳-荷载耦合损伤的应力；εe为弹性应变；为考虑耦合损伤的刚度矩阵；is为四阶对称张量，为克罗内克积；g(dc)为损伤剪切模量；k(dc)为体积模量；其中，(is)
ijkl
＝1/2(δ
ik
δ
jl
+δ
il
δ
jk
)，i、j、k、l均为四阶对称张量的角标；δ
ik
、δ
jl
、δ
il
、δ
jk
均为四阶对称张量的分量；其中，g(dc)、k(dc)可以用初始剪切模量g0和初始体积模量k0表示：
[0117][0118]
式中，g0为初始剪切模量；k0为初始体积模量；μ为泊松比，e(dc)为在损伤变量为dc时的弹性模量；式(17)体现了时间效应-耦合损伤对岩体弹性模量的弱化作用。
[0119]
损伤不仅对弹性模量产生影响，对所述基元的hoek-brown弹塑性损伤模型的中的基元的力学参数也会产生影响。由以下公式得到hoek-brown弹塑性随着时间的变化规律：
[0120]
(mb)
t
＝(1-d(t))m
b0
ꢀꢀꢀꢀꢀ
(18a)
[0121]
(s)
t
＝(1-d(t))s0ꢀꢀꢀꢀꢀ
(18b)
[0122]
式中：(mb)
t
为t时刻的mb值；m
b0
为初始的mb值；(s)
t
为t时刻的s值；s0为初始的s值；
[0123]
对公式(5)的hoek-brown弹塑性损伤模型选用完全隐式应力回映算法，整个求解过程分为弹性预测，塑性修正和损伤修正3个步骤。通常情况下，岩石在开始发生形变时，处于弹性阶段，当形变达到一定程度时，进入塑性形变阶段，因此首先认为岩石处于弹性形变阶段，按照以下步骤求解：
[0124]
s62：获取第t
m+1
时间步下的预测应力，以对岩石是否进入塑性修正阶段进行判断；
[0125]
假设则由给定的应变增量
△
ε获得预测应力表达式为：
[0126][0127]
式中，为第t
m+1
时间步下的预测应力；为考虑耦合损伤的刚度矩阵；εm为第tm时间步下的应变；
△
εm为第tm时间步下的给定的应变增量；为第t
m+1
时间步下的预测耦合损伤变量；(dc)m为第tm时间步下的耦合损伤变量；:表示双点乘，为张量运算中的符号，代表分量按照一定顺序相乘；
[0128]
s63：将第t
m+1
时间步下的预测应力代入公式(5)的屈服函数f中进行判断，若有:
[0129]
则第t
m+1
时间步下的应力σ
m+1
等于第t
m+1
时间步下的预测应力，即执行s62；否则执行s64，进入塑性修正阶段。
[0130]
s64：根据塑性势函数，对hoek-brown弹塑性损伤模型进行塑性修正，获取修正后的等效塑性应变；
[0131]
为避免t
n+1
时间步结束时的计算应力偏离屈服面，本实施例采用完全隐式的应力回映算法，将弹性预测应力修正回屈服面。塑性修正过程中，保持第tm时间步下的应变增量
△
ε
m+1
及损伤变量不变，塑性状态下，应变增量中包含弹性应变增量
△
εe和塑性应变增量
△
ε
p
，塑性应变增量表达式为：
[0132][0133]
式中：
△
ε
p
为塑性应变增量；
△
λb为对应第b个进入塑性曲面的塑性因子，因为公式(5)、(6)在主应力空间是曲面的六棱锥，有6个曲面，式(21)中，q为应力超出曲面的个数；gb是对应第b个进入塑性曲面的塑性是函数，即公式(6)；对应当返回应力位于屈服面、棱线或尖点处时，分别需要求解1个、2个或3个塑性因子；
[0134]
由于hoek-brown弹塑性损伤模型存在数值上的奇异点，无法对其直接进行求导，而采用圆滑处理的方法又会降低求解精度，因此本实施例根据所在应力空间位置分情况进行讨论，判定最终的更新应力位于尖点、屈服面或是棱线处，避免全局求导所带来的困难。此步骤为领域内公知的有限元算法技术，这里不进行详细讨论。
[0135]
s65：根据修正后的等效塑性应变，获取t
m+1
时间步下的名义应力张量σ'
m+1
；将塑性修正步所得的等效塑性应变根据式(7)和(15)中，计算t
m+1
时间步下的耦合损伤变量(dc)
m+1
，并利用式(22)对上一步所得应力再次进行修正。
[0136]
因此，t
m+1
时间步下的名义应力张量σ'
m+1
获取如下：
[0137][0138]
式中：和分别是第tm和t
m+1
时间步下的损伤弹性矩阵；为第t
m+1
时间步的塑性应变增量；(dc)
m+1
是第t
m+1
时间步下的耦合损伤变量；(dc)m为第tm时间步下的耦合损伤变量；
[0139]
s7：根据第t
m+1
时间步下的名义应力张量σ'
m+1
，获取第t
m+1
时间步下的耦合损伤变量值。
[0140]
具体的，反复执行s6进行迭代求解，输出第t
m+1
时间步下的耦合损伤变量值。为保证整体迭代求解过程中的二次收敛速率，在应力求解完成后还应给出该模型的一致切线模量。本实施例采用fortran语言对上述数值计算过程进行编写，实现了对岩石hoek-brown弹塑性损伤模型的有限元求解。并将各个基元求解的损伤值结果输出。此为现有技术，这里不进行展开说明。
[0141]
s8：根据第t
m+1
时间步下的耦合损伤变量值，对岩石的稳定性进行分析。具体的为，
根据每个时间步下的耦合损伤变量值，能够得到耦合损伤变量随时间的变化趋势，当相邻两个时间步下的耦合损伤变量值小于经验阈值时，认为耦合损伤变量值随时间变化趋势趋于平稳，此时岩石的稳定性较好，这是领域内的公知技术，这里不进行详细展开。
[0142]
本节以二维岩石单轴压缩数值试验为例，对考虑时间效应的岩石h-b细观弹塑性损伤模型即hoek-brown弹塑性损伤模型的数值求解过程进行介绍。首先按照标准试件尺寸建立二维平面计算模型如图5所示。模型高0.01m，宽0.005m，其底部施加固定边界条件，顶部施加位移边界条件。利用四边形solid单元对模型进行网格划分，共划分为1250个单元。
[0143]
以mb为例，说明各基元的参数赋值过程。在平均值m
b0
不变的情况下，由于形状参数ζ的改变，mb在所研究空间内的分布方式完全不同，具有细观上的无序性。这种无序性体现了岩石独特的离散性特征。采用monte-carlo方法(蒙特卡洛方法)对在物理空间中具有无序性的强度参数mb进行赋值。首先在(0,1)范围内生成c个(c为总的基元个数)服从均匀分布的随机数u
l
，再利用逆变法将这组随机数映射为服从weibull分布的随机数，映射函数如式(23)：
[0144]
mb＝m
b0
(-ln(u
l
))
1/ζ
ꢀꢀꢀꢀꢀꢀꢀꢀ
(23)
[0145]
对各划分单元材料(基元)进行随机属性赋值，首先根据本技术的方法通过编写程序生成c个服从weibull分布的随机数，然后通过python语言中的循环语句生成c个材料数组(material)，最后对c个基元进行遍历赋值。
[0146]
计算对象为片麻岩，除已给出的计算参数外，利用许宏发提出的参数识别方法，获得时间效应计算参数p＝0.7524e-3，q＝-0.5416e-3，r＝0.028。表1为考虑时间效应的h-b细观弹塑性损伤模型计算参数。
[0147]
表1细观hb-d模型计算参数
[0148][0149]
模型细观参数影响性分析：设m
b0
为固定值6.36，则不同ζ值下mb的分布规律如图6a-图6f所示。可以看出，ζ值越小，mb的离散性越大，其最大值偏离平均值的幅度也就越大。当ζ＝0.5时，试件rve单元(基元)中mb的最大值为1.97，为平均值的9.85倍；当ζ＝10时，试件rve单元中mb的最大值为0.22，为平均值的1.1倍。还也可以看出，ζ值越大，岩石内部各rve单元中的mb值仅在平均值附近有小范围浮动，说明岩石试样的均质性较为明显，而ζ值越小，mb值在平均值附近浮动较大，说明岩石试样具有较强的非均质性。
[0150]
图7a-图7f统计了不同ζ值下，各mb值区域之间的单元个数。由图7可以看出，当ζ＝0.5时，整体模型中mb值接近均值0.2的rve个数仅为133个，占总单元数的10％左右。同时，mb值的浮动范围较大，表明整体模型的非均质性十分明显。ζ值越大，各rve单元上的mb值越接近于均值。当ζ＝4时，mb值小于0.01和大于0.3的rve单元数已经为零。当ζ＝10时，各rve单元上的mb值只分布在0.1和0.3之间，岩石试件整体力学性质呈现出明显的均质性。图7从不同区域间rve单元数量的角度说明了ζ值对岩石均质性的影响，即ζ值越大，岩石的均质性越明显，反之，岩石呈现出强烈的非均质性。
[0151]
在模型顶部施加位移边界条件，记录不同荷载步下mb值及损伤值的演化规律，结果如图8a和图8b所示。可以看出，在加载过程中，损伤值逐渐累积，力学参数较弱的rv e单
元较早进入塑性区，其强度参数mb在损伤的作用下不断发生弱化，最终形成破裂损伤带。由于形状参数ζ的不同，mb的分布集度不同。当ζ较小时(ζ＝0.5)，模拟试样呈现出强烈的非均质性，最终形成的损伤带也与均质条件下有着明显的差异。随着ζ值的增大，试样中各个rve单元中的mb值在平均值m
b0
附近的波动范围减小，模拟试样的非均质性逐渐减弱。此时试样在荷载作用下产生的损伤带形态越趋近于完全均质条件下的损伤带形态。由计算结果可以看出，本技术的考虑时间效应的岩石hoek-brown弹塑性损伤模型能够较好地描述岩石的非均质特性，再现加载岩石的破裂过程。
[0152]
本实施例中的时间及形状参数对局部化影响分析：利用三轴蠕变试验对形状参数ζ＝10条件(试件均质性较强)所建的时效模型进行验证。三轴蠕变试验过程中采用定围压，围压保持5mpa固定不变，轴压初始值为10mpa，当蠕变变形速率小于0.001mm/d，进行下一级加载，每级加载20mpa。数值计算采用同样的加载方式，根据的参数识别方法，以片麻岩蠕变试验数据为例，获得时效计算参数p＝0.8245e-3
，q＝-0.4731e-3
，r＝0.025，h-b模型的参数获取方法前述相同，具体见表1。
[0153]
最终获得计算值与试验值对比结果如图9所示。从图9中可以明显地看出，当蠕变进入稳态蠕变阶段，没有明显的损伤出现。在较低的应力水平下，内部结构的细观损伤积累，还不能导致试样的宏观破坏，因此，逐渐趋于稳定状态；在较高的应力水平下，试样的宏观力学特性逐步劣化。细观损伤积累到一定程度时，加速蠕变阶段开始，试样的宏观力学特性以加速蠕变的形式表现出来。计算结果与试验结果较接近，说明了所建模型的正确性。
[0154]
图10a和图10b，在ζ＝6、8条件下，蠕变加速段岩样损伤值随时间的变化规律。由图可以看出，在最后一级的荷载作用下，岩样内部已经产生荷载损伤。当ζ＝6时，岩样内部最大损伤值为8.457e-2
，随着时间的变化，岩样的力学参数弱化，在荷载不变的情况下，导致损伤值不断累积增长(逐渐增长为1.115e-1
和1.504e-1
)，最终会发生破坏。当ζ＝8时，其变化规律类似，但是损伤带形状发生了变化，均质性越强的试样(ζ值越大)，其损伤带呈现出的对称性越明显。
[0155]
寒区边坡工程长期稳定性评价:
[0156]
本技术对受冻融影响的北方某边坡的长期稳定性展开数值模拟研究。所选断面k400+184位于集安与白山市交界处的丹东至阿勒泰高速公路沿线。边坡临近林峰水库。研究断面主要岩石组成成分为强风化石英岩，节理裂隙较发育，岩体较破碎，结构面结合好，块状结构，岩体基本质量等级为ⅳ类。边坡所在的地理位置如图11所示。
[0157]
根据工程现场情况建立边坡有限元计算模型如图11所示。模型长36m，高22.6m，共划分为4786个节点和4610个单元。模型两侧施加水平约束，底部施加固定约束。根据地勘报告及室内试验初步确定计算参数的均值如表2所示。假设岩石的弹性模量、单轴抗压强度σ
ci
，mb和s值服从weibull分布，考虑形状参数变化和时间变化对边坡安全系数和损伤模式进行研究。其中，当ζ＝10时，对应的参数可以视为等同于均值情况。
[0158]
表2细观hb-d模型计算参数表
[0159]
[0160]
形状参数对边坡稳定性影响分析:图12a-图12c给出了不同形状参数ζ下弹性模量(为简洁起见，只给出了弹性模量的分布规律)以及边坡损伤带的分布云图。由图可见，在选取参数范围内(ζ＝2至ζ＝10)，形状参数ζ对边坡最大损伤值的影响较小，损伤变化幅度的数量级在1
×
10-4
左右。但形状参数ζ对损伤分布形态的影响较大，ζ值越小，边坡损伤带较均质条件下出现越多的“枝杈”状，随着ζ值的增大，边坡损伤带逐渐集中、光滑，与均质条件下的形态更为接近。计算结果表明，在一定范围内，weibull参数取值对边坡整体安全性的影响较小，但对其破坏形态有着较大的影响。
[0161]
时效对边坡工程稳定性影响分析:先不考虑冻融作用，仅进行时效对边坡稳定性影响的分析。图13a-图13f给出了开挖后边坡最大损伤值随时间变化规律。由计算结果可以看出，随着时间的推移，岩体参数不断发生弱化，向长期强度值不断逼近，因此其损伤值也发生明显增长。开挖初期，边坡内部最大损伤值为0.27，随着时间的变化，最大损伤值不断增长。当时间t为20月时，边坡内部最大损伤值达到0.2855，较初始值增长了0.0155。损伤值的变化速率在初期较为明显，在后期逐渐趋近于零。
[0162]
图14为边坡安全系数随时间的变化规律。与位移和损伤值的变化规律相似，在初期，安全系数下降较为明显，在前12个月中，安全系数由1.59下降至1.27，在12个月以后，安全系数值逐渐稳定在1.26，不再发生明显的变化。
[0163]
时效及冻融疲劳对边坡工程稳定性影响分析:据现场调研，所选断面长期受到冻融疲劳循环的扰动作用。因此，分析其在冻融疲劳循环作用下的长期安全储备是十分有必要的。由于岩石在冻融与荷载的共同作用下会产生耦合损伤，而冻融损伤值与冻融次数有关，其关系可以表达为：
[0164][0165]
式中，dr为疲劳损伤变量，在本实施例中为冻融疲劳损伤变量；n在本实施例中为冻融疲劳次数。
[0166]
计算过程中认为模型中的参数e、mb和s值随时间(月)发生变化，时间效应相关参数与前述相同。图15为冻融10次，冻融40次和冻融80次后边坡最大损伤值随时间变化规律，计算时假设整个边坡模型都受冻融疲劳损伤影响。由计算结果可以看出，当冻融次数n＝10时，开挖初期即t＝0月时边坡内部最大耦合损伤值为0.283。随着时间的推移，岩体参数不断发生弱化，向长期强度值不断逼近，因此其损伤值也发生明显增长。当时间t为20月时，边坡内部最大耦合损伤值达到0.3035，较初始值增长了0.0205。当冻融次数n＝40时，在t＝0月时边坡内部最大耦合损伤值为0.3031，当t＝20月时边坡内部最大耦合损伤值增长为0.3136。当冻融次数n＝80时，在t＝0月时边坡内部最大耦合损伤值为0.3214，当t＝20月时边坡内部最大耦合损伤值增长为0.328。耦合损伤值的变化速率在初期较为明显，在后期逐渐趋近与零。
[0167]
图16a-图16f为冻融80次后，边坡位移云图随时间的变化规律。由计算结果可以看出，边坡开挖完成后，其位移也随着时间会发生一定的变化。以最大位移值作为标准，当t＝0月时，边坡最大位移值为0.016m，当t＝20月时，边坡最大位移值为0.023m，增长了约41.8％。图17为坡面各监测点处位移随时间变化规律。由图可以看出，开挖初期边坡位移变化较为明显，随着时间的推移，边坡位移逐渐趋于稳定，不在随时间发生变化。以位移值最
大的2号监测点为例，当t＝0、4、8、12、16和20月时，该点位移为0.015m，0.018m、0.02m、0.021m、0.0214m和0.216m，变化幅度逐渐减小，最终趋近为零。
[0168]
现场采用预应力锚索对边坡进行加固，以保证边坡的长期稳定性。本工程的锚索长度为15～18m，锚固段长度均为8m，每孔锚索6束钢绞线，锚索倾角为20
°
。锚索预应力设计为600kn。为合理计算锚索预应力，在上述模型的基础上，利用embedded region方法计算土体与锚索之间的接触问题。通过温度场法设置锚索预应力，预应力计算公式为：
[0169][0170]
式中，δt为预设降低温度值，单位为℃；fn为锚索的预应力值，单位为牛；e为材料弹性模量，单位为pa；a为截面面积，单位为m2；α
′
为锚索的热膨胀系数，单位为1/℃。
[0171]
分别设置锚索预应力为200kn，300kn，400kn，500kn和600kn，计算得到边坡最终损伤区分布如图18a-图18d所示。
[0172]
由计算结果可以看出，当施加预应力锚索后，边坡损伤带向坡内移动，减小了坡面发生滑落的潜在风险。随着预应力的增加，边坡损伤值逐渐减小。与不施加预应力锚索时相比，当施加了200kn的锚索时，边坡损伤值发生明显的下降，最大损伤值由0.328下降为0.202。当预应力为400kn以上时，边坡损伤值不再发生明显变化。根据计算结果，并考虑到施工成本问题，建议该边坡中锚索预应力设置为400kn为宜。
[0173]
本发明的一种考虑流变和疲劳作用的岩体稳定性计算方法，能够对冻融疲劳扰动下岩土工程的长期稳定性进行评价，为特殊地区工程的施工参数设计提供方法依据。从不同形状参数下，试件压缩数值试验中得到的岩样破坏带，以及边坡工程稳定计算中的损伤带分布情况来看，该方法是具有合理性。
[0174]
最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。

技术特征：
1.一种考虑疲劳和时间效应的岩体稳定性计算方法，其特征在于，包括如下步骤：s1：将岩石划分为若干基元v，以获取所述基元的力学参数；所述基元的力学参数包括反映岩体特征的经验参数m
b
、岩体破碎及完整程度参数s、岩体质量参数a和损伤前的弹性模量e0；s2：根据所述基元的力学参数，获取所述基元的hoek-brown弹塑性损伤模型；s3：获取基元的时效损伤变量d(t)；s4：获取基元的疲劳损伤变量d
r
；s5：根据荷载损伤变量d
m
、所述时效损伤变量d(t)和所述疲劳损伤变量d
r
，获取耦合损伤变量dc：s6：根据所述基元的hoek-brown弹塑性损伤模型和所述耦合损伤变量，获取t
m+1
时间步下的名义应力张量σ'
m+1
；s7：根据第t
m+1
时间步下的名义应力张量σ'
m+1
，获取第t
m+1
时间步下的耦合损伤变量值；s8：根据第t
m+1
时间步下的耦合损伤变量值，对岩石的稳定性进行分析。2.根据权利要求1所述的一种考虑疲劳和时间效应的岩体稳定性计算方法，其特征在于，所述s1中，基元的力学参数获取方法如下：p
z
＝p0·
x
z
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)式中，p0为岩石试样的宏观性质；x
z
为weibull随机数；p
z
为基元的力学参数；其中，weibull随机数x
z
获取方法如下：获取方法如下：式中：w(x)为weibull分布的累积分布函数；w(x)为weibull分布的概率密度函数；x0>0是比例参数；ζ为形状参数；x为随机变量；x
z
＝x0(-ln(u
m
))
1/ζ
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)式中：u
m
为均匀分布的随机数。3.根据权利要求1所述的一种考虑疲劳和时间效应的岩体稳定性计算方法，其特征在于，所述s2中，所述基元的hoek-brown弹塑性损伤模型获取如下：brown弹塑性损伤模型获取如下：式中，f为屈服函数；j2为第二偏应力不变量；θ
σ
为罗德角；d
m
为荷载损伤变量；σ
ci
为完整岩石的单轴抗压强度；p为静水压力；m
bg
为m
b
对应的塑性势函数参数；s
g
为s对应的塑性势函数参数；a
g
为a对应的塑性势函数参数；g为塑性势函数；
式中，α为决定损伤后岩石材料软化曲线的初始斜率，取值范围为(0,+∞)；β为决定岩石最大损伤值的参数，取值范围为[0,1)；为等效塑性应变；式中，ε
p1
为空间坐标系x轴方向上的主塑性应变；ε
p2
为空间坐标系y轴方向上的主塑性应变；ε
p3
为空间坐标系z轴方向上的主塑性应变。4.根据权利要求1所述的一种考虑疲劳和时间效应的岩体稳定性计算方法，其特征在于，所述s3中，获取基元的时效损伤变量d(t)方法如下；式中，ε为应变；σ0为轴向的应力；k1、k2均为弹性元件的力学参数；η为黏壶的粘滞系数；t为流变的时间；化简式(9)可得ε＝σ0[p+qexp(-rt)]
ꢀꢀꢀꢀꢀꢀꢀꢀ
(10)式中，p、q、r均为中间参数，其中，从公式(10)可以得到：e(t)＝p+qexp(-rt)
ꢀꢀꢀꢀꢀꢀꢀꢀ
(11)e(t)为弹性模量随着时间t变化的函数；式中，d(t)为时效损伤变量，即由时间效应引起的损伤变量；e0为损伤前的弹性模量。5.根据权利要求1所述的一种考虑疲劳和时间效应的岩体稳定性计算方法，其特征在于，所述s4中，基元的疲劳损伤变量d
r
获取如下：式中：d
t
为第n
t
次疲劳扰动作用下岩石断面的损伤值；d
r
为疲劳损伤变量d
r
，即第n
t
次疲劳扰动之前的n次疲劳损伤变量；n
t
为疲劳扰动次数的编号；n为第n
t
次疲劳扰动之前的疲劳扰动总次数；w和u均为材料参数；式中，n0为岩石疲劳作用前的孔隙率；n为岩石疲劳作用后的孔隙率；v
p
为岩石疲劳作用后的纵波波速。6.根据权利要求1所述的一种考虑疲劳和时间效应的岩体稳定性计算方法，其特征在于，所述s5中，
d
c
＝d(t)+d
m
+d
r-d(t)d
m-d
r
d
m-d
r
d
t
+d(t)d
m
d
r
ꢀꢀꢀꢀꢀ
(15)式中，dc为耦合损伤变量；d
m
为荷载损伤变量。7.根据权利要求1所述的一种考虑疲劳和时间效应的岩体稳定性计算方法，其特征在于，所述s6中，获取t
m+1
时间步下的名义应力张量σ'
m+1
的方法如下：s61：获取考虑时间效应-疲劳-荷载耦合损伤的应力-应变本构关系为：式中，σ为考虑时间效应-疲劳-荷载耦合损伤的应力；ε
e
为弹性应变；为考虑耦合损伤的刚度矩阵；i
s
为四阶对称张量，为克罗内克积；g(d
c
)为损伤剪切模量；k(d
c
)为体积模量；其中，式中，g0为初始剪切模量；k0为初始体积模量；μ为泊松比，e(d
c
)为在损伤变量为d
c
时的弹性模量；(m
b
)
t
＝(1-d(t))m
b0
ꢀꢀꢀꢀ
(18a)(s)
t
＝(1-d(t))s0ꢀꢀꢀꢀ
(18b)式中：(m
b
)
t
为t时刻的m
b
值；m
b0
为初始的m
b
值；(s)
t
为t时刻的s值；s0为初始的s值；s62：获取第t
m+1
时间步下的预测应力，以对岩石是否进入塑性修正阶段进行判断；假设则由给定的应变增量δε获得预测应力表达式为：式中，为第t
m+1
时间步下的预测应力；为考虑耦合损伤的刚度矩阵；ε
m
为第t
m
时间步下的应变；δε
m
为第t
m
时间步下的给定的应变增量；为第t
m+1
时间步下的预测耦合损伤变量；(d
c
)
m
为第t
m
时间步下的耦合损伤变量；:表示双点乘；s63：根据屈服函数f，若:则第t
m+1
时间步下的应力σ
m+1
等于第t
m+1
时间步下的预测应力，即充分执行s62；否则执行s64；s64：根据塑性势函数，对hoek-brown弹塑性损伤模型进行塑性修正，获取修正后的等效塑性应变；s65：根据修正后的等效塑性应变，获取t
m+1
时间步下的名义应力张量σ'
m+1
如下：式中：和分别是第t
m
和t
m+1
时间步下的损伤弹性矩阵；为第t
m+1
时间步的塑性应变增量；(d
c
)
m+1
是第t
m+1
时间步下的耦合损伤变量；(d
c
)
m
为第t
m
时间步下的耦合损伤变量。

技术总结
本发明提供一种考虑疲劳和时间效应的岩体稳定性计算方法，包括：将岩石划分为若干基元，以获取基元的力学参数；根获取所述基元的Hoek-Brown弹塑性损伤模型；获取基元的时效损伤变量；获取基元的疲劳损伤变量；获取耦合损伤变量：获取t


技术研发人员：姜涛 詹涛 姜谙男 吴招锋 陈登开 姚元 万友生 单生彪 叶帅
受保护的技术使用者：南昌轨道交通集团有限公司地铁项目管理分公司
技术研发日：2023.02.08
技术公布日：2023/7/26