有限元四边形八节点程序

1、精选优质文档-倾情为你奉上_土木工程_学院_15_级_结构工程_专业 姓名_杨尚荣_ 学号_（密）（封）（线）广州大学作业（论文）题目： 平面四边形八节点有限元程序 所修课程名称： 有限元法及其应用 任课教师姓名： 张永山 上课时间： 星期五 5-8节 布置作业（论文）日期：2016年 6 月 28 日任课教师打分： 2016年 6 月 28 日悬臂钢梁，尺寸如图所示；v=0.3。h=1，E=2.1e11. 悬臂钢梁 单元划分与结点编号原始输入数据（bjd.txt)2.1E11 0.3 1 5 28 8 2 41 2 3 13 20 19 18 123 4 5 14 22 21 20 135

2、6 7 15 24 23 22 147 8 9 16 26 25 24 159 10 11 17 28 27 26 160.0 0.00.5 0.01.0 0.01.5 0.02.0 0.02.5 0.03.0 0.03.5 0.04.0 0.04.5 0.05.0 0.00.0 0.51.0 0.52.0 0.53.0 0.54.0 0.55.0 0.50.0 1.00.5 1.01.0 1.01.5 1.02.0 1.02.5 1.0%*% 四边形八节点等参元 matlab计算程序% 主 程 序 %* % 变量说明% E v h% 弹性模量 泊松比 厚度% NPOIN NELEM NVFI

3、X NNODE NFPOIN % 总结点数 , 单元数, 约束结点个数, 单元节点数 ,受力结点数% COORD LNODS % 结构节点整体坐标数组, 单元定义数组, % FPOIN FORCE FIXED% 结点力数组， 总体荷载向量, 约束信息数组% HK DISP FP1% 总体刚度矩阵 结点位移向量 数据指针文件% B D stress% 单元应变矩阵 单元弹性矩阵 单元应力%*clc %清屏clear all %清除内存变量 format short e %设定输出类型FP1=fopen(bjd.txt,rt); %打开数据文件%读入控制数据E=fscanf(FP1,%f,1);

4、%弹性模量 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); %结构节点坐标数组(整体坐标下

5、)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否则为0HK,FORCE=H(E,v,h,NELEM,NPOIN,NFPOIN,NVFIX,LNODS,COORD,FPOIN,FIXED);%形成总纲和总荷载DISP=HKFORCE %求位移stress=STRESS(NELEM,COORD,LNODS,NNODE,

6、DISP,E,v)%求应力%*%生成总刚矩阵和总荷载%* functionHK,FORCE=H(E,v,h,NELEM,NPOIN,NFPION,NVFIX,LNODS,COORD,FPOIN,FIXED)HK=zeros(2*NPOIN,2*NPOIN);for m=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

7、(m,7),1),COORD(LNODS(m,7),2); %调用单元刚度矩阵（读取四个角节点）a=LNODS(m,:); %临时向量,用来记录当前单元的节点编号%对总刚度矩阵的处理 for j=1:8 %对行进行循环（按照节点号循环） for 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)+. Ke(j*2-1:j*2,k*2-1:k*2); end %把单刚2乘2子块加入到总纲矩阵中（其中a(j)*2，a(k)*2为求取j节点和k节点第二

8、个位移的整体编码） endend%对荷载向量进行处理FORCE=zeros(2*NPOIN,1); % 张成总荷载向量并清零for i=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%将约束信息加入总刚-划零置一（即对总纲和总荷载进行边界条件的处理） for i=1:NVFIX %对约束个数进行循环 if

9、 FIXED(i,2)=1 %如果x方向存在约束 c1=2*FIXED(i,1)-1; %零位移所在位置 HK(c1,:)=0; %将一约束序号处的总刚列向量清0 HK(:,c1)=0; %将一约束序号处的总刚行向量清0 HK(c1,c1)=1; %将行列交叉处的元素置为1 FORCE(c1)=0; end if FIXED(i,3)=1 %如果y方向存在约束 c2=2*FIXED(i,1); HK(c2,:)=0; HK(:,c2)=0; HK(c2,c2)=1; FORCE(c2)=0; end End%*%求单元应力%* function stress=STRESS(NELEM,COOR

10、D,LNODS,NNODE,DISP,E,v)stress=zeros(3,NELEM);for m=1:NELEM %对单元个数进行循环 u(1:16)=0; %单元位移列向量清零; d=LNODS(m,:); %临时向量,用来记录当前单元的节点编号 for i=1:NNODE %对节点个数进行循环 u(i*2-1:i*2)=DISP(d(i)*2-1:d(i)*2); %从总位移向量中取出当前单元的节点位移 end D=(E/(1-v*v)*1 v 0;v 1 0;0 0 (1-v)/2;%弹性矩阵 %形成应变矩阵BM BM=zeros(3,16); for i=1:NNODE J=Jac

11、obi(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)=B1i 0;0 B2i;B2i B1i/d

12、et(J); end stressm=D*BM*u; stress(:,m)=stressm;end %输出应力%*%求单元刚度矩阵%* function Ke=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)*1 v 0;v 1 0;0 0 (1-v)/2;%弹性矩阵%高斯积分 采用3x3个积分点W1=5/9;W2=8/9;W3=5/9; %加权系数W

13、=W1 W2 W3;r=15(1/2)/5;x=-r 0 r; %积分点 for i=1:3 for j=1:3 B=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%*%求雅克比矩阵%* function J=Jacobi(x1,y1,x3,y3,x5,y5,x7,y7,s,t)%单元坐标%用1,3,5,7点的坐标表示2,4,6,8点的坐标 x2=(x1+x3)/

14、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=x1 x2 x3 x4 x5 x6 x7 x8;y=y1 y2 y3 y4 y5 y6 y7 y8;%调用形函数对局部坐标的导数N_s,N_t=DHS(s,t);%求Jacobi矩阵的行列式的值x_s=0;y_s=0;x_t=0;y_t=0;for i=1:8 x_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

15、)*y(i);endJ=x_s y_s;x_t y_t;%*%求应变矩阵B%* function B=eleB(x1,y1,x3,y3,x5,y5,x7,y7,s,t)N_s,N_t=DHS(s,t); %调用导函数J=Jacobi(x1,y1,x3,y3,x5,y5,x7,y7,s,t); %求Jacobi矩阵B=zeros(3,16); %求应变矩阵Bfor i=1:8 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); B(1:3,2*i-1:2*i)=B1i 0;0 B2i;B2i B1i;endB=B/

16、det(J);%*%设置形函数%* function N=shape(s,t)%,N(1) = (1-s)*(1-t)*(-s-t-1)/4;N(3) = (1+s)*(1-t)*(s-t-1)/4;N(5) = (1+s)*(1+t)*(s+t-1)/4;N(7) = (1-s)*(1+t)*(-s+t-1)/4;N(2) = (1-t)*(1+s)*(1-s)/2;N(4) = (1+s)*(1+t)*(1-t)/2;N(6) = (1+t)*(1+s)*(1-s)/2;N(8) = (1-s)*(1+t)*(1-t)/2; %*%型函数的导函数%* function N_s,N_t=DH

17、S(s,t)%形函数求导%,N_s(1)=-1/4*(1-t)*(-s-t-1)-1/4*(1-s)*(1-t);N_s(3)=1/4*(1-t)*(s-t-1)+1/4*(1+s)*(1-t);N_s(5)=1/4*(1+t)*(s+t-1)+1/4*(1+s)*(1+t);N_s(7)=-1/4*(1+t)*(-s+t-1)-1/4*(1-s)*(1+t);N_s(2)=1/2*(1-s)*(1-t)-1/2*(1+s)*(1-t);N_s(4)=1/2*(1+t)*(1-t);N_s(6)=1/2*(1-s)*(1+t)-1/2*(1+s)*(1+t);N_s(8)=-1/2*(1+t)

18、*(1-t);N_t(1)=-1/4*(1-s)*(-s-t-1)-1/4*(1-s)*(1-t);N_t(3)=-1/4*(1+s)*(s-t-1)-1/4*(1+s)*(1-t);N_t(5)=1/4*(1+s)*(s+t-1)+1/4*(1+s)*(1+t);N_t(7)=1/4*(1-s)*(-s+t-1)+1/4*(1-s)*(1+t);N_t(2)=-1/2*(1+s)*(1-s);N_t(4)=1/2*(1+s)*(1-t)-1/2*(1+s)*(1+t);N_t(6)=1/2*(1+s)*(1-s);N_t(8)=1/2*(1-s)*(1-t)-1/2*(1-s)*(1+t);结果输出DISP = 0 0 -8.5685e-10 -2.2935e-08 1.7374e-08 -6.6907e-08 4.7180e-08 -1.2425e-07 9.5082e-08 -1.1254e-07 9.6672e-08 -7.9875e-08 0 0 -2.3101e-07 -4.5306e