有限元第二次作业
2-2 图示悬臂板， 属于平面应力问题， 其网格图及单元、 节点编号见图 2-1，E=2.1×1011 ， u=0.28，演算其单刚阵到总刚阵的组集过程，并用
MATLAB 软件计算总刚阵。
图 2-1
答：根据图 2-1 所示列出单元节点列表：
节点
j
k
单元
i
1 3 5 4
2 2 5 3
3 2 6 5
4
1
6
2
（1）计算单元刚度阵
0 0 0 0 0 0 1 1 1 0 0
k
3, 3
k
3,4
k
3 ,5
1
1
1
单元 1 的刚度矩阵 ： k 1
k 4,1
3
k 4,1
4
k 41
,5
， k 1
0 0 k 3 ,3 k 3, 4
k 3,5 0 ；
k 5,1
3 k 5
1
,4
k 5
1,5 0 0 k 41,3 k 41 ,4 k 4,51
0 0 k 51,3
k 51, 4
k 51, 4
0 0 0
0 0 0
0 k 22,2 k 22,3 k 22,5 0 k 22, 2 k 22, 3 0 k 22, 5 0
0 k 3,2 2
k 3,
2
0 k 3,2 5 0
单元 2 的刚度矩阵： k
2
k 3
2
, 2
k 32,3 k 32,5 ， k 2 3
；
0 0 0 0 0 0
k 5,
2
k 5
2
,3
k 5,5
2
2
0 k 52, 2 k 5
2
, 3
0 k 52, 4
0 0
0 0 0 0 0 0 k 23,2 k 23,5 k 23,6 0 k 23, 2 0 0 k 23.5 k 23,6
0 0 0 0 0 0 单元 3 的刚度矩阵： k
3
k 53,2
k 5
3
,5
k 5
3
,6 ， k
3
0 0 0 ；
k 63,2 k 63,5 k 63,6
0 k 53, 2 0 0 k 53,5 k 53,6
0 k 63, 2
0 0 k 63,5
k 63,6
k14,1k14,2000k14,6
k14,1k14,2k14,6
k24.1k24, 2000k24,6
000000
单元 4 的刚度矩阵：k4k24,1k24,2k24,6，k4；
k64,1k64,2k64,6000000
000000
k64,10000k64, 6
4
k 1k 2k 3k 4
总刚度矩阵：K k e
e 1
k14,1k14,2000k14, 6
k24,1k22,2 k23, 2 k24,2k22,30k22,5k23,5k23,6k24,6
0k32,2k31,3k32,3k31, 4k3,51k32,50
K
0k41,3k41,4k4,150
0k52, 2k53,2k51,3k52,3k51,4k51,5 k52,5 k53,5k53,6
k64,1k63, 2k64,200k63,5k63,6k64,6
Matlab程序语言的编写：
function Idex
global gNode gElement gMaterial
gNode=[0.0 0.01
0.5 0.01
1.0 0.01
1.0 0.0
0.5 0.0
0.0 0.0]
%gNode 同样是一个矩阵，每一行表示一个结点，第
点的 y 坐标
1 列是结点的x 坐标，第
2 列是结
gElement=[3 4 5
2 3 5
2 5 6
126];
%gElement 是一个矩阵，每一行表示一个单元，第
单元的第 2 个结点号。
1 行是单元的第 1 个结点号，第
2 行是
Return
function k=StiffnessMatrix(ie)
%计算单元刚度矩阵函数
global gNode gElement
k=zeros(6,6);%6x6 单元刚阵
E=2.1*10^11; u=0.28 ;%材料特性%材料特性
t=0.01;%材料特性
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);
d=[1,xi,yi;1,xj,yj;1,xm,ym];
area=det(d);
B=[bi 0 bj 0 bm 0 ;0 ci 0 cj 0 cm;ci bi cj bj cm bm];
B=B/2/area;
D=[1 u 0;u 1 0;0 0 (1-u)/2];
D=D*E/(1-u^2);
k=transpose(B)*D*B*t*abs(area);
阵
Return
function gK=AssembleStiffnessMatrix
%计算总刚阵
global gElement gK ie
gK=zeros(12,12);
for ie =1:1:4
k=StiffnessMatrix(ie);
for i=1:1:3
for j=1:1:3
for p=1:1:2
for q=1:1:2
m=(i-1)*2+p;%每个节点有（i-1)*2+p
n=(j-1)*2+q;%每个节点有（i-1)*2+p
M=(gElement(ie,i)-1)*2+p;
N=(gElement(ie,j)-1)*2+q;
gK(M,N)=gK(M,N)+k(m,n);
end
%计算节点坐标
%计算单元面积
%计算单元刚度矩
%单元循环
%节点循环
%节点循环
%自由度循环
%自由度循环2 个自由度， i 节点的第p 个自由度为
2 个自由度， i 节点的第p 个自由度为
end
end
end
end
Return
则单元 1 的刚度矩阵为
>>StiffnessMatrix(1)
ans =
1.0e+010 *
2.05080-2.05080.04100-0.0410
0 5.69660.0319-5.6966-0.03190
-2.05080.0319 2.0531-0.0729-0.00230.0410
0.0410-5.6966-0.0729 5.69740.0319-0.0008
0-0.0319-0.00230.03190.00230 -0.041000.0410-0.000800.0008
单元 2 的刚度矩阵
>>StiffnessMatrix(2)
ans =
1.0e+010 *
2.0531-0.0729-2.05080.0319-0.00230.0410
-0.0729 5.69740.0410-5.69660.0319-0.0008 -2.05080.0410 2.050800-0.0410
0.0319-5.69660 5.6966-0.03190
-0.00230.03190-0.03190.00230
0.0410-0.0008-0.0410000.0008
单元 3 的刚度矩阵为
>>StiffnessMatrix(3)
ans =
1.0e+010 *
0.00230-0.00230.03190-0.0319
00.00080.0410-0.0008-0.04100
-0.00230.0410 2.0531-0.0729-2.05080.0319
0.0319-0.0008-0.0729 5.69740.0410-5.6966
0-0.0410-2.05080.0410 2.05080 -0.031900.0319-5.69660 5.6966
单元 4 的刚度矩阵
>>StiffnessMatrix(4)
ans =
1.0e+010 *
2.0531-0.0729-2.05080.0319-0.00230.0410
-0.0729 5.69740.0410-5.69660.0319-0.0008
-2.05080.0410 2.050800-0.0410
0.0319-5.69660 5.6966-0.03190
-0.00230.03190-0.03190.00230
0.0410-0.0008-0.0410000.0008
总刚度矩阵为
ans =
1.0e+011 *
Columns 1 through 8
0.2053-0.0073-0.00020.0041000 0
-0.00730.56970.0032-0.0001000 0
-0.00020.00320.4106-0.0073-0.00020.00410
0.0041-0.0001-0.0073 1.13950.0032-0.00010
00-0.00020.00320.20530-0.2051 0.0041
000.0041-0.000100.56970.0032 -0.5697
0000-0.20510.00320.2053 -0.0073
00000.0041-0.5697-0.0073 0.5697
00-0.41020.00730-0.0073-0.0002 0.0032
000.0073-1.1393-0.007300.0041 -0.0001
-0.20510.00410-0.0073000 0
0.0032-0.5697-0.00730000 0
Columns 9 through 12
00-0.20510.0032
000.0041-0.5697
-0.41020.00730-0.0073
0.0073-1.1393-0.00730
0-0.007300
-0.0073000
-0.00020.004100
0.0032-0.000100
0.4106-0.0073-0.00020.0041
-0.0073 1.13950.0032-0.0001
-0.00020.00320.20530
0.0041-0.000100.5697
2-3 在平面问题有限元分析中，
（1）用到了哪些弹性力学中的基本方程？
答：平衡微分方程、几何方程、相容方程（形变协调方程）。
（2）力的平衡条件是如何满足的？
答：根据能量守恒原理，有外力所作虚功应该等于内力虚功。也就是结构在外载荷作
用下处于平衡状态则在结构上的力在任意虚功位移上所作的虚功之和等于零。
以下是用到的方程：
x xy xz X 0x
u
, xy u v
x y z x y x yx y yz
Y0 x y z
zx zy z Z0 x y z y
z
v,
y
w,
z
yz
zx
v w
z y
w u
x z
（ 3）变形协调条件是如何满足的？
答：对材料进行线弹性和各向同性的假设，用弹性力学中应力-应变之间的关系得到变形协调条件。下面是形变协调方程。
2
2
2
2
x
y
xy
zx
xy
yz
2
x
y 2
x 2
,
y
z
x
y z x y x
2
2
2
2
y
z
yz
xy
yz
zx
2
y
z 2
y 2
,
z
x
y
z x
y z y
2
2
2
yz
xy
2
z
y
zx ,
zx
2
z
x 2
x z x 2 z x z y
x y
2-4 在平面三角形单元中的位移、应变、应力具有什么特征？ 位移特征： (1)必须包含单元的刚体位移；
(2)必须包含单元的常应变状态； (3)必须保证不偏惠各坐标轴； (4)必须保证单元内位移连续。
应力特征： (1)三角形单元其应力仅与单元材料和几何尺寸有关，与节点位移有关，而与单
元内位置坐标无关，也即这类单元内的应力是常量。
(2)三角形单元内应力连续，但在公共边界上应力有突变，密布网格可以减少这
种冲突的不合理性。
应变特征： 由于简单三角形单元取线性位移模式，
其应变矩阵为常数矩阵，
即在这样的位移
模式下，三角形单元内的应变为某一常量。
2-5 在平面三角形单元中，当尺寸逐步缩小，单元中的位移、应变、应力具有什么特征？ 当单元尺寸逐步减小时， 单元各点的应变趋于相等，这时常量应变成为主要成分， 因此，位
移应能反应这种常应变状态， 由于应力矩阵也是常数矩阵， 单元应力也是常量。 但是相邻单
元一般具有不同的力， 在单元的公共边上会有应力突变， 随着单元尺寸的逐步减小，
这种突
变会急剧降低，从而不会妨碍有线单元法的简答收敛于精确解。