土工结构安全系数的有限元计算

第19卷 第2期岩 土 工 程 学 报Vo l.19 N o.2 1997年 3月Chinese Jour nal of Geotechnical Eng ineering M ar. 1997土工结构安全系数的有限元计算*宋二祥(清华大学土木系,北京,100084)文 摘 定义土工结构安全系数为其极限承载力与所需承载力之比,给出了按此定义计算土工结构安全系数的有限元法。

在计算中讨论了弧长控制法的应用。

作为算例,首先计算了一座土坝的安全系数,并与Bishop方法的计算结果相比较,二者相当吻合。

此外,还计算了用土工织物加强路基的安全系数,进一步说明了此法的可靠性及适用性。

关键词 安全系数,非线性分析,有限元法。

1 引 言对于一般地面建筑结构,其承载力的安全系数通常可表示为结构的极限承载力与结构在使用阶段所能承受的最大荷载之比。

但对于一些土工结构,比如路基、堤坝等,上述定义并不适用。

此类结构所受荷载主要是土的自重,为计算其极限承载力而逐渐增大土的容重时,土中正应力也相应增大。

而土是一种摩擦材料,在正应力增大时其抗剪强度也增大。

所以,对此类土工结构,增大土的容重未必会使其达到极限状态。

例如,用砂土筑成的土堤,其坡度角将由砂土的摩擦角唯一决定,而与砂土的容重无关,增大砂土的容重并不会使土堤发生滑坡而破坏。

此外,土的容重易于测定,且变异性不大,而其强度参数的测定要复杂得多。

所以,土工结构的安全系数一般定义为结构所具有的承载力与承受荷载所需要的承载力之比。

土力学中土坡稳定分析所给出的安全系数即按此定义。

土工结构安全系数的计算,特别是土坡稳定安全系数的计算,一般采用基于极限平衡理论的方法,这在一般土力学教科书中均有介绍。

这种方法包括两个方面:一是对设定的破坏滑移面计算安全系数,这方面的计算方法有Bishop方法、Janbu方法、Spencer方法等[1];二是寻找真正的破坏面以给出最小安全系数,这方面最简单的方法是取不同的破坏面进行计算、比较。

近年来又采用运筹学的一些方法,最简便的是非线性规划中的单纯形法,也有人用动态规划方法。

极限平衡方法的一个局限是对较复杂的土层及土工结构,计算较困难。

随着有限元法在岩土工程中的广泛应用,它也被用来分析土坡稳定的安全系数,但这方面的研究还较少。

Donald和Giam[2]曾用有限元法得到的结点位移来确定土坡稳定的安全系数,但他们的方法计算量很大。

最近Zou等人也采用有限元法分析土坡稳定安全系数[3],但他们是用有限元法计算土坡内各点的应力,而破坏面的确定是用动态规划法。

本文完全采用非线性有限元法计算土工结构如上定义的安全系数。

其基本思路是,在计国家教委回国留学人员基金资助项目。

到稿日期:1995-08-22.1算过程中将土的强度参数(内摩擦角和内聚力)逐步降低而使结构达到极限状态,土所具有的强度参数值与相应于该极限状态的强度参数值之比,则是所求的安全系数。

在此计算过程中,外部荷载保持不变,而土的强度逐步降低,这显然不同于一般建筑结构极限荷载的计算方法。

本文将给出用有限元计算此类问题的具体方法,特别要着重讨论在此计算过程中怎样利用弧长控制法来调整强度参数的降低量,以便精确得到结构的极限状态及相应的强度参数值。

本文给出的计算方法简便易行,计算速度快,稳定性好。

随后给出的两个算例也说明了这一点。

2 计算方法简介2 1 土的应力应变模式这里采用土工结构计算中广泛应用的理想弹塑性莫尔库仑模式,因为就土工结构的破坏分析来说该模式是相当简便而又可靠的。

该模式的屈服面方程可以写为1 2( 3- 1)+12( 3+ 1)sin -c cos =0(1)当主应力不按大小排列时,屈服面需用6个与上式类似的方程来描述。

在主应力空间中莫尔-库仑屈服面是一个六棱锥面。

决定塑性流动方向的塑性势面方程和屈服面方程有相同的形式,但需用剪胀角来代替内摩擦角,以给出合理的塑性体积应变。

在用初应力法进行弹塑性分析时,每次迭代我们首先得到试探应力为{ e}={ 0}+[D]{ }(2)当{ e}超出屈服面时,则用下式修正该应力使其回到屈服面:{ p}= i{ g i } { }={ e}-[D]{ p}(3)式中 g i为塑性势面方程。

这里我们采用Koiter一般流动法则, i由一致条件f( )=0决定。

将一致条件在{ e}点泰勒展开并略去高阶项,有i i{f j}T[D]{g i}=f j(e)(4)上式是 i的线性代数方程组,由此可求解 i。

具体计算时,为更为可靠有效,可在进行上述修正之前事先判定修正后的应力将位于屈服面的哪一区域。

详细介绍见文献[4]。

这种应力应变模式的积分方法为隐式方法,它与求解总体平衡方程的初应力法相结合,给出无条件稳定的算法。

如引言所述,土工结构的安全系数应定义为其极限承载力与所需承载力之比。

当采用莫尔库仑模式时安全系数F s可表示为F s=(c a+ tg a)/(c r+ tg r)(5)这里下标 a 和 r 分别表示 所能提供的 和 所需要的 。

为简化计算,这里假定c和 按同一比例变化,从而F s也等于c a与c r之比,或tg a与tg r之比。

2 2 基本计算方法为计算一个土工结构的安全系数,我们首先按所给土的力学性质参数(c a, a等)进行计算,以施加结构的全部荷载。

然后逐步减小土的强度参数直至结构破坏,从而求得c r与 r。

每次减小土强度参数后的迭代计算方法与一般弹塑性有限元类似,该迭代过程中典型的应力2岩 土 工 程 学 报 1997年变化大致如图1中的折线ABC 所示。

图中L 1,L 2分别对应于强度参数降低前后的屈服面。

在每次减小强度参数后,对于超出新屈服面的应力(如图1中A 点应力),将按上节所述方法修正到屈服面(如图1中B 点应力)。

这样修正后的应力将不再符合平衡等条件,所以需进行迭代。

在迭代过程中应力在新的屈服面上调整,以使整个区域的内力和变形符合连续体的平衡、几何和物理三方面的条件。

其数学表述可推导如下。

图1 基本迭代法中的应力变化 图2 使用弧长法时的应力变化由有限元平衡方程[B ]T { }d V ={F}(6)式中 {F }是外部荷载向量。

设在迭代过程中有{ i }={ i -1}+{ i },则可得迭代方程[B ]T { i }d V ={F }- [B ]T { i-1}d V (7)在开始迭代时(i =1)取{ 0}为经过修正而符合当前屈服条件的应力。

随迭代次数的增加,应力的子增量{ i }将趋于0。

采用初应力法时,记[K e ]= [B ]T [D][B ]d V ,则上式可写成[K e ]{ d i }={F }-[B ]T { i-1}d V (8)式中 d i 为位移子增量。

这种迭代方法可避免不平衡误差的积累。

2 3 弧长控制法的应用当结构趋于破坏时,在强度参数不变的情况下变形持续发展。

如分别以强度参数比和位移为纵横坐标画曲线,则该曲线趋于水平。

如采用上述基本算法,当土的强度参数在某一步中减小过多,则结构不可能达到平衡,从而使迭代失败,不能求出准确的安全系数值。

为解决此问题,这里采用在Riks 方法基础上发展的一种弧长控制法[5,6]。

弧长法是计算结构破坏荷载及跟踪结构变形曲线软化段的一种方法。

其基本原理是在迭代过程中调整荷载增量的大小,以得到收敛的解。

荷载调整的原则是使位移向量增量的内积在迭代过程中保持为常量。

但是,在安全系数的计算过程中如直接采用这种方法,则意味着在迭代过程中不断改变强度参数的减小量,这样屈服面随迭代过程不断变化。

如此进行计算,收敛性将很差。

为使计算更为有效,这里进行一些简化处理。

如图2所示,设L 1,L 2和L 3分别代表初始屈服面、试探屈服面和迭代结束后的屈服面,我们将应力点沿AD 的移动看成沿折线ABCD 的移动。

在迭代过程中屈服面用L 2,但用弧长法所求出的荷载乘子m 修正所求出的应力,使其近似对应于最终屈服面L 3上的应力点D 。

这样式(8)右端的应力此时可写为{ i-1}={ i-1c }+(1-m i-1)({ A }-{ B })(9)这里下标A ,B ,C 给出应力所对应的点,m 是荷载乘子,在这里为材料强度参数降低量中所3 第2期宋二祥 土工结构安全系数的有限元计算实现的部分。

上式右端第二项是使用弧长法时对试探屈服面上应力的修正,此种修正可在迭代结束之后进行。

由于在新屈服面上的应力点并非总是对应于原来屈服面上的某一点,式(9)的修正也并非精确,故这里的处理方法是近似的。

但如果我们注意到在结构接近破坏时,土的强度参数几乎不再变化,则不难理解上述处理方法的最终误差将是很小的。

将式(9)代入式(8),有[K e]{ d i}={F}-(1-m i)({F}- [B]T{ B}d V)- [B]T{ i-1c}d V(10)在上式的推导中利用了 [B]T{ A}d V {F}的条件。

记m i=m i-1+ m i,则式(10)可写为[K e]{ d i}={R i-1}+ m i{ R}(11)其中{R i-1}={F}-(1-m i-1)({F}- [B]T{ B}d V)- [B]T{ i-1c}d V{ R}={F}- [B]T{ B}d V这里我们采用一种简化的弧长法,它要求相邻两次迭代所求得的位移向量增量的内积为常数,即{ d i-1}T{ d i}=l2(12)注意到{ d i-1}T{ d i-1} l2,可有{ d i-1}T{ d i}=0(13)将式(11)代入上式可以导出迭代过程中依据不平衡力的大小调整荷载乘子的公式为m i=m i-1-{ d i-1}[K e]-1{R i-1}{ d i-1}[K e]-1{ R}(14)对于材料非线性问题,上述方法的一个重要修正是用部分结点的位移来代替其中的整体结点位移向量,以便更准确地探知结构破坏或软化过程中的应变集中[6]。

计算是分步进行的。

在每一步中将强度参数降低之后,我们保持式(12)中的l不变,而据不平衡力的大小按式(14)对强度参数降低量进行调整,当结构破坏时我们即可得到接近水平的强度参数比-位移曲线。

由于采用理想弹塑性本构模式,再利用上述的弧长法,一般总可得到该曲线的水平段,从而定出结构的安全系数。

3 土坝稳定问题的安全系数(算例1)上述计算方法已用于作者与荷兰几位同事合作开发的土工有限元程序PLAXIS,经多次计算表明所采用的方法是相当可靠和有效的。

作为本文的第一个算例,这里计算一座土坝的安全系数。

图3为所计算土坝的示意图,设其左右两侧的水位高度分别为3m和7m,该水位差引起的坝体内渗流水自由表面用变渗透系数的有限元法计算求得[7]。

取土的内摩擦角 = 30 ,内聚力c=10kPa。