MEF-apunte_CA_v2013

Vista previa del material en texto

El Método de los Elementos Finitos (MEF)
1 Introducción
El Método de los Elementos Finitos (MEF) es un “nuevo” método numérico que se suma al tradicional Método
de Diferencias Finitas (MDF) empleado por ingenieros en análisis estructural, resistencia de materiales, mecánica
de fluı́dos, electromagnetismo, etc.
Las ideas básicas del método surgieron a partir de los avances en el área del análisis estructural relacionado
con la aviación (Hrenikoff en 1941 para problemas de elasticidad y Courant en 1943 para problemas de torsión).
Después, a mediados de los años 50 aparece Turner desarrollando matrices de rigidez para la solución de
problemas de elasticidad en barras y vigas, entre otros elementos.
Con grandes logros y siguiendo los pasos de Turner, las Corporaciones MacNeal-Schwendler and Computer
Sciences elaboran en la NASA el primer código de importancia para el análisis de elementos finitos, llamado
NASTRAN (usado en la industria aeroespacial y en áreas de la ingenierı́a civil).
Pero no fue hasta 1960 cuando Clough utilizó por primera vez el término de elemento finito y en 1967 fue
publicado el primer libro de Elemento Finito por Zienkiewicz y Chung poniendo de manifiesto la versatilidad
del método, sus fuertes bases matemáticas y diversas aplicaciones.
Con los grandes avances tecnológicos que se han logrado en el área de la computación y sobre todo en los
sitemas de diseño asistido por computadoras, ahora es relativamente más fácil la modelización de prototipos,
en los cuales se pueden tener geometrı́as y superficies complicadas e irregulares, con aplicaciones de carga en
forma especı́fica para el estudio preciso de esfuerzos internos y lograr una modelización ajustada a los perfiles
y estructuras que se emplean teniendo en consideración ciertas caracterı́sticas como el cambio de secciones,
estructuras huecas con pared delgada y con caracterı́sticas en secciones transversales muy especı́ficas. Por ejem-
plo, la simulación de deformaciones de vehı́culos por impacto y el análisis de mecanismos de transportación
ósea para reducción de fracturas se pueden citar como dos aplicaciones concretas en áreas bien diferentes.
Actualmente se cuenta con una gran cantidad de paquetes computacionales que permiten realizar cálculos
con elementos finitos. Sin embargo, el manejo correcto de este tipo de herramientas exige un profundo
conocimiento no sólo del material con el que se trabaja, sino también de los principios en los que se basa
el método. Sólo en este caso se estará en condiciones de garantizar que los resultados obtenidos en los análisis
se ajustan a la realidad.
2 Fundamentos de los elementos finitos
Uno de los fundamentos del MEF está basado en la discretización de los cuerpos en estudio. De esta manera,
al igual que en el MDF, el método de los elementos finitos require de un espacio geométrico (o dominio) para
ser dividido en subregiones formando una red o malla.
En el MDF, la malla consiste de filas y columnas de lı́neas ortogonales. En cambio en el MEF cada división
es única y no necesariamente ortogonal, lo cual es una ventaja ante sistemas de geometrı́a irregular.
Cada división o subregión (subintervalos en IR, triángulos o cuadriláteros en IR2, tetraedros o he xaedros
en IR3 ) se llama elemento, y en cada uno de éstos se resuelve la ecuación diferencial, para luego generar la
solución total ensamblando las soluciones individuales.
1
La habilidad del MEF para representar con una colección de elementos finitos a dominios de
geometrı́a irregular, es lo que hace de este método una herramienta de mucho valor para resolver problemas a
valores en el borde (PVB) en diversos campos de la ingenierı́a.
Por otra parte, ası́ como el MDF consiste en discretizaciones convencionales, donde las derivadas se aprox-
iman por cocientes en diferencias (cocientes incrementales), el MEF basa su idea en la aproximación de la
función solución. El factor esencial es que la integral de una función se puede escribir como suma de integrales
en dominios disjuntos cuya unión es el dominio original, entonces se puede hacer un análisis local del problema
y haciendo divisiones del dominio en dominios suficientemente pequeños se pueden elegir funciones sencillas
que sean adecuadas para una buena representación del comportamiento de la solución.
Veamos con un ejemplo unidimensional sencillo la implementación de las ideas expuestas.
3 El MEF para problemas unidimensionales
Consideremos el siguiente problema a valores en el borde (P) que modela los desplazamientos axiales de una
barra de longitud unitaria fija en los extremos y sometida a una carga tangencial, o la temperatura en estado
estacionario en un cilindro unidimensional con fuente interna de calor y temperatura cero en los extremos, o
potencial eléctrico en cualquier punto entre dos placas en un medio homogéneo:
(P )

u′′(x) = f(x) , 0 < x < 1
u(0) = 0
u(1) = 0
Idea : A partir de una partición del intervalo (0, 1), aproximar la solución u = u(x) por una función U(x)
que sea una combinación lineal de funciones “adecuadas”, Tj , llamadas funciones tests (o funciones de base).
La aproximación más simple es considerar funciones lineales a trozos definidas en [0, 1]. De esta manera U(x)
es una poligonal que aproxima a u(x)
0 0.2 0.4 0.6 0.8 1
0
0.2
0.4
0.6
0.8
1
1.2
x
y
exacta
aprox.
2
¿Cómo elegir Tj ?
(1o) Se define una partición 0 = x0 < x1 < x2 < ... < xn+1 = 1 del intervalo (0, 1) en n+ 1 subintervalos
(o elementos finitos) Ij = (xj−1, xj) de longitud hj = xj − xj−1. El número h = max {hj} será la
medida de cuan fina es la partición.
(2o) Se construyen las funciones Tj que permitirán expresar de manera única a cualquier función lineal a
trozos que aproxime a la solución del problema (P) y satisfaga sus condiciones de borde. A partir de
estas consideraciones se definen las siguientes n funciones lineales a trozos (una por cada nodo interior
de la malla) y tales que
Tj(xi) =
{
1 si i = j
0 si i 6= j
Ejercicio 1 : Representar graficamente las funciones Tj para una partición del (0, 1) en 4 subintervalos de
igual longitud (partición equiespaciada).
Observación: Notemos que las Tj asi definidas son linealmente independientes y por lo tanto constituyen
una base del espacio vectorial que generan, Vh = gen {T1, T2, ..., Tn} . De esta manera cualquier función
U ∈ Vh tiene una representación única de la forma
U(x) =
∑n
j=1 ujTj(x) donde uj = U(xj)
¿Cómo calcular los coeficientes uj ?
Estos n coeficientes uj son precisamente los valores de la función “aproximante” evaluada en los n nodos
interiores.
(1o) Se reformula (P) obteniendo la llamada forma débil o variacional (Galerkin):
Se multiplica la ecuación diferencial por una función de base arbitraria Ti,
Ti(x)u
′′(x) = Ti(x)f(x) con Ti(0) = Ti(1) = 0 entonces,∫ 1
0
Ti(x)u
′′(x)dx =
∫ 1
0
Ti(x)f(x)dx , integrando por partes resulta
Ti(x)u
′(x)|10−
∫ 1
0
T ′i (x)u
′(x)dx =
∫ 1
0
Ti(x)f(x)dx , y como Ti satisface las condiciones de borde del
problema (P), el primer término se anula y ası́ se obtiene la
Forma débil o variacional : (D)
∫ 1
0
T ′i (x)u
′(x)dx = −
∫ 1
0
Ti(x)f(x)dx
De esta manera se ha demostrado que si u es solución del problema diferencial (P) entonces, u satisface
la forma débil (D) para toda Ti ∈ Vh.
3
Ejercicio 2: Demostrar que si u′′ y f son continuas también vale el recı́proco, es decir, si u es solución de
(D) para toda Ti ∈ Vh entonces u es solución de (P).
(2o) En (D) se reemplaza u por su aproximación U obteniendo la matriz de rigidez y el vector de cargas:
En efecto,
∫ 1
0
T ′i (x)
 n∑
j=1
ujTj(x)
′ dx = −∫ 1
0
Ti(x)f(x)dx , que se puede escribir en la forma
n∑
j=1
uj
∫ 1
0
T ′i (x)T
′
j(x)dx = −
∫ 1
0
Ti(x)f(x)dx para i = 1, 2, ...n,
que no es más que un sistema lineal de n ecuaciones con n incógnitas u1, u2, ..., un.
En forma matricial ese sistema puede escribirse como KU = F donde K = (Kij) es una matriz
n×ncon elementos Kij =
∫ 1
0
T ′i (x)T
′
j(x)dx , F = (Fi) es un vector columna n× 1 con elementos
Fi = −
∫ 1
0
Ti(x)f(x)dx y U = (u1, u2, ..., un)t es el vector de las incógnitas, que son precisamente
los valores aproximados de la variable de estado u en cada uno de los nodos interiores del dominio (0, 1).
La matriz K es conocida como matriz de rigidez y F como vector de carga terminologı́a que proviene
de las primeras aplicaciones del MEF en mecánica estructural.
Observaciones:
• Al ser las funciones Ti lineales a trozos, sus derivadas T ′i no son más que las pendientes de los corre-
spondientes segmentos de recta que las conforman.
• K es una matriz simétrica pues Kij = Kji (son los productos de las correspondientes pendientes).
• En virtud de la elección hecha de las funciones de base Tj la matriz K es tridiagonal, es decir, Kij = 0
si |i− j| > 1 (pues Ti y Tj tienen soportes disjuntos).
• Puede demostrarse que K es invertible con lo cual el sistema lineal tiene única solución U = K−1F
Estimación del error:
Si u′′ es continua en [0, 1] una estimación del error viene dada por
|u(x)− U(x)| ≤ h max
0≤x≤1
∣∣u′′(x)∣∣ para todo x ∈ [0, 1] ,
que demuestra la convergencia del método porque el error tiende a cero cuando el diámetro h de la partición
tiende a cero (pues al ser u′′ continua en (0, 1) resulta acotada).
4
Ejercicio 3:
(a) Aplicar el MEF con una partición equiespaciada de 4 elementos y funciones de base lineales, para obtener
la solución aproximada del siguiente problema a valores en el borde.
(P )

u′′(x) = −1 , 0 < x < 1
u(0) = 0
u(1) = 0
(b) Hallar la solución exacta de (P) y comparar con los valores aproximados.
(a) Para calcular los elementos de la matriz de rigidez, es útil representar graficamente las correspondientes
funciones tests lineales, T1, T2, T3 para h = 0.25. Luego calculamos las correspondientes integrales:
K11 =
∫ 1
0
[
T ′1(x)
]2
dx =
∫ 0.25
0
dx
(0.25)2
+
∫ 0.5
0.25
dx
(−0.25)2
=
0.5
(0.25)2
= 8
Análogamente, K22 = K33 = 8 ( pues la partición es equiespaciada) .
K12 = K21 =
∫ 1
0
T ′1(x)T
′
2(x)dx =
∫ 0.5
0.25
[
−1
0.25
] [
1
0.25
]
=
−0.25
(0.25)2
= −4
Análogamente, K23 = K32 = −4 ( pues la partición es equiespaciada) .
K13 = K31 = 0 pues T1 y T3 tienen soportes disjuntos (la matriz es tridiagonal).
Las componentes del vector de carga serán iguales pues el término fuente es constante, f(x) = −1, y la
partición es equiespaciada, en efecto
para i = 1, 2, 3 Fi = −
∫ 1
0
(−1)Ti dx = 1/4 ( área bajo Ti para cualquier i).
Por lo tanto los valores aproximados de u en los nodos interiores, x1 = 0.25, x2 = 0.5, y x3 = 0.75,
son
U =
4
 2 −1 0−1 2 −1
0 −1 2
−1 1
4
 11
1
 =
 3/321/8
3/32
 =
 u1u2
u3

(b) Hallamos la solución exacta integrando dos veces
u′′(x) = −1 ⇒ u′(x) = −x+ c ⇒ u(x) = −x2/2 + cx+ d.
Luego por las condiciones de borde se tiene d = 0 y c = 1/2 ⇒ u(x) = −x2/2 + x/2.
Notamos que los valores exactos u(1/4) = 3/32, u(1/2) = 1/8, u(3/4) = 3/32, coinciden con los
hallados aplicando el método de los elementos finitos.
5
A fines comparativos se ilustra en la siguiente figura, la solución analı́tica exacta del problema diferencial
(P), u = u(x), conjuntamente con la solución aproximada obtenida con MEF, la poligonal U = U(x),
0 0.2 0.4 0.6 0.8 1
0
5 · 10−2
0.1
x
y
exacta
aprox.
VER más ejercitación en Guı́a de Problemas complementarios (problemas 11 y 12).
4 El MEF para problemas bidimensionales
Se puede ver en la sección 6.7 (pág 267 a 274) de la bibliografı́a básica “Applied Partial Differential Equations”
de Richard Haberman (Cuarta edición 2004) el procedimiento, totalmente análogo al caso unidimensional, que
se puede seguir para la construcción de las funciones de base lineales y la deducción de la forma débil, matriz
de rigidez y vector de carga correspondiente a un problema del tipo
(P )
{
∇2u = f(x, y) para (x, y) ∈ R
u(x, y) = 0 para (x, y) ∈ ∂R
donde con R se designa a cualquier región poligonal en IR2 y ∂R denota su frontera (curva cerrada)
BIBLIOGRAFIA:
• Reddy J.N. “An introduction of the Finite Element Method” - Mc Graw Hill - 1984
• Haberman R. “Applied Partial Differential Equations” - Prentice Hall - 2004
Material elaborado por A. Frausin - Agosto 2012
6