有限单元法基本思想-原理-数值计算过程

有限单元法学习报告

在对力学问题分析求解过程中，方法可以概括为两种方法，一种为解析法，对具体问题具体分析，通过一定的推导用具体的表达式获得解答，由于实际工程中结构物的复杂性，此方法在处理工程问题是十分困难的；另一种是数值法，有限元法是其中一种方法，其数学逻辑严谨，物理概念清晰，又采用矩阵形式表达基本公式，便于计算机编程，因此在工程问题中获得广泛的应用。

有限元法基本原理是，将复杂的连续体划分为简单的单元体；将无限自由度问题化为有限自由度问题，因为单元体个数是有限的；将偏微分方程求解问题化为有限个代数方程组的求解问题。通常以位移为基本未知量，通过虚功原理和最小势能原理来求解。

基本思想是先化整为零，即离散化整体结构，把整体结构看作是由若干个通过结点相连的单元体组成的整体；再积零为整，通过结点的平衡来建立代数方程组，最后计算出结果。

我将采用最简单的三结点三角形为基本单元体，解决弹性力学中的平面问题为例，解释有限单元法的基本原理、演示数值计算过程和一般性应用结论。

一、离散化

解决平面问题时，主要单元类型包括三角形单元（三结点、六结点）和四边形单元（四结点矩形、四结点四边形、八结点四边形）等。选用不同的单元会有不同的精度，划分的单元数越多，精度越高，但计算量也会越大。因此在边界曲折，应力集中处单元的尺寸要小些，

但最大与最小单元的尺寸倍数不宜过大。在集中力作用点及分布力突变的点宜选为结点，不同厚度，不同材料不能划分在同一单元中。三角形单元以内角接近60°为最好。充分利用对称性与反对称性。

二、单元分析

将一个单元上的所有未知量用结点位移表示，并将分布在单元上的外力等效到结点上。

1、位移函数选取：根据有限元法的基本思路，将连续体离散为有限的单元集合后，此时单元体满足连续性、均匀性、各向同性、完全线弹性假设。单元与单元之间通过结点连接并传递力，位移法

（应用最广）以结点位移8 i= （U i V i）T为基本未知量，以离散位移场代替连续位移场。单元体内的位移变化可以用位移函数（位移模式）来表示，因为有限元分析所得结果是近似结果，为了保证计算精度和收敛性，x 位移函数应尽可能反应物体中的真实位移，即满足完备性和连续性的要求：位移模式必须能反映单元的刚体位移。位移模式必须能反映单元的常量应变。位移模式应尽可能反应位移的连续性。

设三角形单元三个结点编号为i、j、m。平面三角形单元位

u= a 1+ a 2X+ a 3y

V= a 4+ a 5X+ a 6y

可以写成"U 0 7的形式，U 。、V o 反映了单元的刚体平动， 反

v V 。 y

映了单元的刚体转动，满足完备性和连续性的要求

采用插值法由单元结点位移列阵 8 e = U i V i U j V j U m V m T 计 算 a i 、a 2、a 3、a 4、a 5、a 6.,求出位移 d=[u (X , y), V

(x , y ) ] o 6个未知量,6个代数方程，得d e =N 8 e

列。N 为形函数矩阵，形函数 Ni 的性质有: N (X i , y i ) =1 N i (X j , y j ) =0 Ni (xm, ym ) =0

N (x , y ) +N (x , y ) +N (x , y ) =1可推出三个形函数中，两个 是独立的，反映了刚体平移。

令z=Ni ，在直接坐标系中画出 Ni 、Nj 、Nm 的函数图形是 以Ni (xi , yi ) =1为高的四面体，所以结点位移影响单元 的位移场，单元的位移场是线性分布的，相邻单元在公共边 上的位移是连续的，单元相邻边的位移只取决于单元相邻公 共边上的结点而与其他结点无关，无论以哪个单元计算相邻 边的位移，结果一定相同。

移函数选取为

式中 d e = N i 0

N j 0 N m 0 0 N i 0 N j 0 N m U i V i U j V j

T V m b i =-

C i = 1 X j y m 1 X m (i 、j 、

m 轮换)A 为三角形面积，为避免

A<0, i 、j 、m 按逆时针排 N=(a i +b i X+ey)/2A , y j y m

形函数N e 决定了单元内的位移模式， 反映了 i 结点位移对 单元内任意点位移的贡献率。

2、根据几何方程用单元结点位移表示单元应变:

b B=（ B B j B m ）T , B= o q 与x ，y 无关，所以该单元内应变是常量，反映单

元的常量 应变，满足完备性和连续性的要求 ，这是一种常应变单元。

3、根据物理方程用单元结点位移表示单元应力：

D 为弹性矩阵

DB e S e S 为应力矩阵

内x ，y 位置无关，这是一种常应力单元。

因为在三结点三角形单元中，位移函数中含有坐标的一次 项，其误差为o x 2，而应力、

应变是常量，其误差为0 X ，比位移精度低。

4、根据虚功原理用单元结点位移表示单元结点力

单元在结点处受力，单元会发生变形，因此单元在结点处 x v

y

u v B i 1 c 0 b j 0 b m 0 T 0 q 0 C j 0 C m U i V i U j V j U m V m 2A C i b i C j b j C m b m

B 为几何矩阵

c ， B

内所有元素

b S=DB 中， 每一个元素都是常数， 所以e 的每一个分量与单元 B 可写为分块矩阵 e B e

所受到的力与单元结点位移肯定有关系。单元间通过结点的

相互作用成为整体，因此每一单元的受力——位移关系找出

来，整体的受力——位移关系也就出来了

*e *

=i

因为B 、D 中元素都是常数， 单

元刚度矩阵。

K 为6行6列矩阵可写为K xx xy k j B j DB j At 二 j At ,炉表示j 结点处发生y 方向的单位位 k ij k ij

移时所引起的 i 结点处 x 方向的结点力。不同类型不同形式 的单元，只有弹性矩阵 D 和几何矩阵B 不同，计算子块矩阵 的公式相同，平面问题中，影响刚度矩阵

K 的只有几何矩阵

B 。

K 的性质有：

K 中每个元素表示个单元结点沿坐标方向发生单位位移时

所引起的结点力。

K 为对称矩阵。 B i T

k ii k ij ij k im tA B T j D B

i B j B m k ji

k j j jj k jm ， B m T k mi k mj k mm

记单元节点力为 F F i F j F m T ，单元结点虚位移为

单元内应力为

xy 单元内虚应变

x y xy 根据虚功原理，

dxdy ?t ，可得 F e B T DBdxdy t

F e B T DBtA e K e ， K=B T DBtA 为

单元做刚体位移时，单元内不产生应变应力，结点力为0, 所以K中每行每列元素之和为0,所以|K 0,所以只根据F e K e无法求得唯一解。

5、根据虚功等效原则计算等效结点力

根据有限元的基本方法，单元内任意点的位移、应变、应力等最终都要用结点位移来表示，所以作用在物体上的外力也要用结点位移表示。为了计算等效结点力，在任意的虚位移上，使原载荷与等效载荷虚功相等。

设外力为f p，结点虚位移为*e，则任意点虚位移为d*e N *e，等效节点载荷为FJ,有

d*eT f p t *eT F L e F L e N T f p t （集中力）

同理得F L e N T fds t （面力），F L e Nfdxdy t （体力）。

S A

三、整体分析

将结构的所有单元通过结点连接起来，形成一个整体的离散结构以代

替实际的连续体，以形成以结点位移为未知量的整体结构的有限元代数方程组，最后求得结点位移。对结点受力分析：结点受到与之相关的单元给它的反作用力和外载荷的等效结点力，这两组力坐标轴方向相反，

所以应该相等，即F i Fu，设有n个结点，每个结点建立两个方

e e

向的方程，不考虑外界约束时，共2n个方程，2n个未知量（ix, iy, jx...），为了建立这个代数方程组，建立整个弹性体的结点力和结点位移的关系式K F L , K（ 2n x 2n）为整理刚度

矩阵，8为整体结点位移列阵，F L整体结点载荷列阵

为了求整体刚度矩阵，要找到它与已求得的单元刚度矩阵的关系，在整体中对结点编码，设整体刚度矩阵中某元素为Kij ，意为j 个结点在x 或y 方向发生位移引起i 个结点x 或y 方向的结点力，找到同时用到i 与j 结点的单元，并用与之对应的单元刚度矩阵中的元素ksm 相加得到Kij ，整体刚度矩阵也是奇异矩阵，必须考虑边界约束条件，消除K 的奇异性，才能求解结点位移。再由单位结点等效载荷得到整体结点载荷列阵F L。这样K、F L已知，求解代数方程，解出整体结点位移列阵8 ，得到相应的单元结点位移8 e。

8 e得到了，相应的d e、a e、& e等就得到了。