一种基于精确有限元法的土质边坡稳定性分析方法

摘要

本发明公开了一种基于精确有限元法的土质边坡稳定性分析方法，建立二维边坡模型并划分网格；采用修正后的双曲圆弧化摩尔‑库仑屈服本构来判断土体是否发生屈服，若发生屈服将采用一种精确的显示欧拉算法来进行应力重分布，并通过应力回拉始终保证应力点在屈服面上或面内；对非线性方程组的求解，采用牛顿‑拉夫逊迭代进行求解；引入强度折减法，并采用足够精度要求的破坏判据，通过迭代搜索求得边坡的安全系数，对结果进行后处理也可获取边坡的位移和塑性区等数据。本发明针对边坡的稳定性分析合理且准确，并能给出一个高精度的边坡安全系数，同时根据计算给出的边坡位移和塑性区等数据，可以为边坡工程防治措施的实施提供合理的依据或建议。

说明书

技术领域

本发明属于计算机仿真计算领域，涉及计算机有限元建模及计算分析，具体涉及一种基于精确有限元法的土质边坡稳定性分析方法。

背景技术

滑坡灾害往往会对河流、农田、交通沿线的公路、铁路及其他建筑物造成破坏，严重威胁了人类生命及财产安全，如何计算出更为精确的边坡安全系数，如何获取边坡地应力应变、位移或变形等数据，这对边坡的工程防治至关重要。传统的极限平衡法和极限分析法通常需要事先知道或者假定滑面，并假设滑体是刚体，诸多的假设限制了这些方法的使用，难以分析复杂边坡的稳定性，如考虑空间变异性的非均质边坡，同时不能获取边坡的变形和应力应变状态。

然而随着计算机技术的发展，克服了大型数据计算的难题，有限元法发展迅速，有限元法没有极限平衡法和极限分析法诸多假设，适用于各类复杂边坡的稳定性分析，深受广大学者的欢迎。

鉴于此，本文基于有限元法和强度折减法，同时使用了一种高精度的显示欧拉算法进行应力重分布，并对每一次的重分布应力再进行回拉，保证应力点始终保持在屈服面上，为了避免传统摩尔-库仑屈服本构的奇异性，对其进行双曲圆弧化的修正。最终获取更为精确的边坡安全系数和边坡变形和应力应变状态。

发明内容

本发明基于一种精确有限元法，发明了一种土质边坡稳定性分析方法，适用于各类复杂边坡，也可考虑土体的空间变异性，能精确地求解出边坡的安全系数，同时通过程序后处理可以获取边坡各位置处的应力应变、位移或变形状态，可提供边坡潜在的塑性区或滑面位置。为滑坡的工程防治提供可靠的依据。

本文所采用的技术方案是：

一种基于精确有限元法的土质边坡稳定性分析方法，包括以下步骤：

步骤一：获取边坡几何参数和边界条件，建立边坡二维有限元模型并求解，获得单元的初始位移、应变和应力；

步骤二：将初始高斯点的应力当作预测应力，考虑修正后的双曲圆弧化摩尔-库仑屈服本构来判断该高斯点是否发生屈服，如果发生屈服，需要进行应力重分布，没有发生屈服直接进行下一高斯点是否屈服的判断，直至所有单元所有高斯点判断完毕；

步骤三：由于发生应力重分布，弹性状态下的线性方程组求解问题变为弹塑性状态下的非线性方程组求解问题，采用牛顿-拉夫逊迭代进行求解，根据重分布后的应力和塑性流动法则求得弹塑性应力应变矩阵，并构建出新的整体刚度方程；

步骤四：根据塑性流动法则计算得到高斯点的塑性应变，同时随着迭代，单元结点位移、应变、应力和塑性应变在不断累积，直至达到所设置的最大迭代次数，视为边坡已发生破坏，或者满足收敛判据；

步骤五：引入强度折减法进行计算，通过二分法不断变更折减系数重复步骤一至步骤四，直至搜索至边坡极限破坏状态并满足精度要求，此时的折减系数即为边坡的安全系数。

进一步地，根据步骤一中计算得到的位移，进一步求得应变增量Δε

Δε

σ

Δu

采用修正的双曲圆弧化摩尔-库仑屈服本构进行高斯点屈服的判断：

s

s

J

J

F为屈服函数，σ

进一步地，在步骤二中，若是第一次发生屈服需要计算弹性到塑性过渡阶段的比例因子，求得用于后续应力重分布计算的塑性应力增量和重分布前的应力：

其中比例因子计算式为：

z为所求比例因子，Δσ

若为第二次进入塑性状态，比例因子z取零，弹性应变增量全作为塑性应变增量进行应力重分布。

进一步地，在步骤二中，采用Abbo提出的显示欧拉算法进行应力重分布，针对每一发生屈服的高斯点逐一进行应力重分布，分为以下五个步骤：

(1)重分布前应力为：

σ

(2)设置伪时间，初始伪时间T＝0，伪时间增量步ΔT＝1，直至T＝1应力重分布结束，每一伪时间步下的应力增量为：

Δσ

采用理想弹塑性模型，A＝0，Δλ

为了获得准确的应力增量，采用了更高阶显示欧拉方法：

当前子步的误差为：

(3)如果R

ΔT＝max{qΔT，ΔT

其中容许误差STOL设置为0.05，最小相对误差Eps设置为10

T＝T+ΔT

(4)由步骤(3)计算得到满足要求的子步下的应力，进一步判断该应力点是否在屈服面上，如果在屈服面外需要进行应力回拉，计算完成后进行下一子步的计算：

如果应力回拉失败，令q＝min{q，1}，设置下一迭代的新伪时间增量步：

ΔT＝min{ΔT，ΔT

ΔT＝qΔT

其中设置最小的伪时间增量步ΔT

(5)重复步骤(2)-(4)直至满足T＝1时，此时的应力即是重分布后的应力。

进一步地，在步骤(4)的所述重分布后的应力还需要进行是否屈服的判断，如果所述重分布后的应力仍处于屈服面外，还需再对应力进行修正回拉至屈服面上，计算步骤如下：

(1)对于某一高斯点进行重分布后但未修正的应力σ

(2)对初始未修正应力进行修正：

σ＝σ

a

(3)由于在应力修正过程中可能会发生修正后应力比初始应力更远离屈服面，因此需要对修正后应力作进一步判断，即：如果|F(σ)|＞|F(σ

σ＝σ

(4)重复步骤(2)-(3)计算10次，如果此时修正后的应力满足|F(σ)|≤FTOL，则退出循环；

(5)若循环次数小于10次，认为应力修正成功，更新应力σ

进一步地，在步骤三中，应力应变曲线不再呈线性关系，通过塑性流动法则求得弹塑性应力应变矩阵D

重新构建新的刚度矩阵，再根据牛顿-拉夫逊迭代求解该大型非线性方程组：

Δu

(K

P为第n步迭代迭代过程中的内力：

P＝∫∫∫B

进一步地，在步骤四中，采用的牛顿-拉夫逊迭代法，设置最大迭代次数为25次，同时设置位移和不平衡力的迭代收敛判据：

||Δu

其中Tol取0.00001，在迭代过程中，位移、应力、应变和塑性应变都在不断累积，其中塑性应变根据流动法则，塑性应变计算式为：

Q为塑性势函数。

进一步地，在步骤五中，为了考虑边坡安全储备问题，引入强度折减法：

其中FS为折减系数，在边坡安全系数的求解过程中，用于二分法搜索，不断重复步骤一至步骤四的计算，直至FS收敛于边坡刚好处于极限状态，此折减系数FS即为边坡安全系数。

进一步地，根据步骤四获取的每个单元高斯点的应力和塑性应变，可得到边坡各位置处的应力或应变，先通过外推的方法获取每个单元节点的应力或塑性应变：

σ

为含

为了更直观体现边坡塑性区的位置，在平面应变问题中，可将每个节点求得的四个塑性应变分量通过下式计算得到等效的塑性应变；

另一方面，本申请还保护一种电子设备，包括存储器、处理器及存储在所述存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现权利要求1至9任意一项所述的分析方法。

与现有技术相比，本发明具有以下益处：

1、本申请采用了一种高阶的应力重分布方法：显示欧拉算法，比传统有限元法采用的粘弹塑性方法的计算精度更高，可以获取更高准确的边坡安全系数和边坡位移、应力和应变等数据。

2、相比于其它方法，如极限平衡法和极限分析法等边坡稳定性分析方法，本方法不需要事先假设滑面，可以计算各类复杂的非均质边坡的安全系数，同时可以考虑土体空间变异性及其他性质的影响。

3、有些数值软件计算只考虑摩尔-库仑屈服本构，或者由于软件内嵌程序不可见，难以根据需求对本构进行直接修改，本方法可以针对不同需求，对屈服本构进行修改而进一步求解，例如本方法采用了修正的双曲圆弧化摩尔-库仑屈服本构。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作一简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本申请一种基于精确有限元法的土质边坡稳定性分析方法的计算流程图。

图2为本申请一种基于精确有限元法的土质边坡稳定性分析方法的边坡计算模型参考图。

图3为本申请一种基于精确有限元法的土质边坡稳定性分析方法计算采用的四节点矩形单元。

图4为本申请一种基于精确有限元法的土质边坡稳定性分析方法修正后的π平面上的双曲圆弧化摩尔-库仑屈服面，可以消除传统摩尔-库仑屈服面在尖端和边缘处的奇异性。

图5为本申请一种基于精确有限元法的土质边坡稳定性分析方法高斯点应力从弹性阶段过渡到塑性阶段，用于求解比例因子。

图6为本申请一种基于精确有限元法的土质边坡稳定性分析方法显示欧拉算法下的应力重分布计算流程图。

图7为本申请一种基于精确有限元法的土质边坡稳定性分析方法应力重分布过程中，应力回拉至屈服面上的计算流程图。

图8为本申请一种基于精确有限元法的土质边坡稳定性分析方法用于计算非线性方程组的牛顿-拉夫逊迭代方法的计算流程图。

图9为本申请一种基于精确有限元法的土质边坡稳定性分析方法的计算机程序后处理得到的边坡变形图(放大处理)。

图10为本申请一种基于精确有限元法的土质边坡稳定性分析方法的计算机程序后处理得到的边坡塑性区云图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

下面结合图1-10，具体介绍本申请的实施例。

如图1所示为本发明的计算流程图，总结并介绍了计算流程，具体步骤如下：

步骤一：获取边坡的如下参数：边坡几何参数、容重参数、强度参数及其变形参数。

边坡几何参数有坡顶长度w1，坡面水平长度s1，坡底长度w2，边坡高度h1，底座厚度h2，坡顶处划分单元数目nx1，坡底处划分单元数目nx2，坡高划分单元数目ny1，底座厚度划分单元数目ny2；具体参数取值如图2所示。

土体参数包括土的重度γ，强度参数有粘聚力c，摩擦角

步骤二：建立二维边坡模型，对模型进行网格划分，获取各单元编号和节点坐标信息，对单元进行材料参数赋值，设置模型的边界条件，构建出完整的有限元边坡计算模型。

边坡模型的建立和网格划分可通过自编程序、MATLAB或ANYSYS软件生成，获取单元的编号顺序、节点的坐标和边界条件等数据，并为每个单元赋予材料参数，用于后续有限元计算。若要考虑土体的空间变异性，需要分别对每个单元的土体参数进行赋值。

步骤三：根据单元和节点信息，在局部坐标系ξ-η下，通过插值的方式构建单元位移场，引入形函数，采用高斯积分构建单元的刚度矩阵和等效节点力矩阵并对其进行集成，考虑模型的边界条件，形成初始的整体刚度方程。

构建单元位移场：

针对四节点四边形单元，如图3所示，位移场为：

N为形函数矩阵，q为节点位移矩阵。

局部坐标系下某单元n节点的形函数表达式为：

N

构建应变场为：

ε(x，y)＝B·q

B为几何矩阵，写成分块形式：

B＝[B

其中：

构建应力场为：

考虑平面应变问题，处于弹性状态时弹性矩阵D为：

构建单元刚度方程：

Kq＝f

其中单元刚度矩阵K和等效节点力矩阵f均采用高斯积分进行求解：

G(ξ，η)＝B

其中f

将单元刚度方程根据节点的编号集成为整体的刚度方程，在求解该线性方程组的过程中，为避免方程的奇异性，采取乘大数法的方法对边界条件进行处理，对现有刚度方程进行修改。采用乘大数法计算举例如下：现假设边界条件为：

u

β

现有刚度方程为：

采用乘大数法修正为：

步骤四：求解初始的整体刚度方程(大型线性方程组)，可获取单元的初始位移、应变和应力。

采用高斯消去法或高斯-赛德尔迭代法求解该线性方程组，也可采用MATLAB自带的求解方法。

步骤五：将初始高斯点的应力当作预测应力，考虑修正后的双曲圆弧化摩尔-库仑屈服本构来判断该高斯点是否发生屈服，如果发生屈服，需要进行应力重分布，没有发生屈服直接进行下一高斯点是否屈服的判断，直至所有单元所有高斯点计算完毕。

根据第四步计算得到的位移，进一步求得应变增量Δε

Δε

σ

采用修正的双曲圆弧化摩尔-库仑屈服本构，如图所示4，依次对每一高斯点进行是否屈服的判断：

屈服函数表达式为：

其中：

s

s

步骤六：若是第一次发生屈服需要计算弹性到塑性过渡阶段的比例因子，如图5所示。

其中比例因子计算式为

若为第二次进入塑性状态，比例因子z取零，弹性应变增量全作为塑性应变增量进行应力重分布。

步骤七：采用Abbo提出的显示欧拉算法进行应力重分布，计算流程如图6所示。

根据第六步求得的比例因子确定塑性应变增量，然后采用Abbo提出显示欧拉算法进行应力重分，针对每一发生屈服的高斯点都需要进行应力重分布。现针对某一发生屈服的高斯点应力进行分析，主要可分为以下五小步：

(1)重分布前应力为：

σ

(2)设置伪时间，初始T＝0，ΔT＝1，直至T＝1应力重分布结束。每一时间步下的应力增量为：

Δσ

采用理想弹塑性模型，A＝0，b为塑性势函数对应力的偏导数，a为屈服函数对应力的偏导数，考虑相关联流动法则塑性势函数等于屈服函数。

为了获得准确的应力增量，采用了更高阶的计算方法

当前子步的误差为：

(3)如果R

ΔT＝max{qΔT，ΔT

其中最小相对误差Eps设置为10

T＝T+ΔT

(4)由步骤七第(3)步计算得到满足要求的子步下的应力，进一步判断该应力点是否在屈服面上，如果在屈服面外根据步骤八进行应力回拉，计算完成后进行下一子步的计算：

如果应力回拉失败，令q＝min{q，1}，设置下一迭代的新伪时间增量步：

ΔT＝min{ΔT，ΔT

ΔT＝qΔT

其中设置最小的伪时间增量步ΔT

(5)重复步骤七第(2)步-第(4)步直至满足T＝1时，此时的应力即是所需要重分布应力。

步骤八：应力重分布过程中若高斯点应力仍在屈服面外则需要对应力进行修正，回拉至屈服面上，计算流程图如图7所示。

根据步骤七第(4)步中求得的重分布后的应力，还需要对其进行是否屈服的判断，如果该重分布后的应力仍处于屈服面外，还需再对应力进行修正回拉至屈服面上，计算分为以下5个小步：

(1)假设对于某一高斯点初始未修正应力σ

(2)对初始未修正应力进行修正：

σ＝σ

(3)由于在应力修正过程中可能会发生修正后应力比初始应力更远离屈服面，这时需要对修正后应力作进一步判断，如果|F(σ)|＞|F(σ

σ＝σ

(4)重复步骤八第(2)步-第(3)步计算10次，如果此时修正后的应力满足|F(σ)|≤FTOL，则退出循环。

(5)若循环次数小于10次，认为应力修正成功，更新应力σ

步骤九：由于发生应力重分布，弹性状态下的线性方程组求解问题变为弹塑性状态下的非线性方程组求解问题，采用牛顿-拉夫逊迭代进行求解，根据重分布后的应力和塑性流动法则求得弹塑性应力应变矩阵，并构建出新的整体刚度方程。

当高斯点发生屈服，应力应变曲线呈现非线性关系，通过塑性流动法则求得弹塑性应力应变矩阵：

再根步骤三的算法重新构建新的刚度矩阵，再根据牛顿-拉夫逊迭代求解该大型的非线性方程组，具体迭代流程图如图8所示，迭代计算式为：

Δu

ψ(u

P＝∫∫∫B

步骤十：根据塑性流动法则计算得到高斯点的塑性应变，同时随着迭代，单元结点位移、应变、应力和塑性应变在不断累积，直至达到所设置的最大迭代次数(视为边坡已发生破坏)或者满足收敛判据。

采用牛顿-拉夫逊迭代法设置的最大迭代次数为25次，同时设置位移和不平衡力的迭代收敛判据：

||Δu

其中Tol取0.00001。

在迭代过程中，位移、应力、应变和塑性应变会不断累积，其中塑性应变根据流动法则，

计算式为：

步骤十一：引入强度折减法进行上述计算，通过二分法搜索求得边坡的安全系数。

考虑边坡安全储备问题，引入强度折减法，获取边坡安全系数，折减计算式为

通过二分法不断变更折减系数，设置折减系数上下限，当在某一折减系数边坡发生破坏，更改折减系数上限值，反之更改折减系数下限值，不断二分直至折减系数上下限的差值小于误差限，通常设置为0.01。

具体为：在边坡安全系数的求解过程中，事先设置好一大一小两个折减系数用于二分法搜索(FSmin＝0，FSmax＝5，可根据不同边坡进行调整)，初始折减系数FS＝(FSmin+FSmax)/2，然后对根据此FS对强度参数c和φ进行折减，计算步骤一至步骤四，如果计算不收敛，认为此时折减系数FS过大，边坡已经发生破坏，令FSmax＝FS，反之计算收敛，认为此时折减系数FS过小，边坡还没有发生破坏，令FSmin＝FS，然后根据FS＝(FSmin+FSmax)/2计算出新的折减系数，再重复步骤一至步骤四，以此不断二分直至某一折减系数下边坡刚好处于极限状态，此折减系数即为边坡安全系数。本方法具体实现是当(FSmax-FSmin)/2≤err时，(err＝0.01，可根据需求调整计算精度)，认为此时的FS＝(FSmin+FSmax)/2即为边坡安全系数。

步骤十二：进行程序后处理，最终获取模型中各节点的位移，应变，应力和塑性应变等数据，可通过Tecplot绘图得到边坡的变形图和塑性区云图，结果展示如图9-10所示。

步骤十二计算结束之后可以获取每个单元高斯点的应力和塑性应变，为了更直观的体现出边坡各位置出的应力或应变状态，需要获取各节点的状态，先通过外推的方法获取每个单元节点的应力或塑性应变：

σ

然而此时获取的仅仅是每个单元节点的应力，仍需要通过加权平均的方法获得边坡各节点的应力：

求解单元塑性应变方法类似，通常为了更直观体现边坡塑性区的位置，在平面应变问题中，可将每个节点求得的四个塑性应变分量通过下式计算得到等效的塑性应变。

上述实施例阐明的方法，具体可以由计算机芯片或实体实现，或者由具有某种功能的产品来实现。一种典型的实现设备为计算机。具体的，计算机例如可以为个人计算机、膝上型计算机、蜂窝电话、相机电话、智能电话、个人数字助理、媒体播放器、导航设备、电子邮件设备、游戏控制台、平板计算机、可穿戴设备或者这些设备中的任何设备的组合。

为了描述的方便，描述以上装置时以功能分为各种单元分别描述。当然，在实施本说明书时可以把各单元的功能在同一个或多个软件和/或硬件中实现。

本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

在一个典型的配置中，计算设备包括一个或多个处理器(CPU)、输入/输出接口、网络接口和内存。

内存可能包括计算机可读介质中的非永久性存储器，随机存取存储器(RAM)和/或非易失性内存等形式，如只读存储器(ROM)或闪存(flash RAM)。内存是计算机可读介质的示例。

由此，本申请提供一种计算机可读存储介质，所述存储介质存储有计算机程序，所述计算机程序被处理器执行时实现前述之一所述的方法。

本申请还提供一种电子设备，包括存储器、处理器及存储在所述存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现前述之一所述的分析方法。

计算机可读介质包括永久性和非永久性、可移动和非可移动媒体可以由任何方法或技术来实现信息存储。信息可以是计算机可读指令、数据结构、程序的模块或其他数据。计算机的存储介质的例子包括，但不限于相变内存(PRAM)、静态随机存取存储器(SRAM)、动态随机存取存储器(DRAM)、其他类型的随机存取存储器(RAM)、只读存储器(ROM)、电可擦除可编程只读存储器(EEPROM)、快闪记忆体或其他内存技术、只读光盘只读存储器(CD-ROM)、数字多功能光盘(DVD)或其他光学存储、磁盒式磁带，磁带磁磁盘存储或其他磁性存储设备或任何其他非传输介质，可用于存储可以被计算设备访问的信息。按照本文中的界定，计算机可读介质不包括暂存电脑可读媒体(transitory media)，如调制的数据信号和载波。

本领域技术人员应明白，本说明书的实施例可提供为方法、系统或计算机程序产品。因此，本说明书可采用完全硬件实施例、完全软件实施例或结合软件和硬件方面的实施例的形式。而且，本说明书可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本说明书可以在由计算机执行的计算机可执行指令的一般上下文中描述，例如程序模块。一般地，程序模块包括执行特定任务或实现特定抽象数据类型的例程、程序、对象、组件、数据结构等等。也可以在分布式计算环境中实践本说明书，在这些分布式计算环境中，由通过通信网络而被连接的远程处理设备来执行任务。在分布式计算环境中，程序模块可以位于包括存储设备在内的本地和远程计算机存储介质中。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。