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

悬臂钢梁，尺寸如图一所示；v=0.3。h=1，E=2.1e11. 图一悬臂钢梁图二单元划分与结点编号 MatlabMatlab 输出结果输出结果附录Ⅰ：附录Ⅰ：有限元 ANSYS 分析结果采用 PLANE183PLANE183 单元（四边形八节点）单元（四边形八节点）单元得出的结构 Y 向最大位移为-0.216E-04-0.216E-04。。约等于 matlab 平面四边形八节点等参元结点 Y 向最大位移-2.4024E-5-2.4024E-5。。附录Ⅱ：附录Ⅱ： %---------------%---------------四边形八节点等参元四边形八节点等参元 matlabmatlab 计算程序计算程序-------------------------------------------------------- %%——————————————————————主主程程序————序—————————————— %*******************************************************************%******************************************************************* %************************************%************************************ %%20122012 年年 %%本程序只能处理集中荷载作用下的情况本程序只能处理集中荷载作用下的情况 %%只输出了节点位移、单元中心点的应力只输出了节点位移、单元中心点的应力 %*******************************************************************%******************************************************************* %***************%*************** % 变量说明 %Evh %弹性模量泊松比厚度 %NPOINNELEMNVFIXNNODENFPOIN % 总结点数,单元数,约束结点个数,单元节点数 ,受力结点数 %COORDLNODS % 结构节点整体坐标数组,单元定义数组, %FPOINFORCEFIXED % 结点力数组，总体荷载向量, 约束信息数组 %HKDISP %总体刚度矩阵,结点位移向量 %****************************** clear all at short e FP1=( bjd.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=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) 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),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 end end %————————————————————————————————— %对荷载向量进行处理 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;%将一约束序号处的总刚行向量清 0 HK(c1,c1)=1;%将行列交叉处的元素置为 1 FORCE(c1)=