# f2 - Rensselaer at Hartford - Rensselaer Polytechnic Institute

Urban and Civil

Nov 15, 2013 (4 years and 8 months ago)

155 views

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
. 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
-

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

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

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

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
Error-2E
Error-4E
Error-5E
Error-8E
Error-10E
Error-20E
Error-40E
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.

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
-
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

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
20E-Finite 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