Finite element analysis by anakidris

VIEWS: 70 PAGES: 44

• pg 1
```									One-Dimensional Problems
We wish to use FEM for solving the following problems:

x

10 kN            x

δ = 2 x 10-2 mm
Objectives
1. To develop a system of linear equations for one-dimensional
problem.
2. To apply FE method for solving general problems involving
bar structures with different support conditions.

Consider a non-uniform bar subjected
shown.

Note: The bar is constrained by a fix support at
the top and is free at the other end. The positive
x-direction is taken downward.

a) Body force, f
Distributed force per unit volume (N/m3)
Example: self-weight due to gravity

b) Traction force, T
Force per unit area (N/m2)
For a 1-D problem,
 force 
T =        × (perimeter of area )
 area 
Examples: Frictional forces, Viscous drag,
and Surface shear.
Concentrated load (in Newton) acting at any point i.
Finite element Modeling
Element Discretization
The first step is to subdivide the bar into several sections – a
process called discretization.

Note: The bar is discretized
into 4 sections, each has a
uniform cross-sectional area.
The non-uniform bar is
transformed into a stepped
bar.
We will use the stepped bar
as a basis for developing a
finite element model of the
original non-uniform bar.
Numbering Scheme
To analyze the stepped bar systematically, a global numbering
scheme is assigned as shown. The x-direction is considered as
the global coordinate direction.
Note:
F1, …, F5 represent global
forces acting on the points
connecting all sections of
the stepped bar.
Q1, …, Q5 represent global
displacements of the points
resulting from the forces
acting on these points.                     {Q} = [Q1    Q2 Q3 Q4 Q5 ]
T

The stepped bar is transformed into         {F } = [F1   F2   F3   F4   F5 ]
T

a finite element model using 1-D
(line) elements.
Element Connectivity
Consider a single line element. It lies in a local coordinate
system, denoted by^  x.

Note: Node number in local              Element connectivity table
coordinate is denoted by a
number with a hat on top.
ˆ
1   ˆ
2
^               ^
ˆ
x

q1 and q2 are nodal displacements
in the local coordinate direction.

Connectivity between global and local nodes must be established for
each element, as tabulated in the table shown.
Natural Coordinate and Shape Functions
Natural Coordinate
Consider a single element. Local node 1 is at distance x1 from a
datum, and node 2 is at x2, measured from the same datum point.

We define a natural or intrinsic coordinate system, ξ,

ξ=
2
(x − x1 ) − 1
(x2 − x1 )
Note: The ξ-coordinate will be used to define shape functions,
required to establish interpolation function for the displacement
field within the element.
Shape Functions
The displacement field, u(x), within the element is not known.
For simplicity, it is assumed that the displacement varies
linearly from node 1 to node 2 within the element.

We establish a linear interpolation function to represent the linear
displacement field within the element. To implement this, linear shape
functions are defined, given by,
1− ?                    1+ ?
N1 (? ) =          and N 2 (? ) =
2                       2
The linear displacement field, u(x), within the element can now
be expressed in terms of the linear shape functions and the
local nodal displacement q1 and q2 as:
^
u ( x ) = N1q1 + N 2 q2
^
1− ξ       1+ξ 
u( x) =       q1 +      q2
 2          2 
In matrix form:
u( x ) = [N ]{q}
^

where                [N ] = [N1   N2 ]
 q1 
and                  {q} =   = [q1     q2 ]
T

 q2 
Isoparametric Formulation
Coordinate x of any point on the element (measured from the
same datum point as x1 and x2) can be expressed in terms of
the same shape functions, N1 and N2 as

x = N1 x1 + N 2 x2
1−ξ       1+ξ 
x=     x1 +      x2
 2         2 

When the same shape functions N1 and N2 are used to establish
interpolation function for coordinate of a point within an element and
the displacement of that point, the formulation is specifically referred
to as an isoparametric formulation.
Example 1

(a) Evaluate ξ, N1, and N2 at point P.
(b) If q1 = 0.003 in and q2 = -0.005 in, determine the value of
displacement u at point P.
Solution
(a) The ξ coordinate of point P is given by
2
ξP =         ( x − x1 ) − 1
x2 − x1
2
=         ( 24 − 20 ) − 1
36 − 20
ξ P = −0.5
The shape functions are:

1 − ξ 1 + 0.5
N1 =        =        = 0.75
2      2
1 + ξ 1 − 0.5
N2 =      =        = 0.25
2      2

(b) Displacement of point P

uP = N1q1 + N 2 q2
= 0.75 ( 0.003) + 0.25 ( −0.005 )
uP = 0.001 in
Strain-Displacement Relation
Normal strain is related to displacement by
du
ε=
dx
Using the chain rule of differentiation

du d ξ
ε=
d ξ dx
The two terms of the above relation are obtained as follows
2                            dξ      2
ξ=             ( x − x1 ) − 1 ⇒         =
( x2 − x1 )                       dx ( x2 − x1 )

 1− ξ          1+ ξ            du ( − q1 + q2 )
u =         q1 +        q2   ⇒       =
 2             2               dξ       2
Thus the normal strain relation can be written as
1
ε=             [ −q1 + q2 ]
( x2 − x1 )
which can be written in matrix form as

 q1 
ε = [B ] 
 q2 
where [B] is a row matrix called the strain-displacement matrix,
given by
1               1
[ B] =               [ −1 1] = [ −1 1]
( x2 − x1 )          le

since x2 – x1 = element length = le.
Stress-Strain Relation
Normal stress is related to the normal strain by a Hooke’s
law,
σ = Eε
where E is modulus of elasticity.

Substitute for the normal strain ε,
we get,

 q1 
σ = E [ B]  
 q2 

Robert Hooke (1635-1703);
(Experimental Philosopher)
Element Stiffness Matrix
We will use the potential energy approach to derive the
element stiffness matrix [ k] for the 1-D element.
Total potential energy of a body subjected to loads is given by,

πp =U +Ω              U = internal strain energy;
Ω = potential energy of external forces.

For the non-uniform bar, its total potential energy is given by
1
σ T ε A dx − ∫ u T fA dx − ∫ uT T dx − ∑ Qi Pi
2 ∫L
πp =
L             L
i

Since the bar has been discretized into finite elements
1
πp = ∑          σ T ε A dx − ∑ ∫ uT fA dx − ∑ ∫ uT T dx − ∑ Qi Pi
e   2 ∫e              e
e
e
e
i
We will derive the element stiffness matrix of the 1-D element
using the internal strain energy term, U as follows,
1
2 ∫e
Ue =        σ T ε A dx

Recall, the stress and strain are given by

σ = E [ B ] {q}    and         ε = [ B ]{q}

Substitute these into the expression for Ue,

U e = ∫ E ([ B ] {q}) [ B ]{q} A dx
1                T

2 e
1
= ∫ {q} [ B ] E [ B ]{q} A dx
T      T

2 e
1
2
(
U e = {q} ∫ [ B ] E [ B ]A dx {q}
T
e
T
)
dξ 2                le
Recall again,      =       ⇒ dx =       dξ
dx le               2
Substitute and simplifying the expression yields,

1   T                   le +1 
U e = {q} [ B ] Ee [ B ] Ae ∫ dξ  {q}
T

2                       2 −1 
1                      le
= {q} [ B ] Ee [ B ] Ae ( 2 ){q}
T     T

2                      2
1             1  −1 1
= {q} Aele Ee   [ −1 1]{q}
T

2             le  1  le
1            1 −1
= {q} Aele Ee 2   [−1 1]{q}
T

2           le  1 
1   T A E  1    −1
U e = {q} e e           {q}
2      le  −1 1 
The internal strain energy for the 1-D element can now be
written in the form,
1
U e = {q} [ k ] {q}
T     e

2
where [k]e represents the element stiffness matrix for the 1-D
element, i.e.
Ee Ae  1 −1
[k ] =
e

le  −1 1 
     
Note: Ee = elastic modulus;
Ae = cross-sectional area;
le = element length.
Element Force Vector
The forces acting on 1-D structures can be of body force, fb,
traction force, T, and concentrated force, P. They may act
individually in various combination.
The total potential energy of the structure,
1
π p = ∑ ∫ σ T ε A dx − ∑ ∫ uT f b A dx − ∑ ∫ uTT dx − ∑ Qi Pi
e 2
e               e                 e
e                 e            i

a) Due to body force, fb
The potential energy due to body force fb in a single element is
given by the second term, i.e.

Ω f = ∫ uT f b A dx
e

= ∫ ( N1q1 + N 2 q2 ) f b A dx
T
e

Ω f = Ae fb ∫ ( N1q1 + N 2 q2 ) dx
T

e
 Ae fb N1 dx 
Ω f = {q}
T           ∫e        
Rewrite,                                       
 Ae fb ∫e N 2 dx 
                 
Recall that,
le
dx =      dξ
2
le +1 1 − ξ  le
Also,          ∫e N1 dx = 2 ∫−1 2 dξ = 2
le +1 1 + ξ  le
∫e N 2 dx = 2 ∫−1 2 dξ = 2
Substitute and simplifying, yields

 Ae f bl e 
 2 
                 T Al f     1
Ω f = {q}                    = {q} ⋅ e e b
T
                            
 Ae f bl e           2      1
 2 
           
The potential energy due to the body force can now be
expressed in the form,

Ω f = {q}      {f}
T         e

where the force vector due to body force fb is,

Ae le fb       1
{f}      =
e

2           1

Quiz: Can you give the physical interpretation of {f}e?
b) Due to traction force, T

The potential energy due to traction force T is given by,

ΩT = ∫ uT T dx = ∫ ( N1q1 + N 2 q2 ) T dx
T
e             e

Recall,                                   le +1 1 − ξ  le
le
dx = dξ           ∫e N1 dx = 2 ∫−1 2 dξ = 2
2
le +1 1 + ξ  le
∫e N 2 dx = 2 ∫−1 2 dξ = 2
Rearranging and simplifying,
 le 
T N1 dx 
 ∫e                   2 
 
ΩT = {q}                     = {q}
T                       T
                       
T ∫e N 2 dx 
                       le 
2 
 
The last equation is in the form,

ΩT = {q} {T }
T       e

le 1
ΩT = {q}
T
i.e.                               ⋅T ⋅  
2 1

Thus, element traction force vector due to traction T,

Tle    1
{T }e
=        
2     1
Quiz: Can you give the physical interpretation of this?
Summary

We have established, for 1-D problems,

1. Stress-strain relation         3. Element force vector
due to body force, fb
E         q1 
σ = [ −1 1]  
q2                          Al f      1
le
{f}
e
= ee b     
2       1
2. Element stiffness matrix
4. Element force vector
AE      1 −1
[k ]                              due to traction force, T
e
= e e    −1 1 
le                                  Tl   1
{T }     = e
e

2   1
Example 2

A thin steel plate has a uniform
thickness t = 1 in., as shown. Its
elastic modulus, E = 30 x 106
psi, and weight density, ρ =
0.2836 lb/in3.
The plate is subjected to a point
load P = 100 lb at its midpoint and
a traction force T = 36 lb/ft.
Determine:
a) Displacements at the mid-point
and at the free end,
b) Normal stresses in the plate, and
c) Reaction force at the support.
Solution
1. Transform the given plate into 2 sections, each having
uniform cross-sectional area.
Note:
Area at midpoint is
Amid = 4.5 in2.
Average area of section 1 is
A1 = (6 + 4.5)/2 = 5.25 in2.
Average area of section 2 is
A2 = (4.5 + 3)/2 = 3.73 in2.

2. Model each section using 1-D
(line) element.
3. Write the element stiffness matrix for each element

5.25 × 30 ×106    1 −1
[k ]     =
(1)
element 1:                                   −1 1 
12               

3.75 × 30 ×106    1 −1
[k ]     =
(2)
element 2:                                   −1 1 
12               

4. Assemble global stiffness matrix,

 5.25 −5.25   0 
30 × 10 
[ K ] = 12 −5.25 9.00 −3.75
6


 0
      −3.75 3.75 


Note: The main diagonal must contain positive numbers only!
5. Write the element force vector for each element

a) Due to body force, fb = 0.2836 lb/in3

5.25 × 12 × 0.2836 1
{ fb }
(1)
element 1                        =                    
2         1
3.75 × 12 × 0.2836 1
{ fb }
(2)
element 2                        =                    
2         1
Assemble global force vector due to body force,

5.25  8.9 
12 × 0.2836            
{ Fb } =             9.00 = 15.3
2     3.75  6.4 
           
b) Due to traction force, T = 36 lb/ft

 36 
  ×12 1      1
 12 
{T }
(1)
element 1                    =           = 18  
2  1      1

 36 
  × 12 1      1
 12 
{T }
(2)
element 2                    =            = 18  
2   1      1

Assemble global force vector due to traction force,

1  18 
   
{ FT } = 18 2 = 36 
1  18 
   
c) Due to concentrated load, P = 100 lb at node 2

 0 
 
{ FP } = 100
 0 
 
6. Assemble all element force vectors to form the global force
vector for the entire structure.

 8.9 + 18 + 0   26.9 
                        
{ F } = 15.3 + 36 + 100  = 151.3 lb
 6.4 + 18 + 0   24.4 
                        
7. Write system of linear equations (SLEs) for entire model

The SLEs can be written in condensed matrix form as

[ K ]{Q} = {F }
Expanding all terms and substituting values, we get

 5.25 −5.25    0   Q1   26.9 
30 ×10                                
−5.25 9.00 −3.75  Q2  = 151.3
6

12                     
           
 0
       −3.75 3.75  Q3   24.4 

Note:
1. The global force term includes the unknown reaction force R1 at
the support. But it is ignored for now.
2. The SLEs have no solutions since the determinant of [K] = 0;
Physically, the structure moves around as a rigid body.
8. Impose boundary conditions (BCs) on the global SLEs

There are 2 types of BCs:
a) Homogeneous = specified zero displacement;
b) Non-homogeneous = specified non-zero displacement.

In this example, homogeneous BC exists at node 1. How to
impose this BC on the global SLEs?

DELETE ROW AND COLUMN #1 OF THE SLEs!

 5.25 −5.25   0   Q1   26.9 
30 × 10 
6
          
−5.25 9.00 −3.75 Q2  = 151.3
12                      
          
 0
       −3.75 3.75  Q3   24.4 

9. Solve the reduced SLEs for the unknown nodal
displacements
The reduced SLEs are,

30 ×106    9.00 −3.75 Q2  151.3
 −3.75 3.75  Q  =  24.4 
12                   3          

Solve using Gaussian elimination method, yields

Q2  1.339 ×10−5 
 =           −5 
in
Q3  1.599 ×10 

Quiz: Does the answers make sense? Explain…
10. Estimate stresses in each elements

1        q1 
Recall,         σ   (e )
= E [ B ]{q} = E ⋅ [−1 1]  
le       q2 
element 1

1            0      
σ   (1)
= 30 ×10 ⋅ [ −1 1] 
6
−5 
= 33.48 psi
12       1.339 ×10 

element 2
1       1.339 × 10 −5 
σ   (2)
= 30 × 10 ⋅ [ −1 1] 
6
−5 
= 6.5 psi
12       1.599 × 10 
11. Compute the reaction force R1 at node 1

We now include the reaction force term in the global SLEs.
From the 1st. equation we get,

 5.25 −5.25   0         0        26.9 + R1 
30 × 10 
6
 1.339 × 10−5  =  151.3 
12     −5.25 9.00 −3.75                          
 0
       −3.75 3.75  1.599 × 10−5   24.4 
                         

We have,
     0       
30 × 10
6
             
R1 =         [ 5.25 −5.25 0 ] 1.339 × 10−5  − 26.9334
12                     1.599 × 10−5 
             
R1 = −202.68 lb
Example 3
A concentrated load P = 60 kN is
applied at the midpoint of a uniform
bar as shown.
Initially, a gap of 1.2 mm exists
between the right end of the bar
and the support there.
1.2 mm
If the elastic modulus E = 20 x 103
N/mm2, determine the:                  250 mm2
a) displacements field,                           P

b) stresses in the bar, and                                        x

c) reaction force at the support.
150 mm    150 mm
Solution
1. Write the element stiffness matrices and assemble the
global stiffness matrix.
 1 −1 0 
20 × 10 × 250 
3

[ K] =                −1 2 −1

150
 0 −1 1 
        
2. Write the element force vectors and assemble the global force
vector.
{F } = 0,
T
     60 ×10 , 0 
3

3. Write the global system of linear equations.

 500 −500   0   Q1         0
10                              
−500 1000 −500  Q2  = 103 60 
3

15                 
 0
      −500 500   Q3 
 
0
 
4. Impose the boundary conditions.
We have; Q1 = 0; Q3 = 1.2 mm. Using Gaussian elimination
method:
a) Delete 1 st row and column.
b) Delete 3 rd row and column and modify the force term.

 500 −500   0   Q1          0
10                               
−500 1000 −500   Q2  = 103 60 
3

15                 
 0
      −500 500  1.2 
 
0
 
The reduced SLE becomes,                          Modification to
force term
103                       500 (1.2 ) 
[1000]{Q2 } = 10 60 +
3

15                          15 
7. Solve the reduced SLE, we get

Q2 = 1.5 mm
8. Compute stresses in the bar,
1          0
σ 1 = 20 × 10 ×
3
[ −1 1] 1.5
150          
σ 1 = 200 MPa
1          1.5
σ 2 = 20 × 10 ×
3
[ −1 1] 1.2 
150          
σ 2 = −40 MPa

9. Compute reaction forces at supports
Using the 1 st and 3rd equations, we obtain,
R1 = -50 x 103 N;   R3 = -10 x 103 N.
Exercise 1
A composite bar ABC is subjected to axial forces as shown.
Given, the elastic moduli, E1 = 200 GPa and E2 = 70 GPa.
Estimate:
a) Displacement of end C; [Answer: δC = 6.62x10-2 mm]
b) Stress in section 2, and
c) Reaction force at support A.
Verify your results with analytical solution.
Exercise 2
Reconsider Exercise 2-1. Suppose a gap of δ = 2 x 10-2 mm
exists between end C and a fixed support there. Estimate:
a) Displacement of point B;
b) Stress in section 1, and
c) Reaction forces at both supports.

10 kN

2 x 10-2 mm
Assignment 1

q Find a journal paper on the application of finite element method
to model and simulate real world problems, from various journ-
als on the internet.
(e.g. : www.sciencedirect.com).