一种土石混合体边坡的有限元可靠度评价方法

1.本发明涉及土石混合体边坡稳定技术领域，具体地指一种土石混合体边坡的有限元可靠度评价方法。


背景技术：

2.土石混合体作为典型的土石二元介质，是工程界一种常见的不稳定复杂地质体。由于土石混合体中土
‑
石物理力学特性的巨大差异，土石接触间复杂的作用关系，使得土石混合体边坡与传统的均质边坡有着较大的差异，采用传统的边坡安全稳定难以对土石混合体边坡的稳定性准确评价。
3.目前，针对土石混合体边坡的研究主要是通过构建各种非均质模型来模拟力学性质巨大差异的土石结构特征，再结合常用的刚体极限平衡法和有限元强度折减法确定性分析法对其进行边坡稳定性评价。然而，对于土石混合体边坡稳定而言，土石之间力学作用复杂导致危险滑动面复杂多变，而刚体极限平衡法难以确定危险滑动面的位置，无法反映土石之间复杂力学作用对边坡变形破坏特征的影响；有限元强度折减法能通过强度折减迭代求解出最危险滑动面，但迭代计算过程中有限元计算量相对较大。与此同时，由于坡体内块石含量及空间分布的不同，导致其力学特性及变形特征表现出了很强的不确定性，因此，对土石混合体边坡进行稳定性分析时，必须同时考虑岩土体材料的非均质性，可利用可靠度方法用明确的、定量的概率来定性评价边坡稳定性。
4.对于土石混合边坡稳定的可靠度而言，可将土石混合体不确定的力学参数作为基本随机变量进行统计样本分析，采用蒙特卡洛(monte carlo)模拟法作为可靠度模拟方法来计算土石混合体边坡的失效概率，其理论简洁明了，程序结构简单，易于操作，且能不受随机变量分布形式和复杂结构功能函数的影响，通过随机的样本去模拟逼近结构的真实状态的功能函数，但monte carlo模拟法需要随机生成大数据量级的样本，特别在与有限元强度折减相结合时，其计算量大的缺点就显得更为突出，且大数量级的有限元强度折减法不利于程序的自动化编制。


技术实现要素：

5.本发明的目的在于克服上述不足，提供一种土石混合体边坡的有限元可靠度评价方法，有效减低土石混合体边坡可靠度评估过程中的计算工作量，便于程序的自动化编制。
6.本发明为解决上述技术问题，所采用的技术方案是：一种土石混合体边坡的有限元可靠度评价方法，它包括以下步骤：
7.步骤1)：根据土体和块石的不确定性参数及其随机分布特征，构建相应的联合概率密度函数f(x)，并以此随机生成n组不确定土体及块石材料参数的随机样本；
8.步骤2)：结合坡体的实际级配和含石量，利用随机多边形块石模型构建土石混合体边坡有限元模型；
9.步骤3)：根据土与块石不确定性参数的均值和确定性参数，利用有限元强度折减
法对土石混合体边坡有限元模型进行稳定性分析，确定边坡极限失稳条件下的危险滑动面及其贯穿特征点的位移u1；
10.步骤4)：将生成n组不确定土体和块石材料参数的随机样本和确定参数作为输入，代入该土石混合体边坡的有限元模型，进行n次重复的边坡稳定有限元数值模拟，得到n组不同岩土体材料参数条件下边坡滑动面同一贯穿特征点的位移u2；
11.步骤5)：将步骤4)和步骤3)的滑动面贯穿特征点的位移进行比较，计算大于极限失稳条件下边坡滑动面贯穿特征点位移的失效样本数目，计算得到边坡的失效概率。
12.优选地，所述步骤1)中，根据n维不确定性参数随机变量x(x1,x2,l,x
n
)间的相关性矩阵ρ＝(ρ
ij
)
n
×
n
，该联合概率密度函数f(x)分为3种情况：
13.一、若不确定性参数间完全独立，则联合概率密度函数为各单一参数概率密度的积，可表示为：
[0014][0015]
其中，x
i
为第i个不确定性参数；
[0016]
二、若不确定性参数间完全相关，则利用nataf变换法根据标准正态概率密度函数、累计分布函数和相关系数矩阵生成相应的联合概率密度函数；具体如下：
[0017]
首先，将随机变量x(x1,x2,l,x
n
)按等概率变换原则转换为标准随机变量y＝(y1,y2,l,y
n
)：
[0018]
x
i
＝f
‑1[φ(y
i
)]
ꢀꢀꢀꢀ
(2)
[0019]
式中，f
‑1(
·
)为累计分布函数f(
·
)的反函数；
[0020]
然后利用高斯copula函数来生成相应的n维联合概率密度函数：
[0021][0022]
式中，ρ0为标准随机变量y的相关系数矩阵；det(ρ0)为矩阵行列式；
[0023]
最后，根据nataf变换法原理可通过各不确定性参数概率密度函数f(
·
)和标准正态累计分布函数φ(
·
)推求不确定性参数随机变量x的联合概率密度函数：
[0024][0025]
其中随机变量(x
i
,x
j
)的相关系数ρ
ij
和标准随机变量(y
i
,y
j
)相关系数矩阵ρ
0ij
间的关系为：
[0026][0027]
式中，φ2(
·
)为2维联合概率密度函数，和分别为随机变量x
i
和x
j
的均值和标准差；
[0028]
三、不确定性参数间部分相关，如土或块石的不确性参数相关，而土与块石间不确定性参数相互独立，此时先求出q个独立参数随机变量(x1,x2,l,x
q
)的联合概率密度函数再利用nataf变换得到相关不确定性参数随机变量(x
q+1
,l,x
n
)的联合概率密度函
数f(x
q+1
,l,x
n
)，故此时的联合概率密度函数为
[0029]
优选地，所述步骤2)中的随机多边形块石模型通过生成随机多边形来模拟边坡体内的随机块石。
[0030]
优选地，所述步骤2)中，在不同粒径及级配区间生成不同直径的随机外接圆d
max
和随机内接圆d
min
，再通过调整畸变系数μ改变在外接随机圆和内接随机圆间生成不同几何形态的随机多边形块石；以圆心(x0，y0)等分为n份，通过调整每等分角度(θ1,θ2,l,θ
n
)的畸变系数来重构任一直径d
n
多边形形状特征块石的顶点坐标：
[0031]
d
min
＝(1
‑
μ)d
max
ꢀꢀꢀꢀꢀ
(6)
[0032]
d
n
＝(1
‑
μrand(1))d
max
ꢀꢀ
(7)
[0033]
(x0+0.5d
n cosθ
n
,y0+0.5d
n sinθ
n
)
ꢀꢀ
(8)
[0034]
式中，rand(1)能生成0～1的随机数；
[0035]
通过计算任两块石间圆心的直线距离与两个随机外接圆直径大小来判定两块石是否入侵，实现土石混合体边坡的入侵判定；同时，含石量p是通过计算所有块石所占面积与该区域总面积的质量百分比：
[0036][0037]
式中：a
i
为第i个块石所占面积；l为块石总数目，a是区域面积。
[0038]
优选地，所述步骤3)中，有限元强度折减法是通过将抗剪强度参数折减至使临界滑动状态时土石混合体边坡稳定性的有限元分析，土石混合体中土体和石块的抗剪强度参数粘聚力和内摩擦角均采用一样的折减系数：
[0039]
c
f
＝c/f
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(10)
[0040][0041]
式中，c、c
f
分别为实际和折减后的粘聚力；分别为实际和折减后的内摩擦角；f为临界状态下折减系数。
[0042]
优选地，所述步骤3)中的边坡有限元模型包含边坡各土层的几何坐标特征，岩土体材料采用的本构模型、土体和块石参数、相关荷载和边界条件。
[0043]
优选地，所述土体和块石参数包含黏聚力、内摩擦角、密度、饱和渗透系数等，其中部分计算参数根据实际相关统计分析被定义为不确定性参数，服从均匀分布、正态分布和对数正态分布随机分布特征，并得到相应的均值和标准差；所述边坡滑动面贯穿特征点指在边坡稳定有限元强度折减计算后贯穿失效滑动面上的某个有限元节点。
[0044]
优选地，所述步骤4)中，将生成n组不确定土体和块石参数的随机样本和确定参数作为输入，代入该土石混合体边坡的有限元模型，采用mohr
‑
coulomb弹塑性本构模型，进行n次重复的边坡稳定有限元数值模拟，得到n组不同岩土体材料参数条件下边坡滑动面同一贯穿特征点的位移u2；
[0045]
其中，mohr
‑
coulomb弹塑性本构模型表示为：
[0046][0047]
式中，c和φ分别为黏聚力和内摩擦角；σ
m
为平均应力；j2为应力偏量第2不变量；
其中θ为lode角。
[0048]
优选地，所述步骤5)具体为：
[0049]
将步骤4)中滑动面同一贯穿特征点位移u2与步骤3)中滑动面贯穿特征点位移u1进行比较，计算不确定参数x
i
条件下边坡的功能函数g(x
i
)＝|u1|
‑
|u2|；当g(x
i
)≥0时，表明边坡稳定；当g(x
i
)<0时，表明边坡失稳，此时不确定参数x
i
为失效样本，并统计失效样本的数目n，最后计算得到边坡的失效概率：
[0050][0051]
式中i[g(x
i
)]为指示函数，当g(x
i
)<0时，i[g(x
i
)]＝1；当g(x
i
)≥0时，i[g(x
i
)]＝0。
[0052]
本发明的有益效果：与现有技术相比，本发明的土石混合体边坡有限元可靠度评价方法具有以下优势：
[0053]
1.揭示土石结构特征对土石混合体边坡稳定性的影响：本发明利用随机多边形块石模型构建土石混合体的边坡模型，能充分反映不同含石量、块石空间分布及土石复杂力学作用对边坡稳定性的影响，更适用于土石混合体边坡稳定性的影响。
[0054]
2.计算效率高：传统的蒙特卡洛卡洛有限元模拟方法在计算过程中，需结合蒙特卡洛模拟进行大数量级边坡稳定的有限元强度折减法分析，计算效率较低，特别是有限元强度折减法中抗剪强度参数的迭代计算，对计算资源要求高。本发明仅需通过一次有限元强度折减求出土石混合体边坡极限失稳条件下的危险滑动面，通过对比不同力学参数及极限失稳条件下同一贯穿特征点的位移大小，有效减少了蒙特卡洛大数量级有限元计算过程中力学参数折减计算的工作量。
[0055]
3.计算理论清晰，便于程序的自动化编制，为土石混合体边坡工程设计提供参考：本发明基于蒙特卡洛卡洛模拟法，避免了有限元强度折减法中边坡极限状态的人工设定，计算理论清晰、计算过程简单，有利于程序的自动化编制，扩展了有限元在土石混合体边坡稳定可靠度评价的应用。
[0056]
4、本发明能充分反映土石结构特征对土石混合体边坡稳定性的影响，计算理论清晰、计算过程简单，有效减少了蒙特卡洛大数量级有限元计算过程中力学参数折减计算的工作量，便于程序的自动化编制，大大扩展了有限元在土石混合体边坡稳定可靠度评价的应用。
附图说明
[0057]
图1为本发明的具体计算流程图；
[0058]
图2为本发明随机块石的生成示意图；
[0059]
图3为本发明实例的土石混合体边坡尺寸；
[0060]
图4为本发明实例的土石混合体边坡模型；
[0061]
图5为本发明实例边坡滑动面及贯穿特征点；
[0062]
图6为本发明实例不同土体黏聚力和摩擦角条件下滑动面贯穿特征点a的位移。
具体实施方式
[0063]
下面结合附图和具体实施例对本发明作进一步的详细描述。
[0064]
如图1至6所示，一种土石混合体边坡的有限元可靠度评价方法，它包括以下步骤：
[0065]
步骤1)：根据土体和块石的不确定性参数及其随机分布特征，构建相应的联合概率密度函数f(x)，并以此随机生成n组不确定土体及块石材料参数的随机样本；
[0066]
步骤2)：结合坡体的实际级配和含石量，利用随机多边形块石模型构建土石混合体边坡有限元模型；
[0067]
步骤3)：根据土与块石不确定性参数的均值和确定性参数，利用有限元强度折减法对土石混合体边坡有限元模型进行稳定性分析，确定边坡极限失稳条件下的危险滑动面及其贯穿特征点的位移u1；
[0068]
步骤4)：将生成n组不确定土体和块石材料参数的随机样本和确定参数作为输入，代入该土石混合体边坡的有限元模型，进行n次重复的边坡稳定有限元数值模拟，得到n组不同岩土体材料参数条件下边坡滑动面同一贯穿特征点的位移u2；
[0069]
步骤5)：将步骤4)和步骤3)的滑动面贯穿特征点的位移进行比较，计算大于极限失稳条件下边坡滑动面贯穿特征点位移的失效样本数目，计算得到边坡的失效概率。
[0070]
优选地，所述步骤1)中，根据n维不确定性参数随机变量x(x1,x2,l,x
n
)间的相关性矩阵ρ＝(ρ
ij
)
n
×
n
，该联合概率密度函数f(x)分为3种情况：
[0071]
一、若不确定性参数间完全独立，则联合概率密度函数为各单一参数概率密度的积，可表示为：
[0072][0073]
其中，x
i
为第i个不确定性参数；
[0074]
二、若不确定性参数间完全相关，则利用nataf变换法根据标准正态概率密度函数、累计分布函数和相关系数矩阵生成相应的联合概率密度函数；具体如下：
[0075]
首先，将随机变量x(x1,x2,l,x
n
)按等概率变换原则转换为标准随机变量y＝(y1,y2,l,y
n
)：
[0076]
x
i
＝f
‑1[φ(y
i
)]
ꢀꢀꢀꢀ
(2)
[0077]
式中，f
‑1(
·
)为累计分布函数f(
·
)的反函数；
[0078]
然后利用高斯copula函数来生成相应的n维联合概率密度函数：
[0079][0080]
式中，ρ0为标准随机变量y的相关系数矩阵；det(ρ0)为矩阵行列式；
[0081]
最后，根据nataf变换法原理可通过各不确定性参数概率密度函数f(
·
)和标准正态累计分布函数φ(
·
)推求不确定性参数随机变量x的联合概率密度函数：
[0082][0083]
其中随机变量(x
i
,x
j
)的相关系数ρ
ij
和标准随机变量(y
i
,y
j
)相关系数矩阵ρ
0ij
间的关系为：
[0084][0085]
式中，φ2(
·
)为2维联合概率密度函数，和分别为随机变量x
i
和x
j
的均值和标准差；
[0086]
三、不确定性参数间部分相关，如土或块石的不确性参数相关，而土与块石间不确定性参数相互独立，此时先求出q个独立参数随机变量(x1,x2,l,x
q
)的联合概率密度函数再利用nataf变换得到相关不确定性参数随机变量(x
q+1
,l,x
n
)的联合概率密度函数f(x
q+1
,l,x
n
)，故此时的联合概率密度函数为
[0087]
优选地，所述步骤2)中的随机多边形块石模型通过生成随机多边形来模拟边坡体内的随机块石。
[0088]
优选地，所述步骤2)中，在不同粒径及级配区间生成不同直径的随机外接圆d
max
和随机内接圆d
min
，再通过调整畸变系数μ改变在外接随机圆和内接随机圆间生成不同几何形态的随机多边形块石；以圆心(x0，y0)等分为n份，通过调整每等分角度(θ1,θ2,l,θ
n
)的畸变系数来重构任一直径d
n
多边形形状特征块石的顶点坐标：
[0089]
d
min
＝(1
‑
μ)d
max
ꢀꢀꢀꢀꢀ
(6)
[0090]
d
n
＝(1
‑
μrand(1))d
max
ꢀꢀ
(7)
[0091]
(x0+0.5d
n cosθ
n
,y0+0.5d
n sinθ
n
)(8)
[0092]
式中，rand(1)能生成0～1的随机数；
[0093]
通过计算任两块石间圆心的直线距离与两个随机外接圆直径大小来判定两块石是否入侵，实现土石混合体边坡的入侵判定；同时，含石量p是通过计算所有块石所占面积与该区域总面积的质量百分比：
[0094][0095]
式中：a
i
为第i个块石所占面积；l为块石总数目，a是区域面积。
[0096]
优选地，所述步骤3)中，有限元强度折减法是通过将抗剪强度参数折减至使临界滑动状态时土石混合体边坡稳定性的有限元分析，土石混合体中土体和石块的抗剪强度参数粘聚力和内摩擦角均采用一样的折减系数：
[0097]
c
f
＝c/f
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(10)
[0098][0099]
式中，c、c
f
分别为实际和折减后的粘聚力；分别为实际和折减后的内摩擦角；f为临界状态下折减系数。
[0100]
优选地，所述步骤3)中的边坡有限元模型包含边坡各土层的几何坐标特征，岩土体材料采用的本构模型、土体和块石参数、相关荷载和边界条件。
[0101]
优选地，所述土体和块石参数包含黏聚力、内摩擦角、密度、饱和渗透系数等，其中部分计算参数根据实际相关统计分析被定义为不确定性参数，服从均匀分布、正态分布和
对数正态分布随机分布特征，并得到相应的均值和标准差；所述边坡滑动面贯穿特征点指在边坡稳定有限元强度折减计算后贯穿失效滑动面上的某个有限元节点。
[0102]
优选地，所述步骤4)中，将生成n组不确定土体和块石参数的随机样本和确定参数作为输入，代入该土石混合体边坡的有限元模型，采用mohr
‑
coulomb弹塑性本构模型，进行n次重复的边坡稳定有限元数值模拟，得到n组不同岩土体材料参数条件下边坡滑动面同一贯穿特征点的位移u2；
[0103]
其中，mohr
‑
coulomb弹塑性本构模型表示为：
[0104][0105]
式中，c和φ分别为黏聚力和内摩擦角；σ
m
为平均应力；j2为应力偏量第2不变量；其中θ为lode角。
[0106]
优选地，所述步骤5)具体为：
[0107]
将步骤4)中滑动面同一贯穿特征点位移u2与步骤3)中滑动面贯穿特征点位移u1进行比较，计算不确定参数x
i
条件下边坡的功能函数g(x
i
)＝|u1|
‑
|u2|；当g(x
i
)≥0时，表明边坡稳定；当g(x
i
)<0时，表明边坡失稳，此时不确定参数x
i
为失效样本，并统计失效样本的数目n，最后计算得到边坡的失效概率：
[0108][0109]
式中i[g(x
i
)]为指示函数，当g(x
i
)<0时，i[g(x
i
)]＝1；当g(x
i
)≥0时，i[g(x
i
)]＝0。
[0110]
实施例1：
[0111]
选取15％含石量的土石混合体边坡，具体尺寸见图3，其中坡高10.0m，左右长35m，坡角45度。
[0112]
根据medley和lindquist的研究，土
‑
石区分值d
thr
＝(0.05～0.15)lc，其中，lc为工程特征尺度，边坡的工程特征尺度为坡高。所以本实施例中可选0.5m为土
‑
石区分值，即大于0.5m的颗粒为块石，小于0.5m的颗粒为土，具体级配如表1：
[0113]
表1坡体的块石级配
[0114][0115]
在该实施例中仅考虑重力的作用，采用mohr
‑
coulomb弹塑性本构模型，土石混合体内部的“块石”与“土体”具体材料参数详见如表2所示，假定土体和块石的黏聚力和摩擦角均服从正态随机分布特征且相互独立。
[0116]
表2土石混合体的计算参数
[0117]
[0118]
具体实施流程：
[0119]
步骤1)、根据表2)土体和块石抗剪强度参数黏聚力和摩擦角均值、变异系数和随机分布特征，并构建相应的联合概率密度函数f(c
土
,c
石
,φ
土
,φ
石
)，利用蒙特卡洛法随机生成1000组土体及块石抗剪强度参数的随机样本。
[0120]
其中：f(c
土
,c
石
,φ
土
,φ
石
)＝f(c
土
)f(c
石
)f(φ
土
)f(φ
石
)
ꢀꢀꢀꢀ
(14)
[0121][0122][0123][0124][0125]
式中，μ1、μ2、μ3、μ4和σ1、σ1、σ1、σ1分别为土体黏聚力c
土
、块石黏聚力c
石
、土体摩擦角φ
土
、块石摩擦角φ
石
的均值和标准差。
[0126]
步骤2)、根据表1的坡体级配和含石量，利用随机多边形块石模型构建土石混合体边坡有限元模型(见图4)，其中含石量p＝0.15；畸变系数μ＝0.2；n＝8，即将外接圆等分为8分。
[0127]
步骤3)、根据表2土与块石抗剪强度参数均值及相关确定性计算参数，利用有限元强度折减法对土石混合体边坡模型进行稳定性分析，确定边坡极限失稳条件下的危险滑动面，选取坡脚a处为贯穿特征点(见图5)，计算其极限失稳条件下的水平位移u1＝
‑
0.149mm。其中将土体和石块抗剪强度参数黏聚力c和摩擦角φ作为折减系数，且均采用同一折减系数。设定折减系数的区间为[0.5，2]，迭代步长为0.05，最大需迭代30次。
[0128]
由图5可以看出，由于坡体内土与石力学性质的巨大差异，块石阻挡了塑性剪切带的扩展趋势，危险滑动面出现了明显“绕石”现象。
[0129]
步骤4)、将生成1000组不确定土体和块石参数抗剪强度参数的随机样本和确定性参数作为输入，代入该土石混合体边坡的有限元模型，进行1000次重复的边坡稳定有限元数值模拟，得到1000组不同材料参数条件下边坡同一坡脚特征点a处的位移u1，如图6所示。
[0130]
图6为1000组不同土体黏聚力和摩擦角条件下滑动面贯穿特征点a的位移。可以看出，贯穿特征点a处的位移随着土体黏聚力c和摩擦角φ的减小而减小。
[0131]
步骤5)、将步骤4)特征点位移u2与步骤3的位移u1(
‑
0.149mm)进行比较，计算边坡的功能函数g(x
i
)＝|u1|
‑
|u2|，即如果g(x
i
)<0，计为失效样本，统计得到失效样本的数目为95，计算得到边坡的失效概率p
f
＝95/1000＝9.5％。
[0132]
针对实施例1步骤2的有限元强度折减每次最多迭代30次，而实施例1中在岩土体抗剪强度均值时计算的强度折减系数为1.20。根据折减系数的区间[0.5，2]和迭代步长0.05，可以算出均值时有限元强度折减需进行14次有限元计算，表明实施例1共进行了1000+14＝1014次有限元模拟。
[0133]
对于常规有限元可靠度方法时，一次蒙特卡洛抽样需进行一次有限元强度折减，
则可认为1000次蒙特卡洛模拟约需进行14000次有限元模拟，表明本发明方法计算效率了提高90％以上。
[0134]
上述的实施例仅为本发明的优选技术方案，而不应视为对于本发明的限制，本技术中的实施例及实施例中的特征在不冲突的情况下，可以相互任意组合。本发明的保护范围应以权利要求记载的技术方案，包括权利要求记载的技术方案中技术特征的等同替换方案为保护范围。即在此范围内的等同替换改进，也在本发明的保护范围之内。