A Fixed Free Cantilever Beam Deflection and Stress Analysis with Matlab

This tutorial is related to analyzing the deflection and stress distribution of a cantilever beam. A fixed-free cantilever beam, (Young’s Modulus 𝐸, 𝑏×ℎ cross-section, length 𝐿) is supported at its left-hand end. A Force 𝐹 is applied at a distance 𝑎 from the left-hand end. We will first calculate the shape of the deflected beam and plot the result of deflection, shear and bending moment diagram. Next, we will generate a 2D contour map of the stress distribution along the beam.

Length of the beam, L = 1000 mm

Width of the beam, b = 50 mm

Height of the beam, h = 155 mm

Distance of the force location, a = 800 mm

Modulus of elasticity, E = 200 Gpa

Applied force, F = 25 KN

The shape of the deflected beam is calculated as [2-4],

The moment of inertia is calculated as,

I = 1/12 bh^3

The following Matlab codes are used to calculate the beam deflection and plot the results.

close all

clear 

clc

 

% Parameters from the question

E = 200*10^9;

a = 0.8;

b = 0.05;

h = 0.155;

I = (1/12)*(b*h^3);

f = 25000;

 

x = 0:0.01:0.8;

delta_1 = - (f*x.^2)./(6.*E.*I).*(3.*a - x);  % for 0 <= x <= a

delta_2 = - (f*a.^2)./(6.*E.*I).*(3.*x - a);  % for a <= x <= L

 

figure(1)

plot(x, delta_1)

xlabel('Length (m)')

ylabel('Beam Deflection (m)')

title('Cantilever Beam Deflection for 0 <= x <= a')

 

figure(2)

plot(x, delta_2)

xlabel('Length (m)')

ylabel('Beam Deflection (m)')

title('Cantilever Beam Deflection for a <= x <= L')

The shear force and bending moment diagrams for the fixed free cantilever beam is shown below:

The Matlab codes for the calculation of shear and bending moments as well as plotting these diagrams are provided below.

MATLAB Codes

close all

clear

clc

 

L = 1; % Unit in meters

n = 2; % Number of point loads on the beam

 

for i = 1:n

fprintf('Enter the point load and distance where load acts for %d node\n',i)

Wc(i)=input('Enter the load in Newton\n');

Lc(i) = input('Enter the distance of the point load from the fixed end\n');

            

if Lc(i)>L || Lc(i)<0                

    error('Please check the Length')            

end

end

 

NL = zeros(1,n);        

NW = zeros(1,n);

       

for i=1:n

    [a,b]=max(Lc);            

    NL(i)=a;           

    NW(i)=Wc(b);          

    Lc(b)=[];           

    Wc(b)=[];        

end

 

NL(n+1) = 0;

        

% Shear force diagram

figure(1)

Ra = sum(NW);

if NL(1)==L

X1 = L;

F1 = NW(1);

elseif i==1

X1 = [L NL(1) NL(1) 0];

F1 = [0 0 Ra Ra];

else

X1 = [L NL(1)];

F1 = [0 0];

end

S=[]; X=[]; F=[];

        if n>=2

            for i=2:n+1

                x = [NL(i-1) NL(i)];

                S= [S NW(i-1)];

                f = sum(S);

                X= [X x];

                F = [F f f];

            end

        end

subplot(2,1,1)

        

plot([X1 X],[F1 F],'b')

xlim([0 L]);

hline = refline(0,0);

hline.Color = 'k';

legend('Shear Force','Reference')          

title('Shear Force Diagram of Cantilever Beam')

xlabel('Length of the Beam in Meter')

ylabel('Shear Force in Newton')

  

hold on

grid on

fprintf('Reaction Force =%d N\n',Ra)

    

% Bending moment diagram

subplot(2,1,2)

if NL(1)<L

X1=NL(1): L;

M1 = zeros(size(X1));

plot(X1,M1,'r')

hold on 

grid on

end

 

Xm=[]; 

Mm=[];

    

for i=1:n

X = NL(i+1):0.1:NL(i);

M = -(NW(1,1:i)*((NL(1,1:i))'-newX(X,i))) ;

Xm =[X Xm];

Mm =[M Mm];

end    

    

plot(Xm,Mm,'r') 

xlim([0 L]);

hline = refline(0,0);

hline.Color = 'k';

legend('Bending Moment','Reference')  

    

title('Bending Moment Diagram of Cantilever Beam')

xlabel('Length of the Beam in Meter')

ylabel('Bending Moment in Newton-Meter')

hold off

grid on

 

function x = newX(X,i)                                                

[~,d] = size(X);                                                 

x = zeros(i,d);                                                

for j=1:i                                                        

x(j,1:d) = X;                                                   

end                                                             

end    

The bending stress in the beam is calculated as [2-4],

Where, δ is calculated as,

In the following, the stress distributions for the cantilever beam are shown. The Matlab codes are appended after the results.

MATLAB CODES

close all

clear 

clc

 

% Parameters from the question

L = 1;

E = 200*10^9;

a = 0.8;

b = 0.05;

h = 0.155;

I = (1/12)*(b*h^3);

f = 25000;

M1 = f*a;

M2 = f*(L-a);

 

x = 0:0.01:0.8;

y = 0:0.01:0.8;

 

delta_1 = - (f*x.^2)./(6.*E.*I).*(3.*a - x);  % for 0 <= x <= a

delta_2 = - (f*a.^2)./(6.*E.*I).*(3.*x - a);  % for a <= x <= L

 

sigma_1 = (M1.*delta_1)./I;

 

x = rand(100,1)*4-2;

y = rand(100,1)*4-2;

z = (M1.*((f*x.^2)./(6.*E.*I).*(3.*a - x)))./I;

z1 = (M2.*((f*a.^2)./(6.*E.*I).*(3.*x - a)))./I;

 

F = TriScatteredInterp(x,y,z);

ti = -1:.25:1;

[qx,qy] = meshgrid(ti,ti);

qz = F(qx,qy);

 

figure(1)

scatter3(x,y,z)

title('Mesh Plot for 0 <= x <= a')

hold on

mesh(qx,qy,qz)

 

figure(2)

contour(qx,qy,qz)

title('Contour Plot for 0 <= x <= a')

 

F = TriScatteredInterp(x,y,z1);

ti = -1:0.25:1;

[qx,qy] = meshgrid(ti,ti);

qz = F(qx,qy);

 

figure(3)

scatter3(x,y,z)

title('Mesh Plot for a <= x <= L')

hold on

mesh(qx,qy,qz)

 

figure(4)

contour(qx,qy,qz)

title('Contour Plot for a <= x <= L')

REFERENCES

[1] MATLAB 2020b Academic Version from MathWorks (https://www.mathworks.com/)

[2] Budynas-Nisbett, "Shigley's Mechanical Engineering Design," 8th Ed.

[3] Gere, James M., "Mechanics of Materials," 6th Ed.

[4] Lindeburg, Michael R., "Mechanical Engineering Reference Manual for the PE Exam," 13th Ed.

[5] Solid Mechanics Part I: An Introduction to Solid Mechanics by ‪Piaras Kelly‬‬

#FixedFreeBeam   #CantileverBeam   #StressDistribution   #ShearForce   #BendingMoment  #Matlab