C++ Visual Studio programs to plot shear force, bending moment and elastic curve for loaded beams.

This Visual Studio C++ project can be used to plot any function, as long as it can be
programmed using C++. I use it to plot diagrams for the design of loaded beams in
mechanical structures.

It is relatively easy to plot diagrams for shear force, bending moment and elastic curves of loaded
beams, particularly when the reactions at the supports of the beams are determined simply by the methods of statics. Examples using Matlab have been posted in this blog before.

Now there are two applications that work together to produce those diagrams, useful for analysis
and design of loaded beams.

The Visual Studio .NET projects Mechanics1 and Graph are compressed in the file Graphs.zip
in the following page:

http://www.ricardoavila.org/mechanics

There is a Word file (Graph_readme) there in the compressed file, with instructions on how to use
the Visual Studio C++ applications.

The results of using those files are the responsibility of the user.

Design of beam: diagrams using Matlab, Hibbeler's F1.7 problem, page 38.

Matlab program to plot the mechanical diagrams for Hibbeler's
"Mechanics of Materials",
8th. edition, problem F1.7, page 38.

For this problem we assume a maximum distributed load
w = 12000 N/m.

The resulting forces in the vertical supporting rods are:
Rod AB force = 12000 N.
Rod CD force = 24000 N.
_________________________________________________________________
 
% hibbeler_f1_7.m
% Joel Aguinaga / Ricardo Avila, 18 Aug. 2014
% Fundamental Problem F1.7, Hibbeler's
% "Mechanics of Materials", 8th. edition, page 38

clc

clear

disp(' ')

disp('Hibbeler''s Fundamental problem F7.1 ')
disp('"Mechanics of Materials", 8th. edition, p.38 ')
disp('******************************************** ')
disp(' ')

% Data pre-processing

delta_x = .01;  % Define step of discretization for x-axis
x = [ 0: delta_x: 6]';  % Generate x-axis discretization.

N = size(x, 1);  % Calculate size of the discretization

V = zeros(N, 1); % Fill with zeros variable in memory: shear force. 
M = V;           % Memory required to store Bending Moment.
y = V;           % Memory to store elastic curve information.

% Material and Geometric Constants:

E = 200E9;  % Assume a steel beam, for Elastic Modulus
I = 1;      % Moment of inertia is not determined, use unit value
EI = E * I; % Multiply those values

% Data processing section of program: 

% ____________________________________________________________
V = -1000* x.^2 + 12000;    % Equation for Shear Force curve
M = -333.3333* x.^3 + 12000*x;  % Bending Moment equation.

% The elastic curve resulting from double integration,

% following Euler's method: EI * (d2y/dx2) = M(x).
y = -16.6666* x.^5 + 2000 *x.^3 - 50400 *x;

% Divide elastic curve by EI

y = y / EI;

% Post-processing section of program:

%_____________________________________________________________

% Plot the Shear Force Diagram

figure(1)
plot(x, V, 'r', 'linewidth', 3)
grid
title(' Shear Force Diagram')
xlabel('x, m')
ylabel('V, N')

% Plot the Bending Moment Diagram

figure(2)
plot (x, M, 'b', 'linewidth', 3)
grid
title('Bending Moment Diagram')
xlabel('x, m')
ylabel('M, N-m')

% Plot the Elastic Curve Diagram

figure(3)
plot (x, y, 'k', 'linewidth', 3)
grid
title('Elastic Curve Diagram')
xlabel('x, m')
ylabel('y, m')

% Find and display the minimum value of elastic curve

y_min = min(y)

% Find and display the maximum value of bending moment
M_max = max(M)

disp(' ')

disp('**********************')
disp('Successful Program Run')


Figure 1. Shear Force Diagram:

 

 

 

Figure 2. Bending Moment Diagram:

 

Figure 3. Elastic Curve Diagram:

 

2D truss: Solution using Principle of Virtual Work with Matlab program

This Matlab program gives the solution to problem 11.3-12 of the
textbook "Mechanics of Materials", by Roy R. Craig, Jr., 2nd. edition, 2000.

Units are US customary system, Lb and inch.

Areas of elements:
Element (1):  1    in^2
Element (2):  1    in^2
Element (3):  1.5 in^2
Element (4):  1.5 in^2

NODES OF STRUCTURE (TRUSS):
Node A = Node 1. Coordinates: (0, 0, 0)
Node B = Node 2. Coordinates: (40, 0, 0)
Node C = Node 3. Coordinates: (0, -30, 0)
Node D = Node 4. Coordinates: (70, -30, 0)

All elements are straight links of steel, with modulus of elasticity E = 30E6 Lb/in^2.
The structure is a truss, with pins joining the elements at the nodes A, B, C and D.

ANSYS SOLUTION FOR LOAD STEP (LS) 1:

Matlab program:
______________

% craig_11_3_12.m

% Ricardo E. Avila Montoya / Luis Adan Villa Chaparro

% University of Ciudad Juarez, Mexico. 9 May 2014

% e-mail: ricardo_avila@hotmail.com
% Chihuahua, Mexico

% The method of virtual work is applied to solve a 2D

% structure, formed with elastic links pinned at their ends.

clc

clear

A = [1; 1; 1.5; 1.5];      % in^2, Area of link elements.

E = 30e6;                  % Lb/in^2, Elastic modulus of steel. 

L = [40; 42.426; 50; 70];  % Lengths of link elements.

K = A * E./L;              % Calculate element stiffness.

% Initialize the Stiffness Matrix

SM = zeros(4, 4);

% Build the Stiffness Matrix

SM(1, 1) =       K(1) + .5 * K(2) + .64 * K(3);

SM(1, 2) =            - .5 * K(2) + .48 * K(3);

SM(1, 3) =            - .5 * K(2);

SM(1, 4) =              .5 * K(2);

SM(2, 1) = SM(1, 2);

SM(2, 2) =              .5 * K(2) + .36 * K(3);

SM(2, 3) =              .5 * K(2);

SM(2, 4) =            - .5 * K(2);

SM(3, 1) = SM(1, 3);

SM(3, 2) = SM(2, 3);

SM(3, 3) =              .5 * K(2)              + K(4);

SM(3, 4) =            - .5 * K(2);

SM(4, 1) = SM(1, 4);

SM(4, 2) = SM(2, 4);

SM(4, 3) = SM(3, 4);

SM(4, 4) =              .5 * K(2);

% Forces applied on the structure

RHS = [0; 0; 0; -2000];

% Calculate and display displacements of nodes B and D 

disp('Displacement of nodes B and D, inch')

format long

% Solve the system of linear equations using Matlab algorithm
% (use the inverted backlash Matlab operator for matrix solution)

U = SM\RHS

% Calculation of forces

K2_sqrt_2 = K(2)/sqrt(2); % Auxiliary arithmetic operation

Force_Matrix = [ K(1)       0          0          0

                -K2_sqrt_2  K2_sqrt_2  K2_sqrt_2 -K2_sqrt_2

                 .8*K(3)    .6*K(3)    0          0

                 0          0          K(4)       0 ];

% Calculate and display vector of forces for the elements             

Force_vector = Force_Matrix * U

% Calculate stress in link elements of 2D structure:

Stress_in_elements = Force_vector ./ A

disp('________________________________')

disp('Successful run of Matlab program')

_______________________________________________________________
ANSYS ANIMATION OF DEFORMED RESULTS:

 

___________________________________________________________________
ANSYS results for displacements of the nodes:

 PRINT U    NODAL SOLUTION PER NODE

  ***** POST1 NODAL DEGREE OF FREEDOM LISTING *****                           

  LOAD STEP=     2  SUBSTEP=     1                                            
   TIME=    2.0000      LOAD CASE=   0                                        

  THE FOLLOWING DEGREE OF FREEDOM RESULTS ARE IN GLOBAL COORDINATES           

    NODE      UX          UY          UZ          USUM 
       1   0.0000      0.0000      0.0000      0.0000   
       2  0.62222E-02-0.14469E-01  0.0000     0.15750E-01
       3   0.0000      0.0000      0.0000      0.0000   
       4 -0.31111E-02-0.29459E-01  0.0000     0.29623E-01

 MAXIMUM ABSOLUTE VALUES
 NODE          2           4           0           4
 VALUE   0.62222E-02-0.29459E-01  0.0000     0.29623E-01
_________________________________________________________________

ANSYS results of forces and stresses in the elements of the truss:


 PRINT ELEMENT TABLE ITEMS PER ELEMENT

  ***** POST1 ELEMENT TABLE LISTING *****                                     

    STAT    CURRENT     CURRENT
    ELEM    FORCE       STRESS 
       1   4666.7      4666.7   
       2   2828.4      2828.4   
       3  -3333.3     -2222.2   
       4  -2000.0     -1333.3   

 MINIMUM VALUES
 ELEM          3           3
 VALUE   -3333.3     -2222.2   
 MAXIMUM VALUES
 ELEM          1           1
 VALUE    4666.7      4666.7   

Hibbeler Problem 6-35. Shear force, bending moment and elastic curve diagrams.

Problem 6-35, 8th. edition of Hibbeler's "Mechanics of Materials".

A Matlab program to calculate and plot shear force, bending moment and elastic curve
of the loaded beam.

Length of beam is 6 m. Establish the x axis with its origin (x = 0) at the left end of the beam.

The beam is simply supported at its ends (x = 0 at left end, x = 6 at right end).

The load on the beam can be divided to form two loads, and thus simplify the problem:

1. A distributed uniform load of 200 N/m is located on the right half, 3 to 6 m of the beam.

2. A distributed uniformly increasing load starts from 0 N/m at x = 3 m, and it increases
to 200 N/m at the right end of the beam, x = 6 m.

Calculation of reaction forces at the supports: lump the loads, then there is a concentrated
load of 600 N located at x = 4.5, and another concentrated force of 300 N at x = 5.

The reactions are found to be RA = 200N, and RB = 700 N.

The figures below show the calculations to find bending moment M(x), shear force V(x) and
the elastic curve y(x).

Click on the figures to make them larger.

 

MATLAB PROGRAM:

% hibbeler_6_35.m
% Irma Iveth Martinez Huerta / Ricardo E. Avila
% mechanicsofmaterials.blogspot.com / ricardo_avila@hotmail.com

clc

clear

disp(' ')

disp('Problem 6-35, Hibbeler''s ')
disp('"Mechanics of Materials", 8th. ed.')
disp(' ')

% Pre-processing section

L = 6;              % Total length of beam
dx = L/600;         % x-axis division is 1 centimeter
x = [0 : dx : L]';  % x-axis is defined
n = size(x, 1);     % Size of the x-axis column of data

% Initialize data structures

V = zeros(n, 1);   
M = V;             
Y = V;

% Processing section of Matlab program

% Since no values of E (elastic modulus) and I (moment of inertia)
% are given, we use E = 1 and I = 1 to plot the elastic curve.

% First constant of integration, calculated from boundary condition

% x = 0 @ y = 6:
C1 = -1065;

for index = 1 : n

     % First portion of the beam, 0 < x < 3

     V(index) = 200;
     M(index) = 200 * x(index);
     Y(index) = 33.3333 * x(index)^3 + C1 * x(index);

     % Second portion of the beam, 3 < x < 6

      if  x(index) >= 3
          V(index) = V(index) - 200 * (x(index)-3) ...
                              - 33.3333 *((x(index)-3)^2);

          M(index) = M(index) - 100 *((x(index)-3)^2) ...

                              - 11.1111 * ((x(index)-3)^3);

          Y(index) = Y(index) -  8.3333 * ((x(index)-3)^4) ...

                              -  0.5556 * ((x(index)-3)^5);

      end

end

% Material and geometric properties are given unit value:

E = 1;      % Elastic constant is unit value (not calculated)
I = 1;      % Moment of Inertia is unit value (not calculated)

% Post-processing

figure(1)
plot(x, V, 'r', 'linewidth', 3)
xlabel('x, m')
ylabel('V(x), N')
grid
title('Shear Force Diagram')

figure(2)

plot(x, M, 'b', 'linewidth', 3)
grid
xlabel('x, m')
ylabel('M(x), N-m')
title('Bending Moment Diagram')

figure(3)

plot(x, Y, 'k', 'linewidth', 3)
grid
xlabel('x, m')
ylabel('y(x), m')
title('Elastic Curve Diagram')

disp(' ')

disp('_______________________________________')
disp('Successful execution of Matlab program.')

% end of Matlab program

Figure 1.

 

Figure 2.

 

Figure 3.