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

三角形单元程序
%------------------------------------------------%compute element matrices and vectors and assemble %------------------------------------------------for iel=1:nel %loop for the total number of element
nd(1)=nodes(iel,1); nd(2)=nodes(iel,2); nd(3)=nodes(iel,3);
%1st connected node for (iel)-th element %2nd connected node for (iel)-th element %3rd connected node for (iel)-th element
%coord values of 1st node %coord values of 2nd node %coord values of 3rd node %extract system dofs for the element
x1=x0(nd(1),1); y1=x0(nd(1),2); x2=x0(nd(2),1); y2=x0(nd(2),2); x3=x0(nd(3),1); y3=x0(nd(3),2); index=feeldof(nd,nnel,ndof);
Ni,x
area=0.5*(x1*y2+x2*y3+x3*y1-x1*y3-x2*y1-x3*y2); %area of triangula area2=area*2; dhdx=(1/area2)*[(y2-y3) (y3-y1) (y1-y2)]; dhdy=(1/area2)*[(x3-x2) (x1-x3) (x2-x1)]; %derivatives w.r.t x %derivatives w.r.t y %kinematic matrice %element stiffness matrice
三角形单元程序设计
Байду номын сангаас
问题描述
考虑一个平面应力问题如图所示，假设厚度h=1，材料为各项 同性，杨氏模量为 E=1，泊松比为 ν=0，相关力和位移边界条 件如图中所示，问题左端为固定约束。试用两个三角形单元 分析此问题，三角形单元的网格划分如图所示。试求问题各 节点位移u、v和应力σx，σy和σxy。
1 T e de K de 2
16:51
-0.5*lengthy*(1+(lx+1-i)/lx)*(1-(j-1)/ly)];
三角形单元程序
%---------------------------------------------------------------------%input data for nodal connectivity for each element %nodes(i,j) where i->element no. and j->connected nodes %---------------------------------------------------------------------nodes=[]; for i=1:lx for j=1:ly nodes=[nodes; (ly+1)*(i-1)+j (ly+1)*i+j (ly+1)*(i-1)+j+1;]; nodes=[nodes; (ly+1)*i+j (ly+1)*i+j+1 (ly+1)*(i-1)+j+1;]; end end
16:51
三角形单元程序
%---------------------------------%input data for boundary conditions %---------------------------------bcdof=[]; bcval=[]; for i=1:ly+1 bcdof=[bcdof 1+2*(i-1) 2+2*(i-1)]; bcval=[bcval 0 0]; end
A(0, 0) D(4, 0)
F 1
fload=-1;
% the total load
B(0, 2)
C (4, 1)
16:51
三角形单元程序
lx=16; % number of element in x-axis ly=8; % number of element in y-axis nel=2*lx*ly; % number of element nnel=3; %number of nodes per element ndof=2; %number of dofs per node nnode=(lx+1)*(ly+1); %total number of nodes in system sdof=nnode*ndof; %total system dofs edof=nnel*ndof; %degrees of freedom per element %------------------------------------------------------%input data for nodal coordinate values %------------------------------------------------------x0=[]; for i=1:lx+1 for j=1:ly+1 x0=[x0; (i-1)*lengthx/lx end end 16:51
Ni, y
kinmtx2=fekine2d(nnel,dhdx,dhdy); k=kinmtx2'*matmtx*kinmtx2*area;
K e BT DBhdxdy BT DBA
A
kk=feasmb_2(kk,k,index); end
%assemble element matrics
16:51
三角形单元程序
%--------------------------------------%find the derivatives of shape functions %---------------------------------------
1 Ni ai bi x ci y 2A
A(0, 0) D(4, 0)
C (4, 1)
B(0, 2)
16:51
三角形单元程序
%------------------------------------------------%initialization of matrices and vectors %------------------------------------------------ff=sparse(sdof,1); %system force vector k=sparse(edof,edof); %initialization of element matrix kk=sparse(sdof,sdof); %system matrix disp=sparse(sdof,1); %system displacement vector eldisp=sparse(edof,1); %element displacement vector 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); %kinematic matrix matmtx=sparse(3,3); %constitutive matrix %------------------------------------------------%compute material matrices %------------------------------------------------matmtx=fematiso(1,emodule,poisson); %constitutive matrice 16:51
A(0, 0) D(4, 0)
F 1
C (4, 1)
B(0, 2)
16:51
三角形单元程序
clear all first_time=cputime; format long %-------------------------------------------%input data for control parameters %-------------------------------------------lengthx=4; %length of x-axis side of problem lengthy=2; %length of y-axis side of problem emodule=1.0; poisson=0.0; %elastic modulus %Poisson's ratio
16:51
三角形单元程序
kk1=kk; %-------------------------------------------------------------------------% force vector %-------------------------------------------------------------------------ff(sdof,1)=fload; %-----------------------% apply boundary condition %-----------------------[kk,ff]=feaplyc2(kk,ff,bcdof,bcval); %------------------------% solve the matrix equation %------------------------%disp=kk\ff; [LL UU]=lu(kk); utemp=LL\ff; disp=UU\utemp;