等参单元法与滑动面法相结合的土坡地震稳定性分析方法与流程

本发明属于地质工程技术领域，具体涉及一种等参单元法与滑动面法相结合的土坡地震稳定性分析方法。

背景技术：

基于当前对于土坡进行地震稳定性分析时，大多都是采用的等效荷载法，该方法能够有效考虑地震作用时的动力特性以及结构的动力特性，但是其忽略了重要的材料问题，土材料自身的动力特性，即应力应变的滞后性并没有在抗震分析中加以考虑。

大量的试验结果表明，土体材料由于其独特的颗粒特性，在经历动荷载作用下，土体的应力应变之间存在着明显的非线性，且由于土体阻尼的存在表现出强烈的滞后性。伴随着动荷载的强度以及持续时长的改变、增加，土体的强度明显降低，乃至发生失稳液化。因而土材料本身的动力特性就成了影响土坡稳定性的重要因素，如何将地震、土坡结构以及材料的动力特性同时在土坡稳定性分析中加以考虑，是现如今常用的等效荷载法需要突破的问题。

技术实现要素：

本发明的目的是提供一种等参单元法与滑动面法相结合的土坡地震稳定性分析方法，解决了现有技术中存在的在土坡稳定性分析中不能将地震、土坡结构以及材料的动力特性同时考虑的问题。

本发明所采用的技术方案是，一种等参单元法与滑动面法相结合的土坡地震稳定性分析方法，具体按照以下步骤实施：

步骤1、分级施加静荷载，并按每级静荷载施加后土中应力应变的变化情况，修正力学参数，按照实际的应力应变特点构建本构模型作为应力-应变关系的数学模型，当所有静荷载施加后，计算出的累积应力即为震前土坡中的静应力；

步骤2、根据土坡抗震设计要求，确定地震加速度时程关系，然后采用二级时段法输入；

步骤3、静力与动力反应计算，求得土坡内某些点的有效应力。

本发明的特点还在于，

步骤1中本构模型为邓肯-张模型或者弹塑性本构模型。

步骤2中二级时段法中，一级时段是指地震加速度曲线发生转折的相邻2个转折点之间的时段；二级时段也称为计算时段，是指将一级时段按积分稳定性和计算精度要求，均匀等分的每个等份时段。

步骤3具体按照以下步骤实施：

步骤3.1、在一级时段内循环；

步骤3.2、假定单元的初始动剪切模量gi和阻尼比λi，i指第i个单元；

步骤3.3、对该一级时段内的所有二级时段循环，计算出该一级时段内的动剪应变时程关系及平均动剪应变γ；

步骤3.4、确定该一级时段的等效周数及液化次数；

步骤3.5、计算该一级时段结束时的动孔隙压力和有效应力。

步骤3.3具体按照以下步骤实施：

步骤3.3.1、求解每个二级时段内的动力方程：

[m]{δ″}+[c]{σ′}+[k]{δ}＝{r(t)}(1)

式中：[m]、[c]和[k]分别为总体质量、阻尼和刚度矩阵，{r(t)}为结点动荷载向量，{δ″}、{σ′}和{δ}分别为结点相对加速度、速度和位移，[c]用瑞利公式计算；

步骤3.3.2、根据平均剪应变γ和平均有效主应力计算动剪切模量g0和阻尼比λ：

λ＝λmaxγh/(1+γh)(3)

γh＝γ[1+a·exp(-bγ/γr)]/γr(4)

式中，λmax、a、b为试验参数，pa为大气压力，它和g0、应有相同的因次，f(e)为与孔隙比有关的函数，建议关系式ocr为超固结比，响程度与塑性指数的大小有关，n为应力的指数，取0.5，γr为参考应变，γh为计算参数，为平均有效主应力；k为固结应力比的指数，取值与塑性指数的大小有关；

步骤3.3.3、如果该一级时段计算出的动剪切模量gi满足下式(5)，则结束迭代，进入下一个一级时段：

gi,j-gi,j-1≤δ(5)

式中：i指第i个单元，j为该一级时段内第j次迭代，j-1为该一级时段内的第j-1次迭代，δ为计算允许误差；

进入下一个一级时段计算时，该时段内的初始剪切模量g和阻尼比λ则由以下公式进行计算：

其中，k1、k2为试验参数，γc为计算参数；γd为动应变，地震最大剪应力沿深度的减小系数；λmax为上一个一级时段的阻尼比最大值；σ0为初始有效应力。

步骤3.4中，确定一级时段的等效周数采用martin方法，确定一级时段的液化次数由液化曲线查得。

步骤3.5中，孔隙压力计算模式通过动力试验确定。

本发明的有益效果是,一种等参单元法与滑动面法相结合的土坡地震稳定性分析方法，弥补了现有等效荷载分析法上的不足，将等参元法与滑动面法相结合，进而实现了在研究土坡抗震稳定性时，不仅考虑了地震的动力特性以及土坡结构的动力特性，还能兼顾材料的动力特性。其优点是：能够在任意时刻确定土坡最危险滑动面的位置与安全系数，同时能考虑土的应力应变特性。

附图说明

图1(a)是本发明一种等参单元法与滑动面法相结合的土坡地震稳定性分析方法中交点坐标确定图；

图1(b)是本发明一种等参单元法与滑动面法相结合的土坡地震稳定性分析方法中滑弧子段中点坐标的确定图；

图2是本发明一种等参单元法与滑动面法相结合的土坡地震稳定性分析方法中滑弧子段法向和切向应力的确定图。

具体实施方式

下面结合附图和具体实施方式对本发明进行详细说明。

本发明一种等参单元法与滑动面法相结合的土坡地震稳定性分析方法，具体按照以下步骤实施：

步骤1、分级施加静荷载，并按每级静荷载施加后土中应力应变的变化情况，修正力学参数，按照实际的应力应变特点构建本构模型作为应力-应变关系的数学模型，当所有静荷载施加后，计算出的累积应力即为震前土坡中的静应力；

其中，本构模型为邓肯-张模型或者弹塑性本构模型；

步骤2、根据土坡抗震设计要求，确定地震加速度时程关系，然后采用二级时段法输入，

其中，二级时段法中，一级时段是指地震加速度曲线发生转折的相邻2个转折点之间的时段；二级时段也称为计算时段，是指将一级时段按积分稳定性和计算精度要求，均匀等分的每个等份时段；

步骤3、静力与动力反应计算，求得土坡内某些点的有效应力，具体按照以下步骤实施：

步骤3.1、在一级时段内循环；

步骤3.2、假定单元的初始动剪切模量gi和阻尼比λi，i指第i个单元；

步骤3.3、对该一级时段内的所有二级时段循环，计算出该一级时段内的动剪应变时程关系及平均动剪应变γ；

步骤3.4、确定该一级时段的等效周数及液化次数，其中，确定一级时段的等效周数采用martin方法，确定一级时段的液化次数由液化曲线查得；

步骤3.5、计算该一级时段结束时的动孔隙压力和有效应力，其中，孔隙压力计算模式通过动力试验确定。

步骤3.3具体按照以下步骤实施：

步骤3.3.1、求解每个二级时段内的动力方程：

[m]{δ″}+[c]{σ′}+[k]{δ}＝{r(t)}(1)

式中：[m]、[c]和[k]分别为总体质量、阻尼和刚度矩阵，{r(t)}为结点动荷载向量，{δ″}、{σ′}和{δ}分别为结点相对加速度、速度和位移，[c]用瑞利公式计算；

步骤3.3.2、根据平均剪应变γ和平均有效主应力计算动剪切模量g0和阻尼比λ：

λ＝λmaxγh/(1+γh)(3)

γh＝γ[1+a·exp(-bγ/γr)]/γr(4)

式中，λmax、a、b为试验参数，pa为大气压力，它和g0、应有相同的因次，f(e)为与孔隙比有关的函数，建议关系式ocr为超固结比，响程度与塑性指数的大小有关，n为应力的指数，取0.5，γr为参考应变，γh为计算参数，为平均有效主应力；k为固结应力比的指数，取值与塑性指数的大小有关；

步骤3.3.3、如果该一级时段计算出的动剪切模量gi满足下式(5)，则结束迭代，进入下一个一级时段：

gi,j-gi,j-1≤δ(5)

式中：i指第i个单元，j为该一级时段内第j次迭代，j-1为该一级时段内的第j-1次迭代，δ为计算允许误差；

进入下一个一级时段计算时，该时段内的初始剪切模量g和阻尼比λ则由以下公式进行计算：

其中，k1、k2为试验参数，γc为计算参数；γd为动应变，地震最大剪应力沿深度的减小系数；λmax为上一个一级时段的阻尼比最大值；σ0为初始有效应力。

本发明一种等参单元法与滑动面法相结合的土坡地震稳定性分析方法，土坡抗震稳定性验算如下：

通过验算可以确定该时刻土坡的最危险滑动面及相应的最小安全系数。搜寻最危险滑动面的方法有多种，如网格搜寻法和二分法。以下内容以某一滑动面s(圆心坐标为(xs，ys)，半径为rs)为例，给出其安全系数的计算方法。对于其他类型的滑动面，同样可以按照类似方法进行计算。计算步骤如下：

①找出相交单元：

如果滑动面与某单元的至少一个边相交，则称该单元为相交单元，反之为非相交单元。设滑动面s的圆心到某一单元第j个结点(j＝1,2,…,ne,ne为单元结点总数)的距离为lj,且δj＝lj-rs)。这样，如果该单元内至少同时存在一个δj1>0和δj2<0，则表示该单元与滑动面相交，若该单元的δj皆大于零，或皆小于零，则该单元为非相交单元。对滑动面进行地震稳定性验算，只与相交单元有关，与其他单元无关。

②建立相交单元内任意点处的有效应力计算公式：

对于8结点等参元，单元内任意点处的某一应力分量，例如σx，可采用hinton等人提出的如下方法计算:

ni＝(1+ξiξ)(1+ηiη)/4(9)

式中:(ξi，ηi)和σxi分别为单元第i个角点的局部坐标和x向应力。由于式(9)对单元内各点均适用，因此，将单元内2×2个高斯点处的局部坐标及x向应力代入，建立4个线性方程，进而计算出σxi这样，只要知道单元内某点的局部坐标(ξ，η)，即可计算出该点的σx。对其它应力分量σy和τxy,可仿照进行。

对于不同类型的单元，甚至复杂的8节点等参元，采用如下方法建立与滑面相交的单元内任意点处的应力计算公式。各应力分量可表示如下：

σx＝α1+α2x+α3y+α4x²+α5xy+α6y²+…(13)

σy＝β1+β2x+β3y+β4x²+β5xy+β6y²+…(14)

τxy＝γ1-β3x-α2y-0.5β5x²-2(α4+β6)xy-0.5α3y²+…(15)

式中：αi、βi、γi为待定参数。在一般情况下，我们选择前三项计算就可以了，过多的项数可能导致计算出的应力呈现虚假波浪分布形式。对于第i个高斯点，设整体坐标为由有限元计算出的有效应力则单元内ng个高斯点上应力的均方误差之和为建立最小二乘泛函:

对(16)式进行泛函变分，有：即：

式中：

通过(17)式可计算出相应的系数α1，α2和α3.这样，当已知单元内某点的整体坐标，即可计算解出对应的σx，σy，τxy。

③计算滑动面与相交单元各边的整体交点坐标：

以相交单元的某一角点为起点，按一定方向(如逆时针方向)对单元各边循环，用相邻二角点构成的直线方程与假定的滑动面方程联立求解，即可计算出滑动面与相交单元各边的整体交点坐标。显然，整个滑动面被相交单元分割成若干个滑弧段。这里，滑弧段是指相邻二个整体交点之间的滑弧。

④计算滑弧段上各子段的有效应力：

为了比较准确地计算出滑动面的安全系数，应按照弧的长短和应力变化的大小，将各滑弧段再分为若干个子段。每一子段上的应力可用该子段中点上的应力近似代表。

如附图1(a)所示。设一相交单元与滑弧形成的一滑弧段为ab。a和b两点的整体交点坐标为(xa，xa)和(xb，xb)。则滑弧段ab对应的圆心角α为：

式中：假如滑弧段ab按α均分为m个子段，第j个子段的起点aj和终点bj的坐标分别为(xj-1,yj-1)及(xj,yj),中点dj的坐标为则从图1(b)有:

将式(20)中的换为则式(19)中的即为bj点的坐标xj、yj.显然，对第一子段(j＝1),xj-1＝xa,yj-1＝ya；对最后一个子段(j＝m),xj＝xb，yj＝yb.并且,bj点即为下一子段的起点。

求出dj点的整体坐标后,如果采用式(13)～(15)计算该点处的应力σxj、σyj和τxyj,只需将代入。对8结点等参元,如果按式(1)计算dj处的应力,还需计算出dj点对应的局部坐标(ξj,ηj)。当按图1(a)的方式对单元进行节点编号，并且每条边上的角点之间的节点(如图1-图2边之间的节点5位于边中时，由公式：

可推出以下方程：

a1ξjηj+a2ξj+a3ηj-a6＝0(23)

式中：(xi，yi)和ni分别为单元第i个节点的整体坐标和形函数，a1＝x1+x3-x2-x4，a2＝x2+x3-x1-x4，a3＝x3+x4-x1-x2，a4＝x1+x2+x3+x4，a6＝a5-a4.将ak(k＝1～6)中的xi换为yi，换为可计算出bk.求解方程(22)和(23)，可得到dj点的局部坐标(ξj，ηj)。但应注意|ξj|≤1，|ηj|≤1.求出ξj和ηj后，代入式(9)可计算出dj点的σxj、σyj和τxyj.

⑤计算滑弧子段上的有效法向应力和切应力：

如前所述,滑弧子段上的应力可用该子段中点处的应力代表。为此，可通过在该点取微元体，根据力的平衡条件计算出该子段上的有效法向应力σnj和切向应力τj,即：

式中：θj＝2tan^-1[(yj-yj-1)/(xj-xj-1)]。

当在静力条件下计算滑动面的安全系数时，前述应力仅指由静荷载作用产生的静应力；在动力条件下计算某一时刻的安全系数时，前述应力是指该时刻静、动荷载共同作用下产生的应力。此外，所有正应力均指的是有效应力。

⑥计算滑动面的安全系数：

设滑动面s被相交单元分割n个滑弧段，第i个滑弧段又被分为mi个子段，则该滑动面的安全系数为：

式中:δlj和cj、φj分别为第j个子段的弧长和有效抗剪强度指标。当考虑液化问题时，如果子段中点处的平均有效应力等于零，则意味着该子段的土体已液化，故σnj＝0。

按照以上方法步骤，可计算出某一滑动面的安全系数。对其它假定的滑动面，逐个进行计算。从所计算出的这些滑动面的安全系数中寻找出最小者，该最小安全系数所对应的滑动面即为该时刻土坡的最危险滑动面。对地震过程中的其它时刻，按照上述方法继续计算，直到地震结束。