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

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

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

3、scanf(FP1,'%f',1); %厚度NELEM=fscanf(FP1,'%d',1); %单元数NPOIN=fscanf(FP1,'%d',1); % 总结点数 NNODE=fscanf(FP1,'%d',1); %单元节点数NFPOIN=fscanf(FP1,'%d',1); %受力结点数NVFI*=fscanf(FP1,'%d',1); %约束结点个数LNODS=fscanf(FP1,'%f',NNODE,NELEM)' % 单元定义： 单元结点号逆时针COORD

4、=fscanf(FP1,'%f',2,NPOIN)' % 结点号 *,y坐标(整体坐标下)FPOIN=fscanf(FP1,'%f',3,NFPOIN)' % 节点力：结点号、*方向力(向右正)，Y方向力(向上正)FI*ED=fscanf(FP1,'%d',3,NVFI*)' %约束信息数组(n,3) n:受约束节点数目, (n,1):约束点号 %(n,2)与(n,3)分别为约束点*方向和y方向的约束情况,受约束为1否则为0%*%*%=平面应力问题的求解=%*%*% 刚度矩阵的生成%计算刚度矩阵，并对约束条件进展处理Ke=

5、zeros(2*NNODE,2*NNODE); % 单元刚度矩阵并清零 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(m,7),1),COORD(LNODS(m,7),2); %调用单元刚度矩阵a=

6、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 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)=FP

7、OIN(i,2); %FPION(i,2)为*方向的节点力 FORCE(b2)=FPOIN(i,3); %FPION(i,3)为y方向的节点力end%将约束信息参加总刚,总荷载 for i=1:NVFI* if FI*ED(i,2)=1 c1=2*FI*ED(i,1)-1; HK(c1,:)=0; %将一约束序号处的总刚列向量清0 HK(:,c1)=0; %将一约束序号处的总刚行向量清0 HK(c1,c1)=1; %将行列穿插处的元素置为1 FORCE(c1)=0; end if FI*ED(i,3)=1 c2=2*FI*ED(i,1); HK(c2,:)=0; HK(:,c2)=0; HK(

8、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,:); %临时向量,用来记录当前单元的节点编号 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:N

9、NODE 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(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

10、;B2i B1i/det(J); end stressm=D*BM*u' stress(:,m)=stressm;endstress %输出应力function Ke=K(E,v,h,*1,y1,*3,y3,*5,y5,*7,y7)%=单元刚度矩阵=% E 弹性模量 % v 泊松比 % h 厚度 % *1,y1,*3,y3,*5,y5,*7,y7 为4个角结点的坐标%矩阵尺寸为16 * 16Ke=zeros(16,16);D=(E/(1-v*v)*1 v 0;v 1 0;0 0 (1-v)/2;%弹性矩阵%高斯积分 采用 3 * 3 个积分点 书74页W1=5/9;W2=8/9;W3=

11、5/9; %加权系数W=W1 W2 W3;r=15(1/2)/5;*=-r 0 r;%积分点for i=1:3 for j=1:3 B=eleB(*1,y1,*3,y3,*5,y5,*7,y7,*(i),*(j); J=Jacobi(*1,y1,*3,y3,*5,y5,*7,y7,*(i),*(j); Ke=Ke+W(i)*W(j)*B'*D*B*det(J)*h; endendfunction B=eleB(*1,y1,*3,y3,*5,y5,*7,y7,s,t)%调用导函数N_s,N_t=DHS(s,t);%求Jacobi矩阵J=Jacobi(*1,y1,*3,y3,*5,y5,*

12、7,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;B2i B1i;endB=B/det(J);function J=Jacobi(*1,y1,*3,y3,*5,y5,*7,y7,s,t)%-Jacobi-%单元坐标%2,4,6,8点的坐标 *2=(*1+*3)/2;y2=(y1+y3)/2;*4=(*3+*5)/2;y4=(y3+y5)/2;*6=(*5+*7)/2;y6=

13、(y5+y7)/2;*8=(*7+*1)/2;y8=(y7+y1)/2;*=*1 *2 *3 *4 *5 *6 *7 *8;y=y1 y2 y3 y4 y5 y6 y7 y8;%调用形函数对局部坐标的导数N_s,N_t=DHS(s,t);%求Jacobi矩阵的行列式的值*_s=0;y_s=0;*_t=0;y_t=0;for i=1:8 *_s=*_s+N_s(i)*(i);y_s=y_s+N_s(i)*y(i); *_t=*_t+N_t(i)*(i);y_t=y_t+N_t(i)*y(i);endJ=*_s y_s;*_t y_t;function N=shape(s,t)%,N(1) = (

14、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=DHS(s,t)%形函数求导%,N_s(1)=-1/4*(1-t)*(-s-t-1)-1/4*(1-s)*(1-t);N_s(3

15、)=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)*(1-t);N_t(1)=-1/4*(1-s)*(-s-t-1)-1/4*(1-s)*(1-t);N_t(3)=-1/4*

16、(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);bjd.t*t 文件数据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