dynamics of damped cantilever beam -凯发k8网页登录

modal analysis

create a modal analysis model for a plane-stress problem.

modelm = createpde("structural","modal-planestress");

create the geometry and include it in the model. suppose, the beam is 5 inches long and 0.1 inches thick.

width = 5; 
height = 0.1;
gdm = [3;4;0;width;width;0;0;0;height;height];
g = decsg(gdm,'s1',('s1')');
geometryfromedges(modelm,g);

plot the geometry with the edge labels.

figure; 
pdegplot(modelm,"edgelabels","on"); 
axis equal
title("geometry with edge labels displayed")

define a maximum element size so that there are five elements through the beam thickness. generate a mesh.

hmax = height/5;
msh = generatemesh(modelm,"hmax",hmax);

specify young's modulus, poisson's ratio, and the mass density of steel.

e = 3.0e7; 
nu = 0.3; 
rho = 0.3/386;
structuralproperties(modelm,"youngsmodulus",e, ...
                            "poissonsratio",nu, ...
                            "massdensity",rho);

specify that the left edge of the beam is a fixed boundary.

structuralbc(modelm,"edge",4,"constraint","fixed");

solve the problem for the frequency range from 0 to 1e5. the recommended approach is to use a value that is slightly smaller than the expected lowest frequency. thus, use -0.1 instead of 0.

res = solve(modelm,"frequencyrange",[-0.1,1e5]')

res = 
  modalstructuralresults with properties:
    naturalfrequencies: [8x1 double]
            modeshapes: [1x1 festruct]
                  mesh: [1x1 femesh]

by default, the solver returns circular frequencies.

modeid = 1:numel(res.naturalfrequencies);

express the resulting frequencies in hz by dividing them by 2π. display the frequencies in a table.

tmodalresults = table(modeid.',res.naturalfrequencies/(2*pi));
tmodalresults.properties.variablenames = {'mode','frequency'};
disp(tmodalresults)

    mode    frequency
    ____    _________
     1       126.94  
     2       794.05  
     3       2216.8  
     4       4325.3  
     5       7110.7  
     6       9825.9  
     7        10551  
     8        14623  

compute the analytical fundamental frequency (hz) using the beam theory.

i = height^3/12;
freqanalytical = 3.516*sqrt(e*i/(width^4*rho*height))/(2*pi)

freqanalytical = 126.9498

compare the analytical result with the numerical result.

freqnumerical = res.naturalfrequencies(1)/(2*pi)

freqnumerical = 126.9416

compute the period corresponding to the lowest vibration mode.

longestperiod = 1/freqnumerical

longestperiod = 0.0079

plot the y-component of the solution for the lowest beam frequency.

figure;
pdeplot(modelm,"xydata",res.modeshapes.uy(:,1))
title("lowest frequency vibration mode")
axis equal

initial displacement from static solution

the beam is deformed by applying an external load at its tip and then released at time t=0. find the initial condition for the transient analysis by using the static solution of the beam with a vertical load at the tip.

create a static plane-stress model.

models = createpde("structural","static-planestress");

use the same geometry and mesh that you used for the modal analysis.

geometryfromedges(models,g);
models.mesh = msh;

specify the same values for young's modulus, poisson's ratio, and the mass density of the material.

structuralproperties(models,"youngsmodulus",e, ...
                            "poissonsratio",nu, ...
                            "massdensity",rho);
                       

specify the same constraint on the left end of the beam.

structuralbc(models,"edge",4,"constraint","fixed");

apply the static vertical load on the right side of the beam.

structuralboundaryload(models,"edge",2,"surfacetraction",[0;1]);

solve the static model. the resulting static solution serves as an initial condition for transient analysis.

rstatic = solve(models);

transient analysis

perform the transient analysis of the cantilever beam with and without damping. use the modal superposition method to speed up computations.

create a transient plane-stress model.

modelt = createpde("structural","transient-planestress");

use the same geometry and mesh that you used for the modal analysis.

geometryfromedges(modelt,g);
modelt.mesh = msh;

specify the same values for young's modulus, poisson's ratio, and the mass density of the material.

structuralproperties(modelt,"youngsmodulus",e, ...
                            "poissonsratio",nu, ...
                            "massdensity",rho);
                       

specify the same constraint on the left end of the beam.

structuralbc(modelt,"edge",4,"constraint","fixed");

specify the initial condition by using the static solution.

structuralic(modelt,rstatic)

ans = 
  nodalstructuralics with properties:
    initialdisplacement: [6511x2 double]
        initialvelocity: [6511x2 double]

solve the undamped transient model for three full periods corresponding to the lowest vibration mode.

tlist = 0:longestperiod/100:3*longestperiod;
rest = solve(modelt,tlist,"modalresults",res);

interpolate the displacement at the tip of the beam.

intrput = interpolatedisplacement(rest,[5;0.05]);

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

plot(rest.solutiontimes,intrput.uy)
grid on
title("undamped solution")
xlabel("time")
ylabel("tip of beam displacement")

now solve the model with damping equal to 3% of critical damping.

zeta = 0.03;
omega = 2*pi*freqnumerical;
structuraldamping(modelt,"zeta",zeta);
rest = solve(modelt,tlist,"modalresults",res);

interpolate the displacement at the tip of the beam.

intrput = interpolatedisplacement(rest,[5;0.05]);

the y-displacement at the tip is a sinusoidal function of time with amplitude exponentially decreasing with time.

figure
hold on
plot(rest.solutiontimes,intrput.uy)
plot(tlist,intrput.uy(1)*exp(-zeta*omega*tlist),"color","r")
grid on
title("damped solution")
xlabel("time")
ylabel("tip of beam displacement")