allkindofbeam.m

A finite element simply supported beam model for vibration analysis

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%
% A finite element simply supported beam model for 
% vibration analysis
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clc
clear
%
%   Input data
%
% default data!!
L = 1.00;   % beam span(meter)
b = 0.05;   % beam width(meter)
h = 0.01;   % beam thickness(meter)
E = 69.e+9; % young's modulus(Pascal)
Iz = b*h^3/12; % second area moment  
EI = E*Iz;  % bending rigidity
V=b*h*L; 
rho = 2700; % density (kg/meter^3)
m_total = rho*V;

Wmax = 2000;% Maximum frequency Wmax
I = 2; % Station # of Excitation 
J = 5; % Station # of Response
disp('  ');
disp('Welcome to vibration analysis of a beam by the finite element method');
disp('Please provide the beam model data');
disp(' ');
nelem =input('number of elements to model the beam = ');
disp(' ')
disp('Please select the bounday conditions');
disp('free-free beam:   1');
disp('fixed-free beam:  2');
disp('fixed-fixed beam: 3');
disp('fixed-simply supported beam: 4');
disp('simple-simple beam: 5');
bccase =input('boundary condition for this analysis  = ');
disp(' ');
n_mode =input('specify number of modes to be computed  =  ');

%
% begin computations
%
node_g = nelem+1; % total number of global nodes
node_e = 2;   % number of nodes in each element
dof_node  = 2;   % degrees of freedom per node
dof_elem = dof_node*node_e;  % total elemental dofs
dof_g = dof_node*node_g;  % total global dofs
%
%  element stiffness matrix
%
ell = L/nelem;  % elemental beam length
kr = EI/ell^3;
ke=kr*[
    12   6*ell         -12    6*ell;
   6*ell 4*ell^2       -6*ell 2*ell^2; 
   -12  -6*ell          12   -6*ell;
   6*ell 2*ell^2       -6*ell 4*ell^2];

mass=rho*b*h*ell;   
% lumped element mass matrix
me = sparse(zeros(dof_elem,dof_elem));
me(1,1)  = 0.5*mass;
me(3,3)  = 0.5*mass;
me(2,2)  = 0.5*mass*ell^2/12; 
me(4,4)  = 0.5*mass*ell^2/12;
%consistent mass case
me = mass*[
         13/35         11*ell/210        9/70          -13*ell/420;
         11*ell/210    ell^2/105        13*ell/420    -ell^2/140;
         9/70          13*ell/420        13/35         -11*ell/210;
         -13*ell/420   -ell^2/140        -11*ell/210   ell^2/105];	  
%
% generate connectivity table
%
conn=sparse(zeros(nelem, node_e));  % connectivity table 
%
IC=0;
for j = 1: node_g; 
IC = IC+1;
conn(IC, 1) = j;
conn(IC, 2) = j+1;
end;
% 
% generate assembly operator 
% 
Ls = sparse(zeros(nelem*dof_elem, dof_g));
i2 = speye(2);
for j = 1:nelem,
  for k = 1:node_e,
    r = (j-1)*dof_elem + (k-1)*dof_node; r = (r+1:r+dof_node);
    c = (conn(j,k)-1)*dof_node; c = (c+1:c+dof_node);
    Ls(r,c) = i2;
  end
end
%
% generate  the nodal coordinates
%
coord= zeros(node_g,1);
nodes = 0;
for i = 1: node_g  % x-direction sweep
  x =(i-1)*ell;
nodes = nodes+1;
coord(nodes, 1) = x;
end;

% obtain the unassembled stiffness matrices

K_ell = sparse(zeros(nelem*dof_elem, nelem*dof_elem));
M_ell = sparse(zeros(nelem*dof_elem, nelem*dof_elem));

for ne =1:nelem,
	  ig = (ne-1)*dof_elem +1: ne*dof_elem;
	  jg = (ne-1)*dof_elem +1: ne*dof_elem;
      K_ell(ig, jg) = K_ell(ig, jg) + ke;
	  M_ell(ig, jg) = M_ell(ig, jg) + me;
end;

% assemble global mass and stiffness matrices

Kg = Ls' * K_ell * Ls;
Mg = Ls' * M_ell * Ls;
%
%  apply boundary conditions.
%
ng = size(Kg,1);
if     (bccase==1)  % free free case
	wmodes_dof =1:2:(ng-1);
	BeamModel = 'Free - Free Beam';
	% do nothing
elseif (bccase==2)	% fixed - free case
	Kg([1,2],:) = []; Kg(:,[1,2]) = []; 
	Mg([1,2],:) = []; Mg(:,[1,2]) = []; 
	nn=size(Kg,1);
	wmodes_dof =1:2:(nn-1);
	BeamModel = 'Fixed - Free Beam';
elseif (bccase==3)	% fixed - fixed case
	Kg([1,2,(ng-1),ng],:) = []; Kg(:,[1,2,(ng-1),ng]) = []; 
	Mg([1,2,(ng-1),ng],:) = []; Mg(:,[1,2,(ng-1),ng]) = []; 
	nn=size(Kg,1);
	wmodes_dof =1:2:(nn-1);
	BeamModel = 'Fixed - Fixed Beam';
elseif (bccase==4)	% fixed - simple support case
	Kg([1,2,(ng-1)],:) = []; Kg(:,[1,2,(ng-1)]) = []; 
	Mg([1,2,(ng-1)],:) = []; Mg(:,[1,2,(ng-1)]) = []; 
	nn=size(Kg,1);
	wmodes_dof =1:2:(nn-1);
	BeamModel = 'Fixed - Simple Support Beam';
elseif (bccase==5)	% simple - simple support case
	Kg([1,(ng-1)],:) = []; Kg(:,[1,(ng-1)]) = []; 
	Mg([1,(ng-1)],:) = []; Mg(:,[1,(ng-1)]) = []; 
	nn=size(Kg,1);
	wmodes_dof =2:2:(nn-1);
	BeamModel = 'Simple - Simple Support Beam';
elseif (bccase>5)
    disp('your boundary condition is not in the menu. program will terminate');
    exit;	
end;

% set damping matrix
Cg = 0.1*Kg;

% set gobal matrix
Kg=full(Kg);
Mg=full(Mg);
Cg=full(Cg);