Contents

FEM for cantilever beam under constant loading

clear

% data
E=200e+3; I=1e-4; q=200e-3; L=2;
% discretization
Nel = 100; Nodes = Nel+1; node_dof=2; Ndofs = node_dof*Nodes; h=L/Nel;
kinematic_bc=[1 0; 2 0];
coord=0:h:L;
Edof=[[1:Nel]' [1:2:Ndofs-3]' [2:2:Ndofs-2]' [3:2:Ndofs-1]' [4:2:Ndofs]'];
% assembling
KG=zeros(Ndofs,Ndofs);
f=zeros(Ndofs,1);
    Kel = 2*E*I /h^3 * ...
         [6   3*h   -6   3*h;
          3*h 2*h^2 -3*h h^2;
          -6 -3*h    6  -3*h;
          3*h h^2 -3*h 2*h^2]; % element stiffness matrix
for iel=1:Nel
    Qel = [q*h/2; q*h^2/12; q*h/2; -q*h^2/12]; % element load vector
    [KG,f]=assem(Edof(iel,:),KG,Kel,f,Qel);
end
% solution
[u,r]=solveq(KG,f,kinematic_bc);
figure(1);
plot(coord,u(1:2:Ndofs),'b*'), axis ij, hold on % FEM approximation
w_exact=q*L^4/8/E/I; % free end deflection for cantilever beam
plot(L,w_exact,'ro')

% nodal forces
for iel=1:Nel
    uel=u((iel-1)*2+1:iel*2+2);
    fel(1:4,iel)=Kel*uel-Qel;
end

Hermite shape functions

figure(2)
h=2;
x=linspace(0.,h,100);
xi=x/h;
y=[1-3*(xi.^2)+2*(xi.^3);
  (xi-2*(xi.^2)+(xi.^3))*h;
  3*(xi.^2)-2*(xi.^3);
  (-xi.^2+xi.^3)*h];
plot(x,y)
grid

Published with MATLAB® 7.6