一般有限元原理

一般有限元原理

一、基本理论

有限元单元法是数值计算方法中发展较早、应用最广的一种方法。利用有限元法，可以解决经典的传统的方法难以解决或无法求解的许多实际问题。其优点是部分地考虑边坡岩土体的非均质、不连续的介质特征，考虑岩土体的应力应变特征，可以避免将坡体视为刚体，过于简化边界条件的缺点，能够接近实际从应力应变的角度分析边坡的变形破坏机制。对了解边坡的应力分布及应变位移变化很有利。

有限单元法实质是变分法的一种特殊的有效形式，其基本思想是：把连续体离散化为一系列的连接单元，每个单元内可以任意指定各种不同的力学形态，从而可以在一定程度上更好地模拟地质体的实际情况，特殊的节理元，可以有效地模拟岩土体中的结构面。

在大多数情况下岩土体材料应采用非线形模型，其中包括岩体弹塑性、蠕变、不抗拉特性以及结构面性质的影响。下面简要叙述有限元法的求解过程和原理。

有限单元法的基本原理

1.有限单元法的实施步骤

有限元的重要步骤归纳起来，主要有以下几步：

（1）建立离散化的计算模型，包括以一定型式的单元进行离散化，按照求解问题的具体条件确定荷载及边界条件；

（2）建立单元的刚度矩阵；

（3）由单元刚度矩阵组集总体刚度矩阵，并建立系统的整体方程组；

（4）引入边界条件，解方程组，求得节点位移；

（5）求各单元的应变、应力及主应力。

2位移模式与单元类型

在一般的有限单元法问题中，我们常以位移作为未知数，称为位移法。为保证解的收敛性，要求位移模式必须满足以下三条：

（1）位移模式必须能包含单元的刚体位移。即当节点位移是由某个刚体位移所引起时，弹性体内不会有应变。

（2）位移模式必须能包含单元的常应变，即与位置坐标无关的那部分应变。

（3）位移模式在单元内要连续，并使相邻单元间的位移必须协调。 同时，还要求所选的位移模式与局部坐标系的方位无关，即具有几何各向同性。对于线形多项式，各向同性的要求通常就等价于必须包含常应变状态，对于高次模式，就是位移模式不应随局部坐标的更换而改变。对于常应变三角形单元，其位移模式十分简单。这里以常用的四边形等参数单元为例。

等参数单元的概念已普遍地应用于有限单元公式中。“等参数”这个术语意味着单元的未知位移和几何形状有着共同的参数表述，其基本思想使用同样的内插函数N 来表达单元的位移和几何形状。

如果把单元中点的位移表示为：

}]{[}{q N u =

在等参概念中，用同样的函数N 来表示单元中某点的坐标

}]{[}{n x N x =

式中，],[}{y x x T =，],,[}{21 x x x T n =由节点坐标所构成。 先以四节点的四边形等参元为例，为了便于积分运算，采用了自然坐标系

η-

3对于二维应力－应变问题来说，单元上的每个节点有两个自由度，x 方向的位移u 和y 方向位移v ，在矩阵记号中

],[}{v u u T =

用节点位移表示就可以写作

}{u =?

?

????4433221100000000N N N N N N N N }{q 此处i N 有上面方程给出，而，

[]44

3

3

2

2

1

1

}{v u v u v u v u q T =

单元的几何形态可用同样的[N]表示为

T y x

x }{}{==??????][]0[]0[][N N }{n x 此处[]43

2

1

4321

}{y y y y x x x x x T n =，在平面应变条间下的应变位移

关系为{}}]{[////}{}{q B y u x u y u x u T

xy y x =??+??????==τεεε式中[]B 可有

N 的某种适当求导而得

=][B ???

????

?????

?

??????

???????????????????????????????????????x N y N y

N x N x N y N y N x N x N y N y N x N x N y N y N x N 434433332222111100000000 以及][i B =???

????

??

?

??????????????????x Ni y Ni y Ni x Ni 00

下面使用变分法来推导单元的刚度方程。 位移法得变分泛函系统的势能p ∏可以表示为

p ∏＝11

)()(),(ds u T u T du u Y u X v u du S y x V

V

????????+-+-

式中),(v u du ――单位体积的应变能；

X 、Y ――规定的体积力分量，即物体的重量；

x T 、y T ――规定的表面牵引力；

V ――单元体积；

1S ---有规定牵引力作用的曲面。

根据材料是线性的假设，利用弹性理论的结果，将du 表示为

{}{}dV du T

σε2

1=

此处，{}T xy x x τσσσ=}{为应力分量向量。

利用广义虎克定律，应力--应变关系是

}]{[}{εσD =

式中，}{D 为应力－应变矩阵，对于各向同性线弹性材料来说，它有一对参数构成，这对参数为杨氏模量E 和泊松比μ，或为体积模量K 和剪切模量G 或为拉梅常数λ和υ。

考虑各向同性线性弹性性状，方程

p ∏＝11

)()(),(ds u T u T du u Y u X v u du S y x V

V

????????+-+-变为

p ∏＝

11}{}{}){}{2}}{{}({21ds T u dV x u D S T

V

T T ?????--εε 式中，T Z Y

X

x }{}{=，{

}{}z y

x T T T T =分别是体力和表面牵引力向量。

将{}{

}T

xy y x τεεε=={}T

y u x u y u x u ??+??????////=[]{}q B 代入方程上式得： p ∏＝

?????--1

1}{][}{})]{[}{2}]{][[][}({21S T

T V T T T ds T N q dV x N q q B D B q 对节点位移取p ∏得一级变分，以及引用最小势能原理，从而得到

0=∏p δ

[]{}{}Q q K =

此式对于四边形等参单元来说有

[]h K =dsdt J B D B T ??+-+-111

1

]][[][

{}h q =dsdt J x N T ??+-+-111

1

][][+?1

1}{][s T ds x N

式中，J ――雅克比矩阵[]J 的行列式，导出整体坐标与局部坐标之间的关系。

[]????????????????=?

???????????????-ηξi i T

i i N N J y N x N 1，()4,3,2,1=I

[]1-J ＝4

[]J =????

?

?

??????????????ηη

ξξy x

y x

＝i n i n

j j i j i i y N N N N x ∑∑==????-????11)(

ξηηξ h 为单元的等厚度常数。对于平面应变情形来说，取1-h ，因此，对于平面问题，图1给出的四边形等参单元的刚度矩阵就可以写成

[]

e

K 4

＝ηξd d J B D B T ??+-+-111

1

]][[][

上式的积分在程序中用高斯积分实现。

相同的推导，可以得到常应变三角形单元的刚度矩阵

[]e K 3

＝ηξd d J B D B T ??+-+-111

1]][[][

对于[]D 及[]B 为常量的常应变三角形单元有

[][][][]t B D B K T e *?=3

式中 ?――三角形单元面积 t ――单元厚度 4.非线性有限元法的求解

岩土工程问题大都为非线性问题，应力应变关系呈非线性状态，非线性算法是有限元解题步骤中非常重要的一步。求解非线性问题的方法可分为三类：增量法、迭代法和混合法。这里主要介绍迭代法，迭代法在每次迭代过程中都施加全部荷载，但逐步修改位移和应变，使之满足非线性的应力应变关系，即荷载也划分为若干增量，而对每一个荷载增量进行迭代计算。

迭代法又称牛顿法，这种方法的特点是全部荷载一次施加，逐步调整位移进行迭代，最终使方程得到满足。这里主要介绍牛顿法和修正牛顿法。

(1)牛顿法

由非线性方程[]{}{}P u u K =?)(出发，从初始刚度[]0K 求得位移}{1u ?，

{}[][]P K u 101-=? (4-12)

由}{1u ?求得}{u ?，由}{1u 从u P -曲线上求得割线刚度][1K ，再由

[]{}P u u K =?)(1求得}{u ?的第二次}{2u ?，如此重复计算，直到}{1u 与}{1-?i u 充

分接近，使得{}{}[]εε≤=?-??-u i u u 11，即给定精度为止，这一过程可由表示，这种方法收敛快，但每一次迭代都要形成新的刚度矩阵，计算量较大。

(2)修正牛顿法

对上述方法的一种修正是每一步迭代步骤均采用初始刚度[]0K ，迭代方程可写为[]{}{}{}10--=?i i P P u K (4-13)

{}{}{}i i i u u u ?=-1