Term Project:
Pure Elastic 2D Beam Bending

Using Numerical Methods
Scott M. Steffan
December 9, 1999
Numerical Analysis for Engineering MEAE 4960H02
Rensselaer Polytechnic Institute
–
Hartford
ABSTRACT
The objective of this project involves solving
for deflections along a prismatic, linear elastic bar
subject to pure bending loads. The report will begin with a brief discussion of the Strength of
Materials and Differential Equations theory involved in generating closed form solutions for a
specific
problem. Then, a detailed look into three numerical methods that could be used to
generate approximate solutions will begin
, followed by
numerical
deflection
approximatio
ns
for
two different loading situations
. Finally, using the information provided a discussion on the
advantages and disadvantages of the various numerical techniques will
be completed.
THEORY
The prismatic bar is a beam having constant cross

sectional area in which the length is typically
much greater than any of its cross

sectional dimensions. Figure 1 is an example showing a
prismatic bar bent by equal and opposite coup
les (Shigley & Mischke, 1989):
Figure 1
: Beam Bending
–
Equal and Opposite Couples
There are many different aspects to the study of beam bending, however, this report is done
assuming the following:
(1)
The beam is subject to pure bending
(2)
The material is is
otropic and homogeneous
(3)
The deflections are small
(4)
The bending is pure elastic
Assumption (1) allows axial loading to be neglected. Rather than having two loading types
involved in the stress solution, there is only one:
A
P
I
c
M
Equation 1
: Elastic Stress Equa
tion
Where:
=
the stress at the specified location
M
= the bending moment
c
= the height above or below the beam center line
I
= the beam's second moment of inertia
Assumptions (2)

(4) allow nonlinear effects to be neglected, therefore
Hooke's Law
is
obeyed:
Equation 2
: Hooke’s Law for Stress and Strain
Where:
=
the stress at the specified location
= the strain at the specified location
E = Young's Modulus of the Material (uniform in all directions)
The problem used for this research involves a
simply supported beam subject to a distributed,
uniform load. The question is: How much will a beam deflect under a specified load? The
answer involves the manipulation of the following fourth order differential equation:
Equation 3
: 4
th
Order Differe
ntial Equation

Distributed Load
Where:
E = Young's Modulus of the Material
I
= The beam's second moment of inertia
w
= The load per unit length over the beam (including the weight of the beam)
This is done using the diagram and coordinate system shown i
n figure 2 below:
Figure 2
: Beam Geometry, Boundary Conditions, and Coordinate System
The fourth order differential equation can be solved fairly painlessly when both the load and
cross

section are uniform. First,
boundary conditions
for the spe
cified
end conditions
(simply
supported) are determined:
y(0) = y(L) = 0 = y’’(0) = y’’(L)
Equation 4
: Boundary Conditions
Successive simple integration aids in obtaining the exact solution:
w
dx
y
d
I
E
dx
d
2
2
2
2
4
3
2
2
3
1
4
3
2
2
1
3
2
1
2
2
2
1
3
3
2
6
24
2
6
2
C
x
C
x
C
x
C
x
w
y
I
E
C
x
C
x
C
x
w
x
d
dy
I
E
C
x
C
x
w
x
d
dy
I
E
C
x
w
x
d
dy
I
E
y
x
y
z
w
Equations 5

8
: Differential Equation Integration
Where
C
1
, C
2
, C
3
, C
4
are arbitrary constants. Using
y’’(0)=0
implies
C
2
=0
, and allows for
C
1
=

wL/2
to be determined. Also,
y(0)=0
implies
C
4
=0
, and with
y(L)=0
, it follows that
C
3
=wL
3
/24
.
Using this information, an equation for deflection as a function of distance
along the beam,
y(x)
is
generated:
Equation 9
: Deflection as a Function of Beam Location
Symmetry also reveals that at
x=L/2
, the maximum deflection occurs:
Equation 10
: Max Deflection
This derivation provides a baseline for which numerical methods c
an be compared with.
Equation #9 above is commonly referred to as a
singularity function
. Appendix A has a table of
singularity functions for various different loading and supporting conditions. Differential Equation
boundary conditions for the other ty
pes of end conditions are also supplied in Appendix A.
Once these deflections are determined, elastic stresses and strains can be solved for. As
mentioned earlier, when the cross

sectional dimensions of the beam are constant, as well as the
load over the
surface, the closed form solution can be solved by hand. This is the situation that
will be studied. Further, a more complex situation will be discussed to show the true benefits of
numerical approximation methods.
NUMERICAL METHODS
The 4
th
order beam d
ifferential equation 3 discussed in the previous section is an example of a
boundary value problem
. The problem being discussed involve constant second moment of
inertia (
I(x)=I
0
) and load (
q(x)=q
0
), making them easy to solve analytically and provide soli
d
ground for generating numerical approximations and comparing them to the actual values. There
are many numerical techniques that could be used to generate approximate deflection results.
The three that will be investigated in this research are the
Fini
te Element Method
,
the
Linear
“Shooting” Method
, and the
Finite Difference Method
.
Finite Element Method
The first and probably one of the most widely used techniques today is the
Finite Element Method
(FEM). FEM is a numerical approach commonly used to
solve engineering mechanics problems.
Many general finite element code packages have been written over the years with user friendly
windows and menus (GUI) which allow for easy geometry setup, boundary condition manipulation
and evaluation/post

processing
of common structural problems. Some of the most popular
commercial codes in the industry are ANSYS, MSC Nastran, and MARC. ANSYS will be the
code used for the work in this research.
The FEM subdivides the model, in this case a beam, into smaller sectio
ns (elements) and
determines values at each end (nodes) of these smaller beams. The physical, finite element
mesh is converted into a mathematical model. The differential equations involved in the model
are solved through algebraic matrix manipulation. V
arious element types could be used for this
particular problem. One is known as a 2D

beam element. This type of element can be used for
beam structures that have either a constant cross

sectional area, or can be divided into elements,
varying the cross

s
ectional dimensions appropriately. The beam element type described above
x
L
Lx
x
I
E
w
x
y
3
3
4
2
24
)
(
EI
wL
y
384
5
4
max
will be used for this work. The diagrams in figure 3 represent a thin beam subdivided into 2 and 8,
two

noded beam elements respectively:
Figure 3
: 2

Noded Linear Beam Elem
ent Meshes
This element type allows the user to input each element’s cross

sectional properties, material
properties, and loading/boundary conditions.
Linear Shooting Method
The Linear Shooting Method uses the following 2
nd
order differential equatio
n boundary

value
problem:
Equation 11
: 2
nd
Order Boundary Value Problem
The boundary value problem is replaced by two initial

value problems shown below:
Equations 12

13
: 2
nd
Order Initial Value Problems
Then, using the 4
th
order Runge

Kutta Method, ap
proximate initial value solutions for
y
1
(x)
and
y
2
(x)
are generated with the input of the endpoint values, and the boundary conditions at these
endpoints. Once
y
1
(x)
and
y
2
(x)
are known, the following can be solved:
Equation 14
: Solution for y(x)
Now t
he deflection at any given location along the length of the beam can be approximated,
provided the appropriate number of subintervals is specified. The method is shown graphically in
figure 4 below (Burden & Faires, 1997):
Elements
Nodes
)
(
,
)
(
,
)
(
)
(
'
)
(
'
'
b
y
a
y
b
x
a
x
r
y
x
q
y
x
p
y
1
)
(
'
,
0
)
(
,
)
(
'
)
(
'
'
0
)
(
'
,
)
(
,
)
(
)
(
'
)
(
'
'
a
y
a
y
b
x
a
y
x
q
y
x
p
y
a
y
a
y
b
x
a
x
r
y
x
q
y
x
p
y
)
(
)
(
)
(
)
(
)
(
2
2
1
1
x
y
b
y
b
y
x
y
x
y
Figure 4
: Graphical Repres
entation of “Shooting” Method
Finite Difference Method
The Finite Difference Method is very similar to the Finite Element Method. This method uses the
linear second

order boundary

value problem shown in equation 6 of the previous section. Each
of the de
rivatives in the differential equation are replaced by
difference

quotient
approximations.
Equation 6 would require difference

quotient approximations for
y’’
and
y’
. This is done using the
central

difference formula
shown by equation 7 below:
RESUL
TS AND DISCUSSION
In order to do some comparisons between the actual solution and approximations determined by
the various numerical methods some beam specifications are given:
E = 1e+11
Width = 1
Height = 1
Second Moment of Inertia,
I(x)=I
0
=bh
3
/12
= 1/12
Length = 10
Uniform Load,
q(x)
=
q
0
= 1e+5
Results recorded every 0.25 in length along the beam.
Finite Element Method
Using two

noded 2D

beam elements in ANSYS, a parametric mesh density study was performed.
This allowed us to see how close the true beam
deflection was approximated. Meshes of 2, 4, 5,
8, 10, 20 and 40 elements were used. The linear results were extracted form ANSYS and plotted
in Excel. The lower element meshes were capable of determining deflections exact to within the
computer epsilo
n at the ends of each element (node locations), however, the result between
nodes is very inaccurate. This is because the code is simply linearly interpolating between nodes
to generate results. Figure 6 below shows 2, 4, 5, and 8 element mesh densities a
nd their results
relative to the strength of material solution:
Finite Element Method: 2D Simply Supported Beam
Deflection vs. Axial Beam Location
0.0016
0.0014
0.0012
0.0010
0.0008
0.0006
0.0004
0.0002
0.0000
0.0
1.0
2.0
3.0
4.0
5.0
6.0
7.0
8.0
9.0
10.0
Axial Beam Location (in.)
Deflection (in.)
Y disp  2E
Y disp  4E
Y disp  5E
Y disp  8E
Closed Form
Figure 6
: 2, 4, 5, and 8 Beam Element Meshes
–
Deflection vs. Axial Beam Location
The finite element method first solves for displacements, and then stresses and strains. For this
reason
, it is critical to model the true deflected shape. Using too few elements, results in a
geometry model that is too stiff to accurately predict the peak stresses in the beam. This means
a small error in displacement results in higher order error in the s
tress. Figure 7 shows the 10,
20, and 40 element mesh densities and their results relative to the strength of material solution:
Figure 7
: 10, 20, and 40 Beam Element Meshes
–
Deflection vs. Axial Beam Location
It is visibly clear that the higher number
of elements, the more accurate the deflected shape.
Querying the deflections every 0.25” along the beam allowed for some error calculations. The
plot in figure 8 shows the percent error of the various mesh densities vs. axial beam location:
Finite Element Method: 2D Simply Supported Beam
Deflection vs. Axial Beam Location
0.0016
0.0014
0.0012
0.0010
0.0008
0.0006
0.0004
0.0002
0.0000
0.0
1.0
2.0
3.0
4.0
5.0
6.0
7.0
8.0
9.0
10.0
Axial Beam Location (in.)
Deflection (in.)
Y disp  10E
Y disp  20E
Y disp  40E
Closed Form
Finite Element Method Mesh Density Study
% Error in Deflection vs. Beam Location
0.0
5.0
10.0
15.0
20.0
25.0
30.0
35.0
40.0
0.0
1.0
2.0
3.0
4.0
5.0
6.0
7.0
8.0
9.0
10.0
Beam Location (in.)
% Deflection Error
Error2E
Error4E
Error5E
Error8E
Error10E
Error20E
Error40E
Figure 8
: %
Error vs. Axial Beam Location for Different Mesh Densities
Due to the very small deflections, all error plots were generated in percent error from closed form
solutions. The analysis data obviously supports a finer mesh density, with the 40

element mesh
giving approximations at every node to within the computer epsilon. It must be mentioned,
however, that increased accuracy comes with additional
solution time
per iteration
(possibly more
iterations necessary as well)
. For this simpler problem, the difference was not very noticeable (a
matter of se
conds), but when these engineering mechanics problems, become larger, and more
complex (3D, more involved boundary conditions, etc.), engineering judgement to the accuracy of
the approximation comes into play. It is impractical and sometimes impossible to
check
approximations with the exact solutions. A parametric mesh density study to show a converged
stress solution is being
reported is also a good check.
These studies typically require at least
three mesh densities with the same geometry and boundary
conditions.
Linear “Shooting” Method
The second method used to approximate deflections for the problem was the Linear “Shooting”
Method. Algorithm 11.1 from the textbook was used in Maple V (Burden & Faires, 1997).
The spacing was set to output result
s every 0.25” along the beam. The percent error relative to
the closed form solution is plotted in figure 9 below:
Figure 9
: “Shooting” Method: % Error vs. Axial Beam Location
The solution produces good results with very quick solution time. They are
not as accurate as the
40

element finite element method approximation, but the trade for solution time may make the
“shooting” method a valid approach for this particular problem. The computational error is highest
at the edges, near the boundary conditio
ns. The only part that made this method tough was
generating the second order differential equation. This involved solving part of the integration of
the 4
th
order DE as discussed in the theory section of this report. In cases that vary in cross

section
al dimensions or load over the length of the beam, the integration becomes more complex,
as will be shown in the next section of this report.
Finite Difference Method
Linear "Shooting" Method: % Error in Deflection vs. Beam Location
0.0000
0.0005
0.0010
0.0015
0.0020
0.0025
0.0030
0.0035
0.0040
0.0
1.0
2.0
3.0
4.0
5.0
6.0
7.0
8.0
9.0
10.0
Axial Beam Location (in.)
% Deflection Error
Finally, the Finite Difference Method was used. Algorithm 11.3,also provided by our te
xtbook
(Burden & Faires, 1997) was used in Maple V to generate deflection approximations along the
beam. The method’s accuracy relative to the closed form solutions are shown in the figure 10
below:
Figure 10
: Finite Difference Method: % Error vs. Axial
Beam Location
Again, the results are reasonable, however this method produced the highest overall error. The
most difficult part of this method was to perform the two integrations, resulting in a 2
nd
order DE
for input into the algorithm.
All data us
ed to make the plots in the results and discussion section can be found in Appendix B.
VARYING LOAD PROBLEM
The purpose of adding this section is to show the power of using a numerical analysis technique
to approximate beam

bending deflections. As mentio
ned earlier, not all beams have a constant
load and cross

section making the actual 4
th
order differential equation much more difficult to
evaluate. This is when numerical methods and engineering judgement as to the accuracy of the
solution become very va
luable. To illustrate, solutions for the following beam properties will be
generated using the three methods:
E = 1e+11
Width,
b
= 1
Height,
h
= 1
Second Moment of Inertia,
I
= bh
3
/12
=
1/12
Length = 10
Distributed Load,
w(x)=100000x
Results recorded eve
ry 0.5 in length along the beam.
For this problem, it is more mathematically intensive to generate a closed form solution, because
there are more terms to be integrated due to the varying load. The following set of equations
reflect the actual integratio
n solution:
Finite Difference Method: % Error in Deflection vs. Beam Location
0.0000
0.0100
0.0200
0.0300
0.0400
0.0500
0.0600
0.0700
0.0
1.0
2.0
3.0
4.0
5.0
6.0
7.0
8.0
9.0
10.0
Axial Beam Location (in.)
% Deflection Error
Equations 14

18
: Differential Equation Integration
Using the boundary conditions given in equation 4, the coefficients,
C
1
,
C
2
,
C
3
,
C
4
were
determined and are shown in table 1:
Table 1
: Coefficients for Deflection Equation 18
Also, due
to the fact that the beam

bending problem begins with a 4
th
order differential equation,
one would need to reduce this into a 2
nd
order boundary value problem to make the finite
difference method, or linear “shooting” algorithms supplied by the textbook ge
nerate
approximations. The simple geometry and boundary condition set

up involved in using a
commercial finite element code are displayed in figure 11:
Figure 11
:
ANSYS
Finite Element Method Physical Model
4
3
2
2
3
1
5
3
2
2
1
4
2
1
3
2
2
1
2
3
3
4
4
2
6
120
5
1
2
24
5
1
6
5
1
2
5
1
5
1
C
x
C
x
C
x
C
x
e
y
I
E
C
x
C
x
C
x
e
x
d
dy
I
E
C
x
C
x
e
x
d
dy
I
E
C
x
e
x
d
dy
I
E
x
e
x
d
dy
I
E
0
360
)
5
1
(
7
0
6
5
1
4
4
3
2
2
1
C
L
e
C
C
L
e
C
The results of the four methods and their percent errors are shown in figures
12 and 13 below:
Figure 12
: Deflection vs. Axial Beam Location
Figure 13
: % Error vs. Axial Beam Location
Figure 12 shows the Finite

Difference Method generates the least accurate results using the
same refinement as the other two numerical methods.
The Linear “Shooting” Method produces
deflections closest to the actual, however the Finite Element Method ranks a close second. All
errors were within 0.3% of the closed form, Strength of Materials calculation. The speed and
Simply Supported Beam w/ Linearly Increasing Load
% Error vs. Beam Location
0.0000
0.0500
0.1000
0.1500
0.2000
0.2500
0.3000
0.0
1.0
2.0
3.0
4.0
5.0
6.0
7.0
8.0
9.0
10.0
Axial Beam Location (in.)
% Error
% Error FEM
% Error "Shooting"
% Error FDM
Simply Supported Beam w/ Linearly Increasing Load
Deflection vs. Beam Location
0.0080
0.0070
0.0060
0.0050
0.0040
0.0030
0.0020
0.0010
0.0000
0.0
1.0
2.0
3.0
4.0
5.0
6.0
7.0
8.0
9.0
10.0
Axial Beam Location (in.)
Deflection (in.)
Closed Form
20EFinite Element
Linear "Shooting"
Finite Difference
accuracy of the numerical ap
proximation techniques used to solve for beam

bending deflections
support their use as the degree of difficulty of the problem increases.
All data used to make the plots in this section can be found in Appendix C.
CONCLUSIONS
Although there are many d
ifferent types of numerical methods that could be used to generate
approximations for beam bending boundary value problems, the three chosen showed good
results. All methods require a good mesh spacing to obtain a good representation of the true
beam defl
ection. The Finite Element Method is the most universal, but also the most computation
intensive of the three. When the geometry and/or boundary conditions become more complex,
this method still allows for good results if modeled correctly. The Linear “
Shooting” Method
produced good approximations in this case, however it can often have instability problems. For
this reason, the Finite Difference Method, although it requires more work to obtain the specified
accuracy, was also investigated and performed
well. Some other types of numerical
approximation methods that could be studied for a beam deflection problem are Collocation,
Raleigh

Ritz, Subdomain, and the Least Squared Method to name a few. These additional
studies were beyond the scope of this re
search. Table 1 below was created to rank the three
methods studied as to their accuracy and usefulness:
Table 2
: Method Selection Chart
Although mathematical computation time is the longest for the Finite Element Method, the
versatility and accuracy whe
n used correctly, with a good, refined mesh in areas of concern
makes it the first choice for solving linear elastic beam

bending type problems. Linear “Shooting”
and Finite Difference Methods rank 2
nd
and 3
rd
respectively. This may not be the same for a
ll
boundary

value problems. The ranking is based merely on the research and analysis studies
performed specifically for this project.
REFERENCES
(1)
Shigley, J.E., and Mischke, C.R.,
Mechanical Engineering Design 5
th
Edition
, McGraw

Hill,
Inc., New York, NY
, 1989
(2)
Burden, R., and Faires, J.D.,
Numerical Analysis 6
th
Edition
, Brooks/Cole Publishing
Company, Pacific Grove, CA, 1997
(3)
Timoshenko, S.P., Goodier, J.N.,
Theory of Elasticity 3
rd
Edition
, McGraw

Hill, Inc., New
York, NY, 1970
(4)
Edwards, C.H. Jr., and Pen
ney, D.E.,
Elementary Differential Equations With Boundary Value
Problems 3
rd
Edition
, Prentice

Hall, Inc., Englewood Cliffs, New Jersey, 1993
(5)
Bathe
, K.J.,
and Wilson
, E.L.
,
Numerical Methods In Finite E
lement Analysis
,
Prentice

Hall,
Inc., En
glewood Cliffs, New Jersey, 1976
(6)
Desai, C.S.,
Elementary Finite Element Method
,
Prentice

Hall, Inc., En
glewood Cliffs, New
Jersey, 1979
(7)
Craig, R.R. Jr.,
Mechanics of Materials
, John Wiley & Sons, Inc., New York, NY, 1996
Criteria
Finite Element
Linear "Shooting"
Finite Difference
Ease Of Use
1
2
2
Solve Time
3
1
2
Accuracy of Solution
1
1
2
Versatility
1
3
2
Availability
1
1
1
Totals
7
8
9
Comments 0
Log in to post a comment