Introducción Introducció n al Método de los Elementos Finitos
Iñaki Zudaire Ingeniero Industrial MIPRO Enginyeria
Introducción al Método de los Elementos Finitos
Índice 1
Conceptos básicos del cálculo matricial 1.1
Sistemas discretos
1 1
1.1.1
Introducción
1
1.1.2
Conceptos básicos del análisis matricial de estructuras de barras
1
1.1.3
Analogía con el análisis matricial de otros sistemas discretos
3
1.1.4
Etapas básicas del análisis matricial de un sistema discreto
4
1.1.5
El principio de los trabajos virtuales
4
1.1.6
Condiciones de contorno
5
El método de los elementos finitos
7
1.2
1.2.1
Introducción
7
1.2.2
Definición
8
i
Introducción al Método de los Elementos Finitos
1 Conceptos básicos del cálculo matricial 1.1 Sistemas discretos 1.1.1 Introducción En numerosas ocasiones nos enfrentamos con el problema de analizar un sistema tipo malla compuesto de una serie de ‘elementos’ diferentes, físicamente diferenciables, conectados por sus extremidades o ‘nudos’ y sometidos a un conjunto de ‘acciones’, en el sentido más amplio de la palabra, normalmente externas al sistema. A éstos sistemas los llamaremos sistemas discretos, y en el mundo de la ingeniería existen muchos ejemplos. Relacionados con las estructuras, por ejemplo, podemos considerar los pórticos, celosías, forjados, etc. En otras áreas de la ingeniería tenemos ejemplos de este tipo de sistemas en las redes hidráulicas y eléctricas, en los métodos de optimización de la producción, en los sistemas de organización del transporte, etc. La mayoría de los sistemas discretos pueden analizarse utilizando técnicas de cálculo matricial muy similares. Concentrándonos en el problema de cálculo de estructuras, presentaremos seguidamente de forma sucinta las ideas del cálculo matricial de estructuras de barras, que serán de gran utilidad como introducción a la metodología del análisis de estructuras por el método de los elementos finitos.
1.1.2 Conceptos básicos del análisis matricial de estructuras de barras Los métodos de cálculo de estructuras de barras más potentes actuales utilizan técnicas de análisis matricial . No obstante, en algunos casos particulares es posible obtener una representación analítica del comportamiento de la estructura. Las ecuaciones matriciales de una estructura de barras se obtienen a partir del estudio del ‘equilibrio’ de las diferentes barras que la componen. Así pues, consideraremos el caso de una barra de longitud l sometida únicamente a fuerzas axiles como las de la figura. l A: Área de la sección
R2
1’
2’
u1
R1
N=σA
u2
Figura 1. 1
Se deduce de la Resistencia de Materiales que la deformación en cualquier punto de la barra es igual al alargamiento relativo de la misma, es decir
=
ε
Δl l
=
u 2 − u1 l
(1.1)
donde u1 y u2 son los desplazamientos de los extremos 1 y 2 de la barra respectivamente. Por otra parte, la tensión axial σ está relacionada con la deformación ε por la ley de Hooke y u − u1
= E ⋅ ε = E ⋅ 2
σ
l
(1.2) 1
Introducción al Método de los Elementos Finitos
donde E es el módulo de elasticidad del material de la barra. Por integración de las tensiones sobre la sección transversal del área A se obtiene el esfuerzo axil N que se transmite a través de los nudos a las barras adyacentes. Suponiendo que el material es homogéneo se tiene N = A ⋅ σ = ( E ⋅ A)
u 2 − u1 l
(1.3)
Finalmente, estableciendo el equilibrio de las fuerzas axiales R 1 y R2 actuantes en los extremos de la barra, se tiene R1 = − R2 = N = ( E ⋅ A)
k =
u 2 − u1 l
(u 2 − u1 ) = k
(1.4)
E ⋅A l
donde
. La equación anterior puede escribirse en forma matricial como
⎧ R1 ⎫ ⎡ 1 − 1⎤ ⎧u1 ⎫ ⋅a ⎬ = k ⎢ ⎥ ⎨ ⎬ = K ⎩ R2 ⎭ ⎣− 1 1 ⎦ ⎩u 2 ⎭
q=⎨
(1.5)
donde K se denomina matriz de rigidez de la barra y es función únicamente de la geometría de la misma (l, A) y de sus propiedades mecánicas (E ), y a y q son los vectores de desplazamiento y de fuerzas de los nudos de la barra, respectivamente. La ecuación presentada és la expresión matricial de equilibrio de la barra aislada. Si además actuara sobre la barra una fuerza uniformemente distribuida por unidad de longitud de intensidad b , la ecuación del equilibrio se modifica repartiendo el efecto total de dicha fuerza en partes iguales en cada nudo como q
(e)
⎧ R1( e ) ⎫ ( e) ⎡ 1 − 1⎤ ⎧u1(e ) ⎫ (b ⋅ l ) ( e ) = ⎨ ( e ) ⎬ = k ⎢ ⎥ ⎨u (e ) ⎬ + 2 1 1 − R ⎣ ⎦⎩ 2 ⎭ ⎩ 2 ⎭
⎧1⎫ ( e) (e) (e ) ⎨ ⎬ = K a + f ⎩1⎭
(1.6)
(b ⋅ l ) ( e ) ⎧1⎫ f = ⎨⎬ 2 ⎩1⎭ es el vector de fuerzas que actúan en los nudos de la barra debidas a donde (e)
la carga distribuida. En la expresión anterior, y con intención de generalizarla, el índice e indica que los valores se refieren a una barra particular de una determinada estructura. La expresión de equilibrio de una estructura compuesta de barras como la considerada se obtiene a partir de la sencilla regla que expresa que la suma de las fuerzas en un nudo, debidas a las diferentes barras que concurren en el mismo, es igual a la fuerza exterior que actúa en dicho nudo. En forma matemática ne
∑1 R ( ) = R e
i
e=
exterior j
(1.7)
donde el sumatorio se extiende a todas las barras ne que concurren en el nudo de numeración (e) global j . Sustituyendo en dicha expresión los valores de las fuerzas de extremo de cada barra R i en función de los desplazamientos de los nudos a través de la ec. 1.6, se obtiene la ecuación matricial de equilibrio global de la estructura como
2
Introducción al Método de los Elementos Finitos
K ⎡ K 11 12 ⎢ K K 22 ⎢ 21 ⎢ M M ⎢ K n1 n2 ⎣ K
K 1n ⎤ ⎧ u1 ⎫
⎧ f 1⎫ ⎥ ⎪ ⎪ ⎪ ⎪ K 2 n ⎥ ⎪u 2 ⎪ ⎪ f 2⎪ ⎨ ⎬=⎨ ⎬ M ⎥⎪ M ⎪ ⎪ M ⎪ ⎥⎪ ⎪ ⎪ ⎪ K nn ⎦ ⎩u n ⎭ n⎭ ⎩ f
K L O L
K ⋅ a = f
(1.8a) (1.8b)
donde K es la matriz de rigidez de la estructura y a y f son, respectivamente, los vectores de desplazamientos y de fuerzas exteriores de todos los nudos de la estructura. El proceso de obtención de las ecuaciones (1.8) recibe el nombre de ensamblaje . La resolución de las mismas proporciona los valores de los desplazamientos en todos los nudos de la estructura a partir de los cuales se pueden conocer los esfuerzos internos de las barras.
1.1.3 Analogía con el análisis matricial de otros sistemas discretos Los pasos explicados entre las ecuaciones 1.1 y 1.8 son muy similares para la mayoría de los sistemas discretos. Así por ejemplo en el caso de una malla eléctrica, el estudio de un elemento aislado (resistencia) proporciona, de acuerdo con la ley de Ohm, la siguiente relación entre los voltajes y las intensidades que entran por cada nudo (e) (e) I 1 = − I 2 =
1 R
(e )
( e) ( e) (e) (e) ( e) − V − V (V 1 2 ) = k (V 1 2 )
(1.9)
Se observa la analogía de dicha ecuación a la 1.4 para la barra, sin más que intercambiar los conceptos de intensidad y voltaje por fuerza y desplazamiento, y el inverso de la resistencia R por E ⋅A l . La expresión 1.9 puede escribirse en formato matricial como (e) (e) ⎧ I 1 ⎡ 1 − 1⎤ ⎧V 1 ⎫ 1 ⎫ ⎥ ⎨ ( e ) ⎬ = ⎨ I e ⎢ (e) ⎬ R ( ) ⎣− 1 1 ⎦ ⎩V 2 ⎭ ⎩ 2 ⎭
(1.10)
La ‘regla de ensamblaje’ es la conocida ley de kirchhoff, que establece que la suma de las intensidades de corriente que concurren en un nudo es igual a cero: ne
∑1 I ( ) = I e
i
exterior j
e=
(1.11)
exterior donde I es la intensidad que entra en el nudo de numeración global j desde el exterior de la j red. Puede comprobarse la analogía de dicha ecuación con la (1.7) para barras.
En el estudio de redes de tuberías se podría hacer una demostración análoga a las planteadas. En el caso de un tramo de tubería, el sistema de ecuaciones matriciales que se plantea tiene la forma
1 − 1⎤ ⎧h1( e) ⎫ ⎧q1( e ) ⎫ k ⎢ ⎥ ⎨ h ( e ) ⎬ = ⎨q ( e ) ⎬ 1 1 − ⎣ ⎦⎩ 2 ⎭ ⎩ 2 ⎭ (e) ⎡
(1.12)
(e) donde k es un coeficiente que depende de la rugosidad de la tubería y de las alturas piezométricas de los nudos, lo que implica que las matrices K (e) no están formadas por constantes sino por funciones conocidas de a (e).
3
Introducción al Método de los Elementos Finitos
En este caso se puede incluir una aportación de caudal uniforme por unidad de longitud de tubería. En este caso, se podría escribir la ecuació (1.6) de manera idéntica para este caso, (e) siendo la fuerza b la aportación de caudal. La regla de ensamblaje se obtiene por la simple condición de equilibrio entre los caudales que concurren en un nudo y el caudal aportado desde el esterior al nudo, es decir ne
∑1 q ( ) = q e
i
exterior j
e=
(1.13)
Se puede deducir fácilmente la analogía de las expresiones anteriores con las correspondientes para estructuras de barras y mallas eléctricas. Las ecuaciones de equilibrio global de una red hidráulica son por tanto idénticas a las (1.8), teniendo en cuenta que la matriz K es de naturaleza no lineal y para sus solución es necesario utilizar métodos iterativos.
1.1.4 Etapas básicas del análisis matricial de un sistema discreto De todo lo anterior se deduce que en el análisis de un sistema discreto intervienen las siguientes etapas: a) Definición de una malla de elementos discretos conectados entre sí por nudos todos ellos convenientemente numerados. Cada elemento e tiene asignadas unas propiedades geométricas y mecánicas conocidas. Todas estas características constituyen los datos del problema y conviene definirlos de la manera más automática posible b) Cálculo de las matrices de rigidez K(e) y los vectores de fuerzas nodales f (e) de cada elemento del sistema. c) Ensamblaje y resolución de la ecuación matricial de equilibrio global para calcular los valores de las incógnitas en los nudos. d) A partir de los valores de las incógnitas en los nudos obtener información sobre otros parámetros de interés del sistema. Todos los resultados deben presentarse con la mayor claridad, i de forma gráfica si es posible, para facilitar la toma de decisiones sobre el diseño.
1.1.5 El principio de los trabajos virtuales Una de las etapas fundamentales del cálculo matricial de estructuras es la obtención de la ecuación de equilibrio de la barra aislada que relaciona las fuerzas actuantes en los nudos con los desplazamientos de dichos nudos (ec. 1.5). Dado el caso de estructuras sencillas, como una barra a tracción, dicha ecuación se puede obtener de manera directa a partir de conceptos intuitivos de la Resistencia de Materiales. En el caso de construcciones más complejas hay que utilizar procedimientos más generales. Uno de los más populares se basa en la aplicación del Principio de los Trabajos Virtuales (PTV) que se anuncia como sigue: ‘Una estructura está en equilibrio bajo la acción de un sistema de fuerzas exteriores si al imponer a la misma unos desplazamientos arbitrarios (virtuales) compatibles con las condiciones en los apoyos, el trabajo realizado por las fuerzas exteriores sobre los desplazamientos virtuales es igual al trabajo que realizan las tensiones en la barra sobre las deformaciones producidas por los desplazamientos virtuales.’ Como es bien sabido, el PTV es condición necesaria y suficiente para el equilibrio de toda la estructura o de cualquiera de sus partes. Aplicaremos ahora dicha técnica a la barra a tracció de la figura 1.1. El PTV se escribe en dicho caso como
4
Introducción al Método de los Elementos Finitos
∫∫∫ ∂ (e)
⋅ ⋅ dV = ∂u1( e ) ⋅ R1(e ) + ∂u 2( e ) ⋅ R2( e )
ε σ
V
(1.14)
(e) (e ) donde ∂u1 y ∂u 2 son, respectivamente, los desplazamientos virtuales de los extremos 1 y 2 de la barra de volumen V(e), y ∂ε la correspondiente deformación virtual que puede calcularse en ( e) ( e) función de ∂u1 y ∂u 2 por (1.1) como
∂ε =
∂u 2( e ) − ∂u1( e ) l
(e)
(1.15)
Sustituyendo los valores de σ y δεen la ecuaciones (1.2) y (1.15) en (1.14), e integrando las tensiones sobre la sección tranversal de la barra se tiene
1
∫ l ( ) [∂u 2 l
( e)
( e)
e
− ∂u1(e) ]( E ⋅ A) ( e)
1 l
[u 2( ) − u1( ) ]dx = ∂u1( ) R1( ) − ∂u 2( ) R2( ) e
(e)
e
e
e
e
e
(1.16)
(e) e integrando sobre la longitud de la barra, considerando E y A(e) constantes
E ⋅ A⎞ ⎛ ⎜ ⎟ ⎝ l ⎠
( e)
[
e u1( )
− u 2( e)
]
∂u1( e)
E ⋅ A⎞ ⎛ +⎜ ⎟ ⎝ l ⎠
(e)
[u2( ) − u1( ) ]∂u2( ) = ∂u1( ) R1( ) − ∂u2( ) R2( ) e
e
e
e
e
e
e
(1.17)
Como los desplazamientos virtuales son arbitrarios, el cumplimiento de (1.17) para cualquier ∂u1(e ) y ∂u 2(e ) exige que los términos que multiplican a cada desplazamiento virtual en los dos miembros sean iguales, lo que proporciona el sistema de dos ecuaciones con dos incógnitas siguiente: E ⋅ A⎞ ⎛ ∂u1(e ) : ⎜ (e ) ⎟[u1( e ) − u 2( e ) ] = R1(e ) l ⎝ ⎠ Para
(1.17a)
E ⋅ A⎞ ⎛ ∂u 2(e ) : ⎜ (e ) ⎟[u 2( e ) − u1( e ) ] = R2( e ) l ⎝ ⎠ Para
(1.17b)
que son las relaciones de equilibrio buscadas entre las fuerzas y desplazamientos de los extremos de la barra. Como puede apreciarse, dichas ecuaciones, escritas en forma matricial, coinciden con las (1.4) obtenidas de manera directa. El PTV se utilizará constantemente a lo largo del libro para obtener las ecuaciones matriciales de equilibrio de los diferentes elementos finitos para cada tipología estructural.
1.1.6 Condiciones de contorno En este capítulo no se va entrar en detalles sobre el proceso de solución del sistema de ecuaciones Ka=f, pues éste es exclusivamente un problema de cálculo numérico que puede resolverse utilizando cualquiera de los múltiples procedimientos que existen, y de los que incluso está disponible su programación en ordenador. No obstante, sí haremos una breve introducción sobre el tratamiento de las condiciones de contorno, pues es un tema de interés general. En un modelos discreto las condiciones de contorno son todas aquellas restricciones que se imponen debidas a la interacción del modelo con su entorno. En el campo del análisis estructural, las condiciones de contorno son desplazamientos impuestos, es decir, que en un determinado nodo fijamos un valor para su desplazamiento (valor que será 0 para un empotramiento, o distinto de 0 si se produce un desplazamiento determinado). Implícitamente, al imponer una restricción de desplazamiento en un nodo estamos provocando que en éste aparezca también un fuerza resultante, para compensar el desequilibrio en las fuerzas nodales. 5
Introducción al Método de los Elementos Finitos
En otros tipos de análisis, como por ejemplo el térmico, una condición de contorno puede ser la temperatura en un nodo, o la transmisión de calor por convección con el medio. En el análisis de redes de tuberías podemos imponer un determinado caudal, con lo que condicionamos el resto de caudales del problema. Para describir el tratamiento de las condiciones de contorno utilizaremos un modelo de cálculo estructural, aunque debe quedar claro que el procedimiento y conclusiones son las mismas para cualquier otro tipo de análisis. Consideremos el sistema global de ecuaciones de equilibrio de una estructura: k 11u1 + k 12 u 2 + k 13 u 3 + L + k 1n u n = f 1 k 21u1 + k 22 u 2 + k 23 u 3 + L + k 2 n u n = f 2 M
k n1u1 + k n 2 u 2 + k n 3u 3 + L + k nn u n = f n
(1.19)
i son fuerzas exteriores (nulas o no nulas) o reacciones en puntos con desplazamiento donde f restringido.
Supongamos que un desplazamiento cualquiera, por ejemplo u2 , está impuesto al valor u 2 , es decir u2 = u2
(1.20)
existen dos procedimientos clásicos para introducir dicha condición en el sistema de ecuaciones (1.19): i del a) Se eliminan la fila y la columna segunda (puesto que se trata de u2 ) y se sustituyen las f i − k i 2 u 2 , es decir, el sistema de n ecuaciones con n segundo miembro de (1.19) por f incógnitas se reduce en una ecuación y en una incógnita como sigue:
k 11u1 + k 12 u 2 + k 13 u 3 + L + k 1n u n = f 1 − k 12 u 2 k 31u1 + k 32 u 2 + k 33 u 3 + L + k 3n u n = f 3 − k 32 u 2 M
k n1u1 + k n 2 u 2 + k n 3 u 3 + L + k nn u n = f n − k n2 u 2
(1.21)
Una vez calculados los u1, u2 , u3, ... u n, el valor de la reacción f 2 (en el caso de que no exista una fuerza exterior aplicada en el nudo 2) se obtinene por f 2 = k 21u1 + k 22 u 2 + k 23 u 3 + ... + k 2nu n
(1.22)
Si el valor restringido para u2 es cero , el procedimiento es el mismo, pero entonces los valores de f i quedan inalterados y el valor de f 2 se obtiene por (1.22) prescindiendo del término que afecta a u2 . b) Otro procedimiento bastante utilizado y que no precisa modificar apenas el sistema de ecuaciones original, consiste en añadir un coeficiente de valor alto al término de la diagonal principal de la fila correspondiente al desplazamiento prescrito, y reemplazar el segundo miembro de la ecuación de dicha fila por el valor del desplazamiento restringido multiplicado 15 por dicho coeficiente. Es decir, si de nuevou 2 = u 2 , sustituiríamos k por k +10 k (por 22
22
22
10 k 22 * u 2 , quedando el sistema de ecuaciones de la ejemplo), y el valor de f 2 por siguiente forma 15
6
Introducción al Método de los Elementos Finitos
k 11u1 + k 12 u 2 + k 13 u 3 + L + k 1n u n = f 1 k 21u1 + (1 + 10
15
15 )k 22 u 2 + k 23 u3 + L + k 2 n u n = 10 k 22 * u 2
M
k n1u1 + k n 2 u 2 + k n 3u 3 + L + k nn u n = f n
(1.23)
15 De esta manera, la segunda ecuación, al ser 10 k mucho mayor que el resto de los coeficientes, 22 equivale a 15 10 15 k 22 u 2 = 10 k 22 u 2
ó
u2 = u2
(1.24)
que es la condición restringida. Con este procedimiento la condición se impone de forma natural en la solución del sistema de ecuaciones con modificaciones mínimas. El valor de la reacción f se calcula a posteriori con la ecuación (1.22). 2
1.2 El método de los elementos finitos 1.2.1 Introducción Con la excepción de los problemas planteados hasta el momento, a los que hemnos llamado discretos, la mayor parte de sistemas en ingeniería son de naturaleza continua y, por tanto, su comportamiento no puede expresarse en función de un número pequeño de variables discretas. Un análisis riguroso de dichas estructuras precisa la integración de las ecuaciones diferenciales que expresan el equilibrio de un elemento diferencial genérico de las mismas. Aunque los sistemas continuos son inherentes tridimensionalmente en algunos casos su comportamiento puede describirse adecuadamente por modelos matemáticos uni, bi o tridimensionales. Así ocurre, por ejemplo, con los problemas de flexión de placas, en los que el análisis se limita al estudio de la deformación del plano medio de la placa, y con todas las estructuras en las que puede hacerse uso de las hipótesis simplificativas de la elasticidad bidimensional o de revolución. El método de los elementos finitos es hoy en dia el procedimiento más potente para el análisis de sistemas de carácter uni, bi o tridimensional sometidos a las acciones exteriores más diversas. La gran analogía existente entre los conceptos del análisis matricial y los del método de los elementos finitos facilitan en gran manera el estudio de éste a los técnicos con dominio de las ideas de cálculo matricial de estructuras tratadas en apartados anteriores. En este curso se presenta la aplicación del método de los elementos finitos al análisis de sistemas de diferentes tipologías. Inicialmente, y a modo de planteamiento general del métod, nos centraremos en el estudio de estructuras de naturaleza tanto discreta como continua. De este modo, con carácter fundamentalmente didáctico, se introducirán los conceptos básicos del método de los elementos finitos en su aplicación a simples estructuras de barras y vigas. Más adelante, una vez conocido el procedimiento interno del cálculo por el método de los elementos finitos, se plantearan las ecuaciones de equilibrio para el análisis de procesos térmicos y dinámicos, y se apuntaran las bases para el cálculo de sistemas no lineales. Es importante destacar desde un principio las analogías entre las etapas básicas del análisis matricial y el de un sistema cualquiera por el método de lo lementos finitos. Desde el punto de vista práctico del cálculo, la característica más atractiva del método, y quizas también la más peligrosa, estriba en el hecho de que es un método aproximado. En las manos de un técnico cuidadoso y experto es un procedimiento muy útil para obtener información sobre el comportamiento de sistema complejos para los que no existen soluciones analíticas disponibles. 7
Introducción al Método de los Elementos Finitos
No obstante, su carácter aproximado le confiere un cierto riesgo, y su utilización, si no se posee una experiencia previa, debe efectuarse con precaución.
1.2.2 Definición En el caso de los sistemas discretos, como ya se ha demostrado, existen soluciones analíticas para cada tipo de problema que se puede plantear, de modo que su resolución puede considerarse casi directa. En el caso de un sistema complejo (o continuo) no existen este tipo de soluciones por lo que obtener información de éllos resulta más difícil y costoso. Para vencer la infranqueabilidad que supone la solución de problemas continuos reales, como ya se ha comentado, ingenieros y matemáticos han ido proponiendo a través de los años, diversos métodos de discretización, para lograr transformar el problema continuo y desconocido en otro discreto y del que existe una solución conocida. El proceso de discretización del sistema consiste en subdividirlo en un número finito de componentes bien definidos, a los que llamaremos elementos finitos, y que están conectados entre sí a través de sus estremos o nodos. De este modo el planteamiento matemático-analítico de las ecuaciones del equilibrio de dichos elementos da lugar a la formulación discreta del problema. No obstante, dicha simplificación hace necesario efectuar alguna aproximación de tal naturaleza que quepa esperar de la misma se acerque, tan estrechamente como se quiera, a la solución continua verdadera, a medida que crezca el número de elementos ( o variables discretas) definidos. La existencia de una manera única de abordar los problemas discretos nos lleva a la primera definición del método de los elementos finitos como procedimiento de aproximación de problemas continuos, de tal forma que: a) El continuo se divide en un número finitos de partes (elementos), cuyo comportamiento se especifica mediante un número finito de parámetros b) La solución del sistema completo como ensamblaje de los elementos sigue precisamente las mismas reglas que se aplican a los problemas discretos.
8