平面四边形八节点等参元程序

平面四边形八节点等参元程序====================主程序====================%变量说明2016%NPOINNELEMNVFIX%总结点数单元数受约束结点数%NFPOINYOUNGPOISSTHICK%结点力数弹性模量泊松比厚度%LNODSCOORDFPOINFIXED%单元定义数组结点坐标数组结点力数组约束信息数组%HKFORCEDISP%总刚度矩阵总荷载向量结点位移向量clc;clearall;%清空变量formatshorte%设定输出类型FP1=fopen('sj.txt','rt');%打开数据文件FP1数据文件指针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(FP1,'%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%对行进行循环---按结点号循环N1=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*N1-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=HK\FORCE;%计算单元应力EDISP=zeros(16,1);%单元位移列向量清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*EDISP;%根据S=D*B*EDISP求应力矩阵endfprintf(FP2,'%f\n',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,mu);%调用函数求D弹性矩阵%生成单元刚度矩阵fori=1:3%对高斯积分点进行循环forj=1:3B=BM(ie,LNODS,COORD,X(i),X(j));%调用函数求B应变矩阵J=Jacob(ie,LNODS,COORD,X(i),X(j));%调用函数求雅可比矩阵DJ=J(1,1)*J(2,2)-J(1,2)*J(2,1);%求雅克比矩阵行列式的值EK=EK+W(i)*W(j)*B'*D*B*DJ*h;%高斯计分法求单元刚度矩阵endendreturn===================BM.m===================%求应变矩阵BfunctionB=BM(ie,LNODS,COORD,KC,ETA)NXY=nxy(ie,LNODS,COORD,KC,ETA);%形函数N对x、y分别求导。B=[NXY(1,1),0,NXY(1,2),0,NXY(1,3),0,NXY(1,4),0,NXY(1,5),0,NXY(1,6),0,NXY(1,7),0,NXY(1,8),0;...0,NXY(2,1),0,NXY(2,2),0,NXY(2,3),0,NXY(2,4),0,NXY(2,5),0,NXY(2,6),0,NXY(2,7),0,NXY(2,8);...NXY(2,1),NXY(1,1),NXY(2,2),NXY(1,2),NXY(2,3),NXY(1,3),NXY(2,4),NXY(1,4),NXY(2,5),NXY(1,5),NXY(2,6),NXY(1,6),NXY(2,7),NXY(1,7),NXY(2,8),NXY(1,8)];return===================DM.m===================%求弹性矩阵DfunctionD=DM(E,mu)A1=mu;%平面应力的弹性系数A2=(1-mu)/2;A3=E/(1-mu*mu);%平面应变的弹性系数%A1=POISS/(1-POISS);%A2=(1-2*POISS)/2/(1-POISS);%A3=E*(1-POISS)/(1+POISS)/(1-2*POISS);%弹性矩阵D=A3*[1A10;A110;00A2];return===================NXY.m===================%形函数对x、y求导functionNXY=nxy(ie,LNODS,COORD,KC,ETA)J=Jacob(ie,LNODS,COORD,KC,ETA);%调用函数求雅可比矩阵JJN=[J(2,2),-J(2,1);...%求J矩阵的伴随矩阵-J(1,2),J(1,1)];DJ=J(1,1)*J(2,2)-J(1,2)*J(2,1);%求J矩阵行列式的值%借助复合函数求导法则求NXYNKC=nkc(KC,ETA);%调用函数求N对KC导数NETA=neta(KC,EYA);%调用函数求N对ETA导数NXY=1/DJ*JN*[NKC;NETA];return===================xhs.m===================%平面四节点等参元形函数functionN=xhs(KC,ETA)N(1)=(1-KC)*(1-ETA)*(-KC-ETA-1)/4;N(2)=(1+KC)*(1-ETA)*(KC-ETA-1)/4;N(3)=(1+KC)*(1+ETA)*(KC+ETA-1)/4;N(4)=(1-KC)*(1+ETA)*(-KC+ETA-1)/4;N(5)=(1-KC*KC)*(1-ETA)/2;N(6)=(1-ETA*ETA)*(1+KC)/2;N(7)=(1-KC*KC)*(1+ETA)/2;N(8)=(1-ETA*ETA)*(1-KC)/2;return===================NKC.m===================%形函数对KC求导functionNKC=nkc(KC,ETA)NKC(1)=(-(1-ETA)*(-KC-ETA-1)-(1-KC)*(1-ETA))/4;NKC(2)=((1-ETA)*(KC-ETA-1)+(1+KC)*(1-ETA))/4;NKC(3)=((1+ETA)*(KC+ETA-1)+(1+KC)*(1+ETA))/4;NKC(4)=(-(1+ETA)*(-KC+ETA-1)-(1-KC)*(1+ETA))/4;NKC(5)=-KC*(1-ETA);NKC(6)=(1-ETA*ETA)/2;NKC(7)=-KC*(1+ETA);NKC(8)=-(1-ETA*ETA)/2;return===================NETA.m===================%形函数对ETA求导functionNETA=neta(KC,ETA)NETA(1)=(-(1-KC)*(-KC-ETA-1)-(1-ETA)*(1-KC))/4;NETA(2)=(-(1+KC)*(KC-ETA-1)-(1-ETA)*(1+KC))/4;NETA(3)=((1+KC)*(KC+ETA-1)+(1+ETA)*(1+KC))/4;NETA(4)=((1-KC)*(-KC+ETA-1)+(1+ETA)*(1-KC))/4;NETA(5)=-(1-KC*KC)/2;NETA(6)=-ETA*(1+KC);NETA(7)=(1-KC*KC)/2;NETA(8)=-ETA*(1-KC);return===================Jacob.m===================%求雅可比矩阵JfunctionJ=Jacob(ie,LNODS,COORD,KC,ETA)NKC=nkc(KC,ETA);NETA=neta(KC,ETA);N1=LNODS(ie,1:8)';%存储单元的8个节点号的数组；x0=zeros(8,1);y0=zeros(8,1);x0=COORD(N1(1:8,1),1);%单元节点在xy坐标系中的坐标值y0=COORD(N1(1:8,1),2);J=[NKC*x0,NKC*y0;...NETA*x0,NETA*y0];return实例分析：用平面四节点等参元分析：如图所示悬臂梁，弹性模量：111.2eE；泊松比3.0；长度：m5；高度m1；厚度：mh1。图1悬臂梁尺寸及受力图图2单元划分===============sj.txt（数据文件）==============285132.1e110.31.01320222272128351820426192757161862517267914168241525911121410231324000.50101.50202.50303.50404.5050514.51413.51312.51211.51110.510150.540.530.520.510.500.5120-50000011122112811===================结果文件===================大型有限元软件ANSYS结果：ANSYS计算的位移最大值是：3285.1em10.0000E+0020.0000E+003-3.1667E-054-2.0926E-055-6.1874E-056-7.2788E-057-8.9079E