第五章 三角形单元的有限元法
有限元方法
§7. 两点边值问题的有限元方法
本节以两点边值问题为例，并从Ritz法和Galerkin法两 种观点出发来叙述有限元法的基本思想及解题过程.
7.1 基于Ritz法的有限元方程 7.2 基于Galerkin法的有限元方程
这样，我们就得到了单元有限元特征式的一般表示形式：
K(i)u(i) F(i)
第二步：总体合成.总体合成就是将单元上的有限元特征 式进行累加，合成为总体有限元方程. 这一过程实际上是将 单元有限元特征式中的系数矩阵（称为单元刚度矩阵）逐个 累加，合成为总体系数矩阵（称为总刚度矩阵）；同时将右 端单元荷载向量逐个累加，合成为总荷载向量，从而得到关 于的线性代数方程组．为此，记
于是有 u(i) (ui1,ui)TB (i)u
从而式(7.16)右端第一个和式为
1 nu iT K iu i 1 nu T [ ( B i) T K iB i] u 1 u T K u ，
2 i 1
2 i 1
2
其中
（未标明的元素均为0）这就是总刚度矩阵． 对式(7.16)右端第二个和式，有
其中，p x C 1 a , b , p 0 , q C a , b , q 0 , f C a , b
精选版课件ppt
3
1. 写出Ritz形式的变分问题
与边值问题(7.1)、(7.2)等价的变分问题是：
求
u*
H
1，使
E
其中，
Ju*m uH in1 EJu J u 1 a u ,u f,u
u j
便得到确定 u1,u2,
,un的线性代数方程组
有限元分析及工程应用-2016第五章
5.1 轴对称问题有限单元法
机械学院
（1）三角形截面环形单元 1）位移模式
qe ui wi u j wj uk wk T
与平面三角形单元相似，仍选取线 性位移模式，即：
u w
a1 a4
a2r a5r
aa36zz
u Niui N ju j Nkuk
,
A2
1 2 2(1 )
单元中除了剪应力外其 它应力分量也不是常量
在轴对称情况下，由虚功原理可推导出单元刚度矩阵
K e VBT DBddrdz 2 BT DBrdrdz
5.1 轴对称问题有限单元法
机械学院
（1）三角形截面环形单元
2）单元刚度矩阵
K e VBT DBddrdz
Loads>Apply>Structural>Displacement>Symmetry B.C.>On Lines，用鼠标在图形窗口上拾取编号为“1”和“3”的线段 ，单击[OK]，就会在这两条线上显示一个“S”的标记，即 为对称约束条件。
(7)施加面力：Main Menu>Solution>Define Loads>Apply>Structural>Pressure>On Lines，用鼠标在图形 窗口上拾取编号为“4”，单击[OK] 在“VALUE Load PRES value”后面的输入框中输入“10”，然后单击[OK]即可
5.1 轴对称问题有限单元法
机械学院
（3）应用实例 (3)建立几何模型：
MainMenu>Preprocessor>Modeling>Create>Areas>Rectangle>By Dimension，在出现的对话框中分别输入：X1=5，X2=10，Y1=0， Y2=20，单击[OK]。
有限元法原理
有限元法原理
有限元法是一种工程计算方法，主要用于求解连续介质的力学问题。
它的基本原理是将连续介质离散成有限个小单元，然后利用有限元的形状函数对每个小单元进行近似，最终利用这些近似解来求解整个连续介质的力学问题。
有限元法的主要思想是将问题的解表示为一个有限个数的基函数的线性组合。
这些基函数与小单元的形状函数相联系，通过对小单元的形状函数进行合适的选取和调整，可以确保解在小单元内满足边界条件。
然后，通过将所有的小单元的解进行组合，就可以得到整个连续介质的解。
在实际的计算中，有限元法通常分为以下几个步骤：首先，需要根据实际问题确定合适的有限元模型，包括选择适当数量和类型的有限元单元。
然后，需要确定边界条件，即确定整个连续介质的边界约束条件。
接下来，根据小单元的形状函数和基函数，可以建立刚度矩阵和荷载向量。
最后，通过求解线性方程组，可以得到整个连续介质的解。
有限元法具有广泛的应用范围，在工程领域中可以用于求解各种静力学、动力学、热力学、流体力学等问题。
它不仅能够提供精确的解，同时也具有较高的计算效率和灵活性。
因此，有限元法已经成为工程计算领域中一种非常重要的数值分析方法。
有限元分析第五章(第一部分)
第五章 等(Isoparametric Elements)在前面的章节中我们已经认识了三角形单元和矩形单元。
这两种单元的边均为直边，用直边单元离散曲边的求解域势必要用更多的单元数才能较准确地描述实际边界。
本章将要介绍的等参数单元是目前应用最广的一类单元，可用这类单元更精确的描述不规则的边界。
这类单元的出现不仅系统的解决了构造协调位移单元的问题，而且自然坐标系的描述方法也广泛为其他类型的单元所采用。
等参数单元在构造形函数时首先定义一个规则的母体单元（参考单元），在母体单元上构造形函数，再通过等参数变换将实际单元与母体单元联系起来。
变换涉及两个方面：几何图形的变换（坐标变换）和位移场函数的变换，由于两种变换采用了相同的函数关系（形函数）和同一组结点参数，故称其为等参数变换。
§5-1四结点四边形等参数单元１、母体单元 自然坐标和形函数母体单元ê ：边长为２的正方形，自然坐标系ξ，η 示于图5-1。
取四个角点为结点，在单元内的排序为１、２、３、４。
仿照矩形单元，可定义出四个形函数显然有如下特点：（i ）是ξ，η的双线性函数 （ii ）(iii)２、实际单元与母体单元之间的坐标变换（１） 坐标变换设xy 平面上的实际单元e 由母体单元经过变换F 得到，即 且规定结点（ξi ，ηi ）与结点（x i , y i ）对应（i =1~4）。
这样的变换不只一个，利用（5-1-1）定义的形函数即可写出这种变换中的一个１图5-1 ())4~1()1(141),(=++=i N i i i ηηξξηξ),(ηξi N ⎩⎨⎧=≠=i j i i N ij i 当 当 ＝10),(δηξ),(ηξi N 1)1)(1(41)1)(1(41)1)(1(41)1)(1(41),(41≡+-++++-++--=∑=ηξηξηξηξηξi i N e e F →： (5-1-2) (5-1-1) ii i i i i y N y x N x ⋅=⋅=∑∑==4141),(),(ηξηξ(5-1-3)（5-1-3）所定义的变换有如下特点：x , y 是ξ，η的双线性函数。
三角形单元数值积分
三角形单元数值积分一、引言数值积分是数值分析中的一个重要内容，它是利用数值方法来近似计算定积分的过程。
在实际应用中，很多函数都无法求出其解析式，因此需要采用数值积分方法来进行近似计算。
本文将重点介绍三角形单元数值积分的相关知识。
二、三角形单元三角形单元是有限元方法中最基本的单元之一，它由三个节点构成。
在实际应用中，我们通常采用局部坐标系来描述三角形单元。
假设三角形的三个顶点为A、B、C，则可以定义局部坐标系x-y为：以AB边为x轴正方向，以C点到AB边垂线为y轴正方向。
三、三角形单元上的积分对于一个在三角形上定义的函数f(x,y)，我们需要对其进行积分。
根据高斯公式，可以将二维平面上任意闭合曲线内部的积分转化为该曲线上的积分。
因此，在三角形内部进行二重积分时，可以将其转化为对该三角形边界上的积分。
四、高斯公式高斯公式是将一个闭合曲线内部的积分转化为该曲线上的积分的公式。
对于一个在平面区域D上连续可微的函数f(x,y)，高斯公式可以表示为：∬D(∂Q/∂x-∂P/∂y)dxdy=∮C(Pdx+Qdy)其中，P和Q是f(x,y)的偏导数，C为D的边界曲线。
五、三角形单元数值积分在实际应用中，我们需要采用数值方法来进行三角形单元上的积分计算。
常见的数值积分方法有梯形法、辛普森法、高斯积分法等。
其中，高斯积分法是一种比较常用和精确的方法。
六、高斯积分法高斯积分法是一种通过求解一组带权重系数和节点坐标的代数方程组来近似计算定积分的方法。
在三角形单元上进行高斯积分时，我们通常需要将其转化为在标准三角形（即顶点坐标为(0,0)、(1,0)、(0,1)）上进行计算。
七、标准三角形上的高斯积分对于一个定义在标准三角形上的函数f(x,y)，可以采用如下公式进行高斯积分：∫∫f(x,y)dxdy=∑wi*f(xi,yi)其中，wi为权重系数，(xi,yi)为高斯积分点的坐标。
在实际应用中，通常采用2-3-4-5阶高斯积分公式进行计算。
有限元平面问题三角形实例
有限元平面问题三角形实例有限元法是一种常用的计算方法，可以用来解决各种工程问题。
其中，有限元平面问题是有限元法的一种应用，常用于分析三角形结构。
在有限元平面问题中，我们通常会将结构划分成许多小的单元，每个单元由节点和单元刚度矩阵组成。
而三角形结构则是有限元平面问题中常用的一种单元形状。
三角形结构的特点是简单而且易于处理，因此广泛应用于各种领域，如土木工程、机械工程、航空航天等。
下面我们就以一个实际的例子来说明如何应用有限元平面问题分析三角形结构。
假设我们要分析一个三角形钢板在受力作用下的变形情况。
首先，我们需要将钢板划分为许多小的三角形单元。
每个单元由三个节点组成，节点之间通过边连接。
在有限元分析中，我们需要对每个单元进行网格划分，并确定节点的坐标和边的长度。
然后，通过求解节点的位移和应力分布，可以得到钢板在受力作用下的变形情况。
具体来说，我们可以通过求解线性方程组来得到节点的位移。
而节点的应力则可以通过应变-位移关系来计算。
通过这种方式，我们可以得到钢板在受力作用下各个节点的位移和应力分布情况。
有限元平面问题的分析结果可以帮助我们了解结构的强度和刚度情况，为设计和优化提供依据。
例如，在钢板的设计中，我们可以通过有限元分析来确定合适的材料和尺寸，以满足结构的强度和刚度要求。
除了钢板，有限元平面问题还可以应用于其他类型的三角形结构。
例如，在土木工程中，我们可以使用有限元分析来分析三角形桥梁或者三角形支撑结构的变形和应力分布情况。
有限元平面问题是一种常用的分析方法，可以应用于各种三角形结构的分析。
通过对节点的位移和应力分布的求解，我们可以得到结构在受力作用下的变形情况。
这对于工程设计和优化至关重要，可以帮助我们提高结构的强度和刚度，确保其安全可靠。
第五章 三角形单元的有限元法
{ } [ x y xy ]T （1-4）
（1-6）
7、物理方程矩阵式 对于弹性力学的平面应力问题，物理方程的矩阵形 式可表示为：
x E y 2 1 xy 1 0 对 x 1 称 y 1 xy 0 2
qVx T {qV } [qVx qVy ] qVy
（1-2）
3、单元内任意点的位移列阵f
{ f } [u
y
i m
]T
y
m
（1-3） ·
i
·u
v
j
j
x
x
图1-1
4、单元内任意点的应变列阵
{ } [ x y xy ]T
（1-4）
[ N ] [[Ni ] [ N j ] [ Nm ]]
(1-21）
其中子矩阵
Ni 0 [ Ni ] Ni [ I ] 0 Ni [I]是2×2的单位矩阵。
(i, j, m) (1-22）
（2）形函数性质 形函数是有限单元法中的一个重要函数，它具 有以下性质： 性质1 形函数Ni在节点i上的值等于1，在其它节点 上的值等于0。对于本单元，有
用形函数把式(1-16）写成矩阵，有
u N i v 0
缩写为
0 Ni
Nj 0
0 Nj
Nm 0
ui v i 0 u j Nm v j um vm
{ f } [ N ]{ }
(1-20）
[N]为形函数矩阵，写成分块形式：
第五章 三角形单元的有限元法
5.1 基本思想
把整体结构离散为有限个单元，研究单元的平衡和变 形协调；再把这有限个离散单元集合还原成结构，研究离 散结构的平衡和变形协调。划分的单元大小和数目根据计 算精度和计算机能力来确定。
有限元方法-第五章--平面三角形单元
D
E
1 2
1
0
对 1 0
称
1
(i)
2
所以，[S]的子矩阵可记为
Si DBi
E
2 1 2
bi
1
bi
2
ci
ci
1
ci
2
bi
（ i
,
j
,
m轮换） (5-19)
对于平面应变问题，只要将 (i) 式中的E换成E/1-2 ， 换成 /1-，即得到其弹性矩阵
D
1
E1 1 2
1
1
起来，便可近似地表示整个区域的真实位移函数。这种 化繁为简、联合局部逼近整体的思想，正是有限单元法 的绝妙之处。
基于上述思想，我们可以选择一个单元位移模式，
单元内各点的位移可按此位移模式由单元节点位移通过
插值而获得。线性函数是一种最简单的单元位移模式，
故设
u 1 2x 3y
v 4 5x 6y
(b)
0
(b)
Ni xm
,
ym
1 2
ai
bi xm
ci ym
0
(c)
类似地有
N j xi , yi 0 , N j x j , y j 1 , N j xm , ym 0 (d) Nm xi , yi 0 , Nm x j , y j 0 , Nm xm , ym 1
式中 1、2、…6是待定常数。因三角形单元共有六个
自由度，且位移函数u、v在三个节点处的数值应该等于 这些点处的位移分量的数值。假设节点i、j、m的坐标分 别为（xi , yi ）、（xj , yj ）、（xm , ym ），代入 (b) 式， 得：
ui 1 2 xi 3 yi
数值模拟：第五讲 平面问题(二)——三角形单元分析
2) 单元刚度方程和单元刚度矩阵的建立是单元分析的核心内容。
3) 一般情况下，单元应变矩阵是坐标的函数矩阵，所以单元刚度矩 阵的计算需要进行积分运算。
4) 所建立的单元刚度矩阵反映了一般弹性体小单元近似的弹性性质， 是单元特性的核心。
.
• 单元刚度矩阵的计算
➢ 弹性力学平面问题的单元刚度矩阵通式：Fra bibliotekllm
s2llm
3）形函数在单元上的积分：
Ni(x,y)dxdy
A 3
(i l,m,n)
.
5.2.4 单元应变和应力
• 已知节点位移插值形式的单元位移模式：
u v
N
e
• 代入平面问题几何方程（应变~位移关系）得到单元应变：
xxyy 0x
0
y uv 0x
0 yN 0l
0 Nl
（简单三角形单元的形函数只有2个独立）
.
➢ 性质3（推论）：简单三角形单元的形函数在边界上的性质。 某节点的形函数在该点邻边上呈线性分布，取值在0~1之间， 在该点对边上值为零。
简单三角形单元形函数的几何意义
❖ 由形函数表达式和性质1可画出下列形函数几何图形。
.
❖ 根据位移模式表达式及其形函数的性质，可以推断出两个相邻 三角形单元上位移分布形状和公共边界上位移的情况：
.
• 针对三节点三角形单元，可以导出单元形函数 的下列性质。
➢ 性质1：单元上某节点的形函数在该节点的值为1，在其它节点 的值为零。
N l ( xl , yl ) 1 N l (xm , ym ) 0 N l (xn, yn ) 0
(l，m，n)
➢ 性质2：单元上所有形函数之和等于1。
Nl NmNn1
有限元基础第五章-线性三角形单元
x3 ) y
N3
1 2A
x1 y2
x2
y1
( y1
y2 )x
( x2
x1) y
其中A是三角形的面积
1 A 1 1
2
x1 x2
y1 y2
1 x3 y3
17:57
5
平面三角形单元
u x, y N1u1 N2u2 N3u3
同理
v x, y N1v1 N2v2 N3v3
u(x, y) N(x, y)de
17:57
6
平面三角形单元
三角形的形函数可统一表示为：
Ni
1 2A
ai
bi x
ci y
其中
ai
x
j
ym
xm
y
j
xj xm
yj ym
1 bi yi ym 1
yj ym
k
j i
1 ci xm x j 1
xj xm
1 xi yi 2A 1 xj yj
1 xm ym
k k = 3, 1, 2
i = 1, 2, 3 i
j j = 2, 3, 1
17:57
7
形函数的性质
在单元任一点上三个形函数之和等于1(单位分解性)
Ni (x, y) N j (x, y) Nm (x, y)
1 2A
(ai
a
j
am )
(bi
bj
bm )x
(ci
cj
cm
) y
1
第一列与它的 第一列与第二
代数余子式乘 列的代数余子
积之和
式乘积之和
17:57
10
形函数的性质
Ni (x,
弹塑性力学及有限元法_
写成矩阵形式
R11 cos 2 θ x 1 Ry1 EA cos θ sin θ 1 = Rx 2 l1 − cos 2 θ R1 2 − cos θ sin θ y cos θ sin θ sin 2 θ − cos θ sin θ − sin 2 θ − cos 2 θ − cos θ sin θ cos 2 θ cos θ sin θ
单元刚度矩阵的子矩阵 K ij 表示：当单元 e 中节点 j 取单 位位移，且其它节点位移为零时，对应于 i 节点的节点力。
第五章 有限元法简介
单元1的节点力和节点位移的关系可写成
R1 K11 = R2 K 21
1
K12 K 22
1
δ1 δ 2
1 θFx1(u1) 3 Fx3 (u3) Fy1(v1 ) Fy3 (v3) y 2 o x
1
Fy2 (v2) Fx2(u2)
2
图5-1 简例结构图
第五章
分析步骤：
有限元法简介
2
1
1 1 Ry2(v2) 1 1 Rx2(u2)
1. 离散结构物为有限个单元 分为2个单元，第一个单元的节点编号 为1和2，第二个单元的节点编号为2和3。 对于第一单元，在第1、2节点处的节点力 为 R 11 , R 11 , R 1 2 , R 1 2 ，表示节点施加在单元1上 x y x y
1 − cos θ sin θ u1 1 2 − sin θ v1 cos θ sin θ u1 2 1 si成
R11 k x 1 11 Ry1 k21 1 = Rx 2 k31 R1 k41 y2 k12 k22 k32 k42 k13 k23 k33 k43
有限元作业：三角形单元求解
《有限元作业》年级2015级学院机电工程学院专业名称班级学号学生2016年05月如下图所示为一受集中力P作用的结构，弹性模量E为常量，泊松比V=1/6,厚度为I=1。
按平面应力问题计算，运用有限元方法，分别采用三角形及四边形单元求解，求节点位移及单元应力（要求三角形单元数量不少于4个，四边形单元不少于2个）图(一)图（二）三角形单元求解图（三）四边形单元求解（1）如图划分三角形单元，工分成四个分别为④（2）如图分别进行编号1、2、3、4、5、6，并建立坐标系（3）编程进行求解，得出结果，其中假设力P=2000N调用Triangle2D3Node_Stiffness函数，求出单元刚度矩阵k1 =1.0e+06 *7.2857 -3.0000 -2.1429 0.8571 -5.1429 2.1429-3.0000 7.2857 2.1429 -5.1429 0.8571 -2.1429 -2.1429 2.1429 2.1429 0 0 -2.14290.8571 -5.1429 0 5.1429 -0.8571 0-5.1429 0.8571 0 -0.8571 5.1429 02.1429 -2.1429 -2.1429 0 0 2.1429k2 =1.0e+06 *5.1429 0 -5.1429 0.8571 0 -0.85710 2.1429 2.1429 -2.1429 -2.1429 0-5.1429 2.1429 7.2857 -3.0000 -2.1429 0.85710.8571 -2.1429 -3.0000 7.2857 2.1429 -5.14290 -2.1429 -2.1429 2.1429 2.1429 0-0.8571 0 0.8571 -5.1429 0 5.1429 k3 =1.0e+06 *2.1429 0 -2.1429 -2.1429 0 2.14290 5.1429 -0.8571 -5.1429 0.8571 0-2.1429 -0.8571 7.2857 3.0000 -5.1429 -2.1429 -2.1429 -5.1429 3.0000 7.2857 -0.8571 -2.14290 0.8571 -5.1429 -0.8571 5.1429 02.1429 0 -2.1429 -2.1429 0 2.1429 k4 =1.0e+06 *2.1429 0 -2.1429 -2.1429 0 2.14290 5.1429 -0.8571 -5.1429 0.8571 0-2.1429 -0.8571 7.2857 3.0000 -5.1429 -2.1429 -2.1429 -5.1429 3.0000 7.2857 -0.8571 -2.14290 0.8571 -5.1429 -0.8571 5.1429 02.1429 0 -2.1429 -2.1429 0 2.1429 调用Triangle2D3Node_Assembly函数，求出总体刚度矩阵求出的节点位移U =-0.00040.00080.00050.00100.00070.0023-0.00070.0026调用Triangle2D3Node_Stress函数，求出应力,S1、S2、S3、中求出的分别为Sx,Sy,SxyS1 =1.0e+03 *-4.4086-0.73483.5914S2 =1.0e+03 *4.4086-0.64050.4086S3 =1.0e+03 *1.8907-1.06012.1093S4 =1.0e+03 *-1.89072.10931.8907二、（1）如图划分四边形单元，工分成四个分别为（2）如图分别进行编号1、2、3、4、5、6，并建立坐标系（3）编程进行求解，得出结果，其中假设力P=2000N调用Quad2D4Node_Stiffness函数，求出单元刚度矩阵调用Quad2D4Node_Assembly函数，求出求出总体刚度矩阵求出节点位移U =0.00120.0017-0.00120.00170.00160.0049-0.00170.0052调用Quad2D4Node_Stress函数，求出单元应力中的的S1、S2、S3分别为Sx,Sy,Sxy应力分量S1 =1.0e+03 *0.0000-0.24782.0000S2 =1.0e+07 *0.68564.1135-1.7137程序附录一、1、三角形单元总程序：E=1e7;NU=1/6;t=1;ID=1;%调用Triangle2D3Node_Stiffness函数，求出单元刚度矩阵k1=Triangle2D3Node_Stiffness(E,NU,t,0,1,0,0,1,1,ID)k2=Triangle2D3Node_Stiffness(E,NU,t,0,0,1,0,1,1,ID)k3=Triangle2D3Node_Stiffness(E,NU,t,1,1,1,0,2,0,ID)k4=Triangle2D3Node_Stiffness(E,NU,t,2,0,2,1,1,1,ID)%调用Triangle2D3Node_Assembly函数，求出总体刚度矩阵KK = zeros(12,12);KK=Triangle2D3Node_Assembly(KK,k1,1,2,3);KK=Triangle2D3Node_Assembly(KK,k2,2,4,3);KK=Triangle2D3Node_Assembly(KK,k3,3,4,5);KK=Triangle2D3Node_Assembly(KK,k4,5,6,3)% 边界条件的处理及刚度方程求解k=KK(5:12,5:12)p=[0;0;0;0;0;0;0;2000]u=k\p%支反力的计算U=[0;0;0;0;u] %为节点位移P=KK*U%调用Triangle2D3Node_Strain函数，求出应变SN1、SN2、SN3中求出的分别为SNx,SNy,SNxyu1=[U(1);U(2);U(3);U(4);U(5);U(6)];u2=[U(3);U(4);U(7);U(8);U(5);U(6)];u3=[U(5);U(6);U(7);U(8);U(9);U(10)];u4=[U(9);U(10);U(11);U(12);U(5);U(6)];SN1=Triangle2D3Node_Strain(0,1,0,0,1,1,u1)SN2=Triangle2D3Node_Strain(0,0,1,0,1,1,u2)SN3=Triangle2D3Node_Strain(1,1,1,0,2,0,u3)SN4=Triangle2D3Node_Strain(2,0,2,1,1,1,u4)%调用Triangle2D3Node_Stress函数，求出应力,S1、S2、S3、中求出的分别为Sx,Sy,Sxyu1=[U(1);U(2);U(3);U(4);U(5);U(6)];u2=[U(3);U(4);U(7);U(8);U(5);U(6)];u3=[U(5);U(6);U(7);U(8);U(9);U(10)];u4=[U(9);U(10);U(11);U(12);U(5);U(6)];S1=Triangle2D3Node_Stress(E,NU,0,1,0,0,1,1,u1,ID)S2=Triangle2D3Node_Stress(E,NU,0,0,1,0,1,1,u2,ID)S3=Triangle2D3Node_Stress(E,NU,1,1,1,0,2,0,u3,ID)S4=Triangle2D3Node_Stress(E,NU,2,0,2,1,1,1,u4,ID)2、求刚度矩阵程序function k=Triangle2D3Node_Stiffness(E,NU,t,xi,yi,xj,yj,xm,ym,ID)%该函数计算单元的刚度矩阵%输入弹性模量E，泊松比NU，厚度t%输入三个节点i、j、m的坐标xi,yi,xj,yj,xm,ym%输入平面问题性质指示参数ID(1为平面应力，2为平面应变)%输出单元刚度矩阵k(6X6)%---------------------------------------------------------------A = (xi*(yj-ym) + xj*(ym-yi) + xm*(yi-yj))/2;betai = yj-ym;betaj = ym-yi;betam = yi-yj;gammai = xm-xj;gammaj = xi-xm;gammam = xj-xi;B = [betai 0 betaj 0 betam 0 ;0 gammai 0 gammaj 0 gammam ;gammai betai gammaj betaj gammam betam]/(2*A);if ID == 1D = (E/(1-NU*NU))*[1 NU 0 ; NU 1 0 ; 0 0 (1-NU)/2];elseif ID == 2D = (E/(1+NU)/(1-2*NU))*[1-NU NU 0 ; NU 1-NU 0 ; 0 0 (1-2*NU)/2]; endk= t*A*B'*D*B;3、求整体刚度矩阵function z = Triangle2D3Node_Assembly(KK,k,i,j,m)%该函数进行单元刚度矩阵的组装%输入单元刚度矩阵k%输入单元的节点编号I、j、m%输出整体刚度矩阵KK%---------------------------------------------------------------DOF(1)=2*i-1;DOF(2)=2*i;DOF(3)=2*j-1;DOF(4)=2*j;DOF(5)=2*m-1;DOF(6)=2*m;for n1=1:6for n2=1:6KK(DOF(n1),DOF(n2))= KK(DOF(n1),DOF(n2))+k(n1,n2);endendz=KK;4、求应变程序function strain=Triangle2D3Node_Strain(xi,yi,xj,yj,xm,ym,u)%该函数计算单元的应变%输入三个节点i、j、m的坐标xi,yi,xj,yj,xm,ym%输入单元的位移列阵u(6X1)%输出单元的应力strain(3X1)，由于它为常应变单元，则单元的应变分量为SNx,SNy,SNz%---------------------------------------------------------------A = (xi*(yj-ym) + xj*(ym-yi) + xm*(yi-yj))/2;betai = yj-ym;betaj = ym-yi;betam = yi-yj;gammaj = xi-xm;gammam = xj-xi;B = [betai 0 betaj 0 betam 0 ;0 gammai 0 gammaj 0 gammam ;gammai betai gammaj betaj gammam betam]/(2*A);strain = B*u;5、求应力程序function stress=Triangle2D3Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,u,ID)%该函数计算单元的应力%输入弹性模量E，泊松比NU，厚度t%输入三个节点i、j、m的坐标xi,yi,xj,yj,xm,ym%输入平面问题性质指示参数ID(1为平面应力，2为平面应变)，单元的位移列阵u(6X1)%输出单元的应力stress(3X1)，由于它为常应力单元，则单元的应力分量为Sx,Sy,Sxy%---------------------------------------------------------------A = (xi*(yj-ym) + xj*(ym-yi) + xm*(yi-yj))/2;betai = yj-ym;betaj = ym-yi;betam = yi-yj;gammai = xm-xj;gammam = xj-xi;B = [betai 0 betaj 0 betam 0 ;0 gammai 0 gammaj 0 gammam ;gammai betai gammaj betaj gammam betam]/(2*A);if ID == 1D = (E/(1-NU*NU))*[1 NU 0 ; NU 1 0 ; 0 0 (1-NU)/2];elseif ID == 2D = (E/(1+NU)/(1-2*NU))*[1-NU NU 0 ; NU 1-NU 0 ; 0 0 (1-2*NU)/2]; endstress = D*B*u;二、1、四边形单元总程序：E=1e7;NU=1/6;h=1;ID=1;%调用Quad2D4Node_Stiffness函数，求出单元刚度矩阵k1= Quad2D4Node_Stiffness(E,NU,h,0,1,0,0,1,0,1,1,ID)k2= Quad2D4Node_Stiffness(E,NU,h,1,0,2,0,2,1,1,1,ID)%调用Quad2D4Node_Assembly函数，求出求出总体刚度矩阵KK=zeros(12,12);KK= Quad2D4Node_Assembly(KK,k1,1,2,3,4);KK= Quad2D4Node_Assembly(KK,k2,3,5,6,4)% 边界条件的处理及刚度方程求解k=KK(5:12,5:12)p=[0;0;0;0;0;0;0;2000]u=k\p%支反力的计算U=[0;0;0;0;u] %为节点位移P=KK*U%调用Quad2D4Node_Stress函数，求出单元应力中的的S1、S2、S3分别为Sx,Sy,Sxy应力分量u1=[U(1);U(2);U(3);U(4);U(5);U(6);U(7);U(8)];u2=[U(5);U(6);U(9);U(10);U(11);U(12);U(7);(8)];S1= Quad2D4Node_Stress(E,NU,0,1,0,0,1,0,1,1,u1,ID)S2= Quad2D4Node_Stress(E,NU,1,0,2,0,2,1,1,1,u2,ID)2、求刚度矩阵程序function k= Quad2D4Node_Stiffness(E,NU,h,xi,yi,xj,yj,xm,ym,xp,yp,ID) %该函数计算单元的刚度矩阵%输入弹性模量E，泊松比NU，厚度h%输入4个节点i、j、m、p的坐标xi,yi,xj,yj,xm,ym,xp,yp%输入平面问题性质指示参数ID(1为平面应力，2为平面应变)%输出单元刚度矩阵k(8X8)%---------------------------------------------------------------syms s t;a = (yi*(s-1)+yj*(-1-s)+ym*(1+s)+yp*(1-s))/4;b = (yi*(t-1)+yj*(1-t)+ym*(1+t)+yp*(-1-t))/4;c = (xi*(t-1)+xj*(1-t)+xm*(1+t)+xp*(-1-t))/4;d = (xi*(s-1)+xj*(-1-s)+xm*(1+s)+xp*(1-s))/4;B1 = [a*(t-1)/4-b*(s-1)/4 0 ; 0 c*(s-1)/4-d*(t-1)/4 ;c*(s-1)/4-d*(t-1)/4 a*(t-1)/4-b*(s-1)/4];B2 = [a*(1-t)/4-b*(-1-s)/4 0 ; 0 c*(-1-s)/4-d*(1-t)/4 ;c*(-1-s)/4-d*(1-t)/4 a*(1-t)/4-b*(-1-s)/4];B3 = [a*(t+1)/4-b*(s+1)/4 0 ; 0 c*(s+1)/4-d*(t+1)/4 ;c*(s+1)/4-d*(t+1)/4 a*(t+1)/4-b*(s+1)/4];B4 = [a*(-1-t)/4-b*(1-s)/4 0 ; 0 c*(1-s)/4-d*(-1-t)/4 ;c*(1-s)/4-d*(-1-t)/4 a*(-1-t)/4-b*(1-s)/4];Bfirst = [B1 B2 B3 B4];Jfirst = [0 1-t t-s s-1 ; t-1 0 s+1 -s-t ;s-t -s-1 0 t+1 ; 1-s s+t -t-1 0];J = [xi xj xm xp]*Jfirst*[yi ; yj ; ym ; yp]/8;B = Bfirst/J;if ID == 1D = (E/(1-NU*NU))*[1 NU 0 ; NU 1 0 ; 0 0 (1-NU)/2];elseif ID == 2D = (E/(1+NU)/(1-2*NU))*[1-NU NU 0 ; NU 1-NU 0 ; 0 0 (1-2*NU)/2]; endBD = J*transpose(B)*D*B;r = int(int(BD, t, -1, 1), s, -1, 1);z = h*r;k = double(z);3、求总体刚度矩阵程序function z = Quad2D4Node_Assembly(KK,k,i,j,m,p)%该函数进行单元刚度矩阵的组装%输入单元刚度矩阵k，单元的节点编号i、j、m、p%输出整体刚度矩阵KK%---------------------------------------------------------------DOF(1)=2*i-1;DOF(2)=2*i;DOF(3)=2*j-1;DOF(4)=2*j;DOF(5)=2*m-1;DOF(6)=2*m;DOF(7)=2*p-1;DOF(8)=2*p;for n1=1:8for n2=1:8KK(DOF(n1),DOF(n2))= KK(DOF(n1),DOF(n2))+k(n1,n2);endendz=KK;4、求应力程序function stress= Quad2D4Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,xp,yp,u,ID) %该函数计算单元的应力%输入弹性模量E，泊松比NU，厚度h，%输入4个节点i、j、m、p的坐标xi,yi,xj,yj,xm,ym,xp,yp,%输入平面问题性质指示参数ID(1为平面应力，2为平面应变)%输入单元的位移列阵u(8X1)%输出单元的应力stress(3X1)%由于它为常应力单元，则单元的应力分量为Sx,Sy,Sxy%---------------------------------------------------------------syms s t;a = (yi*(s-1)+yj*(-1-s)+ym*(1+s)+yp*(1-s))/4;b = (yi*(t-1)+yj*(1-t)+ym*(1+t)+yp*(-1-t))/4;c = (xi*(t-1)+xj*(1-t)+xm*(1+t)+xp*(-1-t))/4;d = (xi*(s-1)+xj*(-1-s)+xm*(1+s)+xp*(1-s))/4;B1 = [a*(t-1)/4-b*(s-1)/4 0 ; 0 c*(s-1)/4-d*(t-1)/4 ;c*(s-1)/4-d*(t-1)/4 a*(t-1)/4-b*(s-1)/4];B2 = [a*(1-t)/4-b*(-1-s)/4 0 ; 0 c*(-1-s)/4-d*(1-t)/4 ;c*(-1-s)/4-d*(1-t)/4 a*(1-t)/4-b*(-1-s)/4];B3 = [a*(t+1)/4-b*(s+1)/4 0 ; 0 c*(s+1)/4-d*(t+1)/4 ;c*(s+1)/4-d*(t+1)/4 a*(t+1)/4-b*(s+1)/4];B4 = [a*(-1-t)/4-b*(1-s)/4 0 ; 0 c*(1-s)/4-d*(-1-t)/4 ;c*(1-s)/4-d*(-1-t)/4 a*(-1-t)/4-b*(1-s)/4];Bfirst = [B1 B2 B3 B4];Jfirst = [0 1-t t-s s-1 ; t-1 0 s+1 -s-t ;s-t -s-1 0 t+1 ; 1-s s+t -t-1 0];J = [xi xj xm xp]*Jfirst*[yi ; yj ; ym ; yp]/8;B = Bfirst/J;if ID == 1D = (E/(1-NU*NU))*[1 NU 0 ; NU 1 0 ; 0 0 (1-NU)/2];elseif ID == 2D = (E/(1+NU)/(1-2*NU))*[1-NU NU 0 ; NU 1-NU 0 ; 0 0 (1-2*NU)/2]; endstr1 = D*B*u;str2 = subs(str1, {s,t}, {0,0});stress = double(str2);。
有限元法的基本思想及计算步骤
用有限元法求解问题的计算步骤比较繁多，其中最主要的计算步骤为： 1)连续体离散化。首先，应根据连续体的形状选择最能完满地描述连续体形状的单元。常见 的单元有：杆单元，梁单元，三角形单元，矩形单元，四边形单元，曲边四边形单元，四面体单 元，六面体单元以及曲面六面体单元等等。其次，进行单元划分，单元划分完毕后，要将全部单 元和结点按一定顺序编号，每个单元所受的荷载均按静力等效原理移植到结点上，并在位移受约 束的结点上根据实际情况设置约束条件。 2)单元分析。所谓单元分析，就是建立各个单元的结点位移和结点力之间的关系式。现以三 角形单元为例说明单元分析的过程。如图1所示，三角形有三个结点i，j，m。在平面问题中每个 结点有两个位移分量u，v和两个结点力分量Fx，Fy。三个结点共六个结点位移分量可用列阵(δ)e 表示： ,δ-e=*ui vi uj vj um vm+T 同样，可把作用于结点处的六个结点力用列阵{F}e表示： {F}e=[Fix Fiy Fjx Fjy Fmx Fmy]T 应用弹性力学理论和虚功原理可得出结点位移与结点力之间的关系 ,F-e=*k+e,δ-e (1)式中 [k]e——单元刚度矩阵。
有限元语言及编译器finiteelementlanguagecompiler以下简称felac是中国科学院数学与系统科是具有国际独创性的有限元计算软件是pfepg系列软件三十年成果1983年2013年的总结与提升有限元语言语法比pfepg更加简练更加灵活功能更加强大
有限元法的基本思想及计算步骤
有限元法是把要分析的连续体假想地分割成有限个单元所组成的组合体，简称离散 化。这些单元仅在顶角处相互联接，称这些联接点为结点。离散化的组合体与真实弹性 体的区别在于：组合体中单元与单元之间的联接除了结点之外再无任何关联。但是这种 联接要满足变形协调条件，即不能出现裂缝，也不允许发生重叠。显然，单元之间只能 通过结点来传递内力。通过结点来传递的内力称为结点力，作用在结点上的荷载称为结 点荷载。当连续体受到外力作用发生变形时，组成它的各个单元也将发生变形，因而各 个结点要产生不同程度的位移，这种位移称为结点位移。在有限元中，常以结点位移作 为基本未知量。并对每个单元根据分块近似的思想，假设一个简单的函数近似地表示单 元内位移的分布规律，再利用力学理论中的变分原理或其他方法，建立结点力与位移之 间的力学特性关系，得到一组以结点位移为未知量的代数方程，从而求解结点的位移分 量。然后利用插值函数确定单元集合体上的场函数。显然，如果单元满足问题的收敛性 要求，那么随着缩小单元的尺寸，增加求解区域内单元的数目，解的近似程度将不断改 进，近似解最终将收敛于精确解。
有限元第五章 有限元动力学基本原理
第五章 有限元动力学分析基本原理
在前面的介绍中，我们均假设作用在弹性体（或结 构）上的载荷与时间无关，与此相应的，位移、应力 及应变等也都和时间无关，即前面介绍的全部内容皆 称结构静力学有限元方法。但工程实际中还存在着另 外一类载荷与时间有关的动载荷作用于结构或弹性体， 此时，相应的位移、应力、应变等都与时间有关，而 且必须考虑惯性力和加速度等因素，这类分析或问题， 成为动力学分析。 对于质点—弹簧系统的振动，大家比较熟悉，例如 一个自由度为n的质点—弹簧振系，其动平衡方程为
停止迭代 此时为低阶特性
2
1
（ i 1）
(i 1)
三、机械结构固有频率与振型
2.矩阵迭代法
例题：已知一振动系统的质量矩阵、刚度矩阵用迭 代法计算其最高阶固有频率和振型。
1 0 0 3 2 0 M 0 2 0 K 2 5 3 0 0 3 0 3 3 1 1 1 解： 1 1 1.5 1.5 K 1 1.5 11 / 6
& & & M C K P
第五章 有限元动力学分析基本原理
上式中每一项的含义不同
& & M C 为阻尼力
K 为弹性力
对于单元体而言，可以得到类似的上述方程
e T N N dV V
于是，令e T V来自m N N dV
一、单元质量矩阵的计算
1.一致质量矩阵
e
m 的计算式是通式，并因为计算质量矩阵和刚度矩
阵使用的形状函数一致，因此被称为一致质量阵。
有限单元法知识点总结
有限单元法知识点总结1. 有限元法概述有限单元法（Finite Element Method ，简称FEM）是一种数值分析方法，适用于求解工程结构、热传导、流体力学等领域中的强耦合、非线性、三维等问题，是一种求解偏微分方程的数值方法。
有限元法将连续的物理问题抽象为由有限数量的简单几何单元（例如三角形、四边形、四面体、六面体等）组成的离散模型，通过对单元进行适当的数学处理，得到整体问题的近似解。
有限元法广泛应用于工程、材料、地球科学等领域。
2. 有限元法基本原理有限元法的基本原理包括离散化、加权残差法和形函数法。
离散化是将连续问题离散化为由有限数量的简单单元组成的问题，建立有限元模型。
加权残差法是选取适当的残差形式，并通过对残差进行加权平均，得到弱形式。
形函数法是利用一组适当的形函数来表示单元内部的位移场，通过形函数的线性组合来逼近整体位移场。
3. 有限元法的步骤有限元法的求解步骤包括建立有限元模型、建立刚度矩阵和载荷向量、施加边界条件、求解代数方程组和后处理结果。
建立有限元模型是将连续问题离散化为由简单单元组成的问题，并确定单元的连接关系。
建立刚度矩阵和载荷向量是通过单元的应变能量和内力作用，得到整体刚度矩阵和载荷向量。
施加边界条件是通过给定位移或力的边界条件，限制未知自由度的取值范围。
求解代数方程组是将有限元模型的刚度方程和载荷方程组成一个大型代数方程组，通过数值方法求解。
后处理结果是对数值结果进行处理和分析，得到工程应用的有用信息。
4. 有限元法的元素类型有限元法的元素类型包括结构单元、板壳单元、梁单元、壳单元、体单元等。
结构单元包括一维梁单元、二维三角形、四边形单元、三维四面体、六面体单元。
板壳单元包括各种压力单元、弹性单元、混合单元等。
梁单元包括梁单元、横梁单元、大变形梁单元等。
壳单元包括薄壳单元、厚壳单元、折叠单元等。
体单元包括六面体单元、锥体单元、八面体单元等。
5. 有限元法的数学基础有限元法的数学基础包括变分法、能量方法、有限元插值等。
第五章 有限元法
e
单元的节点力和节点位移的关系可写成：
F1 K11 K12 U1 F2 K 21 K 22 U 2
1 1 1
j 其中 Fi 表示 j单元 i节点的节点力， Fi
j
f xij j f yi
1 y1
节点2的x方向
Rx 2 f x12 f x22 EA cos2 u1 cos sin v1 cos2 u2 cos sin v2 l1
节点2的y方向
Ry 2 f y12 f y22 EA cos sin u1 sin 2 v1 cos sin u2 sin 2 v2 EA v2 v3 l1 l2
混合法：未知量为力和位移
以推导方法分类
直接法
变分法 加权残值法
§5.1 有限元的直观方法和基本概念
y
由两根杆件组成的桁架，杆 件的截面积都为A，弹性模量为 E，长度为l1及l2，在节点处受有
o x
1
R y2 (v 2)
2
Rx2 (u 2)
2
外力Rx1，Ry1，Rx2，Ry2，Rx3，
写成矩阵形式：
f x11 cos 2 1 f y1 EA cos sin 1 f x 2 l1 cos 2 f y12 cos sin cos sin sin 2 cos sin sin 2 cos 2 cos sin cos 2 cos sin
2=
y 2
1 2
k 41 k31
u =v =0
1 2
cos θ
有限元三角形单元和四边形单元
有限元三角形单元和四边形单元
有限元分析是工程应用最广泛的方法之一，它可以帮助我们了解结构内部力学行为特征。
对于有限元分析来说，有两种典型的单元：三角形单元和四边形单元。
三角形单元由三个顶点组成，四边形单元由四个顶点组成。
其基本概念是根据有限元分析的原则，将被研究的区域分解成若干个小的连续单元，单元由每个顶点表示，并且形成多边形，以此来模拟物体总体的行为。
三角形单元的特性是它的每个内角都能满足三角函数，可以极大的提高计算质量，避免出现趋势不准确的情况。
在较大范围内，一致性面积越小，他们之间的拉伸应力也就越小，更有利于精确计算。
而四边形单元则更加适合于细粒度的物体，在对细粒度物体进行研究时，可以将其细分成多个正方形小块，从而简化计算难度，提高计算效率。
在有限元分析中，三角形单元和四边形单元可以因应不同的需求而采用，只要能充分构建出更准确的结构行为模型，增加更多的灵活性和应用场景就可以得到更精准的计算结果。
第五章轴对称问题有限元法
轴对称结构的几何模型是一个表示子午面 形状的平面图形，与平面问题相比，轴对 称问题的应力与应变分量各多一个。
第二节 轴对称问题有限元法
一、结构离散 本身是一个三维结构，由于形状和载荷 的特殊性，其网格划分仅在任一子午面上 进行，为平面网格。 几何模型：一个表示子午面形状的平面图 形，用相应的轴对称实体单元划分。 三节点三角形环单元
二、轴对称问题的应力应变特点
轴对称问题的特点是结构的位移、应变和应力都呈轴对称分布。 通常采用柱坐标系（ r、、z ），并以 z 轴为对称轴。结构 中的位移、应变和应力与角度 无关。 可以取出结构的任一子午面进行分析，从而将三维问题转化为二维问题
z z
rz rz
四、总刚集成
求出每个三角形环单元的刚度矩阵后，即可按照第二章介 绍的总体刚度矩阵的集成方法，得出结构的总刚矩阵。
K {q} {R}
五、等效节点载荷的计算
计算轴对称问题的等效节点载荷与平面问题有所不同，因为轴对称结构 的子午面上的一个节点是一个关于对称轴中心对称的圆环，故当计算集 中力、表面力和体积力时，应在整个环上积分。
T
T
rc v 0 N i
0 Nj
0 N m drdz
Arc v 0 1 0 1 0 1T 3
（2）惯性离心力 设结构的转速为n，角速度为ω，其材料密度为ρ，则单元 体积的离心力为ρω2r。 单元的单位体积力为
prv 2 rc Pv pzv 0
e T
k2 A1cr bs f s A2br cs k3 A1cs br f r A2bs cr k4 cr cs A2br bs
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
5、单元内任意点的应力列阵
{ } [ x y xy ]T
y
i m
（1-5）
·
j
6、几何方程
x
u v u v x , y , xy x y y x
将上式代入式(1-4），
u { } x y u y x
图1-2
i ]T
(i, j, m)
（1-10）
一个三角形单元有3个节点（以 i、j、m为 序）， 共有6个节点位移分量。其单元位移或单元节点位移 列阵为：
i T { } j [ui i u j j um m ] m (1-11）
（7）位移函数的形式 一般选为完全多项式。为实现（4）—（6）的要 求，根据Pascal三角形由低阶到高阶按顺序、对称地 选取；多项式的项数等于（或稍大于）单元节点自由 度数。
1 x y
2 3 2
x xy y
2 2 3
x x y xy y
x 4 x 3 y x 2 y 2 xy 3 y 4
（1-7）
x
E ( x y ) 2 1
式中 E、——弹性模量、泊松比。
上式可简写为
{ } [ D]{ }
其中
（1-8）
1 E [ D] 2 1 0
对 1 称 1 0 2
(1-9）
矩阵[D]称为弹性矩阵。对于平面应变问题，将式(1-9） 中的E换为
1 u [( ai bi x ci y)ui (a j b j x c j y)u j (am bm x cm y)u m ] 2A (1-16） 1 [( ai bi x ci y) i (a j b j x c j y) j (am bm x cm y) m ] 2A
①、②、③、④单元的位移函数都是
u a1 a2 x a3 y
a4 a5 x a6 y
可以看出： 位移函数在单元内是连续的； 位移函数在单元之间的边界上也连续吗？是。 以③、④的边界26为例
5 6 5 6
③ ④
2y x 3 2
u6 u
6
③
u2 u2
u6
u
④
2 3
两条直线上有两个点重合，此两条直线必全重合。
1 E ， 换为 2
1
。
{ } [ D]{ }
(1-8）
各种类型结构的弹性物理方程都可用式(1-8）描 述。但结构类型不同，力学性态 (应力分量、应变分 量)有区别， 弹性矩阵[D]的体积和元素是不同的。
5.3 位移函数和形函数
• 1、位移函数概念 由于有限元法采用能量原理进行单元分析，因而 必须事先设定位移函数。 “位移函数”也称 “位移 模式”，是单元内部位移变化的数学表达式，设为坐 标的函数。 一般而论，位移函数选取会影响甚至严重影响 计算结果的精度。在弹性力学中，恰当选取位移函数 不是一件容易的事情；但在有限元中，当单元划分得 足够小时，把位移函数设定为简单的多项式就可以获 得相当好的精确度。这正是有限单元法具有的重要优 势之一。
式中
ai x j ym xm y j
bi y j ym
j
ci x j xm
m (i, j, m) (1-17）
i
式(1-17）中（i、j、m）意指：按i、j、m依次轮换下 标，可得到aj、bj、cj～am、bm、cm。后面出现类似情 况时，照此推理。式(1-17）表明： aj、bj、cj～am、 bm、cm是单元三个节点坐标的函数。
1 u [( ai bi x ci y)ui (a j b j x c j y)u j (am bm x cm y)u m ] 2A (1-16） 1 [( ai bi x ci y) i (a j b j x c j y) j (am bm x cm y) m ] 2A
Ni ( xi , yi ) 1 Ni ( x j , y j ) 0 Ni ( xm , ym ) 0
（i、j、m）
利用 N i 1 (ai bi x ci y )和ai、bi、ci公式证明 2A
y vi
i
vm
m
um
v u j
vj uj x
·
ui
本问题选位移函数（单元中任意一点的位移与节点 位移的关系）为简单多项式：
u a1 a2 x a3 y
a4 a5 x a6 y (1-12）
式中：a1、a2、…、a6——待定常数，由单元位移的 6个分量确定。a1、a4代表刚体位移，a2、 a3 、 a5 、 a6 代表单元中的常应变，而且，位移函数是 连续函数。 u v
第五章 三角形单元的有限元法
5.1 基本思想
把整体结构离散为有限个单元，研究单元的平衡和变 形协调；再把这有限个离散单元集合还原成结构，研究离 散结构的平衡和变形协调。划分的单元大小和数目根据计 算精度和计算机能力来确定。
P
2
○
4
6
8
10
4
4
6 6
6
8
② ①
○
④ ③
3 5
⑥ ⑤
7
⑧ ⑦
9
④
⑥ ⑤
① ②
（4）位移函数中必须包含单元的刚体位移。 （5）位移函数中必须包含单元的常应变。 （6）位移函数在单元内要连续。相邻单元间要尽 量协调。 条件（4）、（5）构成单元的完备性准则。 条件（6）是单元的位移协调性条件。 理论和实践都已证明，完备性准则是有限元解收 敛于真实解的必要条件。单元的位移协调条件构成有 限元解收敛于真实解的充分条件。 容易证明，三角形三节点常应变单元满足以上必 要与充分条件。
用形函数把式(1-16）写成矩阵，有
u N i v 0
缩写为
0 Ni
Nj 0
0 Nj
Nm 0
ui v i 0 u j Nm v j um vm
{ f } [ N ]{ }
(1-20）
[N]为形函数矩阵，写成分块形式：
3
③
5 5 5
⑦ ⑧
7 7
1
单元、节点需编号 弹性悬臂板剖分与集合
5.2 基本力学量矩阵表示
1、单元表面或边界上任意点的表面力列阵qs
qsx {qs } [qsx qsy ]T qsy q
y
j
（1-1）
y
j
·
s
qV
i m m
·
x
i
x
2、单元内任意点的体积力列阵qV
y
m（7） j （1） i （2）
1 A 1 xj 2 1 xm
1
xi
yi yj ym
(1-15）
x
特别指出：为使求得面积的值为正值，本单元节点号 的次序必须是逆时针转向，如图所示。至于将哪个节 点作为起始节点i，则没有关系。
1 u [( ai bi x ci y)ui (a j b j x c j y)u j (am bm x cm y)u m ] 2A 同理 1 [( ai bi x ci y) i (a j b j x c j y) j (am bm x cm y) m ] 2A
u j a1 a2 x j a3 y j
um a1 a2 xm a3 ym
ui a1 a2 xi a3 yi
i a4 a5 xi a6 yi j a4 a5 x j a6 y j m a4 a5 xm a6 ym
(1-13）
从式(1-13）左边3个方程中解出待定系数a1、a2、a3为
1 a1 uj 2A um ui xi xj xm yi yj ym
1 a2 1 uj 2A 1 um
1 ui
yi yj ym
(1-14）
1 a3 1 xj 2A 1 xm
1积，有
x
x a2 , y y a5 , xy a3 a5 a2 a6
选取位移函数应考虑的问题
（1）位移函数的个数 等于单元中任意一点的位移分量个数。本单元中 有u和v，与此相应，有2个位移函数；
（2）位移函数是坐标的函数 本单元的坐标系为：x、y；
（3）位移函数中待定常数个数 待定常数个数应等于单元节点自由度总数，以 便用单元节点位移确定位移函数中的待定常数。本 单元有6个节点自由度，两个位移函数中共包含6个 待定常数。
3、形函数
（1）形函数确定 形函数是用单元节点位移分量来描述位移函数的 u a a x a y a a x a y (1 12) 插值函数。
1 2 3 4 5 6
现在，通过单元节点位移确定位移函数中的待定 常数a1、a2、…、a6 。设节点i、j、m的坐标分别为 （xi、yi）、（ xj、yj ）、（ xm、ym ），节点位移分别 为（ui、vi）、 （uj、vj） 、 （um、vm）。将它们 代入式（1-12），有
[ N ] [[Ni ] [ N j ] [ Nm ]]
(1-21）
其中子矩阵
Ni 0 [ Ni ] Ni [ I ] 0 Ni [I]是2×2的单位矩阵。
(i, j, m) (1-22）
（2）形函数性质 形函数是有限单元法中的一个重要函数，它具 有以下性质： 性质1 形函数Ni在节点i上的值等于1，在其它节点 上的值等于0。对于本单元，有
qVx T {qV } [qVx qVy ] qVy
（1-2）
3、单元内任意点的位移列阵f