24 Implementation of Iso-P Triangular Elements IFEM Ch 24 Slide - PDF document

Introduction to FEM 24 Implementation of Iso-P Triangular Elements IFEM Ch 24 – Slide 1

Introduction to FEM Peculiarities of Implementation of Iso-P Triangular Elements ● Special Gauss quadrature rules needed ● Computation of Jacobian and shape function x-y derivatives complicated by having 3 natural coordinates IFEM Ch 24 – Slide 2

Introduction to FEM Gauss Rules for Straight Sided Triangles Centroid rule: 1 point, degree 1 3 Area A 1 Ω e F ( ζ , ζ , ζ ) d Ω ≈ F ( , , ) 1 1 1 Ω e A 1 2 3 3 3 3 1 2 Three-interior-point rule: 3 points, degree 2 1 Ω e F ( ζ , ζ , ζ ) d Ω ≈ F ( , , ) + F ( , , ) + F ( , , ) 1 1 1 2 1 1 1 2 1 1 1 2 A 3 3 3 6 3 1 2 3 6 6 6 3 6 3 6 Midpoint rule: 3 points, degree 2 1 Ω e F ( ζ , ζ , ζ ) d Ω ≈ F ( , , 0) + F ( 0, , ) + F ( , 0, ) 1 1 1 1 1 1 1 1 1 A 1 2 3 3 2 2 3 2 3 2 2 2 For 6 and 7 point rules see Notes - pictures on next slide IFEM Ch 24 – Slide 3

Introduction to FEM The 5 Simplest Gauss Rules Drawn over Straight Sided Triangles (number annotated near sample point is the weight) (a) rule=1 (b) rule=3 (c) rule= − 3 degree=1 degree=2 degree=2 0.33333 0.33333 0.33333 1.0000 0.33333 0.33333 0.33333 (d) rule=6 (e) rule=7 0.10995 0.12594 degree=4 degree=5 0.13239 0.13239 0.22338 0.22338 0.22500 0.22338 0.12594 0.12594 0.10995 0.10995 0.13239 IFEM Ch 24 – Slide 4

Introduction to FEM Triangles of Arbitrary Geometry (e.g. Curved Sides) The element of area, d Ω , of an iso-P triangular element with n nodes can be expressed as d Ω d Ω = J d ζ d ζ d ζ 3 1 2 where J is the "Jacobian determinant" 1 1 1   n n n ∂ N i ∂ N i ∂ N i � � �   x i x i x i ∂ ζ 1 ∂ ζ 2 ∂ ζ 3   J = 1 2 det   i = 1 i = 1 i = 1     n n n ∂ N i ∂ N i ∂ N i   � � � y i y i y i ∂ ζ 1 ∂ ζ 2 ∂ ζ 3 i = 1 i = 1 i = 1 IFEM Ch 24 – Slide 5

Introduction to FEM Triangles of Arbitrary Geometry (cont'd) Centroid rule Ω e F ( ζ , ζ , ζ ) d Ω ≈ J ( , , ) F ( , , ) 1 1 1 1 1 1 2 1 3 3 3 3 3 3 3 Midpoint rule Ω e F ( ζ , ζ , ζ ) d Ω ≈ J ( , , 0) F ( , , 0 ) 1 1 1 1 1 3 1 2 3 2 2 2 2 1 1 1 1 1 1 1 1 1 + J ( 0, , ) F ( 0 , , ) + J ( , 0, ) F ( , 0 , ) 1 3 2 2 2 3 2 2 2 2 2 etc. This can be compactly expressed by saying that the integration rule is applied to J F - pictures on next slide IFEM Ch 24 – Slide 6

Introduction to FEM The 5 Simplest Gauss Rules Drawn over Arbitrary (Curved Side) Triangles (number annotated near sample point is the weight) (a) rule=1 (b) rule=3 (c) rule= − 3 0.33333 0.33333 1.0000 0.33333 0.33333 0.33333 0.33333 0.12594 0.10995 (d) rule=6 (e) rule=7 0.13239 0.22338 0.22500 0.13239 0.22338 0.22338 0.12594 0.13239 0.10995 0.12594 0.10995 IFEM Ch 24 – Slide 7

Introduction to FEM Triangle Gauss Quadrature Module TrigGaussRuleInfo[{rule_,numer_},point_]:= Module[ {zeta,p=rule,i=point,g1,g2,info=Null}, If [p== 1, info={{1/3,1/3,1/3},1}]; If [p==-3, zeta={1/2,1/2,1/2}; zeta[[i]]=0; info={zeta,1/3}]; If [p== 3, zeta={1/6,1/6,1/6}; zeta[[i]]=2/3; info={zeta,1/3}]; If [p== 6, If [i<=3, g1=(8-Sqrt[10]+Sqrt[38-44*Sqrt[2/5]])/18; zeta={g1,g1,g1}; zeta[[i]]=1-2*g1; info={zeta,(620+Sqrt[213125-53320*Sqrt[10]])/3720}]; If [i>3, g2=(8-Sqrt[10]-Sqrt[38-44*Sqrt[2/5]])/18; zeta={g2,g2,g2}; zeta[[i-3]]=1-2*g2; info={zeta,(620-Sqrt[213125-53320*Sqrt[10]])/3720}]]; If [p== 7, If [i==1,info={{1/3,1/3,1/3},9/40} ]; If [i>1&&i<=4,zeta=Table[(6-Sqrt[15])/21,{3}]; zeta[[i-1]]=(9+2*Sqrt[15])/21; info={zeta,(155-Sqrt[15])/1200}]; If [i>4, zeta=Table[(6+Sqrt[15])/21,{3}]; zeta[[i-4]]=(9-2*Sqrt[15])/21; info={zeta,(155+Sqrt[15])/1200}]]; If [numer, Return[N[info]], Return[Simplify[info]]]; ]; IFEM Ch 24 – Slide 8

Introduction to FEM Triangle Gauss Quadrature (cont'd) Invoked by saying {{ ζ 1, ζ 2, ζ 3},w}=TrigGaussRuleInfo[{rule,numer},point] rule = defines Gauss rule (1, 3, -3, 6 and 7 implemented) numer = True or False to get numeric or exact info, resp. point = index of Gauss point, ranges from 1 through Abs[rule] returns abscissa of sample point and weight IFEM Ch 24 – Slide 9

Introduction to FEM Computation of Shape Function Derivatives 3 5 Illustrated for 6-node 6 quadratic triangle: 2 1 4   N 1 1 1 1 1 1 1 1     N 2   Recall the iso-P x x 1 x 2 x 3 x 4 x 5 x 6   N 3       element definition: y y 1 y 2 y 3 y 4 y 5 y 6   =          N 4   u x u x1 u x2 u x3 u x4 u x5 u x6      N 5   u y u y1 u y2 u y3 u y4 u y5 u y6 N 6 N 1 = ζ 1 ( 2 ζ 1 − 1 ), N 4 = 4 ζ 1 ζ 2 , N 2 = ζ 2 ( 2 ζ 2 − 1 ), N 5 = 4 ζ 2 ζ 3 , N 3 = ζ 3 ( 2 ζ 3 − 1 ), N 6 = 4 ζ 3 ζ 1 . IFEM Ch 24 – Slide 10

Introduction to FEM Partial Derivative Computation (Cont'd)  ζ 1 (2 ζ 1 − 1)  ζ 2 (2 ζ 2 − 1)     ζ 3 (2 ζ 3 − 1)   w = [ w 1 w 2 w 3 w 4 w 5 w 6 ]   4 ζ 1 ζ 2     4 ζ 2 ζ 3   4 ζ 3 ζ 1 where w = w ( ζ ,ζ ,ζ ) can be any quantity interpolated 1 2 3 quadratically over the triangle. Then ∂w ∂ N i ∂ N i ∂ζ 1 ∂ x + ∂ N i ∂ζ 2 ∂ x + ∂ N i ∂ζ 3 � � w i w i ∂ x = ∂ x = ∂ζ 1 ∂ζ 2 ∂ζ 3 ∂ x ∂w ∂ N i ∂ N i ∂ζ 1 ∂ y + ∂ N i ∂ζ 2 ∂ y + ∂ N i ∂ζ 3 � � w i w i ∂ y = ∂ y = ∂ζ 1 ∂ζ 2 ∂ζ 3 ∂ y IFEM Ch 24 – Slide 11

Introduction to FEM Partial Derivative Computation (Cont'd) ∂ζ 1 ∂ζ 1   ∂ x ∂ y   � ∂w � � w i ∂ N i � w i ∂ N i � w i ∂ N i ∂ζ 2 ∂ζ 2 ∂w   � � =   ∂ x ∂ y ∂ζ 1 ∂ζ 2 ∂ζ 3 ∂ x ∂ y     ∂ζ 3 ∂ζ 3 ∂ x ∂ y Make w = 1, x, y and stack results rowwise : � ∂ N i � ∂ N i � ∂ N i   ∂ζ 1 ∂ζ 1 ∂ 1 ∂ 1     ∂ζ 1 ∂ζ 2 ∂ζ 3 ∂ x ∂ y ∂ x ∂ y   � x i ∂ N i � x i ∂ N i � x i ∂ N i ∂ζ 2 ∂ζ 2 ∂ x ∂ x        =      ∂ x ∂ y  ∂ζ 1 ∂ζ 2 ∂ζ 3 ∂ x ∂ y      ∂ y ∂ y ∂ζ 3 ∂ζ 3 � y i ∂ N i � y i ∂ N i � y i ∂ N i   ∂ x ∂ y ∂ x ∂ y ∂ζ 1 ∂ζ 2 ∂ζ 3 IFEM Ch 24 – Slide 12

Introduction to FEM Partial Derivative Computation (Cont'd) As shown in Notes, this system reduces to ∂ζ 1 ∂ζ 1 1 1 1     ∂ x ∂ y   0 0 � x i ∂ N i � x i ∂ N i � x i ∂ N i ∂ζ 2 ∂ζ 2     1 0 ∂ζ 1 ∂ζ 2 ∂ζ 3  =       ∂ x ∂ y    � y i ∂ N i � y i ∂ N i � y i ∂ N i 0 1 ∂ζ 3 ∂ζ 3 ∂ζ 1 ∂ζ 2 ∂ζ 3 ∂ x ∂ y unknowns are grouped here IFEM Ch 24 – Slide 13

Introduction to FEM Partial Derivative Computation (Cont'd) Replacing the shape functions of the 6-node triangle  1 1 x 1 ( 4 ζ 1 − 1 ) + 4 x 4 ζ 2 + 4 x 6 ζ 3 x 2 ( 4 ζ 2 − 1 ) + 4 x 5 ζ 3 + 4 x 4 ζ 1  y 1 ( 4 ζ 1 − 1 ) + 4 y 4 ζ 2 + 4 y 6 ζ 3 y 2 ( 4 ζ 2 − 1 ) + 4 y 5 ζ 3 + 4 y 4 ζ 1 ∂ζ 1 ∂ζ 1   ∂ x ∂ y    1 0 0 ∂ζ 2 ∂ζ 2   x 3 ( 4 ζ 3 − 1 ) + 4 x 6 ζ 1 + 4 x 5 ζ 2 1 0  =      ∂ x ∂ y  y 3 ( 4 ζ 3 − 1 ) + 4 y 6 ζ 1 + 4 y 5 ζ 2 0 1 ∂ζ 3 ∂ζ 3 ∂ x ∂ y or J P = R where J and R are known . IFEM Ch 24 – Slide 14

Introduction to FEM Partial Derivative Computation (Cont'd) Solve J P = R for the partials ∂ζ 1 ∂ζ 1   ∂ x ∂ y ∂ζ 2 ∂ζ 2   P =   ∂ x ∂ y   ∂ζ 3 ∂ζ 3 ∂ x ∂ y Plug these in the chain rule to get the x-y partial derivatives of the shape functions, and use these to form the strain-displacement matrix B IFEM Ch 24 – Slide 15