Documentation Center

  • Trial Software
  • Product Updates

Dynamics of a Damped Cantilever Beam

This example shows how to include damping in the transient analysis of a simple cantilever beam analyzed with the Partial Differential Equation Toolbox™. 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:

$$ D = \alpha M + \beta K $$

It is common to express damping as a percentage of critical damping, $\xi$ , for a selected vibration frequency. For a given frequency, $\omega_i$ , the following expression relates $\xi$ to $\alpha$ and $\beta$ .

$$ \xi = \frac{\alpha}{2 \omega_i} + \frac{\beta\omega_i}{2} $$

In this example, we will define $\xi = .03$ (three percent of critical damping) and $\beta$ equal zero so that $\alpha$ can be calculated as $\alpha = 2 \times .03 \omega_i$

Beam Dimensions and Material Properties

The beam is 5 inches long and .1 inches thick.

width = 5; height = .1;
% The material is steel
E = 3.0e7; nu = .3; rho = .3/386;
% Calculate the coefficient matrix from the material properties
G = E/(2.*(1+nu));
mu = 2.0*G*nu/(1-nu);
c = [2*G+mu; 0; G;   0; G; mu; 0;  G; 0; 2*G+mu];
f = [0 0]';
a = 0;
% Boundary function to clamp (displacements equal zero) the left beam-edge
b = @cantileverBeamBoundaryFile;

Lowest Vibration Frequency From Beam Theory

I = height^3/12; A = height;
% From beam theory, there is a simple expression for the lowest vibration
% frequency of the cantilever beam.
eigValAnalytical = 1.8751^4*E*I/(rho*A*width^4);
freqAnalytical = sqrt(eigValAnalytical)/(2*pi);

Geometry and Mesh

Create a simple rectangular geometry

gdm = [3 4 0 width width 0 0 0 height height]';
g = decsg(gdm, 'S1', ('S1')');
% Define a maximum element size (5 elements through the beam thickness)
hmax = height/5;
[p, e, t] = initmesh(g, 'Hmax', hmax, 'MesherVersion', 'R2013a');

Calculation of Vibration Modes and Frequencies

d = rho;
% Use pdeeig to solve for the lowest-frequency vibration mode
[eigenVec, eigenVal] = pdeeig(b,p,e,t,c,a,d,[0, 1e6]);
freqNumerical = sqrt(eigenVal(1))./(2*pi);
longestPeriod = 1/freqNumerical;
% Color plot of y-displacement of the lowest-frequency vibration mode
u = eigenVec(:,1);
u2 = reshape(u, [], 2); % re-dimension eigenvector for convenience
% Plot the deformed shape of the beam with the displacements scaled
% by an arbitrary factor.
scaleFactor = 20;
figure; pdeplot(p+scaleFactor*u2', e, t, 'xydata', real(u2(:,2)));
title('Lowest Frequency Vibration Mode'); axis equal;
xlabel('Inches'); ylabel('Inches'); drawnow;

fprintf('Lowest beam frequency (Hz). Analytical = %12.3e, Numerical=%12.3e\n', ...
  freqAnalytical, freqNumerical);
              Basis= 10,  Time=   0.17,  New conv eig=  2
End of sweep: Basis= 10,  Time=   0.17,  New conv eig=  2
              Basis= 12,  Time=   0.20,  New conv eig=  1
End of sweep: Basis= 12,  Time=   0.20,  New conv eig=  1
Lowest beam frequency (Hz). Analytical =    1.269e+02, Numerical=   1.297e+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 = u2(2,2); u0TipDisp = .1;
u0 = u0TipDisp/uEigenTip*u;
%
% First solve the undamped system.
%
% Calculate the solution for three full periods
tlist = 0:longestPeriod/100:3*longestPeriod;
u = hyperbolic(u0,0, tlist, b,p,e,t,c,a,f,rho);
titl = 'Initial displacements from lowest eigenvector, un-damped solution';
cantileverBeamTransientPlot(u, tlist, titl);
% As can be seen in the figure below, 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.
%
% Now solve the system with damping equal to 3% of critical for this
% lowest vibration frequency.
%
[K, F, B, ud] = assempde(b,p,e,t,c,a,f);
[~,M] = assema(p,t,0,rho,f);
zeta = .03; omega = 2*pi*freqNumerical; alpha = 2*zeta*omega;
u = hyperbolic(u0, 0, tlist, K, F, B, ud, M, 'DampingMatrix', alpha*M);
titl = 'Initial displacements from lowest eigenvector, damped solution';
cantileverBeamTransientPlot(u, tlist, titl);
hold on;
plot(tlist, u0TipDisp*exp(-zeta*omega*tlist), 'color', 'r');
legend('PDE', 'ODE Amplitude Decay', 'location', 'southeast');
%
% The figure below 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.
154 successful steps
33 failed attempts
249 function evaluations
1 partial derivatives
61 LU decompositions
248 solutions of linear systems
165 successful steps
40 failed attempts
278 function evaluations
1 partial derivatives
72 LU decompositions
277 solutions of linear systems

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 it's 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.

P = 1.0;
bTipLoad = @(p, e, u, time) cantileverBeamBoundaryFileTipLoad(p, e, P/height);
u = assempde(bTipLoad,p,e,t,c,a,f);
% 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.
u2 = reshape(u, [], 2);
uEigenTip = u2(2,2); u0TipDisp = .1;
u0 = u0TipDisp/uEigenTip*u;
%
% Calculate the un-damped solution with the initial displacement from
% the static analysis.
%
u = hyperbolic(u0,0, tlist, b,p,e,t,c,a,f,rho);
titl = 'Initial displacements from static solution, un-damped solution';
cantileverBeamTransientPlot(u, tlist, titl);
% We see in the figure below that the displacement is no longer a pure sin 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.

%
% 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.
%
u = hyperbolic(u0, 0, tlist, K, F, B, ud, M, 'DampingMatrix', alpha*M);
titl = 'Initial displacements from static solution, damped solution';
cantileverBeamTransientPlot(u, tlist, titl);
hold on;
plot(tlist, u0TipDisp*exp(-zeta*omega*tlist), 'color', 'r');
legend('PDE', 'ODE Amplitude Decay', 'location', 'southeast');
%
% In the figure below, we 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.
%
639 successful steps
89 failed attempts
937 function evaluations
1 partial derivatives
208 LU decompositions
936 solutions of linear systems
626 successful steps
93 failed attempts
929 function evaluations
1 partial derivatives
209 LU decompositions
928 solutions of linear systems

Utility Plot Function

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

type cantileverBeamTransientPlot.m
function cantileverBeamTransientPlot( u, tlist, titl )
%CANTILEVERBEAMTRANSIENTPLOT Plot y-displacement at the beam tipe
%   u displacements as a function of time
%   tlist vector of output times
%   titl plot title


% Copyright 2013 The MathWorks, Inc.

uu = reshape(u, [], 2, size(tlist,2));
utip = uu(2,2,:);
figure; plot(tlist, utip(:)); grid on;
title(titl); xlabel('Time, seconds');
ylabel('Y-displacement at beam tip, inches');
drawnow;
end

Basic Boundary Condition Function

Boundary file function that clamps (x and y displacements equal zero) at the left end of the beam.

type cantileverBeamBoundaryFile.m
function [ q, g, h, r ] = cantileverBeamBoundaryFile( p, e, u, time )
%CANTILEVERBEAMBOUNDARYFILE Boundary conditions for cantilever beam example
%   [ q, g, h, r ] = boundaryFileClamped( p, e, u, time ) returns the
%   Neumann BC (q, g) and Dirichlet BC (h, r) matrices for the
%   damped cantilever beam example.
%   p is the point matrix returned from INITMESH
%   e is the edge matrix returned from INITMESH
%   u is the solution vector
%   time is the current time

% Copyright 2013 The MathWorks, Inc.

N = 2; ne = size(e,2);
q = zeros(N^2, ne); g = zeros(N, ne);
h = zeros(N^2, 2*ne); r = zeros(N, 2*ne);
for i=1:ne
  geomEdgeID = e(5,i);
  if(geomEdgeID == 4)
    % clamped (u=v=0) at left geometry edge
    h(1,i) = 1; h(1,i+ne) = 1;
    h(4,i) = 1; h(4,i+ne) = 1;
  end
end

end

Boundary Condition Function With Applied Load

Boundary file function that clamps the left end of the beam and applies a shear stress at the right end of the beam.

type cantileverBeamBoundaryFileTipLoad.m
function [ q, g, h, r ] = cantileverBeamBoundaryFileTipLoad( p, e, s )
%CANTILEVERBEAMBOUNDARYFILETIPLOAD Boundary conditions for cantilever beam example
%   [ q, g, h, r ] = boundaryFileTipLoad( p, e, V ) returns the
%   Neumann BC (q, g) and Dirichlet BC (h, r) matrices for the
%   damped cantilever beam example.
%   p is the point matrix returned from INITMESH
%   e is the edge matrix returned from INITMESH
%   s is shear stress applied to the right beam edge

% Copyright 2013 The MathWorks, Inc.

N = 2; ne = size(e,2);
q = zeros(N^2, ne); g = zeros(N, ne);
h = zeros(N^2, 2*ne); r = zeros(N, 2*ne);
for i=1:ne
  geomEdgeID = e(5,i);
  if(geomEdgeID == 2)
    % right geometry edge
    g(2,i) = s; % apply a shear load
  elseif(geomEdgeID == 4)
    % left geometry edge
    h(1,i) = 1;
    h(1,i+ne) = 1;
    h(4,i) = 1;
    h(4,i+ne) = 1;
  end
end

end

Was this topic helpful?