# The Finite Element Method AP ractical Course by pYJ3t1

VIEWS: 68 PAGES: 56

• pg 1
```									The Finite Element Method
A Practical Course

CHAPTER 5:

FEM FOR 2D SOLIDS
CONTENTS
   INTRODUCTION
   LINEAR TRIANGULAR ELEMENTS
–   Field variable interpolation
–   Shape functions construction
–   Using area coordinates
–   Strain matrix
–   Element matrices
   LINEAR RECTANGULAR ELEMENTS
–   Shape functions construction
–   Strain matrix
–   Element matrices
–   Gauss integration
–   Evaluation of me
CONTENTS
– Coordinate mapping
– Strain matrix
– Element matrices
– Remarks
   HIGHER ORDER ELEMENTS
   CASE STUDY
INTRODUCTION
 2D solid elements are applicable for the analysis
of plane strain and plane stress problems.
 A 2D solid element can have a triangular,
rectangular or quadrilateral shape with straight or
curved edges.
 2D solid element can deform only in the plane of
the 2D solid.
 At any point, there are two components in x and y
directions for the displacement as well as forces.
INTRODUCTION
   For plane strain problems, the thickness of the
element is unit, but for plane stress problems, the
actual thickness must be used.
   In this course, it is assumed that the element has a
uniform thickness h.
   Formulating 2-D elements with given variation of
thickness is also straightforward, as the procedure
is the same as that for a uniform element.
2D solids – plane stress and plane strain

Plane stress                Plane strain
LINEAR TRIANGULAR ELEMENTS
 Less accurate than quadrilateral elements
 Used by most mesh generators for complex
geometry
 Linear triangular element

Nodes   Triangular elements
y, v                3 (x3, y3)
(u3, v3)
fsy
fsx
A
2 (x2, y2)
(u2, v2)
1 (x1, y1)
(u1, v1)                            x, u
Field variable interpolation
U h ( x, y )  N( x, y )de
u1         
v           displaceme nts at node 1
 1         
u 
           
where    de   2          displaceme nts at node 2
v 2        
u 3        
            displaceme nts at node 3
 v3 
           
y, v                3 (x3, y3)
 N1    0     N2     0       N3    0                                  (u3, v3)
N
0      N1    0      N2      0     N3 

fsy
fsx
Node 1      Node 2          Node 3                                A
2 (x2, y2)
(Shape functions)                                                                      (u2, v2)
1 (x1, y1)
(u1, v1)                            x, u
Shape functions construction

Assume,
N1  a1  b1 x  c1 y   Ni  ai  bi x  ci y       i= 1, 2, 3
N2  a2  b2 x  c2 y                          ai 
 
or    Ni  1 x     y bi   pT 
N3  a3  b3 x  c3 y                          c 
pT       i

Shape functions construction

Delta function property:
N1 ( x1 , y1 )  1
1   for i  j
Ni ( x j , y j )                                N1 ( x2 , y2 )  0
0    for i  j
N1 ( x3 , y3 )  0

Therefore, N1 ( x1 , y1 )  a1  b1 x1  c1 y1  1
N1 ( x2 , y2 )  a1  b1 x2  c1 y2  0
N1 ( x3 , y3 )  a1  b1 x3  c1 y3  0

x2 y3  x3 y2        y2  y3        x3  x2
Solving,        a1                , b1          , c1 
2 Ae              2 Ae           2 Ae
Shape functions construction
1 x1      y1
1    1                 1
Ae      P  1 x2         y2  [( x2 y3  x3 y2 )  ( y2  y3 ) x1  ( x3  x2 ) y1 ]
2    2                 2
1 x3        y3

Area of triangle           Moment matrix

Substitute a1, b1 and c1 back into N1 = a1 + b1x + c1y:
1
N1         [( y2  y3 )( x  x2 )  ( x3  x2 )( y  y2 )]
2 Ae
Shape functions construction
Similarly,
1
N 2 ( x1 , y1 )  0   N2          [( x3 y1  x1 y3 )  ( y3  y1 ) x  ( x1  x3 ) y ]
2 Ae
N 2 ( x2 , y2 )  1
1
N 2 ( x3 , y3 )  0             [( y3  y1 )( x  x3 )  ( x1  x3 )( y  y3 )]
2 Ae
1
N3 ( x1 , y1 )  0    N3       [( x1 y2  x1 y1 )  ( y1  y2 ) x  ( x2  x1 ) y ]
2 Ae
N3 ( x2 , y2 )  0
1
N3 ( x3 , y3 )  1                [( y1  y2 )( x  x1 )  ( x2  x1 )( y  y1 )]
2 Ae
Shape functions construction

Ni  ai  bi x  ci y
1
where    ai       ( x j yk  xk y j )    i= 1, 2, 3
2 Ae
1                          J, k determined from cyclic
bi       ( y j  yk )           permutation
2 Ae
i = 1, 2
1                                     i
ci       ( xk  x j )
2 Ae
k          j
k = 3, 1                  j = 2, 3
Using area coordinates
   Alternative method of constructing shape
functions

1 x       y
1              1
 2-3-P:        A1      1 x2     y2  [( x2 y3  x3 y2 )  ( y2  y3 ) x  ( x3  x2 ) y ]
y                                     2              2
1 x3     y3
k, 3
A1
L1 
Ae
A2
A1           Similarly,  3-1-P                 A2            L2 
P                                                                          Ae
A3
i, 1                                               1-2-P             A3            L3 
j, 2                                                   Ae
x
Using area coordinates
Partitions of unity:    L1  L2  L3  1

A1 A2 A3 A1  A2  A3
L1  L2  L3                    1
Ae Ae Ae      Ae

Delta function property: e.g. L1 = 0 at if P at nodes 2 or 3

Therefore,   N1  L1 , N 2  L2 , N3  L3

U h ( x, y )  N( x, y )de
Strain matrix
u                                           
 xx                                             0
 x                                 x        
v                                         
 yy                  LU     where L   0
 y
        y 

u v                              
 xy                                             
 y  x                                       
 y
        x 

  LU  LNd e  Bd e
         
       0
 x                    b1 0 b2 0 b3 0 
       
B  LN   0        N       B   0 c1 0 c2 0 c3 
                      
       y 

c1 b1 c2 b2 c3 b3 
                      

          
 y
       x 
             (constant strain element)
Element matrices
h
k e   B cBdV   (  dz )BT cBdA   hBT cBdA
T
0
Ve            Ae                    Ae

Constant matrix
 k e  hAe BT cB

h
me   N NdV 
T
       dz NT NdA   hNT NdA
0
Ve            Ae                   Ae
Element matrices
For elements with uniform density and thickness,
 N1 N1     0       N1 N 2    0       N1 N 3     0 
 0        N1 N1     0       N1 N 2    0       N1 N 3 
                                                     
 N 2 N1     0      N2 N2      0      N 2 N3     0 
me  h                                                           dA
Ae    0        N 2 N1     0      N2 N2      0      N 2 N3 
 N 3 N1     0      N3 N 2     0      N3 N3      0 
                                                     
 0
          N3 N1      0      N3 N 2     0      N3 N3 

Eisenberg and Malvern (1973):

m! n! p!
    L1 Ln L3 dA 
m     p
2A
(m  n  p  2)!
2
A
y, v                3 (x3, y3)
Element matrices                                          (u3, v3)
fsy
fsx
2 0 1 0 1 0                                   A
  2 0 1 0 1
                                                                 2 (x2, y2)
hA      2 0 1 0                                                       (u2, v2)
me                                      1 (x1, y1)
12         2 0 1                      (u1, v1)                            x, u
 sy.     2 0
            

           2


 f sx 
fe   NT
  dl                     0
l  2 3
 f sy                    0
 
1        fx 
 
Uniform distributed load:   f e  l 2 3  
2        fy
 fx 
 
 fy
 
LINEAR RECTANGULAR
ELEMENTS
 Non-constant strain matrix
 More accurate representation of stress and strain
 Regular shape makes formulation easy
Shape functions construction

Consider a rectangular element

 u1                                 y, v
   displacements at node 1
v                                                             
 1                                               4 (x4, y4)             3 (x3, y3)
u2     
          displacements at node 2                 (u4, v4)               (u3, v3)
v      
de   2                                                                      fsy
u3                                         2b
 v3       displacements at node 3                                               fsx
       
2a
u4     
v         displacements at node 4               1 (x1, y1)             2 (x2, y2)
 4     
(u1, v1)               (u2, v2)

x, u
Shape functions construction
4 (1, +1)                 3 (1, +1)
y, v                                                              (u4, v4)                   (u3, v3)

4 (x4, y4)             3 (x3, y3)
(u4, v4)               (u3, v3)                                    
     fsy
2b                                   fsx                  2

2a
1 (x1, y1)             2 (x2, y2)                                2
(u1, v1)               (u2, v2)                   1 (1, 1)                 2 (1, 1)
(u1, v1)                   (u2, v2)
x, u

x  ( x2  x1 ) / 2      y  ( y2  y1 ) / 2
                        , 
a                        b

 N1   0          N2         0 N3         0         N4   0
U ( x, y )  N( x, y )de where
h
N
0     N1         0          N2 0         N3        0    N4 

Node 1           Node 2            Node 3         Node 4
Shape functions construction
N 1  1 (1   )(1   )
4

N 2  1 (1   )(1   )
4

N 3  1 (1   )(1   )
N j  1 (1   j  )(1   j )
4
4

N 4  1 (1   )(1   )                                                        4 (1, +1)                3 (1, +1)
4                                                                           (u4, v4)                  (u3, v3)

N 3 at node 1  1 (1   )(1   )  1  0
4
 1                                                          
N 3 at node 2  1 (1   )(1   )  1  0
4                                   Delta function
 1                                          2
N 3 at node 3  1 (1   )(1   )  1  1         property                                         
4
 1

N 3 at node 4  1 (1   )(1   )  1  0
4
2
 1                                          1 (1, 1)                2 (1, 1)
4                                                                                (u1, v1)                  (u2, v2)

        Ni  N1  N 2  N 3  N 4
i 1                                                                                 Partition of
 [(1   )(1   )  (1   )(1   )  (1   )(1   )  (1   )(1   )]
1
4                                                                               unity
 1 [2(1   )  2(1   )]  1
4
Strain matrix

  1
a     0     1
a      0      1
a    0      1
 a      0 
1          1            1         1           1 
B  LN   0        b     0      b     0      b     0       b 
4 1
         1
 a
1
 b
1    1   1   1
 a 
1
 b                       a      b     a     b           

Note: No longer a constant matrix!
Element matrices

  x a,   y b        dxdy = ab dd

Therefore,

1 1
k e   hB cBdA  
T
   abhB T cBdd
1 1
A

h                           1       1
me   N NdV    dz N NdA   hN NdA  
T                       T             T
        abh NT Nd d
0                           1   1
V             A                    A
Element matrices

 f sx 
fe   N T
  dl
2 3
 f sy 
l

0                        y, v

0                                    4 (x4, y4)             3 (x3, y3)
                                       (u4, v4)               (u3, v3)
 fx                                                         fsy
                                2b                                   fsx
 fy
f e  b                                                  2a
 fx                                  1 (x1, y1)             2 (x2, y2)
 fy                                    (u1, v1)               (u2, v2)
 
0                                                                   x, u
0 
 
Gauss integration

   For evaluation of integrals in ke and me (in practice)

1                  m
In 1 direction: I           f ( )d   w j f ( j )
1
j 1

m gauss points gives exact solution of
polynomial integrand of n = 2m - 1

nx     ny
1       1
In 2 directions: I                  f ( , )d d            wi w j f (i , j )
1       1
i 1   j 1
Gauss integration
m        Gauss Point j            Gauss Weight wj   Accuracy order
n
1   0                          2                           1

2   -1/3, 1/3                1, 1                        3

3   -0.6, 0, 0.6             5/9, 8/9, 5/9               5

4   -0.861136, -0.339981,      0.347855, 0.652145,         7
0.339981, 0.861136         0.652145, 0.347855
5   -0.906180, -0.538469, 0,   0.236927, 0.478629,         9
0.538469, 0.906180         0.568889, 0.478629,
0.236927
6   -0.932470, -0.661209,      0.171324, 0.360762,        11
-0.238619, 0.238619,       0.467914, 0.467914,
0.661209, 0.932470         0.360762, 0.171324
Evaluation of me

4 0 2 0 1 0 2 0
  4 0 2 0 1 0 2
                
     4 0 2 0 1 0
                
hab        4 0 2 0 1
me 
9            4 0 2 0
                
           4 0 2
 sy.         4 0
                

               4

Evaluation of me
1 1
mij  hab               N i N j dd
1 1

hab       1                                1

16 1           (1   i )(1   j )d 
1
(1  i )(1   j )d

hab
      (1  1  i j )(1  1 i j )
3              3
4
hab                                                 hab
E.g. m33                    (1   1  1)(1   1  1)  4 
1
3
1
3
4                                                 9

Note: In practice, gauss integration is often used
ELEMENTS
 Rectangular elements have limited application
 Quadrilateral elements with unparallel edges are
more useful
 Irregular shape requires coordinate mapping
before using Gauss integration
Coordinate mapping
y     4 (x4, y4)
3 (x3, y3)

1,
4 ( +1)        3 (1, +1)


2 (x2, y2)
1 (x1, y1)
x         1 ( 
1, 1)       2 (1, 
1)

Physical coordinates              Natural coordinates

U h ( , )  N( , )d e (Interpolation of displacements)

X( , )  N( , )xe (Interpolation of coordinates)
Coordinate mapping
 x1           
y                coordinate at node 1

X( , )  N( , )xe            1
 x2           
                 coordinate at node 2
x                y2           
where       X  ,        xe   
 y               x3           
 y3              coordinate at node 3
              
 x4           
N1  1 (1   )(1   )         y                coordinate at node 4
4                           4            
N 2  1 (1   )(1   )
4                                   4
N 3  1 (1   )(1   )
4                           x   N i ( ,  ) xi
i 1
N 4  1 (1   )(1   )
4
4
y   N i ( , ) y i
i 1
Coordinate mapping
4
Substitute 1 into x   N i ( , ) xi
i 1

x  1 (1   ) x2  1 (1   ) x3
2               2                                    x  1 ( x 2  x3 )  1  ( x3  x 2 )
2                2
or
y  (1   ) y2  (1   ) y3
1
2
1
2                                   y  1 ( y 2  y3 )  1 ( y3  y 2 )
2                2

( y3  y2 )
Eliminating ,               y                 {x  1 ( x2  x3 )}  1 ( y2  y3 )
( x3  x2 )      2                2

y   4 (x4, y4)
3 (x3, y3)

4 (1, +1)     3 (1, +1)


2 (x2, y2)
1 (x1, y1)
x             1 (1, 1)   2 (1, 1)
Strain matrix
N i N i    x N i    y               N i     N i 
                                          x 
    x        y    
or           J      
N i N i    x N i     y              N i     N i 
                                          y 
    x      y                                   

 x     y 
       
where      J                 (Jacobian matrix)
 x     y 
       
           
 N1   N 2    N3   N 4   x1   y1 
                   x2    y2 
Since X( , )  N( , )xe ,     J                                    
 N1   N 2    N3   N 4   x3   y3 
                                   
                    x4
       y4 

Strain matrix
 N i       N i 
 x           (Relationship between differentials of shape
Therefore,         J 
1
 functions w.r.t. physical coordinates and
 N i       N i 
 y           differentials w.r.t. natural coordinates)
                 

Replace differentials of Ni w.r.t. x and y
with differentials of Ni w.r.t.  and 

 x  0 
B  LN   0
      y  N

 y  x 
          
Element matrices
Murnaghan (1951) : dA=det |J | dd

1       1
ke                   hBTcBdet J d d
1      1

h
m e   N NdV    dxN T NdA   hN T NdA
T
0
V                     A               A
1 1
           hN T Ndet J dd
1 1
Remarks

   Shape functions used for interpolating the coordinates are
the same as the shape functions used for the interpolation
of the displacement field. Therefore, the element is called
isoparametric element.
   Note that the shape functions for coordinate interpolation
and displacement interpolation do not have to be the same.
   Using the different shape functions for coordinate
interpolation and displacement interpolation, respectively,
will lead to the development of so-called subparametric or
superparametric elements.
HIGHER ORDER ELEMENTS
   Higher order triangular elements
(0,0,p)

nd = (p+1)(p+2)/2
(1,0,p1)                (0,1,p1)
Node i,               I J K  p
(2,0,p2)
Argyris, 1968 :
L1
L2
i (I,J,K)
Ni  lII ( L1 )lJJ ( L2 )lK ( L3 )
K

L3                  (0,p1,1)               ( L  L 0 )( L  L 1 )     ( L  L (  1) )
l ( L ) 
( L   L 0 )( L   L 1 )   ( L   L (  1) )
(p,0,0)   (p1,1,0)                                   (0,p,0)
HIGHER ORDER ELEMENTS
    Higher order triangular elements (Cont’d)
y, v           3           N1  (2 L1  1) L1
N 4  4 L1 L2
5
6

2                                                                  1
y, v               3                    N1      (3L1  1)(3L1  2) L1
4                                                                           2
1                                                                    7                  9
x, u                      8                         N4     L1 L2 (3L1  1)
6          2
9                                N10  27 L1 L2 L3
10
2
5

x, u

Cubic element
HIGHER ORDER ELEMENTS
   Higher order rectangular elements
Lagrange type:
                       (Zienkiewicz et al., 2000)
(0,m)                           (n,m)

Ni  N I1D N J D  lIn ( )lJm ( )
1

i(I,J)
(  0 )(  1 ) (   k 1 )(   k 1 ) (   n )
   lkn ( ) 
( k  0 )( k  1 ) ( k   k 1 )( k   k 1 ) ( k   n )
0

(0,0)                             (n,0)
HIGHER ORDER ELEMENTS
   Higher order rectangular elements(Cont’d)
4    7       3
J=2



6
J=1     8     9

1
J=0                                              N 5  N1 D ( ) N 0 D ( )   (1   )(1   )(1   )
1         1

1           5          2                                                 2
I=0         I=1        I=2
1
N 6  N 2 D ( ) N1 D ( )   (1   )(1   )(1   )
1         1

1                                                      2
N1  N 0 D ( ) N 0 D ( )   (1   ) (1   )
1          1
1
4                         N 7  N1 D ( ) N 2 D ( )  (1   )(1   )(1   )
1         1

1                                                    2
N 2  N 2 D ( ) N 0 D ( )    (1   ) (1   )
1          1
1
4                       N8  N 0 D ( ) N11D ( )    (1   )(1   )
1

1                                                        2
N 3  N 2 D ( ) N 2 D ( )   (1   )(1   )
1          1
N 9  N1 D ( ) N1 D ( )  (1   2 )(1   2 )
1         1
4
1
N 4  N 0 D ( ) N 2 D ( )    (1   )(1   )
1          1

4
HIGHER ORDER ELEMENTS
    Higher order rectangular elements(Cont’d)
Serendipity type:
=1        
4              7       3

                                     N j  1 (1   j )(1   j )( j   j  1)
4                                           j  1, 2, 3, 4
6          N j  1 (1   2 )(1   j )
2                                           j  5, 7
8                  0               
N j  1 (1   j )(1   2 )
2                                           j  6,8

1                  5   2
=1

HIGHER ORDER ELEMENTS
    Higher order rectangular elements(Cont’d)

                    N j  32 (1   j )(1   j )(9 2  9 2  10)
1

4        10       9   3
for corner nodes j  1, 2, 3, 4

8
N j  32 (1   j )(1   2 )(1  9 j )
9
11

for side nodes j  7, 8, 11, 12 where  j  1 and  j   1
3
12                    7
N j  32 (1   j )(1   2 )(1  9 j )
9

6   2
for side nodes j  5, 6, 9, 10 where  j   1 and  j  1
3
1        5

(twelve node cubic element)
ELEMENT WITH CURVED
EDGES
3                               3
5
5
6
6
2
4               2                   4
1
1

4       7           3
3
4   7

8                   6                               6
8

1                   2
5
1           5               2

   When the Gauss integration scheme is used, one has to
decide how many Gauss points should be used.
   Theoretically, for a one-dimensional integral, using m
points can give the exact solution for the integral of a
polynomial integrand of up to an order of (2m1).
   As a general rule of thumb, more points should be used for
higher order of elements.
INTEGRATION)
   Using smaller number of Gauss points tends to counteract
the over-stiff behaviour associated with the displacement-
based method.
   Displacement in an element is assumed using shape
functions. This implies that the deformation of the element
is somehow prescribed in a fashion of the shape function.
This prescription gives a constraint to the element. The so-
constrained element behaves stiffer than it should be. It is
often observed that higher order elements are usually softer
than lower order ones. This is because using higher order
elements gives less constraint to the elements.
INTEGRATION)
   Two gauss points for linear elements, and two or three
points for quadratic elements in each direction should be
sufficient for many cases.
   Most of the explicit FEM codes based on explicit
formulation tend to use one-point integration to achieve the
best performance in saving CPU time.
CASE STUDY
   Side drive micro-motor
CASE STUDY
10N/m
Elastic Properties of Polysilicon
10N/m
Young’s Modulus, E             169GPa

Poisson’s ratio,              0.262

10N/m       Density,                2300kgm-3
CASE STUDY

Analysis no. 1: Von Mises stress
distribution using 24 bilinear
CASE STUDY

Analysis no. 2: Von Mises stress
distribution using 96 bilinear
CASE STUDY

Analysis no. 3: Von Mises stress
distribution using 144 bilinear
CASE STUDY

Analysis no. 4: Von Mises stress
distribution using 24 eight-nodal,
CASE STUDY

Analysis no. 5: Von Mises stress
distribution using 192 three-nodal,
triangular elements (129 nodes)
CASE STUDY
Total number    Maximum
Analysis   Number / type of
of nodes in    Von Mises
no.         elements
model      Stress (GPa)

24 bilinear,
1                              41           0.0139

96 bilinear,
2                              129          0.0180

144 bilinear,
3                              185          0.0197