Integración semi-analítica de matrices de rigidez de elementos finitos en simetría axial

Rev. Int. M´et. Num. C´alc. Dis. Ing. Vol. 24, 4, 345-356 (2008) Revista Internacional de M´etodos Num´ericos para C´alculo y Dise˜no en Ingenier´ıa Integraci´on semi-anal´ıtica de matrices de rigidez de elementos finitos en simetr´ıa axial I.J. Lozada Decanato de Ingenier´ıa Civil, Universidad Centroccidental Lisandro Alvarado Apto. Postal 3001, Departamento de Ciencias B´asicas, Barquisimeto, Venezuela email: [email protected] D.V. Griffiths Geomechanics Research Center, Colorado School of Mines, Golden CO 80401, Colorado, USA email: [email protected] M. Cerrolaza Instituto Nacional de Bioingenier´ıa, Universidad Central de Venezuela, Ciudad Universitaria S/N, Caracas, Venezuela email: [email protected] Resumen El m´etodo de los elementos finitos (MEF) es una t´ecnica de an´alisis num´erico que permite la obtenci´on de soluciones aproximadas de un amplio rango de problemas de ingenier´ıa. Estos problemas pueden ser representados matem´aticamente por medio de ecuaciones diferenciales y en muy pocos casos, es posible obtener la soluci´on anal´ıtica de las mismas. Las alternativas que se presentan son b´asicamente dos: considerar hip´otesis simplificativas que eliminen las dificultades y tornen el problema hacia una soluci´on anal´ıtica posible; o mantener las complejidades del problema y obtener soluciones num´ericas aproximadas haciendo uso del avance de las computadoras y de los algoritmos num´ericos. En este trabajo, se presenta una metodolog´ıa general aplicada a la integraci´on simb´olica de las matrices de rigidez de elementos cuadril´ateros de cuatro y ocho nodos, en problemas de simetr´ıa axial. Esta metodolog´ıa permite calcular los elementos de la matriz de rigidez de forma sencilla, mediante una expresi´on semi-anal´ıtica. Su aplicaci´on mejora los tiempos de CPU, en comparaci´on con la t´ecnica de integraci´on num´erica tipo Gaussiano. El sustituir las subrutinas num´ericas utilizadas en los programas comerciales de elementos finitos por las subrutinas simb´olicas proporcionar´a a los profesionales e investigadores que trabajan con el m´etodo de los elementos finitos una reducci´on sustancial de tiempos de CPU, en el an´alisis de grandes estructuras. Palabras clave: Elemento cuadril´atero axisim´etrico, integraci´on semi-anal´ıtica, matriz de rigidez, manipulaci´on simb´olica. SEMI-ANALYTICAL INTEGRATION OF FINTE ELEMENTS STIFFNESS MATRICES IN AXISYMMETRIC PROBLEMS Summary The finite element method (MEF) is a numerical analysis technique which provide approximated solutions for engineering problems. These problems can be represented mathematically through differential equations and, in a few cases, it is possible to obtain their analytical solution. This work discusses a general methodology based on the symbolic integration of stiffness matrices of quadrilateral elements having four and eight nodes in axisymmetric problems. The terms of the stiffness matrix are calculated through a semi analytical expression, thus leading to a reduction of CPU time, when compared with numerical integration techniques. The semi analytical subroutines developed herein will help FEM researchers and engineers, by providing significant reductions of CPU times in the analysis of large models. Keywords: Quadrilateral axisymmetric finite elements, semi-analytical integration, stiffness matrix, symbolic manipulation. c °Universitat Polit`ecnica de Catalunya (Espa˜na). ISSN: 0213–1315 Recibido: Marzo 2008 Aceptado: Mayo 2008 346 I.J. Lozada, D.V. Griffiths y M. Cerrolaza INTRODUCCI´ ON En el campo de la ingenier´ıa y las ciencias f´ısicas, el permanente desarrollo computacional permite la incorporaci´on de nuevas tecnolog´ıas que mejoren las existentes. El c´alculo de matrices de rigidez y masa en el m´etodo de los elementos finitos, ha evolucionado exigiendo dicha incorporaci´on. Hasta 1980, en la integraci´on de matrices de rigidez de elementos finitos cuadril´ateros, no se contaba con metodolog´ıas para calcular de forma expl´ıcita sus coeficientes, salvo en casos espec´ıficos (rect´angulos y cuadril´ateros convexos). Fue Okabe8, uno de los primeros investigadores que present´o f´ormulas expl´ıcitas de integrales de expresiones racionales sobre un elemento cuadril´atero convexo isoparam´etrico de cuatro nodos. Estas f´ormulas requieren la evaluaci´on de productos de la forma φ(αβ)ln[f(α,β)]/(αmβn), donde φyfson funciones racionales y m≥0, n≥0. Para los elementos cuadril´ateros de lados rectos, que usualmente se encuentran en las aplicaciones pr´acticas, los par´ametros |α| y|β|son peque˜nos y haciendo que la evaluaci´on de los t´erminos que contienen logaritmos se dificulte. Este trabajo no arroj´o comparaciones de tiempos. Autores como Babu y Pinter1, basados en el trabajo de Okabe, presentaron f´ormulas de integraci´on semi-anal´ıtica, para evaluar integrales sobre elementos cuadril´ateros sim´etricos de lados rectos, inscritos en un c´ırculo, mejorando la exactitud de los resultados obtenidos con integraci´on gaussiana. Por otra parte, Mizukami7, trabaj´o con paralelogramos y present´o f´ormulas de integraci´on anal´ıtica para el c´alculo de la matriz de rigidez de este elemento. En este tipo de elemento, el jacobiano de la transformaci´on de coordenadas es una funci´on constante, lo cual ayud´o al autor a deducir las f´ormulas presentadas en su trabajo. No se realiz´o comparaci´on alguna de tiempos de CPU. As´ı mismo, Rathod10, generaliza los resultados obtenidos por Babu y Pinter1y Mizukami7, para presentar f´ormulas de integraci´on anal´ıtica para un elemento finito cuadril´atero isoparam´etrico de cuatro nodos. Para ello, Rathod se bas´o en m´etodos b´asicos de integraci´on (integraci´on por partes), transformando todas las integrales involucradas en el c´alculo de la matriz de rigidez a integrales unidimensionales, que luego son expresadas como combinaci´on lineal de cuatro integrales b´asicas. Usando el sistema de ´algebra computacional REDUCE11 desarrollado en 1987, Kikuchi4obtuvo f´ormulas explicitas para la integraci´on de la matriz de rigidez del elemento cuadril´atero de cuatro nodos. Este autor no report´o reducciones en el costo computacional. Yagawa et al.14 hacen una combinaci´on de m´etodos anal´ıticos y num´ericos para integrar la matriz de rigidez de un elemento cuadril´atero isoparam´etrico de cuatro nodos, en elasticidad plana. Dichos autores, con la ayuda del sistema de ´algebra computacional REDUCE11 expandieron y agruparon convenientemente el integrando. Posteriormente, efectuaron la integraci´on num´erica, logrando una mejora en los tiempos de c´alculo de un 15 % en comparaci´on a la integraci´on num´erica Gauss con 4 puntos. Griffiths3present´o una f´ormula semi-anal´ıtica, para el c´alculo de la matriz de rigidez de un elemento finito cuadril´atero isoparam´etrico de cuatro nodos, en problemas de elasticidad bidimensional. Este autor clasific´o los elementos de la matriz de rigidez en seis grupos, cada uno de ellos caracterizado por una condici´on particular de sus grados de libertad. La expresi´on semi-anal´ıtica fue obtenida manipulando simb´olicamente con el software matem´atico MAPLE6la f´ormula de integraci´on de gauss con 4 puntos. Este autor utiliz´o la f´ormula semi-anal´ıtica y transformaciones de coordenadas, para calcular las componentes de la matriz de rigidez. La t´ecnica desarrollada report´o una mejora sustancial de los tiempos de CPU, en comparaci´on con la obtenida al usar el m´etodo de integraci´on num´erica Gauss con 4 puntos. Otro trabajo elaborado en esta l´ınea, fue el de Videla et al.12. Dichos autores presentaron f´ormulas anal´ıticas para la integraci´on de la matriz de rigidez de un elemento cuadril´atero isoparam´etrico de cuatro nodos, en problemas de elasticidad plana. Estos autores manipularon simb´olicamente todas las expresiones obtenidas en el c´alculo de las derivadas parciales de las funciones de forma con respecto a las coordenadas cartesianas y obtuvieron formas generales para cada una de ellas. Estas expresiones permitieron Integraci´on semi-anal´ıtica de matrices de rigidez de elementos finitos en simetr´ıa axial 347 efectuar la integraci´on de una manera m´as sencilla. Los resultados fueron codificados en Fortran y reportaron una reducci´on de un 50 % en los c´alculos al compararlos con la integraci´on num´erica (Gauss con dos y tres puntos). Lozada et al.5generalizaron los resultados obtenidos por Griffiths, para ser aplicados en la integraci´on semi-anal´ıtica de las matrices de rigidez de los elementos finitos cuadril´ateros de ocho nodos en problemas elasticidad bidimensional, obteniendo mejoras de un 37 % en los tiempos de c´alculo en comparaci´on a la integraci´on num´erica (Gauss con 4 puntos). Videla et al.13 presentaron f´ormulas explicitas para la integraci´on de la matriz de rigidez del elemento finito cuadril´atero de ocho nodos, en problemas de elasticidad plana. Este autor logr´o una reducci´on en los tiempos de c´omputo de un 50 % en comparaci´on a la integraci´on num´erica. En este trabajo se generalizan los resultados obtenidos por Griffiths3y Lozada et al5para ser aplicados en la integraci´on simb´olica del elemento finito cuadril´atero de cuatro y ocho nodos, en problemas de simetr´ıa axial. FORMULACI´ ON Los s´olidos tridimensionales de simetr´ıa axial o s´olidos de revoluci´on, pueden ser analizados como problemas bidimensionales. En este trabajo, se considera el elemento cuadril´atero de revoluci´on de cuatro y ocho nodos, con dos grados de libertad por nodo y enumerado en sentido anti-horario. El elemento, como puede apreciarse en la Figura 1, es un cuadril´atero que genera un toroide al girar alrededor del eje z. 3 2 1 Z (r , z ) 1 1 (r , z ) 2 2 (r , z ) 3 3 (r , z ) 4 4 r 4 3 2 1 4 Figura 1. Elemento s´olido de revoluci´on (cuadril´atero) Mediante la transformaci´on de coordenadas, de la ecuaci´on (1), el elemento es transformado en un elemento cuadril´atero con lados paralelos a los ejes coordenados, como se muestra en la Figura 2                  η  ξ Figura 2. Cambio de coordenadas 348 I.J. Lozada, D.V. Griffiths y M. Cerrolaza La transformaci´on de coordenadas que relaciona el plano rz con el plano ς η est´a dada por: r= n X i=1 Ni(ξ, η)ri;z= n X i=1 Ni(ξ, η)zi(1) donde nrepresenta el n´umero de nodos del elemento y las Nison las funciones que interpolan los desplazamientos y la geometr´ıa del elemento. En este trabajo se consideran las siguientes funciones de interpolaci´on para el elemento cuadril´atero de cuatro y ocho nodos. Cuatro nodos: N1=1 4(1 −η) (1 −ξ) ; N2=1 4(1 −η) (1 + ξ) N3=1 4(1 + η) (1 + ξ) ; N4=1 4(1 + η) (1 −ξ) (2) Ocho nodos: N1=−1 4(1 −ξ)(1 −η)(ξ+η+ 1); N2=−1 4(1 + ξ)(1 −η)(−ξ+η+ 1) N3=−1 4(1 + ξ)(1 + η)(−ξ−η+ 1); N4=−1 4(1 −ς)(1 + η)(ς−η+ 1) N5=1 2(1 −η)(1 −ξ2) ; N6=1 2(1 + ξ)(1 −η2) (3) N7=1 2(1 + η)(1 −ξ2); N8=1 2(1 −ξ)(1 −η2) As´ı: Kij =ZZZ e BT iD BjdV (4) =ZZ A 2π Z 0 BT iD Bjr dθ dA (5) = 2πZZ A BT iD Bjr dA (6) = 2π 1 Z −1 1 Z −1 r BT i(ξ, η)DBj(ξ, η) det Jdξdη (7) donde: Bt i=·∂Ni ∂r 0Ni r ∂Ni ∂z 0∂Ni ∂z 0∂Ni ∂r ¸;D=E(1 −ν) (1 + ν)(1 −2ν)   1−ν ν ν 0 ν1−ν ν 0 ν ν 1−ν0 0001−2ν 2   (8) J,Eyνdenotan el jacobiano de la transformaci´on, el m´odulo de Young y el coeficiente de Poisson, respectivamente. Integraci´on semi-anal´ıtica de matrices de rigidez de elementos finitos en simetr´ıa axial 349 De (7) y (8) se tiene: Kij = 2π 1 Z −1 1 Z −1·cij 11 cij 12 cij 21 cij 22 ¸1 det Jdξdη (9) Siendo: cij 11 =    NiNj n P i=1 riNi +Ãn X i=1 riNi!TiTj    E1+ (TiNj+TjNi)E2+Ãn X i=1 riNi!SiSjE3(10) cij 12 =ÃÃ n X i=1 riNi!SjTi+SjNi!E2+Ãn X i=1 riNi!SiTjE3(11) cij 21 =ÃÃ n X i=1 riNi!SiTj+SiNj!E2+Ãn X i=1 riNi!TiSjE3(12) cij 22 =Ãn X i=1 riNi!(SiSjE1+TiTjE3) (13) donde: K=E (1 + ν) (1 −2ν);E1=K(1 −ν) ; E2=Kν;E3=K(1 −2ν) 2(14) Ti=µ∂z ∂η ∂Ni ∂ξ −∂z ∂ξ ∂Ni ∂η ¶= det Jµ∂Ni ∂r ¶(15) Si=µ−∂r ∂η ∂Ni ∂ξ +∂r ∂ξ ∂Ni ∂η ¶= det Jµ∂Ni ∂z ¶(16) Generalmente, la integraci´on exacta de los t´erminos de la ecuaci´on (9) es compleja y los programas de elementos finitos utilizan integraci´on gaussiana para efectuarlas. Con integraci´on gaussiana de orden 2×2, se obtiene la expresi´on semi-anal´ıtica: kij = 6πnA3(E1S1+E2S2+E3S3)+f1(E1S4+E2S5+E3S6) 3A2 3−f2 1 +A3(E1T1+E2T2+E3T3)+f2(E1T4+E2T5+E3T6) 3A2 3−f2 2o (17) donde: f1= (r1+r3)(z4−z2)−(z1+z3)(r4−r2)−2(r2z4−r4z2) (18) f2= (z2+z4)(r3−r1)−(r2+r4)(z3−z1)−2(r3z1−r1z3) (19) A3=1 8[(z2−z4)(r1−r3)+(z3−z1)(r2−r4)] Las funciones SiyTidependen del n´umero de nodos y de las coordenadas locales. 350 I.J. Lozada, D.V. Griffiths y M. Cerrolaza GENERACI´ ON DE LOS T´ ERMINOS DE LA MATRIZ DE RIGIDEZ En este desarrollo se considera la clasificaci´on presentada por Griffiths3para el elemento de cuatro nodos y la de Lozada et al.5para el elemento de ocho nodos, ver Tabla I y II. Grupos T´erminos Descripci´on Ak11,k22, k33,k44, k55, k66, k77,k88 GDL paralelos en el mismo nodo Bk12,k34, k56,k78 GDL ortogonales en el mismo nodo Ck13, k35, k57, k17, k24, k46, k68, k28 GDL paralelos en nodos adyacentes Dk23, k36, k67, k27, k45, k58, k18, k14 GDL ortogonales en nodos adyacentes Ek15, k37, k26, k48 GDL paralelos en nodos opuestos Fk16, k38, k25, k47 GDL perpendiculares en nodos opuestos Tabla I. Clasificaci´on de los t´erminos de la matriz de rigidez del elemento de cuatro nodos GRUPO T´ ERMINOS DESCRIPCI´ ON GRADO DE ADYACENCIA Ak1,1,k2,2, k3,3,k4,4, k5,5, k6,6, k7,7,k8,8, k9,9,k10,10, k11,11,k12,12, k13,13, k14,14, k15,15,k16,16 GDL paralelos en el mismo nodo 0 Bk1,2,k3,4, k5,6,k7,8, k9,10, k11,12, k13,14, k15,16 GDL ortogonales en el mismo nodo 0 Ck1,3, k3,5, k5,7, k1,7, k9,11, k11,13, k13,15, k9,15, k2,4, k4,6, k6,8, k2,8, k10,12, k12,14, k14,16, k10,16 GDL paralelos de nodos separados por un nodo 2 Dk2,3, k3,6, k6,7, k2,7, k10,11, k11,14, k14,15, k10,15, k4,5, k5,8, k1,8, k1,4, k12,13, k13,16, k9,16, k9,12 GDL ortogonales de nodos separados por un nodo 2 Ek1,5, k9,13, k3,7, k11,15, k2,6, k4,8, k10,14, k12,16 GDL paralelos en nodos opuestos 4 Fk1,6, k9,14, k3,8, k11,16, k2,5, k10,13, k4,7, 12,15 GDL perpendiculares en nodos opuestos 4 Gk1,9, k3,9, k3,11, k5,11, k5,13, k7,13, k7,15, k1,15, k2,10, k4,10, k4,12, k6,12, k6,14, k8,14, k8,16, k2,16 GDL paralelos de nodos adyacentes 1 Hk1,10, k3,10, k3,12, k5,12, k5,14, k7,14, k7,16, k1,16, k2,9, k4,9, k4,11, k6,11, k6,13, k8,13, k8,15, k2,15 GDL perpendiculares de nodos adyacentes 1 Ik1,11, k7,11, k7,9, k5,9, k5,15, k3,15, k3,13, k1,13, k2,12, k8,12, k8,10, k6,10, k6,16, k4,16, k4,14, k2,14 GDL paralelos de nodos separados por dos nodos 3 Jk2,11, k8,11, k8,9, k6,9, k6,15, k4,15, k4,13, k2,13, k1,12, k7,12, k7,10, k5,10, k5,16, k3,16, k3,14, k1,14 GDL paralelos de nodos separados por dos nodos 3 Tabla II. Clasificaci´on de los t´erminos de la matriz de rigidez del elemento de ocho nodos Integraci´on semi-anal´ıtica de matrices de rigidez de elementos finitos en simetr´ıa axial 351 En la clasificaci´on anterior se considera la simetr´ıa de la matriz de rigidez del elemento analizado. El elemento de ocho nodos fue clasificado en 10 grupos, considerando la adyacencia entre los grados de libertad (GDL) de los nodos del elemento y una condici´on particular de paralelismo o perpendicularidad entre sus grados de libertad. En la Figura 3, se muestran las diferentes relaciones de adyacencia del nodo 3, con los otros nodos del elemento y los grados de libertad en cada nodo. 1 0 2 1 3 2 4 3                                         Figura 3. Grados de libertad en cada nodo y relaciones de adyacencia con el nodo3 Para generar los elementos de la matriz de rigidez, se utilizan transformaciones de coordenadas sencillas que no modifiquen la geometr´ıa del elemento sino su posici´on, como se muestra en la Tabla III Nodos Transformaci´on ( R) (r1,z1) (r4,z4) (r2,z2) (r1,z1) (r3,z3) (r2,z2) (r4,z4) (r3,z3) Tabla III. Transformaci´on de coordenadas Dado un elemento cualquiera de un grupo (t´ermino padre), los elementos restantes de este pueden ser generados utilizando las funciones SiyTi(funciones generadoras) y transformaciones de coordenadas. Comparando el n´umero de operaciones algebraicas realizadas por las funciones generadoras de cada elemento de un grupo, resulta computacionalmente m´as eficiente generar cada grupo con varios elementos padre, debido a que el n´umero de operaciones algebraicas realizadas por las funciones generadoras SiyTide cada elemento del grupo var´ıan. A continuaci´on se muestra la forma como se generan los t´erminos de la matriz de rigidez en cada grupo, utilizando la ecuaci´on (17) y la transformaci´on T1. Elementos de ocho nodos Grupo A. (T´erminos padres: k1,1,k2,2,k9,9yk10,10) k1,1 R −→k7,7 R −→k5,5 R −→k3,3 k2,2 R −→ k8,8 R −→ k6,6 R −→4,4(20) k9,9 R −→k15,15 R −→k13,13 R −→k11,11 k10,10 R −→k16,16 R −→k14,14 R −→k12,12 352 I.J. Lozada, D.V. Griffiths y M. Cerrolaza Grupo B.(T´erminos padres k1,2yk9,10) k1,2 R −→k7,8 R −→k6,5 R −→k4,3 k9,10 R −→ k15,16 R −→ k13,14 R −→ k11,12 (21) Grupo C.(T´erminos padres: k1,3,k6,8, k9,11 y k10,12 ) k1,3 R −→k1,7 R −→k5,7 R −→k3,5 k6,8 R −→ k6,4 R −→ k2,4 R −→ k2,8(22) k9,11 R −→k9,15 R −→k13,15 R −→k11,13 k10,12 R −→k10,16 R −→k14,16 R −→k12,14 Grupo D.(T´erminos padres: k1,4,k1,8,k9,12 yk14,15) k1,4 R −→k2,7 R −→k5,8 R −→k3,6 k1,8 R −→ k6,7 R −→ k4,5 R −→ k2,3(23) k9,12 R −→k10,15 R −→k13,16 R −→k11,14 k14,15 R −→k12,13 R −→k10,11 R −→k9,16 Grupo E.(T´erminos padres: k1,5,k2,6,k9,13 yk12,16) k1,5 R −→k3,7 k2,6 R −→ k4,8(24) k9,13 R −→k11,15 k12,16 R −→k10,14 Grupo F.(T´erminos padres: k1,6yk9,14) k1,6 R −→k4,7 R −→k2,5 R −→k3,8 k9,14 R −→ k12,15 R −→ k10,13 R −→ k11,16 (25) Grupo G.(T´erminos padres: k7,15,k7,13,k2,10 yk6,12) k7,15 R −→k5,13 R −→k3,11 R −→k1,9 k7,13 R −→ k5,11 R −→ k3,9 R −→ k1,15 (26) k2,10 R −→k8,16 R −→k6,14 R −→k4,12 k6,12 R −→k4,10 R −→k2,16 R −→k8,14 Integraci´on semi-anal´ıtica de matrices de rigidez de elementos finitos en simetr´ıa axial 353 Grupo H.(T´erminos padres k2,9k6,11 k1,10yk5,12) k2,9 R −→k8,15 R −→k6,13 R −→k4,11 k6,11 R −→ k4,9 R −→ k2,15 R −→ k8,13 (27) k1,10 R −→k7,16 R −→k5,14 R −→k3,12 k5,12 R −→k3,10 R −→k1,16 R −→k7,14 Grupo I.(T´erminos padres k1,11,k5,9,k2,12yk6,10) k1,11 R −→k7,9 R −→k5,15 R −→k3,13 k5,9 R −→ k3,15 R −→ k1,13 R −→ k7,11 ; (28) k2,12 R −→k8,10 R −→k6,16 R −→k4,14 k6,10 R −→k4,16 R −→k2,14 R −→k8,12 Grupo J.(T´erminos padres k2,11,k2,13,k1,12yk1,14) k2,11 R −→k8,9 R −→k6,15 R −→k4,13 k2,13 R −→ k8,11 R −→ k6,9 R −→ k4,15 (29) k1,12 R −→k7,10 R −→k5,16 R −→k3,14 k1,14 R −→k7,12 R −→k5,10 R −→k3,16 Los t´erminos de la matriz del elemento de cuatro nodos son generados de manera an´aloga al caso anterior. Las funciones S1,S2,S3,S4,S5,S6,T1,T2,T3,T4,T5yT6utilizadas en cada grupo son las que generan los t´erminos padres. Por ejemplo, las funciones que generan ak22 en el elemento de ocho nodos son: S1= (1/96) (r42)2(7 r1+2 rs24+r3) S2= 0 S3= (1/96) (z42)2(7 r1+2 rs24 +r3) S4= (1/96) (r42)2(4 r1+rs24) S5= 0 S6= (1/96)(z42)2(4r1+rs24) (30) T1= (1/2592) (2 r3 1+7 r3 2+8 r3 3+12 r2 3r4-18 r3r2 4+7 r3 4+6 r2 1(r2+r43)-3 r2 2(6 r3+r4)+ 3r2(4 r2 3+2 r3r4-r2 4)+6 r1(2 r2 2-r2(3 r3+r4)+r4(-3 r3+2 r4))) T2= 0 T3= (1/2592) (2 (z2 1+4 z3z31-z2z4) (r1+2 rs24+r3)+2 (rs13+5 r4-r2) (z1-2 z3)z4+2 (-r1-5 r2+r43) (2 z3-z1)z2+(2 rs13+r4+7 r2)z2 2+(2 rs13+7 r4+r2)z2 4) T4= (-1/7776) r42 (12 (r2 1+r2 2+r2 3+r2 4)+21 r1rs42+6 r2r4-30 r1r3-33 r3rs42) T5= 0 T6= (1/2592) (2 (3 rs24+rs13) (2 z3z4-z1z4)+2 (6 r4-r2)z2z4+2 z2(z1-2 z3) (rs14+r3+3 r2)+8 r42 z1z3+(rs13+4 r2)z2 2-2 r42 z2 1-8 z2 3r42-z2 4(rs13+4 r4))