四边形八节点等参元matlab程序.docx

1、四边形八节点等参元matlab程序悬臂钢梁，尺寸如图一所示；v=0.3。h=1，E=2.1e11. 图一 悬臂钢梁图二 单元划分与结点编号Matlab 输出结果 附录：有限元ANSYS分析结果采用PLANE183单元（四边形八节点）单元得出的结构Y向最大位移为-0.216E-04。约等于matlab平面四边形八节点等参元结点Y向最大位移-2.4024E-5。附录：%-四边形八节点等参元 matlab计算程序-% 主 程 序 %*%*% 2012年 % 本程序只能处理集中荷载作用下的情况% 只输出了节点位移、单元中心点的应力%*%*% 变量说明% E v h% 弹性模量 泊松比 厚度% NPOI

2、N NELEM NVFIX NNODE NFPOIN % 总结点数 , 单元数, 约束结点个数, 单元节点数 ,受力结点数% COORD LNODS % 结构节点整体坐标数组, 单元定义数组, % FPOIN FORCE FIXED% 结点力数组， 总体荷载向量, 约束信息数组% HK DISP% 总体刚度矩阵,结点位移向量%*clear allformat short e FP1=fopen(bjd.txt,rt); %打开数据文件%读入控制数据E=fscanf(FP1,%f,1); %弹性模量 v=fscanf(FP1,%f,1); % 泊松比h=fscanf(FP1,%f,1); %厚度

3、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); % 结点号 x,y坐标(整体坐标下)FPOIN=fscanf(FP1,%f,3,NFPOIN); % 节点力：结点号、X方向力(向右正

4、)，Y方向力(向上正)FIXED=fscanf(FP1,%d,3,NVFIX); %约束信息数组(n,3) n:受约束节点数目, (n,1):约束点号 %(n,2)与(n,3)分别为约束点x方向和y方向的约束情况,受约束为1否则为0%*%*%=平面应力问题的求解=%*%*% % 刚度矩阵的生成%计算刚度矩阵，并对约束条件进行处理Ke=zeros(2*NNODE,2*NNODE); % 单元刚度矩阵并清零 HK=zeros(2*NPOIN,2*NPOIN); % 张成总刚矩阵并清零%调用子程序 生成单元刚度矩阵for m=1:NELEM %m为单元号 Ke=K(E,v,h,. COORD(LNO

5、DS(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,:); %临时向量,用来记录当前单元的节点编号%对总刚度矩阵的处理 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)

6、:a(k)*2)+. Ke(j*2-1:j*2,k*2-1:k*2); end endend%对荷载向量进行处理FORCE=zeros(2*NPOIN,1); % 张成总荷载向量并清零for i=1:NFPOIN 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 FIXED(i,2)=1 c1=2*FIXED(

7、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 c2=2*FIXED(i,1); HK(c2,:)=0; HK(:,c2)=0; HK(c2,c2)=1; FORCE(c2)=0; end end%=%=DISP=HKFORCE %计算节点位移向量%=%=%求解单元应力stress=zeros(3,NELEM);for m=1:NELEM u(1:16)=0; d=LNODS(m,:); %临

8、时向量,用来记录当前单元的节点编号 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=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(

9、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/det(J); end stressm=D*BM*u; stress(:,m)=stressm;endstress %输出应力function Ke=K(E,v,h,x1,y1,x3,y3,x5,y5,x7,y7)%=单元刚度矩阵=% E 弹

10、性模量 % v 泊松比 % h 厚度 % x1,y1,x3,y3,x5,y5,x7,y7 为4个角结点的坐标%矩阵尺寸为16 x 16Ke=zeros(16,16);D=(E/(1-v*v)*1 v 0;v 1 0;0 0 (1-v)/2;%弹性矩阵%高斯积分 采用 3 x 3 个积分点 书74页W1=5/9;W2=8/9;W3=5/9; %加权系数W=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

11、,x7,y7,x(i),x(j); Ke=Ke+W(i)*W(j)*B*D*B*det(J)*h; endendfunction B=eleB(x1,y1,x3,y3,x5,y5,x7,y7,s,t)%调用导函数N_s,N_t=DHS(s,t);%求Jacobi矩阵J=Jacobi(x1,y1,x3,y3,x5,y5,x7,y7,s,t);%求应变矩阵BB=zeros(3,16);for 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;B2

12、i B1i;endB=B/det(J);function J=Jacobi(x1,y1,x3,y3,x5,y5,x7,y7,s,t)%-Jacobi-%单元坐标%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=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=

13、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)*y(i);endJ=x_s y_s;x_t y_t;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;

14、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=DHS(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)*

15、(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)*(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);

16、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); bjd.txt 文件数据2.1E11 0.3 1 5 28 8 1 31 2 3 13 20 19 18 123 4 5 14 22 21 20 135 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.03.0 1.03.5 1.04.0 1.04.5 1.05.0 1.017 0 -100001 1 112 1 118 1 1