四边形八节点等参元matlab程序全解.docVIP

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