（采矿工程专业论文）基于ansys的基坑开挖位移模拟及位移反分析的研究.pdf

内蒙占科技大学硕士学位论文 摘要 目前，基坑工程的围护结构的变形及位移主要通过数值计算的方式得到。合适的计 算参数选取是影响数值计算成果的决定| 生因素之一。通过位移反分析得到有限元计算的 合适参数是行之有效的方法。深大基坑场地内的多层土体的弹塑性位移反分析由于待反 演参数众多，现有的反演方法需要的计算时间很长。 本文提出一种g a b p 位移反分析法，并编制了反分析程序g e n e - - m e b p 。g a b p 算法的主要原理是：通过使用a n s y s 软件进行数值模拟为b p 神经网络提供训练 样本，并用训练好的神经网络作遗传算法的适应度函数求解器，利用遗传算法搜索到合 适的计算参数。本文主要研究了d r u c k p r a g e r 模型的弹塑性参数位移反分析。 a n s y s 模拟基坑丌挖及支护时围护桩体的位移采用水土分算的总应力法计算。土体非 线性分析采用大变形分析。在取到合适的计算参数时。a n s y s 模拟能正确的反映基坑 开挖的变形规律。g e n e - - m e b p 程序经算例验证，收敛较侠且结果正确。 最后用深圳罗湖地铁站基坑工程的位移实测值验证了g a b p 位移反分析法的正确 性和g e n e - - m e b p 程序的有效性。本法的提出为合理利用基坑位移监测数据提供了一 个新的思路。 关键词：基坑变形弹塑性位移反分析a n s y s 软件b p 神经网络遗传算法 内蒙古科技人学硕士学位论文 r e s e a r c h i n g o n d i s p l a c e m e n t s o f s i m u l a t i n g o fe x c a v a t i o no fp i to h a n s y ss o f t w a r ea n d d i s p l a c e m e n t o fi n v e r s e a n a l y s i s a b s t t a c t i nt h i sp a p e r , an e w a p p r o a c h ，ag e n e t i ca l g o r i t h m sc o m b i n e d b a c k p r o p a g a t i o n n e u r a l n e t w o r ki sp u tf o r w a r dt os o l v ep r o b l e mo f ai n v e r s ea n a l y s i so f d i s p l a c e m e n t sf o r e l a s t i c - p l a s t i c m o d e l o f m u l t i l a y e rs o i lo f p ra n d m a k ea p r o g r a mo f g e n e m e b p f i r s t l y ，u s i n ga n s y s s o f t w a r ec a l c u l a t e sd i s p l a c eo f p i l e so fp i tu s i n gd i f f e r e n tp a r a m e t e ro f e l a s t i c - p l a s t i cm o d e lo f s o i li no r d e rt op r o v i d e l e a n i n gs a m # e t ob p n e t w o r k s e c o n d l y ，u s i n g _ f r a i n e d b pn e t w o r kt o b eas o l v e r o f e v a l u a l i n g f u n c t i o n o f g a t h i r d l y ，u s i n g g as e a r c h e so p t i m i z a t i o n p a r a m e t e r o f s o i lp i ta n dv a l i d a t e st h i si n v e r s em e t h o db yp r a c t i c a l p i te x c a v a t i o ne n g i n e e r i n g i nf e m a n a l y s i s ，d r u c k - p r a g e r m o d e l i sa p p l i e d t o m o d e l o f s o i l p i t b e c a u s e o f c o n v e n i e n c e a n d p r a c t i c a b i l i t yo f d p m o d e l ，n o n l i n e a r a n a l y s i so f s o i li sa p p l i e dl a r g ed e f o r m a t i o n c a l c u l a t e d a e x a m p l e s ，s t i m u l a t i n g o fe x c a v a t i o no fp i tp r o v e sa c c u r a t ea n dc o n v e n i e n t u s i n ga n s y s s o f t w a r e f i n a l l y , t h ed i s p l a c e m e n t si n v e r s ea n a l y s i sm e t h o d o f g a - b pi sa p p l i e dt oa n a l y s i sa p i to f s h e n z h e nl u o h u s u b w a y s t a t i o n i ti sp r o v e dt h a tt h em e t h o do f g a - b pi se f f e c t i v ea n d v a l i d i t yo f p r o g r a mo f o e n e - m e b p i si n d u b i t a b l e t h ec o m b i n e dm e t h o d p r o v i d e s an e wi d e a t ou t i l i z eo b s e r v a t i o nd a t ao f d i s p l a c e m e n t si n p i te x c a v a t i o ne n g i n e e r i n g k e y w o r d s ：d e f o r m a t i o n o f p i td i s p 】a c e m e n t so f d a s t i e - p l a s t i ci n v e r s ea n a y s i s a n s y ss o f t w a r eb pn e u r a ln e t w o r kg e n e t i c a l g o r i t h r n s 一2 独创性说明 本人郑重声明：所呈交的论文是我个人在导师指导下进行的研究工 作及取得研究成果。尽我所知，除了文中特别加以标注和致谢的地方 外，论文中不包含其他人已经发表或撰写的研究成果，也不包含为获得 内蒙古科技大学或其他教育机构的学位或证书所使用过的材料。与我一 同工作的同志对本研究所做的任何贡献均已在论文中做了明确的说明并 表1 了谢意。 签名：毛伏 日期：竺! 整旦兰竖 关于论文使用授权的说明 本人完全了解内蒙古科技大学有关保留、使用学位论文的规定， 即：学校有权保留送交论文的复印件，允许论文被查阅和借阅；学校可 以公布论文的全部或部分内容，可以采用影印、缩印或其他复制手段保 存论文。 ( 保密的论文在解密后应遵循此规定) 签名：z 蓉二导师签名：目期：枷妒；f 野巧嬗签名：移务一导师签名：目期：枷y 宇。毋锄趔 内蒙古科技人学硕士学位论文 1 绪论 1 1 引言 随着城市现代化建设的不断发展，有限的城市地面空间已不能满足人们日益增长的 生活和工作需要，于是人们开始向高空和地下寻求发展空间。各类地下工程诸如越江隧 道、地铁车站和区间隧道等己到处可见。国外著名的地下工程有法国巴黎的中央商场， 美国明尼苏达大学土木与采矿工程系的办公大楼和实验室，日本东京八重洲的地下街 等，这些工程的共同特点之一是都需进行大规模地下开挖。以上海市为例，商层建筑 基坑的最大平面尺寸己达2 7 4 m 1 8 7 m ，面积约5 1 0 0 m 2 ，深度最深达3 2 m 。目前基坑 工程具有基坑开挖深度和面积加大，支护结构布置方式更加多样，对支护结构变形控制 更加严格等特点。因此研究基坑开挖后，土体和支护结构的准确的变形量能指导基坑工 程的设计施工，提高基坑工程设计和施工的安全性，通过优化设计可取得较好经济效 益。 1 2 基坑工程的研究现状 基坑工程涉及到土层性质、挡土结构内力和变形的计算理论、基坑的开挖与稳定 性变形机理及计算以及计算参数的反馈分析确定等诸多方面。本文重点研究基坑开挖时 土体和支护结构变形的数值模拟计算、用优化方法反分析合理计算参数。这些课题的研 究现状叙述如下： 1 2 1 基坑开挖时土体和支护结构内力和变形的数值模拟计算方面 应用于土体基坑过程的数值模拟计算方法主要有有限元法、边界元法、无限元法、 拉格朗日元法、半解析元法等。有限元法在土工问题上应用最多。3 土工有限元分析可分为总应力有限元法和有效应力有限元法。有效应力有限元法可 考虑土的固结问题及地下水渗流问题，主要用于开挖时间较长的并受渗流和固结影响很 大的软土深大基坑工程复杂问题中。1 高俊合研究了土与结构相互作用、土体固结的基坑开挖二维有限元问题。“1 平扬等 考虑了基坑开挖时渗流一应力之间的耦合效应，运用了比奥固结理论的有限单元法扩展 到弹塑性领域。“应宏伟等进行了软土深基坑的土体固结的二维有限元分析。唧俞洪良 等用有限元研究了深基坑的渗流特陛对基坑工程的影响。”1 黄春娥等用有限元研究了渗 流对基坑稳定性的影响。“陈洋做了软土基坑的三维流变有限元分析。嘲蒋明镜研究了 平面大变形b i o t 固结有限元程序。”“ 内蒙古科技大学硕士学位论文 更多的基坑开挖问题是用总应力有限元法分析。由于大多数基坑施工时间并不长， t 叮以不考虑土体的固结作用和地下水渗流作用。总应力法用于基坑工程巾一般可以满足 工程精度的需要。 赵利益等进行了基坑开挖三维有限元的弹塑性分析。土体采用德卢克一普拉格弹塑 性模型。”邱弱等进行了基坑开挖锚杆和土钉支护的三维有限元分析，设定了一种双弹 簧元接触元，进行了非线弹性的有限元分析。“詹必秀等进行了基坑开挖的土钉的三维 有限元分析。考虑了莫尔一库仑弹塑性问题。“4 李斯海进行了基坑开挖及土锚支护的二 维有限元分析。设定了古德曼单元，进行了非线弹性的有限元分析。“”宋二祥等进行了 复合土钉的有限元分析。“”刘学增等研究了接触元初始刚度系数对深基坑变形影响。“” 王广国对深基坑的有限元大变形问题进行了研究。“ 1 2 2 优化方法反分析计算参数的方面 周瑞忠对基坑开挖土体进行了b p 人工神经网络的弹塑性位移反分析。采用3 层网 络，土体弹塑性模型采用剑桥模型。”冯夏庭用遗传算法优化人工神经网络的结构，并 进行了弹性材料位移反分析。“”高广运等用遗传算法搜索b p 网络的权值，用于预测基 坑变形。”“ 邓建辉等用b p n g a 法对结构9 个参数进行了位移反分析，用b p n 法代替有限元用 于遗传算法优化的适应度函数的计算。用训练好的b p 网络代替有限元位移的计算，用 遗传算法进行岩石边坡的多层弹性参数反分析。“1 王登刚等用遗传算法进行了地下隧道 的多层介质的弹性位移反分析。2 1 高玮等用模糊数学结合有限元对地下岩体进行了弹塑 粘性反分析。高玮等改进了遗传算法中的交叉算子的算法，并用于岩石弹粘塑性介质反 分析。”王志亮等用遗传算法结合b p 网络对土石坝进行了三维非线弹性位移反分析。” 肖专文提出了g a a n n 法以优化复杂工程问题。 l _ 3 目前国际流行的数值模拟软件在基坑工程中的应用 目前通用的有限元数值模拟软件主要有a n s y s , a d i n a ，a b u q u s e ，m a r c ，s a p 。 其中a n s y s 软件由于前处理和后处理方便，使用简单，而且留有用户与软件的接口， 用户可编制自己程序与软件连接，可定义自己的土体本构模型，破坏准则等而广泛的用 于结构分析。 2 6 1 a n s y s 用于基坑工程的有限元数值模拟比较适合。 a n s y s 有限元分析软件是一个利用计算机技术进行工程分析的大型、通用软件， 可运行于巨型计算机、小型计算机、工作站以及p c 机上。a n s y s 程序能满足从土 木、汽车、电子到宇航、化学等大多数工业领域的有限元分析( f e a ) 的需要，其功能已 2 内蒙占科技大学硕十学位论文 为全世界所公认。a n s y s 软件在1 9 9 5 年成为分析设计领域中第一个，也是迄今为止唯 一一个通过i s 0 9 0 0 1 质量认证的软件，是美国机械工程师协会( a s m e ) 、美国核安全局 ( n q a ) 认证的标准软件之一。a n s y s 软件在中国率先通过全国压力容器标准化技术委 员会的严格考核，成为唯一与中国压力容器分析设计标准0 8 4 7 3 3 2 9 5 ) 相适应的f e a 软 件。 2 7 】 1 4 本课题的提出及主要内容 综上可知，用数值模拟方法可得到基坑开挖的变形。其对于基坑工程的合理设计和 安全施工都是非常有意义的。使用成熟软件进行基坑工程开挖及支护的数值模拟是方便 准确可行的。由于基坑工程土体的复杂性数值模拟时很难选准计算参数。利用工程现场 的监测数据反分析适合于选定的本构关系的参数以解决数值计算等效计算参数不好得到 的问题。 本文主要研究以下两个问题： ( 1 ) 使用数值模拟软件模拟基坑工程的开挖及支护。以获得基坑开挖时的土体和支护结 构的变形量。使用a n s y s 软件作为数值模拟工具，模拟基坑开挖的二维问题。以获得 基坑丌挖时的变形和位移。 f 2 ) 编制遗传算法结合b p 神经网络反分析的程序，利用基坑开挖的实测位移数据反分 析土体的弹塑性参数。研究遗传算法结合b p 神经网络位移反分析法用于多个土层的基 坑开挖弹塑性位移反分析。用数值计算的结果训练人工神经网络，用实测数据检验反分 析的正确性。 3 内蒙古科技大学硕十学位论文 2a n s y s 软件模拟基坑工程中的开挖和支护 由于基坑开挖较怏，傲不考虑基坑_ 丌挖时的固结和地下水渗流的影响。以下的有 限元分析都是采用总应力有限元法进行的。若要进行有效应力有限元法的计算，可以使 用专门的土体有限元软件如c r i s p 进行分析或对a n s y s 进行二次丌发，编制新单元。 2 1 土体弹塑性力学的常用的土体本构模型口8 】 土体的弹塑性本构关系很多，考虑土体的有效应力的本构模型主要有布尔兰特等提 出的修正剑桥模型，沈珠江等提出的南水模型，殷宗则等提出的椭圆一抛物线双屈服面 模型，郑颖人等提出的后工模型。 研究土体的总应力有限元法常用弹塑性模型有d r u c k e r - p r a g e r 弹塑性本构模型。 f 面分别简单叙述。 2 1 1d r u c k e r p r a g e r 弹塑性本构模型的简介 在做基坑开挖及支护的弹塑生有限元分析时，一般不考虑土体固结和渗流作用。总 应力有限元法的分析精度一般可以满足工程要求的精度。进行弹塑性总应力分析时一般 选用d r u c k e r - p r a g e r 硬化模型作为数值模拟的土体的本构模型。简单介绍如下： 屈服准则： f = n l 一万+ 叫r 3 ( 3 + s i n 2 ) i 这里，c 为土的内聚力 庐为土的内摩擦角 ，为第一主应力不变量 j ：为第二篇偏应力不变量 或f = a i l 一属+ k = q 这里，口= s i n e ， r 3 ( 3 + s i n 2 庐) i k = 4 3 c c o s 4 3 c c o s ( 3 + s i n 2 庐) 2 加载条件为： 庐，以) = 吡一万+ 一矿7 = o d = 0 ( 2 1 ) ( 2 2 ) ( 2 3 ) ( 2 4 ) 堕鍪重型垫查堂堡主兰堡堡塞 这罩，h 。= w p = i d h d ； 破坏准则为：m o h r - c o u l o m b 准则，即 f ( 盯i c 乙) 2 一【2 c c o s 十( 口l + ( _ ) s i n 妒 2 ) ( d i 一l ) 2 一 2 c c o s 矿+ ( 0 2 + o - 3 ) s i n 庐 2 ( c r l c r 2 r - 2 c c o s + ( c r l + 盯2 ) s i n 4 2 = 0 增量形式的弹塑性本构关系表达式为： 拈 = d 。a ( d e 式中d e p 为弹塑性矩阵，可表达为： 0 _ d 卜 d 。】 【d _ d 1 。 d 2 d 3 0 0 0 d l3 000 d 2 3 000 d 3 3 000 0 0 0 00 d 5 5 0 000 d 6 6 d t l = d 2 2 = 皿3 = + 2 g d 1 2 = d i3 = d 2 3 = d 4 4 = d 5 5 = d 6 6 = g q = 石i h ：对 h l h 、h ； h 3 h 1h 3 h lh i h 拉、h | h lh h 3 h 5 h l h 5 h 2h 5 h 3 h 6 h 、hb h 2 h 6 h ， 称 h j h s h4h j h 6 h4 h 6 h5 c = a + 3 a 2 ( 3 丑+ 2 g ) + g a = 。c i 。一甄 h t = 0 t ( 3 1 + 2 g ) + g ( t x x 0 0 | 皈 h l = c t ( 3 l + 2 g ) + g ( c r y 0 3 | 瓯 马= 口( 3 五+ 2 g ) + g ( 吒吒) 历 一5 佗5 ) ( 2 6 ) ( 2 7 ) ( 2 8 ) ( 2 9 ) 日b q 0 o o 内蒙占科技大学硕士学位论文 h 4 = g rx y j 2 h 5 = g r f 4 j 2 h s = g rx ：j 2 d r u c k e r - p r a g e r 弹塑性模型的参数较少，有e 、c 、庐、 e 为杨氏弹性模量 “为泊松比 c 为土的内聚力 为土的内摩擦角 2 1 2 修正的剑桥模型简介 2 9 1 修正的剑桥模型是把体积屈服面作为屈服面，该模型的体积屈服面在p 一q 坐标中 为一个椭圆。该模型有五个参数，即e 、五、k 、m ( 或7 ) 。除了弹性参数外 的三个参数都可以利用常规三轴试验测定。五、足值可以利用不同的盯，等向压缩试验 给出。m 值可以通过三轴排水剪或不排水剪试验得到。这种模型能较好的适用于正常固 结粘土和弱固结粘土。 屈服准则： p ，( 等半卜 式中，m 为p q 平面上破坏线的斜率。 m ：! ! ! 璺生 3 一s i n 击 为土体的有效内摩擦角 剑桥模型的弹塑性矩阵为 土业1 a t 2 r - 南 6 ， ( 2 1 0 ) ( 2 1 1 ) ( 2 1 2 ) 一)一矿 行- i + z一2 一m 一 ， 二y p 竺 = = 、lrj 节 嘶万 一 r，j、l d 式 内蒙占科技大学硕十学位论文 2 1 3 其他多屈服面模型简介 南京水利科学研究院双屈服面模型( 简称“南水模型”) 是采用剪切屈服面和体积 屈服面的双屈服面的模型。体积屈服曲线为一组蛋形线，剪切屈服曲线为一组双曲线。 南水模型适用于软粘土。除弹性参数外，有血个计算参数，a 、b 、五、d 、n 。 拉德一邓肯模型是采用剪切屈服面和体积屈服面的双屈服面的模型。共有1 4 个参 数。拉德一邓肯模型是适合砂土的弹塑性模型。 2 2 土体非线胜问题的分析及其分析方法 2 2 ，1 土体的非线性问题 土体的非线性包括材料非线性和几何非线性。其中材料非线性问题是主要问题。材 料非线性中，最常见的是与加载历史无关的非线性弹性和与加载历史有关的非线性弹塑 性。非线性弹性和弹塑| 生都与时间无关。 非线性弹性问题中，应力一应变关系为非线性，但材料为完全弹性的，应力、应变 互为单值函数，与加载历史无关。 非线性弹塑性问题中，应力一应变关系为非线性的，但材料不是完全弹性的，应 力、应变不是互为单值函数关系。而是与加载历史有关。 非缉l 生中还有粘塑性、蠕变等，它们均与时间有关。 本文主要分析材料非线性中的弹塑性问题。 2 2 2 非线性有限元分析方法【3 0 l 当材料的应力一应变关系是非线性时，其刚度矩阵不是常数，而与应变和变形值有 关，可记为i 足( 巧) l 。这时，结构的整体平衡方程是如下的非线性方程组： 渺 = k ( j ) 】 p j = 0 ( 2 1 3 ) 非线性方程组的解法大致有三类，即全量法、增量法和混合法。全量法又称为迭代 法，在每次迭代过程中都施加全部载荷，逐步修改位移和应变，使之满足非线性的应力 一应变关系。增量法是将载荷划分为许多增量，每次施加一个载荷增量。在一个载荷增 量中，假定刚度矩阵是常数；在不同的载荷增量中，刚度矩阵可以有不同的数值，并与 应力一应变关系相对应。混合法同时采用了迭代和增量法，即将载荷也划分为载荷增 量，但增量的个数较少，并对每一个增量进行迭代计算。 2 2 2 1 全量法( 迭代法) 用全量法求解非线性问题时，一次施加全部载荷，然后逐步调整位移，使基本方程 f 2 1 3 1 得到满足。主要有以下方法： 7 内蒙古科技大学硕士学位论文 2 2 2 1 1 直接迭代法( 牛顿法) 求解非线性方程组( 2 1 3 ) 的一个最简单的方法是直接迭代法。它首先给出一个近似 的解 瓯 ，再由应力一应变关系求出k ( 磊) 】= 【k 】，由式( 2 1 3 ) 就可求得第一个改进的 近似解： 慨 = 慨- 1 p ) ( 2 1 4 ) 重复这样的过程，从第”次近似解求第竹+ 1 次近似解的公式为： k 】2 【k ( 瓯) 】 r 2 1 5 1 瓯+ 。) = 置，r 矿) 计算过程如图2 1 所示： 图2 1 直接迭代法 重复上述计算，直至前后两次计算结果充分接近为止。其判断标准可取为： i j l = q a a i 占) 收敛标准为f f 口f 色l f ( x o ) + f ( ) 一) = 0 ( 2 1 6 ) 设厂( x o ) 0 ，则上式的解为= x 0 一f ( x o ) f ( x o ) ，重复上述过程，得到 f ( x ) = 0 的第+ 1 次近似解为一= x o f ( x o ) f ( ) 。以上是n e w t o n - r a p h s o n 法。 对于方程组： 砂 = 医( 占) p ) 】一 p = 0 设 瓯) 是上式的第n 次近似解。有： ) = 【k ( 瓯) 】 瓯) 一 | p ) 0 ( 2 ，1 7 ) 在 田= 瓯 处将 吼) 作泰勒展开，只保留线性项，得到： 8 一 内蒙占科技人学硕士学位论文 妒 = ) + f ( j 一 瓯) ) = 0( 2 1 8 ) 由此得到第n + 1 次的近似解： 瓯+ ， = 皖) 1 f 】“ 以。)( 2 1 9 ) 式中， 髟一为切线剐度矩阵。 2 2 2 1 2 倦正的n e w t o n 法 在n e w t o n 法中，每一次迭代都要重新建立刚度矩阵并求逆矩阵，对于大型问题来 说，这个过程在时间上是不经济的。如果只在第1 次迭代时建立刚度矩阵并求出其逆矩 阵，在以后的各次迭代中都使用这个逆矩阵进行计算，则第n 步的迭代公式为： 瓯+ ， = 吒) - 一 舛一 奶( 2 ，2 0 ) 这种方法称为修正的n e w t o n 法，如图2 2 所示： 为了提高修正n e w t o n 法的收敛速度，可以在每经过k 次迭代后重新计算一个 【墨】，即在式中将【霹】改为【彰】= 【k 。( 艿，) 】，( ，= ，2 k ， ) 。 图2 2修正的牛顿法 2 2 2 1 3 拟n e w t o n 法 n e w t o n 法也是一种修正的n e w t o n 法，其迭代公式为： 瓯+ ， = 皖) - 群】- i 纯 ， 在每次迭代时都要计算刚度矩阵并求出它的逆。在拟n e w t o n 法中建立一个【群一的递 推公式，由 f 】- l 直接计算 f “】。先令【e 吖 _ l ，则： 。) = ( 一【c 。1 吼)( 2 2 1 ) 由 e 】推算 g + 。】，有两个公式：d e f 公式和b f o s 公式。d e f 公式是一个比较成 功的算法，但数值稳定性较差：b f g s 公式是目前最成功的算法之一，具有良好的稳定 ，胜。 一9 一 内蒙古科技大学硕士学位论文 2 - 2 2 2 增量法 采用增量法分析非线性问题时，把载荷划分为许多载荷增量，每次施加一个增量。 在每步计算中，假定方程是线性的，刚度矩阵为常数，在不同的载荷增量中，刚度矩 阵可以具有不同的数值。增量法实质是用一系列线性问题去近似非线性问题，是用分段 线性的折线去代替非线性曲线。 增量法的具体算法为：把载荷分为m 个增量步：二，2 ，相应的位移也分成若干 步；a l ，a 2 ，a 3 。将( 2 1 3 ) 写成如下形式： 妒( 口) = 尸( d ) ) + 矾= 0 ( 2 2 2 ) 其中五是用以表示荷载变化的参数。进一步有： 粤) = 一k r l ( 口) ) 砺) ( 2 2 3 ) l 用e u l e r 法解上述的微分方程组，有： + 1 ) 一 a 。) = - k r l a 。) 瓴) a x ) = 一( 巧) 1 1 ) a u ) ( 2 2 4 ) 其中： 2 五一) 一 乃) 眠) = f m + 。 - ( l ) 2 。2 2 3 方法比较 全量法的计算量较小，对计算精度也能加以控制，但不能给出载荷一位移曲线，适 用范围也较小；增量法的优点是适用范围广泛，可提供载荷一位移曲线，其缺点是计算 时间较长，另外其近似解与真解差异大小不易确定；混合法在某种程度上包含了全量法 和增量法的优点，避免了两者的缺点，减少了对每一载荷增量的计算，由于进行了迭 代，可以估计近似程度，但计算量仍然较大。 2 3 弹塑性计算原理口1 】 设屈服面为 f ( 打i ，r ) = 0 ( 2 2 5 ) 在本增量步丌始时，已知应力，有效塑性应变，硬化参数和变位。弹塑性本构方程 为 d = ( d l p l ) = 【d l 占 ( 2 - 2 6 ) 2 3 1 变位增量 由上一增量步的 a i r 。计算弹塑性矩阵及各单元切线刚度矩阵，形成整体切线刚 度矩阵后，计算第n + l 步的变位增量如下： 1 0 内蒙古科技大学硕七学位论文 占l + 。= 医，r 1 ( 廿k + ，一概 )( 2 2 7 ) 然后计算应变增量，并对每个高斯积分点，按下述步骤计算其应力 盯 n + 1 。 2 3 2 试探应力 试探应力为 扛。j = 扫l + i d l 豁占l 。 ( 2 - 2 8 ) 把上述试探应力带入( 2 2 5 ) 式中判断单元屈服与否。如果满足，表明此时材料是 弹性的，硬化参数保持不变。第f + 1 增量步末尾的应力即为试探应力。若材料的弹性非 线性比较强，可以把应变增量等分为m 部分，然后顺序进行m 次计算，每次都用当时 相应的弹性模量进行计算。 2 ，3 3 比例因子 如果验算的结果，屈服条件被试探应力破坏了。表明该积分点已进入塑性状态。这 时又分为两种情况：一种是增量步开始时的应力处于屈服面以内，另一种情况是增量步 开始时的应力正好处于屈服面上。 对于第一种情况，假设本步丌始时的应力在a 点，并且符合弹性条件。即 ，( 彩，r ”) = 磊 0 ( 2 2 9 ) 表鹳属服条件已被破坏。应力路径在a 点穿过屈服面上c 点到达屈服面外b 点。 这时可把荷载分为两部分：一部分是弹性的与应力路径a c 相应；另一部分是塑性的， 它控制在c 点达到屈服面以后的行为。为了决定c 点的位置，设 p 。) = p 。) + r 仃。 ( 2 3 0 ) 其中r 是一个比例因子。可采用如下近似计算法： ：一去 ( 2 3 1 ) 1 e r 、7 实际上函数f 是非线性的，所以 f ( 扫。 + 1 盯。 s ? ，r ”) = e 0 ( 2 3 2 ) 最= 一剽时沁 r 2 + 2 1 一i 委i 棚 q 。3 4 内蒙古科技人学硕十学位论文 2 3 4 塑性应力增量 应力增量已被分为两个部分： 塑性应力增量： 盯”) = 亡吼拙) 弹性应变增量和塑性应变增量。可由下式( 2 3 5 ) 计算 ( 2 3 s ) 弹性应力增量为 仃。) = f 8 【d 】船 = p 弘占) ( 2 - 3 6 ) 应力增量满足本构方程，有 盯 = 协盯。 一 盯9 ( 2 3 7 ) 由于塑性刚度矩阵随应力状态而变化，一般利用式( 2 3 8 ) 进行近似计算： 盯j ( 1 一r ) d 】， f ( 2 3 8 ) 如果荷载和应变增量都比较小，可直接利用式( 2 3 8 ) 计算。反之，可将( 1 一r ) s j 戈0 分为m 等份，每次都用上式计算，但每次都采用新的【d 】。值。 2 3 5 应力拉回屈服面 以上已算得第n + l 增量步末尾得应力如下： 扫l + = 扛k + 忸盯。 一讧盯， ( 2 3 9 ) 由于以上步骤中有两处近似计算，必须修正计算中误差的积累，使应力拉回屈服 面。f ( 扛上彤，芷”) = e 0 ，应力修正量为： 协卜蔫 ( 2 4 0 ) 最后得到应力增量如下： p k + = p l + 仃8 ) 一 盯9 ) + 衍 ( 2 ，4 1 ) 2 3 6 计算步骤 ( 1 ) 由本增量步开始时的扫k 计算 d 】。及各单元切线模量，形成整体切线模量矩阵， 再由荷载增量 印l + ，计算变位增 a 6 一1 - i q a p n + l ； ( 2 ) 出 6 n + l 计算各高斯积分点的 a n i ； ( 3 ) 计算弹性应力增量 o 。) = 【d 。 ae ) ”1 及试探应力 e ) = n ) + e ) ； ( 4 ) 计算相应于 o 。 的不变量i i , 3 2 , 1 3 及e 。 ( 5 ) 利用当时的瞬时屈服函数，计算f ( o 。 ) f 1 ； 1 2 内蒙古科技大学硕十学位论文 ( 6 ) 若f - 为负值，表明处于弹性阶段，省去以下步骤，令 o ) n + - = o 。】： ( 7 ) 若f l o 且f o o ，转至第( 9 ) 步： ( 8 ) 若f o = o ，令r = o 及1 - r = l ，并转至第( 1 0 ) 步； ( 9 ) 由式( 2 3 4 ) 计算r ； ( 1 0 ) 计算 o = ( o 。) + r o 。 决定m 并计算增量 de l = ( 1 - o m e 什+ l 及，重复( 1 1 ) - ( 1 7 ) 步m 次： ( 11 ) 计算 。) i 的应力不变量i i ，2 2 ，j s ； ( 1 2 ) 利用e ? ，从单轴应力一应变曲线计算硬化参数h ； ( 1 3 ) 计算 a f c q a j 及【d l ， a f 8 0 , ； ( 1 4 ) 计算五及p 】。，若五 o 令兄= 0 ： ( 1 5 ) 计算 d 盯) ，= d c r 8 ) ，一【d 】， 如) ，更新应力p h = 盯) ，+ 拈) ，； ( 1 6 ) 计算单位体积内的塑性功和塑性应变增量和有效塑性应变增量； ( 17 ) 由式( 2 4 0 ) 计算【面。 ，修j 下应力后转至( 1 1 ) 步，重复( 1 1 ) - - ( 1 7 ) 步： 0 8 ) 计算本步未失衡力妇。 ，并合并到下一步的荷载增量 只+ 1 ) 中去； 2 4 弹塑性大变形的理论口2 】 2 4 1 物体的构形及其描述 在小变形固体力学中，因为变形很小而忽略了物体受力以后空间位置的改变。具体 在求解时，所采用的几何方程是线性的，本构方程和能量原理是相对于变形前的状态作近 似的描述。但在变形很大的有限变形中，必须考虑物体空间位置的改变。求解时，采用 非线性几何方程。 为了描述物体变形前后两种不同的状态，引入了参考构形和变形构形。在某一瞬 时，物体在空间所占据的区域成为该物体的构形。如果已知某物体在t = o , t l ，t 2 ，b ，( m = 1 ，2 ，3 ) 各时刻的力学特性，那么在研究该物体在k 时刻的大变形问题时，就可 以把t m 时刻以前任意时刻的构形作为参考构形。一般，以初始构形为参考构形。 令t = o 时，物体的初始构形为v o ，并参考于一个固定的直角坐标系 x i ) ，物体上的 任一质点p 的位置坐标为x i ( i = l ，2 ，3 ) 。设此后在某一瞬时t ，物体发生了变形，此时物体 的构形称为变形构形。不妨用另一直角坐标系 x i 来描述此构形。初始构形中的p 点， 变形后被移动到q 点，其位置坐标为x i 、( i = 1 ，2 ，3 ) 。如果令坐标系( x 与 x i 、 重合，则变 形后的构形如图2 3 所示。显然，对于同一质点而言，变形前后坐标的关系为 x i l - - - - x i ( x 1 , x 2 , x 3 ，t )i - 1 ，2 ，3( 2 4 2 ) 一1 3 内蒙占科技人学硕士学位论文 假设物体的变形和运动都是连续的，则v o 中每一点x 。仅为v 。e e - - + 质点x i 对应， 反之亦然。于是，可以认为函数是单值、连续、可微的，且雅戈比行列式不等于零，即 j = 肄。一2 s 亿a ， 也有x ，= x i ( x l ，工2 ，工3 ，r ) i = l ，2 ，3 ，：i 熹l 。j = 1 , 2 , 3 j 靠，。 f 2 4 4 ) f 2 4 5 ) 图2 3 物体的构形描述 在连续介质力学中，对物体内质点的运动，可以用两种方式来描述和研究，一种是 以x i ( i = l ，2 ，3 ) 和t 作为独立的变量来描述物体的运动，称为物质描述或拉格朗日描述，其 中x 。和t 称为拉格朗同变量。这种描述是随着运动的质点来研究质点的运动状态。另一 种描述是以x 0 = 1 , 2 ，3 ) , 1 1t 作为独立的变量来描述物体的运动，称为空间描述或欧拉描 述，其中x ，和t 称为欧拉变量。空间描述是指在空间某一定点研究各质点经过该点的状 态。通常，采用物质描述来求解固体力学方面的问题，而空间描述则更适合于求解流体 力学方面的问题。 2 4 2 弹塑性大变形本构方程 由于土体的弹塑性本构关系具有非线性性质，且与其应变过程和加载路径有关。故 本构关系是一种瞬态关系。目前，采用的本构关系可分为两类；一种是塑性力学形变理 论，又称全量理论：另一种是塑性流动理论，又称增量理论。在弹塑性有限元法中，常 用增量理论。在描述本构关系时，是按变形形态构形用欧拉参数描述的，用柯西应力张 量和阿尔曼斯应变张量表达应力与应变之间的关系。针对各向同性材料大变形弹塑性本 一1 4 内蒙占科技大学硕士学位论文 构方程，为了反映与加载历史的相关性，大变形弹塑性问题可以采用速率型的本构方 程。但由于柯西应力张量的物质导数f 。不是客观张量，所以若用变形速率张量和柯西 应力张量的物质导数r 。直接建立本构方程，将不满足本构方程的客观性条件。为此， 在定义本构方程时，应该选用满足客观性条件的久曼应力张量。 久曼应力变化率张量与阿尔曼斯应变速率张量关系的张量形式为 r f 。= ， ( 2 4 6 ) 式中d i j k 厂弹塑性张量 “。久曼应力变化率张量 p 。阿尔曼斯应变速率张量 根据塑性流动理论，对于各向同性硬化材料，按普朗特一劳斯方程来确定d 。 速率型的普朗特一劳斯方程为 一熹 ( 碱+ 1 _ - 7 ”五a , , s k 小k l 嘣叫嘶炯2 ， 式中 日：粤 ( 2 4 9 ) d 占 j 。- 克罗内克尔符号 口载荷性质判断因子。塑性加载过程或中性变载过程口= 1 ；弹性加载过程口= o 卸载过程口= o 。 d o - 。 唧2 言 式( 2 4 7 ) 是在小变形理论基础上得到的，o - g 为应力分量；o - # 为应力偏量；白为 一一p 应变分量；o - 为等效应力；f 为等效塑性应变，e 为弹性模量；v 为泊松比。 1 5 一 南 内蒙占科技大学硕士学位论文 把式( 2 4 7 ) 应用于有限变形有限元分析时，把口。改为r 。”，；，改为；一欧拉应 力o - 改为柯西应力f ，把式( 2 4 8 ) 代入( 2 4 7 ) 后得 j e f i l = 一 。 1 + vw 高咿“叫甄篙面 2 l 3 e 对比式( 2 4 5 ) 和式( 2 5 0 ) ，可得 n 一旦 。u b 一1 + v碱+ 高咿一习硒i 7 9 i k l 3 l 3 e f 2 5 0 ) ( 2 5 1 ) 这样，对于弹塑性变形，从小位移小应变的本构方程式( 2 4 7 ) 出发，就得到了更 普遍意义上的本构方程式( 2 5 0 ) ，它既是用于小位移小应变的情况，也适用于大位移 大变形的情况。 2 4 3 弹塑性大变形有限元列式 如前所述，物体变形虽然有欧拉和拉格朗日两种描述方法。由于欧拉描述中，泛函 积分域是当前构形，它本应是求解的结果，所以在求解过程中必须采用迭代求解。通常 认为采用欧拉描述的有限元法更适合于流体力学问题。而固体力学问题的求解，目前均 采用拉格朗日描述。 2 4 3 1 虚功方程 设t 时刻构形v t 中任意一点的位置坐标为x i ，相应点在变形构形为v t + a t 中的位置坐 标为x ，。a u i 是该点从v 。到v ，t 的位移增量，有 x 。= x ，+ a u ， ( 2 5 2 ) 在t + a t 时刻，按拉格朗r 描述的虚功方程为 p 沁料肌占舢r 件叫：她) 7 斜拟 ( 2 s s ) 式中，v 为参考构形v t 的体积； 五) 和f 分别为t + t 时刻物体单位体积力矩阵 和物体的一部分表面a ，上的表面力向量矩阵；张量 s 为t + a t 时刻以v t 为参考构形 的克希荷夫应力张量；j r 是用拉格朗日参数描述的物体内质点的虚位移增量： 1 6 内蒙占科技大学硕十学位论文 j l a e j 是虚格林戍变增量。以初始构形为参考构形的拉格朗f i 有限元法，称为全拉格朗 f i _ l ( t o t a l l a g r a n g e ) 法，简称t l 法。若以前一个相邻构形为参考构形，则称为修正的 拉格朗同( u p d a t e dl a g r a n g e ) 法，简称u l 法。 2 4 3 2u l 法的有限元列式 设由逐步求解己求得物体在t 时刻的构形，t 时刻构形内任一点的坐标为x j ，相应点 存什t 时刻坐标为五。构形内任意一点的位置坐标及位移场量插值得到 = n k a u ，。 式中，k 是哑标，r n 为单元节点数，n k 为形状函数。 u l 、法是把t 时刻的构形看作具有出位移u i 应变e u 和初应力s n 的参考构形，所以 若用x i 表示t 时亥0 参考构形中任意点的拉格朗日坐标，求解t + t 时亥构形，则不考 虑t 时刻构形的变形，即 e i j = 0( t 时刻) 构形v 。内任意一点的虚应变增量和虚位移增量分别为 艿 盥 = 陋p “) 。，巧 “) = 【p “) 。( 2 5 4 ) 式中， b 】几何矩阵， n 】形状矩阵。 将式( 2 5 4 ) 代入虚功方程式( 2 5 3 ) ，经推导，得平衡方程的刚度矩阵 囟n 。】。陋k y z 虹三 2 西 五。，矿 “ 。+ k 。r 疗= 五) 8 一豫j c z 彤， 式中，k j 及p 。 分别为六维及九维线性几何矩阵，院j 为弹塑性矩阵， 五 为 九维应力矩阵，k 】为初应力刚度矩阵， r 为节点载荷，仁。 为初应力节点力。令 k f = j b l r 【d 】，陋d v ( 小变形弹塑性刚度矩阵) k n f _ 【足。卜z j 三 。 而 五 d 矿 c 初应力刚度矩阵， 一1 7 - 内蒙古科技大学硕士学位论文 k 】。= k 。】+ k v c 沁卜似e 一。 所以，平衡方程为 医】。 “ 。= 猷 。