资源描述：
《平面四边形八节点 等参元程序》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。
1、平面四边形八节点等参元程序====================主程序====================%变量说明2016%NPOINNELEMNVFIX%总结点数单元数受约束结点数%NFPOINYOUNGPOISSTHICK%结点力数弹性模量泊松比厚度%LNODSCOORDFPOINFIXED%单元定义数组结点坐标数组结点力数组约束信息数组%HKFORCEDISP%总刚度矩阵总荷载向量结点位移向量clc;clearall;%清空变量formatshorte%设定输出类型FP1=fopen('sj.txt','rt');%打开数据文件FP1数据
2、文件指针FP2=fopen('jg.txt','wt');%写入结果的文件通道号%读入初始数据NPOIN=fscanf(FP1,'%d',1);NELEM=fscanf(FP1,'%d',1);NFPOIN=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(FP
3、1,'%f',[2,NPOIN])';%结点坐标数组FPOIN=fscanf(FP1,'%f',[3,NFPOIN])';%结点力数组FIXED=fscanf(FP1,'%d',[3,NVFIX])';%约束数组HK=zeros(2*NPOIN,2*NPOIN);FORCE=zeros(2*NPOIN,1);%形成总刚度矩阵fori=1:NELEM%对单元个数循环EK=ele_EK(i,LNODS,COORD,YOUNG,POISS,THICK);%调用函数生成单刚%按2*2子块加入总刚中（共计64块）forj=1:8%对行进行循环---按结点号循环N
4、1=LNODS(i,j);fork=1:8%对列进行循环---按结点号循环N2=LNODS(i,k);%单刚2x2子块叠加到总刚中（对号入座）HK((N1*2-1):(N1*2),(N2*2-1):(N2*2))=HK((N1*2-1):(N1*2),(N2*2-1):(N2*2))+EK((j*2-1):(j*2),(k*2-1):(k*2));endendend%由结点力生成总荷载向量列阵fori=1:NFPOIN%对结点荷载个数进行循环N1=FPOIN(i,1);%作用荷载的结点号FORCE((2*N1-1):2*N1,1)=FORCE((2*N
5、1-1):2*N1)+FPOIN(i,2:3)';end%总刚、总荷载进行边界条件处理fori=1:NVFIX%对约束个数进行循环N1=FIXED(i,1);forj=-1:0ifFIXED(i,j+3)==1HK((N1*2+j),1:2*NPOIN)=0;HK(1:2*NPOIN,(N1*2+j))=0;HK((N1*2+j),(N1*2+j))=1;FORCE((N1*2+j),1)=0;%将零位移约束对应行、列变成零，主元素为1。endendend%计算结点位移DISP=HKFORCE;%计算单元应力EDISP=zeros(16,1);%单元
6、位移列向量清0STRES=zeros(3,NELEM);fori=1:NELEM%对单元个数进行循环forj=1:8%对结点进行循环N1=LNODS(i,j);EDISP((j*2-1):(j*2),1)=DISP((2*N1-1):(2*N1),1);%endforj=1:3%对高斯积分点进行循环fork=1:3X=[-sqrt(0.6),0,sqrt(0.6)];B=BM(i,LNODS,COORD,X(j),X(k));%调用函数求应变矩阵BendendD=DM(YOUNG,POISS);%调用函数求弹性矩阵DSTRES(:,i)=D*B*EDI
7、SP;%根据S=D*B*EDISP求应力矩阵endfprintf(FP2,'%f',DISP);fclose(FP1);fclose(FP2);%关闭文件===================ele_EK.m===================%计算单元刚度矩阵functionEK=ele_EK(ie,LNODS,COORD,E,mu,h)EK=zeros(16,16);%8节点的平面单元有8*2个自由度；%3*3高斯积分点和权系数X=[-sqrt(0.6)0sqrt(0.6)];%高斯积分点W=[5/9,8/9,5/9];%加系数权D=DM(E
8、,mu);%调用函数求D弹性矩阵%生成单元刚度矩阵fori=1:3%对高斯积分点进行循环for
