有限元法在岩土工程中的应用

岩土数值分析

－有限元法在岩土工程中的应用河海大学岩土工程研究所朱俊高§1.概述

朱伯芳《有限单元法原理与应用》：有限单元法是在六七十年代发展起来的强有力的数值分析方法，它使许多复杂的工程问题迎刃而解，而且由于前后处理技术的发展，计算效率非常高，实际应用越来越广泛。有限元（FiniteElementMethod FiniteElementModelling---FEM

FiniteElementAnalysisFEA)

一、有限元简介

§1.概述

有限元方法可以用来求解多种问题：这里应力变形问题有限元方法把系统（结构）看作是由无限多个单元组成的连续体，在解这一连续体时将连续体离散化，然后将物理方程、平衡方程、几何方程结合起来，变换为求解线性方程组问题。单元与单元之间只通过结点连接有限元分析可概括为六个步：有限元基本思想

§1.概述

物理方程有限元基本思想

平衡方程几何方程最终建立外荷载与位移之间的关系§1.概述

有限元分析可概括为六步：1.结构的离散化：将分析结构系统划分成有限个单元体，并在单元体的指定点设置节点，把相邻的单元体在节点处连接起来组成单元的集合体，以代替原来的结构。一般情况，单元划分越细则表述变形情况越精确，即越接近实际变形，但计算量越大。2.选择位移模式：位移模式，就可导出用节点位移表示单元内任意点位移有限元基本思想

§1.概述

3.建立单元结点力和结点位移之间的关系有限元基本思想

§1.概述

4.计算结点荷载(包括集中力、面力、体积力)5.集合所有单元的刚度方程，建立整个结构的平衡方程有限元基本思想

6.方程组求解（1）用离散结构代替原结构：离散结构由许多单元体组成，各单元只在结点处有力的联系。（2）荷载简化到结点上，{R}（3）由s～ε→[D]→[K]（4）解线性方程组[K]{δ}={R}（5）由位移{δ}→应变{ε}→应力{s}（6）对非线性问题，需重复（3）-（6）步。有限元基本步骤

非线性：材料非线性与几何非线性§1.概述

（1）平面（应力、应变）有限元：三结点、四结点、6结点、8结点（2）轴对称有限元(本质上是二维问题)（3）三维有限元有限元类型

（1）总应力法（2）有效应力法二、有限单元法（FEM）的优点及应用情况

可用于非均质问题，多层土、多种材料、多区域；可用于非线性材料，各向异性材料；可适应复杂边界条件；可用于各种类型的问题：应力变形、渗流、固结、流变、湿化变形、动力、温度问题等。

优点：缺点：单元形态对计算结果影响较大；计算比较复杂、麻烦；计算模型、参数对结果影响大；非连续性问题困难；岩土工程正是存在上述问题，因此，有限元得到了广泛地应用。在水利工程中常见的是土石坝，随着筑坝技术的不断提高，土石坝已经向300米级高坝发展。高土石坝一般不是均质坝，如有心墙、砼防渗墙、砼面板，墙、板等的受力与土体变形是联系在一起的，相互作用，相互影响，基坑的支护问题，边坡稳定问题，地下洞室，地基基础，隧道盾构施工等。材料本身不是弹性材料，边界条件十分复杂，解析解是不可能的。只有依靠数值解。可用于各种类型的问题：应力变形、渗流、固结、流变、湿化变形、动力、温度问题等对岩土体进行应力变形（应力应变）有限元分析的目的：了解整个结构的位移场和应力场，从而，位移的分布、大小，应力集中的部位，塑性区的大小、范围，为设计及采取相应措施提供依据。塑性区：强度问题，有限元精度较差；尤其剪切带等

又如边坡稳定问题例一：砼面板坝

土与结构的相互作用问题？面板是受拉？受压？（顺坡向）面板哪些部位受压？面板哪些部位受拉？上游立视例二：土质心墙堆石坝

糯扎渡、双江口拱效应例三：三峡二期围堰一道？塑性混凝土？例四：基坑工程开挖问题例五：基础工程例六：地基沉降路基变形问题例七：边坡工程1、有限元极限平衡法：利用有限元计算得到的应力场，认为沿某假定的滑动面达到极限平衡状态时稳定安全系数为1，并定义：Fs＝tf/t。2、强度折减法：c、tanf作折减后分别进行有限元计算，依据计算得到的位移场等进行稳定性判定3、容重增加法：例八：地下洞室洞室开挖的变形、稳定盾构施工、地下空间但是，有限元对动态过程的模拟能力很有限例八：地下洞室盾构施工例八：地下洞室盾构施工第二章土体应力计算土工有限元需特殊处理的问题：（1）本构关系选用，模型参数的确定；（2）计算方法：非线性分析方法、特殊问题（单元破坏、湿化、分期施工、接触问题）的处理。非常重要：要注意岩土体的特殊性有限元计算分析：总应力法和有效应力法；都有线弹性和非线性分析之分。对透水性强的地基或土工建筑物，可用总应力法进行计算；但由于该法较简单，也常用于分析饱和粘土的应力变形。一般情况下，饱和粘土地基或土工建筑物，较严密的方法为有效应力法。§2.非线性分析方法

（1）用离散结构代替原结构离散结构由许多单元体组成，各单元只在结点处有力的联系。（2）荷载简化到结点上，{R}（3）由s

～ε→[D]→[K]（4）解线性方程组[K]{δ}={R}（5）由位移{δ}→应变{ε}→应力{σ}（6）对非线性问题，需重复（3）-（6）步。一、有限元基本步骤

土体在未达到其极限状态之前的使用荷载作用下，应力水平较低，可以近似当作线性弹性体看待。对有些土，即使应力水平较低，应力应变关系也具有明显非线性，这时，应进行非线性分析。对于非线性材料：试验确定→σ～ε非线性的应力应变关系，即本构关系应力应变关系非线性表现为有限元中荷载与变形的非线性，R～δ非线性；非线性问题包括物理非线性（材料非线性）和几何非线性（大应变）。物理非线性是指土的本构关系是非线性的，而应变与位移的关系是线性的；土体在荷载作用下的位移与其几何尺度相比很小，因而在求出位移场后，可以用某单元原来的尺寸计算应力场。土力学中大多数问题属于物理非线性范畴。几何非线性表示单元体的几何性状的有限变化，将引起位移很大的变化。应变与位移的关系不再是线性的。除非象泥炭或吹填土等的应变量达到30％以上，否则一般不需要作几何非线性分析。在有限元计算中，实现这种非线性的方法是迭代法、增量法或增量迭代法。求解非线性问题的关键在于确定能反映土体变形特性的本构关系，并合理准确地确定模型参数。迭代是指一系列逼近正确解的递推计算。每次迭代中，土体受全部荷载的作用，因此，又称为全量迭代法、直接迭代法。二、迭代法每一步迭代相当于进行一次线性分析。首先，根据弹模、泊松比（割线方法确定）形成[D]→[K]解出位移、应变、应力，再确定新的弹模、泊松比……。由应力求割线模量、泊松比：如三轴试验的应力应变关系直接求。不需要本构模型L次迭代

（1）{ε}L-1

非线性{σ}L-1（2）EsL、VsL→[D]L（3）[K]L（4）[K]L{δ}L={R}（5）{δ}L→{ε}L→{σ}L

（一）、割线迭代法

余量迭代法是先将总荷载施加于结构作一次有限元计算，解得的应变在非线性关系上所对应的应力一般地与外荷载是不平衡的。则从总荷载中扣除计算应力所平衡了的那部分荷载，仅将剩余荷载施加于结构作迭代计算。1、切线迭代；2、常劲度迭代（二）、余量迭代法1、切线迭代；首先，根据假定的初始应力求Et、Vt或者直接由弹塑性模型求[D][K]位移、应变、应力L次迭代（1）{б}L-1→Et、Vt→[D]L

（2）[D]L→[K]L（3）{б}L-1→{F}L-1→{R}L-1（4）{DR}L={R}L-{R}L-1（5）[K]{Dδ}L={DR}L

（6）{Dδ}L→{Dε}L→εL=εL-1+DεL

（7）εL

由非线性关系

бL

（二）、余量迭代法2、常劲度迭代直接用初始模量、或直接假定的弹性常数L次迭代（1）{б}L-1→Et、Vt→[D]L

（2）[D]L→[K]L（3）{б}L-1→{F}L-1→{R}L-1（4）{DR}L={R}L-{R}L-1（5）[K]{Dδ}L={DR}L

（6）{Dδ}L→{Dε}L→εL=εL-1+DεL

（7）εL

由非线性关系

бL

（二）、余量迭代法将非线性的σ～ε关系用初应力{σ0}的变化来实现。

（三）、初应力迭代法

用初应力法迭代求解非线性方程组实际上就是设想有一个初应力{σ0}，使其非线性的应力应变关系与线性的应力应变关系＋{σ0}等效（如图）。其中，为弹性应力。事实上，前一种关系式是真实的，而后一种线性关系式是假想的，只有将{σ0}调整到适当值时，两者才能等效。调整{σ0}使得非线性方程组得到求解的方法即初应力迭代法。（三）、初应力迭代法

将非线性的σ～ε关系用初应力{σ0}的变化来实现。

（三）、初应力迭代法

步骤：（L次）(1)求初应力(2)求，(3)求(4)由可求{ε}L，从而求出和（由非线性的σ～ε求）(5)由求，；(6)比较与，直至两者接近，否则，重复（1）～（6）

弹性应力非线性性应力将非线性的σ～ε关系用初应力{σ0}的变化来实现。

（三）、初应力迭代法弹性应力非线性性应力弹性应力2、初应变法

用初始应力{σ}0确定初始[D]1→[K]1，

[K]{δ}={R}解得δ1{δ}1→{ε}1→{σ*}1＝[D]{ε}（弹性应力）；再根据{ε}1由非线性应力应变关系可求得{σ}1。

2、初应变法

步骤：（L次）(1)求初应力(2)求，(3)求(4)由可求{ε}L，从而求出和（由非线性的σ～ε求）(5)由求，；(6)比较与，直至两者接近，否则，重复（1）～（6）。

前面介绍的方法都是荷载施加到结构，引起应力和变形。实际上，对土工结构，除了荷载因素引起结构变形外，温度、蠕变、湿化（湿陷）等因素也可引起结构的变形。这时，一般采用“初应变”法迭代求解结构变形与应力场。而把这种非荷载因素引起的结构应变称为“初应变”。求初应变：如蠕变，要使结构产生这么大应变，应施加的荷载：

实际外荷载流变问题有限元

1、试验方法

2、理论方法（如邓肯－张模型）

（五）真实应力应变关系的确定迭代法困难：要有全量的σ～ε关系（有唯一性）才能进行，不能反映应力路径影响，试验是在特定的应力路径下进行，σ～ε曲线不能作为迭代收敛的依据。再说试验曲线是轴对称试样σ～ε关系，与平面、三维情况不同。增量法是将全荷载分为若干级微小增量，逐级用有限元法进行计算。对于每一级增量，在计算示（一）基本增量法荷载分成若干级三增量法（一）基本增量法步骤：（1）{σ}L-1→[D]L（2）[D]L→[K]L（3）[K]L{Δδ}L={R}L（4）{Δδ}→{Δε}→{Δσ}（5）求{δ}=Σ{Δδ}，{ε}=Σ{Δε}，{σ}=Σ{Δσ}

三增量法（二）中点增量法对每级荷载先用基本增量作一次试算，得出应力后，求该级平均应力，用平均应力对应的E、v作一次修正计算，作为结果。L级增量的步骤：(1)～(5)同基本增量法(6){σ}L=({σ}L-1+{σ}L)/2(7){σ}L→[D]L→[K]L(8)[K]L{Δδ}L={ΔR}L(9)由Δδ→Δε→Δσ(10){δ}=Σ{Δδ}，{ε}=Σ{Δε}，{σ}=Σ{Δσ}

(三)增量迭代法对每一级增量,用某一种迭代法计算多次,直至收敛。可用常劲度迭代法.也可将中点增量法重复多次。四、弹塑性计算用增量法计算关键先要判断加、卸载加载——[Dep]

卸载——[D]1．判别所依据的应力状态（1）依据前一级到本级的变化依据前一级初始应力变化到本级初始应力，实际上依据前级初始与终了两个应力状态来判别（可用基本增量法）。（2）依据本级初始与终了应力状态的变化必须先作一次试算，近似确定本级终了应力状态，才可确定正式计算用加载公式还是卸载公式。（中点增量法）2．由应力状态变化确定加、卸载存在三种可能：（1）在屈服面以内用[D]（2）在屈服面以外用[Dep]（3）由以内到以外

a.对整个当前荷载增量，对该单元都用[Dep]b.将其分成两段，分别用[D]和[Dep]BF——[D]FE——[Dep]

B：{б}0E：{б}={б}0+{Δб}BF：{Δб}‘BE：{Δб}{Δб}＇=α{Δб}

可以利用三个状态的屈服函数；来确定α

F：{б}={б}0+a{Δб}其中，E点应力{б}由第一次试算得到3．中点增量法弹塑性计算对某一级荷载增量（1）由初始应力状态作一次弹性计算

（2）判别是加载还是卸载？（3）作第二次计算（可以用初应力法或中点增量法）①卸载，用[D]②全加载，③部分加载（<1）加权平均定[D]

§3接触面单元

土与砼接触面上可能会错动、拉开，成为不连续结构。

接触面单元将两不连续部分连接起来。

1．物理关系

Goodman单元2．几何关系线性位移分布

3．虚位移原理

W=σ·w·bdW=σ

wdx

发生虚位移{δ}

Goodman单元可能不在水平方向，这时，就需要作坐标变换，用整体坐标表示的单元劲度矩阵才能直接叠加到整体劲度矩阵中[Q]为转换矩阵模型参数：单剪或直剪试验确定Goodman单元优点：能较好模拟错动、张开；缺点：Kn较大，导致sn失真，且导致周围实体单元应力失真；两边嵌入3、Desai薄单元

它相当于普通四结点实体单元，有厚度。但在考虑本构关系时，将切向与法向分量分开考虑，切向G与E、V无关。

G的变化可以反映接触面抗剪能力弱，产生较大错动，厚度t=（0.01～0.1）BDesai的一个算例：E=69000kPa,v=0.3,G=1.38kPa法向的模量可以用线弹性或非线性弹性。4、殷宗泽接触面薄单元5、摩擦单元等§4土体非线性分析

一、地基边界与初始应力状态1．截断边界假定完全约束，u=v=0

不宜采用部分约束，用滚动支座，因土与岩层的夹面上是有剪应力的，水平向的移动并不自由。截断的边界宜远离受力区。2．用无界元它是一种特殊的等参单元，使一条边可以转换到无限远。一、地基边界与初始应力状态3、地基或新填土的初始应力状态及模量a、假定初始应力状态，地基一般假定K0状态，假定s0＝{rz,K0rz}；碾压填土，可考虑前期固结压力。从而求得弹性矩阵或弹塑性矩阵；b、假定假定初始弹性模量

Eo、voc、地基分级计算（尤其对有斜坡情况），计算得到的应力保留，但位移抛弃；无论是假定初始应力状态、或假定初始弹性模量，其结果不仅影响初次加荷计算，也会对以后的各级计算有影响。

对填土，表层土经过了碾压，超固结；用邓肯模型计算E时，用于求[D]矩阵的应力不宜太小，可假定一个下限如50kPa。一、地基边界与初始应力状态4、接触面单元的初始应力有与之相邻的实体单元确定。取用接触面两边单元应力的平均值若接触面竖直：若接触面水平：若接触面倾斜：二、填土（施工逐级加荷的问题）土石坝要考虑逐步填筑加荷的过程。1．网格随施工过程增加；（商业程序：生死单元）2．新填土的网格单元的自重，形成结点荷载；也可直接作为荷载3．新增单元的初始应力{б}={0}；4．用平均自重应力形成[D]或[Dep]；

自重应力假定为，，，K0=0.5。5．计算结果的应力

（1）用有限元计算所得应力

（2）用自重应力

6．位移

二、填土（施工逐级加荷的问题）A点沉降：S=mv*Dp*h

一次加荷逐级加荷矩形分布应力S=h厚度的土层全部自重应力引起的压缩；S’=Z厚度填土重所引起的h厚土层梯形分布应力的压缩，而假定h厚度内自重引起的压缩填筑到A的前已完成。

S=mv*Dp*h

一次加荷逐级加荷可以用来表示新填土算得的位移，作为修正。

Z深度处的沉降：当填到A点高度时，一次性填土，A点以下土层不引起A点沉降；分层填筑相反附加应力分布不同：

如果分层足够多，

也可以接近曲线

也可以不作处理

修正：v’=v*s’/s位移修正是以位移瞬时完成为前提固结问题不必修正！！！逐级施工，对每一级荷载增量，荷载是一次性施加；因此，对每级荷载下位移增量修正后累加；如果分层无穷多，顶面位移为0；逐级施工两种模拟方法：1、假定初始应力，形成网格，单元数增加，求[D]；2、不形成网格，只作为荷载；三、开挖

开挖，就是使开挖面上的应力解除，1．假定开挖面应力为挖去土的自重应力

第Ⅰ层，

在斜面上，由бz和бx两分量叠加，求开挖面上与自重应力等效的结点荷载。

将等大反向的等效结点荷载作用于开挖面上个结点，作有限元计算。

第Ⅱ、Ⅲ层解除的自重应力，分别为γh2、γh3。

2．除第一层外，其它各层用有限元算得的应力:界面应力内插法

第一层开挖后（荷载近似为rz），有限元可算得第二开挖面上的应力。以结点M为例，周围单元1、2、3、4的应力由有限元求得。

假定：

用四个单元形心坐标代入，分别得出四单元的计算应力，有四个方程可得a1、a2、a3、a4，代入上式，用M点坐标代入，求得M的б，即由单元应力内插得开挖面上的应力，进而求得结点力。

由于第一层开挖，第二开挖面上的应力已不是rz，用有限元算得的应力更合理。

如бx对边界点C或D，可由其下面的结点（E、M点）线性外插：图2－21（a）图2－21（b）

3．由有限元解得的位移求开挖面上的结点力

（1）由各开挖面上的初始自重应力求相应结点力{FⅠ}、{FⅡ}、{FⅢ}{RI}{RⅡ}{RⅡI}

（实际上求得所有结点的结点力，但仅取用开挖面结点力保留后用）

（2）挖第一层{RⅠ}由{FⅠ}有累加而得[K]{Δδ}=－{RⅠ}

由第II开挖面以上的各单元，计算相应由第一层开挖引起的结点力{ΔFⅡ}1=实际求得所用结点的结点力仅对开挖面上的结点求结点荷载{ΔRⅡ}1={ΔFⅡ}1

同理求得{ΔFⅢ}1则：{RⅡ}＇={RⅡ}+Δ{RⅡ}1{RⅢ}＇={RⅢ}+Δ{RⅢ}1

（3）挖第二层[K]{Δδ}２=－{RⅡ}‘{ΔFⅡI}2={ΔRⅡI}1={ΔFⅡI}2

求{RⅢ}´′={RⅢ}＇+{ΔRⅢ}２（4）挖第三层[K]{Δδ}3=－{RⅢ}´′

4．Mana法

对各级开挖的开挖区的所有单元进行积分

保留开挖面上的结点荷载（1）由各开挖面上的初始自重应力{s}0求相应结点力