（水工结构工程专业论文）土石坝动力稳定的强度折减法研究.pdf

大连理工大学硕士学位论文 摘要 随着计算机的飞快发展，有限元强度折减法作为一种边坡稳定分析方法逐渐受到越 来越多的工程人员的关注。由于其概念明确，物理模型直观，工程参数易于获取，不用 对滑裂面进行任阿假设，不受边界条件限制，而且计算结果和传统的极限平衡法有着惊 人的相似。因此，有限元强度折减法开始逐步被研究和应用，但是，如何用有限元强度 折减法分析边坡在地震动力作用下的稳定性一直没有进一步深入研究。 本文在用有限元强度折减法分析边坡在自重作用下稳定性的基础上，提出了一种用 有限元强度折减法分析高边坡在地震动力作用下稳定性的方法。首先通过三维非线性有 限元静力分析得到动力分析的初始应力，再进行动力分析得出边坡在动力作用下每一瞬 时的加速度分布，根据每一瞬时的加速度分布的结果，将地震作用力作为惯性力加在结 构上，用有限元强度折减法进行边坡稳定分析，得出在地震动力作用下每一瞬时的安全 系数及对应滑裂面，拟合可得到一条安全系数随时问变化的曲线，根据这条安全系数时 程曲线及滑裂面，综合得出结构的抗震稳定安全系数及对应的滑裂面。依据这些原理， 编制了有限元分析各个时刻稳定性的程序，通过对验证算例的分析，证明该方法分析地 震动力作用下分析边坡稳定性的合理性。 最后，采用有限元强度折减法对糯扎渡粘土心墙堆石坝进行了动、静力稳定分析， 得到各种工况下的抗滑稳定安全系数和临界滑裂面，通过与极限平衡法、有限元搜索法 进行对比分析，论证了有限元强度折减法在实践应用中的合理性和适用性，得出了糯扎 渡直心墙堆石坝边坡设计是稳定的结论。 关键词：边坡；稳定分析：有限元；强度折减；极限平衡 张彩双：土石坝动力稳定的强度折减法研究 r e s e a r c ho ns t r e n g t hr e d u c t i o nt e c h n i q u ef o re a c h r o c k f i l l d a m d y n a m i cs t a b i l i t ya n a l y s i s a b s t r a c t w i t ht h eh i 曲d e v e l o p m e n to fc o m p u t e rt e c h n o l o g y ，a sas l o p es t a b i l i t ym e m o d ，s t r e n g t h r e d u c t i o na t t r a c t sm o r ea n dm o r ee n g i n e e r s a t t e n t i o n ，t h i sm e m o dh a sm a n ya d v a n t a g e s s u c ha si th a sae x p l i c i tc o n c e p t i o n ，av i s u a lp h y s i c a lm o d e l ，p a r a m e t e r st h a ta r ee a s i l y a c q u i r e d ，n e e d l e s ss h a p eo fs l o p ec r i t i c a la s s u m p t i o na n d i ti sf i t t e dt oc o m p l i c a t e db o u n d a r y m o s to fa l li tc a l lg e tas a f e t yf a c t o rt h a ti sh i g h l yc l o s et ot h es a f e t yf a c t o ro ft r a d i t i o n a l l i m i te q u i l i b r i u mm e t h o d ，s ot h em e t h o do fs t r e n g t hr e d u c t i o ni sb e i n gr e s e a r c h e da n da p p l i e d s t e pb ys t e p h o w e v e r ，t h ea p p l i c a t i o no ft h i sm e t h o d t os l o p es t a b i l i z a t i o nw i t ht h ee f f e c to f d y n a m i c a le a r t h q u a k es t i l lh a s n tb e e nw o r k e do v e rf u r t h e r o nt h eb a s i so fa n a l y z i n gs l o p es t a b i l i z a t i o nu n d e rt h ea c t i o no fg r a v i t y ，t h i sp a p e r a p p l i e ss t r e n g t hr e d u c t i o nm e t h o dt oa n a l y z i n gs l o p es t a b i l i z a t i o nu n d e rt h ee a r t h q u a k e f i r s t ， t h r o u g h3 - dn o n l i n e a rf i n i t ee l e m e n ta n a l y s i s ，i n i t i a ls t r e s s e sa r eg o t t h e nt h r o u g hd y n a m i c a l a n a l y s i s ，e v e r yi n s t a n t a n e o u sa c c e l e r a t i o nd i s t r i b u t i n go fs l o p eu n d e rd y n a m i c a le a r t h q u a k e c a nb eg o tt o o a c c o r d i n gt ot h i sr e s u l t ，a d d i n ge a r t h q u a k ef o r c et o s t r u c t u r eb yw a yo f i n e r t i a lf o r c e ，t h i sp a p e ra n a l y s e st h es l o p es t a b i l i z a t i o nb yu s i n gf e ms t r e n g t hr e d u c t i o n m e t h o da n dg a i n se v e r yi n s t a n t a n e o u ss a f e t yf a c t o ra n di t sc o r r e s p o n d i n gs l i d e c u r v e so f s a f e t yf a c t o r sv e r s u sd u r a t i o nc a nb em a d eb yc u r v e f i t t i n g b ys a f e t yf a c t o r sv e r s t l sd u r a t i o n c u r v e sa n ds l i d e s ，w ec a nv a l u es t r u c t u r e sq u a k e p r o o fs t a b i l i z a t i o na n dp o i n ti t ss l i d eo u t b y t h e s ep r i n c i p l e s ，p r o g r a m so ff e ms t a b i l i z a t i o na n a l y s i sf o re a c hm o m e n ta r eo r g a n i z e d s o m ee x a m p l e sa n a l y z e db yt h i sm e t h o dv a l i d a t et h er a t i o n a l i t yo fa p p l y i n gf e ms t r e n g t h r e d u c t i o nt os l o p es t a b i l i z a t i o na n a l y s i su n d e rt h ea c t i o no fd y n a m i c a le a r t h q u a k e a tl a s t ，t h i sp a p e ru s e st h i sm e t h o dt oa n a l y z et h es t a t i ca n dd y n a m i c a ls t a b i l i z a t i o no f n u o z h a d ud a m as a f e t yf a c t o ra n dac r i t i c a ls l i d eo fe a c hc a s ea r ee d u c e d c o m p a r e dw i t h t h er e s u l to ft h el i m i te q u i l i b r i u mm e t h o da n df e ms e a r c h i n gm e t h o d ，c o n c l u s i o n sc a l lb e m a d et h a tf e m s t r e n g t hr e d u c t i o nm e t h o di sr e a s o n a b l ea n da p p l i c a b l ei np r a c t i c ea n dd e s i g n o f n u o z h a d ur u b b l ed a r nj sr e a s o n a b 】e k e yw o r d s ：s l o p e ；s t a b i l i t ya n a l y s i s ；f i n i t ee l e m e n tm e t h o d ；s t r e n g t hr e d u c t i o n l i m i te q u i l i b r i u mm e t h o d 独创性说明 作者郑重声明：本硕士学位论文是我个人在导师指导下进行的研究工 作及取得研究成果。尽我所知，除了文中特别加以标注和致谢的地方外， 论文中不包含其他人已经发表或撰写的研究成果，也不包含为获得大连理 工大学或者其他单位的学位或证书所使用过的材料。与我一。同工作的同志 对本研究所做的贡献均已在论文中做了明确的说明并表示了谢意。 作者签名：日期： 大连理工大学硕士研究生学位论文 大连理工大学学位论文版权使用授权书 本学位论文作者及指导教师完全了解“大连理工大学硕士、博士学位论文版权使用 规定”，同意大连理工大学保留并向国家有关部门或机构送交学位论文的复印件和电子 版，允许论文被查阅和借阅。本人授权大连理工大学可以将本学位论文的全部或部分内 容编入有关数据库进行检索，也可采用影印、缩印或扫描等复制手段保存和汇编学位论 文。 作者签名 导师签名 张彤双 巡年l 月上竺日 大连理t 大。学硕士学位论文 1 绪论 1 1 课题研究的背景和意义 地震作用诱发的边坡滑动和坍塌是常见的灾害之一，特别是在山区和丘陵地带，地 震诱发的滑坡往往分布广、数量大、危害大。崩塌滑坡因其巨大的致灾力而广泛引起人 们的关注，仅2 0 世纪，就已经造成了数万人丧生和几十亿美元的损失。地震时，地震 诱发的滑坡崩塌灾害，特别是在山区地区，其引起的次生危害比地震时造成的直接危害 还要大。 表1 】国内外由地震诱发的重大滑坡灾害 t a b 1 1d i s a s t r o u s l a n d s l i d ea c c i d e n t se a r t h q u a k e i n d u c e da t h o m ea n da b r o a d 我国是世界上多地震国家之一，地震活动分布范围较广，基本烈度在6 度以上的地 区占全国总面积的6 0 以上。在水库工程的挡水建筑物中，土石坝是最常用的坝型，有 很多土石坝建于强烈地震区。在这些地区，土石坝由于地震作用而部分或完全破坏的事 情时有发生【2 】，主要破坏形式为坝顶及坝坡的纵横裂缝，坝体滑动滑裂、凸起和失稳、 坝基失稳液化、坝顶坝坡垂直沉陷、水平位移、渗流量增加、坝内涵洞或埋管、坝顶或 坝坡上竖井破坏，都对大坝正常运行产生不利影响，甚至威胁大坝正常运行。在我国虽 然尚无地震溃坝事例，但也有一些土石坝因地震而发生滑坡、震陷、裂缝等震害，如北 京密云白河土坝、唐山陡河土坝、辽宁石门土坝等。因此在工程设计中研究地震对土石 坝正常运行可能产生的影响，研究土石坝在地震作用下的破坏形式和破坏机理，对坝基 坝体的材料进行动力试验，对坝基坝体系统进行动力分析，并采取相应的工程措施非常 必要。 张彩双：十石坝动力稳定的强度折减法研究 以往设计中，大都采用拟静力法校核坝体的稳定，随着土的动力分析理论的逐步完 善和计算机技术的不断发展，动力分析方法正逐渐地运用于土石坝抗震设计中我国碾 压式土石坝设计规范( s dj 2 1 8 8 4 ) 中规定，对高烈度地震区的i 、i i 坝，除进行 静力的应力和变形计算外，还应进行动力分析。在选择坝体的材料时，坝体内不应采用 可能液化的土料，坝基中如有可能液化的土层，应对该土层进行核算，检验在设计地震 烈度下是否可能发生液化。水工建筑物抗震设计规范( s d l l 0 - - 7 8 ) 也规定对1 5 0 m 高以上的高坝应进行动力分析，并提出了地震荷载的动力分析方法。 因此，了解有关土石坝动力稳定分析方面的基本理论，掌握动力稳定分析的基本方 法，并能在工程设计中实际应用，具有非常重要的现实意义。 1 2 问题的研究现状 土石坝地震稳定性分析方法是随着科技的进步，特别是计算机技术的快速发展， 以及对大量土石坝震害现象进行调查研究的不断深入而发展起来的。 长期以来，许多国家都把用拟静力法对土石坝进行稳定性分析纳入到规范之中p 】。 我国s d j l 0 - 7 8 中规定对土石坝采用拟静力法进行抗震稳定性计算， 日本大坝委员会 1 9 7 8 年发布坝工设计规范、日本建设省河川开发科1 9 9 1 年颁发土石坝抗震设计 指南，其中土石坝的抗震设计与我国的s d j l 0 7 8 类似。随着土石坝震害现象的不断 出现，目前在美国土石坝抗震计算主要采用动力有限元方法，内容包括建立在有限元动 力分析基础上的滑动稳定计算和变形计算。1 9 7 6 年我国唐山地震中，密云水库白河主 坝因保护层液化而引起的滑坡震害后，我国也对建在高烈度区的土石坝进行了动力分 析，用来研究使用拟静力法无法得出的坝体和坝基内动应力分布及地震引起的坝体变 形。 目前常用的方法有规范推荐的拟静力法( p s e u d o s t a t i ca n a l y s i s ) 、n e w m a r k 滑块分 析法( n e w m a r ks l i d i n gb l o c ka n a l y s i s ) 及动力有限元法等。现在分别介绍如下： 121 拟静力法 在过去6 0 年或更长的一段时间内，标准的评价在地震过程中抗滑稳定安全性方法 是被称为拟静力的评价分析方法。在这种分析方法中，地震对有可能发生滑动的作用被 描绘成等效的静的水平力，其大小由地震系数k 或聍。、质量块的重量g 确定。 如图1 1 所示，具有等效加速度的一次地震产生一水平方向的集中力，单位土重量 上所受力的大小为n ，集中力的合成为胛。，像重力一样穿过滑块a b e 的重心d l 。这 人连理丁大学硕士学位论丈 一集中力的水平力臂为f ，产生使滑块发生转动的力矩的大小为? l g f 。因此，地震产 生的安全因子g 可以表示成： 瓯= 丽i s l r 而 n 。数值的大小取决于地震强度：严重地震，烈度为度，n 。= o 1 0 ；强烈地震、破 坏性地震，烈度为x 度，以。= 0 2 5 ；灾难性地震以。= o 5 0 。1 9 0 0 年s a nf r a n c i s c o 地震 对应的地震系数为月。= 0 2 5 。 方程( 1 1 ) 是建立这样一个简单的假设基础之上，即水平加速度n 。一直作用在潜 在的滑块上，且作用方向不变。事实上，这种作用无论是在方向还是数值大小上都是时 刻变化的，因此这种方法所表达的地震作用的概念是不准确的。理论上讲，当g = 1 0 时会发生滑移，实际上当g s 1 0 却发生严重 震害。这取决于坝体土料的特性。 一7 相当长时问内，这种地震系数分析方法被广泛的采用，并为土石坝的抗震设计提供 了指导。各个国家对于地震系数的取值也做了规定。在美国，的范围是0 1 0 5 0 1 1 5 ； 在日本其值也小于o 2 0 。 在不考虑坝坡几何特性和材料特性的情况下，如果所计算的安全因子大于1 ，我们 通常认为坝体地震稳定性问题得到解决。 张彩双：十石坝动力稳定的强度折减法研究 用拟静力进行设计的填筑式土石坝，虽经长期的运行经验证明，在总体上能经受地 震而绝少失事。但这种方法既没考虑地震的特性，诸如振动频率、次数和地震持续时间 等，又没考虑坝身材料的动力性质，如变形性质和阻尼性质等，因而无法反映坝在地震 时的反应特性。另外，建立在刚塑性模型和极限平衡理论基础上的圆弧滑动法所给出 的稳定安全系数，并不能说明坝体抗震的安全储备实际情况。计算的结果表明，与动力 分析法相比，结构上部的地震荷载偏小，而底部则偏大【4 j 。因此，自7 0 年代以来，土石 坝的各种动力分析理论、方法和坝体抗滑安全性评估方法，有许多进展，并逐步为工程 界所应用。 1 22n e w m a r k 滑块分析法 1 9 6 5 年的r a n k i n e le c t u r e 上，n e w m a r k 首次提出评价由于地震作用引起的坝体有 限滑移的基本原理 5 a 在这次重要的发展过程中，他设想如果作用在潜在滑移质量块上 的惯性力在一定程度上超过屈服抵抗力，滑移震害就会出现。滑移运动也就开始，而当 惯性力反向时运动则停止。通过计算当惯性力足够大从而使屈服现象出现时的加速度 ( 如图1 3 ) ，同时将滑块超过屈服加速度的有效加速度作为时间的函数，滑块的速度 和位移就可以计算出来，如图1 4 。 有限滑移分析方法的一股计算步骤： ( 1 ) 求屈服加速度a 1 ( 2 ) 求有效加速度；( 3 ) 求滑块的永久位移：对某一预 测滑动土体，当地震引起的有效加速度吼( f 1 超过所求得的屈服加速度的时候，就认为 滑动体有位移产生，位移大小可由加速度差值( q d ，1 两次积分得到。 图1 2 滑块受力示意图 f i g 1 2s k e t c ho f a c t i n gf o r c e 皆 一 。 一 w ( ) 一 ， k 一 大连理二r 大学硕士学缱论文 划 煳 划 心 渣 趔 t 图1 3 用有效加速度时程确定速度和位移的示意图 f i g 1 3s k e t c ho fe f f e c t i v ea c c e l e r a t i n gt i m eh i s t o r yt od e t e r m i n ev e l o c i t ya n dd i s p l a c e m e n t n e w m a r k 有限滑移分析方法不仅为客观评价坝堤在地震中表现方面推进了一步，而 且在应用方面建议了一种分析方法。对于特定的土类( 在地震中强度不发生明显降低的 土) ，如压密粘性土、紧密的饱和砂和砂砾石、非饱和土等，这一分析方法可以给出合 理的评价结果。这种分析方法的缺点是屈服加速度不易确定，如何确定屈服加速度还需 要进一步研究。 1 2 3 动力有限元法 近2 0 年来，随着计算技术的迅速发展，动力分析方法也发展很快，从最早的线性 总应力法，逐步发展为基于非线性有限单元法基础上的有效应力动力分析方法和采用复 杂弹塑性模型并考虑水土耦合作用的动力分析方法，从只能分析一维问题发展到能够分 析二维、三维问题，从只能分析饱和土体发展到能够分析多相非饱和土体。本构模型也 从早期的线弹性模型发展为粘弹性模型、弹塑性模型、边界面模型、内时模型和结构性 模型等。动力方程的计算方法也发展很快，出现了逐步积分法、振型叠加法、复反应分 析法、有限差分法、子结构法等方法。目前常用的动力分析方法主要有剪切楔法、集中 质量法和数值分析法( 包括有限单元法、有限差分法和边界元法) 等，其中前两种方法还 可区分为总应力法和有效应力法。 张影双：十石坝动力稳定的强度折减法研究 有限单元法可计算二维问题和三维问题，可以按坝的分区考虑不同材料的客重和剪 切模量、阻尼比。动力分析有限单元法的总体思路和静力情况基本一样，也是首先将计 算域划分为有限个有限大小的单元，单元之间在结点处互相连接，各个单元的质量平均 分配在该单元的结点上，然后分别求出各个单元结点的力与位移的关系，最后根据各单 元结点的力的平衡条件求出所有结点的力与位移的关系，进而形成整个振动体系的动力 方程。不过由于动力荷载与时间有关，相应的位移、应变和应力都是时间的函数，因此 在建立单元体的力学特性时，除静作用力外还需要考虑动荷载以及惯性力和阻尼力的作 用。在引入这些量的影响之后，就可以类似静力有限单元分析过程建立单元体和连续体 的动力方程，然后采用适当的动力计算方法进行求解。 土石坝的地震动力反应的直接结果是求得坝体及其地基土层中地震反应，包括各处 的反应加速度、动剪应力和剪应变。然后在已知的应力场中，确定关于安全系数的目标 函数，假定滑裂面形状，用搜索使安全系数最小的滑动面的方法，来找出最危险滑动面。 但是，如何确定安全系数的目标函数，如何根据应力和位移分析结果评价边坡的稳定性， 目前还没有完善的、为工程界接受的方法【6 j 。 1 3 本文所做的工作 本文主要是针对这些问题，用有限元来分析计算坝体在地震动力作用下的动力反 应，根据动力反应分析的结果，提出了把各个时刻的加速度作为惯性力，进行有限元强 度折减稳定性分析的方法，主要做了以下几个方面的工作： ( 1 ) 本文对动力边坡稳定性分析的拟静力法、n e w m a r k 滑块分析法( n e w m a r k s l i d i n g b l o c ka n a l y s i s ) 、有限元动力分析法等作了较为系统的总结和概括，并分析评价 了这些方法的优缺点。基于此，本文提出了用有限元强度折减法分析动力情况下的边坡 稳定情况。 ( 2 ) 分析介绍了有限元强度折减的基本原理、安全系数的物理意义、弹塑性本构模 型的选取、破坏评判的标准以及滑裂面的确定，并通过a n s y s 的a p d l 语言迸行二次开 发并编制了有限元强度折减边坡稳定性分析的程序，通过算例验证程序的合理性和实用 性，并与传统的极限平衡法进行比较。 f 3 ) 提出了在地震动力作用下如何应用强度折减法分析坝坡的稳定性的方法，详细 分析介绍了该方法的原理和计算方法，并通过算例分析验证该分析方法的合理性。 ( 4 1 对糯扎渡粘土心墙坝进行有限元静、动力反应分析，根据有限元动力反应分析 计算的结果，用有限元强度折减法对糯扎渡在地震作用下的稳定性进行计算分析，得出 结论。 大连理工大学硕士学位论文 2 土石坝动力分析的相关理论 2 1 初始应力的确定 在有限元动力计算中需要首先确定单元的初始应力。本文的初始应力采用有限元静 力计算的结果作为初始应力，静力计算采用我国沈珠江在吸收了邓肯模型和剑桥模型1 7 】 的优点，经过多年研究与不断完善而确立的南水双屈服面弹塑性模型，该模型既反映了 土与堆石的剪胀( 缩) 性、应力路径转折后的应力应变特性同时又可以采用目前通用的 常规三轴试验确定其模型参数，因而，便于推广应用。 为了模拟土石坝的施工过程，更好的反应施工过程中各阶段的应力和变形情况，采 用中点增量法【8 】，更好的体现了材料的非线性，因而更符合实际。 2 2 有限元动力计算方法 2 2 1 1 动力计算模型 目前用于地震反应分析的本构模型可分为两大类：一类是等价线性粘弹性模型，另 一类是( 粘) 弹塑性模型。弹塑性模型根据屈服条件、硬化规律、流动法则建立，在理论 上较为合理。但由于弹塑性模型比较复杂，虽然能够避免一些等价粘弹性模型的缺点， 但其参数的确定较为困难，建立及求解方程工作量大，所需时间长，而且其受应力路径 的影响也很大，目前缺乏合理的计算模型，工程上的实际应用较少。等价线性粘弹性模 型尽管存在一些缺点，但概念明确，应用方便，补充一些相关的计算模式，包括残余变 形模式，就能够全面分析地震反应，而且在参数的确定和应用方面积累了较丰富的试验 资料和工程经验，能为工程界所接受，实用性强。 本文的计算中，采用等价非线性粘弹性模型，根据堆石料动力特性试验得出的参数 和归一化动剪切模量q 。( 6 0 ) 一与动剪应变幅，关系以及等效阻尼比 与动剪应变幅y 的关系曲线进行动力计算。 确定瓯。的经验公式： g 。= 2 2 x 。p 。) 0 5 其中m = 半。 ( 2 1 ) 张彩双：士石坝动力稳定的强度折减法研究 表2 1 材料丘2 。取值 t a b 2 1v a l u eo f k 2 。o f m a t e r i a l 2 2 2 动力计算的运动方程及解法 在土石坝地震反应分析中，坝体及坝基可看作由许多质点组成的多自由度体系，经 过有限元法离散后，其动力方程可表示如下： 阻枷) ) + 【c 枷 + 医) = 一阻枷r ) ) 汜z ， 式中，f 村卜系统的质量矩阵，通常采用集中质量矩阵； c 1 一系统的阻尼矩阵； 陋7 一系统的刚度矩阵； 占( f ) 、8 ( t ) 、j ( f ) 一分别为系统各节点t 时刻的相对加速度、速度、位移向量 p 。( ，) t 时刻地基地震加速度； 系统的阻尼矩阵可采用， c - 口 m i + f l k 】 ( 2 3 ) 式中，口、卢一系数；口= 五国；= z ，。 可取体系的第一振型自振频率，五为阻尼比。各单元的材料不同，可取不同的阻 尼比，五还随剪应变，而变。 解运动方程一般采用数值积分法，一般有线性加速度法、n e w m a r k 常值加速度法及 w i l s o n - 0 法，这些方法对求解多质点体系和有限元的动力平衡方程式都适用。本文采用 w i l s o n - o ( e = 1 4 ) 数值积分法【9 】。 2 2 3 等价线性法的计算要点 ( 1 ) 静力有限元法计算出土体中各单元的震前的平均有效应力盯。 大连理工大学硕士学位论又 ( 2 ) 由本章公式( 2 1 ) 求出土体单元的初始l 动- 剪切模量g ，土体单元的初始阻尼比 根据经验取0 0 5 。 ( 3 1 将整个地震历程划分为若干时段。 ( 4 ) 对每个时段的动剪切模量进行迭代求解。 ( 5 ) 计算各单元的质量矩阵和刚度矩阵，形成总体质量矩阵和刚度矩阵，利用子空 间迭代法求坝体的基频甜，并出公式( 2 3 ) 计算单元阻尼矩阵，最后形成总 体阻尼矩阵。 ( 6 ) 选择地震波，输入地震加速度，并求解。 ( 7 ) 根据节点位移计算各单元的动剪应力和动剪应变。 ( 8 ) 重复( 4 ) 一( 7 ) 在每个时段内的动剪应变y 时程。 ( 9 ) 求出每个单元的动剪应变7 时程的最大值。，根据动剪应变7 。= o 6 5 y 一， 查g g 一，和五y 曲线得到新的g 、a 。 ( 1 0 ) 重复( 4 ) - ( 9 ) ，直到前次所用的g 和新的g 得误差小于0 1 。 ( 11 ) 重复( 3 ) 一0 0 ) ，直到整个地震历程结束。 n 2 ) 输出计算结果。计算结果可以输出各节点的位移、加速度的最大值，单元的 最大应力，给定节点的位移、加速度时程，给定单元的应力时程，给定时刻 的节点的加速度、位移和单元应力。 2 3 程序的实现 张彩双：土石坝动力稳定的强度折减法研究 动力计算程序流程 大连理工大学硕士学位论文 3 有限元强度折减法 3 1 概述 和传统的极限平衡法相比，有限单元法全面满足静力许可、应变相容和应力一应变 之间的本构关系，可作为一种更为严格的理论体系用于土体稳定分槲1 0 l 。同时因为是采 用数值分析方法，可以不受边坡几何形状的不规则和材料的不均匀性的限制，因此，该 方法是比较理想的分析边坡应力、变形和稳定性的手段。 虽然，目前考虑土体实际应力应变关系的坝体有限元应力分析已经十分普遍，但是 在有限元分析结果的基础上，如何评价坝体边坡的整体稳定性依然需要进一步探讨j 。 在已知应力场的基础上，假设滑裂面形状，进行最危险滑裂面搜索，来确定最小安全系 数，这个可描述为广义的数学规划问题，也是当前有限元分析稳定性一种比较普遍的方 法，但是这种方法跟传统的极限平衡法的弊端一样，都不能模拟边坡失稳的过程，也不 能确定边坡最先发生破坏的地方。安全系数计算公式不能充分考虑滑坡体的抗滑潜能， 也不能充分反映滑坡体的方向，安全系数偏大。 有限元强度折减方法最早由z i e n l d e w i c z ( 1 9 7 5 ) 提出并应用于边坡的稳定性分析【l ， 只是由于当时数值计算和计算机水平所限而未能得到发展。直到近几年来，随着数值计 算水平的提高和计算机的技术的发展，才逐渐为人们所用。许多学者在这方面做了大量 的工作，u g a i ( 1 9 8 9 ) 假定土体为理想弹塑性材料，采用强度折减有限元法对直立边坡、 倾斜边坡、非均质边坡及存在孔隙压力的复杂边坡的稳定性进行了较为系统的研究，指 出弹塑性强度折减法具有较强的适用性【j 列；m a t s u i 和s a n ( 1 9 9 2 ) 将d u n c a n - z h a n g 双曲线 非线性有限元与强度折减方法相结合，采用剪应变作为边坡破坏评判的指标，对人工填 筑边坡和开挖边坡分别进行了稳定性分析【l ，并与极限平衡法进行了对比研究；u g a i 和l e s h c h i n s k y ( 1 9 9 5 ) 将强度折减方法引入弹塑性有限元法中进行边坡的三维稳定性分 析，并与极限平衡法的计算结果进行了较全面的比较研究。尽管二者的理论基础、实现 手段完全不同，但是强度折减弹塑性有限元法仍可以得出与极限平衡法近乎“惊人”一 致的效果从而间接说明了有限元强度折减的可信 i 5 1 。 建立在强度折减基础上的有限单元法避免了临界滑裂面的搜寻，但需要进行多次有 限元试算，建议在试算中采用二分法计算方案，这样可以提高安全系数的求解效率。随 着岩土材料非线性弹塑性有限元方法的发展，强度折减法得到了广泛应用。这种方法的 优点就是不用假设滑裂面，可以模拟滑动发生的过程，安全系数也反映了强度储备的意 义。 张彩双：土石坝动力稳定的强度折减法研究 3 2 有限元强度折减系数法的基本原理 3 21 强度折减法的基本理论 在有限元静力稳态计算中，如果模型为不稳定状态，有限元计算将不收敛，基于此 原理，在非线性边坡稳定性分析中，边坡在自重或外荷载作用下失稳破坏时，塑性区会 贯穿整个边坡而形成一个滑裂带，整个边坡会沿滑裂带滑动，如图3 1 所示，从而使得 沿滑裂带以上的土体成为机动结构，因而其有限元计算结果将不收敛，强度折减正是基 于这一点提出来的【1 “。 图3 1 土坡失稳时塑性区分布情况 f i g 3 1p l a s t i c i t ya a r e $ ld i s t r i b u t i n gw h e nt h es l o p ei su n s t a b l e 有限元强度折减系数法的基本原理是：将坡体强度参数粘聚力c 和内摩擦角妒值同 时除以一个折减系数f ，得到一组新的c 、妒值，然后作为新的材料参数输入，再进行 计算。当计算不收敛时，对应的f 被称为坝坡的最小安全系数，此时的坝坡达到极限状 态，发生剪切破坏，同时可得到临界滑动面。 强度参数按下式进行折减： 、 c ，：三 f = a r c t a n ( t a n 州f ) ( 3 1 ) ( 3 2 ) 3 2 2 安全系数的物理意义 传统的边坡稳定极限平衡方法采用莫尔一库仑屈服准则，安全系数定义为沿滑动面 的抗剪强度与滑动面上实际剪力的比值，用公式表示如下： 大连理工大学硕士学位论文 f ( c + 仃t a n 妒) d ， f r d l 为说明有限元强度折减系数法求出的安全系数的物理意义，考察图3 2 所示的脱离 体a b c ，若土坡的应力分布为o vx ，y 1 ，相应滑面上正应力和剪应力为：盯和r ，那么 仃和f 的合力与脱离体所受的所有外力相平衡，即仃和r 的分布反应了脱离体a b c 上全 部外力在滑面a b 上的作用效果。 b c a 图3 2 滑动体受力示意图 f i g 3 2s k e t c ho f a c t i n gf o r c e 图3 3 微元t 的受力 f i g 33a c t i n gf o r c ei n f i n i t e s i m a ld i v i d i n g 边坡失稳时，土体的滑裂一般是渐进的，土体的破坏区域通常是从一点开始开始逐 步扩展，土体的内力不断调整，直到沿滑裂面整体失稳。假定当土体强度折减f 倍系数 时，土坡沿某一滑面滑动，则滑面上土体均达到临界状态，即对滑面上每一点均有： ( 3 4 ) 式中f 为沿滑面的剪应力，。为该点的抗剪强度，f 为强度折减系数。n n n i - - 点i 处于极限平衡状态时，在该点沿滑面曲线方向取微元长度向量龇，并令互= t f j 和 r = 一之争f ，( 如图3 3 ) ，贝。有： 一 l + 丑，= 0 ( 3 5 ) 式中r 为微元长度剪应力对应的内力，称为滑动力，置为与滑动力相平衡的折减抗 剪强度对应的阻滑力。由式( 3 5 ) 可得， 张彩双：土石坝动力稳定的强度折减法研究 r , m ，一 。i , a i ，= 0 ， ( 3 6 ) 式( 3 ，6 ) 表明，j 和( 一詈越) 代数和为零。其中t 为微元曲线的长度，于是对 于整个曲面有： z r , a l , x 等a l , = = 0 当一垒f 0 时，总存在一个f 使 f 4 根据式( 3 7 ) 可得： 进一步可得到， 等m f f j f f = z r i , a l , z r i 越 e = t d , i r 讲 i， ( 3 7 ) ( 3 8 ) ( 3 9 ) ( 3 1 0 ) 上式表明了沿滑动面各点均达到极限平衡时土体的平均强度折减系数，若采用摩尔 一库仑抗剪强度公式时即为安全系数定义式( 3 3 ) ，亦即与传统的安全系数有着相同的 物理意义。 3 2 3 弹塑性有限元的本构模型 尽管在有限元分析中可以考虑更加复杂的本构模型，但目前在工程分析中最普遍的 还是理想弹塑性模型【1 7 ，因为理想弹塑性的结果与工程师们熟悉的极限平衡法的结果最 ：赶连理上大学硕_ l 学位论文 具可比性。j o h n s o n 也证明了“如果对理想塑性结构体系施以简单加载，则极限荷载与 相同强度参数的刚塑性体的极限荷载相等”。另外，对一般的边坡和地下工程，主要的 计算控制量是岩土体的剪切破坏，因而对塑性区分布、大小的计算精度要求较高，而对 位移的计算精度要求较低，因而对上述工程用理想弹塑性模型是比较好的选择，计算结 果能更好的反映土体实际的工作性态。 在本文采用的有限单元中，土体的本构采用理想弹塑性模型，屈服准则选用了常用 的便于计算的广义米赛斯屈服准则，在国外被称为德鲁克普拉格准则( d p 准则) 。广 义米赛斯屈服准则在7 平面上的曲线为圆曲线，选择的参数( o f ，k ) 表达式不同时， 广义米赛斯屈服准则在疗平面上为一系列半径不同的圆。这些广义米赛斯屈服圆同相应 的摩尔一库仑屈服准则的不等角六边形位置关系不同，有限元计算中比较常用的有不等 角六边形的外角外接圆、内角外接圆、内切圆以及与摩尔一库仑屈服曲线围成的不等角 六边形面积相等的圆等。它们之间的相互位置关系可从屈服函数的数学表达式体现出 来，摩尔一库仑屈服准则和广义米赛斯屈服准则可以写成统一的数学表达式，如下式所 不： f = 瑾i i 七一j 2 = k ( 3 1 1 ) 这是一个通用表达式，通过变换o f ，k 的表达式就可以在有限元中实现不同的屈服 准则。o f ，k 是与岩土材料的内摩擦角、粘聚力和应力状态有关的常数【1 8 j 。对于不同的 屈服曲线，有不同的口，k 。 ( 1 ) 当口，k 满足下列表达式时 式中巳：；。 s 1 n 够 o f = 。- - - - 。= = ? 一 3 c o s 以+ 3s i n 以s i n ( o _ j ：；! 丝! ! ! 竺 4 3c o s 6 ，斗s i n 8 ，s i n 甲 ( 0 巳 州3 ) 这是摩尔一库仑屈服准则函数式，屈服曲线在万平面上为不等角六边形。 ( 3 1 2 a ) 张彩双：土石坝动力稳定的强度折减法研究 ( 2 ) 当口，j i 满足下列表达式时 舻万高2 s i n 口 七= ：二 , 3 ( 3 一s i n ) ( 3 1 2 b ) 这是广义米赛斯屈服准则之一，记为m 1 准则，屈服曲线在z 平面上为相应的摩尔 一库仑屈服准则不等角六边形的外顶点外接圆。 ( 3 ) 当口，k 满足下列表达式时 2s i n 口 一万商 _ j ：竺! ! ! 翌 4 3 ( 3 + s i n ( p ) ( 3 1 2 c ) 这是广义米赛斯屈服准则之一，记为m 2 准则，屈服曲线在万平面上为相应的摩尔 一库仑屈服准则不等角六边形的内顶点外接圆。 ( 4 ) 当盯，k 满足下列表达式时 s i n 拈舞 ( 3 1 2 d ) 这是经典德鲁克一普拉格屈服准则，是广义米赛斯屈服准则之一，记为m 3 准则，屈 服曲线在z 平面上为相应的摩尔一库仑屈服准则不等角六边形的内切圆。 大连理工大学硕士学位论文 ( 5 ) 当口，k 满足下列表达式时 2 4 3s i n p c z = ；5 = = = = = = = = = = = = = = = = = ： 2 拈( 9 “n 2p ) 七：型! ! ! ! ! 竺 2 廊( 9 “n 2 妒) ( 3 1 2 e ) 这是摩尔一库仑等面积屈服准则，是广义米赛斯屈服准则的一种形式，记为m 4 准 则，屈服曲线圆围成的面积等于相应摩尔一库仑屈服准则所围成的不等角六边形面积。 上述的摩尔一库仑屈服准则和广义米赛斯屈服准则在丌平面上的屈服曲线( 如图 3 4 ) 所示，从图上可以看出各个屈服曲线的相互对应关系，各个屈服曲线分别是摩尔一 库仑屈服曲线的外接圆、内接圆，内切圆以及等面积圆。而且可以看出外接圆的面积最 大，计算的塑性区最小，计算结果将偏于危险，内切圆面积最小，计算塑性区将是这些 屈服准则当中最大的，计算结果将偏于安全。尽管这几个屈服准则的面积不同，计算出 来的安全系数和塑性区也不相同，但是通过变换t 2 ，k ，可以在几个准则中相互转换。 本文计算中均采用摩尔一库仑等面积圆屈服准则，计算结果相对于其它准则最接近 摩尔库仑屈服准则，同时，不存在尖角，计算精度也不会受到影响。 点外接圆 o 1 图3 4 玎平面上不同口，k 的屈服曲线 f i g 3 4t h ey i e l ds u r f a c eo nt h ed e v i a t o r i cp l a n e 张彩双：土石坝动力稳定的强度折减法研究 3 2 4 确定滑裂面的位置 对于给定滑面的边坡，稳定性计算方法己经十分完善，或用极限平衡法，或用有限 元计算结果沿滑裂面进行积分j 1 9 j 。而实际情况往往是滑裂面的形状事先是不知道的。当 潜在的滑面被假定为圆弧时，临界滑裂面可以不借助复杂的优化技术就能搜索出满意的 结果。当滑裂面是任意形状时，问题就复杂多了。目前还没有普遍有效的方法可以很好 的解决任意形状滑裂面的确定问题。 用有限元强度折减法进行边坡稳定分析时，当达到临界状态时，会出现一个自下向 上贯通的一个塑性区，相当于出现了一个破坏带，滑裂面肯定出现在塑性区范围内，等 效塑性应变等值线应变示意图( 如图3 5 ) 所示。土体发生滑动破坏时，潜在滑裂面上 的点的塑性应变值肯定比该点两侧的塑性值要大，所以滑裂面肯定通过较大塑性应变的 转折点，也就是脊部点，拟合这些点所得到的曲线就是潜在滑裂面。通过脊背点拟合的 曲线即位最危险滑裂面。 图3 5 边坡临界滑裂面 f i g 3 5c r i t i c a ls l i po f d o w n s t r e a m 3 2 5 有限元强度折减法精度分析 如前所述，本构模型的选择合理与否会对有限元折减系数法的计算精度造成较大的 影响 2 0 】，除此之外，有限元本身也是误差的主要来源之一。 f 1 ) 网格的疏密 网格的疏密对单元精度的影响甚至大于单元类型的影响，对于精度较低的单元，可 通过加密网格来达到较高的精度。不同疏密的网格对计算结果具有不同影响。对于强度 折减系数法，有限元网格不能太稀，否则结果误差非常大，甚至可能不可用。但是网格 也不是越密越好，当网格密到一定程度的时候，再加密只会增加计算机时，而计算精度 增加不明显。通过大量算例证实，对于4 节点矩形单元。当单元密度达到每1 0m 2 不少 于3 个节点时，计算精度较为理想。 大连理工大学硕士学何沦文 ( 2 ) 边界范围 边界范围的大小在有限元法中对计算结果的影响比在传统的极限平衡法中表现的 更为敏感，在极限平衡法中只要所求滑动面在边界之内就不会对计算结果有影响，安全 系数只与划分的土条有关，而与土条外的区域无关，有限元则不然，边界的