有限元计算原理与方法

1.有限元计算原理与方法
有限元是将一个连续体结构离散成有限个单元体，这些单元体在节点处相互铰结，把荷载简化到节点上，计算在外荷载作用下各节点的位移，进而计算各单元的应力和应变。

用离散体的解答近似代替原连续体解答，当单元划分得足够密时，它与真实解是接近的。

1.1. 有限元分析的基本理论
有限元单元法的基本过程如下：
1.1.1.连续体的离散化
首先从几何上将分析的工程结构对象离散化为一系列有限个单元组成，相邻单元之间利用单元的节点相互连接
而成为一个整体。

单元可采用各种类
型，对于三维有限元分析，可采用四
面
体单元、五西体单元和六面体
单元等。

在Plaxis 3D Foundation
程序中，土体和桩体主要采用包
含6个高斯点的15节点二次楔
形体单元，该单元由水平面为6
节点的三角形单元和竖直面为四
边形8节点组成的，其局部坐标
下的节点和应力点分布见图3.1，图3.1 15节点楔形体单元节点和应力点分布界面单元采用包含9个高斯点的
8个成对节点四边形单元。

在可能出现应力集中或应力梯度较大的地方，应适当将单元划分得密集些；
若连续体只在有限个点上被约束，则应把约束点也取为节点：若有面约束，则应
把面约束简化到节点上去，以便对单元组合体施加位移边界条件，进行约束处理；
若连续介质体受有集中力和分布荷载，除把集中力作用点取为节点外，应把分布
荷载等效地移置到有关节点上去。

最后，还应建立一个适合所有单元的总体坐标系。

由此看来，有限单元法中的结构已不是原有的物体或结构物，而是同样材料
的由众多单元以一定方式连接成的离散物体。

因此，用有限元法计算获得的结果
只是近似的，单元划分越细且又合理，计算结果精度就越高。

与位移不同，应力
和应变是在Gauss 积分点(或应力点)而不是在节点上计算的，而桩的内力则可通
过对桩截面进行积分褥到。

1.1.
2. 单元位移插值函数的选取
在有限元法中，将连续体划分成许多单元，取每个单元的若干节点的位移
作为未知量，即{}[u ,v ,w ,...]e T i i i δ=，单元体内任一点的位移为{}[,,]T
f u v w =。

引入位移函数N (x,y,z )表示场变量在单元内的分布形态和变化规律，以便用
场变量在节点上的值来描述单元内任一点的场变量。

因此在单元内建立的位移模
式为： {}[]{}e f N δ= （3-1）
其中：12315[][,,......]N IN IN IN IN =，I 为单位矩阵。

按等参元的特性，局部坐标(,,)ξηζ到整体坐标,,x y z （）的坐标转换也采用
与位移插值类似的表达式。

经过坐标变化后子单元与母单元(局部坐标下的规则
单元)之间建立一种映射关系。

不管内部单元或边界附近的单元均可选择相同的
位移函数，则为它们建立单元特性矩阵的方法是相同的。

因此，对于15节点楔
形体单元体内各点位移在整体坐标系,,x y z （）下一般取：
151151151(,,)(,,)(,,)i i i i i i i i i u N u v N v w N w ξηζξηζξηζ===⎫=⎪⎪⎪=⎬
⎪
⎪=⎪
⎭∑∑∑
32-（） 上式中的(,,)i i i u v w 为整体坐标系下节点i 处的位移值，(,,)i N ξηζ为在局部
坐标系下节点相应的形函数。

1.1.3. 单元特性分析
利用几何方程、本构方程、虚功原理或位能变分方程求解单元节点力与节
点位移关系的表达式，即单元刚度矩阵。

根据几何方程可建立单元内的应变矩阵{}{,,,,,}x y z xy yz zx εεεεγγγ=：
{}[]{}e B εδ= （3-3）
其中1215[][,......]B B B B =，
/000/000/[]//00///0/i i i i i i i i i i N x N y N z B N y N x N z N y N z N x ∂∂⎡⎤⎢⎥∂∂⎢⎥⎢⎥∂∂=⎢⎥∂∂∂∂⎢⎥⎢⎥∂∂∂∂⎢⎥∂∂∂∂⎢⎥⎣⎦
(34)- 对于小变形线性弹性问题，根据物理方程建立单元内的应力矩阵：
{}[]{}[][]{}e D D B σεδ== （3-5）
其中,[]B 为几何矩阵，[]D 为弹性矩阵，[]S 为应力矩阵，[][][]S D B =。

根据虚功原理求出单元中的节点力{}e F ：
{}[]{}e e F k δ= （3-6）
其中[]k 为单元的劲度矩阵，[][][][]T e k B D B dxdxdz =⎰⎰⎰
{}R 对于整体结构上的任一点 i ，建立平衡方程：
{}{}i i
e F R =∑ (37)-
{}i R 为i 节点上的外荷。

上式表示{}i R 与围绕i 点的各单元在i 点上的节
点力之和相平衡。

1.1.4. 总体特性分析
对每一个位移未知的节点，都可写出3-7式的方程，利用结构力的平
衡条件和边界条件把各个单元按原来的结构重新联接起来，形成分析对象
的整体有限元平衡方程组：
[]{}{}K R δ= （3-8）
其中， 为整体劲度矩阵， ； 为整个结构的节点位移矩阵，
为整个结构的节点荷载矩阵，是已知的。

由式(3-8)求出节点位移 ,
由式(3-3)、式(3-5)求出各单元的应变和应力。

1.2. 非线性有限元分析
非线性现象是在实际的结构分析中经常遇到的问题。

与线性分析相
比，非线性分析中荷载与位移之间的关系已不是直线关系，而是曲线关系。

土体的非线性分析一般来说采用非线性的分析方法，选用适当的土体本构
系，进行有限元计算。

非线性问题一般有材料非线性和几何非线性两种。

几何非线性即存在大变形，其变化的几何形状可能引起结构的非线性
[]k ij ij K k =∑
{}δ{}δ
响应，即应变与位移的关系不里线性，应变不仅包括位移对坐标的一阶导
数，还要包括高阶导数。

在进行小应变或者小变形分析时，假定位移和变
形总是足够小(这种假定取决于特定分析要求中的精度等级)可以忽略结
构变形对系统刚度的影响，即基于最初几何形状的结构刚度的一次迭代足
以计算出分析结果。

随着变形位移增长，一个有限单元的已移动的坐标可
以多种方式改变结构的刚度，进行多次迭代来获得一个有效的解，这就是
几何非线性。

除了结构大变形引起剐度变化以外。

许多与材料有关的参数同样可以
改变结构刚度。

材料的非线性即是材料的应力—应交关系是非线性的。

主
要有弹性非线性模型和弹塑性模型两大类。

弹性非线性理论是以弹性理论
为基础，在微小的荷载增量范围内，把土看作弹性材料，从一个荷载增量
变化到另一个荷载增量，土体的弹性常数发生变化，以考虑非线性；弹塑
性模型理论认为土体的变形包括弹性和塑性变形两部分，把弹性理论和塑
性理论结合起来建立的本构模型。

土体中的弹塑性本构关系都是用增量形
式表示的，因此，计算方法也宜用增量法。

某级荷载增量
作用下，各单元的应力状态不同。

有些可能处于弹性区，则刚度矩阵要用弹性矩阵
，有些可能产生塑性屈服，则须运用屈服准则、硬化规律和流动法则
建立的弹塑性刚度矩阵
来代替 。

反映到式(3-5)，其中的矩阵 不是常量其随应力或应变改变，由此推导的劲度矩阵 也随应力或
变形而变。

对于相适应流动法则
，则：
[]R ∆[]D []ep D []D []D []K g f =[D]{}{}[D][D ][D]{}[D]{}T ep T f f f f A σσσσ∂∂∂∂=-∂∂+∂∂(38)
-
式中A 为塑性硬化模量,是硬化参数函数。

因此,不管是材料非线性还是几何非线性,推出的劲度矩阵将随位移而变。

因此，不管是材料非线性还是几何非线性，推出的劲度矩阵将随位移而变。

（3-10）
这是位移的非线性方程组。

直接解这样的方程组是困难的，因此简化为一系列的线性问题的解逐步逼近非线性问题的解，非线性问题可以理解为一些线性解进行迭代的结果。

1.3. 有限单元法解比奥固结方程
对于土工问题有限元分析可以采用有效应力法、总应力法和准有效应力法三种。

有效应力法严格区分土体中的有效应力与孔隙水压力。

将土体骨架变形与孔隙水的渗透同步考虑，因而比总应力法更真实反映土体自身特性，能更合理计算土体对荷载的响应。

有效应力法有两个未知量，即土体骨架的变形和孔隙水压力。

对于非饱和土还需要增加一个孔隙气压力这个变量。

有效应力法基本上以Biot 动力固结方程为基础，其计算较为复杂，计算工作量也较大。

土体的总应力有限元法实际上与其他结构有限元分析在计算原理上没有大的区别，主要在材料的本构模型的选择上不同，其实质认为土体是一种由土颗粒和孔隙水组成之间的相互关系，将之合成一个整体，共同一个整体，共同研究其整体的应力与变形状态。

总应力法不能反映土体固结作用。

[()]{}{}
K R δδ=
在有效应力分析中，如果采用与总应力法同样的土性参数并令孔隙水压力为0，则有效应力等于总应力，相应的有效应力法转变为总应力法。

因此，总应力法是有效应力法的一个特例。

在土体材料采用不捧水指标时，总应力法计算出来的是加荷瞬间或短期应力和变形，而采用排水指标进行的总应力分析则得到的是有效应力分析的最终结果，也就是孔压消散完毕，土体固结完成时的应力和交形结果。

在土工问题分析中有时还用总应力和太沙基固结理论相结合的方法来进行有效应力分析(简称准有效应力法)，该法是先用总应力法求得应力和变形，然后根据太沙基固结理论考虑孔压的消散以及有效应力和变形随时间的变化。

这种分析法对于二维和三维渗流而已是近似的，对于只有一个方向渗水的固结问题是精确的。

在Plaxis 3D Foundation程序中，进行最终沉降分析时是材料类型为排水指标的总应力法分析，而进行固结有限元沉降分析时采用的是以Biot固结理论为基础的有效应力法．采用有效应力法可以较为全面地得到桩土的应力、变形和孔压变化的情况。

1.3.1.比奥固结理论
太沙基固结理论只在一维情况下是精确的，对二维、三维问题并不精确。

太沙基一伦杜立克理论(扩散方程)将应力应变关系视为常量(E=常数)的同时，假设三个主应力(总应力)之和不变，不满足变形协调条件。

比奥理论从较严格的固结机理出发推导了准确反映孔隙水压力消散
与土骨架变形相互关系的三维固结方程。

该理论将水流连续条件与弹性理论结合求解了土体受力后的应力、应变、孔隙水压力的生成和消散过程，
一般称为“真三维固结理论”。

两理论均假设土骨架是线弹性体，变形为小变形，土颗粒与孔隙水均不可压缩，孔隙水渗流服从达西定律。

在土工数值计算中，可使用非线性弹塑性模型代替线弹性模型与比奥固结理论耦合求解。

比奥固结理论是严格按照弹性理论，使饱和粘土在固结过程中必须满足应力平衡方程、几何方程及虎克定律，因此对于三维固结问题可导出如下三个平衡方程：
（3-11）
根据饱和土的连续性在一个元素体中，在一定的时间内单元土体积的压缩量等于流进和流出该单元体的流量变化之和，并引进达西定律，从而推导如下连续方程：
（3-12）
式(3一11)和式(3一12)联立就是比奥固结方程。

式中 、 、
— 分别为在x,y 和z 三个轴向的位移； — 孔隙水压力；
— 剪切模量；
— 泊松比； 222()012()012()12y x z x y x z y y x z z w w w G u G w v x x y z x w w w G u G w v y x y z y w w w G u G w v z x y z z γ∂⎫∂∂∂∂-∇++++=⎪-∂∂∂∂∂⎪⎪∂∂∂∂∂⎪-∇++++=⎬-∂∂∂∂∂⎪⎪∂∂∂∂∂-∇++++=-⎪-∂∂∂∂∂⎪⎭ 21()y v x z w
w w w k u t t x y z εγ∂∂∂∂=-++=-∇∂∂∂∂∂z w y w x w u
G ν
— 土的重度；
— 体应变； — 渗透系数，假设土的各向渗透性相同；
— 水的容重； — 拉普拉斯算子， 比奥固结方程中含有 、 、
和 四个未知函数·在一定的边界条件和初始条件下，可以解出任何时间及任何一点的 、 、 和。

但问题远不这么简单，就是二维问题也很难求得该未知函数的解析解。

因此，该理论虽早在1941年就提出来了，但未得到推广使用，直到近年来由于电子计算机的出现，才有人开始用有限元法，把上述理论运用于解决固结问题。

1.3.2. 比奥固结有限元方程
根据有效应力原理，总应力为有效应力和孔隙水压力之和，且孔隙水不承受剪应力。

（3-13）
（3-14）
（3-15） 其中： 为节点孔隙水压力, ， 为单元
的节点超静水压力。

由虚位移原理可推导得出单元节点力与某一时刻已产生的位移所对应的骨架应力以及尚未消散的超静水压力两部分相平衡。

γ
v εk w γ2∇2222
222x y z ∂∂∂∇=++∂∂∂x w y w z w u x w y w z w u
{}{}{}u σσ'=+{}[]{}e u N β'={}[D]{}[D][B]{}e
σεδ'=={}u 1215[][,,]N N N N '=⋯⋯{}e β{F}[k]{}[k ]{}e e e
δβ'=+
（3-16）
式中
— 就是通常单元的劲度矩阵
— 单元节点孔隙压力所对应的那部分节点力； 对于所有位移未知的节点建立整体平衡方程,得有限单元法平衡方程:
（3-17）
将每个节点周围各单元内
的“领域”连在一起形成以节
点为中心的闭合“全领域”，
对节点i 其周围各单元的边界 向外流出的流量总和为0，即 图3-2节点i 的“全领域” 对于一个单元来说是流出，对于相邻的另一单元便是流进，可对l 节点的“全领域”建立连续性方程：
（3-18） 和 分别由单元矩阵 和 中的元素叠加而成。

的元
素为节点位移所对应的“节点领域”的体积改变量：
为节点孔隙压力差所产生的水力坡降在Δt 时间内引起的从“节点领域”边界的排水量。

且 [k][k ]'[K]{}[K ]{}{R}
δβ'+
=[K]{}[K]{}0δβ∆+= [K] [K] []k []k []k []k。