有限元第二次作业.docx

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 0 0 0 0 k3, 3 k3,4 k3 ,5 1 1 1 单元 1 的刚度矩阵 ： k 1 k4,1 3 k4,1 4 k41 ,5 ， k1 0 0 k3 ,3 k3, 4 k3,5 0 ； k5,1 3 k51,4 k51,5 0 0 k41,3 k41 ,4 k4,51 0 0 0 k51,3 k51, 4 k51, 4 0 0 0 0 0 0 0 0 0 0 0 0 0 k22,2 k22,3 k22,5 0 k22, 2 k22, 3 0 k22, 5 0 0 k3,2 2 k3,2 0 k3,2 5 0 单元 2 的刚度矩阵： k 2 k32, 2 k32,3 k32,5 ， k 2 3 ； 0 0 0 0 0 0 k5,2 k52,3 k5,52 2 0 k52, 2 k52, 3 0 k52, 4 0 0 0 0 0 0 0 0 0 0 0 0 0 k23,2 k23,5 k23,6 0 k23, 2 0 0 k23.5 k23,6 0 0 0 0 0 0 单元 3 的刚度矩阵： k 3 k53,2 k53,5 k53,6 ， k 3 0 0 0 0 0 ； k63,2 k63,5 k63,6 0 0 k53, 2 0 0 k53,5 k53,6 0 k63, 2 0 0 k63,5 k63,6 k14,1 k14,2 0 0 0 k14,6 k14,1 k14,2 k14,6 k24.1 k24, 2 0 0 0 k24,6 0 0 0 0 0 0 单元 4 的刚度矩阵： k 4 k24,1 k24,2 k24,6 ， k 4 ； k64,1 k64,2 k64,6 0 0 0 0 0 0 0 0 0 0 0 0 k64,1 0 0 0 0 k64, 6 4 k 1 k 2 k 3 k 4 总刚度矩阵： K k e e 1 k14,1 k14,2 0 0 0 k14, 6 k24,1 k22,2 k23, 2 k24,2 k22,3 0 k22,5 k23,5 k23,6 k24,6 0 k32,2 k31,3 k32,3 k31, 4 k3,51 k32,5 0 K 0 k41,3 k41,4 k4,1 5 0 0 0 k52, 2 k53,2 k51,3 k52,3 k51,4 k51,5 k52,5 k53,5 k53,6 k64,1 k63, 2 k64,2 0 0 k63,5 k63,6 k64,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,