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 modied

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

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

n1

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

xM

0

P < xx

P

>

W

2

(< xx

ws

>

2

< xx

we

>

2

)+F

m

< xx

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 dened 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 buer zone to allow for environmental variables and material imperfections.To

accomplish this,a safety factor is specied.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 modied 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 eective solution,the smallest beamwas found for both

high and low strength steels that satises 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 eective 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 modied beam design.The bending shear function now lies

inside the safety bounds indicated by the dashed lines.

3.2.1 Code Verication

A sanity check can be performed by inspecting the graphs to see if initial conditions are

satised.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 conrm 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 dierence 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 dierence 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 dierence in these two methods is on the order of 5%,

which is well within reasonable tolerance for most applications.The safety factor specied

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

## Σχόλια 0

Συνδεθείτε για να κοινοποιήσετε σχόλιο