（固体力学专业论文）自适应网格加密在弹塑性分析中的应用.pdf

自适应网格加密在弹塑性分析中的应用 摘要 ， 、， 、基于自适应有限元网格加密技术0 本文研究了裂纹扩展的网格生成技术 的算法以及运用应力位移恢复继承技术结合弹塑性自适应有限元跟踪分级加 载条件下塑性区的扩展a 一在四边形自适应网格加密算法的基础上，添加了裂纹生成算法。提出了 通过修改单元的结点信息和单元间的邻接关系来生成裂纹的方法，实现了在 指定裂纹的情况下的网格生成，为进一步研究边坡稳定、地基承载力等问题 打下基础。文中给出了裂纹沿各个方向扩展的算例，算例的结果表明该算法 与有限元计算结合良好，说明了该方法的有效性。 塑性区的确定在岩土工程中有着重要的意义。自适应网格加密算法结合z - z 误差估计方法能有效地跟踪塑性区的扩展，由于自适应网格加密的特点，各级 荷载间的应力位移传递需要运用恢复、继承的手段来实现，构造出最优网格上 的光滑的应力位移益面，以保证下级荷载下弹塑性计算和总应力总位移结果的 正确性。文中采用s p r 方法进行了应力恢复，并设计了有效的算法实现了不同 网格上的应力位移映射，与厚壁圆筒理论解比较表明，该算法精度较高。应力 恢复和继承技术结合网格退化技术大大地提高了计算效率。文中探讨了洞室围 岩、挡土墙等构筑物的塑性区扩展问题，上述的些探讨性的研究结果初步显 示出自适应网格加密方法在工程应用中的潜力 关键词： 一 自适应网格加密：裂纹网格生成；弹塑性有限元；应力恢复、继承；、分级加荷： 塑性区扩展 本文的研究工作获国家自然科学基金( 批准号：5 9 8 7 9 0 1 2 ) 资助。 a d a p t i v em e s hr e f i n e m e n ta n d i t sa p p l i c a t i o n s i ne l a s t 0 一p l a s t i ca n 龇s i s a b s t r a c t b a s e do nt h em e s hr e f i n e m e n tt e c h n i q u eo f a d a p t i v ef i n i t ee l e m e n tm e t h o d ，m a i nw o r k d o n ei nt h i sp a p e ri s ：o n ep a r ti sm e s h g e n e r a t i o nt e c h n i q u eo fc r a c kp r o p a g a t i o n ，t h eo t h e rp a r to f p r o c e d u r ei m p l e m e n t a t i o no fe l a s t o - p l a s t i ca n a l y s i so fs t r u c t u r eu n d e rl o a da d d e db yg r a d ew i t h t h em e a s u r eo fs t r e s sr e c o v e r y a n di n h e r i t i n g m e s hg e n e r a t i o na l g o r i t h mo fc r a c kp r o p a g a t i o ni sa d d e dt o q u a d r i l a t e r a lm e s hr e f i n e m e n t a l g o r i t h mt oc o m p l e t et h ec o m p u t a t i o no fs o m e t h i n gw i t hc r a c k ，s u c ha st h es l i d i n go fb a s eo f a r c h i t e c t u r et h i sa l g o r i t h i nn e e dm o d i f y i n gt h en o d em e s s a g eo fn o d e sa n dt h e n e i g h b o r i n g c o m l e c t i o no fe l e m e n t sw h i c hi su pt ot h em e s hg e n e r a t i o n ，s oi tc h a n g e su n d e rd i f f e r e n tl o a d l e v e lt h ep r o g r a mh a sr e a l i z e dt h ef u n c t i o no fm e s h g e n e r a t i o no fd e s i g n a t e dc r a c k ，w h i c hm a k e t h ef o u n d a t i o no fs o l v i n gp r o b l e ms u c ha ss l o p es t a b i l i t ya n dc a r r y i n gc a p a c i t yo ff o u n d a t i o n c r a c ko fa l ld i r e c t i o na r cg i v e ni nt h i sp a p e rt ov a l i d a t et h ea l g o r i t l m a ，a n di ts h o wt l m ti tw o r k w e l lw h e n c o n n e c t i n gw i t ha f e mt h e n n m e r i c a le x a m p l eo ff o u n d a t i o ns h o w st h ei n i t i a lu s eo f c r a c km e s h g e n e r a t i o n i ti so f g r e a ti m p o r t a n c e o f d e t e r m i n ep l a s t i cz o n e si ng e o t e c l m i c a le n g i n e e r i n ga d a p t i v em e s h g e n e r a t i o nc o m b i n e dw i t hz - z e r r o re s t i m a t i o nc a l lt r a c et h ep l a s t i cz o n ee f f i c i e n t l y b e c a u s eo f t h ec h a r a c t e r i s t i co fa d a p t i v em e s hg e n e r a t i o n ，t h ed a t a ，s u c ha ss t r e s so rd i s p l a c e m e n tc a n t t r a n s f e rd i r e c t l yv i am e s ho fd i f f e r e n tl o a dl e v e l ，s ow eh a v et or e s o r tt ot e c l m i q u eo fs t r e s s r e c o v e r ya n d s t r e s s i n t e r p o l a t i n g t o i m p l e m e n tt h ea l g o r i t h mc o m p a r e dw i t h t h i c k w a l l e d c y l i n d e r st h e o r e t i c a ls o l u t i o n ，t h er e s u l t ss h o wg o o da g r e e m e n t ，w h i c hv a l i d a t et h ea l g o r i t h m t h e nw eu s et h ep r o c e d u r et oc o m p u t eo t h e rc o m p u t a t i o n a lm o d e l st os t u d yt h ed e v e l o p m e n to f t h e i rp l a s t i cz o n e ，s u c ha sw a l lr o c kw i t hc a v i t ya n dr e t a i n i n gw a l lw i t hr e i n f o r c e dc o n c r e t e ，t h e r e s u l t ss h o wt h a tt h em e s h g e n e r a t i o nm e t h o d h a si t sp o t e n t i a li ne n g i n e e r i n g sa p p l i c a t i o n k e yw o r d s ：a d a p t i v em e s hr e f i n e m e n t ；c r a c km e s hg e n e r a t i o n ；e l a s t o - p l a s t i cf e m ； s t r e s sr e c o v e r ya n di n h e r i t i n g ；a d d i n gl o a db yg r a d e ；e x p a n d i n g & p l a s t i cz o n e 上海交通大学硕士学位论文 1 1 研究背景 第一章：前言 有限元方法广泛地应用于工程的各个领域。有限元分析的一个重要环节是形成有限元 网格，有限元网格对于计算结果的可靠性、计算精度有着很大的影响，在很大程度上要依 靠计算人员对其所面临问题的判断来选择合适的有限元网格。对于边界及受荷条件都很复 杂的问题，要形成合适的有限元网格并非易事。其复杂的前处理闯题，也使一般的工程技 术人员难以掌握。这些问题在一定程度阻碍了有限元的应用。 自适应有限元分析方法是解决这一问题的有力工具。它是在初始网格的基础上，根据 计算结果和后误差分析，通过适当增加和减少网格中单元的数量，单元的尺寸相应地发生 变化，从而产生个较佳的网格，经再次计算后判断是否还需增加或减少单元数量。这样 反复循环直到收敛于较精确的解或满足计算精度。该方法的一个特点是，在应力集中的地 方，由于离散误差比较大，网格不断加密，从而提高计算精度而且可以从计算结果输出的 网格图直观地看出应力集中区。 岩土工程中，边坡、地基滑动问题一直是人们关注的问题。目前工程中使用的方法一 般采用极限平衡法 1 】，把土体的破坏看作刚体作滑移运动来研究，而不进行应力分析，是 不符合客观实际的。对于复杂地层和坡顶有复杂荷载的情形，传统的方法很难处理，就得 采用数值分析方法，如有限元方法。传统的有限元，网格一旦生成则不再改变，而分级加 荷下塑性区不断扩大，应力集中区也不断改变，自适应网格加密是解决这问题的有效途径。 当土体含局部破裂面的时候，若有限元网格不反映出这种变化，将会大大地影响结果的正 确性因为土体的滑动是由破裂面的产生、扩展导致的。因此，就提出了在有限元的计算 过程不仅能动态划分网格，还能产生破裂面的要求。土作为弹塑性体，其应力和变形与加 载的路径有很大的依赖性，所以在分级加荷中，有效正确地继承前级荷载下的应力对计算 结果有极其重要的作用。 1 2 研究目的与现状： 本项目的研究目的是用自适应有限元跟踪分级加荷下土体的破坏过程为工程实际提 供可靠的依据。本文将其简化为平面应变问题，以裂纹代替破裂面，用自适应网格加密算 芒、一一一 法中添加裂纹生成算法，为研究滑动问题提供可靠的网格生成技术以及不同网格间的应力 位移的映射算法。 、 上海交通大学硕士学位论文2 对于裂纹跟踪的网格生成问题，国内外学者作了大量的研究，大致可分为两类：一类 采取不离散的方法，对裂纹边的单元进行特殊处理使其矩阵接近于真实情况 2 ；另一类是 通过离散，使裂纹边的边界条件接近于真实【3 】。本文基于工程实际的考虑，采用第二类方 法。文献采用的方法也各不相同，包括无单元法【4 6 】、边界元法及由其发展起来的位移不 连续法【7 】、以及有限元法【8 】，还有以上方法相结合的方法。就自适应有限元网格生成而言， b i t t e n c o u r t 和w a w r z y r e k 9 】提出把裂纹周围的四边形单元划去，然后以三角形单元重新划 分，从而得到光滑的裂纹边。这样，裂纹每扩展一步都要进行一次重新划分部分网格，重 新插值得到三角形单元的应力值。s h c p h a r d 和y e l h a 1 0 1 j l 在每步裂纹扩展后，以裂纹作为 边界的一部分，重新划分整个网格。显然，这两种方法的工作量都较大。这里的裂纹扩展 主要指断裂力学中的i i i 复合型，因为这种裂纹更具一般性。本文基于自适应网格加密的 特点，提出通过修改单元形状和单元间邻接关系的方法来模拟裂纹扩展，具有修改范围小， 可实现性强，无须增加单元类型的优点。 本文的另外一项研究内容是应力的恢复与继承的问题。在用般的有限元求解分级荷 载下的弹塑性计算时，总的应力位移可以用各级荷载的结果直接叠加。用自适应有限元计 算则不一样，不同荷载下的网格不太可能重叠，直接叠加误差很大，特别是进入塑性区以 后，塑性区不断扩展，加密集中带不断移动，两级荷载下的网格更不可能重叠。在进行弹 塑性计算时，需要初始应力来决定采用弹性矩阵还是塑性矩阵，所以应力继承的正确与否 直接关系到非线性自适应分析能否正确实现。在分级加荷条件下，第二级荷载的初始潮格 的选用也是一个关键问题，若直接采用第一级荷载下的最终网格，两层网格间的数据可直 接传递，但是这样给第二级荷载造成很大的计算量，更难保证第三级荷载或更多，造成计 算的理论上效率提高而实际上因为太耗时而降低效率。因此，适当地把前级荷载下最优 网格退化到相对粗的网格作为后级荷载下的初始网格，并在退化后的网格正确继承最优网 格上的数据，将大大地提高计算效率。关于应力恢复和应力继承问题，z i e n k i e w i c z ，z h u 等 作了大量研究 1 1 1 4 】。主要的方法是通过应力恢复在最优网格上由高斯点的应力插值形成 计算区域内光滑的应力场，设计有效的映射算法使要求的网格上的单元从最优网格上插值 得到相应的应力。 上海交通大学硕士学位论文3 第二章：误差估计方法和自适应网格加密方法 误差估计方法和自适应网格加密技术是自适应有限元方法的核心。 有限元离散误差的估计是评估分析精度的依据，也是进行自适应网格加密的基础。许 多学者作了大量而广泛的研究，提出了许多不同的估计方法。z i e n k i e w i e z 和z h u 1 5 1 等a 在误差估计方面作了大量而富有成效的工作。他们提出通过局部应力场的结点插值或最小 一一一 _ 三量苎竺兰竺苎堕堡! 堡至坚垄望竺兰滑! 垄：苎竺宴竺量堇塑垦墨塑堕堕垄堑里自! ! 鎏 ，一一 表示离散误差。该方法也称为z z 法，己被人们所广泛采用。但是该方法是针对线性问题 、一 提出的，不能直接应用于非线性问题。钟【1 6 在此基础上对该方法进行了修正，提出了加 权参数加密准则，该准则能清晰地跟踪出弹塑性交界面。本文的弹塑性研究采用的加密准 则均采用加权参数加密准则。 自适应加密是一种能自动调整算法以改进求解过程的数值方法。目前，应用中常用的 方法有以下几类：( 1 ) h - 加密 1 7 - 1 8 ，此法是不改变插值函数的阶次，仅基于局部网格加 密，提高求解精度。( 2 ) p - 改h 芷f 9 1 ，不改变网格大小布局，仅增加局部插值函数的阶次， 以降低求解误差。( 3 ) 移动结点法【2 0 一2 2 1 ，此法不改变单元的类型和数目，通过改变单元 形状减小离散误差。( 4 ) 组合法：如h - p 法 2 3 1 。本研究采用的是h 加密方法。 2 1 误差估计方法 2 1 1 z - z 法 z z 法采用以能量范数表示的误差估计，并同应力光滑方法相结合。设有限元近似解 位移为疗，应力为子，真实解为”，盯，则误差为： 位移：p = “一万 应力：p 。= 盯一子 ( 21 ) 单元误差定义为： 位移：俐2 = e r e d f l 应力：忙。1 1 2 = l 0 。) 20 。k q ( 22 ) 对整个网格，则有 l l e l l 2 ：氨。盱 ( z ， i = 1 式中，n l 为单元总数。 对于一般的工程实际问题，往往很难找到精确解，对于线弹性问题，可采用插值法， 上海交通大学硕士学位论文4 即用与位移插值相同的形函数进行插值。设结点应力为万，位移为盯，则有 o - = n 厅 又有 p 一子= 0 将( 2 4 ) 代入( 25 ) 得： 孑= a 。1 i n 7 d 凇孑 j n 式中 爿：f r 黼 j o 所以有巳= 盯+ 一子 应力计算的相对误差定义为： 可也| f 2 i 2 若计算误差以7 7 表示，给定的误差以叩。表示，则计算退出的条件为 叩叩。缸 假设误差在各个单元间平均分布，则有 f l e o l l 旷m 州 式中e 。为单元误差最大容许值，定义 茧= l l p 。l l e 一 为加密准则，若董 1 ，则该单元该加密a ( 24 ) ( 25 ) ( 26 ) ( 27 ) ( 28 ) ( 2 9 ) ( 2 1 0 ) ( 21 1 ) ( 2 1 2 ) 2 1 2 加权参数加密准则1 1 6 j 对于弹塑性问题，上述方法则不太适用。z - z 方法是针对线弹性问题而提出的，不能 直接用于弹塑性问题。当用于弹塑性问题时，会出现误差估计不收敛，网格加密紊乱的现 象。同时，z z 方法误差估计方法对弹塑性交界处应力梯度变化的反映不灵敏，不能很明 显地跟踪出弹塑性区交界线。加权参数加密准则对这方面作了补充，其做法：给弹性单元 和塑性单元以不同的权重y ，应力计算的相对误差彳可改为( 1 + y ) 可，以此作为网格加密 细化的标准。y 有如下特点：对于弹性单元，r = o ：对于周围没有弹性单元的塑性单元， 上海交通大学硕士学位论文 而弹塑性交界处的单元，0 y 0 那段四种情况，其余四种可化为这四种情况来求。比如求处 于v 区的裂纹，本文的作法： 上海交通大学硕士学位论文 0 口 口os 口 石2 石2 口 t - t 2 0 7 - g 0 如锄( 3 1 ) 石口 万+ a o t + 口o 口 3 2 z c 3 2 ；r 窿 2 石一口o 2 石一口 2 石 。 j ， 影 签衩7 图3 4 裂纹与x 轴的夹角 f i 9 3 4a n g l eb e t w e e nc r a c ka n da x i sx 1 ) o 点搜索单元至终点。 2 ) 终点搜索单元至o 点。 3 ) 从终点扩展至o 点，扩展的情况与在i 区的单元相同。 4 ) 处理起点终点，以及与已扩展裂纹的衔接问题。 裂纹在搜索过程以及终点处理均会产生新的单元，1 ) 过程与2 ) 过程并不是单纯的相 反过程。 ( 2 ) 终点的处理 由于裂纹的终点不一定会在网格线上所以要对终点进行处理。 2 7 5 l 3 9 8 4 图3 5 终点所在区域示意图 f i 9 3 5d e m o n s t r a t i o no f a r e ai nw h i c h t h ee n d p o i n tl i e s 如图3 5 中虚线所示，把单元分为九个区域。终点落入不同的区域将作不同的处理。 除此以外，还根据裂纹方向的不同的处理，对于每个区域，裂纹来自的方向有八种都得作 出不同的处理。如终点落入1 区域，就细化该单元，并把终点的坐标赋予l 区域中的结点。 并根据裂纹来自的方向判别c r a c k i l l 2 取的值。 上海交通大学硕士学位论文 图36 是从计算结果网格图截取出来的部分网格图。图中终点不在单元的边上，经过 终点处理( 主要是在不破坏自适应有限元的数据结构下调整结点位置) ，得到如图的情况， 这样的处理保证了裂纹的端点在有限元网格的结点上。 不同方向的裂纹需要进行衔接处理，要处理的主要是沿x 轴搜索单元的裂纹和沿y 轴 搜索单元的裂纹的衔接，以及扩展方向处于图3 4 中x 轴的负轴区的裂纹在反向扩展后， 要进行与前后裂纹的衔接，图37 为一种衔接处理的效果。 7 9 0 1 7 7 l 1 7 7 31 7 7 己 1 4 9 41 9 8 7 t 9 8 91 9 8 8 图3 6 终点处理效果图 f i 9 3 6d e m o n s t r a t i o no f r e s u l to fp r o c e s s i n gt h e e n dp o i n t i 滁、 、弋l 挝 1 4 q 3 f 。岁卜 ：十 2 0 0 7 2 0 0 6 1 9 8 6 1 9 8 5 】3 7 92 8 5 4 1 9 5 l 2 8 5 5 图37 两段裂纹的衔接 f i 9 37c o n n e c t i o no f t w os e g m e n to f c r a c k 上海交通大学硕士学位论文 ( 3 ) 单元的分裂处理 单元的分裂处理主要有三个步骤： ( i ) 节点置换 ( i i ) 修改裂纹边邻接关系 ( i i i ) 修改侧面邻接关系 节点置换把裂纹所经单元与裂纹的交点作为该单元的两个节点，也就是把单元的结点 拉到裂纹边上，并赋予该交点两个结点号，分别为c r a c k i 0 】，c r a c k i l 1 ，分别为 裂纹边两个单元的结点。 修改节点信息( 节点信息图如图2 1 ) ： n o d e s ；i n e r t j 3 = c r a c k i 1 ； n o d e s n e l l l 2 = c r a c k i - i j 【ij ： n o d e s 【n e l 儿o 】= c r a c k i 0 b n o d e s ； n e l l i = c r a c k i i - 1 1 1 0 1 ； 要使裂纹边的两个单元脱离关系，使单元边成为自由边，就要使两个单元在裂纹边上 无邻接单元。 修改邻接关系( 单元邻接关系图如图22 ) ： n e l c o n n e l 2 0 = 0 ；令0 方向无邻接单元； n e l c o n n e l i 2 = 0 ；令2 方向无邻接单元。 单元的分裂主要分如图3 6 所示的三种情况： 圉日园圉日圉目日圜 ( b ) ( c 1 图38 单元的分裂( a ) 两端端点在节点上：( b )有一个端点在节点上；( c ) 两端均不在节点 上。 f i 9 38s p l i to f e l e m e n t ( a ) b o t h t w oe n dp o i n t sa r en o d e s ；( b ) o n eo f t h et w oe n dp o i n ti sn o d e ；( c ) n e i t h e r t h et w o e n d p o i n t si sn o d e 裂纹的方向包括八种情况( 如图34 ) ，其开裂方向也分为四种0 ，1 ，2 ，3 方向( 如 图2 2 ) 而且网格是局部加密的，搜索开裂方向的过程比较复杂，但基本思想如上述。 如图3 9 两种情况： ( a ) 图3 9 两种复杂特殊情况( a ) 4 方向有单元：( b ) 5 方向有单元。 f 9 39t w os p e c i a la n dc o m p l i c a t es i t u a t i o n ；( a ) w i t he l e m e n ti n4 “d i r e c t i o n ；( b ) w i t h e l e m e n ti n5 “d i r e c t i o n 网格经过修改以后，与裂纹的形状得到一致，裂纹面的邻接关系也为零。但是就这 样是无法进行有限元计算的，为了保证位移的连续性和与自适应网格数据结构的一致性， 还要修改单元的侧面邻接关系。 78 l 。 一r 6 一 ， 5 34 7 ，一r 8 1 1 ，二，。 6 3 45 图31 0 侧面邻接关系修改前的网格图31 1 侧面邻接关系修改后的网格 f i 9 31 0m e s h b e f o r es i d ec o n n e c t i o n sm o d i f y i n gf i 9 3 1 lm e s ha f t e rs i d ec o n n e c t i o n sm o d i f y i n g 如图31 0 ，若单元l 和单元2 都往方向2 ( 见图2 2 ) 开裂，它们的侧面邻接关系没 有改变。此时单元2 在裂纹的上方，按规定的准则，单元6 与单元8 单元裂开，单元6 将 在裂纹的下方，根据原先的邻接关系单元6 是单元2 在1 方向的邻接单元，这就不对了。 若单元2 往方向0 开裂，单元2 、7 、6 、8 的侧面邻接关系没有改变，此时单元2 在 裂纹的下方，而单元1 在裂纹的上方，单元1 和单元2 的侧面邻接关系又将改变。 所以要改变裂纹所经单元的侧面邻接关系。 为了统一，本程序规定如单元2 这种情况均往方向2 开裂，而在单元3 还未开裂的时 候先预测它的开裂方向，预先修改它的邻接关系，避免在裂开新的单元后还要返回修改已 经过单元的邻接关系，因为开裂过程还是在不断地产生新的单元。 图31 1 为修改邻接关系后的网格。 在修改邻接关系的过程中，可能产生新的单元，原来搜索的裂纹所经单元c r a c k i l l 2 】 被细化掉，所以裂开过程中不断地调整c r a c k i 2 ，举个例子，如图31 2 ： 上海交通大学硕士学位论文 ( a ) 开裂前 n e l ， 一 ( c ) 把n e l 细化 裂纹 每一。 r 江 ( b ) 开裂后( 错误的结果) 土一 图3 1 2 一个开裂过程 f i 9 3 1 2o n ep r o c e s so f c r a c k i n g ( d ) 开裂后( 理想的结果) 如图( a ) 为这样一种情况，裂纹扩展的方向如图3 4 的i 区，小箭头表示n e l 将往下 方开裂，该单元的特殊是5 方向有单元，若按前面所述的一般原则，将会得到如图( b ) 的效果，点2 被移到超过点2 的裂纹边上，这样就破坏了单元a 的单元形态，导致有限元 无法正常实现。对这种情况得先细化单元n e l ，但细化以后裂纹所经单元c r a c k i 2 】 不再是原来的单元，即使单元号没改变单元大小已变化，原先通过n e l 的裂纹现在通过的 是两个单元，本程序的处理是增加一个c r a c k i 2 ，并把裂纹所经单元的总数i c r a c k 加l ，原先的c r a c k数组中c r a c k i + i 2 到c r a c k i c r a c k 2 往后移一位，同时 l 1 n e x c o 数组也得做同样的处理。单元的细化还可能导致前一个开裂单元根据n e l 的侧 面邻接关系的修改已经变化，也得作处理。总之，任何一个环节没考虑到，都会使网格出 现紊乱，导致有限元计算无法正常实现。 还得说明一点，由于裂纹生成过程中不断地产生新的单元和结点，正确地继承或插值 得到新结点的位移和新单元的应力与单元信息如泊松比，弹性模量也都是关键，继承或插 值将在下章中谈到。 上海交通大学硕士学位论文1 6 一圭堂奎望奎兰堡主兰垡堕壅 ! 1 3 2 算例 3 2 1 网格生成算例1 本算例主要表现裂纹沿各个方向( 如图3 4 中i 、i i 、v i i 、v i i i 的情况) 的扩展， 计算的对象是带预裂纹受单向拉力的板，计算模型如图3 1 4 ，采用米塞斯屈服准则，取盯 为1 0 助口，在荷载不变的情况下，让裂纹扩展了四个阶段。图3 1 5 ( a ) ( d ) 均为网格 图( 图中标数字的单元表示该单元已屈服，各图均为部分网格图) 。图3 1 5 ( a ) 表示裂纹 方向图3 4 中i i 情况，这里取式3 1 中= 万4 ，3 1 5 ( b ) 可以看出，裂纹扩展后网格 在新的裂纹尖端附近形成新的塑性区，裂纹边的单元形状经过调接，但没有破坏自适应加 密算法中的相邻单元的层次差不能超过i 的准则，网格加密有序进行，塑性区随着裂纹尖 端的移动而移动。3 1 5 ( b ) ( c ) ( d ) 分别为图3 4 中i 、v i i i 、v i i 情况，每次扩展都在 前次扩展的有限元计算满足精度后进行。3 1 6 ( a ) ( d ) 为相应四个阶段的位移图( 位移 图中的结点经过协调，图为全图) 。位移图的结果表明经过开裂，裂纹边成为自由边，裂纹 间的衔接正常，位移的结果也比较理想( 张开的开口越来越大) ，说明了裂纹两边的单元的 邻接关系以及裂纹周围单元间的邻接关系的处理达到预期的目的，能正常地进行自适应有 限元的计算和网格加密，说明了该裂纹生成算法是有效的。 ttt tt t tt 。 争 o m 、 ，0 上上上盯上0 上 上、 、 l o m 7 图3 1 4 带预裂纹的板的 计算模型 f i 9 3 1 4 c o m p u t a t i o n a l m o d e l o f p l a n ew i 也c r a c k 上海交通大学硕士学位论文1 8 ( d ) 第四步网格图 ( 3 必口 2 厅一) 第一步网格图( 口x 2 ) ( c ) 第三步网格图 ( 2 石- 5 0 口 2 7 r ) 图3 1 0 裂纹扩展过程网格图 f i 9 3 1 0m e s h o f t h ef o u r s t e po f c r a c kp r o p a g a t i o n 上海交通大学硕士学位论文 c a ) 第一步位移 ( c ) 第三步位移图 ( b ) 第二步位移图 ( d ) 第四步位移图 图3 1 6 裂纹扩展的四步位移图 f i 9 31 6d i s p l a c e m e n t o ft h ef o u rs t e po f c r a c kp r o p a g a t i o n 3 2 2 网格生成算例2 本算例研究挡士墙受墙背土压力的作用，以及墙脚地面土的释放压力裂纹扩展的效果。 计算模型及尺寸( 单位为米) 如图3 1 9 ，p l = t 2 0 匈a ，p 2 = 2 0 印| 混凝土的弹性模量 25 1 0 5 m p a 土的弹性模量为2 0 0 0 0 0 k p a 。 上海交通大学硕士学位论文 图3 1 7 挡土墙计算模型 f i 9 3 1 7 c o m p u t a t i o n a l m o d e l o f r e t a i n i n gw a l l 5 0 采用米塞斯屈服准则，经过计算，墙背将会出现塑性区，现人为地让墙背的土扩展一条 长度为3 m 的处于区的裂纹，观察裂纹扩展前后墙体及土体的受力状态，并验证裂纹处 于区的程序实现效果。 l l _ i _ l - - l li l ll i j _ l i i _ l lili l llll ( a ) 开裂前的网格图 l jl l +l l l l l _ j l l - l l 一 - ll | l j _- - ll l l - j - l - - l l l1 l ll ( b ) 开裂后的网格图 图31 8 墙背土开裂前后的挡土墙网格图 f i 9 3 1 8 m e s ho f r e t a i n i n g w a l lb e f o r ea n d a f t e r c r a c k 上海交通大学硕士学位论文 2 l 1 0 1 5 0 02 0 0 02 5 3 0 0 03 5 0 04 0 0 0 ( a ) 开裂前的位移图 诵厶l 国 一目 l1 人l l 户斗芷料牛卜r 广liiii1 l 窜l 稻瑶簿 l 、j 1 i l ( b ) 开裂后的位移图 图3 1 9 墙背土开裂前后的挡土墙位移图 f i 9 31 9d i s p l a c e m e n to fr e t a i n i n gw a l lb e f o r ea n da f t e r c r a c k 从图31 8 和31 9 看出，裂纹方向处于区的算法得到正确地实现，裂纹扩展后，墙体 的位移扩大( 位移图的位移经过放大1 0 0 倍) ，墙体往外倾，说明了该算法与有限元结合 良好。实际上，土是松散介质不能承受拉应力，所以土的破裂面法向的位移不会很大，主 要是滑裂面的切向位移，所以完全地让破裂面成为自由面和实际也不相符，这些将有待于 进一步研究来提出更合理的模型，本章的目的是提供有效的裂纹生成算法，这个目的已得 到了实现。 上海交通大学硕士学位论文 第四章：应力、位移恢复和继承技术的算法 为了更好地模拟岩土的性质，本章采用弹塑性模型来模拟分级加荷下塑性区的扩展。弹 塑性问题并不是弹性问题的简单延伸，除了模型，误差估计方法要修改以外，还提出了网 格应力继承，网格退化的问题。在一般有限元，网格确定以后不再改变，对同单元，在不 同荷载级别下的位移应力都可叠加以获得总的位移和应力。而对自适应有限元，不同荷载 级别下的网格不太可能重合，即使同一荷载级别下不同层次的网格也不同，特别是计算带 裂纹的结构时，网格的加密区域随裂纹尖端变化而变化。计算弹塑性有限元时，在第二级 或以上级荷载下的初始应力应为前级荷载下的应力，且在每层网格开始计算时都要继承前 级应力为初始应力，所以必须有效地继承前级荷载下的应力。 4 1 弹塑性模型 弹塑性模型把总的变形分为弹性与塑性变形两部分。对于塑性变形主要作三方面的假定： ( 1 ) 破坏和屈服准则 ( 2 ) 硬化规律 ( 3 ) 流动法则 对于不同的模型，这三个假定的内容不同。 ( 1 ) 破坏准则 土体达到破坏后，变形会不断发展，与破坏前是截然不同的，判断破坏的标准即破坏 准则。因土体的破坏取决于应力状态，散破坏准则可写为： + ( 盯。) = k f ( 4 1 ) 其中( 盯。) 是破坏函数，k f 是试验确定的常数。若，( ) = k f ，则破坏； ，+ ( 盯。) k f ，不破坏；，( ) 不可能超过k f 。 主要的破坏准则有下面几种： ( i ) 屈雷斯卡( t r e s c a ) 准则 这一准则假定最大剪应力达到某一数值时破坏，即 ! ! 二! ! ：k ， ( 42 ) 2 。 它在主应力空间内是以空间对角线为中心轴的正六角形柱面。 ( i i ) 米塞斯( m i s e s ) 准则 上海交通大学硕士学位论文 该准则假定偏应力达到某一数值时破坏，即 9 2 击厄i 阳i 万丽咆 ( 43 ) 它在主应力空i b 3 内为一圆柱面。 ( i i ) 摩尔一库仑准则 对于土体，摩尔一库仑理论受到广泛应用，其破坏准则为： 学= 半s i n 州唧 ( 4 4 ) 或 旰嘴2 ( 。s 。+ 守z c 喀( 。s 。+ 詈) s ， ( 2 ) 屈服准则 材料的受力变形，在应力较小的时候往往是线弹性的，一旦应力超过某种限度，则开始塑 性变形，口占关系呈非线性，这时称材料屈服了。屈服前后应力应变关系很不相同， 屈服主要取决于应力状态，用公式表达为： b 。) = 七 ( 46 ) 式中，b fj 叫屈服函数，| i 是与应力历史有关的常数。若，b ) 七，材料处于弹性变形 阶段，在应力空间内相应的点落在屈服面以内；当，b g ) = 后，材料屈服。，b 。) 有可能 超过原有的后，但此时七提高了，仍然保持，k ) - - k 。这是由于屈服后材料硬化的缘故。 对于处于屈服状态的体系，即当前应力处于屈服面上，施加一新的应力增量咖f ，就有如 下几种可能，如图4 1 ： 【i ) 慨的方向指向屈服面内部。这将使妙2 瓦o f 峨 。，也就是他f + 拈，) 0 。同时，露增加，发生塑性变形。 上海交通大学硕士学位论文 ( 3 ) 硬化规律 当材料达到屈服后，屈服的标准要改变，即式中七要变化。后随什么因素而变，如何变化， 即为硬化规律。k 的变化有三种：l 、屈服后k 增加，叫硬化：2 、屈服后k 减小，叫软化 t 不变，叫理想塑性变形。 一q 一! 峨9 瑙吾 ( 4 7 ) = 幽 s ， 啪一嚣鬻 h 式中，4 f 别7f 劐0 0 - j l 是反映硬化特征的一个变量，与硬化参数h 的选择有 关。 4 2 非线性计算过程 常用的非线性有限元分析有迭代法和增量法两种。迭代法是用修正劲度的方法( 变劲 度法) ，或保持劲度不变而调整荷载的方法( 常劲度法) ，重复试算逐步逼近真实解。迭代 法包括割线迭代法、初应力迭代法、初应变迭代法等。割线迭代法首先要有全量的本构关 系式，还要假定变形与应力路径无关，因此般不用割线迭代法。初应力存在收敛速度慢， 甚至不收敛的缺点。初应变迭代法比初应力法能更快地逼近非线性的应力应变关系，但当 应力水平较高时，应力应变关系曲线平缓。较小的应力变化会产生很大的应变产生，使初 应变计算误差过大，土体分析中很少采用。为了反映加载历史，了解结构弹塑性变化过程， 跟踪塑性区的发生、发展直至破坏的渐变过程，本文采用增量法作非线性有限元分析。 对于基本增量法和中点增量法都作了尝试，对于弹塑性计算，中点增量法的结果与理 论解更为接近，同时也增大了计算量。基本增量法能清晰地跟踪出弹塑性交界面，而中点 增量法只能跟踪出弹塑性交界带。 4 3 位移、应力继承算法 位移应力的继承问题在自适应有限元中有两种：1 、同级荷载下不同网格间的位移应 力继承；2 、不同级荷载下的位移应力继承。 同级荷载下的应力和位移继承主要是指细化单元时新单元从其父单元继承当前的应 力，位移由其父结点的位移插值而得，这种方法误差很小，因为当前的应力还要经过有限 元计算会被更精确的值代替，如果未被新的值代替说明已达到一定的精度。 但是，如果前一级荷载的应力不能正确地继承，将会影响当前级的有限元计算，进而 又影响下一级，误差不断被放大，所以位移、应力继承算法的重点是不同荷载问的位移、 应力继承。 4 3 1 有限元超收敛理论 2 7 1 有限元的应力佳点定理指出，位移的导数在应力佳点处超收敛且应力佳点集与高斯 分点一致。本文采用的是四边形一次单元，单元的中心点即为应力佳点( 同时也是高斯点) 。 4 3 2 网格映射 作分级加荷计算时，要将前级荷载下满足离散误差的那层网格的应力，位移以及一些 相关变量映射到当前网格，作为初始应力，如式( 41 0 ) ( 41 1 ) 表示，下标。表示旧网格， 上海交通大学硕士学位论文 即前级荷载下满足离散误差的那层网格的应力，下标。表示新网格，即当前荷载下当前层 下的网格。其中( x ，y ) 为当前网格中单元的收敛点的坐标，对位移而言，指的是节点坐标， 对应力而言，指的是高斯点的坐标。 儿0 ，y ) 甘以k y ) ( 41 0 ) 仃。0 ，) 盯。0 ，) ( 4 1 1 ) 由于网格的加密和退化，新旧网格的不可能重合，但是在程序采用将退化的结点置为自 由结点，结点号及结点的位移坐标都保留，以及新产生的节点位移由其父单元节点位移插 值的方法，所以新旧网格的节点的位移是一致的，对同节点号，直接继承即可。 但单元就不一样，无论细化还是退化都无法确定新旧单元应力的关系，比如：图4 3 表 示的是退化过程，左图表示的退化前的单元，右图表示的是退化后的单元，图中加粗斜体 表示屈服或破坏的单元，退化前的一组单元部分屈服，退化后的单元是否屈服，要重新插 值该单元的应力才能确定。所以应力无法直接从同号的单元继承而得，同理，细化过程也 是如此。 j? 43 图4 3 单元退化过程 f i 9 43d e g e n e r a t i o no f e l e m e n t 应力继承过程总结起来有三步： 1 求新网格中单元的高斯点坐标 2 旧网格中高斯点应力外插得节点应力( 应力恢复) 。 3 由旧网格中的节点应力插值得到新网格中高斯点应力。 4 3 3 应力恢复 这三步过程应力恢复是关键。国内外不少学者对这方面做了研究，应力恢复方法主要 分为两种，s p r ( s u p e r - c o n v e r g e n tp a t c h e sr e c o v e r y ) 1 3 1 4 1 和r e p ( r e c o v e r y b y e q u i l i b r i u m i n p a t c h e s ) 2 8 。这里采用的是s p r 方法，下面简要介绍s p r 方法。 如图4 4 ： 上海交通大学硕士学位论文 表示收敛点表示高斯点。表示节点口表示拼接点 有限元解图4 4 一维问题的应力和位移的有限元与精确解 精确解? g 4 4 f e m 竺d e x a c 8 0 1 “石。“o fs 仃e 8 5 锄dd i 8 p l a