Shear force, bending moment and elastic curve: Matlab program for cantilever beam

Problem 6.31, Hibbeler's "Mechanics of Materials", 8th. edition.

Solution of a cantilever beam, supported on its left end.
___________________________________________________________________

% hibbeler_6_31.m

% Luis Adan Villa Chaparro / Ricardo E. Avila, 24 Feb. 2014
% 

% Problem 6.31 of Hibbeler's "Mechanics of Materials", 8th. edition.

% A cantilever beam of unit length, supported only at its left end,
% has a uniformly distributed unit load w0 on the left half
% of its length; then the load decreases uniformly from w0 to 
% become zero load on the right half of the beam.

clc

clear
disp('Problem 6.31, Hibbeler''s "Mechanics of Materials"')
disp(' ')

L = 1;      % The length of the beam is set to be unit length

w0 = 1;     % The unit load, which may be increased at will.
dx = L/100; % Division of the beam length into segments.
x = [0 : dx : L]';  % Generate values of the x-axis.

n = size(x, 1);     % Size of the x-axis column of data

V = zeros(n, 1);    % Initialize column of shear force data.
M = V;              % Initialize column of bending moment data.
Y = V;              % Initialize column of elastic curve data.

% Processing of the data: shear force, bending moment

% and elastic curve equation, obtained by double integration.
% The two constants of integration are found to be zero, for a
% cantilever beam supported on its left end, at x = 0.

for index = 1 : n

  % First half of the beam length
  V(index) = .75*w0*L - w0*x(index);
  M(index) = (-7*w0*L^2)/24 + .75*w0*L*x(index) - .5*w0*x(index)^2;
  Y(index) = (-7*w0*L^2*(x(index))^2)/48 + ...
             (3*w0*L*(x(index))^3)/24 - (w0*x(index)^4)/24;

  % Second half of the beam length:
  length_2 = x(index) - L/2;  % Define this variable to save time
  if x(index) >= L/2
      V(index) = V(index) + w0*length_2^2/L;
      M(index) = M(index) + w0*length_2^3/(3*L);
      Y(index) = Y(index) + w0*length_2^5/(60*L) ;
  end
end

E = 1; % Elastic constant is unit value

I = 1; % Moment of Inertia is unit value

figure(1)

plot(x, V, 'r', 'linewidth', 3)
grid
xlabel('x, length units')
ylabel('V(x), force units')
title('Shear Force Diagram')

figure(2)

plot(x, M, 'b', 'linewidth', 3)
grid
xlabel('x, length units')
ylabel('M(x), force times distance')
title('Bending Moment Diagram')

figure(3)

plot(x, Y, 'b', 'linewidth', 3)
grid
xlabel('x, length units')
ylabel('y(x), length units')
title('Elastic Curve Diagram')

disp('________________________________ ')

disp('Successful run of Matlab program ')
disp(' ')

% end of Matlab program

Figure(1)

Figure(2)

Figure(3)

Shear force and bending moment diagrams using Matlab, Hibbeler's Mechanics.

Problem 6-11 of "Mechanics of Materials", Hibbeler, 8th. edition is solved using a Matlab program.

Matlab program

% hibbeler_6_11.m

% Ricardo E. Avila Montoya, Feb. 18, 2014

clc

clear

disp(' ')

disp('Solution of Prob. 6-11, Hibbeler,')

disp('"Mechanics of Materials", 8th. ed')

disp('_________________________________')

% Solution of equilibrium equations.

% The unknown forces of the system are Cx, Cy, Dx, Dy

% Matrix of coefficients: distance is measured in feet

A = [-1  0  1  0;     % Equation:   -Cx + Dx = 0;

      0 -1  0  1;     % Equation:   -Cy + Dy = 800;

      0  0 -2 -4;     % Equation:   -2Dx -4Dy = -8000;

      0  0  3 -4];    % Equation:   4Dx -3Dy = 0;

  

% Right Hand Side of equilibrium equations:

RHS = [0; 800; -8000; 0];

% Solve the linear system to find the unknown forces:

Force_vector = A \ RHS; % A Matlab algorithm is used here.

disp(' ')

disp('Forces at points C and D, Lb: ')

disp('_________________________________')

format compact

% Assign values of the equilibrium solution to the forces:

Cx = Force_vector(1)

Cy = Force_vector(2)

Dx = Force_vector(3)

Dy = Force_vector(4)

% From here on, distance measurement unit is inches

L = 120;    % length of A-C main beam, inch

dx = .25;   % discrete division of x-axis, along beam A-C

% Generate discrete values for the x-axis of beam A-C

x = [0 : dx: L]';

N = size(x, 1);   % number of rows in column of x-axis values

% Initialize values of variables to be used in beam problem:

V = zeros(N, 1); % Shear force

M = V;            % Bending moment

y = V;            % elastic curve, y(x) displacement

% Solution of linear system to calculate integration constants:

A = [72 1; 120 1];              % From boundary conditions of beam

RHS = [49766400; 164044800];    % Right-hand side of equations

C = A \ RHS;    % Solution of the linear system of equations

disp(' ')   % display empty line

disp('First constant of integration: ')

disp(C(1))

disp('Second constant of integration: ')

disp(C(2))

% y(x) is obtained by the method of double integration

for j = 1 : N

   % First section of beam, from A to B 

   V(j) =  -800;

   M(j) =  -800 * x(j);

   y(j) =  -133.333 * x(j)^3 + C(1) * x(j) + C(2); 

  

   if x(j) >= 72; % 2nd. section of beam, from B to C 

      V(j) = V(j) + Dy;

      M(j) = M(j) + Dy * (x(j) - 72) + 24 * Dx;

      y(j) = y(j) + (Dy/6) * (x(j) - 72)^3 + 12 * Dx * (x(j) - 72)^2;

   end

end

% Material properties of square hollow section 4" x 4" x 1/4" wall

E = 29E6;   % Lb/in2, ASTM A-36 steel

I = 7.8;    % in^4, from tables of structural materials

% Divide by the product E*I, to obtain elastic curve 

y = y/(E * I);

% Post-processing section

disp('Maximum absolute value of bending moment (Lb-in): ')

disp(max(abs(M)))

figure(1)

plot(x, V, 'r', 'linewidth', 2)

axis([0 120 -850 450])

grid

xlabel('x, in')

ylabel('V(x), in')

title('Shear force diagram')

figure(2)

plot(x, M, 'b', 'linewidth', 2)

grid

xlabel('x, in')

ylabel('M(x), Lb-in')

title('Bending moment diagram')

% The beam is considered as simply soported at points B (x = 72)

% and C (x = 120), for the calculation of the elastic curve.

% The effect of the motion of point D under load, due to bending

% of the projected arm BD, can be added later, and that effect

% superimposed to the elastic curve as it has been calculated here.

figure(3)

plot(x, y, 'k', 'linewidth', 2)

grid

xlabel('x, in')

ylabel('y, in')

title('Elastic curve of beam')

disp('_________________________________')

disp('Successful run of Matlab program. ')

disp(' ')

% end of Matlab program

Figures