Dynamics of Damped Cantilever Beam

This example shows how to include damping in the transient analysis of a simple cantilever beam.

The beam is modeled with a plane stress elasticity formulation. The damping model is basic viscous damping distributed uniformly through the volume of the beam. Several transient analyses are performed where the beam is deformed into an initial shape and then released at time, t=0. Analyses with and without damping are considered. Two initial displacement shapes are considered. In the first, the beam is deformed into a shape corresponding to the lowest vibration mode. In the second, the beam is deformed by applying an external load at the tip of the beam. No additional loading is applied in this example so, in the damped cases, the displacement of the beam decays as a function of time due to the damping.

The transient analyses are performed using the PDE Toolbox hyperbolic function. One form of this function allows a transient analysis to be performed with the stiffness, mass, and damping matrices and load vectors as input. Typically these matrices and vectors are calculated using other PDE Toolbox functions. That approach will be demonstrated in this example.

A particularly simple way to construct a damping matrix is by using what is commonly referred to as Rayleigh damping. With Rayleigh damping, the damping matrix is defined as a linear combination of the mass and stiffness matrices:


It is common to express damping as a percentage of critical damping, ξ, for a selected vibration frequency. For a given frequency, ωi, the following expression relates ξ to α and β.


In this example, we will define ξ=0.03 (three percent of critical damping) and β equal zero so that α can be calculated as α=2×.03ωi

This example specifies values of parameters using the imperial system of units. You can replace them with values specified in the metric system. If you do so, ensure that you specify all values throughout the example using the same system.

Beam Parameters

Specify the beam dimensions. It is 5 inches long and 0.1 inches thick.

width = 5; 
height = 0.1;

Specify that the material is steel.

E = 3.0e7; 
nu = 0.3; 
rho = 0.3/386;

Calculate the coefficient matrix from the material properties.

G = E/(2.*(1+nu));
mu = 2.0*G*nu/(1-nu);

From beam theory, there is a simple expression for the lowest vibration frequency of the cantilever beam.

I = height^3/12; 
A = height;

eigValAnalytical = 1.8751^4*E*I/(rho*A*width^4); 
freqAnalytical = sqrt(eigValAnalytical)/(2*pi);

Problem Setup

Create a PDEModel with two independent variables to represent the analysis.

numberOfPDE = 2;
model = createpde(numberOfPDE);

Create a simple rectangular geometry.

gdm = [3;4;0;width;width;0;0;0;height;height];
g = decsg(gdm,'S1',('S1')');

Plot the geometry and display the edge labels for use in the boundary condition definition.

axis equal
title 'Geometry With Edge Labels Displayed'

Provide the model with a definition of the geometry.


Specify the coefficients by deriving them from the material properties.

c = [2*G+mu;0;G;0;G;mu;0;G;0;2*G+mu];
f = [0;0];
a = 0;
m = rho;

Specify the following boundary condition to clamp (displacements equal zero) the left beam-edge.

applyBoundaryCondition(model,'dirichlet','Edge',4,'u',[0 0]);

Define a maximum element size (5 elements through the beam thickness) and generate a mesh.

hmax = height/5;

Calculation of Vibration Modes and Frequencies

Use solvepdeeig and then compute the lowest-frequency vibration mode.

res = solvepdeeig(model, [0,1e6]');
              Basis= 10,  Time=   0.28,  New conv eig=  2
End of sweep: Basis= 10,  Time=   0.28,  New conv eig=  2
              Basis= 12,  Time=   0.36,  New conv eig=  1
End of sweep: Basis= 12,  Time=   0.37,  New conv eig=  1
eigenVec = res.Eigenvectors;
eigenVal = res.Eigenvalues;

Color plot of y-displacement of the lowest-frequency vibration mode.

freqNumerical = sqrt(eigenVal(1))./(2*pi);
longestPeriod = 1/freqNumerical;

Plot the deformed shape of the beam with the displacements scaled by an arbitrary factor.

scaleFactor = 20;
[p,e,t] = meshToPet(msh);
title('Lowest Frequency Vibration Mode');
axis equal


fprintf('Lowest beam frequency (Hz). Analytical = %8.3e, Numerical = %8.3e\n', ...
Lowest beam frequency (Hz). Analytical = 1.269e+02, Numerical = 1.269e+02

Transient Analysis, Initial Displacement From First Mode Shape

In the first two transient analyses, we define an initial displacement. in the shape of the lowest vibration mode. By doing this, we convert the PDE to a single ODE with time as the independent variable. The solution to this ODE is the same as that of the classical spring-mass-damper system with a frequency equal the frequency of this vibration mode. Thus we are able to compare the numerical solution with the analytical solution to this well-known problem.

For convenience, we will scale the eigenvector shape so that y-displacement at the tip is .1 inches. This makes the transient history plots simpler.

uEigenTip = eigenVec(2,2);
u0TipDisp = .1;

u0 = u0TipDisp/uEigenTip*eigenVec;

First solve the undamped system.

Calculate the solution for three full periods.

tlist = 0:longestPeriod/100:3*longestPeriod;

Create a function handle that can be used to provide initial conditions.

R = createPDEResults(model, u0(:));
ice = icEvaluator(R);

Set the initial conditions and solve the system.

setInitialConditions(model, @ice.computeIC, 0);
tres = solvepde(model,tlist);

The displacement at the tip is a sinusoidal function of time with amplitude equal to the initial y-displacement. This agrees with the solution to the simple spring-mass system.

titl = 'Initial Displacements from Lowest Eigenvector, Undamped Solution';

Now solve the system with damping equal to 3% of critical for this lowest vibration frequency.

fem = assembleFEMatrices(model);
zeta = .03; 
omega = 2*pi*freqNumerical; 
alpha = 2*zeta*omega;
dampmat = alpha*fem.M;
tres = solvepde(model,tlist);

This figure shows the y-displacement at the tip as a function of time. Superimposed on this plot is a second curve which shows the envelope of the maximum amplitude as a function of time, calculated from the solution to the single degree-of-freedom ODE. As expected, the PDE solution agrees well with the analytical solution.

titl = 'Initial Displacements from Lowest Eigenvector, Damped Solution';
hold on;
legend('PDE','ODE Amplitude Decay','Location','southeast');

Transient Analysis, Initial Displacement From Static Solution

It would be very unusual for a structure to be loaded such that the displacement is equal to a multiple of one of its vibration modes. In this more realistic example, we solve the transient equations with the initial displacement shape calculated from the static solution of the cantilever beam with a vertical load at the tip.

Perform a static analysis with vertical tip load equal one. Follow the model building steps as before.

P = 1.0;
pdeTipLoad = createpde(2);
pg = geometryFromEdges(pdeTipLoad,g);

Specify the equation to solve.

tipLoad = applyBoundaryCondition(pdeTipLoad,'neumann','Edge',2,'g',[0 P/height]);
clampedEdge = applyBoundaryCondition(pdeTipLoad,'dirichlet','Edge',4,'u',[0,0]);

statres = solvepde(pdeTipLoad);

To make comparison with the eigenvector case clearer, we will also scale this static solution so that the maximum y-displacement at the tip equals .1 inches.

u = statres.NodalSolution;
uEigenTip = u(2,2); 
u0TipDisp = 0.1;
u0 = u0TipDisp/uEigenTip*u;

Calculate the undamped solution with the initial displacement from the static analysis.

specifyCoefficients(model, 'm', m, 'd', 0, 'c', c, 'a', a, 'f', f);

Set the initial conditions and solve the system.

R = createPDEResults(model, u0(:));
ice = icEvaluator(R);
setInitialConditions(model, @ice.computeIC, 0);
tres = solvepde(model,tlist);

The displacement is no longer a pure sinusoidal wave. The static solution that we are using as the initial conditions is similar to the lowest eigenvector but higher-frequency modes are also contributing to the transient solution.

titl = 'Initial Displacements from Static Solution, Undamped Solution';

Calculate the damped solution with the initial displacement from the static analysis. The damping matrix is the same as that used in the eigenvector case, above.

specifyCoefficients(model, 'm', m, 'd', dampmat, 'c', c, 'a', a, 'f', f);
tres = solvepde(model,tlist);

Plot the tip displacement from this solution as a function of time. Again we superimpose a curve of the damped amplitude as a function of time obtained from an analytical solution to the single degree-of-freedom ODE. Because the initial condition differs from the lowest eigenvector, this analytical solution only approximates the amplitude of the PDE solution.

titl = 'Initial Displacements from Static Solution, Damped Solution';
hold on;
legend('PDE','ODE Amplitude Decay','Location','southeast');

Utility Plot Function

Utility function for creating the plots of tip y-displacement as a function of time.

type cantileverBeamTransientPlot.m
function cantileverBeamTransientPlot( tdres, titl )
%CANTILEVERBEAMTRANSIENTPLOT Plot y-displacement at the beam tip
%   tdres - Time-dependent results object representing displacements as a function of time
%   titl plot title

% Copyright 2013-2015 The MathWorks, Inc.

tlist = tdres.SolutionTimes;
uu = tdres.NodalSolution;
utip = uu(2,2,:);
figure; plot(tlist, utip(:)); grid on;
title(titl); xlabel('Time, seconds');
ylabel('Y-displacement at beam tip, inches');

Function Handle for Specifying Initial Condition (IC)

The evaluation function returns the value of the IC at any point within the mesh. A results object represents the values and the interpolation function that it provided is used to compute the ICs.

type icEvaluator.m
classdef icEvaluator
% icEvaluator Evaluates Initial Conditions data at requested locations
% ICE = icEvaluator(R) Creates an initial conditions evaluator from a 
% results object. The evaluator provides a function that can be called at 
% specific locations within the geometry. This class customized to represent
% a stationary solution.
% icEvaluator methods:
%    computeIC - Computes ICs at locations within the geometric domain

% Copyright 2015 The MathWorks, Inc.              

        function obj = icEvaluator(R)                 
            obj.Results = R;                                                                                                                                 
        function ic = computeIC(obj,locations)
            is2d = size(obj.Results.Mesh.Nodes, 1) == 2;
            if is2d
                querypts = [locations.x; locations.y];
                  querypts = [locations.x; locations.y; locations.z];
            numeqns = size(obj.Results.NodalSolution,2);
            if numeqns == 1 
                ic = interpolateSolution(obj.Results, querypts);
                ic = zeros(numeqns, numel(locations.x));
                for i = 1:numeqns
                    ic(i,:) = interpolateSolution(obj.Results, querypts, i);