（水工结构工程专业论文）边坡稳定极限分析有限元法下限解的研究.pdf

摘要 基于土体塑性极限分析理论而产生的边坡极限荷载上、下限解答的问题是目 前工程界广泛考虑的问题。其中，上、下限解答与边坡稳定极限荷载的真实解的 逼近问题更成为了岩土工程中关注的焦点，因此本文就极限分析有限元的下限解 答做了初步探讨。 在分析边坡稳定极限荷载的问题中采用的是s l o a n 等学者提出的极限分析有 限元法，极限分析有限元法将塑性力学中的极限分析方法与传统的有限单元法相 结合，通过对结构物的离散化以及约束条件分析，找到静力容许的应力场。并在 此基础上，利用数学规划方法来确定目标函数的最大值( 或最小值) ，从而得到 边坡稳定极限荷载的下限解答。 本文介绍了极限分析有限元法下限解程序的基本原理，并在已有研究工作成 果的基础上，对极限分析有限元法下限解的实现过程编写了相应的程序。在程序 开发过程中，着重解决了下限解有限元法与传统有限元法的衔接问题，下限解有 限元法转化为标准线性规划问题。在对数学模型中线性优化求解时，选用了功能 强大的数学工具软件来进行线性优化计算，同时通过对典型例题的分析对比，验 证了程序的正确性和数值结果的可靠性。 经过改进后的极限分析有限元法下限解的程序，不仅弥补了目前己丌发程序 中的不足，而且使传统有限元分析技术与极限分析有限元理论结合起来，对工程 研究中的数值计算具有很高实用价值。采用软件进行计算时避开了对求解初始值 的选取，对于简单问题能够得到满意的答案，如何应用此方法解决更为复杂的工 程实际问题将是今后进一步深入研究的重点。 关键词边坡稳定；极限分析；下限解；有限元；线性优化 a b s t r a c t 1 1 1 eu p p e rb o u n da n dl o w e rb o u n dl i m i t e dl o a d i n go f s l o p es t a b i l i z a t i o n w h i c h a r eb a s e do nt h el i m i ta n a l y s i st h e o r yo fs o i l p l a s t i c i t y , h a v eb e e nc o m p r e h e n s i v e l y r e s e a r c h e dr e c e n t l y e n g i n e e r sh a v ea l w a y sb e e nf o c u s i n go nt h ep r o b l e m so fu p p e r b o u n da n dl o w e rb o u n d sc o m p a r a b i l i t yt ot h ee x a c tv a l u eo fs l o p es t a b i l i z a t i o n t h i s p a p e rw i l ld e a lw i t ht h el o w e rb o u n do fl i m i ta n a l y s i sb yf i n i t e e l e m e n tm e t h o d 1 1 1 em e t h o do fl i m i ta n a l y s i sb yf i n i t e e l e m e n t ，w h i c hi sp r o p o s e db ys l o a ne t c ， i su s e di nt h ea n a l y s i so fl i m i tl o a d i n gf o r t h es t a b i l i t yo fs l o p e i ti sm a d eu po fl i m i t a n a l y s i sa n dc l a s s i cf i n i t ee l e m e n t w ec o u l df i n das t a t i c a l l ya d m i s s i b l es t r e s sf i e l db y s c a t t e r i n gt h es t r u c t u r ea n da n a l y z i n gc o n s t r a i n te q u a t i o n ，a n df i n dt h em a x i m u mo r m i n i m u mo fo b j e c tf u n c t i o nv i am a t h e m a t i c a lp r o g r a m m i n g , t h e ni nt h i sw a yw ec a n o b t a i nt 1 1 es o l u t i o nt ot h el o w e rb o a n dv a l u e t 1 1 i sp a p e rc o n c l u d e st h et h e r o mo ff i n i t e e l e m e n tm e t h o do ft h el i m i ta n a l y s i s a n dc o r r e l a t i v ea c h i e v e m e n t s ，a n dd e v e l o p et h ep r o c e d u r eo ft h ef i n i t ee l e m e n t m e t h o do fl o w e rb o u n d t h ep r o b l e m so f j o i n tb e t w e e nf i n i t ee l e m e n tm e t h o do fl i m i t a n a l y s i sa n dc l a s s i c a lf i n i t e - e l e m e n tm e t h o d ，a n dt h et r a n s f o r m a t i o no fs t a n d a r dl i n e a r p r o g r a m m i n gh a v eb e e ns o l v e dw e l l i nt h i sp r o c e d u r e i nt h ec o u r s eo fl i n e a r p r o g r a m m i n g , t h em a t h e m a t i c a ls o f t w a r eh a sb e e nu s e de f f e c t i v e l y a n dt h er e s u l ti s s a t i s f y i n gt h r o u g hs o m ec l a s s i c a le x a m p l e s 1 1 1 er e v i s e dp r o c e d u r eo fl o w e rb o u n df i n i t e e l e m e n tm e t h o do fl i m i ta n a l y s i s b e c o m eu s e f u l n o to n l yc a ni tm a k eu pt h ed i s a d v a n t a g eo f e x i s t i n gp r o c e d u r e ，b u t a l s oi tc o u l dc o u n e c it h et h e o r e mo fl i m i ta n a l y s i st oc l a s s i c a l i t ce l e m e n t n e c h o i c eo fi n i t i a l i z a t i o nh a sb e e na v o i d e db yu s i n gm a t h e m a t i c a ls o f t w a r e ，a n dc a n p r o d u c es a t i s f i e da n s w e r si ns i m p l ee x a m p l e s b u tt h e r e a r es t i l lm a n yp r o b l e m s s h o u l d b es o l v e di nf u t u r er e s e a r c hw h e nt h i sm e t h o di su s e di n c o m p l e x c i r c u m s t a n c e k e y w o r d ss l o p es t a b i l i z a t i o n ；l i m i ta n a l y s i s ；f i n i t ee l e m e n t ；l o w e rb o u n d ； l i n e a ro p t i m i z a t i o n 独创性声明 本人声明所呈交的论文是我个人在导师指导下进行的研究工作及取得的研 究成果。尽我所知，除了文中特别加以标注和致谢的地方外，论文中不包含其他 人已经发表或撰写过的研究成果，也不包含为获得北京工业大学或其它教育机构 的学位或证书而使用过的材料。与我一同工作的同志对本研究所做的任何贡献均 已在论文中作了明确的说明并表示了谢意。 签名： 毒毕 同期：! 生3 关于论文使用授权的说明 本人完全了解北京工业大学有关保留、使用学位论文的规定，即：学校有权 保留送交论文的复印件，允许论文被查阅和借阅；学校可以公布论文的全部或部 分内容，可以采用影印、缩印或其他复制手段保存论文。 ( 保密的论文在解密后应遵守此规定) 签名：强导师签名：纽开期： 导师签名：z ! 二! 兰开期： z o o7 ，6 ， 第1 帚绪论 1 1 极限分析研究概述 第1 章绪论 到目前为止，国内外进行土体极限分析的方法有两种：一种是基于弹塑性模 型的变形分析方法1 2 l 【3 】，另一种是基于刚塑性体模型的极限分析方法。 变形分析法首先要建立土体的弹塑性本构关系模型，结合土体的平衡方程、 几何条件，根据具体问题的边界条件、初始条件及荷载历史，应用解析解法或数 值解法，逐步求解土体应力和变形，从而求得土体的极限承载力。然而，土的变 形特性十分复杂，除了不连续性，非线性的本构关系以及在破坏时呈现的体胀和 软化、大变形等特性影响外，还有天然地基土的各向异性、动力荷载，地震作用 下的特性，与时间有关的特性以及区域性土等等。要想使土的各种复杂特性都在 模型中得到反映，必然会使模型变得十分复杂而无法应用。目前岩土材料的本构 关系还很不成熟，虽然已有的各种模型数以百计，但均处于发展之中，各种弹塑 性模型还有待于进一步验证。另外，应用弹塑性一般理论求解时，在机理及数学 上的严密性未必就能换取计算结果的改善，因为参数选取的关系影响极大。同时 计算方法本身也存在计算误差，常常会遇到无法收敛的问题。因此，用变形分析 法来求解土体极限荷载的问题显得十分困难和复杂，在实际工程中的应用范例也 较少。 基于刚塑性模型的极限分析方法在土体分析中，回避了在工程中最不容易弄 清的本构关系表达式，而是直接从土体达到塑性屈服状态直至破坏的角度去分析 土体的受力。正是因为这一特点，使得多年来大量的研究学者致力于这一理论方 法的研究工作。从2 0 世纪5 0 年代柏林大学的d c d r u c k e r 和w p r a g e r 把静力场和运 动场结合起来，提出极值原理、建立了土体极限分析理论【4 】伊始，很多的学者在 此基础上做了大量研究工作。1 9 7 5 年，美国学者w fc h e n ( 陈惠发) 在其专著【5 l 中系统深入地阐明了极限分析理论在土工问题中的应用。1 9 7 6 年，d a v i d n h u t u l a 将有限元法和极限分析法想结合，针对服从m i s e s s 屈服准则的二维理想塑性平面 结构提出了一套数学规划问题的求解算法【“。在1 9 8 1 年，r y s z a r dj i z b i c h i 研究了 用上限定理来确定土坡安全系数的图解法【”，并与极限平衡理论的若干种方法做 了分析对比。1 9 8 2 年，德国斯图加特大学教授p g u s s m a n n 在e d d m o n t e n 岩石力学 数值方法会议上提出了滑移单元法【8 】i 【9 l ( k i n e m a t i ce l e m e n tm e t h o d ，简称k e m 法) ，采用有限元理论将结构离散化，利用最优化方法求得上限解。1 9 8 5 年r 本 两位学者k a t s u h i k oa r a i 和m i t s u on a k a g a w a 联合运用极限平衡法中的薄片法与极 限分析法的下限原理采用数学规划方法得出了边坡沿己知滑动面的稳定安全系 数i i “。1 9 8 8 年美国学者a d r e s c h e r 和c c h r i s t o p o u l o s 用上限分析方法研究了无限 北京t 业j 、学t 学埘j 卜学位论文 均质自由超载土坡的稳定问题j 。1 9 8 8 年澳大利亚学者s w s l o a n 将土力学中的 下限分析法与有限元法相结合，考虑了单元之间相邻边上的应力间断的条件，在 理想塑性平面应变问题下将下限分析的解答转化为求线性规划的问题1 2 1 。随后 s w s l o a n 叉在1 9 8 9 年提出了有限元塑性极限分析上限解答的算法，从此形成 了完善的极限分析有限元理论的研究方法，为以后的研究奠定了坚实基础。 1 2 课题的研究意义 滑坡是一种常见而且重大的地质灾害，每年的山体滑坡或泥石流都要给人民 的生命，财产带来巨大损失。因此，研究边坡稳定分析的理论和计算方法一直都 是岩土工程师们的一项长期而艰巨的任务。现阶段，用于边坡稳定分析1 1 4 j 的主要 理论基础还是极限平衡理论【1 6 】，基于这一理论基础上的边坡稳定分析方法把 土体划分为连续的垂直条块，通过对条块或者整个滑坡体的力的平衡和力矩平衡 求解边坡的安全系数。由于边坡稳定系统本质上是一个静不定的系统，故而传统 的极限平衡法对条块界面上的力的作用方向和作用点均做了相应假设，从而使问 题变成静定可解。尽管边坡稳定的极限平衡法有其理论上的局限性，但是由于其 理论基础简单，容易被设计人员理解，己在实践中得到了广泛的应用。 寻求边坡工程问题中理论上严密、计算结果合理的真实解答一直是学术界和 工程界的研究目标。近年来，随着计算机技术硬件水平的飞速发展和大型数值分 析软件的出现，用数值方法1 1 7 】1 1 8 1 确定边坡稳定问题的严密解答( 见文献【1 i j ) 就变得可行。一种选择是应用能考虑土的弹塑性模型的有限元方法模拟边坡的整 个加载历史，从而确定边坡稳定的极限荷载，本文中把这一方法称为传统有限元 法1 2 2 】【2 3 j 。另外一种选择就是应用经典塑性力学的极限分析理论【2 4 1 ，即上限原理 和下限原理去直接确定边坡承载力。虽然上、下限原理只适用于考虑相关联流动 的刚塑性模型，但是对于一些服从相关联流动法则的边坡问题，极限分析原理也 能从上限和下限两个方向来逼近问题的真实解答f 2 5 】，而对于那些不服从相关联流 动法则的问题，极限分析方法也能得到非常有益的解答。 尽管我们可以用不同的数值方法进行严格的极限分析，但是有限元极限分析 技术可以用来进一步地提升极限分析方法的可用性和灵活性，这一点已经得到了 岩土工程界的广泛认可。极限分析有限元法与传统有限元法最大的不同在于：前 者能从上、下限解答来逼近问题的真实解答，从而可以提供一个可以参考的真实 解的范围，而后者只能提供所求问题的一个近似解答。 在目前的研究中存在的一个最大问题就是边坡极限荷载的上、下限值与真实 解的逼近程度问题。克服这一问题的关键在于寻找有效的求解最大值或是最小值 的数值优化算法，因而对本课题展丌研究也显得尤为重要。 第1 章绪论 1 3 本文所做的工作 本次课题研究的主要内容如下： 简要论述塑性力学中的下限原理，并在此基础上，结合有限元理论深入 阐述极限分析下限解的基本原理及实现方法。其重点在于形成恰当的目标函数， 将工程问题转化为数学问题。 根据这一求解机理，编写了相应的程序来进行解答。本文就目i j 已有程序 中存在的不足，编写了相应的前处理程序，与现有极限分析有限元结合起来， 完成对数学问题的求解过程。 运用数学规划理论中线性规划的方法来完成对该数学问题的求解过程。本 文中采用了数学工具软件m a t l a b 中提供的优化求解函数来实现这一目标。同时 对m a t l a b 中所采用方法的原理也做了相应说明。 通过实际算例，对边坡稳定的极限荷载进行求解，验证本方法的j 下确性， 并将得到的解答与理论解析解相对比，分析所求解答的准确程度以及产生误差 的原因。 1 4 本章小结 本章回顾了极限分析研究理论发展的概况，指明了这一领域研究的现状， 在此基础上提出了对边坡稳定极限分析有限元法下限解进行研究的必要性，并 对本次研究的具体工作做了系统安排。 第2 章型忡力学下跟定理 2 1 概述 第2 章塑性力学下限定理 在某一边坡表面作用有一个外荷载，不断增加该外荷载值一直到l 边坡最 终达到极限状态2 6 1 。此时的外荷载r 称为极限荷载。对极限荷载，可作如下定 义： s t q c a ) 加载过程和极限荷载( b ) “真实的破坏机构” 图2 - 1边坡加载和破坏 f i g u r e 2 1l o a d i n ga n di n v a l i d a t i o n 如图( 图2 1 ) 所示，为边坡荷载丁与滑动位移s 关系的曲线图。当边坡表 面某点的垂直荷载r 从零增大到某一值时，滑动位移一直为零，此后随着t 的 增大，滑动位移s 也随之变大。在现场进行荷载试验时，即可得到类似曲线。该 曲线最终变为垂直的时刻，即为定量确定地基极限承载力的依据。另一方面，我 们认为，r 达到极限状态时，边坡内出现了一个塑性区0 和一个弹性区e 。在塑 性区n 内，处处达到极限平衡，而在弹性区内，则各点均处于弹性状态。在 弹性区e 内，刚性体假定认为各点的位移值为零。在o 和e 之间，出现了一个 位移问断面r ，通常称其为滑裂面。为了论述方便，称这一处于极限状态的塑性 区0 和位移间断面r 为“真实的破坏机构”。当然，这一提法只是从理论意义上 将其与以后要出现的各种“假设的破坏机构”区分开，并无意与岩土工程中实际发 生的破坏状态直接挂钩。实际上，当土体沿着某一滑裂面失稳破坏时，滑裂面上 的点是处于完全塑性状态的，而滑体内部并不是所有点都处于完全塑性状态，可 能有部分区域还处于弹性状态或是接近于塑性状态。 2 2 下限定理 下限定理认为：在不违反屈服准则的条件下，对应于任何一个静力容许的应 力场都可以求得个相应的外荷载，这个外荷载总是小于或等于极限荷载，是极 限荷载的下界( 或安全) 的估计。一个静力容许的应力场要同时满足应力边 界条件、应力平衡条件及屈服条件( 应力空问的应力必须位于屈服面内) 。 q o 幽2 - 2 证明f 限定理不恿幽 f i g u r e 2 - 2 s k e t c hf o rp r o v i n gt h el o w e rb o u n dt h e o r e m 上图( 图2 2 ) 中，斜纹区域表示。真实的破坏机构”，散点区表示满足静力 平衡条件的应力场。假设处于极限状态的“真实的破坏机构”相应的表面荷载为 l 滑面上和滑体内的应力用表示。下限定理已经指出，对于同样的破坏机构 ( 即同样的区域q ，和厂) ，如果存在另一个静力容许的应力场仃：，在滑体内， 总有厂( 盯；) 蔓o ；在滑面上总有，( 盯；) = o ，那么，与这一应力场相平衡的表面 荷载r o 不比r 大。 假定“真实的破坏机构”对应的应力场为o r # ，设外部荷载的微小增量导致 的塑性应变率为毛和u ，那么，根据虚功原理有 n 勺咖+ r o r f 勺凼= f t u d s + f o w u d v ( 2 - - 1 ) 上式中l 为t 作用的边界。现在，对假设的静力容许的应力场0 及相对应的荷 载r o 也用毛和“作虚功原理的操作， n 0 勺咖+ f r 0 勺d s = t o u d s + f n w u d v ( 2 2 ) 式( 2 2 ) 减去式( 卜1 ) ，可得 j n ( o r 一司勺d v + j r ( o r 口一0 勺d s = ( t t 。) j u d s ( 2 3 ) 上式左侧的第项为滑体内各点的内能耗散，第二项则为滑面上内能耗散，根据 d r u c k e r 公设，此两项均为非负，式子( 2 3 ) 左侧总是证的，故有 t t o ( 2 4 ) 从而下限定理得证。 值得说明的是，在边坡稳定极限分析领域，我们引入下限定理增加了较多的 限制，与极限状态相比较的这个静力容许的应力场盯：和相应“真实的破坏机构” 的应力场盯。具有相同的塑性区口和滑裂面，1 。 2 3 塑性极限分析的基本假定 塑性力学下限定理为边坡稳定的极限分析提供了理论基础，在这一过程中引 入了以下基本假定。 ( 1 ) 材料是理想的刚塑性体。因此，在静力法分析中，屈服函数不随应力的 改变而发生变化。 ( 2 ) 结构物即使在临近破坏状态时，其变形也很小。因此，平衡方程在结构 物变形前后不变。 2 4 下限解的极限分析方法 针对塑性力学中的下限原理，采用的极限分析方法我们称之为静力法。 静力法是以节点应力为未知量，在构造被研究土体静力容许的应力场时，不 必考虑破坏机构条件，而只要满足平衡条件、应力边界条件以及应力屈服条件即 可。依据塑性力学下限定理，从众多静力容许的应力场所对应的荷载值中找到最 大值，就是找到了最接近于真实极限荷载的荷载值。 2 5 本章小结 本章从极限状态的概念出发，说明了极限平衡以及下限定理的含义，对塑性 力学中的下限原理进行了证明，指明了极限分析下限解的分析方法。 第3 章下限解有限单元法的实现原理和基本方法 下限解有限单元法是以塑性力学的下限定理为基础的，且下限定理认为：在 不违反屈服准则的条件下，对应于研究对象中任何一个静力容许的应力场都可以 求得一个相应的外部荷载，这个荷载总是小于或等于所研究对象的极限荷载，并 且是极限荷载的下限估计值。同时，每个静力容许的应力场又要满足应力边界条 件、应力平衡条件以及屈服条件。 在下限解有限单元法1 2 2 l 中，我们先将节点的应力作为未知量，构造同时满足 节点应力平衡条件、应力边界条件以及应力屈服条件的表达式，使得这些应力值 同时满足上述所有条件。然后计算每组应力状念所对应的荷载，从众多荷载中找 到最大的那个荷载，从而就能找到最能接近真实极限荷载的荷载值。 3 1 结构物的离散 众所周知，工程实际中的地质条件是非常复杂的。地基工程或是边坡工程所 研究的土体对象可能将会是砂土和石快的混合体，其力学性质具有不均匀性、复 杂性以及易变性，因而在进行数学模拟的过程中，我们简化认为所要研究的土体 结构都是由性质相同的砂土所构成。 离散时，我们采用常见的三角形3 节点单元对结构物进行离散化，3 个节点 i ，m 的编号顺序按逆时针方向，具体形式见下图( 图3 一1 ) 。 y 图3 - i 三角形3 1 i 点示意i ! i f i g u r e3 - 1 3 - n o d e st r i a n g l em e s h x 北泉t 业人学t 学坝i 掌位论j ( 三角形3 节点单元常被用作模拟平面应变条件下的应力场，因此本文沿用了 这一单元划分形式。单元上所有点的应力状态变量都是线性变化的，每个节点都 只含有3 个未知的应力变量盯，、1 7 ，和f 。，单元上每一点的应力大小都可以通过 下面的式子( 式3 1 ) 计算得到： 仃j = n i f y “ 仃，2 舌f ，( 3 - - 1 ) r w 。i n , r 圳 其中盯。、口，和f ，代表了某点的节点应力，n i 为位移插值函数，它是单元内该 点坐标的线性函数，由以下公式给出： 其中， i = 寺( 最怕x y ) j ：= 击( 善：幌j 坞y ) ( h ) 3 = 击( 最工嵋j ，) f 石= x z y 3 一屯y 2 易： y l 一 ) ，l k ： _ y 2 一工2 一 仇= y 2 一y 1 叩2 = y ) 一y 仉= j ，l y 2 式中的a 代表该三角形单元的面积，并1 t2 a = i 仉“- r 1 2 玉i 。下图( 见下页 图3 2 ) 中给出了三角形单元线性应力的网格。和通常所用的有限元网格的形式 不同的是，每个节点对某一特定单元来说编号都是唯一的，而且多个节点可能会 分享共同的坐标。在两个相邻单元的公共边上，允许出现应力的不连续性，但是 根据力的相互作用原理，公共边上的法向应力和切向应力是大小相等、方向相反 的。 3 2 应力平衡条件 在平面应变状态下的应力平衡条件可用下式表示： )3 3( 屯毛 一 一 一 b 屯 = = = 色白 第3 辛下限舯南限甲，i 法的实现原理和壮4 】= 方法 一0 0 x + 堡：o 叙 匆 一o a r , + 笠：， 砂 西 。 式中拉应力为正值，7 为土体的天然容重，坐标系采用向右旋转的笛卡尔坐 标系统。将公式( 3 一1 ) ，( 3 2 ) 和( 3 3 ) 代入公式( 3 4 ) 中，得到平衡方 程节点应力的约束条件表达式： 【鸽】p - 屹 ( 3 5 ) 式中， 心卜寺心曩：2z 曩：3 篾 盯) 7 = 盯二盯；lr j 。y o r 二2 ：f ；：o a ，) ( 3 - - 6 ) 制 其中a 。表示该三角形单元的面积。这样，对每个三角形单元来说，平衡条 件就能够用包含节点应力变量的两个等效约束方程来表示。 y 图3 2 极限荷载分析法的线性应力三角网格图 f i g u r e3 - 2 m e s ho fl i n e a rs u e s st r i a n g l e sf o rl i m i tl o a d sa n a l y s i s 连续面 北京t 业人学t 学倾f 学位论殳 3 3 应力不连续面平衡条件 为了满足相邻单元在边界上的静态容许的不连续性，在节点应力上加强额外 的约束是很有必要的。这种静态容许的不连续性允许三角形单元的应力是不连续 的，但是它必须保证剪应力和节点的j 下应力分量是连续的。参照下图( 图3 3 ) ， 葛 x 图3 - 3 作用在平面上的止应力和剪鹿力分布图 f i g u r e3 - 3 r e s o l u t i o nf o rs t r e s si n t on o r m a la n ds h e a rc o m p o n e n t sa c t i n go nap l a n e 可以看到，当正应力和剪应力作用的方向是在与工轴呈p 角的斜面上时，其数值 大小可由下面公式( 式3 7 ) 得到： h = s i n 2 0 c r 。+ c o s 2 0 c r y s i n 2 0 x 卜扣胁，巾加町一 帖7 对于相邻的两个单元之白j ，如图3 - 4 ( 见下页) 所示，两个三角形单元a 和 b ，其中边界d 为不连续日j 断面，它是由两组节点对( 1 、2 节点与3 、4 节点) 来定义的。在边界d 上，每个点的应力值可以由下面公式得到： j 盯：= 砖( 3 - - 8 ) l f 4 = f 6 y 凹3 - 4 相邻三角形间静态容许不适续间断面 f i g u r e3 - 4s t a t i c a l l ya d m i s s i b l es t r e s sd i s c o n t i n u i t yb e t w e e na a j a c e n tt r i a n g l e s 因为应力在每个单元的边界上是随着线性变化的，因此上面的条件等价于以 下的约束条件： 仃：l = 盯。b 2 盯二= 盯l ( 3 - - 9 ) f ? = r ； f ；= 巧b 将以上各式代入公式( 3 7 ) ，就可以集合成一个矩阵平衡方程。 【一： p 。 = 6 三) ( 3 1 0 ) 式中， k 】= r o 一； 阱醢岛羞一篡 p 。 7 = 盯：f 刍吒b ：q b ：f ，b ：d 二，巧，盯，d 。盯，d 。r ，d 。) 6 三) 7 = 000o 因此每个静态容许不连续性的平衡条件在单元的边界上产生了等价的节点 应力限制条件。 3 4 边界条件 在岩土工程的许多问题中，产生的应力边界条件都遵从下面的形式： j o r 2q ( 常量) ( 3 一1 2 ) i f = t( 常量) 由于沿着三角形单元边界上各点的三个应力分量盯。o r ，和r 。都是作线性变 化的，因而可以借助于局部坐标的形式来表示每一点的应力状况，各点的应力值 由公式( 3 一1 3 ) 可直接算得( 参见图3 5 ) 。 y 吒= q + ( 9 2 一哂) 舌 f = f l + ( 乞一f 1 ) 善 图3 - 5 应力边界条什 f i g u r e3 - 5 s t r e s sb o u n d a r yc o n d i t i o n 盯：2 吼+ ( g z 咱) 孝( 卜1 3 ) 1 f ，：f l + ( f ：一f ，) 掌 、h “ 式中， ，需要做特殊边界处理的三角形单元e 的边界( 用1 节点和2 节点表示) ； 善沿边界，的局部坐标，且善( o ，1 ) ； q l , q 2 节点1 和节点2 上的正应力( 以受拉为正) ； t l , t 2 节点1 和节点2 上的剪应力( 沿边界面以顺时针方向为正) 。 由方程式( 卜1 3 ) 所确定的边界条件要满足下面的要求： f 盯二= 叮 j 吒e z = g ： l r i = r 。 l f 2 。= t 2 令只为，边与工轴的央角，代入到式( 3 7 ) 中，得到应力边界条件引发的 约束条件，如下： 【一乞】p = ( 叱)( 3 一1 5 ) 其中， 乩啦； l s i n 2 研c o s 2 岛一s i n 2 0 , i 纠2 b n z 只扣帅s 2 叫 盯 7 = 盯二盯品f 孟：盯；：r ；：) i 7 = g l t i9 2t 2 这样，在需要进行边界处理的每条边界，上就会产生以节点应力为变量的四 个等价约束条件。 3 5 屈服条件 假定拉应力为正值，并假设研究的问题属于平面应变问题，m o h r - c o u l o m b 屈服准则就能用下面的式子表示出来： f = ( 盯。一盯j ) 2 + ( 2 r 掣) 2 一【2 c c o s 一( 盯。+ 盯y ) s i n 妒】2 = 0 ( 3 1 7 ) 为了保证应力值不会超出屈服条件，并能够完全满足下限极限定理的需要， 在每个三角形单元中都必须保证f 值大于或等于零。我们需要把此类规划问题转 变成为具体数学公式来实现，所以假设应力屈服函数为线性函数，用线性规划的 思想来近似处理等式( 3 一1 7 ) 是可取的。为了确保得到的方程式的解正好是破 坏时真实荷载的下限值，线性优化的屈服面必须位于应力空间上m o h r - c o u l o m b 屈服面以内。 北京t 业人学t 学硕i 学位论立： 摩尔 令 x = ( 一q ) 伊一( o r ，+ 盯，) s i n 尹 图3 - 6 从内部线性化近似的m o h r - c o u l o m b 屈服面 f i g u r e3 - 6 l i n e a rm o h r - c o u l o m by i e l df u n c t i o n ( p = 6 ) x = o x o r ， y = 2 r 。 ( 3 1 8 ) r = 2 c c o s # 一( 盯。+ 仃，) s i n # m o h r - c o u l o m b 屈服条件可以用公式x 2 + y 2 = r 2 来表示。变量x 和y 的含 义，在上图中( 见图3 6 ) 用一个应力圆来描述。m o b r - c o u l o m b 屈服面用一个有 p 条边界p 个点的内接j 下多边形来近似模拟。 第k 个x 、y 的局部坐标和第七+ l 的最高次项由下面的式子给出： 有： 对于每个应力点来说， x 产r c o s x ( 2 k - 0 p k ：r - s i n 剑 p ( 卜1 9 ) 丘+ l ：r c o s y c ( 2 k + 1 ) k + ：r 咖丝业 p 相关的雒标x 、y 倌都应该位于屈服多边形以内，目 ( x i + 。一x ) g y ) 一( r 一x ) ( 圪+ i y ) o ：k = 1 , 2 ，p ( 3 2 0 ) 将公式( 3 一1 8 ) 、( 3 一1 9 ) 代入到上面式子中，得到 一t 盯，+ b i 仃，+ c i f 咿d ：k = 1 , 2 ，p ( 3 2 1 ) 一1 6 其中 4 ：s 2 k x + s i n 妒c 。s 三 b - ：s i n c o s 一7 - - c o s 2 k z p p(3-22) c k ：2 s i n 2 k x d ：2 c - c o s 西c o s 三 为了保证等式( 3 2 1 ) 对所有的网格单元都成立，必须加强在每个三角形 单元上各个节点的限制条件。另外需要特别注意的是，并不是要求所有单元的凝 聚力都一样。为了证明这一点，我们假设凝聚力随三角形单元的线性变化为： 3 c ：y n i c 其中，n i 为等式( 3 _ 2 ) 确定的线性插值函数，c i 为节点f 处的凝聚力。 等式( 3 2 1 ) 中的常量d 是坐标x 和y 的线性函数，用下式表示： 3 d = q ( 3 2 3 ) 其中口= 2 ge o s # c o s l ，将等式( _ 1 ) 、( 3 2 3 ) 代入( 3 2 1 ) 式中，每个 p 三角形单元节点处所需要的应力约束条件为： 3 3 苫m ( 4 盯，。+ b k o y i + c t f 圳) 茎苫i d i ；七= l ，2 ，p ( 3 - 2 4 ) 如果我们在每个节点上加强约束条件，那么等式( 3 2 4 ) 满足三角形单元 内的所有的点，如下： a k a “+ b k o “+ c k f 斜is d 0k = 1 , 2 ，p 下面来证明此式。我们注意到，如果定义o n i 1 ，那么在单元内的所有 点就会有： m ( 一i c r “+ b t c r 一+ c t f 驯) n i 口；k = 1 , 2 ，p ：f = 1 , 2 ，3 ( 3 2 6 ) 简化后写成： 3 ， 善f ( 一t 盯“+ 口i 盯+ g r 州) s 苫，d i k ：1 , 2 ，p ( 3 - 2 7 ) 北柬t 业人学t 学坝l 掌位论文 倘若考虑的是纯粘性土壤，便可通过将内摩擦角西赋为零值后代入到公式 ( 3 2 2 ) 中得到。由于在每个三角形中的凝聚力是按照线性变化的，所以凝聚 力随深度增加而增大的纯粘性土体的问题也可以得到很好的解决。因此线性化后 的屈服准则施加在节点之上的应力限制条件可以通过矩阵的形式来表示，如下： 【一二】p = 彰。 ( 3 2 8 ) 其中， 盯甲= 屹 7 = 2 q c o s # c o s2 c ic o s 驴c o s l r ，2 q c o s # c o s 三 ppp k 】- 4马g 4 马c 2 4 毋g ( 3 - 2 9 ) 系数4 ，b ，c 都由公式( 3 2 2 ) 给出，q 表示在节点i 处的凝聚力。 因此线性化后的屈服条件能够在每个节点应力向量上产生p 个等价的约束条件。 对于每个网格，如果有e 个三角形单元的话，就会有3 e 个节点，产生的总的不 等式约束条件就有3 p e 个。 3 6 目标函数 对于大多数平面应变的土力学问题，我们希望能找到一个满足静态容许应力 场并且能够反映出整个结构所能承受的某个最大荷载，用下式表示： q = 6 n d a ( 3 - 3 0 ) 其中q 为极限荷载，h 是出平面的厚度，o r 是作用在边界s 上某处的j 下应力 值。图3 7 ( 见下页) 通过一个三角形单元的边缘的节点l 、节点2 说明了盯作 用的最大值。因为假定应力沿各个三角形的分布都是线性变化的，所以q 可以表 达为： 第3 章下限舯仃限单儿注的实现原删手门皋奉方让 q = 譬( + o r n ：) ( 3 3 1 ) 式中l 是三角形单元的边界s 的长度，o r 盯。：是三角形单元节点1 和节点 2 上的j 下应力。 令见为s 边界与x 轴的夹角，假设出平面方向的厚度为单位厚度，即h = l ， 将( 3 7 ) 式代入上式，得到： q = p 7 ( 盯 ( 3 0 2 ) 其中， y p o r j 糍竺等咖c o s “啦卜，) 1 = 一tf ；一盯二c r ；：f ；：j x 矗为出平面厚度 q = 等( 吒+ 吒，) 法向荷载 图3 7 某条边界上的法向荷载 f i g u r e3 - 7 l o a di nad i r e c t i o nn o r m a lt oab o u n d a r ye d g e 在线性规划表达式中， c 是目标函数的系数向量，并且是指单位面积上的 数值。当假设应力分量盯。是拉应力时，其值为正值。如果是压应力，则需要在 系数向量 ，) 前面乘以“一l ”来进行求值。 北京t 业人学t 学坝i 学位论义 3 7 约束方程的集合表达式 经过以上步骤分析，得到了用极限分析下限理论所表示的各个约束条件矩 阵，剩下的步骤就是将所有限定条件矩阵集合汇总成总体约束条件矩阵，并对目 标函数系数进行详细讨论。 联立等式( 3 5 ) ，( 3 - - - 1 0 ) ，( 3 一1 5 ) ，( 扣之8 ) ，各种约束条件可以用如下 式子统一表示，也就得到了整体约束条件矩阵： 悱耋+ 主川+ 圭k 】+ 砉k 】 c s 卅， 式中，各个系数都是根据行和列来编排进去的，代表节点的总数，e 代表单 元总数，d 代表不连续边( 面) 的总数，工代表要进行目标函数求值的边界面的 总数。 类似地，向量 b 和目标函数系数 c 可以按照下面的整体集合表达式给出： 6 = 耋阮】+ 耋k 嵯k 】+ 毫吃】 斜= ( 3 - 3 5 ) ( 3 - 一3 6 ) 其中s 代表需要优化的正应力的边界的总数。 这样的话，就能把找到一个静态容许应力场的的极限荷载的问题就可以转化 成为求下面方程的解的问题。 最大值： c 7 p 约束条件：阻】 盯 6 其中 盯) 为应力数值，由下面的式子计算出： 盯 7 = 盯1 盯，if 掣i 仃，2 盯，2r 叫2 盯，3o r y 3f 叫3 ( 3 3 7 ) 如果求得的j 下应力为压应力的话，就在目标函数系数向量前面乘以“1 ”， 求解表达式如下： 最小值： 一 c p 约束条件： - 】 盯 ( 6 ： ：】 盯 = 6 ： 其中 4 】为屈服条件矩阵，【a 2 】为其它等式约束矩阵。求解这类问题的线性优化 方法将会在下一章节中具体给出。 3 8 本章小结 按照极限分析有限元法下限解的基本原理，本章深入探讨了该方法的实现过 程，采用有限元法将整个结构进行离散划分，设定每个结点的应力，把原工程问 题变成一个以边坡或基础的极限承载力大小为目标函数、以结点应力分量为优化 变量，以静态容许应力场的各种制约为约束条件的线性优化问题，并借助线性规 划方法对极限分析理论的下限解问题做了详细的阐述，明确了总体约束条件矩阵 的形成过程以及目标函数所代表的涵义。 4 1 概述 第4 章线性规划及其优化方法 根据上一章节的叙述，我们已经将求解边坡稳定( 或地基承载力) 极限荷载 的问题转化成为数学中的求函数最大值的问题，从而可以直接通过数值计算方法 求得实际工程问题的解答。运用有限元理论将研究土体进行离散化，以节点应力 分量为优化变量，求得符合约束条件下目标函数的最大值，归跟结底就是线性规 划【2 7 】( l i n e a rp r o g r a m m i n g ) 的问题。 4 1 1 线性规划的发展 线性规划是运筹学的重要分支，它是一门实用性很强的应用数学学科。这门 学科产生于2 0 世纪3 0 年代，是为了解决二战期间美军后勤运输问题而提出来的。 最早提出线性规划问题的，是前苏联数学家康托洛维奇，他于1 9 3 9 年在生产组 织与计划中的数学方法一书中，最早提出并研究了这一问题。随后，在1 9 4 7 年，美国数学家d a n t a i g ( 丹泽格) 提出了一般的线性规划数学模型和求解线性 规划问题的通用方法单纯形法，为这门学科奠定了理论基础。从此，线性规 划的理论和算法逐步得到了完善与发展，它的应用范围r 益扩大，几乎遍及商业 活动、工业生产等一切经济领域。n 2 0 世纪7 0 年代后期线性规划又取得了重大进 展。1 9 7 9 年，前苏联数学家哈奇扬提出运用求解线性不等式组的椭球法【2 8 l 去求解 线性规划问题，并证明该算法是一个多项式时间算法。这一工作具有重要的理论 意义，但实用效果不佳。进入2 0 世纪8 0 年代后，线性规划的另一种思路完全不同 的解法内点算法【2 9 】得到了迅速发展。1 9 8 4 年，在美国工作的印度数学家卡玛 卡尔( k a r m a r k a r )