# Comparing FEM and Analytical Methods for Solving a Beam Loading Design Problem

Πολεοδομικά Έργα

15 Νοε 2013 (πριν από 4 χρόνια και 5 μήνες)

99 εμφανίσεις

ME 128,Project 1
Professor Liwei Lin
UC Berkeley
Comparing FEM and Analytical
Methods for Solving a Beam Loading
Design Problem
Yue (Amy) Fu
October 6,2010
1 Introduction
There are two ways that beam bending problems are typically solved:analytically using
statics,and computationally using the Finite Element Method (FEM) and a computer sim-
ulation.There are inaccuracies associated in both methods.In this project,the same beam
calculations were performed analytically and using FEMand the results were compared and
analyzed.These tools were then used to design a beam that would not fail under the given
loads.
2 Theory
2.1 Problem Statement
The design problemconsists of a steel wide ange I-Beamwith a xed support on one end and
a simple support on the other.The beam is loaded with a distributed load,a concentrated
load,and a concentrated moment as shown in Figure 1.
Figure 1:The beam and external loads used in the design problem.
The following variables are given,and all locations are measured from the xed end.
 Beam Length L,120 in
 Cross-Sectional Area Moment of Inertia I,16.4 in
4
 Distributed Load Intensity W,180 lb/in
 Location Where Distributed Load Begins x
ws
,45 in
 Location Where Distributed Load Ends x
we
,100 in
 Concentrated Load Magnitude P,475 lb
1
 Concentrated Load Location x
P
,30 in
 Concentrated Moment Magnitude F
m
,1800 in-lb clockwise
 Concentrated Moment Location x
M
,30 in
 Young's Modulus E,29  10
6
psi
 Beam Height d,5.9 in
 Yield Strength 
y
,36,000 psi
The following functions need to be computed and graphed as a function of x.
 Shear Force,V (x)
 Bending Moment,M(x)
 Rotation Angle,(x)
 Neutral Axis Displacement,y(x)
 Bending Stress,(x)
The von Mises failure criteria is then applied,and the beam is analyzed for failure based on
a given safety value.If the beam does not meet the safety criteria,the design is modied
until it does.
2.2 Analytical Solution
2.2.1 Free Body Diagram
The rst step to solve this problem is to draw the free body diagram of the system.
Figure 2:Free body diagram of the beam.
The equations of equilibrium are then written based on the free body diagram.
2
X
F
horizontal
= 0:A
x
= 0 (1)
X
F
vertical
= 0:A
y
+F
y
= P +W(x
we
x
ws
) (2)
X
M
A
= 0:M
0
+F
y
L = Px
P
+W(x
we
x
ws
)x
W
(3)
where x
W
= x
ws
+
x
we
x
ws
2
(4)
Since there are four external reaction forces,the system is indeterminate.Bending equations
must be used in order to calculate the value of all four reaction forces.
2.2.2 The Bending Equations
The relationship between the internal moment function M(x) and the neutral axis displace-
ment (bending) function y(x) is described in Equation 5.
EI
d
2
y(x)
dx
2
= M(x) (5)
Equation 5 can be integrated as follows.
d
2
y(x)
dx
2
=
1
EI
M(x) (6)
dy(x)
dx
= (x) =
1
EI
Z
M(x)dx +C
1
(7)
y(x) =
Z
(x)dx +C
2
(8)
The booundary condition for Equation 7 is as follows because there is no rotation at a xed
end.
(0) = 0 (9)
The booundary conditions for Equation 8 are as follows because there is no displacement at
a xed or simple end.
y(0) = 0 (10)
y(L) = 0 (11)
From the bending equations,there appears two additional unknowns C
1
and C
2
,but also
three equations,one for each boundary condition.This gives a grand total of 5 equations
and 5 unknowns,which is enough to solve for all the external reaction forces on the beam.
3
2.2.3 Determining M(x)
M(x),the internal moment in the beam,can either be computed analytically using the
method of sections or by integrating V (x),the internal shear force in the beam.
V (x) =
dM(x)
dx
(12)
V (x) can be determined using the method of sections and superposition.Superposition is
the principle in which V (x) for each of the three loads are analyzed independently and the
results are added together to determine the nal function.Using the method of sections
is to slice the beam at strategic locations and determining the internal forces at the slice,
generating a piecewise function of the internal loads.
2.2.4 Computing V (x) Using Superposition and Sections
Figure 3:Free body diagram of the beam with force P alone.
Figure 4:Free body diagram of the beam with force P alone,sectioned at 0  x < x
P
.
Since there are no horizontal forces on the segment (A
y
= 0),only the equation for force
equilibrium in the x direction (
P
F = 0) is necessary to calculate V
p
(x).
V
p
(0  x < x
P
) = A
y;P
(13)
4
Figure 5:Free body diagram of the beam with force P alone,sectioned at x
P
 x < L.
V
p
(x
P
 x < x
L
) = A
y;P
P (14)
Figure 6:Free body diagram of the beam with moment Fm alone.
Since there are no vertical forces throughout the beamin Figure 6,V is constant throughout.
V
m
(0  x < L) = A
y;m
(15)
Figure 7:Free body diagram of the beam with distributed load alone.
5
Figure 8:Free body diagram of the beam with distributed load alone,sectioned at 0  x < x
ws
.
V
W
(0  x < x
ws
) = A
y;W
(16)
Figure 9:Free body diagram of the beam with distributed load alone,sectioned at x
ws
 x < x
we
.
V
W
(x
ws
 x < x
we
) = A
y;W
W(x x
ws
) (17)
Figure 10:Free body diagram of the beam with distributed load alone,sectioned at x
we
 x < L.
V
W
(x
we
 x < L) = A
y;W
W(x
we
x
ws
) (18)
6
2.2.5 The Singularity Function
Since piecewise functions are a hassle to integrate,a singularity function (shown below) is
dened to aid in computation.
< x x
0
>
n
=
(
0 if x < x
0
(x x
0
)
n
if x  x
0
(19)
The singularity function can be integrated and dierentiated like normal:
Z
< x a >
n
dx =
1
n +1
< x a >
n+1
+C for n  0 (20)
d
dx
< x a >
n
= n < x a >
n1
for n  1 (21)
Using the singularity function,the shear functions V (x) computed in Section 2.2.4 can be
rewritten as follows:
V
P
(x) = A
y;P
P < x x
P
>
0
(22)
V
m
(x) = A
y;m
(23)
V
W
(x) = A
y;W
W(< x x
ws
>  < x x
we
>) (24)
Using superposition and adding these equations,we obtain:
V (x) = A
y
P < x x
P
>
0
W(< x x
ws
>  < x x
we
>) (25)
Integrating Equation 25 as per Equation 12 using the boundary condition M(0) = M
0
and
taking into account the concentrated moment in Figure 6,we get:
M(x) = A
y
xM
0
P < xx
P
> 
W
2
(< xx
ws
>
2
 < xx
we
>
2
)+F
m
< xx
m
>
0
(26)
Using Equation 7 and 8 and the boundary conditions in Equations 9 and 10,we can nd the
rotation angle and neutral axis displacement functions respectively:
(x) =
1
EI
(
1
2
A
y
x
2
M
0
x 
1
2
P < x x
P
>
2
(27)

1
6
W(< x x
ws
>
3
 < x x
we
>
3
) +F
m
< x x
m
>)
y(x) =
1
EI
(
1
6
A
y
x
3
M
0
x 
1
6
P < x x
P
>
3
(28)

1
24
W(< x x
ws
>
4
 < x x
we
>
4
) +
1
2
F
m
< x x
m
>
2
)
Lastly,maximum bending stress (x) is dened as
Mc
I
,where c is the maximum vertical
distance from the centroid to the outside edge of the beam,and I is the area moment inertia
of the beam cross section.
7
2.3 Failure Criteria
The von Mises failure criteria is applied to the maximum bending stress (x) to determine
whether the beam will fail.Since this problem deals with uniaxial stress,the von Mises yield
criteria is reduced to:
j(x)j = 
y
(29)
This means when the maximum bending stress reaches the yield stress 
y
,the beam will fail.
However,this is not a safe criteria to use when designing a beam,since it is highly desirable
to have a buer zone to allow for environmental variables and material imperfections.To
accomplish this,a safety factor is specied.The rejection criteria used is:
j(x)j =

y
Safety Factor
(30)
In the case when this criteria is exceeded,the beam needs to be redesigned either with
stronger materials or larger cross sections while keeping the cost down to a minimum.
2.4 Analytical Solution Using MATLAB
R
A semi-automatic way to solve this design problem was made using MATLAB
R

.A function
was written that takes in the inputs listen in Section 2.1 plus a safety factor and graphs the
output functions,also listed in Section 2.1.
The function uses matrices to solve Equations 2,3,and 11 simultaneously for A
y
,M
0
,and
F
y
.A singularity function was written separately and referenced by the function.It then
uses 1000 equally spaced points between 0 and L and computes the output functions at those
points.These functions are then graphed in a subplot layout.The graph for (x) includes
a dashed line indicating the failure criteria stated in Equation 30,so it is clearly visible
whether or not the beam meets this criteria.In addition,the function outputs a sentence
stating whether or not the beam meets the safety criteria.
The results of this function are printed in Section 3.2.The complete source code is attached
in Section 5.
8
2.5 FEM Using SolidWorks
The cross section sketch shown in Figure 11 was extruded 45 in,55 in,and 20 in to form
a 120 inch long beam.This was done to be able to simulate the distributed load on the
second segment.A point was sketched at 30 in as a location for the concentrated load
and moment.The material used was Cast Carbon Steel with its yield strength and Young's
modulus slightly modied to match the problemstatement.The part features and simulated
loads are visible in Figure 12.
Figure 11:SolidWorks sketch of the beam cross section.Dimensions are per original beam design.
Fillet radius was determined arbitrarily.
9
A xed geometry boundary condition was used at one end of the beam,and two translational
degrees of freedom were constrained on the other end to simulate a simple support.Three
external loads were applied to the part:one distributed load on the second extrude,one
concentrated force on Point 1,and one clockwise moment at Point 1.A mesh was generated
that divided up the beam into around 100 segments,and the simulation was run.The same
output plots were generated for analysis.
Figure 12:Feature tree and applied loads of the beam in SolidWorks Simulator.
10
3 Results
3.1 Summary of Results
Both the analytical method and the nite element method results show that the current
beam does not satisfy the required safety conditions,as shown in Figure 13.The maximum
bending stress exceeds the yield strength of the low strengh steel at a safety factor of 2.
Figure 13:The current beam design does not satisfy the safety criteria,as shown by the red
segments in the FEM graph and the portions of the maximum bending stress curve that lies
outside the dashed lines.
This issue can be xed by either using a bigger beam or by using high strength steel ( =
52;000psi).However,high strength steel costs \$2.85/lb as opposed to \$1.25/lb for low
strength steel.To nd the most cost eective solution,the smallest beamwas found for both
high and low strength steels that satises the safety criteria and their costs were compared.
The analytical method was used for this process.These results are summarized in Table 1.
Steel Strength
Beam Size
Cost
Safety Criteria Met?
Low Strength
W6x9
\$113.36
No
Low Strength
W6x15
\$189.81
Yes
High Strength
W8x15
\$427.25
Yes
Table 1:Summary of beam strengths and costs.Density of steel = 0.282 lb/in
3
From these results,it can be concluded that the best solution would be to use a W6x15 low
strength steel beam,costing \$189.81.
11
3.2 MATLAB
R

Results
The function was run with the original inputs,and again with the most cost eective safe
beam design detailed in Section 3.1.The output plots are shown in Figures 14 and 15.
Figure 14:Output functions for the original beam design.
12
Figure 15:Output functions for the modied beam design.The bending shear function now lies
inside the safety bounds indicated by the dashed lines.
3.2.1 Code Verication
A sanity check can be performed by inspecting the graphs to see if initial conditions are
satised.In the above two gures,the rotation angle at x = 0 is zero,and displacement at
both ends of the beam is also zero.Therefore,the code is presumably correct.
13
3.3 FEM Results
Figure 16:FEMplots for the original beam design,using the same layout as the MATLAB
R
plots
in Figures 14 and 15 with a safety factor plot in place of a bending stress plot for clarity.
Graphically,the analytical and FEM results appear to be very similar.To conrm these
observations,I plotted both data sets in the same graph for a side by side comparison in
Figure 17.
14
Figure 17:FEM and analytical results are very similar.
Futhermore,I graphed the absolute value of the dierence in the two plots,shown in Figure
18.The error is constant for areas with constant shear,but error is accumulated when shear
is changing as a result of the distributed load.However,the error never exceeds more than
5% of the shear.
15
Figure 18:Absolute value of the dierence between FEM and analytical results.
4 Conclusion
The more accurate result is generated by nite element analysis because the bending equa-
tions are derived assuming very small deformations,while the nite element method does
not use these assumptions.This,however,does not make the analytical method any less
useful than FEM.It is a lot faster to run code based o of the beam equations than it is to
run FEMfor the same resolution.The dierence in these two methods is on the order of 5%,
which is well within reasonable tolerance for most applications.The safety factor specied
allows for some variation in terms of which calculation method was used.
I chose to use the analytical method for the design problem presented because it was a lot
faster and more convenient to change a few numbers in a function than it is to change a
SolidWorks model and run the simulation for every beamsize tested.However,high precision
applications such as MEMS systems require a much more accurate failure model,and FEM
simulations must be used.The analysis method used depends on application;it is up to the
engineer to decide which is best on a case by case basis.
16
5 Source Code
The MATLAB
R
functions used are attached to the end of this report.
17