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

1广州大学《有限元方法与程序设计》学院：土木工程学院专业：结构工程姓名：曾一凡学号：21115160**2%平面四边形八节点等参元MATLAB程序%变量说明&（2015级——结构工程——曾一凡）%YOUNGPOISSTHICK%弹性模量泊松比厚度%NPOINNELEMNVFIXNFORCE%总结点数单元数约束结点个数受力结点数%COORDLNODSFORCE%结构节点整体坐标数组单元定义数组结点力数组%ALLFORCEFIXEDHKDISP%总体荷载向量约束信息数组总体刚度矩阵结点位移向量%1本程序计算了节点位移和单元中心应力并输出到nonde8out.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',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坐标(整体坐标下)FORCE=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方向坐标的向量并清零%--------------------------求总刚-------------------------------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);%调用单元刚度矩阵3a=LNODS(i,:);%临时向量,用来记录当前单元的结点编号forj=1:8%对行进行循环---按结点号循环fork=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);%调用函数生成荷载向量%-------------------------处理约束--------------------------------forj=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=HK\ALLFORCE;%方程求解，HK先求逆再与力向量左乘得到位移%-------------------------求应力---------------------------------stress=zeros(3,NELEM);%应力向量并清零fori=1:NELEM%对单元个数进行循环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(X,Y,0,0);%调用单元中心的应变矩阵endstress(:,i)=D*B*U';end%===============计算单元刚度矩阵函数===================functionEK=eKe(X,Y,YOUNG,POISS,THICK)EK=zeros(16,16);%张成16*16矩阵并清零D=(YOUNG/(1-POISS*POISS))*[1POISS0;POISS10;00(1-POISS)/2];%弹性矩阵%高斯积分采用3*3个积分点A1=5/9;A2=8/9;A3=5/9;%对应积分的加权系数A=[A1A2A3];r=(3/5)^(1/2);x=[-r0r];%积分点fori=1:3forj=1:3B=eBe(X,Y,x(i),x(j));%调用应变矩阵BJ=Jacobi(X,Y,x(i),x(j));%调用雅可比矩阵JEK=EK+A(i)*A(j)*B'*D*B*det(J)*THICK;4endend%===============计算雅可比矩阵函数===================functionJ=Jacobi(X,Y,s,t)[N_s,N_t]=DHS(s,t);x_s=0;y_s=0;x_t=0;y_t=0;forj=1:8x_s=x_s+N_s(j)*X(j);y_s=y_s+N_s(j)*Y(j);x_t=x_t+N_t(j)*X(j);y_t=y_t+N_t(j)*Y(j);endJ=[x_sy_s;x_ty_t];%===============计算应变矩阵函数===================functionB=eBe(X,Y,s,t)[N_s,N_t]=DHS(s,t);J=Jacobi(X,Y,s,t);B=zeros(3,16);fori=1:8B1=J(2,2)*N_s(i)-J(1,2)*N_t(i);B2=-J(2,1)*N_s(i)+J(1,1)*N_t(i);B(1:3,2*i-1:2*i)=[B10;0B2;B2B1];endB=B/det(J);%===============计算形函数的求导函数==================function[N_s,N_t]=DHS(s,t)N_s(1)=1/4*(1-t)*(s-t-1)+1/4*(1+s)*(1-t);N_s(2)=1/2*(1+t)*(1-t);N_s(3)=1/4*(1+t)*(s+t-1)+1/4*(1+s)*(1+t);N_s(4)=1/2*(1-s)*(1+t)-1/2*(1+s)*(1+t);N_s(5)=-1/4*(1+t)*(-s+t-1)-1/4*(1-s)*(1+t);N_s(6)=-1/2*(1+t)*(1-t);N_s(7)=-1/4*(1-t)*(-s-t-1)-1/4*(1-s)*(1-t);N_s(8)=1/2*(1-s)*(1-t)-1/2*(1+s)*(1-t);N_t(1)=-1/4*(1+s)*(s-t-1)-1/4*(1+s)*(1-t);N_t(2)=1/2*(1+s)*(1-t)-1/2*(1+s)*(1+t);N_t(3)=1/4*(1+s)*(s+t-1)+1/4*(1+s)*(1+t);N_t(4)=1/2*(1+s)*(1-s);N_t(5)=1/4*(1-s)*(-s+t-1)+1/4*(1-s)*(1+t);N_t(6)=1/2*(1-s)*(1-t)-1/2*(1-s)*(1+t);N_t(7)=-1/4*(1-s)*(-s-t-1)-1/4*(1-s)*(1-t);N_t(8)=-1/2*(1+s)*(1-s);end%===============计算总荷载矩阵函数===================functionALLFORCE=FOECEXL(NPOIN,NFORCE,FORCE)%本函数生成荷载向量ALLFORCE=zeros(2*NPOIN,1);%张成特定大小的向量，并赋值0fori=1:NFORCE5ALLFORCE((FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3);%FORCE(i,1)为作用点,FORCE(i,2:3)为x,y方向的结点力end%-------------------------输出节点位移和单元中心应力-------------------------------fori=1:NPOINfprintf(FP2,'x%d=%d\n',i,DISP(2*i-1));%输出结点x方向的位移fprintf(FP2,'y%d=%d\n',i,DISP(2*i));%输出结点y方向的位移endforj=1:NELEMfprintf(FP2,'%dx=%f\n',j,stress(1,j));%输出单元x方向的应力fprintf(FP2,'%dy=%f\n',j,stress(2,j));%输出单元y方向的应力fprintf(FP2,'%dxy=%f\n',j,stress(3,j));%输出单元切应力end%------------------------实例计算并用ANSYS进行对比结果----------------------------如图所示一个4m*1m悬臂梁，在3节点作用1*N的竖向力，参数如下：弹性模量2.0*，泊松比0.3，划分成四个单元，每个单元八个节点，单元尺寸是1m*1m。1.用MATLAB进行计算在MATLAB当前工作目录下存入初始数据文件，nonde.txt文本数据内容在表第1列，计算结果（位移和应力）在表第2、3、列，第4列为ANSYS计算的结点竖向位移234162E080.3112345678765910111213121110141516171817161519202122234040.5413.513130.5303.50x1=-2.335725e-02y1=-1.298092e-01x2=-8.552411e-05y2=-1.304313e-01x3=2.414924e-02y3=-1.314740e-01x4=2.339276e-02y4=-1.063855e-01x5=2.212796e-02y5=-8.301568e-02x6=1.768304e-05y6=-8.281498e-02x7=-2.217334e-02x16=5.504682e-07y16=-1.135171e-02x17=-1.009374e-02y17=-1.224831e-02x18=-1.428159e-02y18=-2.498899e-02x19=5.239181e-03y19=-3.600960e-03x20=0y20=0x21=0y21=0x22=0ANSYS竖向位移1-0.129512-0.130633-0.132164-0.106345-0.83270E-016-0.82963E-017-0.83027E-018-0.106849-0.61320E-0110-0.41589E-0111-0.41216E-0112-0.41644E-0162.512120.5202.501.51