土工数值分析

1、岩土工程研究所土土 工数值分析工数值分析有限元法在岩土工程中的应用河海大学河海大学 岩土工程研究所岩土工程研究所岩土工程研究所1.1.概述概述 朱伯芳有限单元法原理与应用：有限单元法是在六七十年代发展起来的强有力的数值分析方法，它使许多复杂的工程问题迎刃而解，而且由于前后处理技术的发展，计算效率非常高，实际应用越来越广泛。有限元（有限元（Finite Element Method Finite Element Method Finite Element Finite Element ModellingModelling - FEM) - FEM) 一、有限元简介一、有限元简介 岩土工程研究所1

2、.1.概述概述 有限元方法把系统（结构）看作是由无限多个单元组成的连续体，在解这一连续体时将连续体离散化，然后将物理方程、平衡方程、几何方程结合起来，变换为求解线性方程组问题。单元与单元之间只通过结点连接有限元分析可概括为六个步：有限元基本思想有限元基本思想 yggxyxRRxym(5,6)j(3,4)i(1,2)f(7,8)e1234岩土工程研究所1.1.概述概述 物理方程有限元基本思想有限元基本思想 D平衡方程几何方程gyxzxzyzyxyzxzzxyzyyzxyxx00zuxwzwywzvyvxvyuxuxzzyzyxyx岩土工程研究所1.1.概述概述 有限元分析可概括为六步：1. 结构

3、的离散化：将分析结构系统划分成有限个单元体，并在单元体的指定点设置节点，把相邻的单元体在节点处连接起来组成单元的集合体，以代替原来的结构。一般情况，单元划分越细则表述变形情况越精确，即越接近实际变形，但计算量越大。2. 选择位移模式：位移模式，就可导出用节点位移表示单元内任意点位移 有限元基本思想有限元基本思想 yggxyxRRxym(5,6)j(3,4)i(1,2)f(7,8)e1234332211332211,uyxNuyxNuyxNyxuvyxNvyxNvyxNyxv岩土工程研究所1.1.概述概述 3. 建立单元结点力和结点位移之间的关系有限元基本思想有限元基本思想 yggxyxRRxy

4、m(5,6)j(3,4)i(1,2)f(7,8)e1234 B xNyNxNyNyNxNyNxNyNyNxNxNBmmjmjmjiiiji000000 BD eeeKR dxdydzBDBKTe岩土工程研究所1.1.概述概述 4. 计算结点荷载计算结点荷载(包括集中力、面力、体积力包括集中力、面力、体积力)5. 集合所有单元的刚度方程，建立整个结构的平衡方程集合所有单元的刚度方程，建立整个结构的平衡方程有限元基本思想有限元基本思想 yggxyxRRxym(5,6)j(3,4)i(1,2)f(7,8)e1234 RK6. 方程组求解方程组求解岩土工程研究所（1）用离散结构代替原结构：离散结构由许

5、多单元体组成，各单元只在结点处有力的联系。（2）荷载简化到结点上，R（3）由DK（4）解线性方程组 K=R（5）由位移应变应力（6）对非线性问题，需重复（3）-（6）步。 有限元基本步骤有限元基本步骤 dxdydzBDBKTe非线性：材料非线性与几何非线性岩土工程研究所1.1.概述概述 （1）平面（应力、应变）有限元：三结点、四结点、6结点、8结点（2）轴对称有限元（3）三维有限元有限元类型有限元类型 yggxyxRRxym(5,6)j(3,4)i(1,2)f(7,8)e1234（1）总应力法（2）有效应力法57624318（a）（b）（c）314625岩土工程研究所二、有限单元法（二、有限单

6、元法（FEMFEM）的优点及应用情况的优点及应用情况 可用于非均质问题，多层土、多种材料、多区域；可用于非线性材料，各向异性材料； 可适应复杂边界条件；可用于各种类型的问题：应力变形、渗流、固结、流变、湿化变形、动力、温度问题等。 优点：缺点：单元形态对计算结果影响较大；计算比较复杂、麻烦；计算模型、参数对结果影响大；非连续性问题困难；岩土工程研究所岩土工程正是存在上述问题，因此，有限元得到了广泛地应用。在水利工程中常见的是土石坝，随着筑坝技术的不断提高，土石坝已经向300米级高坝发展。高土石坝一般不是均质坝，如有心墙、砼防渗墙、砼面板，墙、板等的受力与土体变形是联系在一起的，相互作用，相互影

7、响，基坑的支护问题，边坡稳定问题，地下洞室，地基基础，隧道盾构施工等。 材料本身不是弹性材料，边界条件十分复杂，解析解是不可能的。只有依靠数值解。 岩土工程研究所对土工建筑物进行有限元分析的目的：了解整个结构的位移场和应力场，从而，位移的分布、大小，应力集中的部位，塑性区的大小、范围，为设计及采取相应措施提供依据。 稳定性实际上是首先得到应力位移场，再分析处理。岩土工程研究所例一：砼面板坝例一：砼面板坝 面板是受拉？受压？（顺坡向） 面板哪些部位受压？ 面板哪些部位受拉？ 砼面板 堆石 上游立视岩土工程研究所例二：土质心墙堆石坝例二：土质心墙堆石坝 糯扎渡、双江口拱效应 粘土心墙 岩土工程研究

8、所例三：三峡二期围堰例三：三峡二期围堰一道？塑性混凝土？ 岩土工程研究所例四：基坑工程例四：基坑工程 开挖问题岩土工程研究所 例五：边坡工程例五：边坡工程1、极限平衡法；2、强度折减法岩土工程研究所第二章第二章 土体应力计算土体应力计算土工有限元需特殊处理的问题：（1）本构关系选用，模型参数的确定；（2）计算方法：非线性分析方法、特殊问题（单元破坏、湿化、分期施工、接触问题）的处理。 岩土工程研究所有限元计算分析：总应力法和有效应力法；都有线弹性和非线性分析之分。对透水性强的地基或土工建筑物，可用总应力法进行计算；但由于该法较简单，也常用于分析饱和粘土的应力变形。一般情况下，饱和粘土地基或土工

9、建筑物，较严密的方法为有效应力法。 岩土工程研究所2.2.非线性分析方法非线性分析方法 （1）用离散结构代替原结构离散结构由许多单元体组成，各单元只在结点处有力的联系。（2）荷载简化到结点上，R（3）由 DK（4）解线性方程组 K=R（5）由位移应变应力（6）对非线性问题，需重复（3）-（6）步。 一、有限元基本步骤一、有限元基本步骤 dxdydzBDBKTe岩土工程研究所土体在未达到其极限状态之前的使用荷载作用下，应力水平较低，可以近似当作线性弹性体看待。对有些土，即使应力水平较低，应力应变关系也具有明显非线性，这时，应进行非线性分析。对于非线性材料：试验确定非线性的应力应变关系，即本构关系

10、 D D 岩土工程研究所应力应变关系非线性表现为有限元中荷载与变形的非线性，R非线性； R R KK 非线性问题包括物理非线性（材料非线性）和几何非线性（大应变）。物理非线性是指土的本构关系是非线性的，而应变与位移的关系是线性的；土体在荷载作用下的位移与其几何尺度相比很小，因而在求出位移场后，可以用某单元原来的尺寸计算应力场。土力学中大多数问题属于物理非线性范畴。几何非线性表示单元体的几何性状的有限变化，将引起位移很大的变化。应变与位移的关系不再是线性的。除非象泥炭或吹填土等的应变量达到30以上，否则一般不需要作几何非线性分析。岩土工程研究所 在有限元计算中，实现这种非线性的方法是迭代法、增量

11、法或增量迭代法。 求解非线性问题的关键在与确定能反映土体变形特性的本构关系，并合理准确地确定模型参数。 岩土工程研究所迭代是指一系列逼近正确解的递推计算。每次迭代中，土体受全部荷载的作用，因此，又称为全量迭代法。 二、迭代法 岩土工程研究所每一步迭代相当于进行一次线性分析，即根据弹模、泊松比（割线方法确定）形成DK解出位移、应变、应力，再确定新的弹模、泊松比。 L次迭代 （1）L-1 非线性 L-1 （2）EsL、VsLDL （3）KL （4）KLL=R （5）LLL （一）、割线迭代法 R 岩土工程研究所余量迭代法是先将总荷载施加于结构作一次有限元计算，解得的应变在非线性关系上所对应的应力一

12、般地与外荷载是不平衡的。则从总荷载中扣除计算应力所平衡了的那部分荷载，仅将剩余荷载施加于结构作迭代计算。 1、切线迭代； 2、常劲度迭代L次迭代（1）L-1Et、VtDL （2）DLKL （3）L-1FL-1RL-1 （4）DRL=RL-RL-1（5）KDL=DRL （6）DLDLL=L-1+ DL （7） L 由非线性关系 L （二）、余量迭代法 R R * * 1 1 1 1 dxdyBFTe1岩土工程研究所2、常劲度迭代L次迭代（1）L-1Et、VtDL （2）DLKL （3）L-1FL-1RL-1 （4）DRL=RL-RL-1（5）KDL=DRL （6）DLDLL=L-1+ DL （7

13、） L 由非线性关系 L （二）、余量迭代法 dxdyBFTe1 O R O R 岩土工程研究所将非线性的关系用初应力0的变化来实现。 （三）、初应力迭代法 O -D O O RA1 A2 A3 RR0R0用初应力法迭代求解非线性方程组实际上就是设想有一个初应力0，使其非线性的应力应变关系 与线性的应力应变关系 0等效（如图）。其中， 为弹性应力。事实上，前一种关系式是真实的，而后一种线性关系式是假想的，只有将0调整到适当值时，两者才能等效。调整0使得非线性方程组得到求解的方法即初应力迭代法。)(ijijfDDE岩土工程研究所（三）、初应力迭代法 0RRK O O RA1 A2 A3 RR0R

14、0 eeeTeTeRkdxdzBdxdzBDBR00 O -D10*LLLD岩土工程研究所将非线性的关系用初应力0的变化来实现。 （三）、初应力迭代法 步骤：（L次）(1) 求初应力(2) 求 ， (3) 求(4) 由 可求L，从而求出 和 （由非线性的求）(5) 由 求 ， ；(6) 比较 与 ，直至两者接近，否则，重复（1）（6） 11*0LLLL0dxdyBFLTeL00eeLLFR00LLRRK0L10*LLLDLdxdyBFLTeL eeLLFRL LLL O -D O O RA1 A2 A3 RR0R010*LLLD岩土工程研究所2、初应变法 用初始应力0确定初始D1K1， K=R

15、解得111*1D（弹性应力）；再根据1由非线性应力应变关系可求得1。 O e A B O 岩土工程研究所2、初应变法 步骤：（L次）(1) 求初应力(2) 求 ， (3) 求(4) 由 可求L，从而求出 和 （由非线性的求）(5) 由 求 ， ；(6) 比较 与 ，直至两者接近，否则，重复（1）（6）。 11*0LLLL0dxdyBFLTeL00eeLLFR00LLRRK0LLLD*LdxdyBFLTeL eeLLFRLLLL 岩土工程研究所前面介绍的方法都是荷载施加到结构，引起应力和变形。实际上，对土工结构，除了荷载因素引起结构变形外，温度、蠕变、湿化（湿陷）等因素也可引起结构的变形。这时，

16、一般采用“初应变”法迭代求解结构变形与应力场。而把这种非荷载因素引起的结构应变称为“初应变”。求初应变：如蠕变，要使结构产生这么大应变，应施加的荷载： dxdyDBFTe00eeFR000RRK岩土工程研究所1、试验方法 2、理论方法（如邓肯张模型） （五）真实应力应变关系的确定 O -D岩土工程研究所迭代法困难： 要有全量的关系（有唯一性）才能进行，不能反映应力路径影响，试验是在特定的应力路径下进行，曲线不能作为迭代收敛的依据。再说试验曲线是轴对称试样关系，与平面、三维情况不同。岩土工程研究所增量法是将全荷载分为若干级微小增量，逐级用有限元法进行计算。对于每一级增量，在计算示（一）基本增量法

17、 荷载分成若干级 三 增量法 O l1 l N N Nl Nl1 ll l R O l1 l M M Ml Ml1 l Rl 岩土工程研究所（一）基本增量法 步骤： （1）L-1DL （2）DLKL （3）KLL= RL （4） （5）求=，=，= 三 增量法 O l1 l N N Nl Nl1 ll l R O l1 l M M Ml Ml1 l Rl 岩土工程研究所（二）中点增量法 对每级荷载先用基本增量作一次试算，得出应力后，求该级平均应力，用平均应力对应的E、v作一次修正计算，作为结果。L级增量的步骤： (1)(5)同基本增量法 (6)L=(L-1+L)/2 (7)LDLKL (8)K

18、LL=RL (9)由 (10)=，=，= l O Nl l l-1 Nl-1 Dl R O Ml Ml-1 Kl R Rl-1 岩土工程研究所(三)增量迭代法对每一级增量,用某一种迭代法计算多次,直至收敛。可用常劲度迭代法.也可将中点增量法重复多次。 岩土工程研究所四、弹塑性计算 用增量法计算 关键先要判断加、卸载 加载Dep 卸载D 1判别所依据的应力状态 （1）依据前一级到本级的变化 依据前一级初始应力变化到本级初始应力，实际上依据前级初始与终了两个应力状态来判别（可用基本增量法）。 （2）依据本级初始与终了应力状态的变化 必须先作一次试算，近似确定本级终了应力状态，才可确定正式计算用加载

19、公式还是卸载公式。（中点增量法） 岩土工程研究所 2由应力状态变化确定加、卸载存在三种可能： （1）在屈服面以内 用D （2）在屈服面以外 用Dep （3）由以内到以外 a.对整个当前荷载增量，对该单元都用Dep b.将其分成两段，分别用D和Dep BFD FEDep A B C c D E F G 岩土工程研究所 B：0 E：=0+ BF： BE： = 可以利用三个状态的屈服函数；来确定 EBF 0DD )()()()(BEBFffff已知）（）（当前的屈服函数，0)(HFffFF)()()()(0BEBfffHFF：=0+ 岩土工程研究所3中点增量法弹塑性计算 对某一级荷载增量 （1）由初

20、始应力状态作一次弹性计算 （2）判别 是加载还是卸载？ （3）作第二次计算 卸载，用D 全加载， 部分加载 （ 1挤压，甚至达4.0。 K01.0uiCE1600800岩土工程研究所 （二）群桩1.简化为二维问题 桩 板桩 2.半解析元法在y向用级数展开，以反映材料的变化和应力、变形的变化及xz向有限元数值解。3.三维有限元法单元、结点太多：群桩承台模拟：变刚度随着计算机容量、速度的提高，也不成问题。 DLnAEE 岩土工程研究所 五、土坝浸水变形五、土坝浸水变形 蓄水对坝的作用：（1）水压力作用于心墙，面力或渗流力（2）浮托力作用于上游坝壳（3）湿化变形土颗粒重新调整位置，体积缩小 岩土工程

21、研究所 1、 湿化试验（1）直接法（单线法）：干 湿（对于给定的应力状态）优点：反映了由干到湿的过程对变形的影响。缺点：不能反映应力路径的影响，试验工作量大； 由 和 转为应变分量还要作假定； 轴对称试样的湿化变形与平面应变有差异。 ),(31faD),(31fvDaDvD岩土工程研究所 （2）间接法（双线法）：分别作干样和湿样的变形试验。干： 确定邓肯模型参数湿： 确定邓肯模型参数假定湿化变形缺点：作了上假定，没有反映湿化过程。这一假定不一定合理，也许宜作修正。目前看来，大部分试验资料表明，这样得到的湿化变形偏小；优点：试验工作量小；不必进行总应变与应变分量的转换。 va31va31 干湿

22、岩土工程研究所 2、计算（1）假定结点位移限制，求应力改变量用直接试验资料（单线法）由用间接法 分成若干增量湿D wwwDDDD31 n dD wwDww wD假定应变同干态，则用湿态应力应变关系得到的应力比原小需将主应变分量转换为各应变分量，还需假定中主应变增量岩土工程研究所 DDTi将主应变分量转换为各应变分量 zxzxzxzyzyzyyxyxyxzzzyyyxxxnnmmllnnmmllnnmmllnmlnmlnmlT222222222222222222yzxyzxyxDl11zxxyyzxyDl112111xyyxzDl212121zxzyxyzyzxyzxyxzxyyzxDlx，ly

23、, lz 分别为1与x,y,z轴的方向余弦岩土工程研究所 （2）相应结点力：（3）将 与浮力、水压力对应的R相加，作为外载，解有限元方程得：（4）得结果应力 DdxdyBFTewRF wRD DDD岩土工程研究所六、单元破坏的处理六、单元破坏的处理为什么要应力修正？修正原则：莫尔库仑破坏准则 岩土工程研究所六、单元破坏的处理六、单元破坏的处理三维有限元计算中的应力修正：1、修正原则：莫尔库仑破坏准则；2、修正方法：(1)、假定修正前后大主应力不变，同时假定b (23)/ (13)或u (22 13)/ (13)不变，主应力方向不变；1*= 1 ；莫尔库仑准则 3* ； b 2* x 、y、 z

24、、xy、 yz 、 xz l1、m1、n1、 l2、m2、n2 、 l3、m3、n3x * 、y * 、 z * 、xy * 、 yz * 、 xz *(2)、假定修正前后 (13)/2不变，同时假定b (23)/ (13)或u (22 213)/ (13)不变，主应力方向不变；岩土工程研究所 p x z z x xz 六、单元破坏的处理六、单元破坏的处理1.拉裂破坏 令（1）假定 ，且主应力方向不变03p3或0311311313131112xyzx岩土工程研究所 六、单元破坏的处理六、单元破坏的处理（2）假定 不变，且主应力方向不变令 z31xzazxaa121xyxz岩土工程研究所 2.剪

25、切破坏（1）假定 不变，且主应力方向不变zzx3131cos2sin22xzxzzx岩土工程研究所 （2）假定 不变，且主应力方向不变令z31sincossin2xzzcaaxzzxa应力迁移问题： D DDdxdyBRt岩土工程研究所5 5 固结有限元固结有限元 （一）Biot固结方程1.平衡方程 fTxyyxT00 u uMu 011M D w fuMwDTT岩土工程研究所 2.连续性方程 0uktNwv wMMTTv ukwtMNwT岩土工程研究所 （二）有限元方程1.平衡方程 etFdxdyB eTTFdxdyuBdxdyB eeTeTFdxdyNMBdxdyBDB eNw eeBNw eNu uMu eeeFkk岩土工程研究所 2.连续加权余量法取近似解近似解代入分式，不为0，出现余量。令余量乘以其中权函数后，在单元内积分为0，可作出 应服从的关系式，取权函数为Ni（x,y）。miixyxNw1,yiimimiywyxNw,11miiiuyxNu1, eNwiiuw、 e0Ndxdyukwt