（水工结构工程专业论文）岩土边坡的有限元稳定性分析.pdf

摘要 基于强度折减的有限元稳定分析方法是目前应用及研究较多的一种分析方法， 如何根据有限元计算结果来判别边坡是否处于整体破坏状态，是强度折减有限元稳 定分析方法的一个关键性问题。 本文首先从非线性有限元方程的求解入手，介绍了求解非线性方程的迭代法、 增量法及混合法，对各自的基本理论、求解过程及特点作了分析。并对有限元计算 中采用的位移收敛准则、不平衡力收敛准则和能量收敛准则做了阐述。 有限元强度折减法的失稳判别标准主要有三种：第一种以有限元解的收敛性确 定失稳状态，第二种根据计算域内最大节点位移与折减系数之间关系曲线的变化特 征确定失稳状态，第三种通过计算域内塑性区是否贯通确定失稳状态。本文利用 a d i n a 有限元通用软件，对三个不同的算例进行强度折减计算，然后分别采用三种 失稳判别标准进行稳定性分析，对三种失稳判别标准适用性、三者之间的一致性、 各自的适用范围进行了分析研究。对于均质边坡，三种失稳判别标准存在较好的一 致性，但第三种失稳判别标准在应用范围上存在局限性。之后对计算精度进行了分 析，主要分析了网格密度、容许限值对计算精度的影响。网格太稀疏、容许限值设 太大都会使计算结果的精度降低，在稳定性判别时就容易出现偏差。 最后，对一实际边坡工程进行有限元计算，对其应力、变形进行了分析，并利 用前面分析研究的三种失稳判别标准对其稳定性进行了判定。在计算工况下，边坡 最大拉应力出现在坡面，由于边坡表层为强风化岩石，其强度很低，可能出现拉剪 破坏。交通洞四周岩体强度较低，且弹模较小，变形较大，可能发生压剪破坏。结 合三种失稳判别标准判定该边坡在计算工况下可能发生深层滑动破坏，需采取一定 工程措施进行加固。 关键词：岩土边坡强度折减非线性有限元稳定分析失稳判别标准收敛准则 a d i n a a b s t r a c t t h es t a b i l i t ya n a l y s i sm e t h o d su s i n gf e m ，w h i c hb a s e do ns t r e n g h tr e d u c t i o n ，h a s b e e nm a d eal o to fa p p l y i n ga n dr e s e a r c h i n gi ni t sf i e l d s h o wt oj u d g et h ec a l c u l a t e d s l o p e ss t a t e s ，s t a b i l i t yo ri n s t a b i l i t y , i sak e yp r o b l e mf o rs h e a rs t r e n g t hr e d u c t i o n m e t h o d s f i r s t l y , t h i sp a p e ri n t r o d u c e dt h es o l u t i o nm e t h o d so f n o n l i n e a rs t a t i ce q u a t i o n sa n d i t sb a s i ct h e o r y , t h es o l v es t e p sa n dc h a r a c t e r i s t i c s e x c e p tf o rt h o s e ，t h r e ec o n v e r g e n c e c r i t e r i aw e r ei n t r o d u c e d t h e r ea r et h r e em a j o rc r i t e r i o n sf o re v a l u a t i n go ff a i l u r es t a t e t h ei n s t a b i l i t ys t a t ei s j u d g e db yt h ec o n v e r g e n c eo fe q u a t i o n ss o l u t i o ni nt h ef i r s tc r i t e r i o n t h ec u r v ew h i c h e x p r e s s e st h er e l a t i o n s h i po fn o d el a r g e s td i s p l a c e m e n tw i t hr e d u c t i o nc o e f f i c i e n ti st h e j u d g m e n to ft h es e c o n dc r i t e r i o n t h ec o n n e c t i v i t yo fp l a s t i cf l a gi st h ej u d g m e n to ft h e t h r e ec r i t e r i o n s t h r e ed i f f e r e n tc a s ew e r ec a l c u l a t e du s i n gt h eg e n e r a lf e ms o f t w a r e a d i n a b a s e do nt h er e s u l t so ff e m t h i sp a p e ra n a l y z e dt h es t a b i l i t yo fc a l c u l a t e d s l o p es e p a r a t e l yb yt h r e ef a i l u r ee v a l u a t i o nc r i t e r i o n s b e s i d e s ，t h ei n f l u e n c ef a c t o r sf o r a c c u r a c yo fr e s u l t sw e r ea n a l y z e d t h ed e n s i t yo fg r i da n dt o l e r a n c eo fc o n v e r g e n c e c r i t e r i ah a v er e l a t i v e l yb i gi n f l u e n c ef o rt h ea c c u r a c yo fc a l c u l a t e dr e s u l t s f i n a l l y , t h i sp a p e rm a d eu s eo ff e mt oa n a l y z et h es t r e s ss t a t ea n dt r a n s f o r m a t i o n s t a t eo fas l o p ei np r a c t i c e t h e nc o m b i n i n gt h r e ef a i l u r ee v a l u a t i o nc r i t e r i o n st oa n a l y z e t h es l o p e ss t a b i l i t y , t h ea n a l y s i sr e s u l ti st h a tt h es l o p ei si n s t a b i l i t y s ot h e r en e e dt ot a k e s o m ee n g i n e e r i n gm e a s u r e st or e i n f o r c et h es l o p e k e y w o r d ：r o c ka n ds o i ls l o p e s ，s h e a rs t r e n g t hr e d u c t i o n ，n o n l i n e a rf i n i t ee l e m e n t ， s t a b i l i t ya n a l y s i s ，f a i l u r ee v a l u a t i o nc r i t e r i o n ，c o n v e r g e n c ec r i t e r i a ，a d i n a 学位论文独创性声明： 本人所呈交的学位论文是我个人在导师指导下进行的研究工作及 取得的研究成果。尽我所知，除了文中特别加以标注和致谢的地方外， 论文中不包含其他人已经发表或撰写过的研究成果。与我一同工作的同 事对本研究所做的任何贡献均已在论文中作了明确的说明并表示了谢 意。如不实，本人负全部责任。 论文作者( 签名) ：险邀6 4 年5 月2 歹e l 学位论文使用授权说明 河海大学、中国科学技术信息研究所、国家图书馆、中国学术期刊 ( 光盘版) 电子杂志社有权保留本人所送交学位论文的复印件或电子文 档，可以采用影印、缩印或其他复制手段保存论文。本人电子文档的内 容和纸质论文的内容相一致。除在保密期内的保密论文外，允许论文被 查阅和借阅。论文全部或部分内容的公布( 包括刊登) 授权河海大学研究 生院办理。 论文作者( 签名) ：叁丝。6 年占月矽 e l 第一章绪论 第一章绪论 1 1 有限元边坡稳定分析方法研究现状 边坡的稳定性问题一直是土木工程中一个重要的问题，它己变成同地震和火山并 列的全球性三大地质灾害( 源) 之一，也就是广为人知的滑坡、泥石流。这种破坏，轻 者引起边坡的破坏，重者摧毁边坡下的构筑物，引起重大的经济损失，甚至造成人员 伤亡。因此探寻一种能够较准确地分析边坡的稳定性、减少滑坡灾害发生的方法是一 件十分有益的事情。 在边坡稳定性分析中常采用的方法是传统的刚体极限平衡法，如瑞典条分法、 j a n b u 法、b i s h o p 法、s a r m a 法、s p e n c e r 法、通用条分法等n q l ，其理论基础是将滑 动区域可能的滑动体视为刚体，设在滑动面( 带) 上的岩土体处于塑性极限平衡状态， 而后利用刚体力学的观点分别计算滑动体所受的抗滑力( 或力矩) 和滑动力( 或力矩) ， 利用二者的比值作为边坡稳定性的衡量标准，其值即为安全系数。在计算时作了许多 假定陋 1 1 | ： ( 1 ) 假定了一定形状的滑动面； ( 2 ) 滑动面上的所有点同时达到极限状态； ( 3 ) 沿整个滑面上的安全系数是常量； ( 4 ) 土条是一系列刚性体； ( 5 ) 滑面上的抗剪强度采用了摩尔一库仑准则，滑面上的土体为刚塑性材料。 该法的优点在于理论基础简单，便于实现，且得出的安全系数易于为人们所接受， 但其缺点是忽视了边坡岩土体本身是变形体，边坡的破坏过程是其内部应力分布和变 形不断调整的过程。鉴于此，利用近现代发展起来的基于弹塑性理论的各种数值分析 方法计算边坡稳定性已受到了人们的重视。在诸多数值分析方法中有限元得到了极大 的应用，其计算理论已从早期的弹性发展到弹塑性或粘弹塑性。几何非线性有限元方 法能够较准确地模拟岩土体实际地质状况，并能考虑岩土体的非线性本构关系，可以 比较真实地反映边坡的实际变形和应力分布规律。但由于有限元法的分析成果一般是 应力、应变和节点处的位移，故其对边坡稳定性的评价不如极限平衡法直观，为此许 多研究者对有限元在边坡稳定性研究中的应用进行了大量的研究工作。 河海人学硕f ：论文 1 9 7 5 年，z i e n k i e w i c z n 2 1 等首次在土工弹塑性有限元数值分析中提出了抗剪强度 折减系数概念，由此所确定的强度储备安全系数与b is h o p 在极限平衡法中所给出的稳 定安全系数在概念上是一致的，且同时还可考虑边坡的变形和其内部塑性区的发展变 化情况。与传统的刚体极限平衡法和极限分析法相比，有限元边坡稳定分析方法有以 下优点口3 叫引： ( 1 ) 能够全面满足静力平衡、变形协调以及应力、应变之间的本构关系，不必 引入假定条件，保持了严密的理论体系。 ( 2 ) 基本不受边坡几何形状的不规则和材料的不均匀性的限制，可以较真实地模 拟边坡的地形地貌以及边坡内复杂的地质条件。 ( 3 ) 无需事先假定破坏面的形状和位置，破坏很“自然地”发生在边坡岩土体抗 剪强度不能抵抗剪应力的位置。 ( 4 ) 可以分析边坡破坏的发生和发展过程，模拟边坡开挖及加固的施工过程，考 虑岩土体与支挡结构的共同作用及其变形协调。 ( 5 ) 有限单元法可以提供应力、应变、位移等变量的全部信息。 正是由于有限元法上述优点，边坡稳定分析的有限元法近年来受到众多学者的关 注。但由于传统的有限元法无法直接得到工程实践广泛应用的安全系数和滑动面形状 及位置，因此如何将其计算成果与传统的边坡稳定安全系数联系起来，已成为边坡稳 定有限元中的一个重要研究课题，不少学者进行了比较全面的研究。目前根据有限元 法分析结果来评价岩土体稳定性的方法大致可以分为两类： 第一类是采用刚体极限平衡法原理分析有限元计算结果，和刚体极限平衡法的计 算步骤一样，根据一定的安全系数定义方法来寻求最小安全系数和最危险滑动面，这 种方法可以称之为间接法。该方法与极限平衡法计算步骤相似，先确定求解安全系数 的函数关系式厂= 厂 y ( x ) ，然后寻找最小安全系数兀2 呀署厂 少( x ) 】= f y e ( x ) 】和 最危险滑动面位置n 6 2 。z i e n k i e w i c z 等( 1 9 7 5 ) 提出了基于应力水平的安全系数定义 法，安全系数表示为当前应力圆直径与破坏应力圆直径之比，反映了边坡内岩土体强 度的发挥程度，即： 胁糍 m ， d o n a l d 和t a n ( 1 9 8 5 ) 睢圳在以上两种方法的基础上，综合考虑剪应力和应力水平两方面 第一章绪论 的因素，提出了基于应力水平加权强度的安全系数定义法，即： 名。： 兰! 兰芝! 垒呈笙! 竺一( 1 - 2 ) 腰3 一( c 7 + 仃7 t a n # 7 ) ( 叫一) j ( 一) 厂 与极限平衡方法一样，在对安全系数的确定方法进行研究的同时，也有不少学者致力 于研究寻求最小安全系数和最危险滑动面的方法。g i a m 和d o n a l d ( 1 9 8 8 ) 1 提出了由已 知应力场确定最小安全系数和最危险滑动面的方法。z o u 等( 1 9 9 5 ) 瞳们通过有限元法获 得的应力分布规律确定滑动面的范围和初始滑动面，然后利用动态规划的数值方法搜 索最小安全系数及对应滑动面。河海大学殷宗泽、吕擎锋心神等人提出圆弧滑动有限元 分析法，这种基于有限元应力变形计算的边坡稳定分析方法，假定滑动面形状为圆弧 形，有限元网格由一组同心圆作为纬线，一组竖向线为经线构成。对两相邻圆弧线所 夹的弧形带分析滑动力和抗滑力，建立力矩平衡方程，确定安全系数： ( t a n 矽+ c ) i ， 兀= 旦_ o 乞 j = l ( 1 - 3 ) 其中最小安全系数对应的弧形带为可能的滑动带。变化圆心位置用优选方法寻找最小 安全系数对应的圆心，从而得出真正的滑动面。圆弧滑动有限元分析法突出优点是由 有限元计算直接得出滑面上的应力，而不须作近似的插值处理，因而应力更准确。该 方法可以考虑土的非线性变形特性，也更符合土的实际情况。此外，用有限元计算得 出位移，亦可将稳定分析和变形联系起来，为现场通过位移监测来估计边坡的稳定性 提供了可能性，同时也为膨胀土边坡的稳定分析中考虑膨胀性的影响提供了可能性。 但此方法仅仅研究了圆弧滑动的情况，对于非圆滑滑动面还要进一步研究，同时圆弧 滑动有限元方法还要重复进行有限元计算，工作量大。对于填筑和开挖的模拟、膨胀 土和黄土浸水变形的模拟虽然在道理上是可行的，但也有待深入研究。 另一类方法是直接法，这种方法直接使用有限元法，通过不断降低边坡岩土体强 度或增加岩土体的自重使边坡岩土体达到临界状态，得到边坡的安全系数。采用这类 方法无需事先假定滑动面的形状和位置，通过不断降低岩土体的强度或增加岩土体的 自重，破坏将很“自然地”发生在边坡岩土体抗剪强度不能抵抗剪应力的位置，从而 得到最危险滑动面及相应的安全系数。这两种方法分别称之为强度折减有限元法和容 重增加有限元法。 一3 一 河海人学硕一 ：论文 1 强度折减有限元法 d u n c a n ( 1 9 9 6 ) 瞳6 1 指出边坡安全系数可以定义为使边坡刚好达到临界破坏状态时， 对土的剪切强度进行折减的程度。这种强度折减技术应用到有限元中可以表述为：保 持岩土体的重力加速度为常数，通过逐步减小抗剪强度指标，将c 、矽值同时除以折 减系数，得到一组新的强度指标c 、矽进行有限元分析，反复计算直至边坡达到 临界破坏状态，此时采用的强度指标与岩土体原具有的强度指标之比即为该边坡的安 全系数： c ，：一c ( 卜4 ) s r 矽a r c t a n ( t a nc f , ，) ( 卜5 ) 赵尚毅、郑颖人等乜7 3 通过比较毕肖普法和强度折减法的安全系数定义，认为两者 安全系数具有相同的物理意义，强度折减法在本质上与传统方法是一致的。采用强度 折减法分析边坡稳定时，通常将c 、矽值同时除以折减系数进行折减即采用等比例强 度折减的方法。国内外试验研究的结果证明，矽值较稳定，波动较小，而c 值受外界因 素影响较大，不够稳定，波动较大，因此把c 、矽值同等看待按等比例强度折减，显 然是不够合理，也不符合实际的，而采用不等比例强度折减则较合理且符合实际，不 等比例强度折减中的不等比例如何选用又是一个值得研究解决的课题。 陆述远、常晓林瞳8 2 们等在重力坝坝基稳定分析中采用等保证率法，即按等保证率 选取c 、矽值，随着保证率的增大，c 、矽值逐渐降低，得到两个不同的折减系数磁 和巧值，随着k c 和k 值逐渐增大，引起坝体从局部到整体的渐进破坏。但不等比 例强度折减方法分析边坡稳定仍然存在诸如c 、矽值按何种不同比例折减、如何根据 两个不同折减系数k c 和心值来评价边坡的稳定性等一系列的问题，这些问题还有待 进一步研究。 2 容重增加有限元法 容重增加法保持岩土体的抗剪强度指标c 、彩为常数，通过逐步增加重力加速度g 的方式，反复进行有限元分析，直至边坡达到临界破坏状态，此时采用的重力加速度 第一章绪论 g h m 与实际重力加速度c o 之比即为该边坡的安全系数，即： 1 九： - y l i m i t ( 卜6 ) ，a 一 lv g o c h e n 着h m i z u n o ( 1 9 9 0 ) 口叫采用容重增加法研究了粘性土边坡的稳定性，s w a n 和 s e o ( 1 9 9 8 ) 口1 3 2 1 在此基础上进行了比较系统地研究。对于摩擦角较大的平缓边坡应用容 重增加法时应特别注意，重力增加时土体中的平均正应力比剪应力增长得快，因此土 体强度增长会超过剪应力增长速度，土体可能不会破坏。 强度折减法和容重增加法分别通过强度折减系数和加速度比值来确定安全系数， 对于滑动面目前通过某一物理变量的等值线图或是分布云图来确定，或是在可能的滑 动区域内进行最危险滑动面的搜索。文 3 3 利用贯通的广义塑性剪应变的等值图来定 义滑动面。对于人工边坡，文 3 4 建议采用位移增量等值线来确定潜在滑面。文 3 5 通过使用非关联流动法则，将剪胀角沙取为0 ，发现在变形后的网格中会出现一条明 显的畸变带，他们将这条畸变带就定义为潜在滑面。取为0 意味着完全忽略了岩土 材料的剪胀特性，而仅突出其剪切变形。果然，在使用他们的程序时发现，取非零的沙 后这条畸变带并不明显；而且，即使将沙取为0 ，所使用的网格也必须相当规则，否 则，也难以出现畸变带。根据应力计算结果，文 2 4 利用动态规划法，文 3 6 利用人 工智能型优化方法的蚂蚁算法，讨论了基于瑞典法安全系数( 抗滑力比滑动力) 概念的 潜在滑面的搜寻问题。文 3 7 提出的基于能量准则的岩体稳定性分析方法是一个在概 念上有别于经典极限平衡法的分析方法。该方法是通过干扰能量等值线图来确定潜在 的滑面。 强度折减法和容重增加法的思路清晰，原理均比较简单，但是在实际问题中，边 坡所处的受力状态往往很复杂，不只是考虑自重的影响，所以容重法在实际应用中存 在一定局限性，而强度折减法的应用范围则相对较广。 1 2 问题的提出 强度折减法在本质上与传统的极限平衡法是一致的，鉴于它所具有的一些优点， 河海犬学硕一l ：论文 逐渐受到人们的关注，并在实际工程中得到广泛应用口8 4 引。采用有限元强度折减法分 析边坡稳定性的一个关键问题就是如何根据有限元计算结果来判别边坡是否处于整体 破坏状态。目前，岩土边坡破坏的标准主要有： ( 1 ) 根据有限元解的收敛性确定失稳状态，即在给定的非线性迭代次数及限值条 件下，最大位移或不平衡力的残差值不能满足所要求的收敛条件，则认为边坡在所给 定的强度折减系数下失稳破坏。对于材料的非线性，经过有限元离散后，问题归结为 求解一组非线性代数方程。一般的，将非线性方程问题转化为一系列线性问题，通过 迭代法或增量法使一系列线性解收敛于非线性解。在迭代求解时，必须给出迭代的收 敛标准，否则无法终止迭代计算。迭代收敛准则有位移准则、失衡力准则和能量准则 三种。通过强度折减，使边坡达到破坏状态时，滑动面上的位移将产生突变，产生很 大的无限值的塑性流动，有限元程序无法从有限元方程组中找到一个既能满足静力平 衡，又能满足应力一应变关系和强度准则的解，所以此时，不管是从力的收敛标准还 是位移收敛标准，计算都不能收敛。 ( 2 ) 根据计算域内某一部位的位移或最大位移与折减系数之间关系曲线的变化特 征确定失稳状态。边坡的变形破坏总具有一定的位移特性，因此有限元计算的位移结 果是边坡失稳最直观的表达。如当折减系数增大到某一特定值时，某一部位位移突然 迅速增大，则认为边坡发生失稳。 ( 3 ) 通过域内广义剪应变等某些物理量的变化和分布来判断，如当域内的塑性区 ( 或某一幅值的广义剪应变) 连通时，则判断边坡发生破坏。理论上，边坡的变形过 程总是伴随着一些物理量的出现和发展，如塑性区、塑性应变、广义剪应变和应力水 平等，当这些物理量达到一定的值时，边坡失稳。 6 r i f f i t h s 和l a n e ( 1 9 9 9 ) 用强度折减法研究了简单边坡、软弱地基路堤、软弱 夹层边坡和土坝边坡的稳定性，以迭代次数超过1 0 0 0 次仍未收敛为失稳判据，计算的 安全系数和极限平衡法十分接近。d s w s o n 和r o t h n 3 1 ( 1 9 9 9 ) 假定当结点的失衡力与外 荷载的比值超过1 0 。3 来确定安全系数，其中也隐含着必须以某一迭代次数作为收敛准 则。郑颖人，赵尚毅等人口4 4 5 1 认为土体滑动面塑性区贯通是土体破坏的必要条件，但 不是充分条件，还要看是否产生很大的且无限发展的塑性变形和位移，有限元计算中 表现为塑性应变和位移产生突变，在突变前计算收敛，突变之后计算不收敛，表征滑 面上土体无限流动，因此可把有限元静力平衡方程是否有解、有限元计算是否收敛作 为边坡破坏的依据。 6 第一章绪论 z i e n k i e w i c e 最初提出有限元强度折减方法所采用的失稳判据就是最大结点位 移。t a n 和d o n a l d ( 1 9 8 5 ) 用结点的位移和折减系数的关系曲线处理为两段直线，用 两直线的交点对应折减系数作为安全系数；宋二祥( 1 9 9 7 ) 采用坡项位移折减系数关 系曲线的水平段作为失稳判据。m a t s u i 和s a n n 印( 1 9 9 2 ) 在模拟填土时，以剪应变超 过1 5 作为边坡失稳的判据；而模拟开挖时以剪应变增量作为失稳判据。g r i f f i t h sd v 2 4 和张培文、陈祖煜h 7 1 认为当无量纲位移e 瓯觚( y h 2 ) 有一个突然增加时，边坡处 于临界平衡状态，其中e 为材料的弹性模量，瓦觚为计算结果有限元网格节点的最大 位移，y 为材料的比重，日为边坡的高度。 连镇营n 踟( 2 0 0 1 ) 基于强度折减弹塑性有限元分析结果，绘制边坡内广义剪应变 分布，并认为若某一幅值广义剪应变的区域在边坡中相互贯通，则意味边坡已经失稳 破坏。栾茂田h 町( 2 0 0 3 ) 将强度折减概念与弹塑性有限元数值分析及结果的计算机实 时显示技术相结合，并建议采用塑性应变作为失稳判断标准，由计算得到的塑性应变 及其分布，认为塑性开展区相互贯通时边坡失稳。 应该说三种失稳判别标准的应用都有其适用范围，也在不同工程的应用中呈现出 其可行性和合理性。既然在实际应用中它们的判定结果都是合理可用的，那么它们相 互之间必然存在某种联系或是在本质上存在一致性。所以，本文将在前人的工作基础 上，对三种判别标准的适用范围以及它们之间的联系和差异进行探讨研究。 1 3 本文的主要工作 有限元边坡稳定性分析方法的应用越来越广泛，对其研究工作也开展了不少，取 得了不少成果。如何判定边坡何时处于破坏状态是对该方法研究比较热的一个问题， 本文将在前人研究成果的基础上，针对三种失稳判别标准在实际应用中的问题进行研 究分析。主要内容有： 1 对目前有限元边坡稳定分析方法做一个概述，重点分析前人在强度折减法中判 别标准问题上的研究成果。 2 非线性方程的求解是有限元边坡稳定分析计算的一个重要内容，文中将对求解 河海人学硕： ：论文 非线性方程的迭代法、增量法及混合法做简单介绍，对各自的基本理论及求解过程做 一详细叙述。并对迭代计算中采用的三种收敛准则及其各自的适用条件做简单论述。 3 阐述有限元计算中广泛采用的、适用于岩土类材料的弹塑性本构模型，包括屈 服准则、加卸载准则、塑性势理论以及加载硬化规律等原理。 4 论述强度折减法的三种失稳判别标准，并对前人的研究成果进行总结。应用 a d i n a 有限元软件对三个边坡算例进行有限元强度折减计算，分别采用三种失稳判别标 准进行稳定性分析，对三种失稳判别标准适用性、三者之间的一致性、各自的适用范 围进行分析研究。 5 在前面的分析研究基础上，对一实际边坡工程进行有限元计算，对其应力、变 形状态进行分析，并利用三种失稳判别标准对其稳定性进行判定。 8 第二章边坡稳定性分析的非线性有限元理论 第二章边坡稳定性分析的非线性有限元理论 2 1 引言 在有限元线性问题中假定1 ) 材料的应力一应变关系是线性的，即假定材料符合 虎克定律；2 ) 应变一位移关系是线性的，即小位移假定。在实际工程中，有不少问题 是符合上述假定的，因此线性有限元得到了广泛的应用，但是绝大多数问题是不符合 上述假定的，需进行非线性分析咖 5 引。 图2 - 1 给出了实际非线性分析中应用非常方便的一种分类瞄引。根据产生非线性的 原因，非线性问题主要有材料非线性、几何非线性和状态变化非线性三种畸4 5 引。 r r t l l 卜l 叫 ( 2 ) ( 3 ) 9 g = l ，b 0 0 4 河海大学硕i = 论文 ( 4 ) ( 5 ) 图2 - 1 非线性分析的分类 ( 1 ) 线弹性( 无穷小位移)( 2 ) 仅材料非线性( 无穷小位移，但应力一位移关系是非线性的) ( 3 ) 火位移和大转角但小应变( 线性或非线性的材料特性)( 4 ) 大位移、人转角和人应变 ( 5 ) 在位移处边界条件的变化 1 材料非线性 材料的应力应变关系为非线性时即需进行这类分析，该关系通常与加载历史有关， 加载和卸载不是同一途径。但材料非线性问题仍属于小变形问题，位移和应变是微量， 其几何方程是线性的。土、岩石、混凝土等具有典型的材料非线性性质，所以土坝、 岩土地基的稳定性和加固，地下洞室和边坡的稳定性都应当按材料非线性问题处理。 非线性材料包括非线性弹性、超弹性、热弹性、弹塑性、塑性、粘弹塑性及蠕变等材 料。 2 几何非线性 几何非线性属于大变形问题，位移和应变或者它们中的一个是有限量。可能会有 三种情况：大位移( 包括线位移和角位移) 、小应变，小位移、大应变和大位移、大应 变。几何非线性分析研究结构在荷载作用下几何模型发生的改变，如何改变及几何改 变的大小。所有这些均取决于结构受载时的刚性或柔性，且此时应变、位移的关系是 非线性的。工程中的实体结构和板壳结构都存在几何非线性问题，如弹性薄壳的大挠 度分析，压杆或板壳受载弹性屈曲后的稳定性问题。 第二章边坡稳定性分析的非线性有限元理论 3 状态变化非线性 此类非线性是考虑物体在运动时边界条件有变化的这类问题的分析，此外还包括 生死单元的定义。由于接触体的变形和接触边界的摩擦作用，使得部分边界条件随加 载过程而变化，且不可恢复。这种由边界条件的可变性和不可逆性产生的非线性问题， 称为接触非线性。工程中，有许多接触非线性的问题，如混凝土坝纵缝和横缝缝面的 接触，面板堆石坝中钢筋混凝土面板与垫层之间的接触，岩体节理面或裂隙面的工作 状态等。 2 2 非线性有限元的求解方法 一般说来，非线性分析的基本问题是找出物体对应于所加荷载下的平衡状态，假 设外加荷载为时问的函数，则物体的有限元平衡条件可表示为 r 一f = 0( 2 1 ) 其中向量r 表示在时问t 位形中外加结点力，向量f 表示对应于这个位形的单元应力 的结点力。考虑非线性响应的解时应注意到，在整个加载过程中( 即时间变量t 可以从 0 到感兴趣的最大时间中取任何值) 必须满足平衡关系式( 2 一1 ) 。在除加载等级的定 义外无其它时间效应的静态分析中，时间仅是便于表示荷载作用的不同强度及其相应 位形的一个度量。但在动态分析中和有材料时间效应的静态分析中，时间变量是严格 地包括在建立真实物理状况模型之中的一个实际变量。因此，使用时间变量来描述荷 载作用和解的历程，是一个非常普遍的方法，且符合“动态分析本质上是包括惯性效 应的静态分析”瞄3 | 。 在许多非线性分析中仅要求在特定加载等级或特定时间下达到的应力和位移，一 些非线性静态分析中可以算出对应于那些加载等级的平衡位形而无需求解其它的平衡 位形。然而，当分析含有与途径有关的非线性的几何或材料条件，或随时间变化的现 象时，则必须在全部考虑感兴趣的时间范围内解出平衡关系式( 2 1 ) 。这个响应计算 可有效地用逐步增量分析作出。如果在与时间无关的静态解中，总荷载是同时作用的， 并仅需计算对应于该荷载的位形，这个分析就简化为一步分析。然而在实际上为了计 算上的原因，甚至这种分析情况也常常需要增量解，经过很多荷载步长最终达到总的 一1 1 一 河海人学硕 ：论文 作用荷载。 求解非线性问题的方法主要有迭代法、增量法和混合法。下面对各种求解方法做 一介绍。 2 2 1 迭代法，5 7 1 用迭代法求解非线性问题时，一次施加全部荷载，然后逐步调整位移，使基本方 程得到满足，具体的计算方法主要有牛顿一拉菲逊迭代法( n e w t o n - - r a p h s o n 迭代法) 、 拟牛顿法。 2 2 1 1 牛顿一拉菲逊迭代法 牛顿一拉菲逊法是求解非线性方程组 、i ，( 6 ) 三f ( 6 ) 一r ( 2 2 ) 的一个著名方程，简称牛顿法。设v ( 6 ) 为具有一阶导数的连续函数，6 = 6 是方程 ( 2 - 2 ) 第i 次迭代的近似解。若 、l ，1 = v ( 6 。) 三f ( a 1 ) 一r 0 ( 2 3 ) 希望找到一个更好的、方程( 2 2 ) 的近似解为 6 = 6 件1 ：6 + 6 ( 2 - 4 ) 将式( 2 3 ) 带入式( 2 2 ) ，并在6 = 6 附近按一阶泰勒级数展开，则、l ，( 6 ) 在6 处的 线性近似公式为 、- ，川叫+ 尉 5 ) 苴由 、i = 6 甜 ( 2 6 ) 第二章边坡稳定性分析的非线性有限元理论 f ，却、一 l 面厂 a a 匹 a a t a a 瓯 、i ，。v 2 、l ，。 l 、i ，lv 2 、l ，。i ( 2 - 7 ) 在非线性有限兀中，有限兀的平衡要求相当于求f 述方程的解 、i ，( 6 ) = 什出r 件血f ( 6 ) = 0 ( 2 8 ) 式中，h 垃r 为f + f 时刻作用于结点的外荷载向量，件出f ( 6 ) 为t + a t 时刻等价于单 元应力的结点力向量。式( 2 7 ) 中的差值即为结点失衡力。假定在迭代求解中已经算 得什血6 ( n ，而什出6 + 为真实解，i = 0 ，1 ，2 ，n 为迭代步数。则由前述的一阶泰勒 展开式得到 、i ，( t + a t 6 , h ( f + 吲f ) ) + 舢6 i ( f + f 6 + 一舢6 ) ( 2 _ 9 ) 把式( 2 - 8 ) 带入式( 2 - 9 ) 得到 仙6 r 呵一舢6 户地他一) 其中假定作用外荷载与位移无关。定义 a f t ( = h 6 + 一+ 出6 ( )( 2 11 ) 并认为 k = 舢戤) 式中什出k ( 是第i 次迭代的切向刚度矩阵。因此( 2 1 0 ) 可写为 + 出k ( ) 6 ( “1 ) = + 出r 一+ f ( )( 2 1 3 ) 由上式所得位移增量修正值a 6 ( 件，并得到下一个位移近似值 + 血6 ( 洲) t + a t6 ( + a 6 ( 洲 ( 2 1 4 ) 河海大学硕l ：论文 这样，式( 2 一- 1 3 ) 、( 2 1 4 ) 就构成了式( 2 - 8 ) 的牛顿一拉菲逊解。重复上诉迭代过程 直到满足适当的收敛准则为止。牛顿法的迭代公式就可归纳为式( 2 - 1 2 ) ( 2 - 1 4 ) 。 牛顿法每次迭代都需重新计算切向刚度矩阵，因此计算量相当大。迭代中的切向 刚度矩阵均采用初始的切线刚度矩阵o k ，式( 2 - 1 3 ) 就变为 o k 6 ( + 1 ) ：+ 垃r 一+ 血f ( )( 2 1 5 ) 从而避免了多次重新计算和分解( 2 1 2 ) 的系数矩阵。这就是修正的牛顿迭代法的一 种“初始应力法”，而前述的牛顿迭代法称为完全牛顿法。这种初始应力法相当于 有限元系统初始形状响应的线性化，其结果会导致非常缓慢的收敛甚至解的发散。 在修正的牛顿法中有一个介于完全牛顿法和上述初始应力法之问的方法。在这个 方法中用 7 k - 6 ( + 1 ) = + r 一+ 出f ( )( 2 1 6 ) 和初始条件什出f ( o ) = f ，+ 出6 ( o ) = 6 ；丁看作是对应于时刻0 ，a t ，2 a t ，或f 得 到的某个平衡形状。修正的牛顿法对重新形成刚度矩阵的需要比完全的牛顿迭代法要 少，而只根据得到的平衡形状不断地修改刚度矩阵。当必须修改刚度时，时间步长的 选择与系统响应的非线性程度有关，即响应越是非线性就越要经常进行修改。在对系 统的性质事先没有任何了解的情况下，可能在每一时间步长的开始时修改刚度矩阵最 为有效，此时, - = t 。 2 2 1 2 拟牛顿法 作为另一种牛顿迭代形式，是一类为非线性方程组的迭代求解而发展起来的矩阵 修改法，用该方法修改系数矩阵( 或它的逆) 以提供一个从第f 次到z + 1 次迭代的矩 阵的割线逼近。定义位移增量 6 ( 。十1 ) = + 垃5 汁1 一h 出6 ( i )( 2 1 7 ) 及不平衡荷载r n 、不平衡荷载增量丫 a r c ) = h 垃r 一+ 出f ( i )( 2 1 8 ) 丫( “1 ) = a i r ( n a r ( + 1 ) ( 2 1 9 ) 修改的矩阵件k ( 州应满足拟牛顿方程 一 1 4 第二章边坡稳定性分析的非线性有限元理论 。+ 出k ( 件1 ) 6 ( “) = 丫( 1 “) ( 22 0 ) 拟牛顿法提供了一个介于在完全牛顿法中重新形成刚度矩阵与在修正牛顿法中采用前 面形状的刚度矩阵之间的一种折衷办法。在可用的拟牛顿法中，b f g s ( b r o y d e n - - f l e t c h e r - - g o l d f a r b - - s h a n n o ) 方法似乎最有效，并由m a t t h i e s 和s t r a n g 最早建议 用于有限元分析中。 在b f g s 法中，用下述方法对第j 次迭代讨算”出6 ( “1 1 和7 + 出k ( 川) ，其中 “出k ( o ) = k 。 第一步：计算位移向量增量 a 6 = ( + 出k 1 ) 。( t + a t r 一+ f ) ( 2 2 1 ) 对个实际位移增量，这个位移向量确定了一个“方向”。 第二步：沿方向6 作线性搜索以满足该方向的“平衡”。在此线性搜速中我们计 算位移向量 + 6 。6 7 “= + 6 ( ) + 臌6( 22 2 ) 其中是一个标量乘子，并且计算相应于这些位移的不平衡荷载( “r 一”出f ) 。 改变参数直到在6 方向上按内积( 6 ) 7 ( ”出r 一”“f ( o ) 所定义的不平衡荷载分量 很小为止。若收敛限值为s t o l ，当下述方程满足时， a 8 7 ( ”出r 一+ “f ) _ s t o l a 8 7 ( t + a t r 一7 + “f ) ( 2 2 3 ) 该条件也就满足了。 使式( 2 2 3 ) 满足的口的终值决定了式( 2 2 2 ) 中的”6 “。从而可用式( 21 7 ) 和式( 2 1 9 ) 算出a 5 川和1 r ，并进行修改矩阵的计算以满足式( 22 0 ) 。 第三步：计算系数矩阵的修改值。在b f ( ；s 方法中修改矩阵可表示为积的形式： ( 舢k ) = ( a ) 7 ( 舢k ) 一1 a 川 ( 22 4 ) 式中矩阵a ( ”1 是简单形式的( n x n ) 矩阵 a = i + v 川( w 州) 7 ( 22 5 ) 从已知的结点力和位移用下述公式计算向量v ”和w 川) ： 河海大学硕i 二论义 v“+，：一ii予呈；掌)aft(1 2 一m k “，”n 一丫“， 。 i ( 6 “+ 1 ) “k “”j 。 w ”吐蔷铬。i ( 2 2 6 ) ( 2 2 7 ) 式( 2 2 6 ) 中的向量t + a t kc o a 6 “1 等于少( i + a i r 一“f ) 并且已经计算出来。 由于式( 22 4 ) 确定的积是正定和对称的，为避免出现数值上的危险，要计算修 改矩阵a ( 川的条件数c 州) ： c “+ ”_ 茄a 6 ( m ) 篙a 8 i 2 协z s ， 。 i () “k “+ 1 l 、。 用这个条件数和某个预定的限值( 例如1 0 i ) 相比较，如果条件数超过了这一限值，则 中止修改。 从以上迭代法的计算过程可以看到，随着迭代次数的增加，不平衡力逐渐减小， 并趋于平衡。可见迭代法就是用总荷载作用下不平衡的线性解去逼近平衡的非线性解， 迭代的过程就是消除不平衡力的过程。对于不同的迭代方法，这一过程的快慢也就是 收敛速度是不同的。一般来说，牛顿法最快，拟牛顿法次之，修正的牛顿法最慢。不 过，各种方法的效率不仅与收敛速度有关，还与每一步迭代所需的计算量有关。牛顿 法每一步的计算量最大，拟牛顿法次之，修正牛顿法最小。在本文的计算中，将采用 完全牛顿法进行迭代计算。 2 2 2 增量法”，“ 采用增量法分析非线性问题时，把荷载分为许多荷载增量，这些增量可以相等， 也可以不等。每次旌加一个荷载增量。在每一步计算中，假定方程是线性的，刚度矩 阵k 是常数，在不同的荷载增量中，刚度矩阵可以具有不同的数值。增量法是用一系 列线性问题去近似非线性问题，实际上是用分段线性的折线去代替非线性曲线。 把荷载分为m 个增量，f = 1 ，2 ，m 。施加第i 个荷载增量a r ；，得到一个 位移增量a 6 。，累积后得到位移6 。则总荷载为 r = a r ， ( 2 2 9 ) 第二章边坡稳定性分析的非线性有限元理论 在施加第z 个荷载增量后，荷载为 f r f = r ( 2 3 0 ) _ ，= l 每一个荷载增量产生一个位移增量6 ，和应力增量6 ，因此在施加第j 个荷载增量 后，位移和应力分别为 f 5 ；= a 6 ， ( 2 _ 3 1 ) _ ，= l f 6 f = 6 ， ( 2 3 2 ) _ ，= l 在有限元中，常把荷载定义为时间的函数，荷载增量也就对应于时间上的荷载增 量。如果所求问题为与时间无关的静力问题，荷载步即为时间步，即 a r f = r f r f l = + 血r 一r ( 2 3 3 ) 如何由荷载增量计算位移增量和应力增量，主要有欧拉法、修正的欧拉法。 2 2 2 1 欧拉法 设t 时刻( 即加了f 一1 个荷载增量后) 的位移6 和应力6 已求出，根据应力应变 关系就可以确定该时刻的刚度矩阵k 。假定在f 时刻至f + 垃时刻内刚度矩阵保持不 变，并近似等于k 。于是t 时刻到t + a t 时刻( 即加第f 个荷载增量) 的位移增量h & 6 可由下式计算得到： k a t + a t 6 = + 舡r 一r ( 2 3 4 ) 根据位移增量件出6 可进一步求出应力增量a t + a t 6 ，则f + a t 时刻的位移和应力分别 为 + 垃6 = 6 + + 出6( 2 3 5 ) + 出6 = 6 + a t + a t 6 ( 2 3 6 ) 由该方法的求解过程图2 2 可以看出，每步计算都会引起偏差( 不平衡力) ，使折 线偏离曲线，解答产生漂移，随着求解步数的增加，由于偏差的累积使最后的解答离 开真解较远，从而降低了计算精度，因此需对其作些改进。 河海人学硕i j 论文 f r h r r 。 0 2 2 2 2 修正的欧拉法 616m6m + l 图2 2 欧拉法求解过程 修正的欧拉法就是首先施加荷载增量的一部分，即秒( 件垃r 一r ) ，0 0 1 ， 用t 时刻的刚度矩阵k 计算t 时刻到t +