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, . 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, , the following expression relates to and .
In this example, we will define (three percent of critical damping) and equal zero so that can be calculated as
The beam is 5 inches long and 0.1 inches thick.
width = 5; height = 0.1; % 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); c = [2*G+mu;0;G;0;G;mu;0;G;0;2*G+mu]; f = [0;0]; a = 0;
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);
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. figure; pdegplot(g,'EdgeLabels','on'); axis equal title 'Geometry With Edge Labels Displayed'
% Create a pde entity for a PDE with two dependent variables numberOfPDE = 2; model = createpde(numberOfPDE); % Create a geometry entity geometryFromEdges(model,g); % Boundary condition to clamp (displacements equal zero) the left beam-edge applyBoundaryCondition(model,'Edge',4,'u',[0 0]); % Define a maximum element size (5 elements through the beam thickness) hmax = height/5; msh=generateMesh(model,'Hmax',hmax,'MesherVersion','R2013a');
d = rho; % Use pdeeig to solve for the lowest-frequency vibration mode [eigenVec,eigenVal] = pdeeig(model,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; [p,e,t] = meshToPet(msh); 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.29, New conv eig= 2 End of sweep: Basis= 10, Time= 0.29, New conv eig= 2 Basis= 12, Time= 0.67, New conv eig= 1 End of sweep: Basis= 12, Time= 0.67, New conv eig= 1 Lowest beam frequency (Hz). Analytical = 1.269e+02, Numerical = 1.297e+02
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,model,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(model,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
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; pdeTipLoad = createpde(2); pg = geometryFromEdges(pdeTipLoad,g); tipLoad = applyBoundaryCondition(pdeTipLoad,'Edge',2,'g',[0 P/height]); clampedEdge = applyBoundaryCondition(pdeTipLoad,'Edge',4,'u',[0,0]); generateMesh(pdeTipLoad,'Hmax',hmax,'MesherVersion','R2013a'); u = assempde(pdeTipLoad,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 = 0.1; u0 = u0TipDisp/uEigenTip*u; % % Calculate the un-damped solution with the initial displacement from % the static analysis. % u = hyperbolic(u0,0,tlist,model,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 936 function evaluations 1 partial derivatives 207 LU decompositions 935 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 function for creating the plots of tip y-displacement as a function of time.
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