CN110110371A - 基于极限分析下限定理的三维边坡安全系数迭代求解方法 - Google Patents

Abstract

本发明公开了一种基于极限分析下限定理的三维边坡安全系数迭代求解方法，该方法以块体界面作用力的大小、方向和作用点为主要变量，以满足极限分析下限定理许可静力场要求的静力平衡和屈服准则为约束条件，将边坡稳定分析问题转化为求解安全系数最大值的非线性优化问题；同时引入过渡变量临界加速度系数η，降低由于强度折减系数F引起的高度非线性；采用二分法迭代求得当η＝0时对应的折减系数即为所求安全系数。本发明为实际复杂边坡工程的三维稳定性分析提供了高效的手段。

Description

基于极限分析下限定理的三维边坡安全系数迭代求解方法

技术领域

本发明涉及一种评价岩土工程边坡稳定性的方法，具体地说，本发明涉及一种根据极限分析下限定理，通过迭代算法求解三维边坡安全系数的方法。本发明属于岩土工程边坡稳定性极限分析领域。

背景技术

目前在边坡稳定性分析领域，因缺少成熟有效的三维分析方法，通常都将三维边坡稳定性分析简化为平面问题来考虑，而忽视了实际工程的空间效应。这样的做法在多数情况下是可以接受的，但对于某些规模巨大、空间效应十分显著的三维边坡，如仍采用二维分析的方法显然是不合适的。

国内外不少学者曾尝试通过引入更多假定将边坡二维极限平衡法的思想推广应用到三维边坡稳定性分析中，这些假定进一步降低了极限平衡法的理论基础。为克服极限平衡法的不足，近年来已有不少学者探索将塑性力学极限分析理论与有限元结合，提出了基于优化算法的三维边坡稳定性分析方法，但是由于优化方程的高度非线性，使得计算三维边坡安全系数变得十分困难。

因此，进一步深入研究基于优化算法的三维边坡稳定性分析方法，降低约束方程非线性，提高基于优化算法的三维边坡极限分析方法的适用范围和精度是非常必要的。

发明内容

鉴于上述原因，从降低优化方程非线性、提高基于极限分析下限定理的三维边坡稳定性计算方法计算精度和适用范围出发，本发明的目的是提供一种基于极限分析下限定理的三维边坡安全系数迭代求解方法。

为实现上述目的，本发明采用以下技术方案：一种基于极限分析下限定理的三维边坡安全系数迭代求解方法，它包括如下步骤：

S1：根据三维边坡块分模型，引入过渡变量临界加速度系数η，以块体界面作用力的大小、方向和作用点为变量体系；

S2：建立满足极限分析下限定理许可静力场要求的约束方程，该约束方程包括静力平衡方程，不违反屈服准则、岩土体材料不受拉和荷载作用点在界面内部的约束不等式；

S3：在给定折减系数条件下，以临界加速度最大值为目标函数，以满足极限分析下限定理许可静力场要求的约束方程为约束条件，形成边坡临界加速度优化模型，计算该折减系数对应的临界加速度；

S4：采用二分法迭代求得当η＝0时对应的折减系数即为所求安全系数。

本发明根据三维边坡块分模型，引入过渡变量临界加速度系数，以块体界面作用力的大小、方向和作用点为主要变量；以满足极限分析下限定理许可静力场要求的约束方程为约束条件，在给定折减系数条件下，以临界加速度最大值为目标函数，形成三维边坡临界加速度优化模型，计算该折减系数对应的临界加速度，再通过二分法迭代求得当临界加速度系数为零时对应的折减系数，即为所有安全系数。本发明将三维边坡稳定性分析问题转化为求解三维边坡安全系数最大值的非线性优化问题；同时引入过渡变量临界加速度系数，降低由于强度折减系数引起的高度非线性；采用二分法迭代求得当临界加速度系数等于零时对应的折减系数即为所求安全系数。

附图说明

图1是本发明基于极限分析下限定理的三维边坡安全系数迭代求解方法流程图；

图2是本发明实例虚拟的三维边坡模型；

图3是临界加速度系数η定义示意图；

图4是本发明根据极限分析下限定理优化求解三维边坡临界加速度流程图；

图5是本发明根据极限分析下限定理通过二分法迭代优化求解三维边坡安全系数的流程图；

图6是利用本发明评价图的实施例图。

具体实施方式

下面结合附图进一步阐明本发明，应理解这些附图及实例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

为降低目前三维边坡稳定性分析中优化方程的非线性、提高基于极限分析下限定理的三维边坡稳定性计算方法的计算精度和适用范围，如图1所示，本发明提出一种基于极限分析下限定理的三维边坡安全系数迭代求解方法，它包括如下步骤：

S1：根据三维边坡块分模型，引入过渡变量临界加速度系数η，以块体界面作用力的大小、方向和作用点为变量体系；

S2：建立满足极限分析下限定理许可静力场要求的约束方程，该约束方程包括静力平衡方程，不违反屈服准则、岩土体材料不受拉和荷载作用点在界面内部的约束不等式；

S3：在给定折减系数F条件下，以临界加速度最大值为目标函数，以满足极限分析下限定理许可静力场要求的约束方程为约束条件，形成边坡临界加速的优化模型，计算该折减系数对应的临界加速度；

S4：采用二分法迭代求得当η＝0时对应的折减系数即为所求安全系数。

本发明评价三维边坡稳定性的原理是：根据三维边坡块分模型，引入过渡变量临界加速度系数，以块体界面作用力的大小、方向和作用点为主要变量；以满足极限分析下限定理许可静力场要求的约束方程为约束条件，在给定折减系数条件下，以临界加速度最大值为目标函数，形成三维边坡临界加速度优化模型，计算该折减系数对应的临界加速度，再通过二分法迭代求得当临界加速度系数为零时对应的折减系数，即为所有安全系数。本发明将三维边坡稳定性分析问题转化为求解三维边坡安全系数最大值的非线性优化问题；同时引入过渡变量临界加速度系数η，降低由于强度折减系数F引起的高度非线性；采用二分法迭代求得当η＝0时对应的折减系数即为所求安全系数。

本发明步骤S1：根据三维边坡块分模型，引入过渡变量临界加速度系数η，以块体界面作用力的大小、方向和作用点建立变量体系，具体方法如下：

S1.1：如图2所示，根据三维边坡的实际情况，确定该三维边坡的多边形界面，将三维边坡分割成若干个块体，通过多边形界面集合和块体集合定义该三维边坡模型V。

该三维边坡模型V通过构成该三维边坡的多边形界面PO_j集合定义为：

V＝PO₁PO₂…PO_j…PO_p-1PO_p(1)

式中，j为边界面的编号，PO_j为编号为j的多边形界面，p为三维边坡V中多边形界面总数，PO_j的地层岩性、风化程度等特征将通过界面参数(黏聚力c_j和内摩擦角)反映。

该三维边坡模型V通过组成该三维边坡的块体B_k集合定义为：

V＝B₁B₂…B_k…B_q-1B_q(2)

式中，k为边界面的编号，B_k为编号为k的块体，q为三维边坡V中块体总数，B_k的地层岩性、风化程度等特征将通过块体参数(重度d_k)反映。

空间多边形界面PO_j由一系列按顺序排列的顶点来描述，可用下式表示为：

式中i表示多边形界面PO_j中顶点的编号，为编号为j的界面上第i个顶点，h1为多边形PO_j的顶点总数；为表示方便，将多边形界面PO_j最后一个点与第一点设为重合以形成闭合环路，即

多边形界面PO_j的边界为：

式中为多边形PO_j的边界线段。

多边形PO_j的法线向量np_j可由式(5)求得：

当i＝1时

多边形PO_j的边界线段所对应的内法线向量为：

需要指出的是，通过上述公式计算得到的多边形边界线段内法线向量永远指向多边形内部，与顶点旋转方向无关。

根据上述法线向量定义，三维边坡空间任意一点x位于空间多边形界面PO_j内的充要条件为：

式中，ap_j为多边形界面PO_j内任一定点。

S1.2：根据强度储备的概念，引入描述边坡强度的折减系数F，通过该变量描述边坡安全储备，如下式所示；

c_e＝c/F (8)

式中c、为边坡材料实际凝聚力和内摩擦系数，c_e、f_e为边坡材料进行折减后的凝聚力和内摩擦系数，c_e、f_e对应于不同的折减系数F计算得到。

S1.3：为降低非线性求解难度，引入过渡变量临界加速度系数η，如图3所示；通过引入该变量，在主滑方向l(l_x，l_y，0)上对每个块体B_k施加一个虚拟的水平惯性力η_k，使得三维边坡达到极限平衡状态：

η_k＝ηG_k(l_x，l_y，0) (10)

式中G_k为块体B_k的重力；通常，η取值范围在(-1，1)之间。

S1.4为描述三维边坡界面力的方向特性，在每个结构面上建立局部坐标系。

在三维边坡模型V中的任一界面PO_j(除边坡轮廓面外)的作用力包括法向力和位于该界面的切向力(可用两个分量表示)，其量值分别为N_j、Tm_j和Tl_j。

为描述界面力的方向特性，在每个空间多边形界面PO_j上建立临时局部坐标系，以该面的法线向量np_j为N‘轴正向；以空间多边形任一边界线段的内法线向量为Tm‘轴正向，i可取为多边形PO_j的任一边界线段，通常取i＝1；Tl‘轴正向则可根据N‘轴、Tm‘轴通过右手螺旋法则确定对于包含该界面的块体，还需进一步根据该界面在对应块体内的内法线向量的关系，最终确定该界面上局部坐标和轴的正向：

并有

式中i表示块体B_k中边界面的局部编号(区别与全局编号，块体B_k内从1～n的连续编号)，为在块体B_k中编号为i的多边形界面，为界面上作用力的三个分量；j为边界面的编号，PO_j为编号为j的多边形界面；和PO_j表示同一个界面，但局部坐标系可能不同。

通过这样的处理，每条界面上的界面力只需要一组变量N_j、Tl_j和Tm_j，就可以准确描述共用界面上作用力与反作用力的关系。

S1.5在三维问题中，多边形界面PO_j上界面力的作用点r_j可用其可空间坐标(r_xj，r_yj，r_zj)表示。为保证作用点始终位于空间多边形界面PO_j内，必须满足如下条件：

条件1，作用点在空间多边形PO_j所在平面上：

(r_j-ap_j)·np_j＝0(12)

条件2：作用点位于界面凸多边形区域内：

本发明根据三维边坡块分模型，引入过渡变量临界加速度系数，建立以块体界面作用力的大小、方向和作用点为主要变量的变量体系后，如图4所示，执行步骤S2：建立满足极限分析下限定理许可静力场要求的约束方程，该约束方程包括静力平衡方程，不违反屈服准则、岩土体材料不受拉和荷载作用点在界面内部的约束不等式；具体方法如下：

S2.1：根据极限分析下限定理，建立目标函数；

在给定折减系数F的条件下，求临界加速度最大值：

max(η)(14)

S2.2：在三维边坡块分系统中需要保证每一个块体力系平衡，力系平衡的充分必要条件是力系中各力在三个坐标轴上投影的代数和为零，且各力对三个坐标轴的力矩之和也分别为零。

具体到每一个块体B_k而言，力和力矩的平衡方程可以表示为：

式中，i表示块体B_k中边界面或外部荷载的局部编号(区别与全局编号，块体B_k内从1～n的连续编号)；为在块体B_k中编号为i的多边形界面， 为界面上作用力的三个分量，h2为块体B_k内界面总数；为块体B_k中编号为i的外部荷载，为已知量，h3为作用于B_k的外部荷载总数；l_p为投影轴的方向向量；M_k为旋转轴原点；为多边形边界上的作用点到力矩点M_k的力臂矢量；为外部荷载作用点到力矩点M_k的力臂矢量；RG_k＝RG_k-M_k为块体B_k重心RG_k到力矩点M_k的力臂矢量。并有：

S2.3：不违反屈服准则约束条件：每一个界面PO_j上的作用力应不违反摩尔库仑屈服准则，如下式所示：

c_j、分别为多边形界面PO_j的凝聚力和内摩擦角；A_j为界面PO_j的面积。除边坡轮廓面外。

上述不等式方程(15)、(16)、(17)与各界面一一对应。式(17)的几何意义为，切向力T_j(Tm_j，Tl_j)在多边形界面PO_j上位于以界面力作用点r_j为圆心、半径为的圆内。

为降低屈服准则约束条件的非线性，将上述圆用其内接等边多边形代替，则可转化为：

式中，h4为内接等边多边形边的条数，θⁱ为mlⁱ与Tm轴负方向的夹角。

S2.4：为反映岩土体不受拉的特性，采用了式(19)所示的约束条件：

N_j≥0(19)

本发明建立了满足极限分析下限定理许可静力场要求的约束方程后，执行步骤S3：在给定折减系数F条件下，以临界加速度最大值max(η)为目标函数，以满足下限定理许可静力场要求的约束方程为约束条件，形成边坡临界加速度的优化模型，通过输出常用的优化求解算法计算文件，或通过编程软件，计算求出对应于不同折减系数F的临界加速度η，具体方法如下：

所述临界加速度优化模型为：

N_j≥0

根据优化模型生成优化计算文件计算临界加速度。

本发明步骤S4：根据折减系数F与临界加速度η之间的对应关系，采用二分法迭代求得当η＝0时对应的折减系数即为所求安全系数，具体方法，如图5所示：

S4.1：确定初始搜索范围，根据一般工程边坡安全系数初定搜索范围为：F^l＝1.0、F^u＝2.0；

S4.2：对初始范围合理性进行验证，并确定最终的搜索范围；

计算F^l对应的临界加速度η(F^l)，如η(F^l)＞0则令F^l＝F^l/2直到满足η(F^l)＜0时为最终使用的F^l；同理，计算F^u对应的临界加速度η(F^u)，如η(F^u)＜0则令F^u＝2F²直到满足η(F^u)＞0时为最终使用的F^u。

S4.3：确定新的搜索范围，令F^m＝(F^l+F^u)/2，得到η(F^l)、η(F^m)、η(F^u)，以相邻两个η乘积为负的原则，确定新的F^l，F^m。

S4.4：重复S4.3，直到F^m＜1e-3或η(F^m)＜1e-5，此时对应的F^m即为安全系数。

对初始范围合理性进行验证，并确定最终的搜索范围。计算F^l对应的临界加速度η(F^l)，如η(F^l)＞0则令F^l＝F^l/2直到满足η(F^l)＜0时为最终使用的F^l；同理计算F^u对应的临界加速度η(F^u)，如η(F^u)＜0则令F^u＝2F²直到满足η(F^u)＞0时为最终使用的F^u。

图6是利用本发明公开的一种基于极限分析下限定理的三维边坡安全系数迭代求解方法的实施例图。

为便于算例对比分析，对A-A、B-B、C-C和D-D四个剖面采用EMU法进行了二维边坡稳定性分析，计算结果表明剖面A-A、B-B计算得到的安全系数分别为1.34、1.27，说明边坡具有一定安全裕度，而剖面C-C、D-D计算得到的安全系数则分别为0.80、0.53，说明边坡已经失稳破坏。显然，根据这样的结果无法对边坡的整体稳定性做出判断，需要考虑实际边坡的三维效应，通过三维分析才能对边坡的稳定情况给出切合实际的评价。

本发明和EMU上限解法计算得到的三维稳定分析结果表明：本发明计算得到的三维边坡稳定安全系数为1.054，表明边坡基本处于临界状态，与边坡实际现状完全吻合，进一步验证了本发明的合理性和可行性，而EMU给出的计算结果为1.102，与边坡实际状况相比明显偏高，这也从一个侧面说明了上限解会在很多情况下给出偏于危险的结果，在实际应用中应引起高度重视。

本发明通过引入过渡变量临界加速度系数η使得式(17)由原来4次方程变为2次方程，式(18)由二次方程变为线性方程，极大降低了约束方程的非线性，使整个优化计算过程更简化。另外，本发明结合边坡安全系数与过渡变量之间的相关特性，提出基于二分法的迭代求解方法，求得η＝0时对应的折减系数F即为安全系数，提高了计算求解效率。本发明公开的基于极限分析下限定理的三维边坡安全系数迭代求解方法，为实际边坡工程的三维稳定性分析提供了新的手段。

Claims (5)

1.一种基于极限分析下限定理的三维边坡安全系数迭代求解方法，其特征在于：它包括如下步骤：

S1：根据三维边坡块分模型，引入过渡变量临界加速度系数η，以块体界面作用力的大小、方向和作用点为变量体系；

S2：建立满足极限分析下限定理许可静力场要求的约束方程，该约束方程包括静力平衡方程，不违反屈服准则、岩土体材料不受拉和荷载作用点在界面内部的约束不等式；

S3：在给定折减系数F条件下，以临界加速度最大值为目标函数，以满足极限分析下限定理许可静力场要求的约束方程为约束条件，形成三维边坡临界加速度优化模型，计算该折减系数对应的临界加速度；

S4：采用二分法迭代求得当η＝0时对应的折减系数即为所求安全系数。

2.根据权利要求1所述的基于极限分析下限定理的三维边坡安全系数迭代求解方法，其特征在于：所述步骤S1根据三维边坡块分模型，引入过渡变量临界加速度系数η，以块体界面作用力的大小、方向和作用点建立变量体系，具体方法如下：

S1.1：根据三维边坡的实际情况，确定该三维边坡的多边形界面，将三维边坡分割成若干个块体，通过多边形界面集合和块体集合定义该三维边坡模型V，

该三维边坡模型V通过构成该三维边坡的多边形界面PO_j集合定义为：

V＝PO₁PO₂…PO_j…PO_p-1PO_p (1)

式中，j为边界面的编号，PO_j为编号为j的多边形界面，p为三维边坡V中多边形界面总数，PO_j的地层岩性、风化程度、结构面分布特征通过界面参数反映；

该三维边坡模型V通过组成该三维边坡的块体B_k集合定义为：

V＝B₁B₂…B_k…B_q-1B_q (2)

式中，k为边界面的编号，B_k为编号为k的块体，q为三维边坡V中块体总数，B_k的地层岩性、风化程度等特征将通过块体重度参数反映；

空间多边形界面PO_j由一系列按顺序排列的顶点来描述，可用下式表示为：

式中i表示多边形界面PO_j中顶点的编号，为编号为j的界面上第i个顶点，h1为多边形PO_j的顶点总数；为表示方便，将多边形界面PO_j最后一个点与第一点设为重合以形成闭合环路，即

多边形界面PO_j的边界为：

式中为多边形PO_j的边界线段；

多边形PO_j的法线向量np_j可由式(5)求得：

当i＝1时

多边形PO_j的边界线段所对应的内法线向量为：

需要指出的是，通过上述公式计算得到的多边形边界线段内法线向量永远指向多边形内部，与顶点旋转方向无关；

根据上述法线向量定义，三维边坡空间任意一点x位于空间多边形界面PO_j内的充要条件为：

式中，ap_j为多边形界面PO_j内任一定点；

S1.2：根据强度储备的概念，引入描述边坡强度的折减系数F，通过该变量描述边坡安全储备，如下式所示；

c_e＝c/F (8)

式中c、为边坡材料实际凝聚力和内摩擦系数，c_e、f_e为边坡材料进行折减后的凝聚力和内摩擦系数，c_e、f_e对应于不同的折减系数F计算得到；

S1.3：为降低非线性求解难度，引入过渡变量临界加速度系数η，通过引入该过渡变量，在主滑方向l(L_x，L_y，0)上对每个块体B_k施加一个虚拟的水平惯性力η_k，使得三维边坡达到极限平衡状态：

η_k＝ηG_k(l_x，l_y，0) (10)

式中G_k为块体B_k的重力；通常，η取值范围在(-1，1)之间；

S1.4为描述三维边坡界面力的方向特性，在每个结构面上建立局部坐标系；

在三维边坡模型V中的任一界面PO_j(除边坡轮廓面外)的作用力包括法向力和位于该界面的切向力(可用两个分量表示)，其量值分别为N_j、Tm_j和Tl_j；

为描述界面力的方向特性，在每个空间多边形界面PO_j上建立临时局部坐标系，以该面的法线向量np_j为N‘轴正向；以空间多边形任一边界线段的内法线向量为Tm‘轴正向，i可取为多边形PO_j的任一边界线段，通常取i＝1；Tl‘轴正向则可根据N‘轴、Tm‘轴通过右手螺旋法则确定对于包含该界面的块体，还需进一步根据该界面在对应块体内的内法线向量与np_j的关系，最终确定该界面上局部坐标和轴的正向：

并有

式中i表示块体B_k中边界面的局部编号(区别与全局编号，块体B_k内从1～n的连续编号)，为在块体B_k中编号为i的多边形界面，为界面上作用力的三个分量；j为边界面的编号，PO_j为编号为j的多边形界面；和PO_j表示同一个界面，但局部坐标系可能不同；

通过这样的处理，每条界面上的界面力只需要一组变量N_j、Tl_j和Tm_j，就可以准确描述共用界面上作用力与反作用力的关系；

S1.5在三维问题中，多边形界面PO_j上界面力的作用点r_j可用其可空间坐标(r_xj，r_yj，r_zj)表示；为保证作用点始终位于空间多边形界面PO_j内，必须满足如下条件：

条件1，作用点在空间多边形PO_j所在平面上：

(r_j-ap_j)·np_j＝0 (12)

条件2：作用点位于界面凸多边形区域内：

3.根据权利要求2所述的基于极限分析下限定理的三维边坡安全系数迭代求解方法，其特征在于：所述步骤S2：建立满足极限分析下限定理许可静力场要求的约束方程，该约束方程包括静力平衡方程，不违反屈服准则、岩土体材料不受拉和荷载作用点在界面内部的约束不等式；具体方法如下：

S2.1：根据极限分析下限定理，建立目标函数；

在给定折减系数F的条件下，求临界加速度最大值：

max(η) (14)

S2.2：在三维边坡块分系统中需要保证每一个块体力系平衡，力系平衡的充分必要条件是力系中各力在三个坐标轴上投影的代数和为零，且各力对三个坐标轴的力矩之和也分别为零；

具体到每一个块体B_k而言，力和力矩的平衡方程可以表示为：

式中，i表示块体B_k中边界面或外部荷载的局部编号(区别与全局编号，块体B_k内从1～n的连续编号)；为在块体B_k中编号为i的多边形界面， 为界面上作用力的三个分量，h2为块体B_k内界面总数；为块体B_k中编号为i的外部荷载，为已知量，h3为作用于B_k的外部荷载总数；lp为投影轴的方向向量；M_k为旋转轴原点；为多边形边界上的作用点到力矩点M_k的力臂矢量；为外部荷载作用点到力矩点M_k的力臂矢量；RG_k＝RG_k-M_k为块体B_k重心RG_k到力矩点M_k的力臂矢量，并有：

S2.3：不违反屈服准则约束条件：每一个界面PO_j上的作用力应不违反摩尔库仑屈服准则，如下式所示：

c_j、分别为多边形界面PO_j的凝聚力和内摩擦角；A_j为界面PO_j的面积；

式(17)的几何意义为，切向力T_j(Tm_j，Tl_j)在多边形界面PO_j上位于以界面力作用点r_j为圆心、半径为的圆内；

为降低屈服准则约束条件的非线性，将上述圆用其内接等边多边形代替，则可转化为：

式中，h4为内接等边多边形边的条数，θⁱ为mlⁱ与Tm轴负方向的夹角；

S2.4：为反映岩土体不受拉的特性，采用了式(19)所示的约束条件：

N_j≥0 (19)。

4.根据权利要求3所述的基于极限分析下限定理的三维边坡安全系数迭代求解方法，其特征在于：所述步骤S3在给定折减系数F条件下，以临界加速度最大值max(η)为目标函数，以满足下限定理许可静力场要求的约束方程为约束条件，形成边坡临界加速度的优化模型，通过输出常用的优化求解算法计算文件，或通过编程软件，计算求出对应于不同折减系数F的临界加速度η；

所述三维边坡临界加速度优化模型为：

N_j≥0。

5.根据权利要求4所述的基于极限分析下限定理的三维边坡安全系数迭代求解方法，其特征在于：所述步骤S4根据折减系数F与临界加速度η之间的对应关系，采用二分法迭代求得当η＝0时对应的折减系数即为所求安全系数，具体方法为：

S4.1：确定初始搜索范围，根据一般工程边坡安全系数初定搜索范围为：F^l＝1.0、F^u＝2.0；

S4.2：对初始范围合理性进行验证，并确定最终的搜索范围；

计算F^l对应的临界加速度η(F^l)，如η(F^l)＞0则令F^l＝F^l/2直到满足η(F^l)＜0时为最终使用的F^l；同理，计算F^u对应的临界加速度η(F^u)，如η(F^u)＜0则令F^u＝2F²直到满足η(F^u)＞0时为最终使用的F^u；

S4.3：确定新的搜索范围，令F^m＝(F^l+F^u)/2，得到η(F^l)、η(F^m)、η(F^u)，以相邻两个η乘积为负的原则，确定新的F^l，F^m；

S4.4：重复S4.3，直到F^m＜1e-3或η(F^m)＜1e-5，此时对应的F^m即为安全系数。