资源描述：
《平面四边形八节点等参元matlab程序.docx》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。
1、广州大学《有限元方法与程序设计》学院：土木工程学院专业：结构工程姓名：曾一凡学号：**%平面四边形八节点等参元MATLAB程序%变量说明&（2015级——结构工程——曾一凡）%YOUNGPOISSTHICK%弹性模量泊松比厚度%NPOINNELEMNVFIXNFORCE%总结点数单元数约束结点个数受力结点数%COORDLNODSFORCE%结构节点整体坐标数组单元定义数组结点力数组%ALLFORCEFIXEDHKDISP%总体荷载向量约束信息数组总体刚度矩阵结点位移向量%1本程序计算了节点位移和单元中心应力并输出到n
2、onde8out.txt文本里%2在第四页举了一个实例用MATLAB算出结果再用ANSYS算出的结果对比%======================主程序=====================formatshorte%设定输出类型clear%清除内存变量FP1=fopen('nonde8.txt','rt');%打开初始数据文件FP2=fopen('nonde8out.txt','wt');%打开文件存放计算结果NPOIN=fscanf(FP1,'%d',1);%结点数NELEM=fscanf(FP1,'%d'
3、,1);%单元数NFORCE=fscanf(FP1,'%d',1);%作用荷载的结点个数NVFIX=fscanf(FP1,'%d',1);%约束数YOUNG=fscanf(FP1,'%e',1);%弹性模量POISS=fscanf(FP1,'%f',1);%泊松比THICK=fscanf(FP1,'%f',1);%厚度LNODS=fscanf(FP1,'%d',[8,NELEM])';%单元结点号数组（逆时针）COORD=fscanf(FP1,'%f',[2,NPOIN])';%结点号x,y坐标(整体坐标下)FORC
4、E=fscanf(FP1,'%f',[3,NFORCE])';%结点力数组%节点力：结点号、X方向力(向右正)，Y方向力(向上正)FIXED=fscanf(FP1,'%d',NVFIX)';%约束信息：约束对应的位移编码（共计NVFIX组）EK=zeros(2*8,2*8);%单元刚度矩阵并清零HK=zeros(2*NPOIN,2*NPOIN);%张成总刚矩阵并清零X=zeros(1,8);%存放单元8个x方向坐标的向量并清零Y=zeros(1,8);%存放单元8个y方向坐标的向量并清零%--------------
5、------------求总刚-------------------------------fori=1:NELEM%对单元个数循环form=1:8%对单元结点个数循环X(m)=COORD(LNODS(i,m),1);%单元8个x方向坐标的向量Y(m)=COORD(LNODS(i,m),2);%单元8个y方向坐标的向量endEK=eKe(X,Y,YOUNG,POISS,THICK);%调用单元刚度矩阵a=LNODS(i,:);%临时向量,用来记录当前单元的结点编号forj=1:8%对行进行循环---按结点号循环for
6、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)+...EK(j*2-1:j*2,k*2-1:k*2);%单刚子块叠加到总刚中endendendALLFORCE=FOECEXL(NPOIN,NFORCE,FORCE);%调用函数生成荷载向量%-------------------------处理约束--------------------------------for
7、j=1:NVFIX%对约束个数进行循环N1=FIXED(j);HK(1:2*NPOIN,N1)=0;HK(N1,1:2*NPOIN)=0;HK(N1,N1)=1;%将零位移约束对应的行、列变成零，主元变成1ALLFORCE(N1)=0;endDISP=HKALLFORCE;%方程求解，HK先求逆再与力向量左乘得到位移%-------------------------求应力---------------------------------stress=zeros(3,NELEM);%应力向量并清零fori=1:NE
8、LEM%对单元个数进行循环D=(YOUNG/(1-POISS*POISS))*[1POISS0;POISS10;00(1-POISS)/2];%弹性矩阵fork=1:8%对单元结点个数进行循环N2=LNODS(i,k);%取单元结点号U(k*2-1:k*2)=DISP(N2*2-1:N2*2);%从总位移向量中取出当前单元的结点位移B=eBe(
