main.m

应力应变计算的有限元程序。上课的时候用matlab写的一个代码

% % main
% the Finite element codes for strain and stress analysis. 
% Copyright, Yingchun Zhang, 09/13/2006

% material parameter, modulus of elasticity E, Poisson's ratio NU,
E=210e6;
NU=0.3;

% FE model 
Element=[1,1,2,4,6;
    2,1,4,3,7;
    3,6,5,7,1;
    4,6,7,8,4;
    5,1,6,4,7];
Node=[0,0,0;
    0.025,0,0;
    0,0.5,0;
    0.025,0.5,0;
    0,0,0.25;
    0.025,0,0.25;
    0,0.5,0.25;
    0.025,0.5,0.25];
ElementN=length(Element(:,1));
NodeN=length(Node(:,1));

% calculate the global stiffness matrix
K=zeros(NodeN*3,NodeN*3);
for i=1:ElementN
    x1=Node(Element(i,2),1);
    y1=Node(Element(i,2),2);
    z1=Node(Element(i,2),3);
    x2=Node(Element(i,3),1);
    y2=Node(Element(i,3),2);
    z2=Node(Element(i,3),3);
    x3=Node(Element(i,4),1);
    y3=Node(Element(i,4),2);
    z3=Node(Element(i,4),3);
    x4=Node(Element(i,5),1);
    y4=Node(Element(i,5),2);
    z4=Node(Element(i,5),3);
    
    tempK=TetrahedronElementStiffness(E,NU,x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4);
    K=TetrahedronAssemble(K,tempK,Element(i,2),Element(i,3),Element(i,4),Element(i,5));
end

% boundary condition
subK1=[K(7:12,:);K(19:24,:);K(1:6,:);K(13:18,:)];
subK2=[subK1(:,7:12),subK1(:,19:24),subK1(:,1:7),subK1(:,13:18)];
subK=subK2(1:12,1:12);

% load matrix
subF=[0;3.125;0;0;6.25;0;0;06.25;0;0;3.125;0];

% solve the linear equation system
subU=subK\subF;
U=[0;0;0;0;0;0;subU(1:6);0;0;0;0;0;0;subU(7:12)];

F=K*U;

% calculate the element stress vector
sigma=zeros(6,ElementN);
S=zeros(3,ElementN);
for i=1:ElementN
    Elist=[((Element(i,2)-1)*3+1):((Element(i,2)-1)*3+3),((Element(i,3)-1)*3+1):((Element(i,3)-1)*3+3),...
        ((Element(i,4)-1)*3+1):((Element(i,4)-1)*3+3),((Element(i,5)-1)*3+1):((Element(i,5)-1)*3+3)]';
    subU=U(Elist);
    x1=Node(Element(i,2),1);
    y1=Node(Element(i,2),2);
    z1=Node(Element(i,2),3);
    x2=Node(Element(i,3),1);
    y2=Node(Element(i,3),2);
    z2=Node(Element(i,3),3);
    x3=Node(Element(i,4),1);
    y3=Node(Element(i,4),2);
    z3=Node(Element(i,4),3);
    x4=Node(Element(i,5),1);
    y4=Node(Element(i,5),2);
    z4=Node(Element(i,5),3);
    sigma(:,i)=TetrahedronElementStresses(E,NU,x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,subU);
    S(:,i)=TetrahedronElementPStresses(sigma(:,i));
end