Vibration of Supported Cantilever Beam

[ Problem Description ] [ Supported Cantilever Beam ] [ Input File ] [ Output File ]


PROBLEM DESCRIPTION

The purpose of this numerical example is to demonstrate the use of:

  1. Element level stiffness and mass matrices, as would be found in finite element analysis.
  2. Destination arrays as a means of mapping degrees of freedom in the mass and stiffness finite element matrices onto the global stiffness matrix (and mass matrix) degrees of freedom.
  3. The eigenvalue functions to compute the natural periods, and model shapes of vibration.

We have deliberately selected a problem with a small number of unknowns so that step-by-step details of items [1] and [2] may be written down, and solved using functions from the matrix library alone.


SUPPORTED CANTILEVER BEAM

Figure 1 shows uniform cantilever beam with a supported end point.

Figure 1 : Cantilever Beam with End Support

The cantilever has a length of 20 m, modulus of elasticity E=200 GPa, and a constant moment of inertia along its length of I = 15.5 mm^4. We will assume that the cantilever is rigid in its axial direction, and that all deformations are small. The boundary conditions are full-fixity at the base and zero translational displacement at the cantilever end.

Overall behavior of the cantilever will be modeled with five beam finite elements, and nine global degrees of freedom.They are a translation and rotation at each of the internal nodes, and a single rotational degree of freedom at the hinged cantilever support. Numbering for the global degrees of freedom is shown on the left-hand side of figure 1. Lateral displacements along each beam element will be expressed as the superposition of four cubic interpolation functions, each weighted by a nodal displacement factor. The four nodal displacements are two beam-end lateral displacements, and two beam-end rotations, as shown on the right-hand side of the Figure 1.


INPUT FILE

A finite element model for the supported cantilever is defined in parts [a] to [f] of the following input file. The eigenvalue problem is solved in part [g].

/* [a] : Define section/material properties */

   E    = 200000   MPa;
   I    = 15.5E+6 mm^4;
   L    = 4       m;
   mbar = 31.6 kg/m;

/* [b] : Define (4x4) stiffness matrix for beam element */

   stiff = Matrix([4,4]);
   stiff = ColumnUnits( stiff, [N,  N/m,  N,  N/m] );
   stiff =    RowUnits( stiff, [m], [1] );
   stiff =    RowUnits( stiff, [m], [3] );
   stiff[1][1] =  4*E*I/L;
   stiff[1][2] =  6*E*I/(L^2);
   stiff[1][3] =  2*E*I/L;
   stiff[1][4] = -6*E*I/(L^2);

   stiff[2][1] =   6*E*I/(L^2);
   stiff[2][2] =  12*E*I/(L^3);
   stiff[2][3] =   6*E*I/(L^2);
   stiff[2][4] = -12*E*I/(L^3);

   stiff[3][1] =  2*E*I/L;
   stiff[3][2] =  6*E*I/(L^2);
   stiff[3][3] =  4*E*I/L;
   stiff[3][4] = -6*E*I/(L^2);

   stiff[4][1] =  -6*E*I/(L^2);
   stiff[4][2] = -12*E*I/(L^3);
   stiff[4][3] =  -6*E*I/(L^2);  
   stiff[4][4] =  12*E*I/(L^3);

   PrintMatrix(stiff);

/* [c] : Define (4x4) consistent mass matrix for beam element */

   mass = Matrix([4,4]);
   mass = ColumnUnits( mass, [ kg*m, kg, kg*m, kg] );
   mass =    RowUnits( mass, [m], [1] );
   mass =    RowUnits( mass, [m], [3] ); 

   mass[1][1] = (mbar*L/420)*4*L*L;
   mass[2][1] = (mbar*L/420)*22*L;
   mass[3][1] = (mbar*L/420)*-3*L^2;
   mass[4][1] = (mbar*L/420)*13*L;

   mass[1][2] = (mbar*L/420)*22*L;
   mass[2][2] = (mbar*L/420)*156;
   mass[3][2] = (mbar*L/420)*-13*L;
   mass[4][2] = (mbar*L/420)*54;

   mass[1][3] = (mbar*L/420)*-3*L^2;
   mass[2][3] = (mbar*L/420)*-13*L;
   mass[3][3] = (mbar*L/420)*4*L^2;
   mass[4][3] = (mbar*L/420)*-22*L  ;

   mass[1][4] = (mbar*L/420)*13*L;
   mass[2][4] = (mbar*L/420)*54;
   mass[3][4] = (mbar*L/420)*-22*L;   mass[4][4] = (mbar*L/420)*156;

   PrintMatrix(mass);

/* [d] : Destination Array beam element connectivity */

   LD = [ 0, 0, 1, 2  ;
          1, 2, 3, 4  ; 
          3, 4, 5, 6  ; 
          5, 6, 7, 8  ; 
          7, 8, 9, 0 ]; 

/* [e] : Allocate memory for global mass and stiffness matrices */

   GMASS = Matrix([9,9]);
   GMASS = ColumnUnits( GMASS, [ kg*m, kg, kg*m, kg, kg*m, kg, kg*m, kg, kg*m ] );
   GMASS =    RowUnits( GMASS, [m], [1] );
   GMASS =    RowUnits( GMASS, [m], [3] );
   GMASS =    RowUnits( GMASS, [m], [5] );
   GMASS =    RowUnits( GMASS, [m], [7] ); 
   GMASS =    RowUnits( GMASS, [m], [9] ); 

   GSTIFF = Matrix([9,9]);  GSTIFF = ColumnUnits( GSTIFF, [ N, N/m,  N,  N/m, N, N/m, N, N/m, N ] );
   GSTIFF =    RowUnits( GSTIFF, [m], [1] );
   GSTIFF =    RowUnits( GSTIFF, [m], [3] );
   GSTIFF =    RowUnits( GSTIFF, [m], [5] );
   GSTIFF =    RowUnits( GSTIFF, [m], [7] );
   GSTIFF =    RowUnits( GSTIFF, [m], [9] );

/* [f] : Assemble Global Stiffness/Mass Matrices for Two Element Cantilever */

   no_elements = 5;
   for( i = 1; i <= no_elements; i = i + 1) {
   for( j = 1; j <= 4; j = j + 1) {

   row = LD [i][j];
   if( row > 0) {
       for( k = 1; k <= 4; k = k + 1) {
          col = LD [i][k];
          if( col > 0) {
              GMASS [ row ][ col ] = GMASS [ row ][ col ] + mass[j][k];
              GSTIFF[ row ][ col ] = GSTIFF[ row ][ col ] + stiff[j][k];
          }
       }
   } 
   }   }

/* [g] : Compute and Print Eigenvalues and Eigenvectors */

   no_eigen    = 2;
   eigen       = Eigen( GSTIFF, GMASS, [ no_eigen ]);
   eigenvalue  = Eigenvalue(eigen);
   eigenvector = Eigenvector(eigen);

   for(i = 1; i <= no_eigen; i = i + 1) {
       print "Mode", i ," : w^2 = ", eigenvalue[i][1];
       print " : T = ", 2*PI/sqrt(eigenvalue[i][1]) ,"\n";
   }

   PrintMatrix(eigenvector);

Points to note are:


OUTPUT FILE

Details of the beam element stiffness and mass matrices, and ensuing eigenvalues and eigenvectors are as follows:

MATRIX : "stiff"

row/col          1            2            3            4          
      units            N          N/m            N          N/m   
   1        m   3.10000e+06  1.16250e+06  1.55000e+06 -1.16250e+06
   2            1.16250e+06  5.81250e+05  1.16250e+06 -5.81250e+05
   3        m   1.55000e+06  1.16250e+06  3.10000e+06 -1.16250e+06
   4           -1.16250e+06 -5.81250e+05 -1.16250e+06  5.81250e+05

MATRIX : "mass"

row/col          1            2            3            4          
      units         kg.m           kg         kg.m           kg   
   1        m   1.92610e+01  2.64838e+01 -1.44457e+01  1.56495e+01
   2            2.64838e+01  4.69486e+01 -1.56495e+01  1.62514e+01  
   3        m  -1.44457e+01 -1.56495e+01  1.92610e+01 -2.64838e+01
   4            1.56495e+01  1.62514e+01 -2.64838e+01  4.69486e+01

*** SUBSPACE ITERATION CONVERGED IN  7 ITERATIONS 

Mode    1.0000e+00  : w^2 =    145.8 rad.sec^-2.0 : T =   0.5203 sec
Mode    2.0000e+00  : w^2 =     1539 rad.sec^-2.0 : T =   0.1602 sec

MATRIX : "eigenvector"

row/col          1            2          
      units                       
   1            1.24900e-01 -2.28500e-01
   2        m   3.02715e-01 -7.71256e-01
   3            1.03061e-01  1.77895e-01
   4        m   8.01556e-01 -9.46803e-01
   5           -1.29310e-02  3.44764e-01
   6        m   1.00000e+00  3.03399e-01
   7           -1.37281e-01 -5.51890e-02
   8        m   6.87189e-01  1.00000e+00
   9           -1.89640e-01 -3.58198e-01

Points to note are:

Figure 2 shows the mode shapes and natural periods of vibration for the supported cantilever beam.

Figure 2 : Modes of Vibration

We have drawn the mode shapes as connected straight lines -- actually, the mode shapes will be cubic splines, as defined by the beam element shape functions.


Developed in April 1996 by Mark Austin
Last Modified April 17, 1996
Copyright © 1996, Mark Austin, Department of Civil Engineering, University of Maryland