有限元第二次作业.docx

2-2 图示悬臂板，属于平面应力问题，其网格图及单元、节点编号见图2-1，E=2.1×1011，u=0.28，演算其单刚阵到总刚阵的组集过程，并用MATLAB软件计算总刚阵。图2-1答：根据图2-1所示列出单元节点列表：节点 单元ijk1354225332654162（1）计算单元刚度阵单元1的刚度矩阵： ，；单元2的刚度矩阵：，；单元3的刚度矩阵：，；单元4的刚度矩阵：，；总刚度矩阵：Matlab 程序语言的编写：function Idexglobal gNode gElement gMaterialgNode=[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个结点号。Returnfunction 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); %计算单元刚度矩阵Returnfunction 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