《四边形八节点等参元matlab程序》由会员分享，可在线阅读，更多相关《四边形八节点等参元matlab程序（11页珍藏版）》请在金锄头文库上搜索。
1、悬臂钢梁，尺寸如图一所示；v=。h=1，E=. 图一 悬臂钢梁图二 单元划分与结点编号Matlab 输出结果 附录：有限元ANSYS分析结果采用PLANE183单元（四边形八节点）单元得出的结构Y向最大位移为。约等于matlab平面四边形八节点等参元结点Y向最大位移。附录：%-四边形八节点等参元 matlab计算程序-% 主 程 序 %*%*% 2012年 % 本程序只能处理集中荷载作用下的情况% 只输出了节点位移、单元中心点的应力%*%*% 变量说明% E v h% 弹性模量 泊松比 厚度% NPOIN NELEM NVFIX NNODE NFPOIN % 总结点数 , 单元数, 约束结点个
2、数, 单元节点数 ,受力结点数% COORD LNODS % 结构节点整体坐标数组, 单元定义数组, % FPOIN FORCE FIXED% 结点力数组， 总体荷载向量, 约束信息数组% HK DISP% 总体刚度矩阵,结点位移向量%*clear allformat short e FP1=fopen(,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); %
3、 总结点数 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方向力(向右正)，Y方向力(向上正)FIXED=fscanf(FP1,%d,3,NVFIX); %约束信息数组(n,3)
4、 n:受约束节点数目, (n,1):约束点号 %(n,2)与(n,3)分别为约束点x方向和y方向的约束情况,受约束为1否则为0%*%*%=平面应力问题的求解=%*%*% % 刚度矩阵的生成%计算刚度矩阵，并对约束条件进行处理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,3),1),
5、COORD(LNODS(m,3),2),. COORD(LNODS(m,5),1),COORD(LNODS(m,5),2),. COORD(LNODS(m,7),1),COORD(LNODS(m,7),2); %调用单元刚度矩阵a=LNODS(m,:); %临时向量,用来记录当前单元的节点编号%对总刚度矩阵的处理 for j=1:8 for k=1:8 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%对荷载向量进
6、行处理FORCE=zeros(2*NPOIN,1); % 张成总荷载向量并清零for i=1:NFPOIN b1=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 if FIXED(i,2)=1 c1=2*FIXED(i,1)-1; HK(c1,:)=0; %将一约束序号处的总刚列向量清0 HK(:,c1)=0; %将一约
7、束序号处的总刚行向量清0 HK(c1,c1)=1; %将行列交叉处的元素置为1 FORCE(c1)=0; end if FIXED(i,3)=1 c2=2*FIXED(i,1); HK(c2,:)=0; HK(:,c2)=0; HK(c2,c2)=1; FORCE(c2)=0; end end%=%=DISP=HKFORCE %计算节点位移向量%=%=%求解单元应力stress=zeros(3,NELEM);for m=1:NELEM u(1:16)=0; d=LNODS(m,:); %临时向量,用来记录当前单元的节点编号 for i=1:NNODE u(i*2-1:i*2)=DISP(d(i
8、)*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,16); for i=1:NNODE J=Jacobi(COORD(LNODS(m,1),1),COORD(LNODS(m,1),2),. COORD(LNODS(m,3),1),COORD(LNODS(m,3),2),. COORD(LNODS(m,5),1),COORD(LNODS(m,5),2),. COORD(LNODS(m,7),1),COORD(LNODS(m,7),2),0,0
9、); 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); end stressm=D*BM*u; stress(:,m)=stressm;endstress %输出应力function Ke=K(E,v,h,x1,y1,x3,y3,x5,y5,x7,y7)%=单元刚度矩阵=% E 弹性模量 % v 泊松比 % h 厚度 % x1,y1,x3,y3,x5,y5,x7,y7 为4个角结点的坐标%矩阵尺寸为16 x 16Ke=zeros(16,16
