四边形四节点等参元matlab程序

1、悬臂钢梁，尺寸如图一所示； v=0.3 。 h=1， E=2.1e11.图一 悬臂钢梁图二 单元划分与结点编号word 文档 可自由复制编辑附录：% 四边形四节点等参元 matlab 计算程序% 2013 年% 13 级建筑与土木工程 Bruce% B15-405% 本程序只输出了节点位移、单元中心点的应力 %*%变量说明%Evh%弹性模量泊松比厚度%NPOINNELEMNVFIX NNODE NFPOIN%总结点数 ,单元数 ,约束结点个数 , 单元节点数 ,受力结点数%COORDLNODS%结构节点整体坐标数组 ,单元定义数组 ,%FPOINFORCEFIXED%结点力数组，总体荷载向量,

2、 约束信息数组%HKDISP%总体刚度矩阵 ,结点位移向量%*clear all format short e%打开数据文件%弹性模量% 泊松比%厚度%单元数 % 总结点数%单元节点数%受力结点数 %约束结点个数FP1=fopen(' sjd.txt','rt');%读入控制数据E=fscanf(FP1,'%f',1);v=fscanf(FP1,'%f',1);h=fscanf(FP1,'%f',1);NELEM=fscanf(FP1,'%d',1);NPOIN=fscanf(FP1,'%d

3、',1);NNODE=fscanf(FP1,'%d',1);NFPOIN=fscanf(FP1,'%d',1);NVFIX=fscanf(FP1,'%d',1);LNODS=fscanf(FP1,'%f',NNODE,NELEM)' % 单元定义： 单元结点号(逆时针)COORD=fscanf(FP1,'%f',2,NPOIN)' % 结点号 x,y 坐标 (整体坐标下 ) FPOIN=fscanf(FP1,'%f',3,NFPOIN)'% 节点力：结点号、 X 方向

4、力 (向右正 )，Y 方向力 (向上正) FIXED=fscanf(FP1,'%d',3,NVFIX)'%约束信息数组 (n,3) n:受约束节点数目 , (n,1):约束点号word 文档 可自由复制编辑%(n,2)与(n,3)分别为约束点 x 方向和 y 方向的约束情况 ,受约束为 1 否则为 0 % 刚度矩阵的生成 %计算刚度矩阵，并对约束条件进行处理Ke=zeros(2*NNODE,2*NNODE); % 单元刚度矩阵并清零 HK=zeros(2*NPOIN,2*NPOIN); % 张成总刚矩阵并清零 %调用子程序 生成单元刚度矩阵 for m=1:NELEM

5、%m 为单元号Ke=K(E,v,h,.COORD(LNODS(m,1),1),COORD(LNODS(m,1),2),.COORD(LNODS(m,2),1),COORD(LNODS(m,2),2),.COORD(LNODS(m,3),1),COORD(LNODS(m,3),2),.COORD(LNODS(m,4),1),COORD(LNODS(m,4),2);%调用单元刚度矩阵a=LNODS(m,:); %临时向量 ,用来记录当前单元的节点编号 %对总刚度矩阵的处理for j=1:4 for k=1:4HK(a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)=HK(a(j

6、)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)+. Ke(j*2-1:j*2,k*2-1:k*2);endendend % %对荷载向量进行处理FORCE=zeros(2*NPOIN,1); % 张成总荷载向量并清零 for i=1:NFPOINb1=FPOIN(i,1)*2-1;b2=FPOIN(i,1)*2;%FPION(i,1) 为作用点FORCE(b1)=FPOIN(i,2);%FPION(i,2)为x 方向的节点力FORCE(b2)=FPOIN(i,3);%FPION(i,3)为y 方向的节点力end % %将约束信息加入总刚 ,总荷载for i=1:NVFIX i

7、f FIXED(i,2)=1 c1=2*FIXED(i,1)-1;HK(c1,:)=0;%将一约束序号处的总刚列向量清0HK(:,c1)=0;%将一约束序号处的总刚行向量清0HK(c1,c1)=1;%将行列交叉处的元素置为 1word 文档 可自由复制编辑FORCE(c1)=0; end if FIXED(i,3)=1c2=2*FIXED(i,1);HK(c2,:)=0;HK(:,c2)=0;HK(c2,c2)=1;FORCE(c2)=0; end endDISP=HKFORCE%计算节点位移向量%求解单元应力 stress=zeros(3,NELEM);for m=1:NELEMu(1:8)

8、=0;d=LNODS(m,:);%临时向量 ,用来记录当前单元的节点编号for i=1:NNODEu(i*2-1:i*2)=DISP(d(i)*2-1:d(i)*2); %从总位移向量中取出当前单元的节点位移 endD=(E/(1-v*v)*1 v 0;v 1 0;0 0 (1-v)/2;% 弹性矩阵%形成应变矩阵 BMBM=zeros(3,8);for i=1:NNODEJ=Jacobi(COORD(LNODS(m,1),1),COORD(LNODS(m,1),2),.COORD(LNODS(m,2),1),COORD(LNODS(m,2),2),.COORD(LNODS(m,3),1),C

9、OORD(LNODS(m,3),2),.COORD(LNODS(m,4),1),COORD(LNODS(m,4),2),0,0);N_s,N_t=DHS(0,0);B1i=J(2,2)*N_s(i)-J(1,2)*N_t(i);B2i=-J(2,1)*N_s(i)+J(1,1)*N_t(i);BM(1:3,2*i-1:2*i)=B1i 0;0 B2i;B2i B1i/det(J);endstressm=D*BM*u'stress(:,m)=stressm; end stress%输出应力word 文档 可自由复制编辑function Ke=K(E,v,h,x1,y1,x2,y2,x3,

10、y3,x4,y4) %= 单元刚度矩阵 = % E 弹性模量% v 泊松比% h 厚度% x1,y1,x2,y2,x3,y3,x4,y4 为 4 个角结点的坐标%矩阵尺寸为 8 x 8Ke=zeros(8,8);D=(E/(1-v*v)*1 v 0;v 1 0;0 0 (1-v)/2;% 弹性矩阵%高斯积分 采用 3 x 3 个积分点W1=5/9;W2=8/9;W3=5/9; %加权系数W=W1 W2 W3;r=15(1/2)/5;x=-r 0 r;% 积分点for i=1:3for j=1:3B=eleB(x1,y1,x2,y2,x3,y3,x4,y4,x(i),x(j);J=Jacobi(

11、x1,y1,x2,y2,x3,y3,x4,y4,x(i),x(j); Ke=Ke+W(i)*W(j)*B'*D*B*det(J)*h;endendfunction B=eleB(x1,y1,x2,y2,x3,y3,x4,y4,s,t)%调用导函数N_s,N_t=DHS(s,t);%求 Jacobi 矩阵J=Jacobi(x1,y1,x2,y2,x3,y3,x4,y4,s,t);%求应变矩阵 BB=zeros(3,8);for i=1:4B1i=J(2,2)*N_s(i)-J(1,2)*N_t(i);B2i=-J(2,1)*N_s(i)+J(1,1)*N_t(i);B(1:3,2*i-1

12、:2*i)=B1i 0;0 B2i;B2i B1i; endword 文档 可自由复制编辑B=B/det(J);function J=Jacobi(x1,y1,x2,y2,x3,y3,x4,y4,s,t)%Jacobi%单元坐标 x=x1 x2 x3 x4;y=y1 y2 y3 y4;%调用形函数对局部坐标的导数N_s,N_t=DHS(s,t);%求 Jacobi 矩阵的行列式的值 x_s=0;y_s=0;x_t=0;y_t=0;for i=1:4x_s=x_s+N_s(i)*x(i);y_s=y_s+N_s(i)*y(i); x_t=x_t+N_t(i)*x(i);y_t=y_t+N_t(i

13、)*y(i);endJ=x_s y_s;x_t y_t;function N=shape(s,t)% , N(1)=(1-s)*(1-t)/4;N(2)=(1+s)*(1-t)/4;N(3)=(1+s)*(1+t)/4;N(4)=(1-s)*(1+t)/4;function N_s,N_t=DHS(s,t)%形函数求导% , N_s(1)=-1/4*(1-t);N_s(2)=1/4*(1-t);N_s(3)=1/4*(1+t);N_s(4)=-1/4*(1+t);N_t(1)=-1/4*(1-s);N_t(2)=-1/4*(1+s);N_t(3)=1/4*(1+s);N_t(4)=1/4*(1-s);word