有限元四边形八节点程序

广州大学作业（论文）题目：平面四边形八节点有限元程序所修课程名称：有限元法及其应用任课教师姓名：张永山上课时间：星期五5-8节布置作业（论文）日期：2016年6月28日任课教师打分：2016年6月28日悬臂钢梁，尺寸如图所示；v=0.3。h=1，E=2.1e11._______土木工程______学院___15_____级__结构工程____专业姓名__杨尚荣______学号__2111516045____………………………………（密）………………………………（封）………………………………（线）………………………………1m1m1m1m1m20kN10kN悬臂钢梁1812119132214236258271020143221552416726179281112345单元划分与结点编号原始输入数据（bjd.txt)2.1E110.3152882412313201918123451422212013567152423221478916262524159101117282726160.00.00.50.01.00.01.50.02.00.02.50.03.00.03.50.04.00.04.50.05.00.00.00.51.00.52.00.53.00.54.00.55.00.50.01.00.51.01.01.01.51.02.01.02.51.0%******************************************************************%四边形八节点等参元matlab计算程序%主程序%******************************************************************%变量说明%Evh%弹性模量泊松比厚度%NPOINNELEMNVFIXNNODENFPOIN%总结点数,单元数,约束结点个数,单元节点数,受力结点数%COORDLNODS%结构节点整体坐标数组,单元定义数组,%FPOINFORCEFIXED%结点力数组，总体荷载向量,约束信息数组%HKDISPFP1%总体刚度矩阵结点位移向量数据指针文件%BDstress%单元应变矩阵单元弹性矩阵单元应力%******************************************************************clc%清屏clearall%清除内存变量formatshorte%设定输出类型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])';%结构节点坐标数组(整体坐标下)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[HK,FORCE]=H(E,v,h,NELEM,NPOIN,NFPOIN,NVFIX,LNODS,COORD,FPOIN,FIXED);%形成总纲和总荷载DISP=HK\FORCE%求位移stress=STRESS(NELEM,COORD,LNODS,NNODE,DISP,E,v)%求应力%********************************************************************%生成总刚矩阵和总荷载%********************************************************************function[HK,FORCE]=H(E,v,h,NELEM,NPOIN,NFPION,NVFIX,LNODS,COORD,FPOIN,FIXED)HK=zeros(2*NPOIN,2*NPOIN);form=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(LNODS(m,7),1),COORD(LNODS(m,7),2));%调用单元刚度矩阵（读取四个角节点）a=LNODS(m,:);%临时向量,用来记录当前单元的节点编号%对总刚度矩阵的处理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)+...Ke(j*2-1:j*2,k*2-1:k*2);end%把单刚2乘2子块加入到总纲矩阵中（其中a(j)*2，a(k)*2为求取j节点和k节点第二个位移的整体编码）endend%—————————————————————————————————%对荷载向量进行处理FORCE=zeros(2*NPOIN,1);%张成总荷载向量并清零fori=1:NFPION%对节点荷载个数进行循环b1=FPOIN(i,1)*2-1;b2=FPOIN(i,1)*2;%FPION(i,1)为作用点%荷载对应的节点号FORCE(b1)=FPOIN(i,2);%FPION(i,2)为x方向的节点力FORCE(b2)=FPOIN(i,3);%FPION(i,3)为y方向的节点力end%—————————————————————————————————%将约束信息加入总刚---划零置一（即对总纲和总荷载进行边界条件的处理）fori=1:NVFIX%对约束个数进行循环ifFIXED(i,2)==1%如果x方向存在约束c1=2*FIXED(i,1)-1;%零位移所在位置HK(c1,:)=0;%将一约束序号处的总刚列向量清0HK(:,c1)=0;%将一约束序号处的总刚行向量清0HK(c1,c1)=1;%将行列交叉处的元素置为1FORCE(c1)=0;endifFIXED(i,3)==1%如果y方向存在约束c2=2*FIXED(i,1);HK(c2,:)=0;HK(:,c2)=0;HK(c2,c2)=1;FORCE(c2)=0;endEnd%********************************************************************%求单元应力%********************************************************************functionstress=STRESS(NELEM,COORD,LNODS,NNODE,DISP,E,v)stress=zeros(3,NELEM);form=1:NELEM%对单元个数进行循环u(1:16)=0;%单元位移列向量清零;d=LNODS(m,:);%临时向量,用来记录当前单元的节点编号fori=1:NNODE%对节点个数进行循环u(i*2-1:i*2)=DISP(d(i)*2-1:d(i)*2);%从总位移向量中取出当前单元的节点位移endD=(E/(1-v*v))*[1v0;v10;00(1-v)/2];%弹性矩阵%形成应变矩阵BMBM=zeros(3,16);fori=1:NNODEJ=Jacobi(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(LNODS(m,7),1),COORD(LNODS(m,7),2),0,0);[N_s,N_t]=DHS(0,0);B1i=J(2,2)*N_s(i)-J(1,2)*N_t(i);B2i=-J(2,1)*N_s(i)+J(1,1)*N_t(i);BM(1:3,2*i-1:2*i)=[B1i0;0B2i;B2iB1i]/det(J);endstressm=D*BM*u';stress(:,m)=stressm;end%输出应力%********************************************************************%求单元刚度矩阵%********************************************************************functionKe=K(E,v,h,x1,y1,x3,y3,x5,y5,x7,y7)%E弹性模量v泊松比h厚度%x1,y1,x3,y3,x5,y5,x7,y7为4个角结点的坐标%矩阵尺寸为16乘16Ke=zeros(16,16);%即(2*NNODE,2*NNODE)D=(E/(1-v*v))*[1v0;v10;00(1-v)/2];%弹性矩阵%高斯积分采用3x3个积分点W1=5/9;W2=8/9;W3=5/9;%加权系数W=[W1W2W3];r=15^(1/2)/5;x=[-r0r];%积分点fori=1:3forj=1:3B=eleB(x1,y1,x3,y3,x5,y5,x7,y7,x(i),x(j));J=Jacobi(x1,y1,x3,y3,x5,y5,x7,y7,x(i),x(j));Ke=Ke+W(i)*W(j)*B'*D*B*det(J)*h;end%求解单元刚度矩阵（根据虚功原理）End%********************************************************************%求雅克比矩阵%********************************************************************functionJ=Jacobi(x1,y1,x3,y3,x5,y5,x7,y7,s,t)%单元坐标%用1,3,5,7点的坐标表示2,4,6,8点的坐标x2=(x1+x3)/2;y2=(y1+y3)/2;x4=(x3+x5)/2;y4=(y3+y5)/2;x6=(x5+x7)/2;y6=(y5+y7)/2;x8=(x7+x1)/2;y8=(y7+y1)/2;x=[x1x2x3x4x5x6x7x8];y=[y1y2y3y4y5y6y7y8];%调用形函数对局部坐标的导数[N_s,N_t]=DHS(s,t);%求Jacobi矩阵的行列式的值x_s=0;y_s=0;x_t=0;y_t=0;fori=1:8x_s=x_s+N_s(i)*x(i);y_s=y_s+N_s(i)*y(i);x_t=x_t+N_t(i)*x(i);y_t=y_t+N_t(i)*y(i);endJ=[x_sy_s;x_ty_t];%*************************************************