8 第四章 用常应变三角形单元解力学平面问题 (2)
三角形常应变单元程序的编制与使用11页
三角形常应变单元程序的编制与使用有限元法是求解微分方程边值问题的一种通用数值方法，该方法是一种基于变分法（或变分里兹法）而发展起来的求解微分方程的数值计算方法，以计算机为手段，采用分片近似，进而逼近整体的研究思想求解物理问题。
有限元分析的基本步骤可归纳为三大步：结构离散、单元分析和整体分析。
对于平面问题，结构离散常用的网格形状有三角形、矩形、任意四边形，以三个顶点为节点的三角形单元是最简单的平面单元，它较矩形或四边形对曲边边界有更好的适应性，而矩形或四边形单元较三节点三角形有更高的计算精度。
Matlab语言是进行矩阵运算的强大工具，因此，用Matlab语言编写有限元中平面问题的程序有优越性。
本章将详细介绍如何利用Matlab语言编制三角形常应变单元的计算程序，程序流程图见图1。
有限元法中三节点三角形分析结构的步骤如下：1）整理原始数据，如材料性质、荷载条件、约束条件等，离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。
2）单元分析，建立单元刚度矩阵。
3）整体分析，建立总刚矩阵。
4）建立整体结构的等效节点荷载和总荷载矩阵5）边界条件处理。
6）解方程，求出节点位移。
7）求出各单元的单元应力。
8）计算结果整理。
计算结果整理包括位移和应力两个方面；位移计算结果一般不需要特别的处理，利用计算出的节点位移分量，就可画出结构任意方向的位移云图；而应力解的误差表现在单元内部不满足平衡方程，单元与单元边界处应力一般不连续，在边界上应力解一般与力的边界条件不相符合。
图1 程序流程图1.1 程序说明% 三角形常应变单元求解结构主程序●功能：运用有限元法中三角形常应变单元解平面问题的计算主程序。
●基本思想：单元结点按右手法则顺序编号。
●荷载类型：可计算结点荷载。
●说明：主程序的作用是通过赋值语句、读取和写入文件、函数调用等完成算法的全过程，即实现程序流程图的程序表达。
1 程序准备format short e %设定输出类型clear all %清除所有已定义变量clc %清屏●说明：format short e －设定计算过程中显示在屏幕上的数字类型为短格式、科学计数法；clear all －清除所有已定义变量，目的是在本程序的运行过程中，不会发生变量名相同等可能使计算出错的情况；clc －清屏，使屏幕在本程序运行开始时2 全局变量定义global NNODE NPION NELEM NVFIX NFORCE COORD LNODS YOUNG POISS THICKglobal FORCE FIXEDglobal BMATX DMATX SMATX AREAglobal ASTIF ASLOD ASDISPglobal FP1●说明：NNODE—单元结点数，NPION—总结点数，NELEM—单元数，NVFIX—受约束边界点数，NFORCE—结点力数，COORD—结构结点坐标数组，LNODS —单元定义数组，YOUNG—弹性模量，POISS—泊松比，THICK—厚度FORCE —节点力数组(n,3) n:受力节点数目,(n,1):作用点,(n,2):x方向,(n,3):y 方向；FIXED—约束信息数组(n,3) n:受约束节点数目, (n,1):约束点(n,2)与(n,3)分别为约束点x方向和y方向的约束情况,受约束为1否则为0BMATX—单元应变矩阵(3*6)，DMATX—单元弹性矩阵(3*3)，SMATX—单元应力矩阵(3*6)，AREA—单元面积ASTIF—总体刚度矩阵，ASLOD—总体荷载向量，ASDISP—结点位移向量FP1—数据文件指针3 打开文件FP1=fopen('input.txt','rt'); %打开输入数据文件存放初始数据●说明：FP1=fopen('input.txt','rt'); －打开已存在的输入数据文件input.txt，且设置其为只读格式，使程序在执行过程中不能改变输入文件中的数值，并用文件句柄FP1来执行FP2=fopen('output.txt','wt'); －打开输出数据文件，该文件不存在时，通过此命令创建新文件，该文件存在时则将原有内容全部删除。
有限元经典PPT第4章
Pii Kiiui
Ki1u1 Ki2u2 Kiiui K u i,i1 i1
ui
n
Kiiui Kiiui
Kiju j
4.1.2 平面应力问题有限元的基本思想和瑞雷－里兹法
v3 f3y
3
u3
f3x
f1y v1 u1
1 f1x
v2 f2y u2
2 f2x
给定一个三角形单元和作用在角点上 的六个力，要求得六个角点的位移。 或者是要求三角形角点发生指定的位 移，在三角形三个角点如何加力？
很显然，问题的精确解很困难。采用 瑞雷－里兹法求近似式解
e号单元的三个节点I,j,k的力对应的 力的平衡方程是第2i-1,2i;2j-1,2j;2k1,2k个平衡方程
e号单元的三个节点I,j,k的位移是第 2i-1,2i;2j-1,2j;2k-1,2k个未知数
弹性模量：E 横截面积：A
1
1 L
2
2L
3
局部系单元刚度阵：
k
1
EA L
1 -1
-1
1
2 集成总刚：
0 1
解得：
ux uy
L EA
3.8284L
EA
i
j
第一类位移条件：
Ki1u1 Ki2u2 Kiiui Ki1ui1
ui 0
令： Kij 0 i j
m
vi 0
Kii 1
um 0
Pi 0
ui 0
第二类位移条件：um um
大数
充大数法： Kii Kii
第一步：求转换矩阵
k2
EA 1 2L -1
-1
1
P
cos 0
T sin
常应变三角形单元
四.应力结果的整理
位移的计算结果一般比应力、内力结果精度高。位移达到满意结果， 由几何方程求应变，再由物理方程求应力，结果的精度较差。上述三角形 单元为常应力，矩形单元应力线性变化，而工程问题的应力是比较复杂的。 为更好地反应实际应力情况，需要对计算结果进行整理。常用处理方法有 两种:绕结点平均法和两单元平均法。
四.应力结果的整理
位移的计算结果一般比应力、内力结果精度高。位移达到满意结果， 由几何方程求应变，再由物理方程求应力，结果的精度较差。上述三角形 单元为常应力，矩形单元应力线性变化，而工程问题的应力是比较复杂的。 为更好地反应实际应力情况，需要对计算结果进行整理。常用处理方法有 两种:绕结点平均法和两单元平均法。
宜
不宜
5.相邻单元的尺寸尽可能接近。 6.结点所连接的单元个数尽可能一致。 宜 二.结点编码 尽可能使相关结点的结点编码差值最小. 1 2 3 4 5 6 7 1 3 5 7 9 11 13 不宜
8
9
10
11 12
13
14
2
4
6
8
10
12
14
总刚半带宽=(相关结点最大差值+1)*结点位移数 总刚半带宽=(7+1)*2=16 总刚半带宽=(2+1)*2=6 总刚需占用的存贮空间为: 6 *14*2=168
宜
不宜
5.相邻单元的尺寸尽可能接近。 6.结点所连接的单元个数尽可能一致。 宜 不宜
§3.3 有限元分析应注意的问题和结果整理
一.结点的选择和单元划分 1.集中力作用点、分布力突变点、支承点应选作结点。 2.不同厚度、不同材料的部分不应划在同一个单元。 3.应力变化大处单元应密集一些。结点的多少与疏密要考虑计算 机的容量和计算精度。 4.单元边界的边长之比应尽可能靠近1。
初应变
4.4
平面问题程序(一)
起点号，终点号，生成 的点数，相邻点号差值， “相邻两点间距”。 2-3）读入结点荷载值 有荷载的结点数 结点号，X方向荷载值， Y方向荷载值。 2-4）单元的整体结点码 如果规则无孔：不需要 输入 如果规则有孔： 待修改
2-2) 结点坐标 如果不规则：按结点号 顺序读入全部结点的坐 标值。 如果规则无孔：孔标志, X方向单元数，Y方向单 元数，X方向单元长度， Y方向单元长度。 如果规则有孔： 控制结点数，生成结点 类数。 结点号，X，Y
e
第一项体积力、第二项结 点力、第三项表面力的外 F ( t N dl)]d 力势。
T e T l ij e
4.2
TT
常应变三角形单元
T ee TT T
5) 令总势能一阶变分等于零，推导单元刚度方程
t S B F e N F dA tB t tSS Bdd Fd t N Fb bdA b
4.2.3 单元列式 1) 微分算子矩阵 x 0 A 0 y
T
常应变三角形单元
y x
平面应力问题 B B B
2) 应变、应力矩阵
A d Bd e
D S d e
4.3
矩形双线性单元
4.3
矩形双线性单元
4 1 3 [d]=[u，v]T
3) 位移模式 u=Niui ； v=Nivi。 或以矩阵表示为
2
2
N N1
Ni N i 0
N4 0 d N d e Ni
2
单元结点 位移矩阵
4.2
常应变三角形单元
有限元分析第四章
19
4）形函数的性质
形函数是有限单元法中的一个重要函数，它具 有以下性质： 性质1 形函数Ni在节点i上的值等于1，在其它节点 上的值等于0。对于本单元，有
20
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
对于一个具体问题进行分析，不管采用什么样的单元， 分析过程与思路是一样的，所不同的只是各种单元的位移模 式和单元刚度矩阵不一样，其他的包括整体刚度矩阵的组装 过程都完全一样，所以我们仅仅对矩形单元位移模式的求取 和单元刚度矩阵的求解加以介绍。
4.7 收敛准则
可以证明，对于一个给定的位移模式，其刚度系统的数 值要比精确值大。所以，在给定载荷的作用下，有限元计算 模型的变形要比实际结构的变形小。因而，当单元网格分得 越来越细时，位移的近似解将由下方收敛于精确解，即得到 真实解的下界。 为了保证解答的收敛性，要求选取的位移模式必须满足 以下三个条件： 1）位移模式必须包含单元的刚体位移 也就是说，当节点位移是某个刚体位移所引起时，弹 性体内将不会产生应变。所以位移模式不但要具有描述单元 本身形变的能力，而且还要具有描述由其他变形而通过节点 位移引起单元刚体位移的能力。例如，三角形三节点位移模 式中，常数项就是用于提供刚体位移的。
Ni(x、y)
1 i(xi，yi) x xi
x xi N i ( x, y ) 1 x j xi
N m ( x, y ) 0
证
N
y j (xj，yj)
m (xm，ym)
xj
x
N i ( x, y )
有限元方法课件 第四章 平面三角形单元
第四章 平面三角形单元
§4–1 有限元法的基本思想 §4–2 三角形常应变单元 §4–3 形函数的性质 §4–4 刚度矩阵 §4–5 等效节点力载荷列阵 §4–6 有限元分析的实施步骤 §4–7 计算实例
§4-1 有限元法的基本思想
一、有限元法的基本思想 假想的把一连续体分割成数目有限的小体（单元），
vi (Vi )
i ui (Ui )
m
um (Um )
o
x
图4-2 平面三角形单元
将 (d) 式代入 (b) 式的第一式，经整理后得到
u 1 2ai源自bi x ci yuiaj
bjx cj y
uj
am bm x cm yum
(e)
其中 同理可得
ai
xj xm
yj ym
x j ym xm y j
这样，位移模式 (e) 和 (f) 就可以写为
u Ni ui N j u j N mum v Nivi N jv j Nmvm
(4-11)
也可写成矩阵形式
f
u v
Ni I
NjI
NmI e N e
(4-12)
式中 I是二阶单位矩阵；Ni 、Nj 、Nm 是坐标的函数， 它们反映了单元的位移状态，所以一般称之为形状函数，简 称形函数。矩阵 [N] 叫做形函数矩阵。三节点三角形单元的 形函数是坐标的线性函数。单元中任一条直线发生位移后仍 为一条直线，即只要两单元在公共节点处保持位移相等。则 公共边线变形后仍为密合。
f N e
（4-1）
f ——单元内任一点的位移列阵； e——单元的结点位移列阵；
N ——单元的形函数矩阵；（它的元素是任一点位置坐
平面问题三角形单元有限元课件
(i, j, m)
(1-26）
由于 A, bi , ci , b j , c j , bm , cm 与x、y无关，都是常量，因此 [B]矩阵也是常量。单元中任一点的应变分量是[B]矩阵
与单元位移的乘积，因而也都是常量。因此，这种单元
被称为常应变单元。
2、单元应力
{} [B]{ }
j
bj
x
c
j
y)
j
(am
bm
x
cm
y)
m
]
1
u
2 A [(ai
bi
x
ci
y)ui
(a
j
bj
x
c
j
y)u
j
(am
bm x
cm y)um ]
(1-16）
1 2A
[(ai
bi
x
ci
y)
i
(a
j
b
j
x
c
j
y)
j
(am
bm
x
cm
y)
m
]
j
式中
ai x j ym xm y j
mi
bi y j ym
(i, j, m) (1-17）
bi
x
ci
y)ui
(a
j
bj
x
c
j
y)u
j
(am
bm
x
cm y)um ]
(1-16）
1 2A
[(ai
bi
x
ci
y)
i
(a
j
bj
x
c
j
y)
j
(am
三角形常应变单元程序的编制与使用共18页文档
三角形常应变单元程序的编制与使用有限元法是求解微分方程边值问题的一种通用数值方法，该方法是一种基于变分法（或变分里兹法）而发展起来的求解微分方程的数值计算方法，以计算机为手段，采用分片近似，进而逼近整体的研究思想求解物理问题。
有限元分析的基本步骤可归纳为三大步：结构离散、单元分析和整体分析。
对于平面问题，结构离散常用的网格形状有三角形、矩形、任意四边形，以三个顶点为节点的三角形单元是最简单的平面单元，它较矩形或四边形对曲边边界有更好Array的适应性，而矩形或四边形单元较三节点三角形有更高的计算精度。
Matlab语言是进行矩阵运算的强大工具，因此，用Matlab语言编写有限元中平面问题的程序有优越性。
本章将详细介绍如何利用Matlab语言编制三角形常应变单元的计算程序，程序流程图见图1。
有限元法中三节点三角形分析结构的步骤如下：1）整理原始数据，如材料性质、荷载条件、约束条件等，离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。
2）单元分析，建立单元刚度矩阵。
3）整体分析，建立总刚矩阵。
4）建立整体结构的等效节点荷载和总荷载矩阵5）边界条件处理。
6）解方程，求出节点位移。
7）求出各单元的单元应力。
8）计算结果整理。
计算结果整理包括位移和应力两个方面；位移计算结果一般不需要特别的处理，利用计算出的节点位移分量，就可画出结构任意方向的位移云图；而应力解的误差表现在单元内部不满足平衡方程，图1 程序流程图单元与单元边界处应力一般不连续，在边界上应力解一般与力的边界条件不相符合。
1.1 程序说明% 三角形常应变单元求解结构主程序●功能：运用有限元法中三角形常应变单元解平面问题的计算主程序。
●基本思想：单元结点按右手法则顺序编号。
●荷载类型：可计算结点荷载。
●说明：主程序的作用是通过赋值语句、读取和写入文件、函数调用等完成算法的全过程，即实现程序流程图的程序表达。
1 程序准备format short e %设定输出类型clear all %清除所有已定义变量clc %清屏●说明：format short e －设定计算过程中显示在屏幕上的数字类型为短格式、科学计数法；clear all －清除所有已定义变量，目的是在本程序的运行过程中，不会发生变量名相同等可能使计算出错的情况；clc －清屏，使屏幕在本程序运行开始时2 全局变量定义global NNODE NPION NELEM NVFIX NFORCE COORD LNODS YOUNG POISS THICKglobal FORCE FIXEDglobal BMATX DMATX SMATX AREAglobal ASTIF ASLOD ASDISPglobal FP1说明：NNODE—单元结点数，NPION—总结点数， NELEM—单元数，NVFIX—受约束边界点数，NFORCE—结点力数，COORD—结构结点坐标数组，LNODS —单元定义数组，YOUNG—弹性模量，POISS—泊松比，THICK—厚度FORCE —节点力数组(n,3) n:受力节点数目,(n,1):作用点,(n,2):x方向,(n,3):y方向； FIXED—约束信息数组(n,3) n:受约束节点数目, (n,1):约束点 (n,2)与(n,3)分别为约束点x方向和y方向的约束情况,受约束为1否则为0BMATX—单元应变矩阵(3*6)， DMATX—单元弹性矩阵(3*3)，SMATX—单元应力矩阵(3*6)，AREA—单元面积ASTIF—总体刚度矩阵，ASLOD—总体荷载向量，ASDISP—结点位移向量FP1—数据文件指针3 打开文件FP1=fopen('input.txt','rt'); %打开输入数据文件存放初始数据●说明：FP1=fopen('input.txt','rt'); －打开已存在的输入数据文件input.txt，且设置其为只读格式，使程序在执行过程中不能改变输入文件中的数值，并用文件句柄FP1来执行FP2=fopen('output.txt','wt'); －打开输出数据文件，该文件不存在时，通过此命令创建新文件，该文件存在时则将原有内容全部删除。
用常应变三角形单元解弹性力学平面问题的程序
用常应变三角形单元解弹性力学平面问题的程序******************************************************************** ANALYSIS PROGTAM OF FINITE ELEMENT METHOD ** FOR PLANE STRESS/STRAIN OF TRIANGULAR ELEMENT ** ----- FEMT3.FOR ----- **------------------------------------------------------------- ** Subroutines: 1-SDATA, 2-STE, 3-ATE, 4-DTE, 5-BTE, 6-STIFF ** 7-EQUPE, 8-INSCD, 9-BGSMT, 10-SIGME ********************************************************************DIMENSION LND(50,3),X(100),Y(100),JR(20,3),PJ(20,3),P(200)REAL KS(200,100)OPEN(5,FILE='FEMT3.DAT')OPEN(6,FILE='FEMT3.OUT',STATUS='NEW')READ(5,*) NJ,NE,NS,NPJ,IPS(结点、单元、支承、荷载、类型)WRITE(6,*)' FINITE ELEMENT ANALYSIS IN PLANE PROBLEM'WRITE(6,*)' SOURCE DATA OUTPUT'WRITE(6,20) NJ,NE,NS,NPJ,IPS20 FORMAT(4X,'NJ',3X,'NE',3X,'NS',3X,'NPJ',2X,'IPS'/1X,5I5)IF(IPS.EQ.0) WRITE(6,*)' PLANE STRESS PROBLEM'IF(IPS.EQ.1) WRITE(6,*)' PLANE STRAIN PROBLEM'CALL SDATA(NJ,NE,NS,NW,NPJ,IPS,E,PR,T,V,LND,X,Y,JR,PJ)NJ2=2*NJWRITE(6,50) NJ250 FORMAT(/1X,'DEGREES OF FREEDOM=',I5)WRITE(6,60) NW60 FORMAT(1X,'BAND WIDTH=',I5)CALL STIFF(NJ,NE,NJ2,NW,LND,X,Y,E,PR,T,KS)(总刚6)CALL EQUPE(NJ,NE,NPJ,NJ2,T,V,LND,X,Y,PJ,P)({P}7)CALL INSCD(NS,NW,NJ2,JR,KS,P)(引入支承条件8)CALL BGSMT(NJ,NJ2,NW,KS,P)（解方程9）CALL SIGME(NE,NJ,NJ2,E,PR,LND,X,Y,P)（求应力10）CLOSE(5)CLOSE(6)END*--------------------------------------------------------C SUBPROGRAM-1C INPUT STRUCTURAL DATASUBROUTINE SDATA(NJ,NE,NS,NW,NPJ,IPS,E,PR,* T,V,LND,X,Y,JR,PJ)DIMENSION LND(NE,3),X(NJ),Y(NJ),JR(NS,3),PJ(NPJ,3)READ(5,*) E,PR,T,V(弹性模量、泊松比、厚度、容重)WRITE(6,10) E,PR,T,V10 FORMAT(/6X,'E',10X,'PR',9X,'T',9X,'V'/,4F10.2)READ(5,*)((LND(I,J),J=1,3),I=1,NE)（结点编码）WRITE(6,20)20 FORMAT(/1X,'ELEMENT INFORMATION'/3X,'ELEM',3X,* 'I J K'/)WRITE(6,30)(I,(LND(I,J),J=1,3),I=1,NE)30 FORMAT(1X,4I5)READ(5,*)(X(I),Y(I),I=1,NJ)（结点坐标）WRITE(6,40)40 FORMAT(/1X,'COORDINATES OF NODES'/3X,'NODES',* 8X,'X',13X,'Y')WRITE(6,50)(I,X(I),Y(I),I=1,NJ)50 FORMAT(1X,I5,2E15.6)READ(5,*)((JR(I,J),J=1,3),I=1,NS)（约束信息）WRITE(6,60)60 FORMAT(/1X,'CONSTRAINED NODES'/3X,'NODE',3X,'X',4X,'Y') WRITE(6,70)((JR(I,J),J=1,3),I=1,NS)70 FORMAT(1X,3I5)READ(5,*)((PJ(I,J),J=1,3),I=1,NPJ)（荷载信息）WRITE(6,80)80 FORMAT(/1X,'LOAD CASES'/3X,'NODE',8X,'X',13X,'Y')WRITE(6,90)((PJ(I,J),J=1,3),I=1,NPJ)90 FORMAT(1X,F5.0,2E15.6)100 NW=0（半带宽）DO 110 IE=1,NEDO 110 I=1,3DO 110 J=1,3IW=IABS(LND(IE,I)-LND(IE,J))IF(NW.LT.IW) THENNW=IWENDIF110 CONTINUENW=(NW+1)*2IF(IPS.NE.0) THENE=E/(1.0-PR*PR)PR=PR/(1.0-PR)ENDIFEND*---------------------------------------------------------C SUBPROGRAM-2C CALCULATE ELEMENT STIFFNESS MATRIXSUBROUTINE STE(IE,NJ,NE,LND,X,Y,E,PR,T,KE)DIMENSION LND(NE,3),X(NJ),Y(NJ),B(3,6),D(3,3)REAL KE(6,6)CALL ATE(IE,NJ,NE,LND,X,Y,AE)CALL DTE(E,PR,D)CALL BTE(IE,NJ,NE,LND,X,Y,AE,B)DO 10 I=1,6DO 10 J=1,6KE(I,J)=0.DO 10 K=1,3DO 10 K1=1,310 KE(I,J)=KE(I,J)+B(K,I)*D(K,K1)*B(K1,J)C=AE*TDO 30 I=1,6DO 30 J=1,630 KE(I,J)=KE(I,J)*CEND*------------------------------------------------ C SUBPROGRAM-3C CALCULATE ELEMENT AREASUBROUTINE ATE(IE,NJ,NE,LND,X,Y,AE)DIMENSION LND(NE,3),X(NJ),Y(NJ)I=LND(IE,1)J=LND(IE,2)K=LND(IE,3)XIJ=X(J)-X(I)YIJ=Y(J)-Y(I)XIK=X(K)-X(I)YIK=Y(K)-Y(I)AE=.5*(XIJ*YIK-XIK*YIJ)END*---------------------------------------------- C SUBPROGRAM-4C CALCULATE ELASTICITY MATRIXSUBROUTINE DTE(E,PR,D)DIMENSION D(3,3)DO 10 I=1,3DO 10 J=1,310 D(I,J)=0.D(1,1)=E/(1.-PR*PR)D(1,2)=E*PR/(1.-PR*PR)D(2,1)=D(1,2)D(2,2)=D(1,1)D(3,3)=.5*E/(1.+PR)END*------------------------------------------------ C SUBPROGRAM-5C CALCULATE MATRIX [B]SUBROUTINE BTE(IE,NJ,NE,LND,X,Y,AE,B)DIMENSION LND(NE,3),X(NJ),Y(NJ),B(3,6)I=LND(IE,1)J=LND(IE,2)K=LND(IE,3)DO 10 II=1,3DO 10 JJ=1,610 B(II,JJ)=0.B(1,1)=Y(J)-Y(K)B(1,3)=Y(K)-Y(I)B(1,5)=Y(I)-Y(J)B(2,2)=X(K)-X(J)B(2,4)=X(I)-X(K)B(2,6)=X(J)-X(I)B(3,1)=B(2,2)B(3,2)=B(1,1)B(3,3)=B(2,4)B(3,4)=B(1,3)B(3,5)=B(2,6)B(3,6)=B(1,5)DO 20 I1=1,3DO 20 J1=1,620 B(I1,J1)=.5/AE*B(I1,J1)END*------------------------------------------------------- C SUBPROGRAM-6C CALCULATE GLOBAL STIFFNESS MATRIXSUBROUTINE STIFF(NJ,NE,NJ2,NW,LND,X,Y,E,PR,T,KS) DIMENSION LND(NE,3),X(NJ),Y(NJ)REAL KS(NJ2,NW),KE(6,6)DO 5 I=1,NJ2DO 5 J=1,NW5 KS(I,J)=0.DO 10 IE=1,NECALL STE(IE,NJ,NE,LND,X,Y,E,PR,T,KE)DO 10 I=1,3IZ=LND(IE,I)DO 10 II=1,2IH =2*(I -1)+IIIDH=2*(IZ-1)+IIDO 10 J=1,3JZ=LND(IE,J)DO 10 JJ=1,2L =2*(J -1)+JJIL=2*(JZ-1)+JJIF(IL.GE.IDH) THENIDL=IL-IDH+1KS(IDH,IDL)=KS(IDH,IDL)+KE(IH,L)ENDIF10 CONTINUEEND*-------------------------------------------------------- C SUBPROGRAM-7C CALCULATE NODAL LOAD VECTORSUBROUTINE EQUPE(NJ,NE,NPJ,NJ2,T,V,LND,X,Y,PJ,P) DIMENSION LND(NE,3),X(NJ),Y(NJ),PJ(NPJ,3),P(NJ2) DO 10 I=1,NJ210 P(I)=0.DO 20 I=1,NPJII=PJ(I,1)P(2*II-1)=PJ(I,2)20 P(2*II)=PJ(I,3)30 IF(V.EQ.0.) GOTO 50DO 40 IE=1,NECALL ATE(IE,NJ,NE,LND,X,Y,AE)PE=-V*AE*T/3.DO 40 I=1,3II=LND(IE,I)40 P(2*II)=P(2*II)+PE50 RETURNEND*---------------------------------------------C SUBPROGRAM-8C INTRODUCE BOUNDARY CONDITIONSUBROUTINE INSCD(NS,NW,NJ2,JR,KS,P)DIMENSION P(NJ2),JR(NS,3)REAL KS(NJ2,NW)DO 30 I=1,NSIR=JR(I,1)DO 30 J=2,3IF(JR(I,J).EQ.0) GOTO 30II=2*IR+J-3KS(II,1)=1.DO 10 JJ=2,NW10 KS(II,JJ)=0.IF(II.GT.NW) JO=NWIF(II.LE.NW) JO=IIDO 20 JJ=2,JO20 KS(II-JJ+1,JJ)=0.P(II)=0.30 CONTINUEEND*-------------------------------------------C SUBPROGRAM-9C SOLVE EQUATIONS BY GS METHODSUBROUTINE BGSMT(NJ,NJ2,NW,KS,P)DIMENSION P(NJ2)REAL KS(NJ2,NW)NJ1=NJ2-1DO 20 K=1,NJ1IF(NJ2.GT.K+NW-1) IM=K+NW-1IF(NJ2.LE.K+NW-1) IM=NJ2K1=K+1DO 20 I=K1,IML=I-K+1C=KS(K,L)/KS(K,1)IW=NW-L+1DO 10 J=1,IWM=J+I-K10 KS(I,J)=KS(I,J)-C*KS(K,M)20 P(I)=P(I)-C*P(K)P(NJ2)=P(NJ2)/KS(NJ2,1)DO 40 I1=1,NJ1I=NJ2-I1IF(NW.GT.NJ2-I+1) JO=NJ2-I+1IF(NW.LE.NJ2-I+1) JO=NWDO 30 J=2,JOK=J+I-130 P(I)=P(I)-KS(I,J)*P(K)40 P(I)=P(I)/KS(I,1)WRITE(6,50)50 FORMAT(/1X,'NODAL DISPLACEMENTS'/3X,* 'NODE',5X,'X-DISP.',8X,'Y-DISP.')DO 60 I=1,NJ60 WRITE(6,70) I,P(2*I-1),P(2*I)70 FORMAT(1X,I5,2E15.6)END*--------------------------------------------------- C SUBPROGRAM-10C CALCULATE ELEMENT STRESS MATRIXSUBROUTINE SIGME(NE,NJ,NJ2,E,PR,LND,X,Y,P)DIMENSION LND(NE,3),X(NJ),Y(NJ),D(3,3),B(3,6), * S(3,6),ST(3),P(NJ2),DE(6)WRITE(6,*)WRITE(6,*)' ELEMENT STRESSES'CALL DTE(E,PR,D)DO 50 IE=1,NECALL ATE(IE,NJ,NE,LND,X,Y,AE)CALL BTE(IE,NJ,NE,LND,X,Y,AE,B)DO 10 I=1,3DO 10 J=1,6S(I,J)=0.DO 10 K=1,310 S(I,J)=S(I,J)+D(I,K)*B(K,J)DO 20 I=1,3DO 20 J=1,2IH=2*(I-1)+JIW=2*(LND(IE,I)-1)+J20 DE(IH)=P(IW)DO 30 I=1,3ST(I)=0.DO 30 J=1,630 ST(I)=ST(I)+S(I,J)*DE(J)SGX=ST(1)SGY=ST(2)TXY=ST(3)ASG=(SGX+SGY)*.5RSG=SQRT(.25*(SGX-SGY)**2+TXY*TXY)SGMA=ASG+RSGSGMI=ASG-RSGIF(SGY.EQ.SGMI) CETA=0.IF(SGY.NE.SGMI) CETA=90.-57.29578*ATAN* (TXY/(SGY-SGMI))50 WRITE(6,60) IE,SGX,SGY,TXY,SGMA,SGMI,CETA60 FORMAT(1X,'ELEMENT NO.=',I4/2X,'SIGX=',E10.4, * 2X,'SIGY=',E10.4,2X,'TXY =',E10.4/2X,'SGMA=', * E10.4,2X,'SGMI=',E10.4,2X,'CETA=',E10.4)END。
8 第四章 用常应变三角形单元解力学平面问题 (2)解析
um xm ym
1 um ym
1 xm um
其中
1 xi yi
2 1 x j y j
1 xm ym
(c) (d) (1)
从解析几何可知，式中的 就是三角形i、j、m的面积。
为保证求得的面积为正值，节点i、j、m的编排次序必须是逆 时针方向，如图1所示。
7. 由单元的节点位移列阵计算单元应力
解出整体结构的节点位移列阵 后，再根据单元节点的 编号找出对应于单元的位移列阵 e，将 e代入（3-3）式就
可求出各单元的应力分量值。
8. 计算结果输出
求解出整体结构的位移和应力后，可有选择 地整理输出某些关键点的位移值和应力值，特别 要输出结构的 变形图、应力图、应变图、结构仿 真变形过程动画图及整体结构的弯矩、剪力图等 等。
平面问题可用三角元，四边元等。
例如：
3. 选择单元的位移模式
结构离散化后，要用单元内结点的位移通过插值来获得 单元内各点的位移。在有限元法中，通常都是假定单元的位 移模式是多项式，一般来说，单元位移多项式的项数应与单 元的自由度数相等。它的阶数至少包含常数项和一次项。至 于高次项要选取多少项，则应视单元的类型而定。
有限元法的实质是：把有无限个自由度的连续体， 理想化为只有有限个自由度的单元集合体，使问题简化 为适合于数值解法的结构型问题。
二、经典解与有限元解的区别：
微分 经 典 解 法 —— （解析法）
数目增到∞ 大小趋于 0
建立一个描述连续体 性质的偏微分方程
有限单元 离散化 集合
总体分析解
有限元法——连续体——单元——代替原连续体
式中：
Re ke e
（3-4）
——单元刚度矩阵
ke BT DBdxdydz
弹性力学平面问题2
i
ui
j
uj
o
x
（1）
单元及其位移表示
为待定常数
ui = a1 + a2xi + a3 yi vi = a4 + a5xi + a6 yi
（2）
位移(1)在结点上有： 位移 在结点上有： 在结点上有
uj = a1 + a2xj + a3 yj vj = a4 + a5xj + a6 yj uk = a1 + a2xm + a3 ym vm = a4 + a5xm + a6 ym
Ω SσΒιβλιοθήκη 对于平面问题：{ε * }T {σ }dxdy = ∫∫ {u * }T { f }dxdy + ∫ {u * }{ f }dS ∫∫
Ω Ω Sσ
相容位移：即为满足位移边界条件的位移。 相容位移：即为满足位移边界条件的位移。
最小势能原理 在一切可能位移和形变中,真正的位移和形变使总势能取最小值； 在一切可能位移和形变中,真正的位移和形变使总势能取最小值；反 使总势能取最小值者也必是真正的位移和形变。 之,使总势能取最小值者也必是真正的位移和形变。
弹性力学平面问题的有限元法
一 弹性力学的基本方程 二 弹性力学的变分原理 三 平面问题的三角形单元 四 平面问题的四边形等参单元
一 弹性力学基本方程
1、基本物理量
位移
{u} = {u( x, y, z), v( x, y, z), w( x, y, z )}
ui (i = 1, 2,3)
张量表示： 张量表示：
, ci =
1 xm ym
1 xj 1 xm
有限元平面问题
平面应力 H =
（5）单元刚度方程
K e ⋅ δ e = Pe
讨论1：平面三节点三角形单元的节点位移和 坐标变换
由于该单元的节点位移是以整体坐标系中的X方向位移（ui）和Y 方向位移（vi）来定义的，所以没有坐标变换的问题。
讨论2：平面三节点三角形单元的应变矩阵和应力矩 阵为常系数矩阵
单元的位移场为线性关系，由几何函数矩阵Be可知，由于△ 是常系数，因而Be、Se为常系数矩阵，不随X、Y的变化， 即这种单元在单元内任意一点的应变和应力都相同，因此， 三节点三角形单元称为常应变单元。在应变梯度较大的部 位，单元划分应适当密集，否则将不能真实反映应变的变化 而导致误差较大。
由节点位移条件可求得待定系数：
1 a1 = uj xj yj 2Δ um xm ym 1 a3 = 1 xj uj 2Δ 1 xm um 1 xi ui
ui xi yi
1 a2 = 1 uj yj 2Δ 1 u m ym 1 xi yi 2Δ = 1 x j y j 1 xm ym
1 ui
yi
1 a4 = vj xj yj 2Δ vm xm ym 1 a6 = 1 xj vj 2Δ 1 xm vm 1 xi vi
第四章
连续体平面问题
杆梁结构系统由于本身存在有自然的连接关系 即自然节点，所以他们的离散化均叫做自然离 散，这样的计算模型对原始结构具有很好的描 述，而连续体结构不同，它本身内部不存在有 自然的连接关系，而是以连续介质的形式进行 物质间的相互关联，所以，必须人为地在连续 体内部和边界上划分节点，以分片（单元）连 续的形式来逼近原来复杂的几何形状，这种离 散过程叫做逼近性离散。
N（x,y）为形状函数：
⎡ Ni 0 N j 0 N m 0 ⎤ N ( x, y ) = ⎢ ⎥ ⎢ ⎣ 0 Ni 0 N j 0 N m ⎥ ⎦
常应变三角形单元
宜
不宜
5.相邻单元的尺寸尽可能接近。 6.结点所连接的单元个数尽可能一致。 宜 不宜
§3.3 有限元分析应注意的问题和结果整理
一.结点的选择和单元划分 1.集中力作用点、分布力突变点、支承点应选作结点。 2.不同厚度、不同材料的部分不应划在同一个单元。 3.应力变化大处单元应密集一些。结点的多少与疏密要考虑计算 机的容量和计算精度。 4.单元边界的边长之比应尽可能靠近1。
宜
不宜
5.相邻单的尺寸尽可能接近。 6.结点所连接的单元个数尽可能一致。 宜 二.结点编码 尽可能使相关结点的结点编码差值最小. 1 2 3 4 5 6 7 1 3 5 7 9 11 13 不宜
8
9
10
11 12
13
14
2
4
6
8
10
12
14
总刚半带宽=(相关结点最大差值+1)*结点位移数 总刚半带宽=(7+1)*2=16 总刚半带宽=(2+1)*2=6 总刚需占用的存贮空间为: 6 *14*2=168
结点4的应力由结点1、2、3的应力 外插得到
F
E
2.两单元平均法 三角形单元时，以两相邻单元应力平均值作为边中点的应力近似值。 矩形单元时，以两相邻单元公共边两端结点四个应力的平均值作为边中点 的应力近似值。对于边界处的结点，同样由内结点结果的外插得到。
总刚需占用的存贮空间为: 16 *14*2=448
三.充分利用结构的对称性 P P P P P
四.应力结果的整理
位移的计算结果一般比应力、内力结果精度高。位移达到满意结果， 由几何方程求应变，再由物理方程求应力，结果的精度较差。上述三角形 单元为常应力，矩形单元应力线性变化，而工程问题的应力是比较复杂的。 为更好地反应实际应力情况，需要对计算结果进行整理。常用处理方法有 两种:绕结点平均法和两单元平均法。
三角形平面问题
ci xm x j 0
b j ym yi 0
ai x j ym xm y j 0 bi y j ym a
a j xm yi xi ym 0
c j xi xm a
am xi y j x j yi a 2 bm yi y j a
• 4、应力、应变矩阵 • 将位移函数代入平面问题几何方程，得应变矩阵：
u x x x v y 0 y xy u v x y y 0 Ni y 0 x ui v i 0 u j Nm vj um vm
平面问题的有限单元法
1 有限单元法的计算步骤
弹性力学平面问题的有限单元法包括五个主要步骤： 1、所分析问题的数学建模 2、离散化 3、单元分析 4、整体分析与求解 5、结果分析
图 1
2 平面问题的常应变(三角形)单元
有限单元法的基础是用所谓有限个单元的集合体 来代替原来的连续体，因而必须将连续体简化为由 有限个单元组成的离散体。对于平面问题，最简单， 因而最常用的单元是三角形单元。因平面问题的变 形主要为平面变形，故平面上所有的节点都可视为 平面铰，即每个节点有两个自由度。单元与单元在 节点处用铰相连，作用在连续体荷载也移置到节点 上，成为节点荷载。如节点位移或其某一分量可以 不计之处，就在该节点上安置一个铰支座或相应的 连杆支座。如图1
3)单元内任一点的三个形函数之和恒等于1。
Ni ( x, y) N j ( x, y) Nm ( x, y) 1
4)形函数的值在0—1间变化。
数值模拟：第五讲平面问题(二)——三角形单元分析
u(x, y) a1 a2x a3 y v(x, y) a4 a5x a6 y
a1 ~ a6 为待定系数，称为广义坐标。
由于有限元法中未知量是节点位移，所以上面单 元位移模式需要转换为以节点位移分量为待定参量的 形式。过程如下：
位移多项式写成矩阵形式：
代入各节点取值条件后：
u 1 x
y
aa12
坐标取节点值
a3
uuml
1 1
xl xm
un 1 xn
yl ym
aa12
yn a3
v 1 x
y
aa54
a6
坐标取节点值
vvml vn
1 1 1
xl xm xn
yl ym yn
aaa654
分别解出6个待定系数：
2) 理论和数值试验均可证明，用该单元求解问题时，误差随单元 尺寸的减小和单元数目的增加而减小，这就是有限元解的收敛 性。事实上，当单元变得越来越小时，结构在一个小单元区域 上的应力趋于均匀，而大量小单元拼接在一起就可能很精确地 模拟结构中变化的应力。
5.2.5 单元刚度方程与刚度矩阵
• 前面已经得到了用单元节点 位移表示的单元内部位移、 应变、应力的分布：
•
三角形顶点设为节点，其局部编号为l，m，n（逆时针）。
每节点有总体坐标x，y方向两个待求位移分量：u，v。单元共
有6个位移分量——6个自由度。
（1）单元节点位移列阵
l
e
m
ul
vl
um
vm
un
vn T
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
1
1 uj 2 um
ui
xi xj xm
1 y j , 2 1 uj 2 ym 1 um
yi
1
ui
1 y j , 3 1 xj 2 ym 1 xm
yi
1
xi
ui uj um
(d)
u 1 2x 3y v 4 5x 6 y
(b)
图1 平面三角形单元 将 (d) 式代入 (b) 式的第一式，经整理后得到
图1 弹性体和有限元计算模型
图1 平面三角形单元
第二节 单元分析
一、由结点位移求内部任一点的位移
首先，建立以单元节点位移表示单元内各点位移的 关系式。设单元 e 的节点编号为 i 、 j 、 m 。由弹性力学平 面问题可知，每个节点在其单元平面内的位移可以有两 个分量，所以整个三角形单元将有六个节点位移分量， 即六个自由度。用列阵可表示为：
6. 求解修改后的整体结构刚度方程 考虑整体结构的约束情况，修改整体刚度方程之后，（36）式就变成以节点位移为未知数的代数方程组。解此方程组 可求出节点位移。 解出整体结构的节点位移列阵 后，再根据单元节点的 e e 编号找出对应于单元的位移列阵 ，将 代入（3-3）式就 可求出各单元的应力分量值。 7. 由单元的节点位移列阵计算单元应力
(g)
bm 1 N j x, y a j b j x c j x xi yi 2 cm bm c j 1 x xi a j b j xi c j yi b j x xi 2 cm
有限单元 离散化 集合 总体分析解 有限元法——连续体——单元——代替原连续体 （近似法） （单元分析） 线性方程组
三、有限元法算题的基本步骤 1. 力学模型的选取
（平面问题，平面应变问题，平面应力问题，轴对称问题， 空间问题，板，梁，杆或组合体等，对称或反对称等） y 例如：
Байду номын сангаас
x
为平面应力问题 ，由于结构的对 称性可取结构的 1/4来研究，故 所取的力学模型
T T e iT T j m
ui
vi
uj
vj
um
vm
T
(3-7)
其中的子矩阵
i ui
T vi
（i，j，m 轮换） (a)
式中 ui、vi 是节点i在x轴和y轴方向的位移。
选择一个单元位移模式，单元内各点的位移可按此位移 模式由单元节点位移通过插值而获得。线性函数是一种 最简单的单元位移模式，故设
1 ai bi x j ci y j 0 2 1 N i xm , ym ai bi xm ci ym 0 2 Ni x j , y j
类似地有
N m xi , yi 0 , N m
N j xi , yi 0 , N j x j , y j 1 , N j xm , ym 0
u
1 ai bi x ci y ui a j b j x c j y u j am bm x cm y um 2
(e)
其中
ai
xj xm 1
yj ym yj
x j y m xm y j y j ym
bi
1 ym
（i , j , m轮换） (3-9)
f N e
（3-1）
f ——单元内任一点的位移列阵； e ——单元的结点位移列阵； N ——单元的形函数矩阵；（它的元素是任一点位置坐
标的函数）
4. 单元的力学特性分析 把（3-1）式代入几何方程可推倒出用单元结点位移表示 的单元应变表达式：
B
同理可得
v
1 xj ci x j xm 1 xm
1 ai bi x ci yvi a j b j x c j y v j am bm x cm yvm 2
(f)
若令
Ni 1 ai bi x ci y 2
（i , j , m轮换） (3-10)
ui 1 2 xi 3 yi uj 1 2xj 3yj um 1 2 x m 3 y m
ui xi xj xm yi
, , ,
v j 4 5 xi 6 yi v j 4 5x j 6 y j vm 4 5 xm 6 ym
j j
x
, y 0 ,
N m xm , ym 1
(d)
⒉ 在单元的任一节点上，三个形函数之和等于1，即
N i x , y N j x , y N m x , y 1 ai bi x ci y a j b j x c j y a m bm x cm y 2 1 ai a m a m bi b j bm x ci c j cm y 2 1
8. 计算结果输出 求解出整体结构的位移和应力后，可有选择 地整理输出某些关键点的位移值和应力值，特别
要输出结构的 变形图、应力图、应变图、结构仿
真变形过程动画图及整体结构的弯矩、剪力图等
等。
第一节 结构的离散化
在运用有限单元法分析弹性力学平面问题时，第一步就 是要对弹性体进行离散化，把一个连续的弹性体变换为一个 离散的结构物。对于平面问题，三角形单元是最简单、也是 最常用的单元，在平面应力问题中，单元为三角形板，而在 平面应变问题中，则是三棱柱。 假设采用三角形单元，把弹性体划分为有限个互不重叠 的三角形。这些三角形在其顶点（即节点）处互相连接，组 成一个单元集合体，以替代原来的弹性体。同时，将所有作 用在单元上的载荷（包括集中载荷、表面载荷和体积载荷）， 都按虚功等效的原则移置到节点上，成为等效节点载荷。由 此便得到了平面问题的有限元计算模型，如图1所示。
有限元法
Finite Element Method
第四章
用常应变三角形单元解 弹性力学问题
第一节 结构的离散化 第二节 单元分析 第三节 外力、应力、应变和位移
有限元法基本思想和解题步骤 一、有限元法的基本思想 假想的把一连续体分割成数目有限的小体（单元）， 彼此间只在数目有限的指定点（节点）相互连结，组成 一个单元的集合体以代替原来的连续体，再在节点上 引进等效力以代替实际作用于单元上的外力。选择一个
e
（3-2）
式中：
——单元内任一点应变列阵； B ——单元的应变矩阵；（它的元素仍为位置坐标的
函数）
再把（３－２）式代入物理方程，可导出用单元结点位 移列阵表示的单元应力表达式：
DB
e
（3-3）
式中：
——单元内任一点的应力列阵； D ——单元的弹性矩阵，（它与材料的特性有关）
(e)
简记为
Ni N j N m 1
(3-22)
这说明，三个形函数中只有二个是独立的。 ⒊三角形单元任意一条边上的形函数，仅与该边的两端节 点坐标有关、而与其它节点坐标无关。例如，在i j 边上， 有
x xi N i x, y 1 x j xi x xi N j x, y x j xi N j x, y 0
这样，位移模式 (e) 和 (f) 就可以写为
u N i ui N j u j N m u m v N i vi N j v j N m v m
也可写成矩阵形式 u f Ni I v (3-11)
N jI
N m I N
u 1 2x 3y v 4 5x 6 y
(b)
式中 1、2、…6是待定常数。因三角形单元共有六个 自由度，且位移函数u、v在三个节点处的数值应该等于 这些点处的位移分量的数值。假设节点i、j、m的坐标分 别为（xi , yi ）、（xj , yj ）、（xm , ym ），代入 (b) 式， 得：
(3-23)
事实上，因i j 边的直线方程方程为
y yi y j xi x j
x xi
bm x xi yi cm
(f)
代入（3-10）式中的Nm (x , y) 和Nj (x , y)，有
bm 1 N m x, y a m bm x cm x xi yi 2 cm 1 am bm xi cm yi 0 2
2. 单元的选取、结构的离散化 根据题目的要求，可选择适当的单元把结构离散化。对于 平面问题可用三角元，四边元等。 例如：
3. 选择单元的位移模式 结构离散化后，要用单元内结点的位移通过插值来获得 单元内各点的位移。在有限元法中，通常都是假定单元的位 移模式是多项式，一般来说，单元位移多项式的项数应与单 元的自由度数相等。它的阶数至少包含常数项和一次项。至 于高次项要选取多少项，则应视单元的类型而定。
简单的函数来近似地表示位移分量的分布规律，建立位
移和节点力之间的关系。 有限元法的实质是：把有无限个自由度的连续体， 理想化为只有有限个自由度的单元集合体，使问题简化
为适合于数值解法的结构型问题。
二、经典解与有限元解的区别：
微分 经 典 解 法 —— （解析法）
数目增到∞ 大小趋于 0
建立一个描述连续体 性质的偏微分方程
cj 和am 、bm 、cm 分别是行列式2的第一行、第二行和第 三行各元素的代数余子式，我们有 ⒈ 形函数在各单元节点上的值，具有“本点是1、它 点为零”的性质，即 在节点i上，
N i xi , yi