可靠度分析法采用可靠度指标(失效概率)代替安全系数进行边坡稳定性分析[1-4]，是一种非确定性方法，更加符合边坡岩土体的非均匀性及失稳破坏不确定性。但由于可靠度分析法的计算相对复杂，目前在岩土工程分析和设计中的应用尚处于研究和探索阶段。
目前的可靠度方法研究和应用中，蒙特卡洛模拟[5-6]应用较为广泛。已有学者基于随机场理论，使用蒙特卡洛模拟进行边坡稳定性研究。宋永东[7]运用Matlab离散随机场，利用Excel作为衔接手段，将离散后的土体强度参数导入有限差分软件Flac3D，计算边坡的稳定性；胡金政等[8]用Flac3D建模，利用Fish语言将离散场与网格单元一一对应，反复N次，进行边坡稳定性计算，最后使用Matlab读取计算结果；曹少刚[9]使用Matlab编写程序，得到能够表现土体参数空间变异性的一系列随机变量，然后使用Flac3D计算边坡的安全系数；蒋水华[10]利用有限元软件Abaqus和GeoStudio编写接口程序，计算边坡的可靠度指标；Griffths等[11]使用Fortran语言编写耦合随机场理论与边坡可靠度分析软件；王新[12]使用Matlab获取离散随机场，然后与Abaqus模型相结合进行边坡的可靠度计算；袁葳等[13]以随机场理论为基础，使用Abaqus提供的用户子程序接口编写随机有限元程序，使用Python脚本进行后期处理。
上述研究在“随机有限元程序”应用方面取得了一定的进展，但仍然存在不足。首先，使用Flac3D与Matlab计算时，对不同的土体参数进行敏感性或影响程度分析时，需要在两者之间进行数以万次反复转换，计算量较大且耗时较长；其次，调用Abaqus内核进行批处理时，并没有涉及地应力迭代过程，这将导致计算结果存在一定的误差；最后，使用GeoStudio、Abaqus与Matlab相结合时，有限元软件与编写程序所使用的语言不一致，会降低原程序的计算效率。
本文利用Abaqus脚本建模使用的Python语言编写程序，将有限元建模、随机场赋值和强度折减计算有机结合起来，进行批量自动化运算，实现高效精确的边坡可靠度分析。
Phoon等[14]指出，土体的强度是各向异性的，可以从两个方向对土体的强度进行分析。对正常固结土来说，在垂直方向上，从地表开始，土体的强度随深度加深逐渐增加；在水平方向上，土体的强度是与深度无关的随机波动量。由于非平稳随机场能够合理地模拟土体强度参数随埋深逐渐增加的特征，因此，许多学者进行了相关研究。Li等[15]通过式(1)研究了土体黏聚力随有效应力及固结比的变化关系。
式中：su、σv′和OCR分别表示土体的黏聚力、有效应力及固结比。Jiang等[16]考虑土体抗剪强度参数的随机波动量服从对数正态分布，将其随深度的变化关系表示为式(2)。
式中:su(x, z)为某位置处的土体黏聚力；t为比例因子；w(x, z)服从正态分布。Griffiths等[11, 17]、Der Kiureghian等[18]等利用式(3)建立非平稳随机场，研究边坡的失效概率。
式中:cz为某深度处的土体黏聚力；c0为由非平稳随机场得到的土体黏聚力值；μcu0为地表处黏聚力均值；ρ为比例因子；z为特定的土体深度。
相对于式(1)和式(2)，式(3)的认可度较高，因此，使用式(3)将平稳随机场转化为非平稳随机场，研究土体参数的空间变异性对边坡可靠度影响。
非平稳随机场的形成通过以下4步来实现，如图 1所示。
1) Abaqus平台模型网格划分。首先，给定边坡，在Abaqus平台下划分网格，得到各个单元所对应的初始节点及节点坐标，并将其导出。
2) Python读取数据并对单元排序。将步骤1)导出的单元重新排序，保持节点序号不变，目的是使离散后的随机场变量能够批量赋值给对应的边坡单元。
3) 平稳随机场。主要包括随机场的离散和有限元的结合，利用中心点法离散随机场，得到一系列随机变量，然后按照边坡的实际空间位置，将随机变量映射到步骤2)得到的有限元单元中。
4) 非平稳随机场。非平稳随机场与土体参数实际分布比较接近，将步骤3)得到的平稳随机场转化为非平稳随机场。
求解过程包括7步，如图 2所示。
1) Python形成初始Inp文件。将边坡的几何参数、材料信息及随机场数据写入Inp文件，该文件称为初始Inp文件。
2) Inp文件进行初次运算。该步骤的主要目的是平衡地应力。在初始Inp文件中，有施加土体重力的分析步骤。边坡在初始状态下，由于自重作用，存在与重力相平衡的应力状态，因此，在进行数值模拟时，需要在边坡开始运算之前建立相应的应力场。
3) 提取初始应力生成Rpt文件。经过步骤2)的初始运算，得到一系列包含各个单元应力的Job文件，然后使用Python编写的脚本文件，提取各个单元的内力，并生成包含各个应力提取代码的Rpt文件。
4) 地应力平衡的Csv文件。在Abaqus平台上运行步骤3)得到Rpt文件，得到与每种情况相对应的Csv文件，以便平衡地应力。
5) 地应力平衡。模型的地应力平衡结果满足要求后，程序自动调用提前编写的命令读取Csv文件。
6) 得到最终的Inp文件。在初始Inp文件中加入强度折减法的分析步，得到最终的Inp文件。
7) Abaqus强度折减法运算。调用Abaqus求解器得到最终包含边坡变形、应力和场变量等信息的Job文件。
8) 计算边坡的失效概率。根据Pf=Nfs < 1/N(Pf表示边坡的失效概率; Nfs < 1表示安全系数小于1的数量; N表示总的计算次数)输出边坡的失效概率。
为验证编写的自动化算法程序的精度，采用经典边坡算例。边坡尺寸如图 3所示，坡高比为1:2。黏聚力均值为15 kPa，标准差为4 kPa，水平相关距离为38 m，竖向相关距离为3.8 m。在随机场理论中，相关距离是指土体中任意两点性质不相关的最小距离，是土体的天然特性。土体天然密度ρ为2 000 kg/m3，变形模量为10 MPa，泊松比v=0.3[17]。为简化计算，只考虑黏聚力生成的随机场，内摩擦角为0°。
边坡采用平面应变单元CPE4，共划分910个单元，971个单元节点，土体失效模式采用Mohr-Coulomb屈服准则。边界条件为约束边界的侧向位移及底部的水平及竖向位移。Der Kiureghian等[18]和Huang等[19]指出单元尺寸与相关距离之比应小于0.25。单元水平长度为2 m，高度为0.5 m，其中，单元水平长度/水平相关距离=2/38=0.05 < 0.25，单元高度/竖向相关距离=0.5/3.8=0.13 < 0.25，单元尺寸符合要求。
地应力平衡是岩土工程数值模拟过程中的重要内容，根据一般岩土工程对地应力平衡的要求，土体变形小于10-4 m即可满足工程实际要求[20]，非均匀随机场下自动化程序计算的边坡地应力平衡结果如图 4(a)所示，土体变形最大值所在的量级为10-5 m，满足要求。在非均匀随机场下边坡的失效变形模式如图 4(b)所示，为圆弧形面，符合规律。
图 5为黏聚力在非平稳随机场下的变化规律及均值线性趋势图。由于图 5可知，本程序计算的结果与Jiang等[17]的非平稳随机场下土体的黏聚力值均在其均值线性趋势线的右边。这是因为将平稳随机场转化为非平稳随机场，并没有将土体强度参数随深度增加的趋势分量与波动分量分开，波动分量并不明显。而假定的土体强度参数分布为对数正态分布，由对数正态函数的频率分布图可知，均值右侧的随机变量远多于左侧数据，因此，由随机场得到的随机变量大多浮动在线性趋势的右侧。
根据自动化程序得到10 000组安全系数的散点图，如图 6所示，安全系数大多分布在1~3之间。经统计得出，安全系数的均值为1.967，标准差为0.479，最小值为0.716，最大值为3.86。根据安全系数的分布直方图以及拟合的正态分布，可知安全系数服从正态分布。
图 7为不同模拟次数下对应的边坡失效概率，可知，当模拟次数在1 000~10 000之间时，由自动化程序得到的该边坡失效概率曲线趋于稳定，此时，所对应的边坡的失效概率为7.2‰。Jiang等[17]采用同样土体参数计算得到边坡失效概率为5.28‰，误差来源主要为失效概率计算方法的偏差，Jiang等[17]采用子集模拟法计算边坡的失效概率，而本文使用蒙特卡洛模拟法。子集模拟是一种求解失效概率的近似方法，所得到的结果是近似结果，而蒙特卡洛模拟是检验其他方法的依据，并且两者结果仅相差1.92‰，可以认为本文的计算结果准确。
基于Abaqus开放接口，使用Python语言进行二次开发，编写随机有限元脚本文件，用以计算边坡的可靠度。当边坡的几何形状确定后，只需运行几个特定的脚本文件，便可利用该程序求解基于随机场理论的边坡可靠度，使用方便。主要结论如下：
1) 该程序能够使用随机场理论自动计算边坡的失效概率。
2) 当土体黏聚力服从对数正态分布时，边坡的安全系数分布比较集中，服从正态分布。使用随机场理论计算边坡的稳定性时，边坡失效时的滑动面为圆弧面，符合规律。
3) 计算边坡的失效概率时，蒙特卡洛模拟次数较大时，计算得到的失效概率逐渐趋于稳定，在计算未知边坡的失效概率时，不必过多设置模拟次数，以免耗时过长，可以近似认为边坡失效概率曲线稳定时，对应的失效概率为边坡的实际失效概率。
2019, Vol. 41 Issue (5): 67-72
PDF
2019, Vol. 41 Issue (5): 67-72
PDF