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

Department of Engineering Mechanics PhD. TRUONG Tich Thien Introduction to FEM 24 Implementation of Iso-P Triangular Elements IFEM Ch 24 – Slide 1

Department of Engineering Mechanics PhD. TRUONG Tich Thien Introduction to FEM Peculiarities of Implementation of Iso-P Triangular Elements ● Special Gauss quadrature rules ● Computation of Jacobian and shape function x-y derivatives IFEM Ch 24 – Slide 2

Department of Engineering Mechanics PhD. TRUONG Tich Thien 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 3

Department of Engineering Mechanics PhD. TRUONG Tich Thien 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.22338 0.13239 0.22338 0.13239 0.12594 0.10995 0.12594 0.10995 IFEM Ch 24 – Slide 4

Department of Engineering Mechanics PhD. TRUONG Tich Thien Introduction to FEM Gauss Rules for Straight Sided Triangles 1-point (centroid) rule 1 � � ( e ) F (ζ 1 , ζ 2 , ζ 3 ) d � ( e ) ≈ F ( 1 3 , 1 3 , 1 3 ), A 3-internal-point rule � 1 � ( e ) F (ζ 1 , ζ 2 , ζ 3 ) d � ( e ) ≈ 1 3 F ( 2 3 , 1 6 , 1 6 ) + 1 3 F ( 1 6 , 2 3 , 1 6 ) + 1 3 F ( 1 6 , 1 6 , 2 3 ). A 3 midpoint rule � 1 � ( e ) F (ζ 1 , ζ 2 , ζ 3 ) d � ( e ) ≈ 1 3 F ( 1 2 , 1 2 , 0 ) + 1 3 F ( 0 , 1 2 , 1 2 ) + 1 3 F ( 1 2 , 0 , 1 2 ). A In the above A is the triangle area   1 1 1 � � ( e ) d � ( e ) = 1 A 2 det x 1 x 2 x 3 =   y 1 y 2 y 3 For 6 and 7 point rules see notes IFEM Ch 24 – Slide 5

Department of Engineering Mechanics PhD. TRUONG Tich Thien Introduction to FEM Triangles of Arbitrary Geometry (e.g. Curved Sides) Introduce the Jacobian determinant 1 1 1   n n n ∂ N i ∂ N i ∂ N i � � �   x i x i x i   d � ( e ) = J d ζ 1 d ζ 2 d ζ 3 , ∂ζ 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 Centroid rule � � ( e ) F (ζ 1 , ζ 2 , ζ 3 ) d � ( e ) ≈ J ( 1 3 , 1 3 , 1 3 ) F ( 1 3 , 1 3 , 1 3 ). Midpoint rule � � ( e ) F (ζ 1 , ζ 2 , ζ 3 ) d � ( e ) ≈ 1 3 J ( 1 2 , 1 2 , 0 ) F ( 1 2 , 1 2 , 0 ) + 1 3 J ( 0 , 1 2 , 1 2 ) F ( 0 , 1 2 , 1 2 ) + 1 3 J ( 1 2 , 0 , 1 2 ) F ( 1 2 , 0 , 1 2 ). etc IFEM Ch 24 – Slide 6

Department of Engineering Mechanics PhD. TRUONG Tich Thien 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]]]; ]; rule = defines Gauss rule (1, 3, -3, 6 and 7 implemented) numer = True or False to get numeric or exact info point = index of Gauss point, ranges from 1 through Abs[rule] returns abscissa { ζ , ζ , ζ } and weight 1 3 2 IFEM Ch 24 – Slide 7

Department of Engineering Mechanics PhD. TRUONG Tich Thien Introduction to FEM Computation of Shape Function Derivatives 3 5 Illustrated for 6-node 6 quadratic triangle: 2 4 1  N 1  1 1 1 1 1 1 1     N 2   x x 1 x 2 x 3 x 4 x 5 x 6   Recall the iso-P N 3       y y 1 y 2 y 3 y 4 y 5 y 6   =      element definition: 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 8

Department of Engineering Mechanics PhD. TRUONG Tich Thien 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 can be any interpolated quantity. Then � ∂ N i � ∂w ∂ N i ∂ζ 1 ∂ x + ∂ N i ∂ζ 2 ∂ x + ∂ N i ∂ζ 3 � � w i w i ∂ x = ∂ x = ∂ζ 1 ∂ζ 2 ∂ζ 3 ∂ x � ∂ N i � ∂w ∂ 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 9

Department of Engineering Mechanics PhD. TRUONG Tich Thien 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 and stack the results row-wise w ≡ 1 , x , y � ∂ 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 10

Department of Engineering Mechanics PhD. TRUONG Tich Thien 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 ∂ζ 2 ∂ζ 3 1 0  =       ∂ 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 11

Department of Engineering Mechanics PhD. TRUONG Tich Thien 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 JP = R where J and R are known. IFEM Ch 24 – Slide 12

Department of Engineering Mechanics PhD. TRUONG Tich Thien 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 13