平面问题的有限元法实例

有限单元法是求解连续区域内的边值问题和初值问题的数值方法。其实质是将

e

分析区连续的域V和边界S离散成为有限个、只在结点连接的子域Ve和面域 S（每个子域称为一个有限单元，各个子域连接的点称为结点），用全部有限单元的集合等价于连续域的近似分析方法。计算中只求结点上的待定函数（即：要求的函数）值，用结点上的待定函数值反映整个域待定函数的分布和变化规律。

有限单元法的分析步骤 1. 连续体的离散 2. 选择形函数 3. 单元特性分析

4. 总体特性分析（单刚－总刚） 5. 引入边界条件

6. 求解线性方程组，得到求解函数值。

下面我们就以平面应变问题的三角形单元有限元为例，阐明有限元解题的过程，从位移模式的建立，到形函数的形成以及单元刚度矩阵、总体刚度矩阵，求解等全过程。

例题：如下图是一块三角形薄板，荷载沿厚度均匀分布，荷载集度为10，板的各部分尺寸如图所示。为了计算方便，我们去板厚t=1，弹性模量E=1，泊松比

0.25，容重0。

首先我们先把下图中的薄板划分成四个三角形单元，并进行标号，标出单元号和节点码，选取坐标。

计算过程如下：

由上图可知坐标1(0,2),2(0,1),3(1,1),4(0,0),5(1,0),6(2,0). 位移边界条件：x轴上，v=0； y轴上，u=0；

荷载边界条件：由计算得1上的荷载为5，2上的荷载为10，3上的荷载为5。 所要建立的刚度矩阵为：

第 1 页 共 9 页 1

K11K 61K16K6616FL1FL6 节点荷载为FL1(5,0),FL2(10,0),FL3(0,0),FL4(0,0),FL5(0,0),FL6(0,0)

TTF（F F FFFF）（）LLLLLLL荷载阵列： T位移阵列：()(uvuvuvuvuvuv)

先求劲度矩阵k；

对于Ⅰ、Ⅱ、Ⅳ，可求得A=0.5㎡。

Axjxmxiyj(xiyixjyixmyixiymxjyixmyj)ymyi 其中

同理对于Ⅲ，A=0.5㎡。

为简单起见，我们取,tm。

设三角形单元的形函数为：(x,y)xy （1）

i(xi,yi)j(xj,yj)(x,y)mmm (2)

)ma())amija(A(aiiajjam(biibjjbm)mb()jb()bmiA)mc()jc()cmiA(ciicjjcm  (3)

Axjxmxiyiyjym得

由

aixjymxmyjajxmyixiymajxiyjxjyibiyjymbjymyibmyiyjccixmxjcjxixmcm （4） 将(3)代入(1)，

(x,y)Ni(x,y)iNj(x,y)jNm(x,y)m[N(x,y)]{(e)}

得形函数矩阵N(x,y)[Ni(x,y)Nj(x,y)Nm(x,y)]

第 2 页 共 9 页 2

N(aibixciy)xiAN(ajbjxcjy)yjANmA(ambmxcmy)xy  i{(e)}jm

设ui,j,mxi,j,myi,j,mvi,j,mxi,j,myi,j,m

uNiuiNjujNmumxuiyuj(xy)um vNiviNjvjNmvmxviyvj(xy)vm uNiuiNjujNmumdNvNvNvvjjmm ii 位移模式矩阵

e dN

eT()u(iviujvjumvm )ijm

NiN应变：

NiNjNjNmxyxyNmxyxy

uiviujumucmjvjvmvjbmuvivmujummvm

dudxbidvdyAcidvdudxdycibibjcjcjbjbmcmBe 由B*，则

应力：由D

第 3 页 共 9 页 3

ED 其中



D*BeuiEviujvjumvmuiviEujvjumvm

(()e)TFeT劲度矩阵：由虚功原理

A()dxdyt

其中

()T(B()e)T(()e)TBT 整理得

(()e)TFe(()T)(ABTDBdxdyt)e

FeABTDBdxdyteFeke(adxdyA.m)kBTDBtABTDBt

代入得：

kiikijkimkkkkEtjijjjmkk()mikmjmm取,tm时，得

第4

4 页 共 9 页 

....E....k........

eTeeF(FFF)ijmFk知道了k，我们由，即可求出节点力

Krskrse，这里整体编码与局部编码之间存在下列关系：

kiikkjikmikijkjjkmjkimkjmkmmFFkFFFkF

FFkFFFkFⅠⅠⅠⅠⅠⅠⅠⅠ KkkKkkKkkKkk Fjjjmji先看，有

ⅠⅠⅠⅢⅡⅠⅡⅢKkkKkkkkkkFmjmmjjii ，有

同理有F,F,F,F的矩阵元素，如下：

第 5 页 共 9 页 5

kⅠjjⅠkmjkⅠKijkⅠjmⅡⅢkⅠmmkjjkiiⅢkⅠimkmikⅡmjⅢkⅡijkjikⅠjiⅢkⅠmikimⅢⅣkⅠiikmmkjjⅣkⅢjmkmjⅣkijkⅡjmkⅡmmkⅡjmⅢkⅡjikijⅢkmjkⅣjmⅡkmiⅡⅣkiikⅢjjkmmⅣkimkⅣjiⅣkmiⅣkii

...............E...........

整体平衡方程组为：

..................

uv.uv..uv.u.vu...vuv...................E.................

可求得节点位移：

第 6 页 共 9 页 6

再由位移求应力 对于Ⅰ、Ⅱ、Ⅳ单元

u.v.u.v.u.v.uEvu.v.u.v e由应力转换矩阵S，及DDB有

D*BSE....

同理对于Ⅲ单元

D*BSE....

则各单元的应力为

uvxEyv....xyv Ⅰ：

uxyE....vxy Ⅱ：





第 7 页 共 9 页 7

vxuEy....xyuv Ⅲ：

uxuyE....vuxy Ⅳ：





我们编程要输入的数据如下： 1. 基本参数

1) 单元数NE=； 2) 节点数NJ=； 3) 支承数NZ=；

4) 节点荷载数NPJ=； 5) 半带宽DD=；

6) 节点位移数NJ2=NJ*2=； 2. 其他参数

1) 问题类型码LXM，平面应力LXM=0；

2) 弹性常数E，，弹性模量EO=1，泊松比MU=0.25； 3) 容重，LOU=0； 4) 板厚t，TE=1； 5) 节点坐标数组AJZ

021012113AJZ(N,J2)004105206

AJZ(NJ,2)中的元素按节点整体码顺序输入，数组行号为节点整体码，每个节点的坐标值存一行，第一列存x值，第二列存y值。 6) 节点码数组JM

第 8 页 共 9 页 8

123Ⅰ245ⅡJM(NE,3)253Ⅲ356Ⅳ

JM（NE，3）中的元素按单元输入，单元号为行号。每个单元的整体码存一行。局部码为列号，局部码对应的整体码放在相应的列。 7) 支承数组NZC

7NZC(NZ)812

NZC(NZ)中的元素按支承对应的位移数从小到大排列。 8) 节点荷载数组PJ

57PJ(NPJ,2)108512

节点荷载数 荷载对应的位移数

PJ(NPJ，2)中的元素排列原则为：一个节点荷载存一行，同一行中的第一列是节点荷载值，第二列是荷载对应的位移数。 3. 程序框图

程序如附件所示。

第 9 页 共 9 页 9