有限元第二次作业

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 个自由度为