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

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

2、量, 约束信息数组% HK DISP% 总体刚度矩阵,结点位移向量%*clear allformat 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,1); % 总结点数 NNODE=fscanf(FP1,%d,1); %单元节点数NFPOIN=fscanf(FP1,%d,1); %受力结点数NVFIX=fsc

3、anf(FP1,%d,1); %约束结点个数LNODS=fscanf(FP1,%f,NNODE,NELEM); % 单元定义： 单元结点号（逆时针）COORD=fscanf(FP1,%f,2,NPOIN); % 结点号 x,y坐标(整体坐标下)FPOIN=fscanf(FP1,%f,3,NFPOIN); % 节点力：结点号、X方向力(向右正)，Y方向力(向上正)FIXED=fscanf(FP1,%d,3,NVFIX); %约束信息数组(n,3) n:受约束节点数目, (n,1):约束点号 %(n,2)与(n,3)分别为约束点x方向和y方向的约束情况,受约束为1否则为0% 刚度矩阵的生成%计算刚

4、度矩阵，并对约束条件进行处理Ke=zeros(2*NNODE,2*NNODE); % 单元刚度矩阵并清零 HK=zeros(2*NPOIN,2*NPOIN); % 张成总刚矩阵并清零%调用子程序 生成单元刚度矩阵for m=1:NELEM %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,

5、4),2); %调用单元刚度矩阵a=LNODS(m,:); %临时向量,用来记录当前单元的节点编号%对总刚度矩阵的处理 for j=1:4 for k=1:4 HK(a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)=HK(a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)+. Ke(j*2-1:j*2,k*2-1:k*2); end endend%对荷载向量进行处理FORCE=zeros(2*NPOIN,1); % 张成总荷载向量并清零for i=1:NFPOIN b1=FPOIN(i,1)*2-1;b2=FPOIN(i,1)*2; %FPION(i,

6、1)为作用点 FORCE(b1)=FPOIN(i,2); %FPION(i,2)为x方向的节点力 FORCE(b2)=FPOIN(i,3); %FPION(i,3)为y方向的节点力end%将约束信息加入总刚,总荷载 for i=1:NVFIX if FIXED(i,2)=1 c1=2*FIXED(i,1)-1; HK(c1,:)=0; %将一约束序号处的总刚列向量清0 HK(:,c1)=0; %将一约束序号处的总刚行向量清0 HK(c1,c1)=1; %将行列交叉处的元素置为1 FORCE(c1)=0; end if FIXED(i,3)=1 c2=2*FIXED(i,1); HK(c2,:)

7、=0; HK(:,c2)=0; HK(c2,c2)=1; FORCE(c2)=0; end endDISP=HKFORCE %计算节点位移向量%求解单元应力stress=zeros(3,NELEM);for m=1:NELEM u(1:8)=0; d=LNODS(m,:); %临时向量,用来记录当前单元的节点编号 for i=1:NNODE u(i*2-1:i*2)=DISP(d(i)*2-1:d(i)*2); %从总位移向量中取出当前单元的节点位移 end D=(E/(1-v*v)*1 v 0;v 1 0;0 0 (1-v)/2;%弹性矩阵 %形成应变矩阵BM BM=zeros(3,8);

8、for i=1:NNODE J=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),COORD(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)=B1

9、i 0;0 B2i;B2i B1i/det(J); end stressm=D*BM*u; stress(:,m)=stressm;endstress %输出应力function Ke=K(E,v,h,x1,y1,x2,y2,x3,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

10、; %加权系数W=W1 W2 W3;r=15(1/2)/5;x=-r 0 r;%积分点for i=1:3 for j=1:3 B=eleB(x1,y1,x2,y2,x3,y3,x4,y4,x(i),x(j); J=Jacobi(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

11、);%求应变矩阵BB=zeros(3,8);for i=1:4 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); B(1:3,2*i-1:2*i)=B1i 0;0 B2i;B2i B1i;endB=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

12、=0;y_t=0;for i=1:4 x_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)*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(