第7章——有限元三角形单元程序设计

湖南大学 机械与运载工程学院
Hunan University
College of Mechanical & Vehicle Engineering
崔向阳
第7讲: 三角形单元程序设计
问题描述
考虑一个平面应力问题如图所示，假设厚度h=1，材料为各项 同性，杨氏模量为 E=1，泊松比为ν=0，相关力和位移边界条 件如图中所示，问题左端为固定约束。试用两个三角形单元 分析此问题，三角形单元的网格划分如图所示。试求问题各 节点位移u、v和应力σx，σy和σxy。
A(0, (0 0)
D (4, (4 0)
C ( 4, 1)
B(0, 2)
12:45
有限单元法
崔向阳
三角形单元程序
%------------------------------------------------%initialization of matrices and vectors %------------------------------------------------ff ff=sparse(sdof,1); ( d f 1) % t %system force f vector t k=sparse(edof,edof k=sparse( edof,edof); ); %initialization of element matrix kk=sparse( kk =sparse(sdof,sdof sdof,sdof); ); %system matrix disp=sparse(sdof,1); disp =sparse(sdof,1); %system displacement vector eldisp=sparse(edof,1); eldisp =sparse(edof,1); %element displacement vector stress=zeros(nel,3); stress=zeros(nel 3); %matrix containing stress components strain=zeros(nel,3); %matrix containing strain components index=sparse(edof,1); %index vector kinmtx=sparse(3,edof); kinmtx =sparse(3,edof); %kinematic matrix matmtx=sparse(3,3); %constitutive matrix %------------------------------------------------%compute material matrices %------------------------------------------------matmtx=fematiso matmtx= fematiso(1 fematiso(1,emodule,poisson); (1 emodule poisson); %constitutive matrice (1,emodule,poisson); 12:45
-0.5*lengthy*(1+(lx+1 0.5*lengthy*(1+(lx+1-i)/lx)*(1 i)/lx)*(1-(j (j-1)/ 1)/ly ly)]; )];
有限单元法
崔向阳
三角形单元程序
%---------------------------------------------------------------------%input data for nodal connectivity for each element %nodes(i,j %nodes( i,j) ) where i->element no. and j j->connected nodes %---------------------------------------------------------------------nodes=[]; for i=1:lx for j=1:ly nodes=[nodes; (ly+1)*(i(ly+1)*(i-1)+j (ly+1)*i+j (ly+1)*i+j (ly+1)*(i(ly+1)*(i-1)+j+1;]; nodes=[nodes; (ly+1)*i+j (ly+1)*i+j (ly+1)*i+j+1 (ly+1)*(i(ly+1)*(i-1)+j+1;]; end end
A(0, 0)
D (4, 0)
F 1
C ( 4, 1)
B(0, 2)
12:45
有限单元法
崔向阳
三角形单元程序
clear all first_time= first_time =cputime cputime; ; format long g %-------------------------------------------%input data for control parameters %-------------------------------------------lengthx=4; lengthx =4; %length of xx-axis side of problem lengthy=2; %length of yy-axis side of problem emodule=1.0; emodule=1.0; ; poisson=0.0; poisson =0.0; fload= fload =-1; %elastic modulus %Poisson's ratio % the total load
12:45单元程序
%---------------------------------%input data for boundary conditions %---------------------------------bcdof=[]; bcdof =[]; bcval=[]; bcval =[]; for i=1:ly+1 bcdof=[ bcdof =[bcdof bcdof 1+2*(i 1+2*(i-1) 2+2*(i2+2*(i-1)]; bcval=[ bcval =[bcval bcval 0 0]; end
Ni,x
area=0.5*(x1*y2+x2*y3+x3*y1-x1*y3 area=0.5*(x1*y2+x2*y3+x3*y1x1*y3-x2*y1 x2*y1-x3*y2); %area of triangula area2=area*2; dhdx=(1/area2)*[(y2dhdx=(1/area2)*[(y2 -y3) (y3 (y3-y1) (y1 (y1-y2)]; dhdy=(1/area2)*[(x3 dhdy=(1/area2) dhdy =(1/area2) =(1/area2)*[(x3[(x3-x2) (x1[(x3(x1-x3) (x2(x2-x1)]; %derivatives w.r.t x %derivatives w.r.t wrt y %kinematic matrice %element stiffness matrice
B(0, 2)
A(0, 0) D (4, 0)
F 1
C ( 4, 4 1)
12:45
有限单元法
崔向阳
三角形单元程序
lx=16; % number of element in x x-axis ly=8; ly =8; % number of element in yy-axis nel=2*lx* nel =2*lx*ly ly; ; % number of element nnel=3; nnel =3; %number of nodes per element ndof=2; ndof =2; %number of dofs per node nnode=(lx+1)*(ly+1); nnode =(lx+1)*(ly+1); (lx 1) (ly 1); %total number of nodes in system sdof= sdof =nnode* nnode*ndof; ndof; %total system dofs edof= edof =nnel* nnel*ndof; ndof; %degrees of freedom per element %------------------------------------------------------%input data for nodal coordinate values %------------------------------------------------------x0=[]; for i=1:lx+1 for j j=1:ly+1 1:ly+1 x0=[x0; (i(i-1)* 1)*lengthx lengthx/lx /lx end end 12:45
12:45
有限单元法
崔向阳
三角形单元程序
%--------------------------------------%find the derivatives of shape functions %---------------------------------------
1 Ni ai bi x ci y 2A
Ni, y
kinmtx2=fekine2d( kinmtx2=fekine2d (nnel,dhdx,dhdy nnel,dhdx,dhdy); ); k=kinmtx2'*matmtx k=kinmtx2'* matmtx*kinmtx2*area; *kinmtx2*area;