一用三结点三角形平面单元计算平面结构的应力和位移讲解
平面问题的有限元分析及三角形单元的应用
第八章 平面问题的有限元分析及三角形单元的应用第一节 概述分析弹性力学平面问题时，最简单的单元式由三个结点组成的三角形单元。
当用以分析平面应力问题时，可将其视为三角板；当用以分析平面应变问题时，则可式为三棱柱。
各单元在结点处为铰结。
图8-1所示位移悬臂梁离散为三角形单元的组合体以矩阵形式列出弹性力学平面问题的基本量和基本方程。
谈形体所受体力分量可表示为[]Tyxy x p p p p p =⎥⎥⎦⎤⎢⎢⎣⎡= (8-1)所受面力分量可表示为[]Tyxy x p p p p p =⎥⎥⎦⎤⎢⎢⎣⎡= （8-2）体内任一点应力分量可表示为[]T xy y x τδδδ= （8-3）任一点的应变分量可表示为[]T xy y x γεεε= （8-4）任一点的位移分量可表示为[]Tv u =δ （8-5）弹性力学平面问题的几何方程的矩阵表达式为⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎣⎡∂∂+∂∂∂∂∂∂=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=x u y v y v x u xy y x εεεε （8-6） 平面应力问题的物理方程的矩阵表达式为⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡--=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡xy y x xy y x E γεεμμμμτσσ2100010112 （8-7）或简写成为 εσD = （8-8）式中⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡--=2100010112μμμμE D （8-9）称为弹性矩阵。
平面应变问题的物理方程也可写成式（8-8），但须将式（8-9）中的E 换成21μ-E，μ换成21μμ-，因此得出⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎣⎡-----+-=)1(22100011011)21)(1()1(22μμμμμμμμμE D （8-10）平衡微分方程及边界条件也可以用矩阵表示，但弹性力学有限元位移法中，通常用虚功方程代替平衡微分方程和应力边界条件。
虚功方程的矩阵表达式为⎰⎰⎰⎰⎰***=+tdxdy tds p f ptdxdy f T T σε （8-11）式中：[]Tv u f ***=，表示虚位移；[]Txyx x ****=γεεε，表示与虚位移相对应的虚应变。
平面问题有限元法
第三章平面问题有限元法重庆大学机械工程学院一、平面单元一、平面单元矩形单元正方形单元二、三角形三节点单元2.1 单元位移模式xy{}(,,)Ti ii u v i j m δ ={}T TeTT T i j mi i j j m m u v u v u v δδδδ ==节点数:3；自由度自由度（（DOF ): 6节点位移节点位移：：单元位移单元位移：：二、三角形三节点单元三角形三节点单元位移模式123456u x y v x y αααααα=++=++(3-1)节点：i ()i i y x ,()i i v u ,节点：j ()j j y x ,()j j v u ,()m m y x ,()m m v u ,节点：m二、三角形三节点单元将三个节点的坐标和位移代入将三个节点的坐标和位移代入（（3-1），），得得ii i y x u 321ααα++=jj j y x u 321ααα++=mm m y x u 321ααα++=321ααα，，ii i y x v 654ααα++=j j j y x v 654ααα++=mm m y x v 654ααα++=654ααα，，二、三角形三节点单元mm m j j ji i iy x v y x v y x v A214=αmmmj j ji i i y x u y x u y x u A211=αm mj j i i y u y u y u A111212=αmm j ji i u x u x u x A111213=αmm j ji i y v y v y v A111215=αmm j j i i v x v x v x A111216=αmmj ji iy x y x y x A 1112=(3-2)二、三角形三节点单元将（3-2）代入代入（（3-1），），并整理并整理i i j j m m i i j j m m u N u N u N u v N v N v N v =++=++(3-3))(2111121y c x b a A y x y x yxAN i i i m mj ji ++==)(2111121y c x b a A y x y xy x AN j j j mmii j ++==)(2111121y c x b a A yxy x y x AN m m m j jiim ++==m j i N N N ,,称为形函数二、三角形三节点单元jm m j mmj j i y x y x y x y x a −==m i i m ii mmj y x y x y x y x a −==ij j i jji i m y x y x y x y x a −==mj mji y y y y b −=−=11im im j y y y y b −=−=11ji jim y y y y b −=−=11二、三角形三节点单元)(m j mj i x x x x c −−==)(11i m im j x x x x c −−==)(11j i jim x x x x c −−==二、三角形三节点单元将（3-3）写成矩阵形式：{}}}ee v uf =(3-4)形函数的性质1）形函数形函数在节点处的值为处的值为11，在其余节点处之值为零i N i≠==ij i j y x N j j i 01),((3-5)mj N N ,？？形函数的性质2）在单元内任一点的三个形函数之和等于在单元内任一点的三个形函数之和等于在单元内任一点的三个形函数之和等于11(3-6)1i j m N N N ++=3）在单元某一边上的形函数与第三个顶点的坐标无关形函数的性质0),(=y x N m )/()()/()(),(i j j i j j i y y y y x x x x y x N −−=−−=)/()()/()(),(i j i i j i j y y y y x x x x y x N −−=−−=在边上ij形函数的性质4）形函数在单元面积形函数在单元面积A A 上的二重积分之值上的二重积分之值，，等于高为等于高为11、底为底为A A 的三角锥的体积的三角锥的体积。
有限元 2-弹性力学平面问题(2.1三角形单元,2.2问题讨论)
有限单元法
土木工程学院
P-7/78
平面应变问题
(2)平面应变问题：图示一圆形隧洞的横截面。由于 隧洞的长度比直径大得多．而荷载又都与0xy平面平 行，且沿z轴为均匀分布，因此可以认为沿z轴方向 的位移分量等于零。这种问题称为平面变形问题。
隧道
挡土墙
有限单元法
土木工程学院
P-8/78
位移与应变之间的几何方程
E xy 2(1 )
1
x 1 y
xy G xy
用矩阵表示为：
x 1 E y 2 xy 1 0 0 x x 1 0 y [ D] y 1 xy xy 0 2
有限单元法
土木工程学院
P-13/78
1 x (1 2 ) x (1 ) y E E E 1 1 令 2 1 1
1 2 E
y x 1
1 x 1 y 则 E1 1 y 1 x 同理 y E1 xy 2(1 ) 2(1 2 ) 1 xy xy xy G E E (1 )
由此可以看出：两类平面问题可以同样的方法进 行分析，只须将相应的E与换成E1或1，或者将 弹性矩阵[D]换成
有限单元法
土木工程学院
P-15/78
常用的平面单元
三角形单元
六结点三角形单元
四结点四边形单元
八结点曲线四边形元
有限单元法
土木工程学院
P-16/78
第2章 弹性力学平面问题有限单元法
如三节点三角形单元
1 Ai 1 x j 2 1 xm
1
x
1 y j ( ai bi x ci y ) 2 ym
y
（7-23）
ai x j ym xm y j bi y j ym ci x j xm
由式（7-23），式（7-21）化为 1 Li (ai bi x ci y ) (i, j , m) 2A
第7章
平面问题高阶单元
7.1 位移模式阶次的选择 在前面两章中讨论了平面问题三结点三角形单元，其位 移模式的最高阶是坐标 x、 y的一次项。这种位移模式导致 单元常应变、常应力特性，单元应变矩阵、应力矩阵、刚 度矩阵均为常数矩阵，因此计算非常简单。但这种单元难 以反映应力梯度的迅速变化。要想提高计算精度，必须细 分网格，增加单元数和点数，因而加大输入数据的工作量 提高计算精度的另一条有效途径是采用高阶单元。由于 高阶单元的应变、应力不再是常数，因此采用少量单元就
Pmi, Pijd的面积。这三个比值Li, Lj, Lm称为P点的面积坐标。 由于 则
Ai Ai Am A
Li L j Lm 1
（7-22）
由此可见，P点的三个面积坐标不是独立的。同时，面积
坐标只是用以确定三角形内部某点的位置，因而是一种局 部坐标。下面进一步给出面积坐标的几个性质。 （2）面积坐标与直角坐标的关系 在图7-5中，三角形Pjm的面积为
各待定系数(a1… a8)。将这些系数再代入式（7-1），可得：
u N i ui N j u j N l ul N mum v N i ui N j u j N l ul N mum
式中形函数为 ：
1 x y N i 1 1 4 a b 1 x y N j 1 1 4 a b 1 x y N l 1 1 4 a b 1 x y N m 1 1 4 a b
第六讲 平面问题(三)——三角形单元综合举例、收敛准则
5.3 简单三角形单元综合举例
• • 图示一平面应力问题。结构为直角三角形薄片，厚度为h。承受 集中载荷P。 有限元离散化结构如图所示。直角三角形薄片的每边中点取为 节点，共划分4个单元、6个节点，编号如图。 各单元节点定义如下表：
•
单元 1 2 3 4
l 1 2 2 3
1 kll 1 kml 1 knl 0 0 0 1 klm 1 kln 1 1 kmm kmn 1 1 knm knn
0 0 0
0 0 0
0 0 0 0 0 0 2 2 0 0 0 0 kll kln 2 2 0 0 0 + 0 knl knn 0 0 0 0 0 0 2 2 0 0 0 0 kml kmn 0 0 0 0 0 0
由于单元3、4跟单元1的几何形 状和局部节点编号顺序完全相同， 因此单元刚度矩阵相等：
[k ]3 = [k ]4 = [k ]1
这里将单元刚度矩阵子块的局部 编号和整体编号对照后，可以方 便总刚度矩阵的叠加！
二、整体刚度矩阵的叠加 1) 单元刚度矩阵扩大成整体规模： 先以单元2为例。 单元 1 2 3 4 l 1 2 2 3 m 2 5 4 5 n 3 3 5 6
整体刚度矩阵列式中各子块的局部编号改为整体编号：
1 1 1 kll klm kln 0 0 0 1 1 2 3 1 2 3 2 3 kml kmm + kll + kll kmn + kln klm klm + kln 0 1 1 2 1 2 4 2 4 4 knl knm + knl knn + knn + kll 0 knm + klm kln [K] = 3 3 3 kml 0 kmm kmn 0 0 2 3 2 4 3 2 3 4 4 0 kml + knl kmn + kml knm kmm + knn + kmm kmn 4 4 4 0 knl 0 knm knn 0
有限元程序设计--第五章 线性三角形单元
16
单元应变和应力矩阵
由于同一单元中的D、B矩阵都是常数矩阵，所以S矩阵也是常 数矩阵。也就是说，三角形三节点单元内的应力分量也是常 量。 当然，相邻单元的E， ， A和bi、ci(i，j， m)一般不完全相同， 因而具有不同的应力，这就造成在相邻单元的公共边上存在 着应力突变现象。但是随着网格的细分，这种突变将会迅速 减小。
b2 B2 0 c2
b1 , c1 , b2 , c2 , b3 , c3 与x、y无关，都是常量，因此B矩
阵也是常量。单元中任一点的应变分量是 B矩阵与单元节点位
移的乘积，因而也都是常量。因此，这种单元被称为常应变单
元。
03:25
15
单元应变和应力矩阵
平面应力： x σ ( x, y ) y
03:25
12
形函数的性质
当单元的位移函数满足完备性要求时，称单元是完备的（通
常较容易满足）。当单元的位移函数满足协调性要求时，称 单元是协调的。
当势能泛函中位移函数的导数是2阶时，要求位移函数在单元
的交界面上具有C1或更高的连续性，这时构造单元的插值函 数往往比较困难。在某些情况下，可以放松对协调性的要求 ，只要单元能够通过分片试验 (Patch test)，有限元分析的解 答仍然可以收敛于正确的解。这种单元称为非协调单元。
3
平面三角形单元
显然，三角形三个节点的的位移可由下列方程给出，在各节点上 的水平位移方程为： u1=1+2 x1 +3 y1 u2=1+2 x2 +3 y2 u3=1+2 x3 +3 y3
1 1 x1 解出 2 1 x2 1 x 3 3 y1 y2 y3
3-6 解题步骤
节点 坐标
１ ０ ０
２ １ ０
３ １ １
４ ０ １
x y
表4-1
1，2，3）
b2 y 3 y1 1 b3 y1 y 2 0
c1 x 2 x 3 0
1
c 2 x 3 x1 1
c 3 x1 x 2 1
K11 1 K K 21 K 31
K12 K 22 K 32
单元的刚度矩阵 K 表示为
2
K
2
K 22 K 42 K 32
K 24 K 44 K 34
K 23 K 43 K 33
2
对于单元
b2 y 4 y 3 1 0 1 c 2 x3 x 4 0 0 0 b4 y 3 y 2 0 1 1 c 4 x 2 x3 2 0 2 b3 y 2 y 4 1 1 0 c3 x 4 x 2 0 2 2
第四章 平面问题的有限元分析
例4-1 §4-7 计算实例
图示为一厚度t的均质长方形薄板，上下受均匀拉力q=106N/m ， 左端固定。材料弹性模量为E，泊松比 u，不记自重，试用有限元 法求其应力分量。
1、离散结构物
为了计算简单，划分为2个单元， 单元号和节点编号如图所示。
y 4 (0,1) ② 1m ① o 3 (0,0) 2m 1 (2,0) x t 2 (2,1)q (KN/m)
2
S
2
2
S2
S4
u 2 v 2 u 4 S3 v4 u 3 v3
一用三结点三角形平面单元计算平面结构的应力和位移讲解
一：用三结点三角形平面单元计算平面结构的应力和位移。
1，设计说明书计算简图，网格划分，单元及结点的编号如下图所示。
由于结构对称，去四分之一结构分析。
其中E=2e10pa,mu=0.167,h=1m.变量注释：Node ------- 节点定义gElement ---- 单元定义gMaterial --- 材料定义，包括弹性模量，泊松比和厚度gBC1 -------- 约束条件gNF --------- 集中力gk------------总刚gDelta-------结点位移子程序注释：PlaneStructualModel ———定义有限元模型SolveModel ———————求解有限元模型DisplayResults ——————显示计算结果k = StiffnessMatrix( ie )———计算单元刚度AssembleStiffnessMatrix( ie, k )—形成总刚es = ElementStress( ie )————计算单元应力function exam1% 输入参数：无% 输出结果：节点位移和单元应力PlaneStructualModel ; % 定义有限元模型SolveModel ; % 求解有限元模型DisplayResults ; % 显示计算结果return ;function PlaneStructualModel% 定义平面结构的有限元模型% 输入参数：无% 说明：% 该函数定义平面结构的有限元模型数据：% gNode ------- 节点定义% gElement ---- 单元定义% gMaterial --- 材料定义，包括弹性模量，泊松比和厚度% gBC1 -------- 约束条件% gNF --------- 集中力global gNode gElement gMaterial gBC1 gNF% 节点坐标% x ygNode = [0.0, 2.0 % 节点10.0, 1.0 % 节点21.0, 1.0 % 节点30.0, 0.0 % 节点41.0, 0.0 % 节点52.0, 0.0] ; % 节点6% 单元定义% 节点1 节点2 节点3 材料号gElement = [3, 1, 2, 1 % 单元15, 2, 4, 1 % 单元22, 5, 3, 1 % 单元36, 3, 5, 1]; % 单元4 % 材料性质% 弹性模量泊松比厚度gMaterial = [1e0, 0, 1] ; % 材料1% 第一类约束条件% 节点号自由度号约束值gBC1 = [ 1, 1, 0.02, 1, 0.04, 1, 0.04, 2, 0.05, 2, 0.06, 2, 0.0] ;% 集中力% 节点号自由度号集中力值gNF = [ 1, 2, -1] ;returnfunction SolveModel% 求解有限元模型% 输入参数：无% 说明：% 该函数求解有限元模型，过程如下% 1. 计算单元刚度矩阵，集成整体刚度矩阵% 2. 计算单元的等效节点力，集成整体节点力向量% 3. 处理约束条件，修改整体刚度矩阵和节点力向量% 4. 求解方程组，得到整体节点位移向量global gNode gElement gMaterial gBC1 gNF gK gDelta% step1. 定义整体刚度矩阵和节点力向量[node_number,dummy] = size( gNode ) ;gK = sparse( node_number * 2, node_number * 2 ) ;f = sparse( node_number * 2, 1 ) ;% step2. 计算单元刚度矩阵，并集成到整体刚度矩阵中[element_number,dummy] = size( gElement ) ;for ie=1:1:element_numberk = StiffnessMatrix( ie ) ;AssembleStiffnessMatrix( ie, k ) ;end% step3. 把集中力直接集成到整体节点力向量中[nf_number, dummy] = size( gNF ) ;for inf=1:1:nf_numbern = gNF( inf, 1 ) ;d = gNF( inf, 2 ) ;f( (n-1)*2 + d ) = gNF( inf, 3 ) ;end% step4. 处理约束条件，修改刚度矩阵和节点力向量。
第七章 平面问题的有限单元法(Q4)
8
4节点四边形单元
y, v
u1 v 1 u2 u de 2 u3 u3 u4 u 4 displacements at node 1 displacements at node 2 displacements at node 3 displacements at node 4
x 1 2 3 4 N1 x1 N 2 x2 N 3 x3 N 4 x4 y 1 2 3 4 N1 y1 N 2 y2 N 3 y3 N 4 y4
1 N (1 )(1 ) 1 4 N 1 (1 )(1 ) 2 4 1 N (1 )(1 ) 3 4 N 1 (1 )(1 ) 4 4
1 4
Nj 1 4 (1 j )(1 j )
4 ( 1, +1) ( u4, v4)
1
N3 1 4 (1 )(1 ) N4 1 4 (1 )(1 )
N 3 at node 1 1 4 (1 )(1 ) 1 0 N 3 at node 2 1 4 (1 )(1 ) 1 0
同理：
1 1 1 1 1 y1 2 1 1 1 1 1 y2 1 1 1 1 4 3 y3 1 1 1 1 y4 4
K e B DBtd
e
T
11
等参单元
对于一般的四边形单元，在总体坐标系下构造 位移插值函数，则计算形状函数矩阵、单元刚 度矩阵及等效节点载荷列阵时十分冗繁；而对 于矩形单元，相应的计算要简单的多。 矩形单元明显的缺点是不能很好的符合曲线边 界，因此可以采用矩形单元和三角形单元混合 使用（网格划分困难）。更为一般的方法是通 过等参变换将局部自然坐标系内的规格化矩形 单元变换为总体坐标系内的任意四边形单元（ 包括高次曲边四边形单元）。 等参单元的提出为有限元法成为现代工程实
结构力学练习题及答案
一．是非题（将判断结果填入括弧：以O 表示正确，X 表示错误）（本大题分4小题，共11分）1 . (本小题 3分)图示结构中DE 杆的轴力F NDE =F P /3。
（ ）.2 . (本小题 4分)用力法解超静定结构时，只能采用多余约束力作为基本未知量。
（ ）3 . (本小题 2分)力矩分配中的传递系数等于传递弯矩与分配弯矩之比，它与外因无关。
（ ）4 . (本小题 2分)用位移法解超静定结构时，基本结构超静定次数一定比原结构高。
（ ）二．选择题（将选中答案的字母填入括弧内）（本大题分5小题，共21分） 1 （本小题6分）图示结构EI=常数，截面A 右侧的弯矩为：（ ）A ．2/M ；B ．M ；C ．0； D. )2/(EI M 。
2. （本小题4分）图示桁架下弦承载，下面画出的杆件内力影响线，此杆件是：（ ） Ａ.ch; B.ci; C.dj; D.cj.F p /2M2a2a a aa aA F p /2F p /2 F p /2F p F pa a aa F PED3. (本小题 4分)图a 结构的最后弯矩图为：A. 图b;B. 图c;C. 图d;D.都不对。
（ ）( a) (b) (c) (d)4. (本小题 4分)用图乘法求位移的必要条件之一是： A.单位荷载下的弯矩图为一直线； B.结构可分为等截面直杆段； C.所有杆件EI 为常数且相同； D.结构必须是静定的。
( ) 5. （本小题3分）图示梁A 点的竖向位移为（向下为正）：( ) Ａ．F P l 3/(24EI); B. F P l 3/(!6EI); C. 5F P l 3/(96EI); D. 5F P l 3/(48EI).三（本大题 5分）对图示体系进行几何组成分析。
A l /2l /2EI 2EIF Pa d c eb fgh iklF P =11j llM /4 3M /4M /43M /43M /4M /4M /8 M /2EIEIM四（本大题 9分）图示结构B 支座下沉4 mm ，各杆EI=2.0×105 kN ·m 2，用力法计算并作M 图。
一 三节点三角形单元
有限元课程总结一 三节点三角形单元 1位移函数移函数写成矩阵形式为：确定六个待定系数矩阵形式如下：[]{}em m j j i i m jim j iN v u v u v u N N N N N N v u δ=⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧⎥⎦⎤⎢⎣⎡=⎭⎬⎫⎩⎨⎧000002 单元刚度矩阵的计算1）单元应变和节点位移的关系由几何方程可以得到单元的应变表达式，{}⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎪⎪⎪⎭⎪⎪⎪⎬⎫⎪⎪⎪⎩⎪⎪⎪⎨⎧∂∂+∂∂∂∂∂∂=m m j j i i m mjjiim j i m j i v u v u v u b c b c b c c c c b b b A x v y u y v x u 00000021ε⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧321111a a a y x y x y x u u u m m j j i i m j i ⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧m j i m j i m j i m j i u u u c c c b b b a a a A a a a 21321⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧m j i m j i m jim ji 654v v v c c c b b b a a a A 21a a a2）单元应力与单元节点位移的关系[][][]⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡---==i i i i i ii i b c c b c b A E B D S 2121)1(22μμμμμ),,;,,(21212121)1(4]][[][][2m j i s m j i r b b c c cb bc b c c b c c b b A Et B D B K s r s r sr s r s r s r s r s r s T r rs ==⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-+-+-+-+-==μμμμμμμ3)单元刚度矩阵[][][][][][]Tm j iT m T j T i mm mj mi jm jj ji im ij ii e B B B D B B B tA K K K K K K K K K K ⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=][][][][][][][][][3载荷移置1)集中力的移置如图3所示，在单元内任意一点作用集中力 {}⎭⎬⎫⎩⎨⎧=y x P P P图3由虚功相等可得，()()}{][}{}{}{**P N R TTe eTe δδ=由于虚位移是任意的，则 }{][}{P N R Te =2)体力的移置令单元所受的均匀分布体力为⎭⎬⎫⎩⎨⎧=y x p ρρ}{由虚功相等可得，()()⎰⎰=tdxdy p N R T TeTe }{][}{}{}{**δδ⎰⎰=tdxdy p N R T e }{][}{ 3)分布面力的移置设在单元的边上分布有面力{}TY X P ],[=，同样可以得到结点载荷，⎰=s T e tdsP N R }{][}{4. 引入约束条件，修改刚度方程并求解1)乘大数法处理边界条件图3-4所示的结构的约束和载荷情况，如图3-7所示。
三角形平面问题
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间变化。
现代设计方法4-3 三角形三节点平面单元概要
单元分析的步骤：
节点 (1) 位移
单元内部 各点位移
(2)
单元 (3) 应变
单元 应力
(4)
节点 力
单元分析
2.1 由节点位移求单元内部任一点位移
(1) 单元位移模式
有限元法中，通常采用位移法进行计算，即取节点的位移分量
{a }=[L]{d }e
为三角形单元面积。
1 1 A 1 xj yj 2 2 1 xm ym 1 xi yi
1 ( x j ym xm yi xi y j xm y j xi ym x j yi ) 2
将a写成矩阵形式，有{a }=[L]{d }e
ai b i 1 ci L 0 2 0 0 0 a j 0 am 0 0 b j 0 bm 0 0 c j 0 cm 0 ai 0 a j 0 am bi 0 b j 0 bm ci 0 c j 0 cm
ui xi xj xm yi yj ym
ui 1 xi u j 1 x j u 1 x m m
1 ui yi yj ym
yi a1 yj a 2 a ym 3
1 xi ui
(a)
1 a1 uj A um
v ( x, y ) 0
a1 a 2 x y 0 0 0 a 3 (4.23) m ( x , y ) a a 0 0 1 x y 4 a 5 a 6
(2) 由单元节点位移{d }e求位移参数{a }
4.3 平面问题的有限单元法
数值模拟第五讲 平面问题——三角形单元分析 ppt课件
1 xl yl
其中： 2 1 xm ym 1 xn yn
为三角形数面 值模拟积第五讲 平面问题——三 角形单元分析
节点坐标行列式
9
由第二组方程求解 a4 ~ a6 :
vl 1 xl
vm
1
xm
vn 1 xn
yl a4
ym
a5
yn a6
a a5 4 1 1
xl xm
a6 1 xn
25
角形单元分析
根据上面推导，有：
外力虚功： W(*e)Tpe
虚应变能： U (* e )T V eB T D B d V e
❖由虚功原理得到：
WU
(* e ) Tp e (* e ) T B T D B d V e V e
考虑到虚位移 *e的任意性，由上式立即得到下列方程：
n
•
单元平衡时要在节点处受到节点力（节点对单元的作用
力），每节点有2个节点力分量，单元有6个节点力分量。
（2） 单元节点力列阵
pe p pm l pxl pyl pxm pym pxn pynT
pn
• 下面要研究的问题是该三角形薄片弹性体在保持平衡时所受节
点力和节点位移的关系。
数值模拟第五讲 平面问题——三
数值模拟第五讲 平面问题——三
15
角形单元分析
❖ 根据位移模式表达式及其形函数的性质，可以推断出两个相邻 三角形单元上位移分布形状和公共边界上位移的情况：
uNlul NmumNnun vNlvl NmvmNnvn
两个单元上位移线性连续分布，各单
元在公共边界上位移线性分布，数值
相同——边界位移协调！
由第一组方程求解 a1 ~ a3 :
现代设计方法4-3 三角形三节点平面单元
[m(x,y)][L]10
x 0
y 0
0 1
0 x
0y21c0i
0 ai
cj 0
0 aj
cm 0
0 am
0
bi
0 ci
0 bj 0 cj
0 bm 0 cm
N i( 0 x ,y )N i( 0 x ,y )N j( 0 x ,y )N j( 0 x ,y )N m ( 0 x ,y )N m ( 0 x ,y )
N i(x,y)2 1 (a ib ixciy)1b x
1
y
N j(x,y)2 (aj b jx cjy) 1 a
N m (x ,y ) 2 1 (a m b m x c m y ) 1 b x a y
精品课件
[N ] 1b x 0
0 1x
1y 0 a
0 1y
(1xy) ba 0
0
1
1 A
ui uj
xi xj
yi yj
1
2
1 A
1
ui uj
yi yj
1
3
1 A
1
xi xj
ui uj
um xm ym
1 um ym
1 xm um
(b)
1 xi yi A 1 xj yj
其中，
1 x y精品课件
m
m
对v同理可列出a4、a5、a6的方程 ：
vi = a4+a5 xi+a6 yi
ci bi
(i, j, m)
2.3 由单元节点位移求单元的应力
物理方程
{s }=[D]{e}
而
{e }=[B]{d}e
{s } = [D][B]{d }e
简支梁的有限单元法分析-三角形三节点单元
u x x v y y xy u v y x
B
e
1 Ni ( x, y ) (ai bi x ci y ) 2A
（i , j , m）
y j (0,2 )
x
m (0,0 ) i (2,0 )
三角形三节点单元
代入[D],[B]得三角形单元的单元刚度矩阵：
1 4 0 0 2 Et e [k ] 1 m2 m 4 1 4 m 4 对 1 m 8 1 m 8 0 1 m 8 1 m 8 称 1 m 8 0 1 m 8 1 m 8 1 4 3 m 8
简支梁的有限单元法分析
三角形三节点平面单元
王 峰
有限元分析的基本步骤：
结构离散化
单元分析
整体分析
1 结构离散化
图示为简支梁，梁的厚度为t，泊松比m =0.3，弹性 模量为E=2e+5Mpa，用三节点三角形单元进行离散， 直角三角形边长为2。
2 单元分析
单元分析的主要内容：由节点位移求内部任一点的
边长度，t为单元厚度。 单元等效节点荷载为 P
e
ql0t T 0 1 0 0 0 1 2
j
结点 编号
1 51 ... 2 102 53 52 ... ... 561 ... 511
将单元节点荷载集成为结构的节点荷载列阵[P]
结点载荷列阵
2 4 6 102 104 1122
同理
k11 K 2 n i 1 , 2 n i 1
k32 K 2 n j 1 , 2 n i
分别叙述位移法以及应力法求解步骤
分别叙述位移法以及应力法求解步骤位移法呀，就像是给结构的变形玩一场找规律的游戏。
对于一个结构，咱先得确定基本未知量。
这就好比是找出这个结构里那些关键的会动的点或者方向，一般就是节点的位移啦。
这些位移包括线位移和角位移哦。
比如说一个框架结构，那些节点是怎么歪呀、怎么扭呀的，这就是我们要找出来的基本未知量。
接着呢，要建立位移法的基本结构。
这就像是给这个结构搭一个特殊的架子，把那些有位移的地方给它加上约束，让它先“老实”下来。
然后就是列出位移法的典型方程啦。
这个方程呢，就像是这个结构的变形和受力之间的一个小秘密协议。
方程的左边是关于那些位移未知量的系数，右边呢，就是这个结构受到的荷载在这些约束上产生的反力。
这一步就像是在给结构的变形和受力牵红线，让它们之间的关系明明白白的。
最后就是求解这些方程啦。
求出那些位移未知量之后呢，就可以根据位移和内力的关系，算出结构各个杆件的内力啦。
这样，整个结构在荷载作用下的情况就被我们摸得清清楚楚啦。
应力法呢，也很有趣味的哦。
一开始呀，要确定应力分量的函数形式。
这就像是在猜这个结构里面应力是怎么分布的，是像小波浪一样呢，还是像规规矩矩的小方块一样。
这一步要根据结构的形状、受力情况这些来猜个大概的形式。
然后呢，要根据平衡条件、几何条件和物理条件来建立方程。
平衡条件就像是结构在受力的时候不能自己跑掉或者翻跟头，要安安稳稳的。
几何条件呢，就是结构变形的时候，各个部分之间的形状关系。
物理条件就是应力和应变之间的那种亲密关系啦。
把这些条件都用上，建立起方程，这就像是给应力这个调皮的小家伙定规矩。
接下来就是求解这个方程啦。
这个方程可能有点难搞，不过没关系，就像解开一个神秘的小盒子一样，慢慢算，总能算出应力分量的。
最后呀，根据算出来的应力，还可以进一步算出结构的变形之类的其他情况呢。
这样，整个结构在应力方面的情况就被我们掌握在手中啦。
有限元平面问题
平面应力 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 ⎥ ⎦
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
一：用三结点三角形平面单元计算平面结构的应力和位移。
1，设计说明书计算简图，网格划分，单元及结点的编号如下图所示。
由于结构对称，去四分之一结构分析。
其中E=2e10pa,mu=0.167,h=1m.变量注释：Node ------- 节点定义gElement ---- 单元定义gMaterial --- 材料定义，包括弹性模量，泊松比和厚度gBC1 -------- 约束条件gNF --------- 集中力gk------------总刚gDelta-------结点位移子程序注释：PlaneStructualModel ———定义有限元模型SolveModel ———————求解有限元模型DisplayResults ——————显示计算结果k = StiffnessMatrix( ie )———计算单元刚度AssembleStiffnessMatrix( ie, k )—形成总刚es = ElementStress( ie )————计算单元应力function exam1% 输入参数：无% 输出结果：节点位移和单元应力PlaneStructualModel ; % 定义有限元模型SolveModel ; % 求解有限元模型DisplayResults ; % 显示计算结果return ;function PlaneStructualModel% 定义平面结构的有限元模型% 输入参数：无% 说明：% 该函数定义平面结构的有限元模型数据：% gNode ------- 节点定义% gElement ---- 单元定义% gMaterial --- 材料定义，包括弹性模量，泊松比和厚度% gBC1 -------- 约束条件% gNF --------- 集中力global gNode gElement gMaterial gBC1 gNF% 节点坐标% x ygNode = [0.0, 2.0 % 节点10.0, 1.0 % 节点21.0, 1.0 % 节点30.0, 0.0 % 节点41.0, 0.0 % 节点52.0, 0.0] ; % 节点6% 单元定义% 节点1 节点2 节点3 材料号gElement = [3, 1, 2, 1 % 单元15, 2, 4, 1 % 单元22, 5, 3, 1 % 单元36, 3, 5, 1]; % 单元4 % 材料性质% 弹性模量泊松比厚度gMaterial = [1e0, 0, 1] ; % 材料1% 第一类约束条件% 节点号自由度号约束值gBC1 = [ 1, 1, 0.02, 1, 0.04, 1, 0.04, 2, 0.05, 2, 0.06, 2, 0.0] ;% 集中力% 节点号自由度号集中力值gNF = [ 1, 2, -1] ;returnfunction SolveModel% 求解有限元模型% 输入参数：无% 说明：% 该函数求解有限元模型，过程如下% 1. 计算单元刚度矩阵，集成整体刚度矩阵% 2. 计算单元的等效节点力，集成整体节点力向量% 3. 处理约束条件，修改整体刚度矩阵和节点力向量% 4. 求解方程组，得到整体节点位移向量global gNode gElement gMaterial gBC1 gNF gK gDelta% step1. 定义整体刚度矩阵和节点力向量[node_number,dummy] = size( gNode ) ;gK = sparse( node_number * 2, node_number * 2 ) ;f = sparse( node_number * 2, 1 ) ;% step2. 计算单元刚度矩阵，并集成到整体刚度矩阵中[element_number,dummy] = size( gElement ) ;for ie=1:1:element_numberk = StiffnessMatrix( ie ) ;AssembleStiffnessMatrix( ie, k ) ;end% step3. 把集中力直接集成到整体节点力向量中[nf_number, dummy] = size( gNF ) ;for inf=1:1:nf_numbern = gNF( inf, 1 ) ;d = gNF( inf, 2 ) ;f( (n-1)*2 + d ) = gNF( inf, 3 ) ;end% step4. 处理约束条件，修改刚度矩阵和节点力向量。
采用乘大数法[bc_number,dummy] = size( gBC1 ) ;for ibc=1:1:bc_numbern = gBC1(ibc, 1 ) ;d = gBC1(ibc, 2 ) ;m = (n-1)*2 + d ;f(m) = gBC1(ibc, 3)* gK(m,m) * 1e15 ;gK(m,m) = gK(m,m) * 1e15 ;end% step 5. 求解方程组，得到节点位移向量gDelta = gK \ f ;returnfunction DisplayResults% 显示计算结果% 输入参数：无global gNode gElement gMaterial gBC1 gNF gK gDeltafprintf( '节点位移\n' ) ;fprintf( ' 节点号x方向位移y方向位移\n' ) ;[node_number,dummy] = size( gNode ) ;for i=1:node_numberfprintf( '%6d %16.8e %16.8e\n',...i, gDelta((i-1)*2+1), gDelta((i-1)*2+2)) ;endfprintf( '\n\n单元应力\n' ) ;fprintf( ' X-STR Y-STR XY-STR\n' ) ;[element_number, dummy] = size( gElement ) ;for ie = 1:element_numberes = ElementStress( ie ) ;fprintf( '单元号%6d %16.8e %16.8e %16.8e\n', ...ie, es(1), es(2), es(3));endreturnfunction k = StiffnessMatrix( ie )% 计算单元刚度矩阵% 输入参数:% ie ------- 单元号global gNode gElement gMaterialk = zeros( 6, 6 ) ;E = gMaterial( gElement(ie, 4), 1 ) ;mu = gMaterial( gElement(ie, 4), 2 ) ;h = gMaterial( gElement(ie, 4), 3 ) ;xi = gNode( gElement( ie, 1 ), 1 ) ;yi = gNode( gElement( ie, 1 ), 2 ) ;xj = gNode( gElement( ie, 2 ), 1 ) ;yj = gNode( gElement( ie, 2 ), 2 ) ;xm = gNode( gElement( ie, 3 ), 1 ) ;ym = gNode( gElement( ie, 3 ), 2 ) ;ai=xj*ym-xm*yj;aj=xm*yi-xi*ym;am=xi*yj-xj*yi;bi=yj-ym;bj=ym-yi;bm=yi-yj;ci=-(xj-xm);cj=-(xm-xi);cm=-(xi-xj);A=(ai+aj+am)/2;B=[bi 0 bj 0 bm 00 ci 0 cj 0 cmci bi cj bj cm bm];B=B/2/A;D=[1 mu 0mu 1 00 0 (1-mu)/2];D=D*E/(1-mu^2);k=transpose(B)*D*B*h*abs(A);returnfunction AssembleStiffnessMatrix( ie, k )% 把单元刚度矩阵集成到整体刚度矩阵% 输入参数:% ie --- 单元号% k --- 单元刚度矩阵global gElement gKfor i=1:1:3for j=1:1:3for p=1:1:2for q =1:1:2m = (i-1)*2+p ;n = (j-1)*2+q ;M = (gElement(ie,i)-1)*2+p ;N = (gElement(ie,j)-1)*2+q ;gK(M,N) = gK(M,N) + k(m,n) ;endendendendreturnfunction es = ElementStress( ie )% 计算单元的应力% 输入参数% ie ----- 节点号% es ----- 单元应力global gElement gNode gDelta gMateriales=zeros(1,6);de=zeros(6,1);for j=1:1:3de(2*j-1)=gDelta(2*gElement(ie,j)-1);de(2*j)=gDelta(2*gElement(ie,j));endE = gMaterial( gElement(ie, 4), 1 ) ;mu = gMaterial( gElement(ie, 4), 2 ) ;h = gMaterial( gElement(ie, 4), 3 ) ;xi = gNode( gElement( ie, 1 ), 1 ) ;yi = gNode( gElement( ie, 1 ), 2 ) ;xj = gNode( gElement( ie, 2 ), 1 ) ;yj = gNode( gElement( ie, 2 ), 2 ) ;xm = gNode( gElement( ie, 3 ), 1 ) ;ym = gNode( gElement( ie, 3 ), 2 ) ;ai=xj*ym-xm*yj;aj=xm*yi-xi*ym;am=xi*yj-xj*yi;bi=yj-ym;bj=ym-yi;bm=yi-yj;ci=-(xj-xm);cj=-(xm-xi);cm=-(xi-xj);A=(ai+aj+am)/2;B=[bi 0 bj 0 bm 00 ci 0 cj 0 cmci bi cj bj cm bm];B=B/2/A;D=[1 mu 0mu 1 00 0 (1-mu)/2];D=D*E/(1-mu^2);S=D*B;es(1:3)=S*de;es(6)=0.5*sqrt((es(1)-es(2))^2+4*es(3)^2);es(4)=0.5*(es(1)+es(2))+es(6);es(5)=0.5*(es(1)+es(2))-es(6);return3，数据文件：输入数据：gNode = [0.0, 2.00.0, 1.01.0, 1.00.0, 0.01.0, 0.02.0, 0.0] ;gElement = [3, 1, 2, 15, 2, 4, 12, 5, 3, 16, 3, 5, 1];gMaterial = [1e0, 0, 1] ;gBC1 = [ 1, 1, 0.02, 1, 0.04, 1, 0.04, 2, 0.05, 2, 0.06, 2, 0.0] ;gNF = [ 1, 2, -1] ;输出数据：节点位移节点号x方向位移y方向位移1 -8.79120879e-016 -3.25274725e+0002 8.79120879e-017 -1.25274725e+0003 -8.79120879e-002 -3.73626374e-0014 1.17216117e-016 -8.35164835e-0165 1.75824176e-001 -2.93040293e-0166 1.75824176e-001 2.63736264e-016单元应力X-STR Y-STRXY-STR单元号 1 -8.79120879e-002 -2.00000000e+000 4.39560440e-001单元号 2 1.75824176e-001 -1.25274725e+000 2.56410256e-016单元号 3 -8.79120879e-002 -3.73626374e-001 3.07692308e-001单元号 4 0.00000000e+000 -3.73626374e-001 -1.31868132e-001。