• Author / Uploaded
• Anonymous QUJMWr

Citation preview

8QLYHUVLGD )DFXOWD 

7HFQROyJLF 5HJLRQD 




1DFLRQDO

%XHQR 

$LUHV

CURSO DE ESPECIALIZACIÓN

,QWURGXFFLy

D 0pWRG 

)RUPXODFLy



G 

9DULDFLRQD G 

(OHPHQWR 



(OHPHQWR



)LQLWRV

)LQLWRV

Dr. Ing. Claudio E. Jouglard Notas del Curso dictado en el 2º cuatrimestre de 2002

,QGLFH

FORMULACIÓN VARIACIONAL DE ELEMENTOS FINITOS ...............................

1

1 Metodos aproximados de solución .................................................................... 1.1 Método de Rayleigh-Ritz ............................................................................... 1.1.1 Ejemplo de aplicación del método de Rayleigh-Ritz ...................... 1.2 Métodos de residuos ponderados ................................................................. 1.2.1 Aproximación mediante residuos ponderados ......... ...................... 1.2.2 Aplicación a la ecuación diferencial de Poisson .............................. 2 Funciones de forma ...................................... ..................................................... 2.1 Orden de continuidad ................................................................................... 2.2 Polinomios de Lagrange en una y dos dimensiones .................................... 2.3 Polinomios de Hermite en una y dos dimensiones ....................................... 2.4 Convergencia de la aproximación. Test de la parcela .................................. 3 Implementación matricial del método de elementos finitos ................................. 3.1 Expresión matricial de la energía potencial total .......................................... 3.2 Aproximación por elementos finitos .............................................................. 3.3 Matriz gradiente ................................. .......................................................... 3.4 Matriz de rigidez y vector de cargas nodales equivalentes ......................... 3.5 Ejemplo de ensamblaje ................................................................................. 4 Elemento finito triangular .................................................. ................................. 4.1 Coordenadas de área ................................................................................... 4.2 Campos de interpolación para triángulos ...................................................... 5 Formulación isoparamétrica ................................................................................ 5.1 Variación de funciones en los elementos curvilíneos ................................... 5.1.1 Evaluación de matrices del elemento .............................................. 5.2 Integración numérica ..................................................................................... 5.2.1 Cuadratura de Gauss-Legendre ...................................................... 5.2.2 integración completa y reducida....................................................... Referencias ............................................................................................................

1 1 3 6 6 7 9 12 13 15 17 18 18 19 21 22 23 25 26 27 28 31 31 33 33 34 35

)RUPXODFLy

9DULDFLRQD G




(OHPHQWR 

)LQLWRV

La formulación de elementos finitos puede deducirse para ciertos problemas, como por ejemplo el análisis de estructuras, como una extensión de los métodos matriciales utilizados para calcular estructuras de vigas y reticulados. Sin embargo, dicha deducción encuentra serias limitaciones cuando se quiere extender la formulación a problemas no estructurales. Por ello se mostrarán en este apunte algunos conceptos básicos de la formulación variacional del método de elementos finitos que pueden aplicarse a una gran variedad de problemas. En primer lugar describiremos algunos conceptos sobre métodos aproximados de solución para ecuaciones diferenciales, en particular veremos el método de Rayleigh-Ritz y el método de residuos ponderados. Luego veremos la utilización de estos métodos con elementos finitos y se describirá la implementación matricial y los elementos más utilizados. 1 MÉTODOS APROXIMADOS DE SOLUCIÓN Describiremos los métodos más usuales empleados para la resolución aproximada de ecuaciones diferenciales en dos y tres dimensiones, sobre los que basan la mayoría de implementaciones de elementos finitos. En general los métodos empleados son de dos tipos, por un lado están los métodos basados en principios variacionales, generalmente asociados a la minimización de algún funcional, y por otro lado tenemos los métodos del tipo de residuos ponderados que se aplican directamente sobre la ecuación diferencial y sus condiciones de contorno, no precisando de ningún funcional asociado. 1.1 MÉTODO DE RAYLEIGH-RITZ Con este método es posible obtener soluciones aproximadas de ecuaciones diferenciales mediante principios variacionales. La idea básica consiste en aproximar a las soluciones u, v que hacen estacionario un funcional mediante una suma ponderada de funciones m

u~ = ∑ ai N i ( x, y ), i =1

v~ =

n

∑a

i

N i ( x, y )

(1)

i = m +1

donde ai son constantes a determinar llamadas coordenadas generalizadas. Las funciones Ni (x, y) son llamadas funciones de prueba y pueden ser elegidas arbitrariamente, pero deben ser admisibles, esto es, deben satisfacer las condiciones esenciales de contorno y las condiciones de compatibilidad. En general, se utilizan polinomios, aunque pueden utilizarse funciones trigonométricas ú otro tipo de funciones. Este método fue utilizado por primera vez por Lord Rayleigh en 1870 usando un campo de aproximación con una única función de prueba. Posteriormente, en 1909, Ritz generalizó el método construyendo un campo de aproximación con varias funciones.

2

Formulación de Elementos Finitos

Consideremos un sistema de ecuaciones diferenciales cuya solución es equivalente a hacer estacionaria la primera variación de un funcional asociado Π que es función de u, v y sus derivadas primeras

Π = ∫ F (x, y, u , v, u,x , u,y , v ,x , v ,y )dx dy Ω

Las derivadas de las funciones de aproximación u son m m ∂u~ ∂N i ( x , y ) ∂u~ ∂N i ( x , y ) , = ∑ ai = ∑ ai ∂x i =1 ∂x ∂y i =1 ∂y

(2)

(3)

y en forma análoga para la aproximación de v. Si substituimos las funciones de aproximación u, v y sus derivadas en el funcional Π este se transformará en una función de las coordenadas generalizadas ai, cuyos valores por ahora desconocemos, esto es

Π = Π (ai )

i = 1,2,  , n

(4)

Luego deseamos conocer cuales son los mejores valores de las constantes ai, tal que reemplazados en las expresiones (1) nos brinden la mejor aproximación a la solución del sistema de ecuaciones diferenciales asociado al funcional Π. Para ello aplicamos la condición de estacionaridad a este funcional, que sabemos que debe ser satisfecha por la solución exacta, resultando n

δΠ = ∑ i =1

∂Π δai = 0 ∂a i

(5)

Como esta ecuación debe ser válida para variaciones arbitrarias δai, luego deben cumplirse las siguientes n ecuaciones algebraicas: ∂Π =0 ∂a i

i = 1,2,  , n

(6)

Si ahora consideramos el caso particular, pero muy común en problemas físicos, donde el funcional Π es una función cuadrática de las funciones u, v y sus derivadas primeras, entonces al substituir las aproximaciones dadas por la ec. (1) este funcional será una función cuadrática de las coordenadas generalizadas ai. Por lo tanto, las derivadas de este funcional serán funciones lineales en las coordenadas generalizadas ai que se pueden expresar como: ∂Π = (k i1a1 + k i 2 a 2 +  + k in a n − f i ) = 0 ∂a i

i = 1, 2,  , n

(7)

En forma matricial este sistema de ecuaciones se puede escribir como  k11 k  21     k n1

k12 k 22  kn2

 k1n   a1   f1   k 2 n  a 2   f 2   =            k nn  an   f n 

(8)

y en forma abreviada usando notación matricial Ka = f

(9)

Introducción al Método de Elementos Finitos

3

La solución de este sistema de ecuaciones nos dá los valores de las constantes ai. Para un número grande ecuaciones la solución a mano de este sistema es prácticamente imposible y debe ser resuelto por computadora. Notemos que las derivadas segundas del funcional Π respecto de las constantes ai valen ∂Π = k ij ∂a i ∂a j

(10)

y dado que Π es una función continua en las coordenadas generalizadas ai, por ser un polinomio de segundo grado en estas coordenadas, entonces deben ser iguales las derivadas segundas cruzadas ∂Π ∂Π = ∂a i ∂a j ∂a j ∂ a i

⇒

k ij = k ji

(11)

Por lo tanto la matriz K es simétrica. Esto tiene implicancias desde el punto de vista computacional, ya que para almacenar esta matriz solo precisamos guardar la mitad de sus coeficientes. Una vez conocidas las constantes ai reemplazándolas en las expresiones (1) tenemos las aproximaciones buscadas a la solución del problema variacional y equivalentemente al sistema de ecuaciones diferenciales asociado. 1.1.1

Ejemplo de aplicación del método de Rayleigh-Ritz.

Consideremos una viga de eje recto de sección constante y longitud L, simplemente apoyada en sus extremos y sometida a un cargamento distribuido p de forma senoidal p = p0 sen πx/L. p = p 0 sen

πx L x v

L y, v

Figura 1. Viga de eje recto y sección constante.

El funcional a minimizar es el de la energía potencial total V, que se compone de la suma de la energía de deformación U y el potencial de las fuerzas externas Ve V = U + Ve

(12)

La energía de deformación U para la viga, despreciando las deformaciones por corte, es 2

1 L  d 2v  U = ∫ E I  2  dx 2 0  dx 

(13)

Siendo v el desplazamiento del eje centroidal de la viga, E es el módulo de elasticidad e I es el momento de inercia de la sección transversal.

4

Formulación de Elementos Finitos El potencial de fuerzas externas es L

Ve = − ∫ p v dx 0

(14)

Luego la energía potencial total es 2

L 1 L  d 2v  V = ∫ E I  2  dx − ∫ p v dx 0 2 0  dx 

(15)

Si adoptamos la siguiente aproximación polinómica v~ para los desplazamientos verticales v~ = a (L x − x 2 ) (16) Donde a es un coeficiente a determinar. Notemos que esta aproximación satisface las condiciones esenciales de contorno, esto es v~ (0 ) = 0 (17) ~ v (L ) = 0 Además, las derivadas de la aproximación v~ valen d v~ = a (L − 2 x ) dx d 2 v~ = −2 a dx 2

(18)

Substituyendo en la expresión de la energía potencial total tenemos V (a ) =

L 1 L πx  2 2 E I a dx 2 ( ) −  p0 sen  a (Lx − x )dx ∫ ∫ 0 0 L  2 

(19)

Calculando las integrales resulta  4 L3  V ( a ) = 2 E Ia 2 L − a  3  p0  π 

(20)

Minimizando la energía potencial total  4 L3  dV = 4 E I a L −  3  p0 = 0 da  π 

(21)

resultando a=

p0 L2 π3 E I

(22)

Por lo tanto la solución aproximada v~ queda p L4 v~ = 30 π EI Si comparamos con la solución exacta

 x   x  2    −     L   L  

(23)

Introducción al Método de Elementos Finitos

v=

p0 L4 πx sen 4 π EI L

5

(24)

en el punto medio x = L/2 tenemos 1 p0 L4 p0 L4 v~ ( L2 ) = = 0 . 008063 4 π3 E I EI v ( L2 ) =

p0 L4 π4 E I

= 0.010266

(25)

p0 L4 EI

Notemos que la solución aproximada nos da un error del 21 % por debajo para la flecha en el punto medio. Esto es, la solución aproximada es más rígida que la solución exacta dando deformaciones menores. Los momentos flectores se relacionan con las derivadas segundas de los desplazamientos M = − EI

d 2v dx 2

(26)

Si ahora comparamos los momentos en el punto medio p0 L2 ~ L M ( 2 ) = 2 3 = 0.0645 p0 L2 π p L2 M ( L2 ) = 0 2 = 0.1013 p0 L2 π

(27)

Notemos que en este caso el momento dado por la solución aproximada nos da un error del 36 % por debajo para el momento flector en el punto medio. Esto es, el error en el momento es mayor que el de la aproximación de los desplazamientos. Además notemos que la solución aproximada nos da un momento constante para toda la viga ~ M ( x ) = EI 2 a (28) Mientras la solución exacta varia sinusoidalmente M ( x) =

p0 L2 πx sen 2 π L

(29)

En este caso hemos aplicado en realidad el método de Rayleigh, pues hemos utilizado una única función de aproximación. Para aplicar Rayleigh-Ritz deberíamos utilizar al menos dos funciones de aproximación. Por ejemplo, podríamos utilizar además de la parábola cuadrática usada en el ejemplo anterior, un polinomio de cuarto grado como v~1 = a1 (L x − x 2 ) (30) v~2 = a2 (x 4 − 2 x 3 L + x 2 L ) Notemos que ambas aproximaciones satisfacen las condiciones esenciales de contorno y son simétricas respecto del punto medio. Se deja como ejercicio para el lector la determinación de los coeficientes a1, a2.

6

Formulación de Elementos Finitos

1.2 MÉTODOS DE RESIDUOS PONDERADOS Con estos métodos es posible obtener soluciones aproximadas de ecuaciones diferenciales que no tienen un funcional asociado. Consideremos una ecuación diferencial definida sobre un dominio Ω como A(u ) − f = 0

en Ω

(31)

y estando sujeta a condiciones naturales de contorno sobre la parte ΓN de su frontera de la forma M (u ) − g = 0

en ΓN

(32)

siendo A(u), M(u) operadores diferenciales. Algunos ejemplos de los operadores diferenciales A(u) son: A(u ) = − A(u ) =

d  du  a  + cu dx  dx 

d2 dx 2

(33)

 d 2u   b 2   dx 

 ∂  ∂u  ∂  ∂ u     A(u ) = −   k x  +  k y  ∂ x  ∂x  ∂ y  ∂y   Si reemplazamos la solución exacta u por una solución aproximada u~ en la ecuación diferencial (31) y en sus condiciones naturales de contorno (32), estas no serán satisfechas exactamente, generando un residuo RΩ en el dominio y un residuo RΓ en el contorno, esto es (34) R (u~ ) = A(u~ ) − f ≠ 0 en Ω Ω

RΓ (u~ ) = M (u~ ) − g ≠ 0

en ΓN

La idea básica de los métodos de residuos ponderados es imponer la condición

∫

Ω

W RΩ (u ) dΩ + ∫ W RΓ (u ) dΓ = 0 ΓN

(35)

para cualquier par de funciones arbitrarias W , W que sean integrables y no idénticamente nulas. Se puede demostrar que si esto se cumple entonces la función u debe ser la solución exacta de la ecuación diferencial y de sus condiciones de contorno naturales. 1.2.1

Aproximación mediante residuos ponderados.

A partir de esta idea es posible generar soluciones aproximadas de la forma n

u~ = ∑ ai N i ( x, y )

(36)

i =1

donde las funciones de prueba Ni (x, y) satisfacen las condiciones esenciales de contorno y las n constantes ai son coeficientes a determinar imponiendo las n condiciones

∫

Ω

Wi RΩ (u~ ) dΩ + ∫ Wi RΓ (u~ ) dΓ = 0 ΓN

i = 1,2,  , n

donde las funciones Wi , Wi son llamadas funciones de ponderación.

(37)

Introducción al Método de Elementos Finitos

7

Dependiendo del tipo de funciones de ponderación empleadas tenemos diferentes métodos de aproximación: a) Método de Colocación por puntos: en este método se impone que el residuo sea nulo en n puntos (xi, yi) del dominio y de la parte del contorno donde se hayan impuesto condiciones naturales de contorno. (38) R (u~ , x , y ) = 0 i = 1,2,  , p Ω

i

i

RΓ (u~ , xi , yi ) = 0

i = p + 1, p + 2,  , n

Esto es equivalente a adoptar a las funciones de ponderación como las funciones delta de Dirac δ(x-xi, y-yi) que se definen como:

∫

Ω

f ( x , y ) δ ( x − xi , y − y i ) d Ω = f ( xi , y i )

(39)

b) Método de Colocación por subdominios: en este método se impone que el residuo sea nulo en n subregiones Ωi del dominio y Γi de la parte del contorno donde se hayan impuesto condiciones naturales de contorno.

∫ ∫

Ωi Γi

RΩ (u~ ) = 0

i = 1,2,  , p

RΓ (u~ ) = 0

i = p + 1, p + 2,  , n

(40)

En este caso las funciones de ponderación valen 1 en cada subdominio y 0 en el resto. c) Método de los cuadrados mínimos: en este método las constantes ai se determinan de la condición de hacer mínimo un funcional I ∂I =0 ∂a i

(41)

i = 1,2,  , n

siendo el funcional I definido como 2 I = ∫ (RΩ (u~ ) ) dΩ + ∫ Ω

ΓN

(RΓ (u~ ) )2

dΓ

(42)

En este caso las funciones de ponderación son iguales a sus respectivos residuos. d) Método de Galerkin: en este método se adoptan las funciones de ponderación iguales a las funciones de prueba resultando

∫

Ω

N i RΩ (u~ ) dΩ + ∫ N i RΓ (u~ ) dΓ = 0 ΓN

i = 1,2,  , n

(43)

Nótese que en todas estos métodos las funciones de prueba Ni deben tener derivadas definidas hasta el máximo orden que aparece en el residuo, esto es, del mismo orden que la ecuación diferencial. 1.2.2

Aplicación a la ecuación diferencial de Poisson.

La ecuación diferencial de Poisson aparece en un gran número de aplicaciones de ingeniería, incluyendo problemas de transferencia de calor, electrostática, mecánica de los fluidos y mecánica de los sólidos, entre otros. Esta ecuación viene definida como

8

Formulación de Elementos Finitos  ∂  ∂u  ∂  ∂u    k y   = Q −  kx + ∂ x ∂ x ∂ y ∂ y     

(44)

donde kx, ky son llamadas constantes del material y Q es llamado el término fuente. Si Q es nulo esta ecuación es llamada ecuación de Laplace. Las condiciones naturales de contorno son del tipo kx

∂u ∂u nx + k y ny = q ∂x ∂y

(45)

Aplicando el método de residuos ponderados tenemos:  ∂  ∂u  ∂  ∂ u   ∫Ω W − ∂x  k x ∂x  − ∂y  k y ∂y  − Q  dΩ + ∫ΓN W

 ∂u  ∂u  k x ∂x n x + k y ∂y n y − q  d Γ = 0  

(46)

Si integramos por partes la integral del residuo en el dominio resulta  ∂W ∂ u  ∂W ∂u ∫Ω k x ∂x ∂x + k y ∂y ∂y − W Q  dΩ  ∂u ∂u  nx + k y n y dΓ − ∫ W q dΓ = 0 + ∫ (W − W ) k x ΓN ΓN ∂y   ∂x

(47)

La primera integral es conocida como la forma variacional de la ecuación diferencial y notemos que solo tiene derivadas primeras, a diferencia de la forma original de residuos ponderados (46) que contenía derivadas segundas. También se suele llamar a esta formulación como formulación débil, por los menores requerimientos que deben cumplir las funciones de aproximación en cuanto al orden de las derivadas definidas. Por otro lado, si utilizamos esta formulación como base para una aproximación por el método de Galerkin, para cada función de prueba Ni tenemos la siguiente ecuación:  ∂N i ∂u~  ∂N i ∂u~ (48) k + k i = 1, 2,  , n ∫Ω  x ∂x ∂x y ∂y ∂y − N i Q  dΩ − ∫ΓN N i q dΓ = 0 Si ahora reemplazamos la aproximación u~ dada por la ec.(36), resulta un sistema de n ecuaciones algebraicas, cuyas incógnitas son las coordenadas generalizadas ai: k i1a1 + k i 2 a 2 +  + k in an − f i = 0

i = 1, 2,  , n

(49)

siendo los coeficientes kij definidos como  ∂N i ∂N j ∂N i ∂N j  k ij = ∫  k x + ky dΩ Ω ∂y ∂y   ∂x ∂x

(50)

además, los términos fi valen: f i = ∫ N i Q dΩ + ∫ N i q dΓ Ω

ΓN

En forma matricial este sistema de ecuaciones se puede escribir como

(51)

Introducción al Método de Elementos Finitos

 k11 k  21     k n1

k12 k 22  kn2

k1n   a1   f1   k 2 n  a 2   f 2   =            k nn  an   f n  

9

(52)

y en forma abreviada usando notación matricial Ka = f

(53)

La solución de este sistema de ecuaciones nos da los valores de las constantes ai. Observemos, que para el caso particular de la ecuación de Poisson se cumple k ij = k ji

(54)

Por lo tanto, la matriz K resulta simétrica. Esto no sucede con todas las ecuaciones diferenciales, sino solo con aquellas que cumplen con la propiedad de ser autoadjuntas[1]. Esta propiedad está relacionada con la simetría de la forma variacional respecto de las funciones de ponderación y la función solución. La ecuación de Poisson es autoadjunta. 2 FUNCIONES DE FORMA Describiremos los conceptos de funciones de forma y su continuidad para elementos finitos. Hasta ahora expresamos a las soluciones aproximadas como: n

u~ = ∑ ai N i ( x, y )

(55)

i =1

donde hemos asumido implícitamente que las funciones de prueba Ni estaban definidas por una expresión simple válida en todo el dominio Ω. Esto puede ser posible en el caso de dominios con geometría sencilla como rectángulos, círculos, elipses, pero no en el caso de geometrías más complejas. Una forma alternativa de definir las funciones de prueba consiste en subdividir el dominio Ω en una serie de subdominios ó elementos Ωe que no se superpongan, y luego las aproximaciones u~ se construyen por trozos usando definiciones simples de las funciones de prueba sobre estos subdominios. Si los subdominios son de forma relativamente simple y la definición de las funciones de prueba sobre estos subdominios puede ser hecha de manera repetitiva, es posible aproximar dominios complejos de forma bastante directa. Esta es la idea básica del método de elementos finitos el cual puede interpretarse como un método de aproximación donde las funciones de prueba se definen en forma local en cada elemento y son llamadas funciones de forma, estas funciones de forma se combinan para dar lugar a una aproximación por trozos. Considere, por ejemplo, un dominio unidimensional, esto es una recta de longitud L, sobre la cual queremos aproximar una función arbitraria u(x) mediante una aproximación lineal por trozos. Para ello subdividimos esta recta en m elementos de recta vinculados por sus extremos Nótese que hemos asociado los valores de las coordenadas generalizadas ai a los valores de la aproximación en los extremos de los elementos. Estos puntos de cada elemento que tienen asociados valores de las coordenadas generalizadas son llamados nodos del elemento.

10

Formulación de Elementos Finitos

u

a2

a4

a1 a3 e1

e2

x

e3

0

L Figura 2. Aproximación de u(x) con elementos lineales.

Si observamos la definición (55) de las funciones de aproximación es evidente que en un nodo i asociado a la coordenada generalizada ai la función de prueba Ni debe valer 1 y el resto de las funciones de prueba Nj ( j ≠ i ) deben ser nulas en ese punto. Luego, una vez identificados los puntos nodales, es muy sencillo definir las funciones de prueba Ni ya que deben valer 1 en el punto nodal asociado y 0 en el resto de los puntos nodales.

e1

e2

e3

N1 1 N2 1 N3 1 N4 1

Figura 3. Funciones de forma globales Ni lineales por trozos.

También puede observarse que las únicas funciones de forma globales que son diferentes de cero en cada elemento son aquellas asociadas con los nodos de ese elemento. Luego, en cada elemento e con nodos i,j la aproximación puede ser expresada simplemente en función de dos funciones de forma lineales del elemento, Nei, Nej y de los valores nodales ai, aj como

Introducción al Método de Elementos Finitos

11

u~ e = ai N ie + a j N ej

(56) Nej

Nei 1

1

j

i e

Figura 4. Funciones de forma locales N i lineales en cada elemento.

La aproximación global se genera a partir de la combinación de las funciones de forma locales en cada elemento. Por otro lado, observemos que si tenemos dos valores nodales por elemento podemos reproducir cualquier variación lineal sobre este elemento. En particular, si tenemos un valor constante de la aproximación u~ sobre el elemento e esto implica que los valores nodales deben ser iguales a este valor, esto es: (57) u~ e = c = a N e + a N e = c (N e + N e ) = cte i

i

j

j

i

j

Luego, sobre cada elemento se debe verificar

(N

e i

+ N ej ) = 1

(58)

Esto es, la suma de las funciones de forma de cada elemento debe ser igual a uno. La extensión de estos conceptos a dos dimensiones es bastante directa. En este caso, la subdivisión del dominio se efectúa utilizando triángulos ó cuadriláteros.

Figura 5. Aproximación en dos dimensiones.

12

Formulación de Elementos Finitos

En este caso los puntos nodales quedan asociados, en general, a los vértices de la malla y para el caso de aproximación lineal sobre cada triángulo las funciones de forma del elemento son planos. 2.1 Orden de continuidad. Nótese que en el ejemplo anterior la función de aproximación u~ era una función continua pero con discontinuidad en la derivada primera. Luego estas funciones se dicen que poseen continuidad de tipo C0. En general, si una función y sus m primeras derivadas son continuas, se dice que la función posee continuidad de tipo Cm. Consideremos una función f con apenas continuidad C0 formada por la unión de dos funciones como se muestra en la figura 6. Podemos imaginar que en el límite la derivada primera tiene una variación grande en torno del punto de unión pero se mantiene continua e integrable. Pero la derivada segunda tiende a infinito y no podemos asegurar que la derivada segunda sea integrable alrededor del punto de unión.

f

x

∂f ∂x

x

∂2 f ∂x 2

x

∞

Figura 6. Diferenciación de una función con discontinuidad de la derivada primera.

Luego se puede demostrar que si la máxima derivada de una función u que aparece en las integrales provenientes de métodos de residuos ponderados ó principios variacionales es de orden s entonces una condición suficiente para que estas integrales sean bien definidas, es que

Introducción al Método de Elementos Finitos

13

la función u posea como mínimo continuidad de tipo Cs-1. Luego, los elementos finitos cuyas funciones de forma posean este orden mínimo de continuidad se denominan compatibles. Notemos, que tanto para el caso de análisis de tensiones como para los problemas gobernados por la ecuación de Poisson, las máximas derivadas involucradas en las integrales son las derivadas primeras. Luego las funciones de prueba Ni solo deben poseer como mínimo continuidad de tipo C0. 2.2 Polinomios de Lagrange en una y dos dimensiones. Analizaremos interpolaciones del tipo polinómicas sobre elementos finitos en una dimensión y con continuidad C0.

Lineal

1

0

Cuadrática

0

2

1

Cúbica 0

1

2

3

Figura 7. Interpolación polinómica en un elemento unidimensional.

En estos elementos la continuidad de tipo C0 se consigue sencillamente ubicando puntos nodales en los extremos del elemento. Además con p+1 valores es posible definir un polinomio completo de grado p, luego si consideramos un elemento con nodos numerados de 0 a p distribuidos sobre el elemento con coordenadas x0, x1, x2, ... , xp-1, xp entonces las funciones de forma del elemento Nei serán polinomios de grado p que deben valer 1 en el nodo i asociado y cero en todos los otros nodos. Luego podríamos derivar una expresión para tales funciones de forma como N ie = α 0 + α 1 x + α 2 x 2 +  + α p x p

(59)

donde para encontrar las constantes αi imponemos las condiciones N ie ( xi ) = 1 N ej ( x j ) = 0

(60) j≠i

Sin embargo, este proceso es innecesario ya que por inspección el siguiente polinomio cumple con estas propiedades

14

Formulación de Elementos Finitos

N ie =

(x − x0 )(x − x1 ) (x − xi −1 )(x − xi +1 ) (x − x p ) (xi − x0 )(xi − x1 ) (xi − xi −1 )(xi − xi +1 ) (xi − x p )

(61)

Este polinomio se llama polinomio de interpolación de Lagrange de grado p. Si los nodos se encuentran equiespaciados entonces es posible expresar las funciones de forma en función de una coordenada local ξ del elemento definida como

ξ=

2(x − xce ) he

(62)

Donde xec es la coordenada del centro del elemento, he es la longitud del elemento. Luego la coordenada ξ está definida en el elemento en el rango –1 ≤ ξ ≤ 1 y las funciones de forma son para interpolación lineal: 1−ξ 2 1+ ξ e N1 = 2

(63)

N 0e =

para interpolación cuadrática las funciones de forma son

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

(64)

N 0e =

y para interpolación cúbica las funciones de forma son N 0e N1e N 2e N 3e

= − 169 (ξ + 13 )(ξ − 13 )(ξ − 1) 27 = 16 (ξ + 1)(ξ − 13 )(ξ − 1) 27 = − 16 (ξ + 1)(ξ + 13 )(ξ − 1) = 169 (ξ + 1)(ξ + 13 )(ξ − 13 )

(65)

La extensión a dos dimensiones es bastante directa para elementos finitos rectangulares si disponemos los nodos en una grilla rectangular sobre el elemento. y s=3 s=2

s=1 x

s=0 r=0

r=1

r=2

r=3

Figura 8. Elemento rectangular con interpolación Lagrangiana.

Introducción al Método de Elementos Finitos

15

Luego si rotulamos cada nodo con dos índices r,s (r = 0,1,2; s = 0,1,2) para el elemento cuadrático, entonces podemos escribir las funciones de forma para el nodo (r,s) como: N rse (x, y ) = N rp (x ) N sp ( y )

(66)

donde Npr es el polinomio de Lagrange de grado p asociado al nodo r en la dirección x. Notemos que el elemento rectangular es de aplicación limitada, debido a que no es posible aproximar un dominio complicado solo con rectángulos, sin embargo como veremos más adelante es posible extender esta formulación a elementos de forma cuadrilateral que si son ampliamente utilizados. 2.3 Polinomios de Hermite en una y dos dimensiones. En algunos problemas, como por ejemplo el análisis de la flexión de placas y láminas, se requieren funciones con continuidad C1 para interpolar los desplazamientos transversales. En una dimensión las funciones tipo “viga” que tienen como variables nodales a los desplazamientos y rotaciones en los extremos del elemento tienen esta propiedad. θ1

θ2

w2

w1

Figura 9. Variables nodales para elemento tipo “viga”.

Luego los desplazamientos transversales se pueden interpolar como: (67)

2

w = ∑ H wi wi + H θi θ i i =1

donde Hwi, Hθi son las funciones de forma asociadas a los desplazamientos y rotaciones nodales. Estas funciones son polinomios cúbicos conocidos como polinomios de Hermite y se muestran en la figura 10.

1

H w1

1

H w2

H θ2

1

H θ1

1

Figura 10. Funciones de forma con continuidad C1 para elementos tipo viga.

16

Formulación de Elementos Finitos

Estas funciones de forma se pueden expresar como N w1 = 1 − 3ξ 2 + 2ξ 3

(68)

N w 2 = 3ξ 2 − 2ξ 3

N θ1 = L (ξ − 2ξ 2 + ξ 3 ) N θ 2 = L (− ξ 2 + ξ 3 )

donde ξ = x/L. Uno de los primeros elementos finitos para placas fue un elemento rectangular que poseía como variables nodales a los desplazamientos y rotaciones en los vértices. θ x4 z

x

θ x3 θ y4

w4

θ x1

w3

θ y3

θ x2

y

θ y1 w1

θ y2 w2

Figura 11. Variables nodales para elemento de placa rectangular.

Obsérvese que hemos representado las rotaciones mediante vectores y el sentido positivo adoptado para las rotaciones no coincide con el de los ejes. Es importante verificar cuando se emplea un programa de elementos finitos de placas el sentido positivo adoptado para las rotaciones. Luego, para este elemento los desplazamientos transversales se aproximan como 4

w = ∑ N wi wi + N θxi θxi + N θyi θyi

(69)

i =1

Las funciones de forma de este elemento se obtienen mediante el producto de funciones de tipo viga. Así, para el nodo 1 las funciones de forma son N w1 = H w1 ( x ) H w1 ( y ) N θx1 = H θx1 ( x ) H w1 ( y ) N θy1 = H w1 ( x ) H θy1 ( y )

(70)

La desventajas de estos elementos es que solamente poseen continuidad C1 cuando el elemento tiene forma rectangular pero si se lo distorsiona, usando una formulación isoparamétrica, como veremos más adelante, para transformarlo en un cuadrilátero que posee mayor habilidad para representar geometrías complejas se pierde la continuidad C1 y queda solo con continuidad C0. Por ello este elemento si bien tiene gran precisión es poco utilizado por las restricciones que le impone la forma rectangular.

Introducción al Método de Elementos Finitos

1

1

Nw1 = H w1 (x) H w1 (y)

17

Nθ x1 = H θx1 (x) H w1 (y)

Figura 12. Funciones de forma con continuidad C1 para elementos rectangulares.

2.1.4

Convergencia de la aproximación. Test de la parcela.

Usando elementos finitos obtenemos soluciones aproximadas, idealmente si incrementamos el número de elementos empleados en un análisis, la solución debería converger a la solución exacta del problema. Hasta ahora solo se ha requerido que de las funciones de aproximación sean integrables, pero esto no es suficiente para garantizar la convergencia. Notemos que los integrandos de las integrales de los métodos de residuos ponderados y de las integrales que provienen de principios variacionales deben ser constantes en el límite cuando el tamaño de los elementos tiende a cero. Esto es, a medida que el tamaño de los elementos disminuye estos integrandos deben tender a un valor constante sobre cada elemento. Luego una condición necesaria para la convergencia de las funciones de aproximación es que para un tamaño finito los elementos finitos deben ser capaces de reproducir cualquier valor constante de los integrandos que pueda ocurrir en el limite del tamaño de los elementos tendiendo a cero. Así, por ejemplo, para un análisis de tensiones, la condición de convergencia implica que un elemento finito debe ser capaz de reproducir cualquier estado de deformación constante. Ya que los integrandos, en este caso, son funciones de las derivadas primeras de los desplazamientos, la condición de convergencia implica que las funciones de aproximación deben ser capaces de reproducir cualquier variación lineal de desplazamientos sobre el elemento. La construcción de funciones de prueba con continuidad mayor que C0 puede ser bastante complicada, en particular para dos dimensiones. Esto ha hecho que en algunos casos se utilicen elementos con menor continuidad que la mínima requerida, a estos elementos se los denomina incompatibles. La utilización de elementos incompatibles requiere que estos elementos cumplan requisitos adicionales para poder ser utilizados con seguridad. Uno de estos requisitos adicionales es el llamado test de la parcela[2] (patch test), que consiste en que el elemento debe poder reproducir a escala finita cualquier comportamiento básico infinitesimal para asegurar la convergencia. En general, este test consiste en testear con la computadora estados básicos a escala finita, y no sirven solamente para verificar la convergencia sino también el correcto funcionamiento del programa.

18 3

Formulación de Elementos Finitos IMPLEMENTACIÓN MATRICIAL DEL MÉTODO DE ELEMENTOS FINITOS.

Describiremos los procedimientos matriciales usados en el método de elementos finitos aplicado a problemas de análisis de tensiones. Estos métodos incluyen el ensamblaje de elementos, la imposición de condiciones de contorno, la solución del sistema de ecuaciones para obtener las cantidades nodales y el procesamiento de elementos para obtener cantidades tales como las tensiones. 3.1 Expresión matricial de la energía potencial total. Consideremos un cuerpo plano que puede representarse mediante un dominio bidimensional Ω discretizado mediante elementos finitos. La energía potencial total V de un cuerpo elástico lineal viene dada por la suma de la energía potencial de deformación U y de la energía potencial Ve asociada al trabajo de las fuerzas externas. V = U + Ve

(71)

La energía potencial de deformación U se puede expresar como U=

1 2

∫

Ω

ε T D 0 h dΩ

(72)

donde ε es el vector de deformaciones, D es la matriz constitutiva y h es el espesor. εx    0 = ε y  γ   xy 

,

1 ν E  ν 1 D= 1 −ν 2   0 0

0  0  (1 − ν )  2 

(73)

La energía potencial Ve asociada al trabajo de las fuerzas externas es Ve = − ∫ b T u h dΩ − ∫ t T u h dΓ Ω

(74)

Γ

donde b es el vector de fuerzas de volumen, t es el vector de fuerzas de superficie y u es el vector de desplazamientos bx  b=  b y 

,

t x  t=  t y 

,

u  u=  v 

(75)

Luego la energía potencial total se puede expresar como V = U + Ve =

1 2

∫

Ω

ε T D 0 h dΩ − ∫ b T u h dΩ − ∫ t T u h dΓ Ω

Γ

(76)

Si usamos una aproximación por elementos finitos es necesario dividir el dominio Ω en elementos y podemos expresar la energía potencial total como nel

V = U + Ve = ∑ 12 ∫ ε T D 0 h dΩ − ∫ u T b h dΩ − ∫ u T t h dΓ e =1

Ωe

Ωe

Γe

(77)

siendo nel el número de elementos, Ωe la región ocupada por cada elemento y Γe su contorno cargado

Introducción al Método de Elementos Finitos

19

Para poder obtener una aproximación por elementos finitos debemos aplicar el método de Rayleigh-Ritz utilizando los campos de desplazamientos u formados por las funciones de forma de los elementos finitos. 3.2 Aproximación por elementos finitos. El primer paso para obtener una aproximación por elementos finitos es realizar una discretización del dominio. Esto es, debemos generar una malla de elementos finitos que cubra todo el dominio. Además debemos numerar los nodos de la malla, que son aquellos puntos que tienen asociadas coordenadas generalizadas. Para el caso particular de análisis de tensiones las coordenadas generalizadas son los desplazamientos nodales. Así para el nodo i sus desplazamientos nodales serán ui y vi. Consideremos una aproximación por elementos finitos de los desplazamientos u, como 3

2

⇒

Ω

4

7 11

5

8 9 10

1 13

12

Figura 13. Discretización de un dominio Ω. n

n

u = ∑ u i N i ( x , y ),

v = ∑ vi N i ( x, y )

i =1

(78)

i =1

donde ui, vi son los desplazamientos nodales. Cada función de prueba Ni se compone de las funciones de forma asociadas al nodo i de todos los elementos que contienen ese nodo, esto es nel

N i ( x, y ) = ∑ N ie ( x, y )

(79)

e =1

donde nel es el número de elementos de la malla. Además, los desplazamientos ue, ve en cada elemento se pueden expresar como ue =

nnod

∑ u ej N ej ( x, y ),

ve =

j =1

nnod

∑v j =1

e j

N ej ( x, y )

(80)

donde nnod es el número de nodos del elemento. Nótese que en este caso el índice j se refiere a la numeración local del nodo en el elemento. Así, por ejemplo, para un triángulo de 3 nodos los desplazamientos sobre el elemento son u e = N 1e u1e + N 2e u 2e + N 3e u 3e v =N v +N v +N v e

e 1

e 1

e 2

e 2

esta ecuación puede escribirse matricialmente como

e 3

e 3

(81)

20

Formulación de Elementos Finitos ue = Ne de

(82)

v3e u3e

3

v1e

u1e

1

v2e u2e 2

Figura 14. Desplazamientos nodales para el triángulo de 3 nodos.

siendo Ne la matriz de funciones de forma del elemento  N1e N = 0

0 N1e

e

N 2e 0

N 3e

0 N 2e

0

0   N 3e 

(83)

y de es el vector de desplazamientos nodales del elemento. d e = {u1e T

v1e

u 2e

v2e

u3e

v3e }

(84)

En forma particionada la matriz de funciones de forma se puede escribir como

[

N e = N1e

N e2

N 3e

]

(85)

donde las submatrices Nie que están asociadas a cada nodo del elemento son 

e i

N e = i 0

0   N ie 

(86)

El vector de de desplazamientos nodales del elemento se puede también expresar en forma particionada como d1e    d e = d e2  d e   3

(87)

donde los vectores die que están asociados a los desplazamientos de cada nodo i del elemento son

Introducción al Método de Elementos Finitos

21

u e  d ie =  ie   vi 

(88)

siendo uie, vie los desplazamientos del nodo i del elemento. 3.3 Matriz gradiente. Si reemplazamos los campos de desplazamientos aproximados por elementos finitos en las expresiones de las deformaciones, en cada elemento tenemos ∂u e ∂N 1e e ∂N 2e e ∂N 3e e u1 + u2 + u3 = ∂x ∂x ∂x ∂x ∂v e ∂N 1e e ∂N 2e e ∂N 3e e εy = v1 + v2 + v3 = ∂y ∂y ∂y ∂y

εx =

γ xy =

(89)

∂u e ∂v e ∂N 1e e ∂N 1e e ∂N 2e e ∂N 2e e ∂N 3e e ∂N 3e e u1 + v1 + u2 + v2 + u3 + v3 + = ∂y ∂x ∂y ∂x ∂y ∂x ∂y ∂x

y en forma matricial  ∂u e   ∂N 1e    ∂x    ∂x  ∂v e   ε= = 0  e ∂y e   e ∂v   ∂N 1  ∂u  ∂ y + ∂x   ∂y   

∂N 2e ∂x

0 ∂N 1e ∂y ∂N 1e ∂x

∂N 3e ∂x

0 ∂N 2e ∂y ∂N 2e ∂x

0 ∂N 2e ∂y

0 ∂N 3e ∂y

e  u1  0  v e    1e  ∂N 3e  u 2    ∂y   v2e   ∂N 3e  u 3e  ∂x   v3e 

(90)

y en forma abreviada 0

= Be de

(91)

donde Be es la matriz gradiente del elemento y de es el vector de desplazamientos nodales del elemento. En forma particionada la matriz gradiente se puede escribir como

[

B e = B1e

B e2

B 3e

]

(92)

donde las submatrices Bie que están asociadas a cada nodo del elemento son  ∂N ie   ∂x e Βi =  0   e  ∂N i  ∂y

 0   ∂N ie  ∂y   ∂N ie  ∂x 

(93)

Observemos que un caso general la matriz gradiente del elemento Be estará compuesta de tantas submatrices Bie como nodos tenga el elemento.

22

Formulación de Elementos Finitos

3.4 Matriz de rigidez y vector de cargas nodales equivalentes. Si reemplazamos los campos de desplazamientos aproximados por elementos finitos en la expresión (2.31) de la energía potencial total tenemos nel

V = U + Ve = ∑ 12 ∫ d e B e D B e d e h dΩ − ∫ d e N e b h dΩ − ∫ d e N e t h dΓ e =1

T

T

T

Ωe

T

Ωe

T

T

(94)

Γe

Definiendo a la matriz de rigidez del elemento como e

= ∫ B e D B e h dΩ T

(95)

Ωe

esta matriz es una matriz cuadrada de dimensión igual a la cantidad de desplazamientos nodales del elemento y definiendo además al vector de cargas nodales equivalentes del elemento como f e = ∫ N e b h dΩ + ∫ N e t h dΓ T

T

Ωe

Γe

(96)

luego la energía potencial total se puede expresar como nel

V = ∑ 12 d e K e d e − d e f e T

T

(97)

e =1

Si empleamos la forma particionada (2.46) para las matrices gradiente del elemento Be entonces la matriz de rigidez del elemento Ke se puede expresar en forma particionada como e  K 11  K e = K e21 e  K 31 

e K 12

K e22 K e32

e  K 13 e  K 23  e  K 33 

(98)

siendo e ij

= ∫ B ie D B ej h dΩ T

Ωe

(99)

la submatriz de rigidez del elemento que relaciona los nodos numerados localmente como i, j en el elemento. Si definimos al vector d de desplazamientos de la malla con n nodos, como  d1  d   2 d=   d n 

(100)

entonces la energía potencial de deformación se puede expresar como nel

U = ∑ 12 d e K e d e = 12 d T K d T

e =1

siendo K la matriz de rigidez global formada por las submatrices Kij valen

(101)

Introducción al Método de Elementos Finitos

23

nel

K ij = ∑ K ije

(102)

e =1

Esto es, si dos nudos están vinculados por un elemento, entonces dicho elemento debe contribuir con una submatriz a la matriz de rigidez global. Por otro lado, la energía potencial de las fuerzas externas se puede expresar como: Ve = −d T f

(103)

siendo f el vector de fuerzas externas global cuyas componentes fi valen nel

fi = ∑ fi e

(104)

e =1

Finalmente, la energía potencial total queda V = 12 d T K d − d T f

(105)

Aplicando Rayleigh-Ritz debemos minimizar esta expresión respecto de las coordenadas generalizadas, que en este caso son los desplazamientos nodales d, esto es ∂V = Kd −f = 0 ∂d

(106)

Resultando el siguiente sistema de ecuaciones Kd =f

(107)

que tiene por incógnitas a los desplazamientos nodales de toda la malla. En general, algunos de estos desplazamientos tendrán valores prescritos por lo que no serán incógnitas, en este caso deberíamos eliminar la línea correspondiente a este desplazamiento de la matriz de rigidez global. Obsérvese que el primer paso para resolver este sistema de ecuaciones es el montaje de la matriz K y del vector f a partir de las contribuciones de los elementos, este proceso se denomina ensamblaje. 3.5 Ejemplo de ensamblaje. Considere una malla formada únicamente por dos elementos y cuatro nodos 4

3 2 1

1

4 Figura 15. Ensamblaje de dos elementos.

24

Formulación de Elementos Finitos

Notemos que los nodos 1 y 3 son compartidos por ambos elementos, luego habrá aportes de ambos elementos a las posiciones K11, K13 y K31 de la matriz global de rigidez. Luego la matriz de rigidez global ensamblada queda (1) ( 2) ( 2) (K 11 + K 11 ) K 12  K (2) K (222 ) K =  (1) 21 ( 2 ) ( 2) (K 31 + K 31 ) K 32  K (412 ) 0 

(1) ( 2) (2)  (K 13 + K 13 ) K 14  K (232 ) 0  (1) (2)  (K 33 + K (332 ) ) K 34  K (432 ) K (442 ) 

(108)

donde hemos indicado los aportes de cada elemento por el índice superior entre paréntesis. Notando que el único lado cargado está sobre el elemento 1, entonces el vector de fuerzas nodales queda:  0   0    f =  (1)  f 3  f 4(1) 

(109)

Para poder resolver el sistema de ecuaciones debemos eliminar previamente los desplazamientos prescriptos. Esto es debemos eliminar la fila y columna de la matriz asociadas a un desplazamiento de valor conocido. Notemos que los nodos deben ser compartidos por todos los elementos adjacentes a ese nodo. Además, no deben mezclarse elementos con diferentes números de nodos por lado. Esto es, no son permitidas situaciones como las mostradas en la figura 16.

cuadrático lineal

nodos desconectados Figura 16. Ensamblaje no permitido de elementos.

Sin embargo si está permitido ensamblar triángulos con cuadriláteros, siempre que posean el mismo número de nodos por lado.

cuadrático cuadrático

Figura 17. Ensamblaje válido de elementos.

Introducción al Método de Elementos Finitos 4

25

ELEMENTO FINITO TRIANGULAR.

Consideremos un triángulo de tres nodos. Los desplazamientos de un punto cualquiera del elemento se pueden expresar en función de los desplazamientos nodales como: u e = N 1e u1e + N 2e u 2e + N 3e u 3e v =N v +N v +N v e

e 1

e 1

e 2

e 2

e 3

(110)

e 3

Los tres nodos del elemento definen una variación lineal del campo de desplazamientos que puede escribirse como u e = α1 + α 2 x + α 3 y

(111)

v = α 4 + α5 x + α6 y e

Estos campos de desplazamientos deben coincidir en los nodos con las correspondientes incógnitas nodales u1e = α1 + α 2 x1 + α 3 y1

(112)

u = α 1 + α 2 x2 + α 3 y 2 e 2

u3e = α1 + α 2 x3 + α 3 y3 donde xi, yi son las coordenadas del nodo i. Luego tenemos tres ecuaciones con tres incógnitas, α1, α2, α3. Resolviendo este sistema de ecuaciones se obtiene la siguiente expresión para u u1e =

[

]

(113)

i , j , k = 1,2,3

(114)

1 (a1 + b1 x + c1 y )u1e + (a 2 + b2 x + c2 y )u 2e + (a3 + b3 x + c3 y )u3e e 2A

donde ai = x j y k − x k y j ,

bi = y j − y k ,

ci = x k − x j ,

y Ae es el área del elemento definida como Ae =

1 (c3b2 − c2b3 ) 2

(115)

Comparando la expresión (110) con la (113) se deduce que las funciones de forma son N ie =

1 (ai + bi x + ci y ) , 2 Ae

i = 1,2,3

(116)

Luego las matriz gradiente del nodo i es  ∂N ie   ∂x Β ie =  0  e  ∂N i  ∂y

 0   ∂N ie  1 =  ∂y 2 Ae  ∂N ie  ∂x 

(117) bi 0  c i

0 c i  bi 

Con la matriz gradiente ya podemos calcular la matriz de rigidez del elemento.

26

Formulación de Elementos Finitos

4.1 Coordenadas de área. Si se une un punto interior P de un triángulo de área A con los tres vértices se obtienen tres subáreas A1, A2, A3. Las coordenadas de área ξi se definen como los cocientes ξ1 =

A1 A

ξ2 =

A2 A

ξ3 =

A3 A

(118)

Como el área total A = A1 + A2 + A3, esto implica que las coordenadas de área no son independientes ya que deben cumplir la restricción ξ1 + ξ 2 + ξ 3 = 1

(119)

3

y

Lado 1 A2 Lado 2

2

A1 P A3 Lado 3

1 x Figura 18. Coordenadas de área para triángulos.

Cada coordenada de área ξi varía linealmente sobre el triángulo, ya que cada área Ai es proporcional a la distancia al lado i. Luego existe una relación lineal entre las coordenadas cartesianas x, y de un punto y las coordenadas de área que juntamente con la condición de restricción (119) se pueden escribir como x = ξ1 x1 + ξ 2 x2 + ξ 3 x3 y = ξ1 y1 + ξ 2 y 2 + ξ 3 y3

(120)

1 = ξ1 + ξ 2 + ξ 3 Tenemos un sistema de tres ecuaciones con tres incógnitas ξ1, ξ2, ξ3 resolviendo resultan ξi =

1 (ai + bi x + ci y ) 2 Ae

(121)

Nótese que las coordenadas de área coinciden con las funciones de forma lineales (116). La ventaja de usar coordenadas de área para definir las funciones de forma reside en la facilidad para integrar funciones polinómicas en estas coordenadas, ya que es posible utilizar la siguiente fórmula

Introducción al Método de Elementos Finitos

∫

e

ξ1k ξ l2 ξ 3m dA = 2 A e

27

k!l! m! (2 + k + l + m )

(122)

donde k, l, m son enteros no negativos. La formulación de matrices de elementos requiere que las funciones de forma Ni(ξ1, ξ2, ξ3) sean derivadas respecto de la coordenadas cartesianas, aplicando la regla de la cadena tenemos ∂N i ∂N i ∂ξ1 ∂ N i ∂ ξ 2 ∂N i ∂ ξ 3 = + + ∂x ∂ξ1 ∂x ∂ξ 2 ∂x ∂ξ 3 ∂x

(123)

∂N i ∂N i ∂ξ1 ∂ N i ∂ ξ 2 ∂N i ∂ ξ 3 = + + ∂y ∂ξ1 ∂y ∂ξ 2 ∂ y ∂ξ 3 ∂ y Siendo b ∂ξ i = ie ∂x 2 A ∂ξ i c = ie ∂y 2 A

(124)

4.2 Campos de interpolación para triángulos. Los elementos triangulares permiten que se usen polinomios completos para definir las funciones de forma. Los coeficientes polinomiales que intervienen en un polinomio completo vienen dados por el triángulo de Pascal. Así, un polinomio de grado p debe incluir todos los coeficientes de las primeras (p+1) filas. El número n de términos de un polinomio completo de grado p es n = (p+1)( p+2)/2. Triángulo de Pascal

Grado del polinomio

Número de términos n = (p+1) (p+2) / 2

1

0 1 2 3 4

1 3 6 10 15

x

y

x2 xy y2 x3 x2y xy2 y3 4 x x3y x2y2 xy3 y4

Una forma de generar funciones de forma consiste en utilizar un polinomio completo en coordenadas de área cuyas constantes deben ser determinadas. Así, por ejemplo, para un polinomio completo de segundo grado las funciones de forma serian N i (ξ1 , ξ 2 ) = α 0 + α1 ξ1 + α 2 ξ 2 + α 3 ξ1ξ 2 + α 4 ξ12 + α 5 ξ 22

(125)

Luego se definen tanto puntos nodales como constantes haya para determinar y se impone la condición de que cada función de forma valga uno en un nodo y cero en el resto, lo que genera suficientes ecuaciones para hallar las constantes. La ubicación de los puntos nodales no puede ser totalmente arbitraria, ya que en general las funciones de forma deben ser continuas entre elementos. Esto implica que debe ser impuesta la misma variación polinómica sobre los lados de los elementos. Esto es, sabemos que un polinomio de grado p en una variable tiene p + 1 coeficientes, luego si sobre el triángulo se

28

Formulación de Elementos Finitos

define una función de forma polinómica completa de grado p entonces debemos tener exactamente p + 1 nodos sobre cada lado. En la figura 19 tenemos los puntos nodales correspondientes a funciones de forma polinómicas completas de primer, segundo y tercer grado.

3

3

3

7

5

6 8

2

9

4

4

2 5

1

1

1

10

2

6

Figura 19. Elemento triangular lineal, cuadrático y cúbico.

Las funciones de forma para el triángulo lineal son N1 = ξ1

N 2 = ξ2

N 3 = ξ3

(126)

Las funciones de forma para el triángulo cuadrático son N1 = ξ1 (2ξ1 − 1)

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

N 4 = 4ξ1ξ 2

N 5 = 4ξ 2 ξ 3

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

(127)

Las funciones de forma para el triángulo cúbico son N1 = 12 ξ1 (3ξ1 − 1)(3ξ1 − 2 )

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

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

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

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

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

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

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

(128)

9 2

N 9 = 92 ξ 3 ξ1 (3ξ1 − 1)

N10 = 27 ξ1ξ 2 ξ 3

5

FORMULACIÓN ISOPARAMÉTRICA

Presentaremos una descripción de la formulación isoparamétrica, la cuál es utilizada para generar muchos tipos útiles de elementos. Actualmente, la gran mayoría de los programas comerciales de elementos finitos utiliza elementos basados en esta formulación debido a su mayor precisión. Esta formulación está basada en realizar un cambio de coordenadas debido a lo cual es necesario replantear las integrales sobre los elementos, las cuales se tornan complicadas de evaluar analíticamente, por lo que se debe recurrir a integrarlas numéricamente en forma aproximada. Por ello se presentará también una descripción de los métodos de integración numérica utilizados. Los elementos rectangulares son más eficientes que los elementos triangulares para la misma cantidad de grados de libertad, sin embargo tienen escasa capacidad para modelar geometrías complejas. La formulación isoparamétrica permite generar elementos con lados

Introducción al Método de Elementos Finitos

29

curvos ó no rectangulares lo que tiene obvias ventajas para el modelado de geometrías de formas arbitrarias y bordes curvos. Para formular elementos isoparamétricos deben usarse sistemas de coordenadas naturales. Consideremos la generalización de un elemento cuadrado de cuatro nodos a una forma cuadrilateral arbitraria. η η =1

4

y

3

3 4

x

ξ 1

η =-1

1

2

ξ =-1

2

ξ =1

Figura 20. Mapeo del elemento cuadrilateral de cuatro nodos.

Una forma posible de definir la geometría del cuadrilátero consiste en emplear funciones de forma definidas sobre el cuadrado para interpolar las coordenadas de los vértices del cuadrilátero. Empleando un sistema de coordenadas paramétrico ξ, η sobre el cuadrado, donde cada una de estas coordenadas varía entre –1 y 1, y si xi, yi son las coordenadas cartesianas de cada vértice del cuadrilátero entonces es posible establecer la siguiente transformación de coordenadas 4

x = ∑ N i* (ξ ,η ) xi

(129)

i =1 4

y = ∑ N i* (ξ , η ) yi i =1

Donde Ni*(ξ,η) son las funciones de forma Lagrangianas, que expresadas usando las coordenadas paramétricas ξ,η valen N1* = 14 (1 − ξ )(1 − η )

(130)

N = (1 + ξ )(1 − η ) * 2

1 4

N 3* = 14 (1 + ξ )(1 + η ) N 4* = 14 (1 − ξ )(1 + η ) También es posible utilizar funciones cuadráticas para reproducir bordes curvos η η =1

y 7

4

7

4 6

8

η =-1

1

ξ =-1

5

3

3

2

ξ =1

ξ

6

x

8 1

5 2

Figura 21. Mapeo del elemento de lados cuadráticos de ocho nodos.

30

Formulación de Elementos Finitos

En este caso tenemos dos posibilidades de elementos cuadráticos, de 8 ó de 9 nodos, según si se considera ó no el nodo central. El elemento de 9 nodos es llamado Lagrangiano pues sus funciones de forma son productos de polinomios de Lagrange

½ η (1-η )

η

7

4 8

6

1 5

½ ξ(1-ξ)

8

ξ

9

η

7

4

(1-η 2)

3

3 6

5

½ ξ(1-ξ)

2

ξ

9

1 2

N 8 =½ ξ(1-ξ)(1-η 2 )

N 1 =¼ ξ(1-ξ)η (1-η ) (1-ξ2 )

η

7

4 8

3

ξ

9 6

1 5

(1-η 2)

2

N 9 = (1-ξ2)(1-η 2 )

Figura 22. Elemento Lagrangiano de nueve nodos.

El elemento de ocho nodos se obtiene combinando las funciones de forma del elemento de cuatro nodos con funciones de nodos intermedios modificados y es llamado Serendípito. η

7

4

(1-η 2)

3

8 6

1 5

½ (1-ξ)

η

7

4

½(1-η )

5

½ (1-ξ)

½(1-η )

ξ

6

1

N 8 =½(1-ξ)(1-η 2)

2

3

8

ξ

-½

2

-½

7

4

η

5

2

N 5 = ½(1-ξ2)(1-η )

η

3 6

1

ξ

7

4

=

η

N 8 =½(1-ξ)(1-η 2)

2

3

8 6

1 5

ξ

6

1

(1-ξ2)

8

½ (1-ξ)

3

8

N 1L =¼ ξ(1-ξ)η (1-η )

(1-η 2)

7

4

5

N 1 = N 1L -½ N 2-½ N 5

Figura 23. Elemento Serendípito de ocho nodos.

2

ξ

Introducción al Método de Elementos Finitos

31

En tres dimensiones también es posible generar elementos en coordenadas curvilíneas. Uno de los elementos más usados es el hexaedro trilineal ó elemento tipo ladrillo que posee ocho nodos y una variación lineal de los lados. ζ

8 5

8

5

1

7

6

x3

η 1 2

ξ

4

6

2

3

4

7 3 x2

x1

Figura 24. Mapeo del elemento hexaédrico de ocho nodos

Nótese que en todos los casos la variación de la geometría de cada lado solo depende de la posición de los nodos ubicados sobre ese lado. Esto permite asegurar que no habrá discontinuidad de la geometría entre elementos adyacentes. 5.1 VARIACIÓN DE FUNCIONES EN LOS ELEMENTOS CURVILÍNEOS Con la geometría del elemento definida mediante funciones de forma Ni* debemos considerar como especificar la variación de una función u sobre el elemento distorsionado. Una forma conveniente de definir esta variación es mediante el uso de funciones de forma también definidas en coordenadas paramétricas, de la manera usual nen

u = ∑ N i (ξ , η ) u i

(131)

i =1

Donde nen es el número de nodos del elemento donde se han especificado valores nodales ui. Estos valores nodales pueden estar asociados ó no a los mismos nodos empleados para especificar la geometría. Si se usan los mismos puntos para definir la geometría y la variación de la función u el elemento tenemos las mismas funciones de forma N i = N i*

(132)

en este caso el elemento es llamado isoparamétrico. Si se introducen más nodos para definir la variación de u que la geometría el elemento es llamado subparamétrico, y si por el contrario se utilizan más nodos para definir la geometría que la función u el elemento es llamado superparamétrico. 5.1.1

Evaluación de matrices del elemento.

Para calcular los coeficientes de las matrices de elemento debemos evaluar integrales sobre el dominio del elemento, por ejemplo, la matriz de rigidez de elemento para un problema de estados planos es e

= ∫ B e D B e h dx dy T

Ωe

donde Be es la matriz gradiente del elemento

(133)

32

Formulación de Elementos Finitos

[

B e = B1e

B e2

B 3e  B enen

]

(134)

y Bie es la submatriz gradiente del nodo i formada por derivadas de las funciones de forma respecto de las coordenadas cartesianas  ∂N i   ∂x e Βi =  0   ∂N  i  ∂y

 0   ∂N i  ∂y  ∂N i   ∂x 

(135)

Aplicando la regla de diferenciación tenemos ∂N i ∂N i = ∂ξ ∂x ∂N i ∂N i = ∂η ∂x

∂ x ∂N i + ∂ξ ∂y ∂x ∂N i + ∂η ∂ y

∂y ∂ξ ∂y ∂η

(136)

ó en forma matricial  ∂N i   ∂x  ∂ξ   ∂ξ  ∂N  =  ∂x  i   ∂η   ∂η

∂y   ∂N i   ∂N i   ∂x  ∂ξ   ∂x    ∂N  = J  ∂N  ∂y   i  i ∂η   ∂y   ∂y 

(137)

donde J es la matriz Jacobiana de la transformación cuyos coeficientes son ∂x nnge ∂N i* xi =∑ ∂ξ i =1 ∂ξ

∂x nnge ∂N i* xi =∑ ∂η i =1 ∂η

(138)

siendo nnge el número de nodos empleados para definir la geometría. Luego, es posible hallar la relación inversa  ∂N i   ∂N i     ∂x  −1  ∂ξ   ∂N  = J  ∂ N   i  i  ∂y   ∂η 

(139)

Donde J-1 es la inversa de la matriz Jacobiana

J

−1

 ∂y 1  ∂η =  J  − ∂x  ∂η

∂y  ∂ξ   ∂x  ∂ξ 

−

(140)

siendo J = J(ξ,η) el determinante de la matriz Jacobiana, usualmente llamado el Jacobiano

Introducción al Método de Elementos Finitos

J=

33

∂x ∂y ∂x ∂y − ∂ξ ∂η ∂η ∂ ξ

(141)

El Jacobiano también permite expresar el diferencial de área en coordenadas paramétricas dx dy = J dξ dη

(142)

Luego, la submatriz gradiente del nodo i se puede expresar como Gxi 1 e 1 Β = Gi =  0 J J Gy i e i

0  Gy i  Gxi 

(143)

donde Gxi, Gyi son polinomios en ξ,η que valen ∂N i ∂ξ ∂N i Gy i = ∂η Gx i =

∂y ∂ N i − ∂η ∂η ∂x ∂N i − ∂ξ ∂ξ

∂y ∂ξ ∂x ∂η

Luego definiendo a la matriz G del elemento como

[

G e = G 1e

G e2

G 3e  G enen

(144)

]

(145)

podemos escribir a la matriz de rigidez como e

T

Ge DGe h dξ dη =∫ ∫ −1 −1 J 1

1

(146)

Debe observarse que el integrando es un cociente de polinomios en ξ,η y la integración analítica de estas integrales se torna bastante complicada, por lo que en general se utiliza integración numérica. 5.2 INTEGRACIÓN NUMÉRICA Existen varias fórmulas para obtener el valor de una integral definida en forma aproximada. La mayoría de las expresiones son del tipo n

I = ∫ f (ξ ) dξ = ∑ wi f (ξ i ) 1

−1

(147)

i =1

Esto es, se transforma el cálculo de la integral en una suma ponderada por pesos wi de evaluaciones del integrando en n puntos ξi. Dado que los integrandos que debemos evaluar son en general complicados es conveniente utilizar un método que utilice el menor número de evaluaciones posibles. En particular, es ampliamente utilizado la cuadratura de Gauss-Legendre que requiere un mínimo de evaluaciones del integrando. 5.2.1

Cuadratura de Gauss-Legendre.

En este método la posición de los puntos ξi está especificada y tenemos diferentes pesos según la cantidad de puntos utilizadas. Si se utilizan n puntos de integración se puede integrar

34

Formulación de Elementos Finitos

exactamente un polinomio de grado (2n-1). En la tabla 3.1 se muestran las coordenadas de los puntos y sus respectivos pesos para diferentes número de puntos, para evaluar una integral del tipo de la ec.(147). ± ξi 0 0.577350269 0.774596669 0.000000000 0.861136312 0.339981044 0.906179846 0.538469310 0.000000000

n 1 2 3 4 5

wi 2 1 0.555555555 0.888888888 0.347854845 0.652145154 0.236926885 0.478628671 0.568888889

Tabla 3.1. Coordenadas y pesos para cuadratura de Gauss-Legendre

En dos dimensiones aplicamos la fórmula para cada dimensión, esto es 1

I =∫

1

∫

−1 −1

1  n  f (ξ , η ) dξ dη = ∫ dξ ∑ wi f (ξ , ηi )  = −1  i =1  n

(148)

n

= ∑ ∑ wi w j f (ξ i ,η j ) i =1 j =1

En forma análoga se puede extender la fórmula a tres dimensiones. ξ =-1

η

η

ξ =1

η

η =1 ξ

ξ

ξ

η =-1

2x2

1x1

3x3

Figura 25. Posición de puntos de Gauss para cuadrilátero.

Los puntos de Gauss tienen, en general, la propiedad de ser puntos óptimos para el cálculo de tensiones. Esto es, las tensiones en estos puntos suelen ser más precisas y más cercanas a los valores exactos. 5.2.2

Integración completa y reducida.

Una cuestión importante es determinar el mínimo número de puntos de integración para las integrales de elementos finitos. Estas integrales se pueden expresar como I =∫

Ωe

donde Ωe es el volumen de cada elemento.

f dΩ e

(149)

Introducción al Método de Elementos Finitos

35

Se desea obtener el mínimo número de puntos de manera de asegurar la convergencia, esto es, que cuando el tamaño del elemento tienda a cero se obtenga la solución exacta. Sabemos que cuando el tamaño del elemento tiende a cero el integrando f tiende a ser constante, ya que depende de derivadas que en el límite tienden a ser constantes. Luego debe emplearse un regla de integración que sea capaz de integrar exactamente el volumen del elemento, esto es Ω e = ∫ e dΩ e = ∫

1

Ω

∫

1

−1 −1

J dξ dη

(150)

Notese que haciendo el cambio de coordenadas aparece el Jacobiano en el integrando, esto implica que se debe utilizar un número mínimo de puntos de integración de manera de integrar exactamente al Jacobiano en las coordenadas paramétricas. A esta cuadratura se la denomina integración completa. Se puede demostrar que para cumplir esta condición para un elemento cuadrilateral ó para un hexaedro trilineal se deben utilizar al menos dos puntos de integración en cada dirección, luego el número total de puntos de integración es de cuatro puntos para el elemento cuadrilátero y de ocho puntos para el hexaedro lineal. Debemos notar que el costo computacional para integrar numéricamente las matrices de elemento es directamente proporcional al número total de puntos de integración, por lo tanto, podría pensarse que usando un número menor de puntos que los mínimos requeridos, a lo cual se lo denomina integración reducida, igual se obtendría una buena aproximación a mucho menor costo. De hecho para el elemento cuadrilateral y el hexaedro la integración reducida solo requiere un punto de integración, con los que los costos de integración se reducen a la cuarta parte en cuadriláteros y a la octava parte en hexaedros. Sin embargo este tipo de integración presenta inestabilidades de malla, esto es, la solución aparece contaminada por una función oscilatoria que oculta totalmente la solución correcta. Esto es debido a la presencia de singularidades en la matriz de rigidez. Por lo tanto, no es recomendable el uso de elementos con integración reducida. Algunos programas comerciales dan la opción de usar elementos con integración reducida estabilizada, esto es, usan integración reducida pero mediante alguna técnica eliminan las inestabilidades. REFERENCIAS [1] J.N Reddy y M.L. Rasmussen, Análisis Matemático Avanzado, Limusa, 1992. [2] O.C. Zienkiewicz and R.L. Taylor, The Finite Element Method, McGraw Hill, Vol. I., 4 edn., McGraw-Hill, London, 1989.