有限元 第二次作业

2-2 图示悬臂板,属于平面应力问题,其网格图及单元、节点编号见图2-1,E=2、1×1011,u=0、28,演算其单刚阵到总刚阵得组集过程,并用MATLAB 软件计算总刚阵。

图2-1

答:根据图2-1所示列出单元节点列表:

i j k 1 3 5 4 2 2 5 3 3 2 6 5 4

1

6

2

(1)计算单元刚度阵 单元1得刚度矩阵: ,; 单元2得刚度矩阵:,; 单元3得刚度矩阵:,; 单元4得刚度矩阵:,; 总刚度矩阵:

[]⎥⎥⎥⎥⎥⎥⎥

⎥⎦

⎤

⎢⎢⎢⎢⎢

⎢⎢⎢⎣⎡++++++++++++=4

6,636,635

,642

,632,641

,636

,535

,525,515,514,523

,513,53

2,522,515,414,413,425

,315,314,323

,313,322

,346,236,235

,225,223

,242

,232,222,241,246,142

,141

,10

00000

00

000k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k K 节点

单元

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 同样就是一个矩阵,每一行表示一个结点,第1 列就是结点得x 坐标,第2 列就是结点得y坐标

gElement=[3 4 5

2 3 5

2 5 6

1 2 6 ];

%gElement 就是一个矩阵,每一行表示一个单元,第1 行就是单元得第1 个结点号,第2 行就是单元得第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; %每个节点有2个自由度,i节点得第p个自由度为(i-1)*2+p

n=(j-1)*2+q; %每个节点有2个自由度,i节点得第p个自由度为(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

end

end

end

end

Return

则单元1得刚度矩阵为

>> StiffnessMatrix(1)

ans =

1、0e+010 *

2、0508 0 -2、0508 0、0410 0 -0、0410

0 5、6966 0、0319 -5、6966 -0、0319 0

-2、0508 0、0319 2、0531 -0、0729 -0、0023 0、0410

0、0410 -5、6966 -0、0729 5、6974 0、0319 -0、0008

0 -0、0319 -0、0023 0、0319 0、0023 0

-0、0410 0 0、0410 -0、0008 0 0、0008 单元2得刚度矩阵

>> StiffnessMatrix(2)

ans =

1、0e+010 *

2、0531 -0、0729 -2、0508 0、0319 -0、0023 0、0410

-0、0729 5、6974 0、0410 -5、6966 0、0319 -0、0008 -2、0508 0、0410 2、0508 0 0 -0、0410

0、0319 -5、6966 0 5、6966 -0、0319 0

-0、0023 0、0319 0 -0、0319 0、0023 0

0、0410 -0、0008 -0、0410 0 0 0、0008 单元3得刚度矩阵为

>> StiffnessMatrix(3)

ans =

1、0e+010 *

0、0023 0 -0、0023 0、0319 0 -0、0319

0 0、0008 0、0410 -0、0008 -0、0410 0

-0、0023 0、0410 2、0531 -0、0729 -2、0508 0、0319

0、0319 -0、0008 -0、0729 5、6974 0、0410 -5、6966

0 -0、0410 -2、0508 0、0410 2、0508 0

-0、0319 0 0、0319 -5、6966 0 5、6966