Finite element method (continued)

Document Sample
Finite element method (continued) Powered By Docstoc
					                    EN0175                                                        11 / 30/ 06


Finite element method (continued)

Summary of essential concepts:

FEM analysis begins with calculations on the element level.

Element

                                                   #3




                        #1                                                    #2
Element nodal displacement:

                                                      ⎧ u1#1 ⎫
                                                      ⎪ #1 ⎪
                                                      ⎪u2 ⎪
                                                      ⎪u1# 2 ⎪
                                                      ⎪ ⎪
                                               u el = ⎨ # 2 ⎬
                                                      ⎪u 2 ⎪
                                                      ⎪u1#3 ⎪
                                                      ⎪ #3 ⎪
                                                      ⎪u2 ⎪
                                                      ⎩ ⎭


Element interpolation function:

                              N 1 ( x1 , x2 ) , N 2 (x1 , x2 ) , N 3 (x1 , x2 )


Element displacement field:

                        ⎧ u1 ⎫ ⎡ N 1        0       N2      0      N3      0⎤
                        ⎨ ⎬=⎢                  1               2             ⎥ u el
                        ⎩u 2 ⎭ ⎣ 0         N        0      N        0      N3⎦


Element strain field:




                                                     1
                     EN0175                                                              11 / 30/ 06



                              ⎡ ∂N 1                       ∂N 2              ∂N 3            ⎤
                              ⎢                 0                        0               0 ⎥
                                ∂x                         ∂x1               ∂x1
                   ⎡ ε 11 ⎤ ⎢ 1                                                              ⎥
                                             ∂N 1                    ∂N 2               ∂N 3 ⎥
               ε = ⎢ ε 22 ⎥ = ⎢ 0
                          ⎥ ⎢                                0                0                u el = Bu el
                   ⎢                         ∂x2                     ∂x2                ∂x2 ⎥
                   ⎢2ε 12 ⎥ ⎢ 1
                   ⎣      ⎦                                                                  ⎥
                              ⎢ ∂N           ∂N 1          ∂N 2      ∂N 2    ∂N 3       ∂N 3 ⎥
                              ⎢ ∂x2
                              ⎣              ∂x1           ∂x2       ∂x1     ∂x2        ∂x1 ⎥⎦


Element stress field:

                                                     ⎡          ⎤
                        ⎡σ 11 ⎤                      ⎢1 ν    0 ⎥ ⎡ ε 11 ⎤
                    σ = ⎢σ 22 ⎥ =                            0 ⎥ ⎢ ε 22 ⎥ = Dε = D Bu el
                                  E
                                                     ⎢ν 1
                        ⎢ ⎥ 1 −ν 2                   ⎢     1 −ν ⎥ ⎢      ⎥
                        ⎢σ 12 ⎥
                        ⎣ ⎦                          ⎢ 0 0      ⎥⎣⎢2ε 12 ⎥
                                                                         ⎦
                                                     ⎣       2 ⎦


Element strain energy:
                                                1                    1 T
                                     U el =
                                                2 ∫Vel σ ijε ij dV = 2 u el K el u el
                                                     K el = Ael B D B
                                                                     T




Element nodal force:

                                      ∫
                                      Vel
                                             f i ui dV + ∫ ti ui dS = F el u el
                                                              Ael



After the element quantities are calculated, the next step is to assemble the global stiffness matrix.

Global strain energy:
                                                                  1 T               1 T
                                U=    ∑U
                                     elements
                                                el   =    ∑
                                                         elements 2
                                                                    u el K el u el = u K u
                                                                                    2

                                                ⎧ u1(1) ⎫
                     ⎧u 1
                         #1
                              ⎫                 ⎪ (1) ⎪
                     ⎪   #1   ⎪                 ⎪u2 ⎪
                     ⎪u 2     ⎪                 ⎪u1(2 ) ⎪
                     ⎪u
                     ⎪
                        #2    ⎪
                              ⎪                 ⎪ ( ⎪
              u el = ⎨  1
                        #2    ⎬             u = ⎨u22 ) ⎬                            (K el )6×6         K 2 n×2 n
                     ⎪u 2     ⎪                 ⎪u (3 ) ⎪
                     ⎪u #3    ⎪                 ⎪ 1 ⎪
                     ⎪  1
                              ⎪                     (
                                                ⎪u 23 ) ⎪
                     ⎪u
                     ⎩
                        #3
                        2     ⎪6×1
                              ⎭                 ⎪ ⎪
                                                ⎩ M ⎭ 2 n×1


Element connectivity

                                             (#1, #2, #3) ⇒ (a, b, c )
                   α (1, 2, 3, 4, 5, 6) ⇒ zα (2a − 1, 2a, 2b − 1, 2b, 2c − 1, 2c )

                                                                 2
                     EN0175                                           11 / 30/ 06



           K zα z β = K zα zβ + Kαβ (assembly of global stiffness matrix over all elements)
                                 el




Nodal force

                                          Fzα = Fzα + Fαel


FEM equation:

                                       K 2 n×2 n u 2 n×1 = F 2 n×1


For n nodes, there will be 2n unknown displacement components (2 degrees of freedom each
node) to be solved. For displacement boundary value problems, there are further restrictions on the
degrees of freedom.



                5                               4                               3

                           1                                         3
                                                 2
               1                                                                    2
At node 1, the displacements at both directions are fixed; at node 2, the displacement in 2-
direction is fixed. The global displacement vector becomes
                                                     ⎡0 ⎤
                                                     ⎢0 ⎥
                                                     ⎢ ⎥
                                                     ⎢×⎥
                                            u 2n×1 = ⎢ ⎥
                                                     ⎢0 ⎥
                                                     ⎢×⎥
                                                     ⎢ ⎥
                                                     ⎢M⎥
                                                     ⎣ ⎦

                     (1)
Example: suppose u2 = Δ

For this prescribed displacement, we can simply replace the equation for the appropriate degree of
freedom with the displacement constraint, but the resulting stiffness matrix is no longer
symmetric. This can be remedied by placing the prescribed displacement to the right side of
equation. The resulting FEM equation becomes



                                                   3
                         EN0175                                                    11 / 30/ 06



             ⎡ K11      0     K13         L L L K1, 2 n ⎤ ⎡ u1(1) ⎤ ⎡ F1(1) − K12 Δ ⎤
                                                                        ⎢                      ⎥
             ⎢ 0        1      0          0 L L    0 ⎥ ⎢ u 21) ⎥ ⎢
                                                                 (
                                                                                  Δ            ⎥
             ⎢                                              ⎥⎢ ⎥
             ⎢ K 31     0     K 33        L L L K 3, 2 n ⎥ ⎢u1(2 ) ⎥ ⎢ F1(2 ) − K 32 Δ ⎥
             ⎢                                              ⎥⎢ ( ⎥ ⎢                           ⎥
             ⎢ M        0       M         M M M    M ⎥ ⎢u 22 ) ⎥ = ⎢ F2( 2 ) − K 42 Δ ⎥
             ⎢ M        M       M         M M M    M ⎥⎢ M ⎥ ⎢                     M            ⎥
             ⎢                                              ⎥ ⎢ (n ) ⎥ ⎢ (n )                  ⎥
             ⎢ M        M   M             M M M    M ⎥ ⎢u1 ⎥ ⎢ F1 − K 2 n −1, 2 Δ ⎥
             ⎢K
             ⎣ 2 n ,1   0 K 2 n ,3        L L L K 2 n , 2 n ⎥ ⎢u 2n ) ⎥ ⎢ F2(n ) − K 2 n , 2 Δ ⎥
                                                                 (
                                                            ⎦⎣ ⎦ ⎣                             ⎦ 3n×1

The modified K then remains symmetric, as shown above.


Post-processing

                    Extract nodal displacement using element connectivity: u ⇒ u el

                                                       ε = Bu el

                                                      σ = D Bu el


Isoparametric elements:


The actual nodal coordinates can be mapped to a normalized domain ξ ⊂ [− 1, 1] , which is

particularly convenient because displacements and positions for different elements can be
interpolated using the same shape functions.


                                     ξ                                                                  x
 −1           0          1                                            #1 #2

u (x ) = u #1 N 1 (ξ ) + u # 2 N 2 (ξ )

x = x #1 N 1 (ξ ) + x # 2 N 2 (ξ )


           ⎧1, ξ = −1                           ⎧0, ξ = −1
N 1 (ξ ) = ⎨          ,              N 2 (ξ ) = ⎨
           ⎩ 0, ξ = 1                           ⎩ 1, ξ = 1

Interpolation functions:

                                                  1− ξ                          1+ ξ
                                     N 1 (ξ ) =        ,           N 2 (ξ ) =
                                                   2                             2


                                                           4
                           EN0175                                                        11 / 30/ 06




2D linear triangular element:

                   ξ2
                                                                                     c
                      #2
                  1



                      #3                     #1
                                                       ξ1               a                                 b
                  0                         1
v v                       v                     v
u = u a N 1 (ξ1 , ξ 2 ) + u b N 2 (ξ1 , ξ 2 ) + u c N 3 (ξ1 , ξ 2 )
v v                       v                     v
x = x a N 1 (ξ1 , ξ 2 ) + x b N 2 (ξ1 , ξ 2 ) + x c N 3 (ξ1 , ξ 2 )


Interpolation functions:

                       N 1 (ξ ) = ξ1 ,               N 2 (ξ ) = ξ 2 ,           N 3 (ξ ) = 1 − ξ1 − ξ 2


2D quadratic triangular element:

                      ξ2
                                                                                      3
                      1
                        #2
                                                                                6               5
                  #5                 #4

                          #3                    #1
                                                         ξ1                 1             4               2
                      0          #6             1

Interpolation functions:

                                                       N 1 = (2ξ1 − 1)ξ1

                                                      N 2 = (2ξ 2 − 1)ξ 2

                                         N 3 = (2(1 − ξ1 − ξ 2 ) − 1)(1 − ξ1 − ξ 2 )

                                                          N 4 = 4ξ1ξ 2

                                                    N 5 = 4ξ 2 (1 − ξ1 − ξ 2 )

                                                                5
                     EN0175                                                 11 / 30/ 06



                                        N 6 = 4ξ1 (1 − ξ1 − ξ 2 )


3D linear tetrahedral element:

                                           ξ3

                                         1 #3



                                              #4                    #2
                                                                            ξ2
                                              0                     1
                            #1
                                    1
                           ξ1
Interpolation functions:

                                                  N 1 = ξ1

                                                N 2 = ξ2

                                                N 3 = ξ3

                                        N 4 = 1 − ξ1 − ξ 2 − ξ 3


3D brick element:

                                              ξ3
                                        5
                                                                        8
                            6                        7
                                                                               ξ2

                                          1                         4
                      ξ1
                                2                            3
                                         − 1 ≤ ξ1 , ξ 2 , ξ3 ≤ 1

                                                     6
                          EN0175                                                             11 / 30/ 06




Interpolation functions:

                     N1 =
                              1
                                (1 − ξ1 )(1 − ξ 2 )(1 − ξ3 ) ,        N2 =
                                                                                1
                                                                                  (1 + ξ1 )(1 − ξ 2 )(1 − ξ3 )
                              8                                                 8


                     N3 =
                              1
                                (1 + ξ1 )(1 + ξ 2 )(1 − ξ3 ) ,        N4 =
                                                                                1
                                                                                  (1 − ξ1 )(1 + ξ 2 )(1 − ξ3 )
                              8                                                 8


                    N5 =
                              1
                                (1 − ξ1 )(1 − ξ 2 )(1 + ξ3 ) ,        N6 =
                                                                                1
                                                                                  (1 + ξ1 )(1 − ξ 2 )(1 + ξ3 )
                              8                                                 8


                    N7 =
                              1
                                (1 + ξ1 )(1 + ξ 2 )(1 + ξ3 ) ,         N8 =
                                                                                1
                                                                                  (1 − ξ1 )(1 + ξ 2 )(1 + ξ3 )
                              8                                                 8


Element stiffness matrix:


                     ∂N a ∂N b       1 1 1      ∂N a (ξ1 , ξ 2 , ξ 3 ) ∂N b (ξ1 , ξ 2 , ξ 3 )
K aibk = ∫ Cijkl
  el
                               dV = ∫ ∫ ∫ Cijkl                                               J dξ1dξ 2 dξ 3
          Vel        ∂x j ∂xl        −1 −1 −1         ∂x j                    ∂xl

       ⎛ ∂x         ⎞
J = det⎜ i          ⎟ is the Jacobian associated with the mapping, which can be computed by
       ⎜ ∂ξ         ⎟
       ⎝ j          ⎠

                                          ∂xi   ∂              ∂N a a
                                              =
                                          ∂ξ j ∂ξ j
                                                     a
                                                           (
                                                    N (ξ )xi =
                                                           a

                                                               ∂ξ j
                                                                    xi      )
                                                  −1
∂N a ∂N a       ⎛ ∂ξ l   ⎞ ∂N a ⎛ ∂x j        ⎞
     =          ⎜        ⎟=     ⎜             ⎟
∂x j   ∂ξ l     ⎜ ∂x     ⎟ ∂ξ ⎜ ∂ξ            ⎟
                ⎝ j      ⎠   l ⎝     l        ⎠

             ∂x j
Calculate:           = A jl
             ∂ξ l

        ∂ξ l
Then:        = A−1
        ∂x j
                jl




Integration scheme for isoparametric elements:

Gaussian integration:
                                                                     NI
                                                       f (ξ )dξ = ∑ WI f (ξ I )
                                                  1
                                               ∫  −1
                                                                     I =1


ξ I : Gaussian integration points


                                                                 7
                          EN0175                                         11 / 30/ 06



WI : weighting coefficients


Gaussian integration ensures that polynomials of order of (2 N I − 1) are exactly integrated.



For N I = 1 ,    ξ1 = 0 , W1 = 2 .


                                                  f (ξ )

                                                                     ξ
                               −1                 0        1
                                            ∫ f (ξ )dξ = 2 f (0)
                                             1

                                             −1

Check:

f (ξ ) = C0 + C1ξ

 ∫ (C       + C1ξ )dξ = 2C0 = 2 f (0 )
  1
        0
 −1




                           1         1
For N I = 2 ,    ξ1 = −       , ξ2 =    , W1 = W2 = 1 .
                            3         3

                                                    ⎛     1 ⎞    ⎛ 1 ⎞
                                     ∫ f (ξ )dξ = f ⎜ −
                                      1
                                                            ⎟+ f ⎜   ⎟
                                     −1
                                                    ⎝      3⎠    ⎝ 3⎠
Check:

f (ξ ) = C0 + C1ξ + C2ξ 2 + C3ξ 3

                                                 ⎛ 1 ⎞     ⎛ 1 ⎞
     f (ξ )dξ = ∫ (C0 + C2ξ 2 )dξ = 2C0 + C2 = f ⎜ −
  1              1                       2
 ∫−1             −1                      3       ⎝
                                                      ⎟+ f ⎜
                                                     3⎠    ⎝ 3⎠
                                                               ⎟



For N I = 3 ,    ξ1 = −0.775 , ξ 2 = 0 , ξ3 = 0.775 , W1 = W3 = 0.556 , W2 = 0.889 .

                    ∫ f (ξ )dξ = 0.556 f (− 0.775) + 0.889 f (0) + 0.556 f (0.775)
                     1

                     −1




                                                      8
                     EN0175                                          11 / 30/ 06




h-method: refine mesh, i.e. reduce the element size to achieve sufficient accuracy
p-method: increase order of polynomial interpolation function



Summary about integration schemes for 2D and 3D elements

Linear triangular element: 1 integration point
Quadratic triangular element: 3 integration points

Linear quadrilateral element: 4 integration points ( N I = 2 )

Quadratic quadrilateral element: 4 integration points ( N I = 2 )

Linear brick element: 8 integration points ( N I = 2 )

Quadratic brick element: 27 integration points ( N I = 3 )




                                                  9