Wednesday, February 19, 2014

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


No comments:

Post a Comment