FEATool  v1.7.1
Finite Element Analysis Toolbox
 All Files Functions Pages
ex_euler_beam2.m File Reference

Description

EX_EULER_BEAM2 1D Euler-Bernoulli beam vibration example.

[ FEA, OUT ] = EX_EULER_BEAM2( VARARGIN ) 1D Euler-Bernoulli beam model example for. calculating vibration modes and eigenfrequencies. Accepts the following property/value pairs.

Input       Value/{Default}        Description
-----------------------------------------------------------------------------------
L           scalar {1}             Beam length
E           scalar {1}             Elastic modulus
I           expression {1}         Cross section moment of intertia
rho         expression {1}         Material density
A           expression {1}         Cross section area
nx          scalar {100}           Number of grid cells
iplot       scalar 0/{1}           Plot solution (=1)
                                                                                  .
Output      Value/(Size)           Description
-----------------------------------------------------------------------------------
fea         struct                 Problem definition struct
out         struct                 Output stuct

Code listing

 cOptDef = { 'L',        1;
             'E',        1;
             'I',        1;
             'rho',      1;
             'A',        1;
             'nx'        100;
             'iplot',    1;
             'tol',      1e-3;
             'fid',      1 };
 [got,opt] = parseopt( cOptDef, varargin{:} );
 fid       = opt.fid;


% Grid generation.
 fea.sdim = {'x'};
 fea.grid = linegrid( opt.nx, 0, opt.L );


% Problem and equation definitions.
% Problem and equation definitions.
 fea = addphys( fea, @eulerbeam );
 fea.phys.eb.eqn.coef{1,end}  = { opt.rho };
 fea.phys.eb.eqn.coef{2,end}  = { opt.A };
 fea.phys.eb.eqn.coef{3,end}  = { opt.E };
 fea.phys.eb.eqn.coef{4,end}  = { opt.I };
 fea.phys.eb.bdr.coef{1,5}{2} = 0;
 fea = parsephys( fea );


% Coefficients and equation/postprocessing expressions.
 fea.expr = { 'L',    opt.L ;
              'E',    opt.E ;
              'I',    opt.I ;
              'rho',  opt.rho ;
              'A',    opt.A };


% Parse and assemble matrices.
 fea = parseprob( fea );
 [M,A] = assembleprob( fea, 'f_m', 1, 'f_a', 1, 'f_sparse', 1, 'icub', 4 );


% Eliminate Dirichlet boundary conditions from matrices.
 ix_b = ones(size(M,1),1);
 [tmp,ix_b] = bdrsetd( [], ix_b, fea );
 ix_rem = find(ix_b==0);
 M(:,ix_rem) = [];
 A(:,ix_rem) = [];
 M = M';
 A = A';
 M(:,ix_rem) = [];
 A(:,ix_rem) = [];
 M = M';
 A = A';


% Solve for eigenvalues and eigenvectors.
 n_eigs = 6;
 [V,D] = eigs( A, M, n_eigs, 0 );
 efq = sqrt(diag(D))/(2*pi);


% Postprocessing.
 if( opt.iplot )
   figure, hold on, grid on
   cols = {'b' 'r' 'y' 'g' 'c' 'm'};
   x = fea.grid.p;
   ind_dof_v = find(ix_b(1:length(x)));
   for i=1:4
     y = zeros(size(x));
     y(ind_dof_v) = V(1:length(ind_dof_v),n_eigs-i+1);
     plot( x, y, 'color', cols{mod(i-1,6)+1}, 'linewidth', 2 )
   end
   legend( {'Mode 1' 'Mode 2' 'Mode 3' 'Mode 4'} )
   title( 'Vibration modes' )
   xlabel( 'x' )
   ylabel( 'Displacement' )
 end

% Error checking.
 efq_ref = [0.5596 3.5069 9.8194 19.2421 31.8086 47.5166]';
 err = abs(sort(efq(1:6))-efq_ref)./efq_ref;


 out.efq  = sort(efq);
 out.err  = err;
 out.pass = all(out.err<opt.tol);
 if ( nargout==0 )
   clear fea out
 end