INTRODUCCIÓN AL

MÉTODO DE LOS
ELEMENTOS FINITOS
Resumen de clases
Actualización 2014
 ÍNDICE GENERAL

Prólogo v

1. INTRODUCCIÓN 1
1.1. El Método de los Elementos Finitos . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2. El Principio de la Mínima Energía Potencial Total . . . . . . . . . . . . . . . . . 5

2. CÁLCULO VARIACIONAL 11
2.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2. Cálculo variacional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2.1. Multiplicadores de Lagrange . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2.2. Cálculo de máximos o mínimos . . . . . . . . . . . . . . . . . . . . . . . . 14
2.2.3. Condiciones naturales de contorno . . . . . . . . . . . . . . . . . . . . . . 16
2.2.4. La notación variacional . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.3. Ejemplos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.3.1. El problema de la braquistócrona . . . . . . . . . . . . . . . . . . . . . . . 18
2.3.2. El principio de Hamilton . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.3.3. Transmisión del calor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.3.4. Cálculo estructural . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

3. TEOREMAS ENERGÉTICOS 25
3.1. Definiciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.1.1. Tipos de materiales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.1.2. Trabajo de las fuerzas exteriores y energía de deformación . . . . . . . . . 25
3.1.3. Expresiones para solicitación axil . . . . . . . . . . . . . . . . . . . . . . . 26
3.2. Principio de equivalencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.3. Teorema de los desplazamientos virtuales . . . . . . . . . . . . . . . . . . . . . . . 29
3.4. Teorema de las tensiones virtuales . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.5. El principio de la mínima energía potencial total . . . . . . . . . . . . . . . . . . 33
3.6. Ejemplos de aplicación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.6.1. Barra solicitada axilmente . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.6.2. Viga solicitada a flexión . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

4. El MÉTODO DE RITZ 41
4.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
4.2. El método de Ritz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
4.3. Ejemplos de aplicación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

i
ÍNDICE GENERAL El Método de los Elementos Finitos

4.3.1. Resolución en primer orden de una viga sometida a flexión . . . . . . . . . 43

4.3.2. Barra solicitada axilmente . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.3.3. Viga a flexión en teoría de segundo orden . . . . . . . . . . . . . . . . . . 47
4.3.4. Transmisión de calor: caso unidimensional . . . . . . . . . . . . . . . . . . 49

5. ELEMENTO DE BARRA 53
5.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
5.2. Elemento de barra con dos nodos . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
5.3. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

6. ELEMENTO DE VIGA 59
6.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
6.2. Cálculo de la matriz de rigidez . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
6.3. La hipótesis de Timoshenko . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
6.4. Determinación del coeficiente κ . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

7. ELEMENTO DE ESTADO PLANO 69

7.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
7.2. Estados planos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
7.2.1. Estado plano de tensión . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
7.2.2. Estado plano de deformación . . . . . . . . . . . . . . . . . . . . . . . . . 70
7.3. Elementos de estados planos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
7.3.1. Elemento cuadrilátero de estado plano de tensión . . . . . . . . . . . . . . 70
7.3.2. Elemento triangular de estado plano . . . . . . . . . . . . . . . . . . . . . 73
7.4. Simetría de revolución . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

8. ELEMENTO DE PLACA 79
8.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
8.2. Placa plana con Hipótesis de Kirchhoff . . . . . . . . . . . . . . . . . . . . . . . . 80
8.2.1. Obtención de las deformaciones . . . . . . . . . . . . . . . . . . . . . . . . 80
8.2.2. Ecuaciones de equivalencia . . . . . . . . . . . . . . . . . . . . . . . . . . 81
8.2.3. Planteo del Teorema de los Desplazamientos Virtuales . . . . . . . . . . . 82
8.3. Placa plana con Hipótesis de Mindlin . . . . . . . . . . . . . . . . . . . . . . . . . 83
8.4. Elementos isoparamétricos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

9. TEOREMAS VARIACIONALES MIXTOS 89

9.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
9.2. Teorema de Hellinger-Reissner . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
9.2.1. Ecuaciones para vigas sometidas a flexión con teoría de primer orden . . . 93
9.2.2. Ecuaciones de la teoría de segundo orden para vigas . . . . . . . . . . . . 96
9.3. Teorema de Hu-Washizu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98

A. EL USO DE LA COMPUTADORA 99
A.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
A.2. Los métodos numéricos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
A.3. Sistemas de ecuaciones lineales . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
A.3.1. Métodos tradicionales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
A.3.2. Métodos más modernos . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
A.4. Integración numérica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
A.5. Los programas de análisis por elementos finitos . . . . . . . . . . . . . . . . . . . 105
A.5.1. Conceptos generales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
A.5.2. Programas de análisis estructural por elementos finitos . . . . . . . . . . . 106

- ii - Preliminar Versión 2014

 ÍNDICE DE FIGURAS

1.1. Modelo real (a) y modelo discreto (b). . . . . . . . . . . . . . . . . . . . . . . . . 4

1.2. Trabajo de la fuerza P y caminos posibles. . . . . . . . . . . . . . . . . . . . . . . 5
1.3. Configuraciones admisibles y no admisibles. . . . . . . . . . . . . . . . . . . . . . 5
1.4. Energía Potencial Total. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.5. Principio de la Mínima Energía Potencial Total. . . . . . . . . . . . . . . . . . . . 7
1.6. Sistema con tres grados de libertad. . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.7. Diagrama de cuerpo libre. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.8. Sistema estructural continuo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.9. Sistema sometido a cargas axiles. . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.1. Máximo para la función y(x) = sen x en x0 = π2 . . . . . . . . . . . . . . . . . . . . 12

2.2. Funciones y(x), η(x) y y(x) + ε η(x). . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.3. El problema de la braquistócrona. . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.4. Transmisión del calor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.5. Viga empotrada–apoyada. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

3.1. Barra solicitada axilmente. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

3.2. Trabajo de las fuerzas exteriores y trabajo complementario de las fuerzas exteriores. 27
3.3. Energía interna de deformación y energía interna complementaria de deformación. 27
3.4. Energías. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.5. Barra solicitada axilmente. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.6. Energía interna de deformación y trabajo de las fuerza exteriores. . . . . . . . . . 29
3.7. Barra con desplazamientos axiles. . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.8. Energía interna de deformación y trabajo de las fuerzas exteriores complementarios. 32
3.9. Barra empotrada–libre con carga distribuida axil. . . . . . . . . . . . . . . . . . . 34
3.10. Viga sometida a flexión. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

4.1. Viga sometida a flexión. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

4.2. Viga sometida a carga concentrada. . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.3. Diagrama de momento flector. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.4. Viga con carga horizontal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
4.5. Transmisión del calor unidimensional. . . . . . . . . . . . . . . . . . . . . . . . . 49

5.1. Elemento de barra. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

5.2. Sistema con barras axilmente cargadas. . . . . . . . . . . . . . . . . . . . . . . . . 55
5.3. Incógnitas globales. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
5.4. Discretización del sistema e incógnitas locales. . . . . . . . . . . . . . . . . . . . . 56

iii
ÍNDICE DE FIGURAS El Método de los Elementos Finitos

6.1. Elemento de viga. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

6.2. Incógnitas del elemento de viga. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
6.3. Hipótesis de Timoshenko. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
6.4. Elemento de barra. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
6.5. Solución real y proyección ortogonal. . . . . . . . . . . . . . . . . . . . . . . . . . 66

7.1. Elemento de estado plano de tensiones. . . . . . . . . . . . . . . . . . . . . . . . . 69

7.2. Elemento de estado plano de deformaciones. . . . . . . . . . . . . . . . . . . . . . 70
7.3. Estado plano. Elementos rectangulares. . . . . . . . . . . . . . . . . . . . . . . . . 71
7.4. Estado plano. Elementos triangulares. . . . . . . . . . . . . . . . . . . . . . . . . 74
7.5. Funciones A1 , A2 y A3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
7.6. Cargas distribuidas. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76

8.1. Hipótesis de Kirchhoff. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

8.2. Elemento isoparamétrico rectangular. . . . . . . . . . . . . . . . . . . . . . . . . . 85

9.1. Viga sometida a carga concentrada. . . . . . . . . . . . . . . . . . . . . . . . . . . 90

9.2. Desplazamiento propuesto. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
9.3. Momento flector propuesto. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
9.4. Desplazamiento virtual propuesto. . . . . . . . . . . . . . . . . . . . . . . . . . . 91
9.5. Variación de la función momento flector. . . . . . . . . . . . . . . . . . . . . . . . 93
9.6. Diagramas para un elemento de viga genérico. . . . . . . . . . . . . . . . . . . . . 94
9.7. Carga concentrada de compresión. . . . . . . . . . . . . . . . . . . . . . . . . . . 96
9.8. Nodo genérico. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97

A.1. Malla gruesa. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

A.2. Malla más densa. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
A.3. Malla más densa sin modificar la numeración original. . . . . . . . . . . . . . . . 103

- iv - Preliminar Versión 2014

El Método de los Elementos Finitos 0. Prólogo

PRÓLOGO

Este apunte no es otra cosa que un resumen de las clases dictadas durante el cuatrimestre
de la materia MÉTODO DE LOS ELEMENTOS FINITOS en la Facultad de Ingeniería
de la Universidad de Buenos Aires. No pretende ser un libro de texto sino simplemente una guía
ordenada de los distintos temas que se abordan en la materia, de manera que el alumno pueda
luego profundizar los conocimientos mediante la bibliografía existente.

La introducción en el método se hace a partir de la resolución de problemas estructurales

conocidos, comenzando por los casos más sencillos y avanzando hacia problemas de mayor com-
plejidad. Para ello, sin embargo, debe introducirse una rama del análisis matemático que no se
estudia en las parte inicial de la carrera como es el cálculo variacional.

Puesto que el método de los elementos finitos es un método numérico para la resolución
de ecuaciones diferenciales, los resultados que se obtendrán serán solamente una aproximación.
Aún así, en la mayoría de los casos, sólo es posible obtener esta solución aproximada, porque no
existe una solución analítica del problema. El cálculo variacional, junto con el Análisis Funcional,
son las herramientas matemáticas que permiten estudiar la calidad de esta aproximación.

Versión 2014 Preliminar -v-

El Método de los Elementos Finitos 1. INTRODUCCIÓN

CAPÍTULO 1

INTRODUCCIÓN

1.1. El Método de los Elementos Finitos

En los últimos años, sobre todo con la generalización del uso de las computadoras, los
ingenieros y científicos han pasado de aplicar modelos matemáticos y físicos puramente analíticos
a utilizar cada vez más modelos matemático-numéricos. Esto se ha incrementado en los últimos
15 a 20 años, principalmente por el abaratamiento y el exponencial incremento en la capacidad
de procesamiento de las computadoras personales.
De los muchos modelos numéricos usados hoy en día, uno de los más importantes es
el Método de los Elementos Finitos. Este método nació para analizar las tensiones en las
piezas, es decir, como parte de un análisis estructural. Luego se extendió a otros campos, como
la transferencia de calor, la mecánica de fluidos y el electromagnetismo, entre otros campos.
Hoy está en constante desarrollo y sirve como base para una disciplina más amplia que es la
Mecánica Computacional.
Uno de los conceptos básicos del método es la «discretización», es decir, la de convertir
un problema que matemáticamente se puede representar con un modelo continuo en otro de
similares características, pero «discreto». De ahí el concepto de «finito» por contraposición al
«infinito» o, en este caso, al «infinitésimo». Por lo tanto, el objetivo del método no es obte-
ner una solución «completa» al problema sino una aproximación mediante un modelo discreto
(finito) que se acerque lo más posible al modelo tradicional continuo. Esta forma de enfrentar
el problema trae aparejada la imposibilidad de encontrar una solución general al problema, al
modo matemático tradicional, pero permite encontrar soluciones particulares a problemas que
no cuentan con soluciones generales, lo que en cierta forma le da una gran versatilidad.
A partir de de lo anterior, podemos ensayar una primera definición:

El método de los Elementos Finitos es un método por aproximación discreta que

utiliza funciones sencillas como base para aproximar la función solución del problema.
Las funciones sencillas suelen ser funciones de interpolación entre puntos llamados
«nodos».

En el campo del análisis estructural, la utilización del método de los elementos finitos
suele incluir los siguientes pasos:

1. Definir el problema (estructura);

2. Dividir el problema (estructura) continuo en pequeñas partes (discretización en elementos);

3. Definir los puntos iniciales y finales de cada parte (nodos);

Versión 2014 Preliminar -1-

1.1. El Método de los Elementos Finitos El Método de los Elementos Finitos

4. Formular las propiedades de cada elemento;

5. Definir las cargas, tanto en los elementos como en los nodos;
6. Ensamblar los elementos y las cargas;
7. Determinar las condiciones de borde del problema (estructura);
8. Resolver el sistema de ecuaciones lineales;
9. Obtener las deformaciones y las tensiones en cada elemento.
En el caso particular del análisis estructural, los resultados obtenidos sirven para el dimen-
sionamiento y la verificación de los elementos estructurales que componen la estructura principal,
o sea, el problema en estudio. En muchos casos, una vez obtenidos los resultados numéricos se
necesitan programas adicionales para poder visualizar esos resultados de forma amigable para el
ingeniero. Estos programas se conocen como «post–procesadores». En otros casos, para poder dis-
cretizar el problema es necesario un programa adicional, conocido como «pre–procesador». En la
actualidad, la mayoría de los programas de análisis estructural incluyen tanto el pre–procesador
como el postprocesador, que generalmente está asociado a una pantalla gráfica (similar a un
programa de dibujo), que ayuda en la definición de los nodos y los elementos.
Debe tenerse bien en cuenta que el método de los elementos finitos sirve para resolver
problemas mediante procedimientos numéricos, por lo que no pueden obtenerse soluciones «ana-
líticas». De todos modos, la capacidad para resolver problemas prácticos es mucho mayor que
las posibilidades de un método analítico. Esta característica del método hace necesario que para
usar un programa comercial (o no) que lo aplique, se deban cumplir estos requisitos:
Apoyarse en la documentación del programa;
Conocer las bases del método y de los principios que rigen el problema a encarar;
Un criterio ingenieril para analizar y entender los datos a ingresar y los resultados obtenidos.
En consecuencia, es importante conocer las bases del método para entenderlo y de esa
manera usarlo de la forma más eficaz y eficiente posible. Y la única forma de conocerlo es
aprender las bases matemáticas y numéricas junto con las bases del análisis estructural. De allí
que el bagaje teórico es importante como conocimiento fundamental.
Justamente, uno de los puntos más importantes en estos principios teóricos está en los
métodos numéricos necesarios para resolver sistemas de ecuaciones lineales, el cálculo de los
coeficientes de la matriz del sistema de ecuaciones lineales mediante métodos de integración
numérica y el uso de funciones interpolantes para representar las características físicas de los
elementos que constituyen el modelo.
Además de los fundamentos numéricos, también hay que agregar para el caso del análisis
estructural, el aprendizaje de otra de las ramas del cálculo o análisis matemático: el cálculo de
variaciones, sobre el que volveremos más adelante.
El análisis estructural cuenta con modelos matemáticos que pueden resolver problemas
con ayuda de herramientas numéricos. De hecho, los métodos de las incógnitas estáticas y cinemá-
ticas, usados para resolver estructuras hiperestáticas, terminan generando sistemas de ecuaciones
lineales que pueden ser resueltos con suma facilidad por cualquier programa que cuente con mó-
dulos de álgebra lineal numérica. Es usual que en los cursos de grado se inste a los alumnos a
hacer uso de estos programas.
De ahí que para algunos casos, existen métodos tradicionales que son posibles de im-
plementar en modelaciones numéricas y obtener resultados tradicionalmente identificados como
«exactos». Pero estos métodos suelen limitarse a modelos sencillos generalmente planos. La im-
plementación de modelos espaciales es más compleja y es ahí cuando el método de los elementos
finitos es altamente eficiente.

-2- Preliminar Versión 2014

El Método de los Elementos Finitos 1. INTRODUCCIÓN

No debe perderse de vista, de todas formas, algo importante: los resultados que presenta
una programa pueden estar equivocados. El usar un elemento no apropiado, una malla pobre-
mente definida, subestimar el pandeo o la fluencia, una mala definición de los apoyos (condiciones
de borde), entre otros, puede llevar a obtener malos resultados al aplicar un análisis por elemen-
tos finitos. Tampoco es el método universal para resolver numéricamente cualquier problema.
Existen otros métodos igual de calificados para efectuar una análisis estructural, como el método
de las diferencias finitas (que son aplicados con gran efectividad en problemas de cáscaras de
revolución) y el método de los elementos de contorno. Hay una frase muy conocida entre los
ingenieros:
El método de los elementos finitos puede hacer mejor a un buen ingeniero, pero
puede volver peligroso a uno malo.
Es por eso que un programa de análisis por elementos finitos no debería ser usado sin un entre-
namiento adecuado, pues sin una comprensión cabal de sus bases teórica (físicas, matemáticas y
numéricas), además de su funcionamiento, puede ser mal utilizado, justificando la frase anterior.
Como primera aproximación al método, diremos que consiste en una transformación de
un problema continuo en otro discreto, que puede ser resuelto mediante algoritmos numéricos.
De esta forma se pueden resolver problemas que eventualmente no tengan solución analítica (o
sea, solución del problema continuo).
Desde un punto de vista estrictamente matemático–numérico, se transforma la ecuación
diferencial que define el problema en un esquema de tipo algebraico discreto. Existe un modelo
previo para esta transformación y es el Método de las Diferencias Finitas, muy utilizado para
resolver cualquier tipo de ecuación diferencial. Muchos problemas y estructuras han sido resueltos
con ayuda de este método en el pasado y aún se sigue utilizando para algunos problemas en
particular. Pero la formulación de la solución por diferencias finitas es estrictamente matemática,
pues no incorpora ninguna consideración física, más allá de algunas características físicas (el
módulo de Young) y geométricas (el momento de inercia estático).
Básicamente, el Método de las Diferencias Finitas se concentra en la obtención
de los valores de los desplazamientos en puntos predeterminados denominados nodos, que son el
resultado de la discretización del modelo continuo, proceso que suele llamarse como creación de la
malla (en inglés, meshing) del modelo numérico. El método resulta en transformar una ecuación
diferencial en un sistema de ecuaciones lineales, para cuya solución se cuenta con varios algoritmos
numéricos muy conocidos. Pero obtener los parámetros estructurales para el dimensionamiento
o verificación de la estructura, lleva a agregar información o condiciones por fuera del modelo
matemático–numérico que no inciden en el resultado original, y en consecuencia, pueden modificar
las hipótesis que se usaron para la definición de los nodos. Más aún, la definición de la malla
para la resolución de las ecuaciones diferenciales mediante diferencias finitas no siempre resulta
del todo adaptada a las condiciones del problema, por lo tanto muchas veces deben realizarse
modelos sucesivos de aproximación, lo que incide en la rapidez y la convergencia de los resultados.
No siempre se consiguen aproximaciones aceptables. Una de las ventajas del método es que a
medida que esa malla es más densa (mayor cantidad de puntos), la solución numérica converge
a una eventual solución analítica (continua).
A causa de estos problemas y por necesidades de los ingenieros aeronáuticos que interve-
nían en el proyecto y desarrollo de los aviones a reacción, se comenzó a estudiar la posibilidad de
analizar las estructuras con otro enfoque. La idea que privó fue la discretizar la estructura como
en el caso de las diferencias finitas pero ya no desde un punto de vista matemático–numérico
sino que el modelo numérico incluyera completamente las características físicas. En una primera
aproximación tomaron como referencia los modelos de cálculo matricial de estructuras, parti-
cularmente del método conocido cono de las incógnitas geométricas (o de las rigideces), que
justamente se basaba en el cálculo de los desplazamientos en un estructura plana (un pórtico,
por ejemplo), y que podía ser resuelto mediante algoritmos de álgebra lineal numérica, que eran
de «fácil» codificación en las computadoras de la década del 50.

Versión 2014 Preliminar -3-

1.1. El Método de los Elementos Finitos El Método de los Elementos Finitos

La idea que tomaron de este método fue la siguiente: en lugar de discretizar el dominio del
problema estructural mediante puntos o nodos solamente, hacerlo con elementos de dimensiones
pequeñas pero no infinitésimos (de ahí el concepto de elemento finito), y que esos elementos
incluyan todas las características físicas y geométricas propias de la estructura en análisis: ma-
teriales y geometría. Como ejemplo, supongamos una columna soportada desde arriba y cuya
sección transversal varíe según una función ´cuadrática, cúbica o cualquier otra, algo que es
difícil de representar en un modelo por diferencias finitas, pero que mediante la discretización
por elementos pequeños puede hacer sin demasiada complicación, como se ve en la figura 1.1.

Figura 1.1: Modelo real (a) y modelo discreto (b).

El resultado de esta discretización es un sistema de ecuaciones lineales, similar al que se

obtiene por diferencias finitas, pero con ciertas particularidades que lo asemejan más a los modelos
matriciales de resolución de estructuras, porque al discretizar mediante elementos finitos y no
por puntos, el modelo discreto incorpora las características físicas del problema. Y al igual que
en el caso de diferencias finitas, cuanto más densa la malla, mayor convergencia de la solución
numérico con la eventual solución continua.
Otro punto de contacto entre ambos métodos son los nodos. El método de los elementos
finitos también utiliza nodos para la definición del modelo numérico. En el caso del método de
los elementos finitos los nodos sirven para dos cosas: evidenciar los desplazamientos del sistema
estructural, al igual que en el método de las diferencias finitas,; y para definir el «elemento finito»,
es decir, permite la caracterización del elemento, introduciendo un modelo de aproximación de
los desplazamientos entre nodos, e incluye las características mecánicas y geométricas de la
estructura, para luego calcular la matriz de las rigideces, o matriz de rigidez, que será parte
fundamental del sistema de ecuaciones lineales. Esto es muy importante, porque el elemento
incorpora las hipótesis mecánicas de comportamiento estructural en los elementos, algo que el
modelo por diferencias finitas ignora completamente, pues no es necesario para la dicretización
de la ecuación diferencial.
De ahí que la aplicación del Método de los Elementos Finitos requiere que los ingenieros
estén capacitados y tengan los conocimientos fundamentales de matemática y mecánica de los
materiales (la parte física) que incorpora la formulación numérica. Estos conocimientos son la
guía para llevar adelante los pasos que requiere en análisis estructural por elementos finitos:

1. El pre–procesamiento;

2. La solución numérica, y;

3. El post–procesamiento.

En todos los pasos es importante la intervención de ingenieros bien capacitados en el uso

del método, más que en la carga de datos en una computadora. Es por eso que este trabajo apunta
a esto, a explicar las bases teóricas del método, tanto desde la mecánica de materiales (física)

-4- Preliminar Versión 2014

El Método de los Elementos Finitos 1. INTRODUCCIÓN

como de la matemática y la numérica, que capacite a los ingenieros para usar apropiadamente
el método. Y uno de esos conceptos fundamentales es el principio de la mínima energía
potencial total.

1.2. El Principio de la Mínima Energía Potencial Total

Uno de los puntos más importantes para poder aplicar el método de los elementos finitos
es conocer sus bases matemáticas. Para ello, analizaremos primero un principio muy importante
para el cálculo estructural: el Principio de la Mínima Energía Potencial Total.
Definamos como sistema a una estructura y las cargas aplicadas en ella. La configuración
de un sistema es el conjunto de posiciones de todas las partículas de esa estructura. Sabemos
que un sistema es conservativo si todos loas trabajos de las fuerzas internas y externas son
independientes del camino seguido por dichas fuerzas (ver caminos A y B en la figura 1.2).
En un estructura elástica, el trabajo externo es igual en magnitud al cambio de la energía de
deformación. Además, tenemos dos tipos de condiciones de borde: las esenciales (a veces llamadas
principales o forzadas) y las no esenciales (o naturales).

Figura 1.2: Trabajo de la fuerza P y caminos posibles.

Generalmente cunado aplicamos el método de los elementos finitos, las condiciones de

borde esenciales son los desplazamientos y los giros, de ahí que suelen conocerse como condi-
ciones de borde cinemáticas, en tanto que las condiciones de borde no esenciales suelen ser
momentos flectores, esfuerzo de corte o esfuerzo normal, lo que lleva a que se las conozca como
condiciones de borde estáticas. Desde un punto de vista matemático, las condiciones de
borde esenciales involucran a la función y a su derivada primera, en tanto que las condiciones de
borde no esenciales involucran a derivadas de orden superior.

Figura 1.3: Configuraciones admisibles y no admisibles.

Para cualquier sistema es necesario definir las configuraciones admisibles. Una configu-
ración admisible es aquella que cumple con la compatibilidad interna y las condiciones de borde
esenciales. Por ejemplo, en la figura 1.3, la configuración 1 no es admisible por varios motivos:

1. No cumple con las condiciones de desplazamientos y giros nulos en el empotramiento;

2. Viola la compatibilidad interna en la discontinuidad en A y en el vértice B.

Los otros dos desplazamientos son configuraciones admisibles, aún cunado sólo una pa-
rece ser razonable (la 3 ). Es importante recalcar que una configuración admisible no necesita
satisfacer las condiciones de borde no esenciales (o naturales).

Versión 2014 Preliminar -5-

1.2. El Principio de la Mínima Energía Potencial Total El Método de los Elementos Finitos

Todo sistema mecánico conservativo tiene una energía potencial. Ésta, que también se
denomina energía potencial total, incluye a la energía interna de deformación y al trabajo
de las fuerzas externas. Surge, entonces, un principio denominado Principio de la Mínima
Energía Potencial Total:
De todas las configuraciones admisibles de un sistema conservativo, aquellas que
satisfacen la ecuación diferencial de equilibrio hace mínima la Energía Potencial Total
respecto de pequeñas variaciones de los desplazamientos admisibles. Y dado que es
mínimo, es estable.
Para ello tomemos el siguiente ejemplo:

Figura 1.4: Energía Potencial Total.

La energía potencial total la definimos de la siguiente forma:

ΠP = U + W, (1.1)

donde U es la energía interna de deformación y W es el trabajo de las fuerza externas. Para

nuestro caso de la figura 1.4, la energía interna de deformación es:
1
U = k D2 , (1.2)
2
en tanto que el trabajo de las fuerzas externas es:

W = −P D. (1.3)

El trabajo de las fuerzas externas lo consideramos negativo porque al desplazarse la

distancia D pierde energía potencial. En consecuencia, la energía potencial total del sistema de
la figura 1.3 resulta ser:
1
ΠP = k D2 − P D. (1.4)
2
Para obtener la configuración de equilibrio basta con obtener D cunado la energía po-
tencial total sea mínima. Para obtener este mínimo basta con derivar la expresión de la energía
potencial total respecto a D e igualarla a cero:
P
dΠP = (k D − P ) dD = 0 → Deq = . (1.5)
k
La ecuación (k D − P )dD = 0 es una forma del principio de los trabajos virtuales:
«el trabajo de todas las fuerzas intervinientes que están en equilibrio es nulo para pequeños des-
plazamientos admisibles o compatibles con los vínculos». En la figura 1.5 está una interpretación
gráfica del principio de la mínima energía potencial.
Hasta aquí nos hemos ocupado de un sistema de un solo grado de libertad. ¿Pero qué
pasa si aumentamos esos grados de libertad? Analicemos ahora un sistema representado mediante
resortes pero con más grados de libertad, como en la figura 1.6.

-6- Preliminar Versión 2014

El Método de los Elementos Finitos 1. INTRODUCCIÓN

Figura 1.5: Principio de la Mínima Energía Potencial Total.

Figura 1.6: Sistema con tres grados de libertad.

Para este sistema la expresión de la energía potencial total es:

1 1 1
ΠP = k1 D12 + k2 (D2 − D1 )2 + k3 (D3 − D2 )2 − P1 D1 − P2 D2 − P3 D3 . (1.6)
2 2 2

Al tener tres grados de libertad (D1 ; D2 y D3 ), el equilibrio se dará cuando los tres
desplazamientos hagan que ΠP sea mínima. Para obtener los valores de los tres desplazamientos,
hallemos la derivada total de ΠP e igualémosla a cero:

∂ΠP ∂ΠP ∂ΠP

dΠP = dD1 + dD2 + dD3 = 0. (1.7)
∂D1 ∂D2 ∂D3

Para que esta derivada total sea siempre nula, para cualquier desplazamiento pequeño,
se debe cumplir que:

∂ΠP
= 0; (1.8a)
∂D1
∂ΠP
= 0; (1.8b)
∂D2
∂ΠP
= 0. (1.8c)
∂D3

Si calculamos las derivadas parciales, las ecuaciones de arriba quedan así:

k1 D1 − k2 (D2 − D1 ) − P1 = 0; (1.9a)
k2 (D2 − D1 ) − k3 (D3 − D2 ) − P2 = 0; (1.9b)
k3 (D3 − D2 ) − P3 = 0. (1.9c)

Versión 2014 Preliminar -7-

1.2. El Principio de la Mínima Energía Potencial Total El Método de los Elementos Finitos

Estas ecuaciones se pueden reescribir así:

(k1 + k2 ) D1 − k2 D2 = P1 ; (1.10a)
− k2 D1 + (k2 + k3 ) D2 − k3 D3 = P2 ; (1.10b)
− k3 D2 + k3 D3 = P3 ; (1.10c)

y pueden escribirse de forma matricial:

    
k1 + k2 −k2 0 D1 P1
 −k2 k2 + k3 −k3  D2  = P2  . (1.11)
0 −k3 k3 D3 P3

Ahora, analicemos el problema pero aplicando el análisis por cuerpo libre. Para el primer
«carrito», la ecuación resulta ser:

P1 − k1 D1 − k2 (D1 − D2 ) = 0; (1.12)

como se ve en la figura 1.7.

Figura 1.7: Diagrama de cuerpo libre.

Lo mismo puede hacerse con los otros dos «carritos»:

P2 − k2 (D2 − D1 ) − k3 (D2 − D3 ) = 0 (1.13a)

P3 − k3 (D3 − D2 ) = 0 (1.13b)

Si juntamos las tres ecuaciones y las reordenamos convenientemente, nos queda:

(k1 + k2 ) D1 − k2 D2 = P1 ; (1.14a)
− k2 D1 + (k2 + k3 ) D2 − k3 D3 = P2 ; (1.14b)
− k3 D2 + k3 D3 = P3 ; (1.14c)

; que al ser escrita en forma matricial resulta ser:

(k1 + k2 ) D1 − k2 D2 = P1 ; (1.15a)
− k2 D1 + (k2 + k3 ) D2 − k3 D3 = P2 ; (1.15b)
− k3 D2 + k3 D3 = P3 ; (1.15c)

y pueden escribirse de forma matricial:

    
k1 + k2 −k2 0 D1 P1
 −k2 k2 + k3 −k3  D2  = P2  ; (1.16)
0 −k3 k3 D3 P3

es decir, el mismo sistema al que llegamos planteando el Principio de la Mínima Energía

Potencial Total.

-8- Preliminar Versión 2014

El Método de los Elementos Finitos 1. INTRODUCCIÓN

Vemos que con planteo puramente matemático, obtenemos un sistema de ecuaciones

lineales cuya solución es la configuración de equilibrio. La forma matricial sencilla es K D = P,
donde      
k1 + k2 −k2 0 D1 P1
K=  −k2 k2 + k3 −k3 , D = D2 y P = P2  ;
   
0 −k3 k3 D3 P3
y de la cual podemos establecer que la Energía Potencial Total puede ser escrita como una forma
cuadrática:
1
ΠP = DT K D − DT P. (1.17)
2
Dado que la matriz K es simétrica y definida positiva, al anular la derivada de la Energía
Potencial Total respecto de D y obtener el Deq , la energía potencial total sea mínima, pues la
segunda derivada es mayor que cero.
Hasta aquí lo que hemos hecho es plantear la solución de problemas con grados de libertad
(desplazamientos) discretos. ¿Qué pasa cunado el problema está representado por un sistema en
el cual los grados son «infinitos», o sea, mediante una función continua? Por ejemplo, podríamos
tomar el caso de una viga simplemente apoyada con carga uniformemente distribuida (figura 1.8).

Figura 1.8: Sistema estructural continuo.

Ahora lo que tenemos como incógnita no son desplazamientos sino una función, la elástica
de deformación. En consecuencia, si bien podríamos utilizar el principio visto, ya no podemos
discretizar el sistema como tal. Por lo tanto ahora nuestras incógnitas discretas dejaron de ser
valores para pasar a ser una función: f (x) (o w(x)). Nuestro problema es cómo obtener una
función que haga mínima la energía potencial total de nuestro sistema estructural. En este
sentido, es importante tener en cuenta que cunado nos referimos al principio de la mínima
energía potencial total, nos estamos refiriendo a un valor numérico, un valor que es el menor de
todos los que pueden calcularse para distintas funciones desplazamientos. Es entonces que queda
en evidencia que cada sistema tiene una expresión particular de la energía potencial total, si bien
podemos hallar alguna forma general para expresarla.

Figura 1.9: Sistema sometido a cargas axiles.

Supongamos ahora un sistema sometido a cargas axiales como el de la figura 1.9. Para
este sistema, la expresión de la energía potencial total es:
Z   Z Z
1 T T T T
ΠP = ε Eε − ε Eε0 + ε σ0 dV − u F dV − uT ΦdS − DT P, (1.18)
V 2 V S

Versión 2014 Preliminar -9-

1.2. El Principio de la Mínima Energía Potencial Total El Método de los Elementos Finitos

donde:
Energía Interna de Deformación: V 12 εT Eε − εT Eε0 + εT σ0 dV ;
R  

Fuerzas de volumen: V uT F dV ;
R

Fuerzas de superficie: S uT ΦdS, y;

R

Fuerzas concentradas: DT P .
Para el caso de la energía interna de deformación, analicemos el incremento de la misma,
o sea, d U0 :

d U0 = σx dεx + σy dεy + σz dεz + τxy dγxy + τyz dγyz + τzx dγzx , (1.19)

expresión en la que se despreciado todos los términos diferenciales de orden superior. d U0 es una
derivada total, podemos inferir algo muy importante:
∂U0 ∂U0
=σ o = Eε − Eε0 + σ0 . (1.20)
∂ε ∂ε
Si integramos la expresión de la derecha, obtenemos la expresión de la energía interna de defor-
mación por unidad de volumen:
1 T
ε Eε − εT Eε0 + εT σ0 , (1.21)
2
escrita en forma matricial y vectorial.
A partir de esto, podemos definir para cada sistema estructural la función que representa
a la energía potencial total del mismo y resolverlo mediante el Principio de la Mínima Ener-
gía Potencial Total. Hemos conseguido armar un modelo físico–matemático que nos permite
resolver nuestro problema estructural desde otra visión matemática y que no requiere el análisis
de la ecuación diferencial de equilibrio en forma explícita. sin embargo, a diferencia del caso
discreto, ahora nuestra incógnita es una función y no un valor (o variable), y es algo que está
fuera de lo estudiado en los cursos tradicionales de análisis matemático: obtener una función que
haga mínima una cantidad determinada.
Para este nuevo análisis, supongamos que nos quedamos con la definición de la energía
potencial total de una estructura sometida solamente a fuerzas de volumen. La expresión entonces
se reduce a la siguiente:
Z  
1 T
ΠP = ε Eε − εT Eε0 + εT σ0 − uT F dV. (1.22)
V 2

Tenemos una expresión que nos indica que la energía potencial total del sistema es el
resultado de una integral definida, donde la función principal es el desplazamiento axil de la
estructura. Esto nos lleva a una rama de la matemática que se ocupa justamente de esto: obtener
una cantidad que para una determinada función es mínima. El ejemplo universal para este tipo
de problemas es el siguiente:
Dados dos puntos A y B en el plano, determinar qué función hace mínima la distancia
entre los dos.
En rigor, el problema que dio origen este tipo de análisis no fue la determinación de
la función que hace mínima la distancia entre dos puntos de un plano sino el problema de la
braquistócrona o la determinación de la función o curva que hace mínimo el tiempo que tarda
una partícula para ir desde un punto B, ubicado a una altura h, hasta el punto A, ubicado en
el suelo, por acción únicamente de la gravedad y sin considerar la fricción.
La rama de la matemática que se ocupa de estos problemas es el Cálculo de Variaciones
o Cálculo Variacional, que veremos en el capítulo siguiente.

- 10 - Preliminar Versión 2014

El Método de los Elementos Finitos 2. CÁLCULO VARIACIONAL

CAPÍTULO 2

CÁLCULO VARIACIONAL

2.1. Introducción
¿Qué tienen en común la recta, una línea geodésica, la elástica de deformación, la carga
crítica de Euler, y la ley de Snell (óptica)?
Esta pregunta parece no tener sentido o una respuesta inmediata. Sin embargo, todas
tienen algo en común. Y lo que es más importante, todas surgen de un razonamiento matemático,
sin la necesidad aparente de una interpretación física del problema, aunque cada una de ellas
fueron obtenidas a partir de problemas concretos.
¿Y cuál es la respuesta? Para ello debemos empezar recordando algunos conceptos ya
vistos. Empecemos por uno de los puntos importantes que se ven en los cursos de análisis ma-
temático que es la determinación de los valores máximos y mínimos de una función y(x), para
una x perteneciente a un intervalo (a; b). Si la y(x) tiene derivada continua en ese intervalo, la
condición necesaria para que exista un máximo o un mínimo, es que la derivada de y(x) respecto
de x sea nula para un x0 , perteneciente al intervalo (a; b). Es decir:
dy
=0 para x0 ∈ (a, b).
dx
Por ejemplo, sea la función y(x) = sen x en el intervalo (a, b). Para obtener un máximo
(o un mínimo) debemos hallar el valor de x tal que se cumpla que:
dy d(sen x)
= =0 para x0 ∈ (0, π).
dx dx
Si resolvemos, obtenemos lo siguiente:
d(sen x) π
= cos x = 0 → x0 = .
dx 2
Por lo tanto, cuando x = x0 el valor de y(x0 ) = sen x0 es un máximo (o un mínimo).
Para que y(x) sea un máximo (o un mínimo), se debe verificar que la derivada segunda de y(x)
respecto de x dos veces para x0 sea menor (o mayor) que cero, es decir:

d2 y d2 y
< 0 ⇒ Máximo para x0 ∈ (a, b) > 0 ⇒ Mínimo para x0 ∈ (a, b)
dx2 dx2
En nuestro ejemplo, verifiquemos si y(x0 es un máximo o un mínimo. La segunda derivada
de y(x) en x0 es:
d2 sen x π
2
= − sen x → − sen = −1 < 0,
dx 2

Versión 2014 Preliminar - 11 -

2.2. Cálculo variacional El Método de los Elementos Finitos

Figura 2.1: Máximo para la función y(x) = sen x en x0 = π2 .

es decir, es un máximo, como se ve en la figura .

Con este procedimiento hemos obtenido un valor de x para el cual la función de nuestro
ejemplo y(x) = sen x alcanza un valor máximo en el intervalo dado.
Hemos visto que existe otro tipo de problemas que consisten en encontrar una función
que haga mínimo (o máximo) una cantidad o valor determinado en un intervalo dado. En este
capítulo veremos que existe otra rama del análisis matemático que se ocupa precisamente de lo
que acabamos de mencionar y que permite responder la pregunta inicial. Esta rama es el cálculo
de variaciones o cálculo variacional 1 .
Justamente, la respuesta a la pregunta inicial surge del cálculo variacional, puesto que
todas las funciones allí indicadas son aquellas que minimizan una cantidad o valor determinado.
Por ejemplo, la recta es la función que hace mínima la distancia entre dos puntos de un mismo
plano. Y en la aplicación del método de los elementos finitos haremos uso de la propiedad de
que la solución dela ecuación diferencial de equilibrio, la elástica de deformación es la función
que hace mínima la energía potencial total de un sistema estático, algo que vimos en el capítulo
anterior.

2.2. Cálculo variacional

El cálculo de variaciones es un caso análogo al cálculo diferencial, pero que se refiere
a la determinación de máximos y mínimos de expresiones que incluyen funciones inicialmente
desconocidas. Su importancia es tal que muchos de los principios físico–matemáticos pueden
deducirse a partir de él.
Tomemos inicialmente el caso general para una función de n variables f (x1 , x2 ; . . . ; xn ).
La condición necesaria para que exista un valor máximo o mínimo relativo es que se cumpla lo
siguiente:
∂f ∂f ∂f
df = dx1 + dx2 + · · · + dxn = 0, (2.1)
∂x1 ∂x2 ∂xn
en un punto interior al intervalo dado, y para todos los valores de las diferenciales dx1 ; dx2 ; . . .;
dxn .
1
El cálculo variacional fue desarrollado como respuesta al desafío que realizó Johann Bernouilli en 1696 a sus
colegas de aquel tiempo. El desafío consistió en resolver el problema de la braquistócrona (ver ejemplo 1). La
historia dice que en 1697 Isaac Newton recibió el desafío un día a las cuatro de la tarde y envió la respuesta a las
cuatro de la mañana del día siguiente. Además de la respuesta de Newton, Johann Bernoulli publicó las soluciones
halladas por su hermano Jakob, Leibniz y L’Hôpital, junto con la suya propia. Para más detalles ver [6] y [13].

- 12 - Preliminar Versión 2014

El Método de los Elementos Finitos 2. CÁLCULO VARIACIONAL

Si la función f (x1 , x2 , . . . , xn cumple con (2.1) en un punto x1 , x2 , . . . , xn , se denomina

estacionaria en dicho punto. Si las n variables son todas independientes, es posible asignar
arbitrariamente las n diferenciales y por lo tanto, (2.1) será equivalente a lo siguiente:
∂f ∂f ∂f
df = = = ··· = = 0. (2.2)
∂x1 ∂x2 ∂xn
Supongamos ahora que las n variables no son independientes y que están relacionadas
por m condiciones de la forma:
Φk (x1 , x2 , . . . , xn ) = 0 con k ∈ (1, m) ∧ m < n. (2.3)
En este caso es posible, al menos en forma teórica, resolver las m ecuaciones para expresar
las m variables en términos de las n − m variables restantes, y expresar f y df en términos de
las n − m variables independientes y sus diferenciales.
También es posible hallar las m relaciones lineales entre la n diferenciales y, de esta mane-
ra, expresar a m diferenciales como combinaciones lineales de las n − m variables independientes.

2.2.1. Multiplicadores de Lagrange

Existe, sin embargo, un procedimiento mucho más conveniente para resolver lo expues-
to anteriormente. Consiste en la introducción de los denominados multiplicadores de Lagrange.
Tomemos como ejemplo la función f (x; y; z) y obtengamos los valores estacionarios haciendo
df = fx dx + fy dy + fz dz, (2.4)
que están sujetos a las siguientes restricciones:
Φ1 (x, y, z) = 0 (2.5a)
Φ2 (x, y, z) = 0. (2.5b)
Si aplicamos las relaciones diferenciales:
Φ1x dx + Φ1y dy + Φ1z dz = 0 (2.6a)
Φ2x dx + Φ2y dy + Φ2z dz = 0, (2.6b)
debemos resolver las ecuaciones (2.6) para hallar dx y dy en términos de dz, e introducir esos
resultados en la ecuación (2.4). De esa manera obtenemos una expresión de la diferencial df sólo
en función de dz , es decir, (. . .)dz = 0, y como dz es arbitrario, entonces solo queda que se
cumpla que (. . .) = 0.
El método alternativo es multiplicar las ecuaciones (2.6) por las cantidades λ1 y λ2
respectivamente y sumarlos a (2.4). Así, el esquema queda:
   
∂Φ1 ∂Φ2 ∂Φ1 ∂Φ2
fx + λ1 + λ2 dx + fy + λ1 + λ2 dy+
∂x ∂x ∂y ∂y
  (2.7)
∂Φ1 ∂Φ2
+ fz + λ1 + λ2 dz = 0,
∂z ∂z
donde λ1 y λ2 son arbitrarios.
Para resolver el esquema, obtengamos λ1 y λ2 de manera de que los paréntesis de (2.7)
sean nulos:
∂Φ1 ∂Φ2
fx + λ1 + λ2 =0 (2.8a)
∂x ∂x
∂Φ1 ∂Φ2
fy + λ1 + λ2 =0 (2.8b)
∂y ∂y
∂Φ1 ∂Φ2
fz + λ1 + λ2 = 0. (2.8c)
∂z ∂z

Versión 2014 Preliminar - 13 -

2.2. Cálculo variacional El Método de los Elementos Finitos

Con las ecuaciones (2.5) y (2.8) disponemos de cinco ecuaciones para obtener x, y, z, λ1
y λ2 . Estas dos últimas cantidades son los denominados multiplicadores de Lagrange.
En el cálculo variacional el principal problema consiste en determinar una función tal que
una integral definida que involucre a esa función y a algunas de sus derivadas adopte un valor
máximo o mínimo. Así, lo que se busca es obtener las funciones estacionarias.
En este caso estudiaremos solamente la parte elemental de la teoría, que se refiere a una
condición necesaria que debe ser satisfecha por la función requerida. Es decir, nos ocuparemos
de estudiar las condiciones que llevan a maximizar o minimizar una función determinada, sin
estudiar si se trata de un máximo o un mínimo. La demostración matemática de esto último es
mucho más difícil que en el cálculo diferencial.
Esta condición necesaria suele estar en forma de una ecuación diferencial con condiciones
de contorno. Un ejemplo sencillo de este tipo de problemas es hallar la mínima superficie de
revolución que pasa por dos puntos dados y que gira alrededor del eje x. La función y(x) debe
cumplir que la integral Z x2 p
I =2π y 1 + y 02 dx,
x1

sea mínima, y que y(x1 ) = y1 y y(x2 ) = y2 . (Suponemos además que tanto y1 como y2 son
mayores que cero.)

2.2.2. Cálculo de máximos o mínimos

Analizaremos ahora cómo determinar las funciones que minimizan (o maximizan) una
expresión dada. Consideremos una función continua diferenciable cualquiera y(x) que cumpla
con las siguientes condiciones:

1. Que la integral Z x2
I= F (x, y, y 0 ) dx, (2.9)
x1

tenga un valor mínimo;

2. Que y(x1 ) = y1 y y(x2 ) = y2 .

Elegiremos otra función cualquiera también diferenciable η(x) que se anule en los puntos
extremos del intervalo (x1 , x2 ), es decir, que η(x1 ) = 0 y η(x2 ) = 0. En este caso, la función
y(x) + ε η(x) cumplirá con las condiciones extremas para cualquier constante ε. Por lo tanto la
integral: Z x2
I(ε) = F (x, y + ε η, y 0 + ε η 0 ) dx, (2.10)
x1

dI
tendrá un valor mínimo cuando ε = 0 si = 0, como se ve en la figura 2.2.
dε
Designemos como F1 al integrando F (x, y + ε η, y 0 + ε η 0 ). Notemos que:

dF1 ∂F1 ∂F1 0

= η+ η, (2.11)
dε ∂y ∂y 0

y por lo tanto queda: Z x2  

dI(ε) ∂F1 ∂F1 0
= η+ η dx = 0. (2.12)
dε x1 ∂y ∂y 0
Como F1 → F cuando ε → 0, la expresión (2.12) toma la forma:
Z x2  
dI(ε) ∂F ∂F 0
= η + 0 η dx = 0. (2.13)
dε x1 ∂y ∂y

- 14 - Preliminar Versión 2014

El Método de los Elementos Finitos 2. CÁLCULO VARIACIONAL

Figura 2.2: Funciones y(x), η(x) y y(x) + ε η(x).

Si integramos el segundo término de la derecha por partes (para unificar en η) obtenemos:

Z x2       x2
∂F d ∂F ∂F
η− η dx + η(x) = 0, (2.14)
x1 ∂y dx ∂y 0 ∂y 0 x1

pero como hemos definido que η(x1 ) = 0 y η(x2 ) = 0, entonces la ecuación (2.14) queda como:
Z x2   
∂F d ∂F
− η(x) dx = 0. (2.15)
x1 ∂y dx ∂y 0

Para que la integral sea siempre nula, debe cumplirse que el integrando sea nulo. Como
la función η(x) es arbitraria, la expresión entre corchetes en (2.15) debe anularse en el intervalo
(x1 , x2 ) para que se cumpla la condición dada. Luego, para que la función y(x) haga mínima
(o máxima) la integral de (2.15)(es decir, sea estacionaria), se debe cumplir la ecuación de
Euler–Lagrange:
 
d ∂F ∂F
0
− = 0. (2.16)
dx ∂y ∂y
Las soluciones que se obtienen al resolver de la ecuación de Euler–Lagrange se denominan
extremas del problema considerado. La ecuación (2.16) puede escribirse también de las siguientes
formas:

∂F dy 0 ∂F
     
∂ ∂F ∂ ∂F dy ∂
+ + − = 0; (2.17a)
∂x ∂y 0 ∂y ∂y 0 dx ∂y 0 ∂y 0 dx ∂y

d2 y dy

Fy0 y0 2
+ Fy0 y + Fy0 x − Fy = 0, o; (2.17b)
dx dx

   
1 d ∂F dy ∂F
F− 0 − = 0. (2.17c)
y 0 dx ∂y dx ∂x

∂F
Como generalmente la F no es función explícita de x, entonces ∂x = 0, y la expre-
sión (2.17c) toma una forma particularmente interesante:

∂F dy
F− = C, (2.18)
∂y 0 dx

conocida como Identidad de Beltrami, donde C es una constante de integración.

Versión 2014 Preliminar - 15 -

2.2. Cálculo variacional El Método de los Elementos Finitos

2.2.3. Condiciones naturales de contorno

Analicemos la ecuación (2.14) nuevamente. Si inicialmente no fijamos los valores para
y(x1 ) y y(x2 ), entonces no es necesario que la diferencia entre y(x) y y(x) + ε η(x) se anule
en el intervalo (x1 , x2 ). Sin embargo, el primer término debe anularse si y(x) es la función que
minimiza (o maximiza) nuestra integral, cualquiera sea η(x), y por lo tanto, también lo hará
para el caso particular de todas las η(x) que se anulen en ambos extremos. Valen, entonces,
las consideraciones ya vistas, y el primer término debe ser nulo. Así, sigue siendo válida la
ecuación (2.15) y la ecuación de Euler–Lagrange (ecuación (2.16)).
Podemos deducir, entonces, que el segundo término también debe anularse para cualquier
η(x), es decir:
∂F ∂F
η(x) − η(x) =0 (2.19)
∂y 0 x=x2 ∂y 0
x=x1

Como hemos supuesto que y(x1 ) y y(x2 ) no han sido prefijadas, la funciones η(x1 ) y η(x2 )
son arbitrarias, por lo tanto, para que se cumpla la expresión (2.19), se debe dar necesariamente
que:

∂F
=0 (2.20a)
∂y 0 x=x2

∂F
= 0. (2.20b)
∂y 0 x=x1

Las condiciones que se deben cumplir en (2.20) se denominan condiciones naturales de

contorno, y veremos que son de gran importancia.

2.2.4. La notación variacional

Con el fin de establecer más claramente la analogía entre el cálculo variacional y el cálculo
diferencial, introduciremos la notación variacional.
Tomemos un conjunto S de funciones que satisfagan ciertas condiciones, por ejemplo, que
dichas funciones sean de una sola variable y que su primera derivada sea continua en el intervalo
(a, b). Diremos que cualquier cantidad que adopte un valor numérico específico en correspondencia
con cada función de S es una funcional de dicho conjunto. Así:
Z b Z b
I1 = y(x)dx I2 = y(x) y 00 (x) − [y 0 (x)]2 dx,
a a

son cantidades que llamaremos funcionales, que son las F (x, y, y 0 ) que analizamos antes, porque
para cualquier función y(x) para las cuales están definidas las operaciones indicadas, dichas
cantidades tienen un valor numérico definido (que I1 e I2 , por ejemplo).
Ya habíamos definido a F = F (x; y; y 0 ), que para un valor fijo de x depende de y(x) y de
su derivada y 0 (x). También hemos convertido a la función y(x), que es nuestra incógnita, en la
función y(x) + ε η(x). Este cambio en x se denomina variación de y . Lo simbolizaremos como
δy 2 , es decir:
δy ≡ ε η(x). (2.21)
Por lo tanto, para un valor fijo de x, la funcional F cambia en una cantidad ∆F , la cual
se puede expresar como:

∆F = F [x, y(x) + ε η(x), y 0 (x) + ε η 0 (x)] − F [x, y(x), y 0 (x)]. (2.22)

2
Esta notación fue introducida por Lagrange para remarcar las similitudes y analogías con el cálculo diferencial.

- 16 - Preliminar Versión 2014

El Método de los Elementos Finitos 2. CÁLCULO VARIACIONAL

Si desarrollamos F [x, y(x) + ε η(x), y 0 (x) + ε η 0 (x)] como potencias de ε, tenemos (en
notación simplificada):

∂F ∂F
F (x, y + ε η, y 0 + ε η 0 ) = F (x, y, y 0 ) + ε η + 0 ε η 0 + T (ε), (2.23)
∂y ∂y

donde T (ε) son los términos de orden superior que se desprecian, con lo cual podemos expresar
∆F así:
∂F ∂F
∆F = ε η + 0 ε η0. (2.24)
∂y ∂y
Por analogía con el cálculo diferencial, la ecuación (2.24) se define como la variación de
F, o sea
∂F ∂F
δF = ε η + 0 ε η0. (2.25)
∂y ∂y
Como δy ≡ ε η, y por consiguiente δy 0 ≡ ε η 0 , podemos rescribir la ecuación (2.25) como:

∂F ∂F
δF = δy + 0 δy 0 . (2.26)
∂y ∂y

Ahora, como x no varía, (δx = 0), podemos expresar la (2.26) como:

∂F ∂F ∂F
δF = δx + δy + 0 δy 0 , (2.27)
∂x ∂y ∂y
y la analogía se cumple en forma completa.
Así, mientras la diferencial de una función es una aproximación de primer orden al in-
cremento de esta función según una curva específica, la variación de una funcional es una
aproximación de primer orden al incremento de curva a curva.
Esta analogía se cumple totalmente, pues:

δ(F1 + F2 ) = δF1 + δF2 (2.28a)

δ(F1 F2 ) = δF1 F2 + F1 δF2 (2.28b)

 
F1 δF1 F2 − F1 δF2
δ = . (2.28c)
F2 F22

(Caso particular de (2.28b): δF 2 = 2 F δF .)

Tomemos el caso en el cual x e y son variables independientes, y u y v son funciones de
x y de y, o sea u(x, y) y v(x, y). Además, consideremos la funcional:

F (x, y, u, v, ux , vy , uy , vy ). (2.29)

Si variamos u y v, tendremos las funciones u + ε ξ y v + ε η, siendo las variaciones de u

y v:

δu = ε ξ(x, y), (2.30a)

δv = ε η(x, y). (2.30b)

Al desarrollar según potencias de ε análogamente al caso anterior, hallemos ∆F , de forma

análoga a la ecuación (2.24):

∂F ∂F ∂F ∂F ∂F ∂F
∆F = εξ+ εη+ ε ξx + ε ηx + ε ξy + ε ηy . (2.31)
∂u ∂v ∂u ∂v ∂u ∂v

Versión 2014 Preliminar - 17 -

2.3. Ejemplos El Método de los Elementos Finitos

Ahora bien, dado que x y y son independientes, entonces δx = δy = 0 y, por las defini-
ciones en (2.30a) y (2.30b), podemos escribir la (2.31) así:
∂F ∂F ∂F ∂F ∂F ∂F ∂F ∂F
δF = δx + δy + δu + δv + δux + δvx + δuy + δvy . (2.32)
∂x ∂y ∂u ∂v ∂u ∂v ∂u ∂v
Hay otra propiedad interesante entre el cálculo diferencial y el cálculo variacional. Hemos
visto que δy = ε η (por (2.21)), por lo que entonces δy 0 = ε η 0 , entonces:
 
d dη 0 dy
δy = ε =εη =δ , (2.33)
dx dx dx
d
y como δx = 0, entonces los «operadores» δ y son conmutativos y podemos definir que la
dx
derivada de la variación de una función con respecto a una variable independiente
es igual a la variación de la derivada de esa función con respecto a esa variable
independiente.

2.3. Ejemplos
Para entender más la aplicación del cálculo variacional, nos ocuparemos a continuación de
algunos problemas conocidos. En el primer lugar, aplicaremos el método para demostrar algunos
ejemplos de orden teórico, como el problema que dio origen al cálculo variacional, el problema
de la braquistócrona, y el principio de Hamilton.
Luego plantearemos dos casos de aplicación práctica, que generalmente son resueltos
aplicando el método de los elementos finitos. Se trata de un problema de transmisión del calor
y de un típico problema de cálculo estructural. En estos dos últimos ejemplos nos limitaremos a
plantear el funcional y su variación.

2.3.1. El problema de la braquistócrona

Supongamos dos puntos A y B en un plano e imaginemos que un cable flexible y fino
une ambos puntos Consideremos que A está por encima de B y que una bolita (o sea, una
partícula) se puede desplazar por dicho cable, despreciando la fricción. Modificando la forma del
cable podemos aumentar o disminuir el tiempo que emplea la bolita en ir de A hacia B sólo bajo
acción gravitatoria (peso propio), sin considerar la fricción. El problema de la braquistócrona
consiste en encontrar la forma (curva) de ese cable para que el tiempo que tarda la bolita en ir
desde A hasta B sea mínimo, como se ve en la figura 2.3.

Figura 2.3: El problema de la braquistócrona.

Partamos de lo siguiente. Definamos como y(x) a la curva que une A con B y asumamos
que A = (a, y(a)) y B = (b, y(b)). Para seleccionar la curva nos restringiremos a las funciones

- 18 - Preliminar Versión 2014

El Método de los Elementos Finitos 2. CÁLCULO VARIACIONAL

y(x) tales que pertenezcan a C 2 (a, b), o sea, las funciones continuas cuyas derivadas primera y
segunda también son continuas. Por lo tanto, el tiempo que le insume a la bolita en ir de A a B
está dado por la funcional:
Z b
ds
I= , (2.34)
a v
donde ds es el diferencial de arco de la trayectoria en un instante, y v es la velocidad en ese mismo
instante. Podemos expresar el diferencial ds en función de x e y, mediante la transformación
p
ds = dx2 + dy 2 . (2.35)

Por conservación de la energía tenemos:

1 p
m v 2 = m g[y(x) − y(a)] ⇒ v = 2 g[y(x) − y(a)]. (2.36)
2
Con (2.35) y (2.36) podemos escribir la (2.34) de la siguiente forma:
Z bs s
1 + y 02 1 + y 02
I= dx con F (x, y, y 0 ) = . (2.37)
a 2 g[y(x) − y(a)] 2 g[y(x) − y(a)]

Si asumimos que el mínimo buscado existe, entonces podemos hallar δI. Por simplicidad
adoptaremos A = [0; 0] (a = 0, y(0) = 0), y aplicando la Identidad de Beltrami tenemos:
s
1 + y 02 1 0 2 g y y0
r
0 ∂F
F −y = C → − y = C, (2.38)
∂y 0 2gy 2 1 + y 02 g y

que, simplificada, queda de la siguiente manera:

1
(1 + y 02 )y = = k2 . (2.39)
2 g C2
La solución analítica de esta ecuación diferencial corresponde a una cicloide, cuya expre-
sión en forma paramétrica es:
1
x(θ) = k 2 (θ − sen θ), (2.40a)
2
1
y(θ) = k 2 (1 − cos θ). (2.40b)
2
Como resultado obtenemos que la curva que hace mínimo el tiempo para ir desde A hasta
B bajo acción gravitatoria y sin considerar la fricción es un arco de cicloide.

2.3.2. El principio de Hamilton

Como segundo ejemplo veremos uno de los principios más importantes de la física, el
principio de Hamilton. Lo haremos a través de considerar la dinámica de una partícula. La
trayectoria real de una partícula según las leyes del movimiento de Newton está dada por la
ecuación vectorial:
d r2
m 2 − f = 0. (2.41)
dt
La energía de la partícula al recorrer la trayectoria está dada por la expresión:
Z t2  2 
dr
I= m − f r dt, (2.42)
t1 dt2

que no es otra cosa que una funcional.

Versión 2014 Preliminar - 19 -

2.3. Ejemplos El Método de los Elementos Finitos

Para hallar la trayectoria, planteemos la variación de dicha funcional e igualémosla a cero.

Con esto estamos definiendo que la trayectoria buscada es la que hace mínima la energía de la
partícula. Por lo tanto tendremos:
Z t2  2 Z t2 
d r2
 
dr
δI = m − f δr dt = m δr − f δr dt, (2.43)
t1 dt2 t1 dt2

además que en los extremos se cumple que δrt1 = δrt2 = 0. Si integramos por partes el primer
término de la integral, nos queda:
"   #
dr t2
Z t2 Z t2 Z t2  2
d r2 dr dr 1 dr
m 2 δr dt = m δr − δ dt = −m δ dt. (2.44)
t1 dt dt t1 t1 dt dt t1 2 dt

Es decir, la ecuación (2.43) puede ser escrita como:

Z t2 "  2 #
1 dr
δI = − mδ dt − f δr dt = 0. (2.45)
t1 2 dt

Podemos ver que el primer término de la integral no es otra cosa que la variación de la
energía cinética de la partícula (δEc ), con lo cual, la expresión se puede escribir de la siguiente
forma: Z t2
δI = (δEc + f δr)] dt = 0. (2.46)
t1

Esta expresión es general, pero si el campo de fuerzas aplicado es conservativo, significa

que existe una función Φ tal que:

f δr = X δx + Y δy + Z δz = δΦ, (2.47)

por analogía al cálculo diferencial (recordemos que f dr = X dx + Y dy + Z dz = dΦ). Esta es

la conocida función potencial. Entonce, si reemplazamos (2.47) en (2.46) tenemos que:
Z t2
δI = (δEc + δΦ) dt = 0. (2.48)
t1

En lugar de la función potencial Φ se suele usar la función energía potencial tal que
Φ = −Ep , de modo que finalmente tenemos:
Z t2 Z t2 Z t2 
δI = (δEc − δEp ) dt = δ (Ec − Ep ) dt = δ (Ec − Ep ) dt = 0. (2.49)
t1 t1 t1

El principio de Hamilton enuncia que la integral de la diferencia entre las energías

cinética y potencial es estacionaria para la trayectoria real de una partícula, pudiéndose demostrar
que es un mínimo.

2.3.3. Transmisión del calor

Una placa de longitud infinita, como se muestra en la figura 2.4, tiene una temperatura
inicial θt , cuando la superficie en x = 0 repentinamente recibe un flujo de calor uniforme. La
otra cara de la placa (x = L) mantiene la temperatura θi , y las superficies paralelas al plano xz
están aisladas.
Consideremos un diferencial de la placa. El equilibrio calórico está dado por:
 
∂q ∂θ
qA−A q+ dx = ρ A dx c , (2.50)
∂x ∂t

- 20 - Preliminar Versión 2014

El Método de los Elementos Finitos 2. CÁLCULO VARIACIONAL

Figura 2.4: Transmisión del calor.

y según la ley de Fourier de transmisión del calor tenemos que:

∂θ
q = −k . (2.51)
∂x
Si reemplazando la (2.51) en la (2.50) nos queda lo siguiente:

∂2θ ∂θ
k 2
= ρc , (2.52)
∂x ∂t
que es la ecuación diferencial en una dimensión de la transmisión del calor. Las condiciones de
borde para este problema en particular son:

∂θ(0, t) q0 (t)
=− ; (2.53a)
∂x k

t > 0; (2.53b)

θ(L, t) = θi . (2.53c)

Ahora planteemos la funcional que gobierna este problema:

Z  2 Z
1 ∂θ
ΠT = k dx − θ q B dx − θ0 q0 . (2.54)
L 2 ∂x L

Al igual que en los casos anteriores, variemos la funcional para obtener un máximo o
mínimo: Z Z
∂θ ∂θ
δΠT = k δ dx − δθ q B dx − δθ0 q0 (2.55)
L ∂x ∂x L
Como hemos visto antes, la variación y la diferenciación son conmutables, y si además
integramos por parte la primera integral, podemos escribir la (2.55) de esta forma:
Z  2   
∂ θ B ∂θ ∂θ
δΠT = k 2 +q δθdx + k δθL − k + q0 δθ0 = 0. (2.56)
L ∂x ∂x x=L ∂x x=0

Nuevamente, para que la funcional se anule, deben ser nulos todos los términos. Así, se
debe cumplir que
∂2θ
k 2 + q B = 0; (2.57)
∂x

Versión 2014 Preliminar - 21 -

2.3. Ejemplos El Método de los Elementos Finitos

dado que δθ0 es arbitrario, se debe cumplir que

∂θ
k + q0 = 0, (2.58)
∂x x=0

que es la condición de borde vista en (2.53a). Como

∂θ
q B = −ρ c , (2.59)
∂t
si reemplazamos esta última en la ecuación (2.56), nos queda

∂2θ ∂θ
k 2
= ρc , (2.60)
∂x ∂t
que es la ecuación diferencial hallada antes (la ecuación (2.52)).

2.3.4. Cálculo estructural

Ahora nos ocuparemos del caso de una viga empotrada–apoyada, aunque este apoyo en
realidad es un resorte cuya constante es k, y sometida a una carga axil H como se ve en la
figura 2.5.

Figura 2.5: Viga empotrada–apoyada.

La funcional del problema es:

 2 2
1 dw 2
Z Z  
1 d w 1 2
ΠP = EI 2
dx − H dx + kR wL , (2.61)
L 2 dx L 2 dx 2

donde wL = w(x = L) y las condiciones esenciales de borde son:

w(x = 0) = 0, (2.62a)

dw
= 0. (2.62b)
dx x=0

Para hallar la condición estacionaria debe ser δΠP = 0, es decir:

d2 w
Z  2  Z  
d w dw dw
δΠP = EI δ dx − H δ dx + kR wL δwL = 0. (2.63)
L dx2 dx2 L dx dx
Integremos por parte los dos primeros términos y así obtendremos una nueva expresión:

d4 w d2 w d2 w dw d2 w dw
Z  
E I 4 − H 2 δw dx + E I 2 δ − EI 2 δ +
L dx dx dx dx x=L dx dx x=0
(2.64)
d3 w d3 w
   
dw dw
− EI 3 +H δw + EI 3 +H δw + k wL δwL = 0.
dx dx x=L dx dx x=0

- 22 - Preliminar Versión 2014

El Método de los Elementos Finitos 2. CÁLCULO VARIACIONAL

Analicemos los términos de la misma. Como δw debe cumplir con las condiciones de
2
borde, tenemos que δw|x=0 = δw0 |x=0 = 0 y entonces son nulos los términos E I ddxw2 δ dw
dx x=0 y
 3

E I ddxw3 + H dw
dx δw . Entonces la ecuación (2.64) queda así:
x=0

d4 w d2 w d2 w dw
Z  
E I 4 − H 2 δw dx + E I 2 δ
L dx dx dx dx x=L
 3
 (2.65)
d w dw
− EI 3 +H δw + k wL δwL = 0.
dx dx x=L

Como δw es arbitraria, para que la ecuación (2.65) se anule, debe cumplirse que:

d4 w d2 w
EI − H =0 Ecuación de equilibrio, (2.66a)
dx4 dx2

d2 w
EI =0 Momento flector en x = L y; (2.66b)
dx2 x=L

d3 w dw
EI 3
+H − k wL = 0 Corte en x = L. (2.66c)
dx dx x=L

Podemos notar que las dos últimas expresiones corresponden a las condiciones de borde
del problema.

Versión 2014 Preliminar - 23 -

El Método de los Elementos Finitos 3. TEOREMAS ENERGÉTICOS

CAPÍTULO 3

TEOREMAS ENERGÉTICOS

3.1. Definiciones
3.1.1. Tipos de materiales
Empezaremos con una serie de definiciones sobre los distintos tipos de materiales que
existen.

Material con efecto viscoso

Cuando en un material se cumple que la relación σ–ε es dependiente del tiempo, es decir,
el material sigue un diagrama σ–ε–t, entonces se dice que el material posee efecto viscoso.
Si en cambio la carga es independiente del tiempo, entonces se habla de un material sin
efecto viscoso, y la representación del diagrama σ–ε–t se reduce al diagrama σ–ε, bidimensional

Material con efecto plástico

Cuando el proceso de descarga en un material no coincide con el proceso de carga, se dice
que el material posee efecto plástico. En particular, si la curva de descarga es una recta paralela
al eje σ, entonces se trata de un material plástico.

Material elástico
Cuando ambos efectos, el plástico y el viscoso, son nulos se dice que el material es elástico.

3.1.2. Trabajo de las fuerzas exteriores y energía de deformación

También definiremos las expresiones de trabajo de las fuerzas exteriores, trabajo comple-
mentario de las fuerzas exteriores, energía de deformación y energía complementaria de defor-
mación.

Trabajo de las fuerzas exteriores

El trabajo de las fuerzas exteriores esta dado por la expresión:
n
"Z Z Z Z #
X
W = Ri dui dV + Ti dui dS . (3.1)
i=1 V u∗i S u∗i

Versión 2014 Preliminar - 25 -

3.1. Definiciones El Método de los Elementos Finitos

Trabajo complementario de las fuerzas exteriores

El trabajo complementario de las fuerzas exteriores esta dado por la siguiente expresión:
n
"Z Z Z Z #
X
W = ui dRi dV + ui dTi dS . (3.2)
i=1 V Ri∗ S Ti∗

Energía de deformación
La energía de deformación está dada por la expresión:
Xn Z Z
U= σij dεij dV . (3.3)
i=1 V ε∗ij

Energía complementaria de deformación

La energía de deformación está dada por la expresión:
Xn Z Z
U= εij dσij dV . (3.4)
V ∗
σij
i=1

3.1.3. Expresiones para solicitación axil

Analicemos cada una de estas definiciones para el caso de una barra de sección transversal
A y longitud L solicitada por una carga axil P en el extremo x = L, como se observa en la
figura 3.1. Para este ejemplo, cada una de las definiciones anteriores se reduce a:

Figura 3.1: Barra solicitada axilmente.

1. Trabajo de las fuerzas exteriores:

Z
W = P dui . (3.5)
u∗i

2. Trabajo complementario de las fuerzas exteriores:

Z
Wc = u∗i dP. (3.6)
P∗

- 26 - Preliminar Versión 2014

El Método de los Elementos Finitos 3. TEOREMAS ENERGÉTICOS

Z
U =A·L σ dε. (3.7)
ε∗

3. Energía complementaria de deformación:

Z
Uc = A · L ε dσ. (3.8)
σ∗

Las figuras 3.2 y 3.3 muestran las representaciones gráficas de las expresiones anteriores.

Figura 3.2: Trabajo de las fuerzas exteriores y

trabajo complementario de las fuerzas exteriores.

Figura 3.3: Energía interna de deformación y

energía interna complementaria de deformación.

Notemos que las expresiones (3.5) a (3.8) son igualmente válidas para el caso de un
material viscoso y plástico. En todas estas expresiones hemos asumido las siguientes expresiones:

P uL
σ= ; ε= . (3.9)
A L

3.2. Principio de equivalencia

El Principio de Equivalencia estipula que en un proceso de carga arbitrario, el trabajo
W desarrollado por las fuerzas exteriores (3.1) es igual a la energía de deformación U (3.3). Es
decir este principio se expresa como:
W = U. (3.10)
A partir de este principio y teniendo en cuenta el equilibrio, se puede deducir inmediatamente
que:
W c = Uc , (3.11)

Versión 2014 Preliminar - 27 -

3.2. Principio de equivalencia El Método de los Elementos Finitos

es decir, que el trabajo complementario de las fuerzas exteriores es igual a la energía comple-
mentaria de deformación. La validez de este principio no está relacionada con las propiedades
del material.
Vamos a relacionar ahora el principio de equivalencia con el principio de conservación
de la energía. Para que resulte más fácil de entender, nos vamos a valer de la barra solicitada
axilmente de la figura 3.1. Fijémonos primero en la figura 3.4.

Figura 3.4: Energías.

La curva OMA representa el proceso de carga de la barra. El área debajo de esta curva,
como vimos antes, es la energía de deformación U dividida por el volumen de la barra (A · L).
Supongamos ahora que el proceso de descarga está representado por la curva ANB. Como
podemos observar, la barra no ha retornado a su longitud inicial, sino que tiene una deformación
residual representada por el segmento OB. Podemos suponer, entonces, que el área bajo la curva
ANB es la energía potencial elástica, pues es la energía restitutiva o recuperable una vez que la
barra ha sido descargada.
Planteemos el Primer Principio de la Termodinámica o principio de la conservación de la
energía. La energía potencial clásica, para una transformación isoterma, está dada por:

U = W − Q, (3.12)

donde W es el trabajo realizado por las fuerzas exteriores, y Q es la cantidad de calor transferida
por la barra al medio ambiente. Por lo tanto, de las ecuaciones (3.10) y (3.12) podemos deducir
que
Q = U −U. (3.13)

Así,la cantidad de calor disipada está dada por el área delimitada por las curvas OMA
y ANB, de la figura 3.4, multiplicada por el volumen de la barra.
Si analizamos un poco más las definiciones vistas en el punto anterior podemos arribar a
las siguientes conclusiones:

1. Para un material elástico, la cantidad de calor disipada en los procesos de carga y descarga
es nula y por consiguiente las energías de deformación y potencial son iguales. Se trata, en
consecuencia, de un sistema conservativo.

2. Para un material plástico, la cantidad de calor disipada es igual a la energía de deformación,

siendo la energía potencial nula. En este caso, estamos ante un sistema disipativo.

3. Para un material viscoso pero sin efecto plástico, la cantidad de calor disipada estará dada
por el área encerrada por las curvas OMA y ANB y el segmento BO, multiplicada por
el volumen de la barra.

- 28 - Preliminar Versión 2014

El Método de los Elementos Finitos 3. TEOREMAS ENERGÉTICOS

3.3. Teorema de los desplazamientos virtuales

Tomemos nuevamente una barra de sección transversal constante (A = cte) solicitada
axilmente con una carga distribuida axil t(x) y una carga concentrada P en el extremo de la
barra.
Supongamos además, que existe una función desplazamiento u(x) virtual, que no conoce-
mos, que cumple con las condiciones de borde cinemáticas, es decir u(0) = 0, tal como se observa
en la figura 3.5. Como la función u(x) es desconocida, le aplicamos una variación de u (o sea δu)
para poder ajustar la curva, que también debe cumplir con las condiciones de borde, y entonces,
δu(0) = 0.

Figura 3.5: Barra solicitada axilmente.

De acuerdo con el principio de equivalencia tenemos que U = W 1 . En nuestro caso queda:

U + δU = W + δW. (3.14)

Una representación gráfica de la ecuación (3.14) se ve en la figura 3.6.

Figura 3.6: Energía interna de deformación y trabajo de las fuerza exteriores.

1
No todos los autores están de acuerdo con esta forma de expresar el principio de equivalencia. La objeción
reside en la caracterización de W como el área bajo la curva de carga. Hoy la mayoría prefiere definir el trabajo de
las fuerzas exteriores como W = P · u, y plantear el principio de la mínima energía potencial total por aplicación
del cálculo variacional, como se vio en su momento en el capítulo 2.

Versión 2014 Preliminar - 29 -

3.3. Teorema de los desplazamientos virtuales El Método de los Elementos Finitos

La ecuación (3.14) se puede escribir de la siguiente forma:

δu U = δu W ⇒ δu (U − W ) = 0. (3.15)

De esta manera, lo que nos queda planteado en la (3.15) es la variación de la funcional,

que en nuestro caso minimiza la energía (pero podría maximizarla) que desarrolla el sistema y
por lo tanto, la función u(x) que tratamos de hallar debe cumplir con esa condición.
Si definimos δu U como Z
δu U = σ δε dV, (3.16)
V
dado que en la barra A = cte y que, por hipótesis, adoptamos que se cumple la relación cinemática
(ε = du
dx ), podemos escribir la (3.16) como:
Z  
du
δu U = A σδ dx, (3.17)
V dx
Como hemos visto, la variación y la diferenciación son conmutativas, de modo que la (3.17)
podemos escribirla de esta forma:
Z Z
d
δu U = A σ δu dx = A σ dδu. (3.18)
V dx V

Ahora, integremos la (3.18) por partes. Nos queda lo siguiente:

Z Z
L dσ
δu U = A [σ δu(x)]0 − A δu dσ = A σL δuL − A δu dx, (3.19)
L L dx
pues por hipótesis hemos establecido que δu(0) = 0.
Si definimos la variación del trabajo de las fuerzas exteriores como
Z
δu W = t(x) δu dx + P δuL , (3.20)
L

podemos reemplazar en la (3.15) las ecuaciones (3.19) y (3.20). La funcional nos queda de esta
forma: Z Z
dσ
A σL δuL − A udxdx − t(x) δu dx − P δuL = 0. (3.21)
L δ L
Al reagrupar los términos, la ecuación nos queda de esta forma:
Z  
dσ
− A + t(x) δu dx + (A σL − P ) δuL = 0. (3.22)
L dx
Como hemos visto, para que la funcional se anule para un δu(x) arbitrario, se debe
cumplir que:
dσ
A + t(x) = 0 Ecuación diferencial de equilibrio, y; (3.23a)
dx
A σL = P Condición de borde estática. (3.23b)

que son las condiciones necesarias para que el sistema esté en equilibrio.
A partir de este resultado, podemos enunciar el Teorema de los Desplazamientos
Virtuales de la siguiente manera:
Si a un sistema sometido a un estado de carga, se le aplica un campo de despla-
zamientos virtuales que sea compatible con los vínculos, y que cumpla con la relación
cinemática, es decir, que δε = ddx
δu
; si se cumple que la variación de la energía interna
de deformación es igual a la variación del trabajo de las fuerzas exteriores, entonces
se puede asegurar que se cumple el equilibrio y las condiciones de borde estáticas.

- 30 - Preliminar Versión 2014

El Método de los Elementos Finitos 3. TEOREMAS ENERGÉTICOS

3.4. Teorema de las tensiones virtuales

En este caso vamos a tomar la misma barra de la figura 3.5, pero en vez de estar solicitada
por un estado de carga, está afectada por un campo de desplazamientos. Partamos de los teoremas
complementarios y apliquemos a esta barra un sistema de fuerzas virtuales en equilibrio que
cumpla con las condiciones de borde estáticas (figura 3.7).

Figura 3.7: Barra con desplazamientos axiles.

Si definamos la energía complementaria interna como:

Z Z
Uc = ε∗ σ ∗ (x) dV = A ε∗ σ ∗ (x) dx, (3.24)
V L

entonces la variación de la energía complementaria es:

Z Z
δσ Uc = ε δσ (x) dV = A ε∗ δσ ∗ (x) dx,
∗ ∗
(3.25)
V L

Hagamos lo mismo pero con el trabajo complementario de las fuerza exteriores:

Z
Wc = u∗ (x) t(x)dx + u∗L P. (3.26)
L

Al igual que para la energía complementaria, la variación del trabajo complementario es:
Z
δ σ Wc = u∗ (x) δt(x)dx + u∗L δP. (3.27)
L

En la figura 3.8 se pueden las representaciones gráficas de las ecuaciones (3.24), (3.25), (3.26)
y (3.27).
En forma análoga al teorema anterior, armemos la funcional δσ (Uc − Wc ) = 0:
Z Z
∗ ∗
A ε δσ (x) dx − u∗ (x) δt(x)dx − u∗L δP = 0. (3.28)
L L

Como dijimos al principio, el sistema de fuerzas virtuales está en equilibrio; deben cumplir
con las condiciones de borde estáticas. En consecuencia, podemos expresar la (3.28) así:
d δσ ∗ (x)
Z Z   Z
∗ ∗ ∗
A ε δσ (x) dx + A + δt(x) u (x) dx − u∗ (x) δt(x)dx − u∗L δP = 0, (3.29)
L L dx L

Versión 2014 Preliminar - 31 -

3.4. Teorema de las tensiones virtuales El Método de los Elementos Finitos

Figura 3.8: Energía interna de deformación y trabajo de las fuerzas exteriores complementarios.

que podemos reescribir de esta manera:

d δσ ∗ (x) ∗
Z Z Z Z
A ε∗ δσ ∗ (x) dx + A u (x) dx + u∗ (x) δt(x)dx − u∗ (x) δt(x)dx − u∗L δP = 0,
L L dx L L
(3.30)
y reducir a:
d δσ ∗ (x) ∗
Z Z
A ε∗ δσ ∗ (x) dx + A u (x) dx − u∗L δP = 0. (3.31)
L L dx
∗
Ahora, integremos por partes L A d δσdx(x) u∗ (x) dx:
R

d δσ ∗ (x) ∗ d u∗ (x) ∗
Z Z
A u (x) dx = [A δσ ∗ (x) u∗ (x)]L
0 −A δσ (x) dx. (3.32)
L dx L dx

Si reemplazamos la (3.32) en la (3.31), nos queda:

d u∗ (x) ∗
Z Z
∗ ∗ ∗ ∗ L
A ε δσ (x) dx + [A δσ (x) u (x)]0 − A δσ (x) dx − u∗L δP = 0. (3.33)
L L dx

Si desarrollamos el términos entre corchetes, nos queda:

d u∗ (x) ∗
Z Z
∗ ∗ ∗ ∗ ∗ ∗
A ε δσ (x) dx + A δσL uL − A δσ0 u0 − A δσ (x) dx − u∗L δP = 0, (3.34)
L L dx

y que podemos agrupar de esta forma:

d u∗ (x)
Z  
∗
A ε − δσ ∗ (x) dx + [A δσL∗ − δP ] u∗L − A δσ0∗ u∗0 = 0, (3.35)
L dx

y reducir a
d u∗ (x)
Z  
∗
A ε − δσ ∗ (x) dx − A δσ0∗ u∗0 = 0, (3.36)
L dx
pues por definición A δσL∗ − δP = 0.
Nuevamente, para que la funcional sea mínima, se debe cumplir que:

d u∗ (x)
ε∗ − =0 Relación cinemática, y; (3.37a)
dx
u∗0 = 0 Condición de borde cinemática, (3.37b)

dado que δσ ∗ (x) y δσ0∗ son arbitrarios.

Nuevamente, con lo que hemos obtenido, podemos enunciar el Teorema de los Des-
plazamientos Virtuales Complementario o Teorema de las Tensiones Virtuales:

- 32 - Preliminar Versión 2014

El Método de los Elementos Finitos 3. TEOREMAS ENERGÉTICOS

Si a un sistema que se ha sometido a un campo un desplazamiento y a un campo de

deformaciones, le aplicamos un sistema de fuerzas virtual en equilibrio con un campo
de tensiones, también virtual, que cumpla con las condiciones de borde estáticas, y
se cumple que variación de la energía interna de deformación es igual a la variación
del trabajo de las fuerzas exteriores, entonces podemos asegurar que se cumplen la
relación cinemática y las condiciones de borde cinemáticas.

3.5. El principio de la mínima energía potencial total

Este principio lo hemos tratado al principio como una de las bases del método de los
elementos finitos, asociado al cálculo variacional. Ahora nos ocuparemos de ver que relación tiene
con los dos teoremas anteriores. Para ello, tomemos nuevamente nuestra barra de la figura 3.5,
y definamos la energía potencial total de dicha barra:

du 2
Z   Z
1
ΠP = EA dx − t(x) u(x) dx − P uL . (3.38)
L 2 dx L

Como vimos, para hallar la función u(x) que haga mínima la funcional (ecuación (3.38)),
debemos variar la misma e igualarla a cero. Eso nos lleva a lo siguiente:

du 2
Z   Z
1
δΠP = EAδ dx − t(x) δu(x) dx − P δuL = 0, (3.39)
L 2 dx L

que podemos escribir como

Z   Z
du du
δΠP = EA δ dx − t(x) δu(x) dx − P δuL = 0. (3.40)
L dx dx L

Como suponemos que se cumple la relación cinemática, entonces podemos reemplazar

du(x)
dx por ε. Así, la ecuación (3.40) queda:
Z Z
δΠP = E ε δε Adx − t(x) δu(x) dx − P δuL = 0. (3.41)
L L

Por la ley de Hooke, σ = E ε. Si reemplazamos en la anterior, nos queda:

Z Z
δΠP = σ δε A dx − t(x) δu(x) dx − P δuL = 0, (3.42)
L L

pero como también se cumple que dV = A dx, entonces tenemos que

Z Z
δΠP = σ δε dV − t(x) δu(x) dx − P δuL = 0, (3.43)
V L

es decir, tenemos la siguiente ecuación:

Z Z
σ δε dV = t(x) δu(x) dx + P δuL , (3.44)
V L

que no es otra cosa que el Teorema de los Desplazamientos Virtuales. Queda en evidencia,
pues, que la aplicación del Principio de la Mínima Energía Potencial Total lleva a la
aplicación de alguno de los teoremas vistos –si uno es válido, el otro lo es también.
Veremos más adelante que este principio es importante porque nos permitirá determinar
cuales de las soluciones aproximadas obtenidas mediante el método de los elementos finitos es la
mejor aproximación a la solución del problema, y consecuentemente, de la calidad de la malla de
elementos usada.

Versión 2014 Preliminar - 33 -

3.6. Ejemplos de aplicación El Método de los Elementos Finitos

3.6. Ejemplos de aplicación

3.6.1. Barra solicitada axilmente
Veamos un caso práctico de cómo aplicar el Teorema de los Desplazamientos Vir-
tuales. Tomemos una barra de sección constante solicitada por una carga distribuida axil, como
en la figura 3.9:

Figura 3.9: Barra empotrada–libre con carga distribuida axil.

La ecuación diferencial del problema es:

d2
A·E· u(x) + t(x) = 0. (3.45)
dx2
Para resolver la ecuación diferencial, basta con integrarla dos veces:
Z
d u(x) t(x) t
=− dx + C1 = − x + C1 ; (3.46a)
dx AE AE

Z  
t 1 t
u(x) = − x + C1 dx + C2 = − x2 + C1 x + C2 . (3.46b)
AE 2 AE

Para obtener la ecuación de los desplazamientos, debemos tomar en cuenta las condiciones
de borde. Éstas son:

1. El desplazamiento en el extremo empotrado es nulo: u(0) = 0, y;

2. El esfuerzo normal en el extremo libre es nulo: N (x = L) = 0, que también se puede escribir

como A σL = 0 y como A d u(x)dx = 0.
L

Con la primera condición de borde obtenemos que C2 = 0, en tanto que con la segunda
condición de borde obtenemos que C1 = At LE . Con estos dos nuevos valores tenemos finalmente
la ecuación de los desplazamientos:

1 t tL tx  x
u(x) = − x2 + x= L− . (3.47)
2 AE AE AE 2
Calculemos el desplazamiento en el extremo libre de la barra:

t L2
 
tL L
u(L) = L− = .
AE 2 2AE

Si queremos hallar la ecuación del esfuerzo normal, nos basta con derivar la función u(x)
respecto a x, y multiplicarla por A E. Así obtenemos:

N (x) = t (L − x). (3.48)

- 34 - Preliminar Versión 2014

El Método de los Elementos Finitos 3. TEOREMAS ENERGÉTICOS

Pero lo que en realidad queremos es aplicar el Teorema de los Desplazamientos Virtuales.

Para ello, propongamos una función desplazamiento que cumpla con las condiciones de borde
cinemáticas:
x
u1 (x) = uL , (3.49)
L
donde uL es el desplazamiento en el extremo libre de la barra, que tomaremos como incógnita.
Propongamos una variación de los desplazamientos con las mismas características:
x
δu1 (x) = δuL , (3.50)
L
donde δuL es análogo a uL .
Notemos que hemos propuesto una función lineal en ambos casos, aún cuando sabemos
que la solución de la ecuación diferencial es una función cuadrática. Pero de acuerdo con el
teorema, la función desplazamiento y su variación, sólo deben cumplir con las condiciones de
borde cinemáticas, algo que la función lineal propuesta cumple sin problemas.
Para ambas funciones tendremos los siguiente:
d u1 (x) uL
= , (3.51)
dx L

d δu1 (x) δuL

= . (3.52)
dx L
Ahora, planteemos el teorema con las funciones (3.50), (3.51) y (3.52):
Z Z
uL δuL δuL
δU = EA dx = t(x) x dx = δW. (3.53)
L L L L L
Como todos los términos son constantes excepto la función δu(x), la ecuación (3.53)
queda así: Z Z
uL δuL δuL
EA dx = t x dx, (3.54)
L L L L L
que al integrar y simplificar queda
uL L
EA =t , (3.55)
L 2
y de donde despejamos uL :
t L2
uL = . (3.56)
2E A
Resulta al aplicar el Teorema de los Desplazamientos Virtuales obtuvimos el mismo valor
para el desplazamiento en L que al resolver la ecuación diferencial, pero con una función u(x)
lineal, en tanto que el resultado obtenido al resolver la ecuación diferencial es una función cua-
drática. ¿Cómo se explica esta diferencia? Antes de ensayar una explicación, armemos la función
desplazamiento y a partir de ésta, la de esfuerzo normal:
tL
u1 (x) = x, (3.57)
2E A

d uD (x) tL tL
N1 (x) = E A =EA = . (3.58)
dx 2E A 2
Vemos que la función del esfuerzo normal es una constante, en cambio, la solución hallada
resolviendo la ecuación diferencial nos dio una función lineal. Si calculamos el valor de N para
x = L2 con la ecuación (3.48), obtenemos lo siguiente:
   
L L tL
N =t L− = , (3.59)
2 2 2

Versión 2014 Preliminar - 35 -

3.6. Ejemplos de aplicación El Método de los Elementos Finitos

es decir, que el esfuerzo normal obtenido al aplicar el teorema de los desplazamientos virtuales
corresponde al esfuerzo normal en la sección media de la barra.
Veamos qué pasa si cambiamos la función δuD (x) por esta otra:
uL 2
δu2 (x) = x . (3.60)
L2
Esta variación también cumple con las condiciones de borde cinemáticas. Entonces cal-
culemos su derivada:
d u2 (x) uL
= 2 2 x. (3.61)
dx L
Ahora reemplacemos (3.60) y la (3.61) en la (3.53):
Z Z
uL δuL δuL
δU = EA 2
2 x dx = t(x) 2 x2 dx = δW. (3.62)
L L L L L
Si integramos nos queda
uL δuL L2 δuL L3
EA 2 = t(x) . (3.63)
L L2 2 L2 3
Si simplificamos δuL y L, nos queda un nuevo valor de uL :
t L2
uL = , (3.64)
3E A
y una nueva función desplazamiento:
tL
u2 (x) = x (3.65)
3E A
Lo primero que notamos es que al cambiar la variación de la función, cambió el valor
de uL . Por lo tanto, para obtener una aproximación de la función desplazamiento no sólo es
importante elegir la función sino también la variación de la función.
Con esta nueva función podemos obtener también la función del esfuerzo normal
tL
N2 (x) = , (3.66)
3
que también es una constante pero diferente a la hallada con u1 (x).
Los resultados al aplicar el teorema son algo desconcertantes porque hemos obtenido dos
soluciones aproximadas. Nos resta saber cuál de las dos es la mejor. Es aquí donde aplicamos
el principio de la mínima energía potencial total. Para cada una de las soluciones tenemos lo
siguiente:

tL 2
Z   Z
1 tL
Π1P = EA dx − t x dx (3.67a)
L 2 2E A L 2E A

Z  2 Z
1 tL tL
Π2P = EA dx − t x dx. (3.67b)
L 2 3E A L 3E A
Al integrar las dos nos queda
t2 L3 t2 L3
 
1
Π1P =− = − , (3.68a)
8E A EA 8

5 t2 L3 t 2 L3
 
5
Π2P =− = − . (3.68b)
54 E A EA 54

- 36 - Preliminar Versión 2014

El Método de los Elementos Finitos 3. TEOREMAS ENERGÉTICOS

Como − 18 es menor que − 54

5
, la mejor aproximación es la función u1 (x).
Nos queda explicar el porqué de las diferencias entre una y otra solución aproximada.
Para ello analicemos qué ocurre si proponemos como función de desplazamiento a la siguiente
función:
x x 
u(x) = uL 1− , (3.69)
L 2L
similar a la función que obtuvimos al resolver la ecuación diferencial.
Tomemos como variación de u(x) la primera, es decir, δu1 (x) y planteemos el teorema de
los desplazamientos virtuales:
Z Z
uL  x  δuL δuL
δU = EA 1− dx = t(x) x dx = δW. (3.70)
L L L L L L
Rearmemos la expresión para facilitar los cálculos:
Z  Z
uL x
EA 1− dx = t x dx. (3.71)
L L L L

Al integrar ambas expresiones nos queda:

uL x2 x2
EA x− =t , (3.72)
L 2L L 2 L

que termina siendo

uL L2
EA =t , (3.73)
2 2
del que obtenemos
t L2
uL = , (3.74)
EA
y que nos da la nueva función u(x):
tLx  x 
u(x) = 1− , (3.75)
EA 2L
que resulta ser la misma función que obtuvimos resolviendo la ecuación diferencial.
Ahora veamos qué pasa si tomamos la otra variación, δu2 (x), y planteamos nuevamente
el teorema de los desplazamientos virtuales. En este caso tendremos los siguiente:
Z Z
uL  x  δuL δuL
δU = EA 1− 2
2 x dx = t(x) 2 x2 dx = δW. (3.76)
L L L L L L
Nuevamente, reduzcamos la ecuación para que sea fácil operar con ella. Entonces nos
queda: Z  Z
uL δuL x δuL
EA 1− 2 x dx = t 2 x2 dx. (3.77)
L L2 L L L L
Al igual que el caso anterior, simplifiquemos la expresión e integremos. Así obtendremos
que
uL 2 2 x3 x3
EA x − =t , (3.78)
L 3L L 3 L
que se reduce a:
uL L L3
EA =t . (3.79)
3 3
Finalmente, si despejamos uL obtenemos que

L2
uL = t , (3.80)
EA

Versión 2014 Preliminar - 37 -

3.6. Ejemplos de aplicación El Método de los Elementos Finitos

que es el mismo resultado que para el caso anterior. Por lo tanto, volvemos a obtener la misma
ecuación que obtuvimos al variar con δu1 (x);

tLx  x 
u(x) = 1− , (3.81)
EA 2L
y que es también la que resuelve la ecuación diferencial.
Verifiquemos que ocurre con la energía potencial total del sistema con esta función. Al
plantearla nos queda:

(t L)2 
Z Z
1 x 2 tLx  x 
ΠP = EA 1 − dx − t(x) 1 − dx. (3.82)
L 2 (E A)2 L L EA 2L

Al integrar todos los términos nos queda:

1 t2 L3 1 t2 L3 1 t2 L3
ΠP = − =− , (3.83)
6 EA 3 EA 6 EA
resultado que es menor a las otras dos energía potenciales totales que hallamos con las soluciones
aproximadas. Al ser la menor de toda, por definición, es la solución del problema.
Podemos decir, entonces, que al plantear como función desplazamiento aproximada una
que esté incluida en la familia de funciones de la solución de la ecuación diferencial, obtendremos,
al aplicar el Teorema de los Desplazamientos Virtuales, la misma función que es solución de la
ecuación diferencial, cualquier fuere la variación que apliquemos.

3.6.2. Viga solicitada a flexión

En este caso, tomemos una viga simplemente apoyada solicitada por una carga distribuida
uniforme, como en la figura .

Figura 3.10: Viga sometida a flexión.

Vimos en capítulos anteriores la formulación de la energía potencial total para el caso de

flexión, o sea, la funcional del problema. Para la estructura en análisis, la funcional es:
2
d2 w(x)
Z  Z
1
ΠF = EI dx − q(x) w(x) dx. (3.84)
L 2 dx2 L

La función desplazamiento que resuelve el problema es aquella que haga mínima la fun-
cional. Por eso busquemos eso variando a ésta e igualando a cero:
2
d2 w(x)
Z  Z
1
δΠF = EI δ dx − q(x) δw(x) dx = 0, (3.85)
L 2 dx2 L

que podemos escribir como

d2 w(x) d2 δw(x)
Z Z
EI dx = q(x) δw(x) dx, (3.86)
L dx2 dx2 L

- 38 - Preliminar Versión 2014

El Método de los Elementos Finitos 3. TEOREMAS ENERGÉTICOS

y, en notación más sencilla:

Z Z
E I w00 (x) δw00 (x) dx = q(x) δw(x) dx. (3.87)
L L

Ahora debemos proponer una función desplazamiento real aproximada y una variación
de dicha función, ambas que cumplan con las condiciones de borde cinemáticas, es decir, que
w(0) = w(L) = 0 y δw(0) = δw(L) = 0. Para esto propongamos como función aproximada real:
x x
w(x) = c1 1− , (3.88)
L L
y como variación de la función:
 x
δw(x) = δc1 sen π . (3.89)
L
Ahora, calculemos las derivadas segundas de (3.88) y de (3.89):

2 c1
w00 (x) = − (3.90)
L2
π2  x
δw00 (x) = −δc1 2 sen π , (3.91)
L L
y reemplacemos en (3.87):

π2
Z    x  Z
2 c1  x
EI − 2 −δc1 2 sen π dx = q(x) δc1 sen π dx. (3.92)
L L L L L L

Saquemos fuera de las integrales todas las constantes. Entonces, tenemos:

π2
  Z Z
2 c1  x  x
EI − 2 −δc1 2 sen π dx = q δc1 sen π dx. (3.93)
L L L L L L

Al simplificar todo lo posible, nos queda:

2 c1 π 2
EI = q, (3.94)
L2 L2
y por lo tanto, podemos despejar c1 :

q L4
c1 = . (3.95)
2 π2 E I
Con este c1 hemos aproximado la función desplazamiento propuesta, que llamaremos
w1 (x), cuya expresión es:
q L3  x
w1 (x) = x 1 − . (3.96)
2 π2 E I L
Ahora, supongamos que tomamos otra variación de w(x). Para facilitar los procedimien-
tos, esta nueva variación de w(x) es:
x x
δw(x) = δc1 1− . (3.97)
L L
Ahora reemplazamos otra vez en (3.87):
Z    Z
2 c1 2 x x
EI − 2 −δc1 2 dx = q(x) δc1 1− dx. (3.98)
L L L L L L

Versión 2014 Preliminar - 39 -

3.6. Ejemplos de aplicación El Método de los Elementos Finitos

Si repetimos lo que hicimos antes nos queda:

Z Z
4 x x
E I c1 4 δc1 1 dx = q δc1 1− dx, (3.99)
L L L L L
L
4 x2 x3 L
E I c1 4
L δc1 = q − δc1 = q δc1 , (3.100)
L 2 L 3 L2 0 6
q L4
c1 = . (3.101)
24 E I
Con este nuevo valor de c1 armamos una nueva función w(x) que llamaremos w2 (x):

q L3  x
w2 (x) = x 1− . (3.102)
24 E I L
Nos resta averiguar cual de las dos funciones propuestas aproxima mejor el problema,
como hicimos en el caso anterior. Para ello, hagamos lo mismo que para el ejemplo anterior,
calculemos la energía potencial total del sistema con cada una de las funciones. Los resultados
son los siguientes:
Para la primera función, w1 (x), tenemos:
Z L  2 2 Z L
1 d w1 (x)
ΠP1 = EI dx − q w1 (x) dx,
0 2 dx2 0
Z L 2 Z L
q L2 q L3 x2
 
1
= EI dx − q x− dx,
2 0 π2 E I 0 2 π2 E I L
Z L Z L
q 2 L4 q L3 x2

1 (3.103)
= EI 4 dx − q x− dx,
2 π (E I)2 0 2 π2 E I 0 L
L
1 q 2 L5 q L3 x2 x3
= − q − ,
2 π4 E I 2 π2 E I 2 3L 0
1 q 2 L5 q 2 L5 π 2 − 6 q 2 L5
ΠP1 = − = − .
2 π4 E I 12 π 2 E I 12 π 4 E I
Y para la segunda función, w( 8x):
Z L  2 2 Z L
1 d w2 (x)
ΠP2 = EI dx − q w2 (x) dx,
0 2 dx2 0
Z L 2 Z L
q L2 q L3 x2
 
1
= EI dx − q x− dx,
2 0 12 E I 0 24 E I L
Z L Z L
1 q 2 L4 q 2 L3 x2

= dx − x− dx, (3.104)
2 122 E I 0 24 E I 0 L
L
1 q 2 L5 q 2 L3 x2 x3
= − − ,
2 122 E I 24 E I 2 3L 0
1 q 2 L5 q 2 L5 1 q 2 L5
ΠP2 = 2
− 2 =− .
2 12 E I 12 E I 288 E I
Si comparamos ΠP1 con ΠP2 , queda en evidencia que ΠP2 < ΠP1 , por lo tanto, la función
w2 (x) es la mejor aproximación.

- 40 - Preliminar Versión 2014

El Método de los Elementos Finitos 4. El MÉTODO DE RITZ

CAPÍTULO 4

EL MÉTODO DE RITZ

4.1. Introducción
En los capítulos anteriores nos hemos ocupado de estudiar el cálculo de variaciones y la
revisión de los distintos teoremas energéticos que se aplican en al campo de la ingeniería civil. En
el capítulo dos nos adentramos en el cálculo de variaciones y vimos que la principal aplicación
es la de minimizar (o maximizar) un funcional. También estudiamos algunas de las propiedades
fundamentales y sus operaciones más usuales.
En el capítulo tres hicimos una revisión de los teoremas energéticos más usuales en el
campo de la ingeniería civil. Luego analizamos tres teoremas, el Teorema de los Desplaza-
mientos Virtuales, el Teorema de las Tensiones Virtuales (de los Desplazamientos Virtua-
les Complementario) y el Teorema de la Mínima Energía Potencial, todos ellos mediante
la aplicación de los conocimientos adquiridos del cálculo variacional.
En todos los casos vistos, se armaba un funcional que luego se minimizaba, y partir de
allí se aseguraba que las funciones que cumplían con dichos teoremas resultaban ser la solución
del problema.
Todavía nos falta ver como utilizar estos conocimientos para hallar dicha funciones.
Hemos visto que la mayoría de los problemas que se presentan en el campo de la ingeniería
consisten en la resolución de ecuaciones diferenciales con sus condiciones de contorno o de borde
que son las que definen las características del problema. Pero también es claro que raramente se
pueden hallar las soluciones exactas de estas ecuaciones diferenciales.
Generalmente debemos abocarnos a resolver la mayoría de estos problemas mediante
métodos aproximados. Precisamente, dentro de los métodos aproximados se encuentra el Método
de los Elementos Finitos. En este capítulo veremos uno de los primeros métodos para la resolución
de problemas mediante elementos finitos. Este método es el Método de Rayleigh–Ritz (o
simplemente Método de Ritz).

4.2. El método de Ritz

El método de Ritz se basa en la existencia un funcional, como puede ser el de la energía
potencial total, que mediante la aplicación del cálculo variacional, pueda ser minimizado. Para
resolver este funcional, el método propone una función solución para dicho funcional que debe
cumplir con una serie de condiciones.

Versión 2014 Preliminar - 41 -

4.2. El método de Ritz El Método de los Elementos Finitos

Para explicar y mostrar cómo aplicar el método, supongamos por un momento un pro-
blema cualquiera cuyo funcional, por ejemplo la energía potencial total, tiene la forma:
Z  
due
Πp = F ue ; ; x dx. (4.1)
dx

y cuyas condiciones de borde son ue (x1 ) = ue (x2 ) = 0, si ue es la solución exacta, y F el

integrando, en este caso integrado en un dominio unidimensional de x.
El método propone que la solución del problema se puede aproximar por una serie de
funciones que satisfagan la condiciones de borde y que se pueden expresar de la forma:
n
X
ue (x) = ci φi (x) = c1 φ1 (x) + c2 φ2 (x) + · · · + cn φn (x)φ, (4.2)
i=1

donde ci es un coeficiente no conocido y φi (x) una función cualquiera (xi ; sen(i π x); etc.) lineal-
mente independiente. Además vamos a imponer que estas funciones deben satisfacer cada una
de ellas las condiciones de borde, es decir, se debe cumplir que:

φi (x1 ) = 0; φi (x2 ) = 0; para i = 1; 2; . . . ; n. (4.3)

De acuerdo con la (4.1), las funciones φi (x) deben ser continuas, no así las primeras
derivadas, que pueden ser discontinuas.
Reemplacemos u en el funcional de (4.1) y teniendo en cuenta que δΠP = 0, llegamos a
que:
n
X ∂Πp ∂Πp ∂Πp ∂Πp
δΠp = δci = δc1 + δc2 + . . . + δcn = 0, (4.4)
∂ci ∂c1 ∂c2 ∂cn
i=1

pues es una función de los ci . Como a su vez los δci son arbitrarios (y por lo tanto distintos de
cero), de (4.4) inferimos que:

∂Πp
δci 6= 0 ⇒ = 0; i = 1; 2; . . . ; n, (4.5)
∂ci

con lo cual obtenemos un sistema de ecuaciones lineales cuando Πp sea una función cuadrática
de u y de du
dx .
Para asegurar la convergencia del método se debe cumplir que:

1. Las funciones de aproximación deben ser continuas hasta un orden menos que la derivada
más alta del integrando. Esta condición significa tener un integrando definido.

2. Las funciones deben satisfacer en forma individual las condiciones de borde esenciales (lo
que hemos impuesto en (4.3)).

Podemos agregar como una tercera condición que la serie de funciones debe ser completa.
Se entiende por completa que el error cuadrático se anule en el límite, o sea:

n
!2
Z x2 X
lı́m = ue − ci ϕi dx = 0. (4.6)
n→∞ x1 i=1

Por todo esto es fácil ver que las funciones polinómicas y las trigonométricas son las
funciones que mejor se ajustan para aplicar el método de Ritz. Veamos algunos casos conocidos
dentro de la ingeniería civil.

- 42 - Preliminar Versión 2014

El Método de los Elementos Finitos 4. El MÉTODO DE RITZ

4.3. Ejemplos de aplicación

4.3.1. Resolución en primer orden de una viga sometida a flexión
Planteo de la funcional
Analicemos el Teorema de los Desplazamientos Virtuales. Hemos visto que el funcional
es: Z Z
σ δε dV = p(x) δw(x)dx. (4.7)
V
L

Tomemos ahora una viga como la de la figura 4.1. donde w(x) son los desplazamientos
del eje de la viga en dirección del eje z.

Figura 4.1: Viga sometida a flexión.

La ecuación de equivalencia para una viga sometida a flexión es:

Z
M= σ z dA. (4.8)
A

Según la figura 4.1 tenemos las siguientes relaciones:

ε = −z w00 ⇒ δε = −z δw00 , (4.9)

de modo que la variación de la energía interna la podemos expresar de la forma:

Z Z Z
δε U = σ δε dV = − σ z δw00 (x)dA dx. (4.10)
V L A

De acuerdo con la (4.8) podemos escribir la (4.10) como:

Z
δε U = − M δw00 (x) dx. (4.11)
l

Siguiendo con este razonamiento, el trabajo de las fuerzas exteriores está definido como:
Z
δw W = q δw(x) dx. (4.12)
L

Planteemos el funcional correspondiente:

Z Z
00
− M δw (x)dx − q(x) δw(x) dx = 0. (4.13)
L L

Versión 2014 Preliminar - 43 -

4.3. Ejemplos de aplicación El Método de los Elementos Finitos

Si integramos por partes el primer término de la expresión, obtenemos lo siguiente:

Z Z
L L
− M δw0 0 + M 0 δw 0 − M 00 δw(x) dx − q(x) δw(x) dx = 0; (4.14)
L L

y si reagrupamos los términos podemos expresar la (4.14) de esta manera:

Z
 00 L L
M + q(x) δw(x) dx + M δw0 0 − M 0 δw 0 = 0.

(4.15)
L

Para que la funcional sea nula, se deben anular todos los términos. En el caso del primero,
lo que se debe anular es el corchete, es decir:

M 00 + q(x) = 0, (4.16)

que es la ecuación diferencial de equilibrio de una viga sometida a flexión, mientras que los tér-
minos restantes corresponden a las condiciones de borde del problema, que pueden ser esenciales
o naturales:

M0 δw00 = 0, (4.17a)
0
ML δwL = 0, (4.17b)
0


Q0 δw0 = 0 M =Q , (4.17c)
QL δwL = 0. (4.17d)

Resolución de la funcional mediante el método de Ritz

Regresemos a la ecuación (4.13)y analicemos como resolver este funcional. En primer
lugar, recordemos que M = −EI w00 , por lo que la ecuación se puede escribir de la siguiente
forma: Z Z
00 00
EI w δw (x) dx − q(x) δw(x) dx = 0 (4.18)
L L

Ritz propone para resolver esta funcional adoptar una función w(x) de la forma:
X
w(x) = ci ϕi (x), (4.19)
i

que cumpla con las condiciones de borde cinemáticas del problema. A partir de esta función se
deducen la primera derivada: X
w0 (x) = ci ϕ0i (x), (4.20)
i

y la segunda: X
w00 (x) = ci ϕ00i (x). (4.21)
i

Hallemos las variaciones de todas estas funciones:

X ∂ci
δcj w(x) = ϕi (x) δcj = ϕj (x) δcj , (4.22)
∂cj
i

∂ci ∂ci
pues ∂c j
= 0 para i 6= j y ∂c j
= 1 para i = j. Esto asegura que los ci sean linealmente
independientes.
Análogamente tenemos la variación de la derivada primera:

δcj w0 (x) = ϕ0j (x) δcj , (4.23)

- 44 - Preliminar Versión 2014

El Método de los Elementos Finitos 4. El MÉTODO DE RITZ

y la variación de la derivada segunda:

δcj w00 (x) = ϕ00j (x) δcj . (4.24)
Para simplificar un poco la notación, definamos B = E I y reemplacemos las ecuacio-
nes (4.21), (4.22) y (4.24) en la (4.18):
Z X Z
00 00
B ci ϕi (x) ϕj (x) δcj dx = q(x)ϕj (x) δcj dx. (4.25)
L i L

Como δcj 6= 0, y por definición es independiente de x, lo mismo que los coeficientes ci , si

permutamos la sumatoria y la integral, podemos reordenar la (4.25) así:
" # Z 
X Z
00 00
ci B ϕi (x) ϕj (x) dx δcj = q(x)ϕj (x) dx δcj , (4.26)
i L L

con lo cual los términos que debemos igualar son los que están entre corchetes. Así:
Xn Z Z
00 00
ci B ϕi (x) ϕj (x) dx = q(x) ϕj (x) dx para j = 1; 2; . . . ; n. (4.27)
i=1 L L

Esta expresión nos muestra que nuestras incógnitas son los coeficientes cj de la función
w(x) propuesta. Lo que obtendremos es un sistema de n ecuaciones con n incógnitas siendo n el
cantidad de coeficientes cj , es decir, obtendremos los cj para j = 1; 2; . . . ; n.
Una última consideración importante. Cualquier función w(x) que se proponga como
solución debe cumplir solamente con las condiciones de borde cinemáticas, en razón que las
condiciones de borde estáticas están implícitas en la funcional, como vimos en la introducción al
cálculo de variaciones.

Ejemplo
Con los conceptos desarrollados el punto anterior y las expresiones encontradas, resolva-
mos en forma práctica un ejemplo de viga. Tomemos una viga simplemente apoyada de longitud
L con una carga concentrada P en el centro de la luz, como se ve en la figura 4.2.

Figura 4.2: Viga sometida a carga concentrada.

Propongamos la siguiente función desplazamiento real w(x):

x x
w(x) = c1 1− . (4.28)
L L
Primero, verifiquemos que la función cumple con las condiciones de borde cinemáticas,

es
decir, que w(x = 0) = w(x = L) = 0. En efecto, para x = 0 tenemos que w(0) = c1 L0 1 − L0 = 0
y para x = L, w(L) = c1 L L


L 1 − L = 0.
En segundo lugar, hallemos la primera y la segunda derivadas de la función w(x):
 
0 c1 2x
w (x) = 1− (4.29)
L L
c1
w00 (x) = −2 2 . (4.30)
L

Versión 2014 Preliminar - 45 -

4.3. Ejemplos de aplicación El Método de los Elementos Finitos

Finalmente, hallemos las variaciones de w(x), w0 (x) y w00 (x). Como las funciones de las
variaciones son las mismas que las propuestas como reales, se verifica que δw(0) = 0 y δw(L) = 0.
Entonces:
x x
δc1 w(x) = 1− δc1 (4.31)
L L
1  x 
δc1 w0 (x) = 1−2 δc1 (4.32)
L L
2
δc1 w00 (x) = − 2 δc1 . (4.33)
L

Reemplacemos estas expresiones en la ecuación (4.27), e integrando el segundo miembro,

nos queda:
Z
c1 P
EJ 2 2 dx δc1 = δc1 , (4.34)
L4 L 4

donde φ00 (x) = − L22 . Por último, integremos el primer miembro de la ecuación, y hallemos c1 :

c1 P P L3
4
E J 4L = ⇒ c1 = (4.35)
L 4 16EJ

Así la función w(x) queda:

P L2  x
w(x) = x 1− (4.36)
16 E J L
L L



Calculemos con esta función w x = 2 yM x= 2 . En el primer caso tenemos:

P L2 L P L3
     
L L L
w = 1− ⇒ w = (4.37)
2 16 E J 2 2 2 64 E J

PL 3
y comparada con la solución exacta, que es fmáx = 48 E J , vemos que el desplazamiento calculado
mediante el Método de Ritz es un 25 % menor, por lo tanto el sistema resulta ser más rígido que
el exacto.
Para calcular M x = L2 , hallemos M (x):

PL
M (x) = , (4.38)
8

una función constante.

Figura 4.3: Diagrama de momento flector.

Observemos que esta función no es la solución para una viga simplemente apoyada, pues
el diagrama de momentos para este caso es el que se ve en la figura 4.3. Lo que hemos obtenido
al aplicar el método de Ritz es el momento flector medio en la viga.

- 46 - Preliminar Versión 2014

El Método de los Elementos Finitos 4. El MÉTODO DE RITZ

4.3.2. Barra solicitada axilmente

Podemos hacer el mismo planteo para una barra solicitada axilmente. Sabemos que en
este caso la funcional es: Z Z
σ δε dV = p(x) δu(x) dx. (4.39)
V L

Adoptemos una función u(x) tal que sea de la forma:

X
u(x) = ai φi (x); (4.40)
i

du
y como además sabemos que σ = E ε y ε = dx , entonces podemos expresar la (4.39) así:
Z   Z
du du
AE δ dx = p(x) δu(x) dx. (4.41)
L dx dx L

De la (4.40) podemos deducir que:

X
u0 (x) = ai φ0i (x), (4.42)
i

y a partir de ambas expresiones podemos hallar las variaciones de las dos:

δaj u(x) = aj φj (x) (4.43)

0
δaj u (x) = aj φ0j (x). (4.44)

Reemplazando las expresiones (4.42), (4.43) y (4.44) en la (4.41), nos queda:

Z X Z
AE ai φ0i (x) φ0j (x) δaj dx = p(x)φj (x) δaj dx. (4.45)
L i L

Aquí también podemos permutar la sumatoria y la integral, y si A y E no son funciones

de x, que podemos escribir la expresión (4.45) de la forma:
" Z # Z 
X
0 0
ai A E φi (x) φj (x) dx δaj = p(x)φj (x) dx δaj . (4.46)
i L L

Como los δaj 6= 0, necesariamente las expresiones en los corchetes deben ser iguales, es
decir: Z Z
X
0 0
ai A E φi (x) φj (x) dx = p(x) φj (x) dx para j = 1; 2; . . . ; n, (4.47)
i L L

que es similar a la expresión (4.27). En este caso podemos proponer una función para u(x) que
sea lineal como solución aproximada del problema.

4.3.3. Viga a flexión en teoría de segundo orden

Vimos anteriormente el funcional de un sistema analizado por la teoría de segundo orden.
La funcional para ese caso es:
Z Z
δΠ = E I w00 δw00 dx − P w0 δw0 dx = 0. (4.48)
L L

Tomemos el caso de una barra simplemente apoyada con una carga horizontal axil, como
se ve en la figura 4.4.

Versión 2014 Preliminar - 47 -

4.3. Ejemplos de aplicación El Método de los Elementos Finitos

Figura 4.4: Viga con carga horizontal.

Adoptemos como función w(x) la siguiente expresión:

n  
X iπ
w(x) = ci sen x , (4.49)
L
i=1

y hallemos las derivadas primera y segunda de w(x):

n  
0
X iπ iπ
w (x) = ci cos x , (4.50)
L L
i=1
n  2  
00
X iπ iπ
w (x) = − ci sen x . (4.51)
L L
i=1

El paso siguiente es hallar las variaciones de estas últimas funciones. Estas variaciones
son:
 
0 jπ jπ
δcj w (x) = cos x δcj , (4.52)
L L
 2  
00 jπ jπ
δcj w (x) = − cos x δcj . (4.53)
L L

Reemplacemos (4.50), (4.51), (4.52) y (4.53) en la (4.43), de modo que quede:

Z n  2     2   
X iπ iπ jπ jπ
EI ci − sen x − sen x δcj dx+
L L L L L
i=1
Z X n     (4.54)
iπ iπ jπ jπ
−P ci cos x cos x δcj dx = 0
L L L L L
i=1

Si agrupamos los términos, nos queda:

( n  2  2 Z       )
X iπ jπ iπ jπ
ci E I − sen x − sen x dx δcj +
L L L L L
i=1
" n     # (4.55)
X iπ jπ Z iπ jπ
−P ci cos x cos x dx δcj = 0
L L L L L
i=1

Resolvamos la ecuación, teniendo en cuenta:

Que las integrales para este caso sólo tienen sentido cuando i = j,
L
En ese caso valen 2.

Con estas premisas, la (4.55)queda:

(  4  2 )
jπ L jπ L
cj E I −P δcj = 0. (4.56)
L 2 L 2

- 48 - Preliminar Versión 2014

El Método de los Elementos Finitos 4. El MÉTODO DE RITZ

Para que la ecuación (4.56) se anule para cualquier δcj , resulta evidente que se debe anular
el término entre llaves. Observemos que al anularse éste, implica que la ecuación se hace nula
también para cualquier cj , es decir, que la solución del problema es independiente de la función
desplazamiento específicamente propuesta, esto es, la solución es independiente del coeficiente,
pero no lo es de la familia de funciones a la que pertenece, que es la que define la expresión.
Por lo tanto podemos hallar P en función de E I y L, y resulta:

E Iπ 2
P = . (4.57)
L2

La solución a la que hemos arribado es la carga crítica de Euler para una barra sim-
plemente apoyada con una carga axil de compresión. Es claro que hemos hallado la solución
«exacta» de este problema, pues partimos suponiendo como solución aproximada para la función
desplazamiento una función que pertenece a la misma familia de funciones a la que pertenece la
solución exacta del problema. Si, en cambio, hubiésemos propuesto como solución aproximada
una función perteneciente a otra familia de funciones, (por ejemplo, una función polinómica), la
solución hubiera sido simplemente una aproximación.

4.3.4. Transmisión de calor: caso unidimensional

Veremos a continuación como resolver mediante el Método de Ritz un problema de trans-
misión del calor unidimensional.

Ecuación unidimensional de transmisión de calor

Determinaremos la distribución de temperatura debido a la conducción de calor a través

de una placa gruesa con un coeficiente de conductividad térmica k(x), como se ve en la figura 4.5.
La placa está cargada por un flujo de calor q, aplicado en la superficie en x = a, mientras que en
la otra superficie está sometida a una temperatura θ = Tb . Por último, asumamos que s(x) = q B
define la distribución de la generación interna de calor.

Figura 4.5: Transmisión del calor unidimensional.

La ecuación diferencial que describe la distribución estacionaria de temperatura es:

 
d dθ
k − q B = 0 con a < x < b. (4.58)
dx dx

Versión 2014 Preliminar - 49 -

4.3. Ejemplos de aplicación El Método de los Elementos Finitos

Las condiciones de borde para el problemas son:

dθ
−k − q = 0 con x = 0, y (4.59)
dx
θ = Tb para x = b. (4.60)

La dirección positiva del flujo de calor es paralela a x; el flujo de calor ingresa a la placa
en x=a, que toma en cuenta el signo de q en (4.59).
La solución para (4.58) a (4.60) se puede hallar analíticamente pues k(x) y s(x) (en este
caso q B ) son funciones simples. Para k y q B constantes, uno integra directamente (4.58) dos veces
y determina las constantes de integración imponiendo las condiciones de borde (4.59) y (4.60).
Esta solución es independiente del origen del eje x; así para L = b−a y definiendo a = 0 en (4.58)
y (4.59), la solución del problema es:

q B L2
  x 2  q L  x
θ (x) = 1− + 1− + Tb . (4.61)
2k L k L

Planteo Variacional

Tomemos la ecuación diferencial (4.58) y definamos el funcional, tal como hacemos en

análisis estático. La funcional nos queda:
Z    
d dθ
Π= − k − q B θ dx. (4.62)
L dx dx

Esta funcional que armamos es muy similar al funcional que ya conocemos para solicita-
ción axil. Es: Z    
d du
Π= E A(x) + t u dx. (4.63)
L dx dx
Variemos la funcional para minimizarlo y operemos de modo de reducir el orden de la
ecuación. Así, integrando por partes, nos queda:
Z Z L
d δθ d θ dθ
δΠ = k dx − q B δθ dx − δθ k = 0. (4.64)
L dx dx L dx 0

Comparemos otra vez con el mismo caso de (4.63) una vez desarrollado:
Z Z L
d u d δu du
δΠ = E A(x) dx − tδu dx − δu E A(x) = 0. (4.65)
L dx dx L dx 0

Esta ecuación la hemos resuelto mediante el método de Ritz para el caso de la barra
solicitada axilmente. Por lo tanto, como ambas ecuaciones son análogas, podemos concluir di-
ciendo que el último término de la expresión (4.64) es la condición de borde (4.59), y entonces
reescribimos la ecuación como:
Z Z
d δθ d θ
δΠ = k dx − q B δθ dx − δθ q|0 = 0. (4.66)
L dx dx L

Para resolver este sistema, podemos aplicar la misma mecánica utilizada en los sistemas
estructurales y definir una función de θ que cumpla con las condiciones de borde de temperatura
vistas en (4.60). Estas condiciones son análogas a las condiciones de borde cinemáticas que se
deben cumplir en análisis de estructuras.

- 50 - Preliminar Versión 2014

El Método de los Elementos Finitos 4. El MÉTODO DE RITZ

Resolución del sistema

Planteemos como función θ(x) la siguiente expresión:
X
θ(x) = Ti φi (x). (4.67)
i

Si tomamos un solo elemento para resolver la ecuación, entonces nuestra función queda
de la siguiente manera:  x x
θ(x) = T1 1 − + Tb , (4.68)
L L
donde nuestra única incógnita es T1 , pues Tb es dato [expresión (4.60)], y las funciones interpo-
lantes son: φ1 (x) = 1 − Lx y φ2 (x) = Lx .
Desarrollemos la expresión (4.66) por el método de Ritz. La expresión queda entonces
como: " #
XZ Z
0 0 B
φj k Ti φi dx − q φj dx − q|0 δTj = 0. (4.69)
i L L

Hallemos las derivadas de φ1 (x) y φ2 (x), y la variación de θ0 (x):

1
φ01 (x) = − (4.70)
L
1
φ02 (x) = (4.71)
L
1
δθ0 (x) = − δT1 . (4.72)
L
Como observamos, la variación de θ0 (x) respecto de Tb es nula, pues es un valor ya definido.
Análogamente pasa con la variación de θ(x). Por lo tanto, sólo queda la variación respecto a T1 ,
que la expresamos de la siguiente manera:
 x
δθ(x) = 1 − δT1 . (4.73)
L
Reemplacemos las expresiones (4.70), (4.71), (4.72) y (4.73) en (4.69)):
Z   
Tb − T1
Z
1 B
 x
− k dx − q 1− dx − q|0 δT1 = 0. (4.74)
L L L L L

Como la solución que buscamos debe ser distinta de la trivial (T1 = 0), el corchete debe
anularse. La expresión para T1 una vez resuelto el corchete, es:

q B L2 q L
T1 = + + Tb , (4.75)
2k k
y reemplazando en (4.68), nos queda la expresión de la función θ(x):
 B 2 
q L qL  x
θ(x) = + 1− + Tb . (4.76)
2k k L

Si comparamos las funciones halladas por los métodos, es decir, mediante la resolución
de la ecuación diferencial y mediante el método de Ritz, veremos claramente la diferencia entre
ambas. La solución de la ecuación diferencial (4.61) es una parábola de segundo grado, mientras
que la solución por Ritz (4.76) es una función lineal.

Versión 2014 Preliminar - 51 -

El Método de los Elementos Finitos 5. ELEMENTO DE BARRA

CAPÍTULO 5

ELEMENTO DE BARRA

5.1. Introducción
Hemos visto en el capítulo anterior la forma de resolver un sistema de ecuaciones diferen-
ciales mediante la aplicación del método de los elementos finitos. Utilizamos para ello el método
de Rayleigh - Ritz, con el cual proponíamos como solución aproximada una función (u(x) o w(x)
en nuestros ejemplos).
Vimos, además, que proponiendo funciones adecuadas como solución del problema, la
aproximación convergía rápidamente a la solución exacta, asegurando de este modo que los
resultados obtenidos son confiables, y por lo tanto, aplicables para el análisis de los diferentes
problemas encarados.
Aquí aparece una nueva cuestión: ¿cuán fácil es poder determinar que la función que
uno propone como solución aproximada de ese problema es una función adecuada, si en primera
instancia no se conoce dicha función? O dicho de otro modo, para poder resolver el problema,
¿necesitamos saber primero la solución? En los ejemplos anteriores la cuestión era sencilla, debido
a que ya conocíamos con bastante certeza la familia de funciones a la cual pertenecía la solución
del problema, por haber resuelto la ecuación diferencial en forma exacta.
¿Qué hacemos entonces cuando nos encontramos ante un problema cuya solución igno-
ramos por completo? Aparece entonces el concepto de elemento finito en toda su dimensión. En
este capítulo vamos a desarrollar en forma completa, la formulación de un elemento en particular,
con el objetivo de sistematizar la utilización del método de los elementos finitos.
Hasta ahora nos hemos ocupado de estudiar la base teórica del método, resolviendo una
serie de ejemplos prácticos, pero siempre tomando el sistema completo. De ahora en más nos
concentraremos en el estudio de diferentes elementos, para los cuales desarrollaremos el algoritmo
de resolución y la manera de vincularse entre ellos.
En primer lugar estudiaremos el elemento de barra. Para ello comenzaremos analizando
el elemento a partir del teorema de los desplazamientos virtuales.

5.2. Elemento de barra con dos nodos

Analicemos primeramente un elemento de barra genérico, de acuerdo con la figura 5.1.
En este elemento tenemos como incógnitas los desplazamientos UA y UB . Adoptemos la
siguiente función de u(x):  x x
u(x) = 1 − UA + UB , (5.1)
L L
que cumple con las condiciones de borde: u(x = 0) = UA , y u(x = L) = UB .

Versión 2014 Preliminar - 53 -

5.2. Elemento de barra con dos nodos El Método de los Elementos Finitos

Figura 5.1: Elemento de barra.

Esa misma función u(x) la podemos expresar también de modo matricial:

 
 x x
 UA
u(x) = 1 − L L → u(x) = HL (x) UL , (5.2)
UB
donde HL (x) es la matriz de interpolación y UL es la matriz de los desplazamientos nodales.
Planteemos el término de la variación de la energía interna de deformación, considerando
que tanto σ y ε son vectores. Entonces nos queda lo siguiente:
Z
δU = δεεT σ dV. (5.3)
V

En esta expresión tenemos, en primer lugar, que:

 
d u(x) d [HL (x)UL ]  1 1  UA
ε= = = −L L ⇒ ε = BL (x) UL , (5.4)
dx dx UB
donde
d [HL (x)UL ]
= BL (x) UL . (5.5)
dx
Por otro lado, tenemos que:
σ = E ε = E BL (x) UL , (5.6)
donde E es el módulo de elasticidad del material de la barra.
Reemplacemos las expresiones (5.4) y (5.6) en la (5.3), con lo cual obtendremos lo si-
guiente:
Z
δU = δUL T BL (x)T E BL (x) UL dV
ZV Z
= δUL T BL (x)T E BL (x) UL dA dx
L A
Z Z  (5.7)
= δUL T BL (x)T E BL (x) dA dx UL ⇒
L A
T
δU = δUL KL UL ,
donde KL es la matriz de rigidez local y cuya expresión es:
Z Z
KL = BL (x)T E BL (x) dA dx. (5.8)
L A

Esta formulación parte de aplicar el Método de Ritz para una barra solicitada axilmente,
pues hemos utilizado la misma función para definir la variación de los desplazamientos.
Ahora planteemos la expresión del trabajo de las fuerzas exteriores:
Z Z
δW = p(x) δu(x) dx = δUL T HL (x)T HL (x) PL dx
L Z L

= δUL T T
HL (x) HL (x) PL dx ⇒ (5.9)
L
T
δW = δUL RL ,

- 54 - Preliminar Versión 2014

El Método de los Elementos Finitos 5. ELEMENTO DE BARRA

donde Z
RL = HL (x)T HL (x) PL dx, (5.10)
L

es el vector de cargas nodales equivalentes local.

Con estas dos expresiones podemos plantear el Teorema de los Desplazamientos Virtuales:

δU = δW. (5.11)

Si ahora reemplazamos las ecuaciones (5.7) y (5.9) en la (5.11), obtenemos la siguiente

expresión:
δUL T KL UL = δUL T ⇒ δUL T (KL UL − RL ) = 0. (5.12)

Resulta evidente que para que sea nula la expresión para cualquier variación de U, es
necesario que el término entre paréntesis sea nulo. Por lo tanto:

KL UL − RL = 0 ⇒ KL UL = RL . (5.13)

Hasta aquí hemos analizado simplemente un elemento. Para el caso de un sistema de

varias barras, utilizamos este mismo razonamiento, ensamblando las distintas matrices KL <i> y
los vectores RL <i> , con i = 1; 2; . . . ; n, siendo n es la cantidad de elementos en que se dividió
el sistema, de manera de obtener la matriz de rigidez global, KG , y el vector de cargas nodales
equivalente global del sistema. Así, lo que obtenemos es la misma ecuación anterior pero de la
forma:
KG UG − RG = 0 ⇒ KG UG = RG , (5.14)
donde la g indica que el sistema es el global y no el local (l).

5.3. Ejemplo
Veamos un ejemplo sencillo como aplicación del método. Tomemos el sistema de la figu-
ra 5.2.

Figura 5.2: Sistema con barras axilmente cargadas.

Vemos que en este ejemplo tenemos un tramo de de la barra con sección constante y el
otro, con sección variable. Fijemos primero las incógnitas globales del sistema, que son U1 y U2 ,
como vemos en la figura 5.3.
Por las características del sistema, vamos a discretizarlo en dos elementos de barra para
analizar el sistema. Como primer elemento (e<1> ) definimos el tramo de barra de sección cons-
tante, de longitud L1 . El segundo y último elemento (e<2> ) es el de sección variable linealmente.
Ahora tomemos la discretización e identifiquemos las incógnitas locales para cada uno
de los elementos (figura 5.4). En cada elemento tenemos dos incógnitas locales (UA y UB ). En
consecuencia, la discretizar tenemos cuatro incógnitas, que son: UA <1> , UB <1> , UA <2> y UB <2> .

Versión 2014 Preliminar - 55 -

5.3. Ejemplo El Método de los Elementos Finitos

Figura 5.3: Incógnitas globales.

Figura 5.4: Discretización del sistema e incógnitas locales.

Armemos la matriz de rigidez KL<1> del primer elemento. Adoptemos las mismas funciones
ya vistas en los puntos anteriores como funciones desplazamiento, por lo tanto, la matriz de rigidez
del primer elemento es:
 1
− L11

L1
Z
T
KL<1> = A1 B<1> E B<1> dx = A1 E  , (5.15)
L1 1 1
− L1 L1

donde A1 es una constante y B <1> es una matriz con coeficientes constantes.

Tomemos ahora el segundo elemento y hallemos la matriz de rigidez, KL<2> . Adoptemos
también las mismas funciones desplazamientos que en el caso anterior. Sin embargo, en este caso
la sección ya no es constante y no puede salir de la integral. La integral resulta ser, entonces:
 1
− L12

L2
Z
KL <2 =
T
B<2> E B<2> A2 (x) dx = E   A1 + A2 . (5.16)
2
L2 − L12 1
L2

Por último, hallemos los vectores de cargas nodales de cada elemento. En el caso del
elemento 1, la matriz de cargas es nula, pues no hay cargas en los nodos A y B. En cambio, en
el elemento 2 tenemos una carga concentrada en el nodo B. Los vectores son:
   
<1> 0 <2> 0
R = y R = . (5.17)
0 P
Ahora debemos ensamblar las matrices locales de rigidez para hallar la matriz global
de rigidez. Lo mismo haremos con los vectores de cargas nodales. Primero, relacionemos las
incógnitas locales con las incógnitas globales mediante la matriz topológica:
 
0 U1
MT = , (5.18)
U1 U2
donde la columna 1 representa al elemento 1, y contiene los desplazamientos de ese elemento, y
la columna 2 representa al elemento 2 y contiene los desplazamientos de ese elemento; en ambos
casos, en términos de los desplazamientos globales. Así, la matriz la leemos de esta forma:
 <1>     <2>   
UA 0 UA U1
<1> = y <2> = .
UB U1 UB U2

- 56 - Preliminar Versión 2014

El Método de los Elementos Finitos 5. ELEMENTO DE BARRA

Con esto, ensamblemos las matrices de rigidez locales y armemos la matriz de rigidez
global, que queda así:
 A 
L1
1
−A L1
1
0
 
 
 A1 A1 A1 +A2 A1 +A2 
KG = E − L1 L1 + 2 L2 − 2 L2  , (5.19)
 
 
0 − A21 +A
L2
2 A1 +A2
2 L2

que resulta ser simétrica, como en los casos ya vistos.

De la misma forma, hallemos la matriz de cargas nodales equivalentes global:
 
0
RG = 0  .
 (5.20)
P

Lo que falta es armar el sistema KG UG = RG :

 A 
1
−A 1
0
   
L1 L1 U0 0
     
 
 A1 A1 A1 +A2 A1 +A2 
   
E − L1 L1 + 2 L2 − 2 L2  U1  =  0  .
    (5.21)
     
 
0 − A21 +A
L2
2 A1 +A2
2 L2
U2 P

Ahora de debemos imponer las condiciones de borde del sistema. En nuestro ejemplo es
U0 = 0, entonces podemos eliminar la primera fila y la primera columna. Nuestro sistema se
reduce al siguiente:
 A1 A1 +A2
− A21 +A
   
L1 + 2 L2
2
L2 U1 0
E    =  . (5.22)
− A21 +A
L2
2 A1 +A2
2 L2
U2 P
Si resolvemos el sistema de ecuaciones, obtenemos los valores de U1 y U2 :
 P L1 
  E A1
U1
=   . (5.23)
 
U2 P L1 2 L2
E A1 + A1 +A2

Con estos dos valores podemos calcular la deformación ε y las tensiones σ para cada
elemento. La deformación ε<1> es:
 
0
ε<1> = B1 U1 = − L11 L11  = P ,
 
(5.24a)
P L1 E A1
E A1

ε<2> es:  P L1 
E A1
2P
ε<2> = B2 U2 = − L12 1
 
 = E (A + A ) . (5.24b)

L2  
P L1 2 L2 1 2
E A1 + A1 +A2

Con ambas deformaciones, calculemos las tensiones σ <1> y σ <2> :

P P
σ <1> = E ε<1> = E ⇒ σ <1> = (5.24c)
E A1 A1
2P 2P
σ <2> = E ε<2> =E ⇒ σ <2> = . (5.24d)
E (A1 + A2 ) A1 + A2

Versión 2014 Preliminar - 57 -

5.3. Ejemplo El Método de los Elementos Finitos

Observemos que para el primer tramo de la barra la solución obtenida es la solución

«exacta», mientras que para el segundo tramo la solución es el promedio de ambos extremos.
En realidad, el método calcula la deformaciones y las tensiones en la sección media de
cada elemento. En el primer elemento, como la sección es constante, tanto la deformación como
la tensión halladas coinciden con la tensión y la deformación «exactas». Pero en el segundo
elemento, la sección de la barra no es constante, y en consecuencia, al calcular los valores en la
sección media, los resultados obtenidos arrastran un error con respecto a la solución «exacta».

- 58 - Preliminar Versión 2014

El Método de los Elementos Finitos 6. ELEMENTO DE VIGA

CAPÍTULO 6

ELEMENTO DE VIGA

6.1. Introducción
En el capítulo anterior estudiamos a los elementos de barra, aquellos que sólo pueden ser
solicitados por cargas axiles. En este capítulo nos concentraremos en el análisis de los elementos
de vigas, elementos que pueden ser solicitados a flexión, tanto simple como compuesta.
Primero estudiaremos el elemento despreciando la influencia en las deformaciones del
esfuerzo de corte y hallaremos la matriz de rigidez. Más adelante nos ocuparemos de introducir
la deformación por corte (distorsión) y veremos que dificultades nos agrega al cálculo de las
deformaciones.

6.2. Cálculo de la matriz de rigidez

Consideraremos que cada elemento tiene tres incógnitas desplazamiento por nodo: u, w
y θ. Analicemos, entonces, en un punto cualquiera de un elemento de viga los desplazamientos
de los puntos en una sección, según la figura 6.1.

Figura 6.1: Elemento de viga.

Definamos la energía interna de deformación para el elemento:

Z Z Z Z Z Z Z
00
UM = σ ε dV dε = σ dx dA d(−z w ) = − M d(w00 ) dx. (6.1)
ε∗ V w∗ A L w∗ L

Versión 2014 Preliminar - 59 -

6.2. Cálculo de la matriz de rigidez El Método de los Elementos Finitos

Tomemos un elemento de barra e identifiquemos todas sus incógnitas desplazamiento.

(Figura 6.2).

Figura 6.2: Incógnitas del elemento de viga.

Fijemos las funciones para u(x), w(x), y w0 (x). Empecemos con u(x):

u(x) = α0 + α1 x, (6.2)

y finalmente, con la función desplazamiento para w(x):

w(x) = α2 + α3 x + α4 x2 + α5 x3 . (6.3)

La función de θ(x) la definimos como w0 (x), por lo que nos queda:

θ(x) = w0 (x) = α3 + α4 2 x + α5 3 x2 . (6.4)

Expresemos las funciones u(x) y w(x) en forma matricial:

 
α0
α1 
    
u(x) 1 x 0 0 0 0  α2  .

= (6.5)
w(x) 0 0 1 x x2 3
x  α3 

α4 
α5

que también podemos expresar como:

U(x)<n> = Φ (x)<n> α <n> , (6.6)

donde el supraíndice n indica el elemento n en que se discretizó la viga.

Armemos la matriz para los desplazamientos de cada uno de los nodos del elemento de
la figura 6.2, es decir uA , wA , θA , uB , wB y θB , suponiendo que el elemento tiene una longitud
L:     
uA 1 0 0 0 0 0 α0
wA  0 0 1 0 0 0  α1 
    
 θA  0 0 0 1 0 0  α2 
 =    ⇒ UL <n> = A<n> α <n> . (6.7)
 uB  1 L 0 0 0 0  α3 
    
wB  0 0 1 L L2 L3  α4 
θB 0 0 0 1 2 L 3 L2 α5

Despejemos α en función de UL <n> y expresemos U(x)<n> también en función de

UL <n> :

α = A−1 UL <n> ⇒ U(x)<n> = Φ L (x)<n> A1 UL <n> = HL (x)<n> UL <n> . (6.8)

- 60 - Preliminar Versión 2014

El Método de los Elementos Finitos 6. ELEMENTO DE VIGA

De esta forma expresamos la función desplazamiento en función de los desplazamiento no-

dales de la barra, y por consiguiente, de las incógnitas de nuestro problema. La matriz HL (x)<n>
tiene la siguiente forma:
1 − Lx x
 
0 0 L 0 0
HL (x)<n> =  . (6.9)
x2 x3 x2 x3 x2 x3 x2 x3
0 1 − 3 L2 + 2 L3 x − 2 L + L2 0 3 L2 − 2 L3 − L + L2
La función propuesta w(x) está obtenida mediante interpolación por el método de Her-
mite, en tanto que la función u(x) es una función obtenida mediante interpolación de Lagrange 1 .
Hallemos la energía interna de deformación para la viga. La vamos a separar en energía
debida al esfuerzo axil y en energía debida a la flexión.
La energía debida al esfuerzo axil es:
Z
δUN = N (x) δε(x) dx, (6.10)
L

y la debida a flexión: Z
δUM = M (x) δw00 (x) dx. (6.11)
L
Definamos ahora el vector ε (x)<n> :
   
du(x) d() 
dx dx 0 u(x)
ε (x)<n> =  = . (6.12)
   

2
d w(x) d2 ()
−z dx2 0 −z dx 2
w(x)

Reemplacemos en (6.12) las funciones u(x) y w(x) por su forma matricial:

d
ε (x)<n> = [HL (x)<n> ] UL <n> = BL (x)<n> UL <n> , (6.13)
dx
donde BL (x)<n> es la matriz derivada primera de HL (x)<n> respecto de x.
Ahora definamos al vector de tensiones σ (x)<n> :

σ (x)<n> = E ε (x)<n> = E BL (x)<n> UL <n> . (6.14)

Por último, nos queda hallar la expresión para la energía interna. Para ello, tomemos
como base las expresiones (6.13) y (6.14):
Z Z Z
<n> T <n> T T
δε U = δεε σ dV = δUL <n> BL (x)<n> E BL (x)<n> UL <n> dA dx, (6.15)
V L A

y, si sacamos fuera de la integral los términos con UL <n> que no son funciones de x, nos queda:
Z Z 
<n> T <n> T <n>
δε U = δUL BL (x) E BL (x) dA dx UL <n> , (6.16)
L A

donde el término entre corchetes resulta ser la rigidez KL <n> :

Z Z
<n> T
KL = BL (x)<n> E BL (x)<n> dA dx. (6.17)
L A

Nos falta todavía encontrar el vector de cargas del trabajo de las fuerzas exteriores.
Tomemos el vector de las cargas del elemento:
 
<n> px (x)
pL (x) = , (6.18)
pz (x)
1
Para mayores detalles sobre funciones de interpolación ver Análisis Numérico de R. Burden y J.D. Faires,
International Thomson Editions. 6a Edición, o cualquier otro libro de análisis numérico.

Versión 2014 Preliminar - 61 -

6.3. La hipótesis de Timoshenko El Método de los Elementos Finitos

a partir del cual obtenemos el término del trabajo de las fuerzas exteriores:
Z Z  
<n> T <n>
  px (x)
δU W = δU(x) pL dx = δu(x) δw(x) dx. (6.19)
L L pz (x)

Como hemos definido U(x)<n> en función de los desplazamientos nodales UL <n> , el

trabajo de las fuerzas exteriores queda expresado de la siguiente forma:
Z
T T
δU W = δUL (x)<n> HL (x)<n> pL (x)<n> dx
L Z
= δUL (x) <n> T T
HL (x)<n> pL (x)<n> dx ⇒ (6.20)
L
<n> T
δU W = δUL (x) RL <n> ,
donde Z
T
RL <n> = HL (x)<n> pL (x)<n> dx. (6.21)
L
Ahora sí, al igualar la variación de la energía interna de deformación con la variación del
trabajo de las fuerzas exteriores nos queda
T T
δU U = δU W ⇒ δUL <n> KL <n> UL <n> = δUL <n> RL <n> ⇒
T
(6.22)
δUL <n> KL <n> UL <n> − RL <n> = 0 ⇒ KL <n> UL <n> = RL <n> ,

que es, nuevamente, la ecuación tradicional para resolver nuestro problema (K U = R).

6.3. La hipótesis de Timoshenko

Hasta aquí hemos visto el caso del elemento de viga suponiendo solamente el cumplimiento
de la hipótesis de Bernoulli-Navier, es decir, que las secciones transversales se mantienen planas
y normales al eje del elemento luego de la deformación.
Esta hipótesis no tiene en cuenta las deformaciones debidas al corte, deformaciones que
introducen en la sección una distorsión. Esta nueva hipótesis, conocida como la hipótesis de
Timoshenko, sostiene que la sección transversal se mantiene plana pero no forma un ángulo de
90o respecto al eje del mismo, como se puede ver en la figura 6.3.
Calculemos en primer lugar la energía interna de deformación por corte según la ecuación
clásica: Z Z Z
1 1
δUQ = γ τ dz = b γ τ dz. (6.23)
2 b h 2 h
τ
Si reemplazamos γ por G, e integramos en z, nos queda:

1 τ2
δUQ = A. (6.24)
2 G
Según ya se vio en la teoría de Jouravsky para flexión y corte en vigas, la distribución
del τ no es constante sino variable, lo que se traduce en que la sección en realidad se alabea.
Por lo tanto, es necesario introducir un factor de corrección para aquellas secciones no circulares,
a fin de mantener la relación τ = G γ. En el caso de la sección rectangular, ese coeficiente de
corrección de la sección, que podemos llamar κ, es 5/6.
Planteemos ahora las funciones que intervienen:

ε T = −z β 0 w0 − β
 
(6.25)

que son las deformaciones, tanto a flexión como por corte, y:

σ T = σx τxy
 
(6.26)

- 62 - Preliminar Versión 2014

El Método de los Elementos Finitos 6. ELEMENTO DE VIGA

Figura 6.3: Hipótesis de Timoshenko.

Ahora introduzcamos estas dos matrices en la expresión de la energía interna. Nos queda:
Z Z Z
0
δε U = T
δε σ dV = δ(−z β ) σx dV + δ(w0 − β) τxy dV, (6.27)
V V V

y si reemplazamos a σx por E (−z β 0 ) y a τxy por G (w0 − β), podemos escribir la (6.27) como:
Z Z
δε U = δ(−z β 0 ) E (−z β 0 ) dV + δ(w0 − β) G (w0 − β) dV. (6.28)
V V

Ahora integremos ambos términos respecto de z y de y, con lo cual la expresión (6.28)

queda como: Z Z
δε U = E I δβ 0 β 0 dx + G κ A δ(w0 − β) (w0 − β) dx. (6.29)
L L

Propongamos ahora las siguientes incógnitas con las que vamos a definir las dos funciones

Figura 6.4: Elemento de barra.

desplazamientos propuestas, como se ve en la figura 6.4, es decir, vamos a proponer una función
para los desplazamientos w:
w(x) = α0 + α1 x, (6.30)

Versión 2014 Preliminar - 63 -

6.3. La hipótesis de Timoshenko El Método de los Elementos Finitos

y otra para los giros β:

β(x) = α2 + α3 x. (6.31)
Así, la matriz de incógnitas desplazamiento queda de la siguiente forma:
 
w1
w2 
U = β1  ;
 (6.32)
β2

y si expresamos las funciones desplazamiento en forma matricial, tenemos:

 
    α0  
w(x) 1 x 0 0  α1  Φ w (x)
u(x) = =   = α. (6.33)
β(x) 0 0 1 x α2  Φ β (x)
α3

Hallemos la matriz A y consecuentemente A−1 :

    
w(0) 1 0 0 0 α0
w(Le ) 1 Le 0 0  α1 
 β(0)  = 0 0 1 0  α2  = A α
     ⇒ α = A−1 U. (6.34)
β(Le ) 0 0 1 Le α3

Con la matriz recién hallada hallemos ahora la expresión matricial de u(x):

 
w1
   x x 
Φ w (x) 1 − 0 0 w2 
A−1 U
 
u(x) = Le Le (6.35)
x x
Φ β (x) 0 0 1 − Le Le  β1 
β2

Si analizamos la matriz H(x), podemos diferenciar dos partes: la fila superior, que corres-
ponde a la función de w(x), la fila inferior, que corresponde a la función de β(x), es decir, que
las podemos identificar como Hw (x) y Hβ (x) , respectivamente. Por lo tanto, podemos hallar
Bw (x) y Bβ (x) en forma separada, que son las matrices derivadas de cada una de las filas de
H(x). Estas matrices resultan:

∂w(x) ∂Hw (x)

= U = Bw (x) U, (6.36)
∂x ∂x
y
∂β(x) ∂Hβ (x)
= U = Bβ (x) U. (6.37)
∂x ∂x
Una vez si reemplazamos en la expresión de la energía interna, obtenemos:
Z
T
δε U =δU Bβ (x)T E J Bβ (x) dx U+
Le
Z (6.38a)
+ δUT [Bw (x) − Hβ (x)]T G A κ [Bw (x) − Hβ (x)] dx U,
Le

que al agrupar convenientemente nos queda:

Z
T
δε U = δU Bβ (x)T E J Bβ (x) dx+
Le
Z  (6.38b)
T
+ [Bw (x) − Hβ (x)] G A κ [Bw (x) − Hβ (x)] dx U,
Le

- 64 - Preliminar Versión 2014

El Método de los Elementos Finitos 6. ELEMENTO DE VIGA

y si definimos que
Z
KF = Bβ (x)T E J Bβ (x) dx, y; (6.38c)
Le
Z
KC = [Bw (x) − Hβ (x)]T G κ A [Bw (x) − Hβ (x)] dx, (6.38d)
Le

donde KF es la rigidez a flexión y KC es la rigidez a corte, la expresión de la variación de la

energía interna de deformación queda:
δε U = δUT (KF + KC ) U = δUT K U. (6.39)
Finalmente, obtenemos el sistema
K U = R, (6.40)
que al resolverlo, define los valores de w1 ; β1 , w2 y β2 .
Sin embargo, este modelo tiene algunos inconvenientes. Uno de ellos es que si tomamos la
energía interna dividida la altura al cubo de una sección rectangular, obtendremos que a medida
que la altura es cada vez menor, la sección se vuelve más resistente, lo que es incompatible con
lo observado en la realidad. Este efecto se denomina «locking» o bloqueo. Hay varias maneras
de corregir esta situación:
Mediante el método de integración «falseada» o reducida;
Proponiendo una relación entre β y w, y;
Utilizando elementos con más nodos.
La integración «falseada» o reducida consiste en integrar numéricamente la matriz de
rigidez aplicando cuadratura de Gauss pero utilizando menos puntos que los necesarios para
obtener una integración «exacta». 2
Respecto a este caso, G. Prathap y y S. Mukherjee (véase [15]) demuestran que el elemento
de viga de tres nodos con las funciones de forma convencionales es inconsistente, y por lo tanto,
no es eficiente para representar las tensiones del elemento, y en consecuencia, sólo disminuye el
efecto de bloqueo, no lo elimina.
Al analizar la viga de Timoshenko de longitud L con un elemento de tres nodos con una
carga concentrada en un extremo y empotrada en el otro, obtienen las siguientes deformaciones
y tensiones expresadas como momentos y esfuerzo de corte:
   
PL 5
− 1 + ξ
 
ε̄ 2E J e+5
= , (6.41)
γ̄ P 5 P L2 3 ξ 2 −1
− G κ A − 12 E J 2 (e+5)
    "   #
σ̄ M̄ − P2L 1 + e+5 5
ξ
= = , (6.42)
τ̄ V̄ −P − 5 P e (3 ξ 2 − 1)
2 e+5
G κ A L2
donde e = 12 E J .
Cuando la viga se vuelve más esbelta, la altura h es cada vez más chica y en consecuencia,
el coeficiente e se hace cada vez más grande. El límite lo tendremos cuando e → ∞, y las
deformaciones y tensiones serán:
   PL 
ε̄ −2E J
= , (6.43)
γ̄ − G Pκ A
− P2L
     
σ̄ M̄
= = , (6.44)
τ̄ V̄ −P − 52P (3 ξ 2 − 1)
2
Para integrar numéricamente un polinomio de grado g en forma «exacta» se deben usar n puntos de Gauss-
Legendre que cumpla la siguiente relación: g ≤ 2 n − 1.

Versión 2014 Preliminar - 65 -

6.4. Determinación del coeficiente κ El Método de los Elementos Finitos

Podemos ver que a medida que la viga se vuelve esbelta, los términos lineales del momento
y las deformaciones por flexión tienden a anularse. Sin embargo, para el corte y la distorsión no
pasa la mismo. Mientras que el término cuadrático de la distorsión tiende a anularse también, el
término cuadrático del corte permanece, lo que introduce oscilaciones espúreas en las tensiones.
Si bien son menos importantes que para el caso de un elemento de dos nodos, no ha desaparecido
el efecto de bloqueo.
La explicación de este fenómeno la obtienen a través del concepto de proyección de la
solución, es decir, asumiendo que la solución por elementos finitos es una proyección ortogonal
de la solución real, como podemos ver en la figura 6.5.

Figura 6.5: Solución real y proyección ortogonal.

Así, cuando la solución aproximada no puede ser expresada como combinación lineal de
polinomios de Legendre, contiene oscilaciones espúreas y no puede representar correctamente las
tensiones del elemento.
Para obtener una solución correcta con un elemento de tres nodos, también debería usarse
la integración reducida.

6.4. Determinación del coeficiente κ

El coeficiente κ se determina a partir de la teoría de Jouravski Según esta teoría, el τxy
en una viga de sección rectangular está dado por la siguiente expresión:
"  #
Q h 2 2
τxy = −z , (6.45)
2 Iz 2
3
donde Iz = b12
h
.
τxy
La distorsión es γxy = G . Así, la energía interna de deformación por corte, por unidad
de longitud, se expresa como:
h
b
3 Q2
Z Z
1 2
UJ = τxy γxy dy dz = . (6.46)
2 0 −h 5bhG
2

La hipótesis de Timoshenko considera una distribución de tensiones constante, es decir

que:
Q
τxy =
, (6.47)
bh
con la misma relación entre tauxy y γxy . Por lo tanto la energía por unidad de longitud para este
caso se expresa de esta manera:
h
b
Q2
Z Z
1 2
UT = , τxy γxy dy dz = . (6.48)
2 0 −h 2bhG
2

- 66 - Preliminar Versión 2014

El Método de los Elementos Finitos 6. ELEMENTO DE VIGA

Para resolver nuestro problema, ambas energías internas de deformación deberían ser
iguales. Al hacer esto nos queda lo siguiente:

3 Q2 Q2
UJ = 6= = UT . (6.49)
5bhG 2bhG
Como podemos ver, no son iguales. Para que sean iguales hay que afectar a la energía in-
terna de deformación por Timoshenko por un coeficiente κ relacionado con la sección transversal.
Entonces tenemos:
3 Q2 Q2 3 1 1 5
= ⇒ = ⇒ κ= . (6.50)
5bhG 2κbhG 5 2κ 6
En consecuencia, la energía interna de deformación por corte para una sección rectangular
según la hipótesis de Timoshenko queda así:

Q2
UT = 5 , (6.51)
2 bhG
6

donde la sección 56 b h es una «sección equivalente» para mantener la energía interna de defor-
mación por la teoría de Jouravski.

Versión 2014 Preliminar - 67 -

El Método de los Elementos Finitos 7. ELEMENTO DE ESTADO PLANO

CAPÍTULO 7

ELEMENTO DE ESTADO PLANO

7.1. Introducción
en este capítulo vamos a analizar y estudiar los elementos de estado plano. Para ello
encararemos el armado de las matrices de rigidez de dichos elementos, tanto de estado plano de
tensión como de deformación.
En primer lugar, recordaremos las características de cada uno.

7.2. Estados planos

7.2.1. Estado plano de tensión
El estado plano de tensiones es aquel que tiene dos tensiones principales preponderantes
sobre la tercera, es decir, se cumple que:

σx 6= 0, σy 6= 0, τxy 6= 0, y σz = τxz = τzy ≈= 0. (7.1)

Este estado se da en estructuras que presentan la particularidad de tener dos dimensiones

preponderantes sobre la tercera y que las cargas están aplicadas en su contorno o en el plano
medio (por ejemplo, una viga de gran altura o viga pared, figura 7.1).

Figura 7.1: Elemento de estado plano de tensiones.

Versión 2014 Preliminar - 69 -

7.3. Elementos de estados planos El Método de los Elementos Finitos

Para este caso veamos las matrices que definen la ley de Hooke:
1 −ν −ν
    
εx 0 0 0 σx
 εy  −ν 1 −ν 0 0 0   σy 
    
  = 1 −ν −ν 1
 εz   0 0 0   σz 
γxy  E  0
  (7.2)
   0 0 2(1 + ν) 0 0  τxy 
 
γyz   0 0 0 0 2(1 + ν) 0  τyz 
γzx 0 0 0 0 0 2(1 + ν) τzx
Como hemos dicho, para un estado plano de tensiones se cumple que σz = τyz = τzx ≈ 0,
de manera que la matriz de la ecuación (7.2) se reduce a una matriz de 3x3. Si invertimos el
sistema, nos queda:     
σx 1 −ν 0 εx
 σy  = E −ν 1 0   εy  (7.3)
1 − ν2 1−ν
τxy 0 0 2 γxy

7.2.2. Estado plano de deformación

Este estado corresponde a aquél en el que dos deformaciones principales son preponde-
rantes sobre la tercera (al igual que en el caso anterior, si ε1 y ε2 6= 0, entonces ε3 ≈ 0). Como
ejemplo de este caso podemos tomar un cuerpo infinitamente largo y en el cual las cargas son
sólo función de x y y, como en la figura 7.2.

Figura 7.2: Elemento de estado plano de deformaciones.

Para este segundo caso, la ecuación constitutiva queda así:

    
σx 1−ν ν 0 εx
E
 σy  =  ν 1−ν 0   εy  (7.4)
(1 + ν)(1 − 2ν) 1−2ν
τxy 0 0 2 γxy

7.3. Elementos de estados planos

7.3.1. Elemento cuadrilátero de estado plano de tensión
Después de recordar las matrices constitutivas de ambos estados planos, nos vamos a
ocupar de desarrollar la matriz de rigidez para un elemento de estado plano de tensión. Para ello
vamos a tomar como ejemplo una ménsula corta de forma cuadrangular, con un empotramiento
en un borde y una carga concentrada en el extremo superior derecho.
Vamos a dividir a esta ménsula en cuatro elementos cuadrangulares 1 de estado plano, tal
como se ve en la figura 7.3, cada uno con ocho incógnitas desplazamiento, dos en cada extremo,
uno horizontal y el otro vertical, como se ve en el elemento genérico a la derecha.
1
Lo usual es empezar por el elemento más sencillo, el triangular. En este caso se verá a continuación del
elemento rectangular.

- 70 - Preliminar Versión 2014

El Método de los Elementos Finitos 7. ELEMENTO DE ESTADO PLANO

Figura 7.3: Estado plano. Elementos rectangulares.

Como en los elementos anteriores, proponemos las siguientes funciones desplazamiento

para un elemento genérico:
u(x, y) = α1 + α2 x + α3 y + α4 xy
(7.5)
v(x, y) = β1 + β2 x + β3 y + β4 xy
Observemos que la cantidad de coeficientes αi y βi es igual a la cantidad de incógnitas
uLi y vLi respectivamente. Con las funciones definidas, armemos la matriz α , que queda de la
siguiente manera:
    
uL1 1 x1 y1 x1 y1 0 0 0 0 α1
uL2  1 x2 y2 x2 y2 0 0 0 0  α2 
 
   
uL3  1 x3 y3 x3 y3 0 0 0 0   α3 
 
  
uL4  1 x4 y4 x4 y4 0 0 0 0 
Uhni =   α4  ,
 
 vL1  = 0 0 0 (7.6)
 
   0 1 x1 y1 x1 y1   
 β1 

 vL2  0 0 0 0 1 x2 y2 x2 y2    β2 
 
  
 vL3  0 0 0 0 1 x3 y3 x3 y3   β3 
vL4 0 0 0 0 1 x4 y4 x4 y4 β4

donde Ahni es la matriz que relaciona a Uhni y α hni .

Si utilizamos notación matricial, podemos escribir las funciones desplazamiento como:

uhni (x, y) = Φ hni (x, y) α hni , (7.7)

donde  
hni 1 x y xy 0 0 0 0
Φ (x, y) .
0 0 0 0 1 x y xy
Al igual que para los elementos vistos en los capítulos precedentes, podemos expresar
α hni en función de Uhni :

Uhni = Ahni · α hni ⇒ α hni = [Ahni ]−1 · Uhni (7.8)

Versión 2014 Preliminar - 71 -

7.3. Elementos de estados planos El Método de los Elementos Finitos

Si reemplazamos la ecuación (7.8) en la ecuación (7.7), obtenemos la siguiente expresión

para el un elemento genérico:

uhni (x, y) = Φ hni (x, y) · [Ahni ]−1 · Uhni = Hhni (x, y) · Uhni , (7.9)

con
Hhni (x, y) = Φ hni (x, y) · [Ahni ]−1 ,
que es la matriz de la función interpolante de los desplazamientos Uhni . Esta matriz la podemos
expresar también de esta forma:
 
h(x, y) [0]
Hhni (x, y) = , (7.10)
[0] h(x, y)
donde  
h(x, y) = h1 (x, y) h2 (x, y) h3 (x, y) h4 (x, y) ,
y  
[0] = 0 0 0 0 .
Si definimos las coordenadas de cada uno de los puntos de esta manera, P1 = (−1, −1),
P2 = (1, −1), P3 = (1, 1) y P4 = (−1, 1), para un elemento genérico cualquiera 2 , cada una de
estas hi (x, y) quedan definidas así:
1
h1 (x, y) = (1 − x)(1 − y)
4
1
h2 (x, y) = (1 + x)(1 − y)
4 (7.11)
1
h3 (x, y) = (1 + x)(1 + y)
4
1
h4 (x, y) = (1 − x)(1 + y)
4
Ahora armemos las matrices que relacionan los desplazamientos u(x, y) y v(x, y), con las
deformaciones εx , εy y γxy :
  ∂  
∂

εx ∂x 0 0
 ∂x ∂  hni
 
∂  u(x, y) hni hni hni
 0 ∂y 
 εy  =  =  0 ∂y  H (x, y) U = B (x, y) U , (7.12)
∂ ∂ v(x, y) ∂ ∂
γxy ∂y ∂x ∂y ∂x

con  
∂
∂x 0
∂ 
Bhni (x, y) =  0 Hhni (x, y).

∂y 
∂ ∂
∂y ∂x

Con este último paso hemos expresado las deformaciones en función de las incógnitas
desplazamiento del elemento. Nos falta aún relacionar las tensiones en el elemento con las defor-
maciones. Esto es muy sencillo porque sabemos que:

σ hni = E ε hni = E Bhni (x, y) Uhni , (7.13)

con  hni 
σx
 hni 
σ hni = σy  ,
hni
τxy
2
Esta forma de definir las coordenadas de un elemento genérico se verá más adelante en la denominada
formulación isoparamétrica.

- 72 - Preliminar Versión 2014

El Método de los Elementos Finitos 7. ELEMENTO DE ESTADO PLANO

y E es una matriz que depende del estado plano en análisis (de tensión o de deformación).
Como ya hemos hallado todos los términos que componen la expresión (7.13), sólo nos
resta plantear los términos de la variación de la energía interna de deformación y la variación del
trabajo de las fuerzas exteriores.
La variación de la energía interna de deformación está dada por:
ZZZ
hni
δU = δεεhni σ hni dz dy dx. (7.14)

Para integrar esta expresión asumiremos primero que el espesor del elemento es constante
y vale h es decir, z ∈ (0, h). En segundo lugar, reemplacemos la integral en z por la sumatoria
de los hhni . Con estas simplificaciones el término de la energía interna nos queda:
X ZZ
hni hni
δU U = h δεεhni σ hni dy dx
n
ZZ (7.15)
hni T
X
hni
= δU h Bhni (x, y)T E Bhni (x, y) dy dx Uhni ,
n

que podemos resumir en:

T
δU U hni = δUhni Khni Uhni . (7.16)
En el caso de la variación del trabajo de las fuerzas exteriores tenemos que:
T
δU W hni = δUhni Rhni , (7.17)

donde Rhni son las cargas nodales equivalentes para cada elemento. Igualando las ecuacio-
nes (7.16) y (7.17) obtenemos:
T T T
 
δUhni Khni Uhni = δUhni Rhni ⇒ δUhni Khni Uhni − Rhni = 0, (7.18)

que resulta en la siguiente ecuación:

Khni Uhni = Rhni , (7.19)

cuya solución son las incógnitas Uhni .

7.3.2. Elemento triangular de estado plano

En el punto anterior hemos desarrollado el elemento de estado plano de forma cuadrangu-
lar. Una de las características más salientes de ese elemento es que las funciones de interpolación
de los desplazamientos son de segundo orden (aparece el término x y, que es un término rectan-
gular), por lo tanto se «pierde» la sencillez de la interpolación lineal. De modo que la resolución
de ese tipo de elemento requiere más recursos, tanto matemáticos como de computación.
Existe, sin embargo, la posibilidad de modelar el sistema a resolver mediante otro elemen-
to, mucho más sencillo, que es el elemento triangular. En este caso volvemos a tener funciones
de interpolación lineales. Par ello, utilicemos el mismo ejemplo anterior pero en esta caso ar-
memos una discretización con elementos triangulares. Vemos que la cantidad de elementos para
representar la malla es el doble que en el caso de elementos cuadrangulares, como se ve en la
figura 7.4.
En este caso, vamos a definir nuevas funciones de interpolación para los desplazamientos
u(x, y) y v(x, y):

u(x, y) = α1 + α2 x + α3 y,
(7.20)
v(x, y) = β1 + β2 x + β3 y,

Versión 2014 Preliminar - 73 -

7.3. Elementos de estados planos El Método de los Elementos Finitos

Figura 7.4: Estado plano. Elementos triangulares.

que sólo tienen términos lineales. Podemos expresar esto de forma matricial:

uhni (x, y) = Φ hni (x, y) α hni , (7.21)

con  
α1
α2 
   
1 x y 0 0 0 α3 
Φhni (x, y) y αhni =
 β1  .

0 0 0 1 x y  
 β2 
β3
Como ya vimos para el caso anterior, tenemos que:

Uhni = Ahni · α hni ⇒ α hni = [Ahni ]−1 · Uhni , (7.22)

con    
1 x1 y1 0 0 0 u1
1 x2 y2 0 0 0 u2 
   
1 x3 y3 0 0 0 u3 
Ahni =
0 0 0
 y Uhni =
 v1  .

 1 x 1 y1 
  
0 0 0 1 x 2 y2   v2 
0 0 0 1 x 3 y3 v3
Finalmente, podemos escribir la función uhni (x, y) como:

uhni (x, y) = Φ hni (x, y) [Ahni ]−1 · Uhni = Hhni (x, y) · Uhni , (7.23)

con Hhni (x, y) = Φ hni (x, y) [Ahni ]−1 .

Si expresamos las funciones de u(x;y) y v(x;y) en forma independiente (no matricial),
obtenemos lo siguiente:

u(x, y) = N1 (x, y) u1 + N2 (x, y) u2 + N3 (x, y) u3 ,

(7.24)
v(x, y) = N1 (x, y) v1 + N2 (x, y) v2 + N3 (x, y) v3 .

Las funciones Ni (x, y) son:

1
N1 (x, y) = [(x2 y3 − x3 y2 ) + (y2 − y3 ) x + (x3 − x2 ) y]
2 AT
1
N2 (x, y) = [(x3 y1 − x1 y3 ) + (y3 − y1 ) x + (x1 − x3 ) y] (7.25)
2 AT
1
N3 (x, y) = [(x1 y2 − x2 y1 ) + (y1 − y2 ) x + (x2 − x1 ) y]
2 AT

- 74 - Preliminar Versión 2014

El Método de los Elementos Finitos 7. ELEMENTO DE ESTADO PLANO

donde AT es el área del triángulo.

Existe otra forma de escribir lo mismo:

A1 A2 A3
N1 = , N2 = y N3 = , (7.26)
AT AT AT

donde A1 , A2 y A3 son funciones de x, y que están representadas en la figura 7.5.

Figura 7.5: Funciones A1 , A2 y A3 .

Ahora reemplacemos las (7.26) en la ecuación (7.24). Nos queda:

1
u(x, y) = (A1 u1 + A2 u2 + A3 u3 ) ,
AT
(7.27)
1
v(x, y) = (A1 v1 + A2 v2 + A3 v3 ) .
AT

Las funciones desplazamientos que obtenemos para este caso resultan ser funciones linea-
les, que son mucho más sencillas de manejar que las funciones interpolantes que obtenemos con
al elemento cuadrilátero, tal como se indicó al principio. Como las tensiones surgen de derivar
esas funciones desplazamiento, que son constantes, al elemento triangular se lo conoce también
como elemento de tensión o deformación constante.

7.4. Simetría de revolución

Un caso particular de los estados planos puede encontrarse con las estructuras o sólidos
con simetría de revolución, como pueden ser una cúpula o un tanque (especialmente en los
sectores alejados de los bordes).
En estos casos se puede aplicar un análisis por elemento finitos utilizando elementos de
estados planos, introduciendo la relación constitutiva para estados con simetría de revolución.
Al igual que en los casos anteriores, se proponen dos funciones para los desplazamientos.
Suponiendo un estado plano de deformaciones y un elemento triangular, las funciones son las
definidas en (7.20), y expresándolas en forma matricial, nos queda como en (7.21).
Ahora, definamos las deformaciones de este elemento en función de los desplazamientos.
Estas deformaciones están dadas por:
 ∂  
∂

0 0

εx ∂x ∂x
∂   ∂ 
 =0
 εy   0
∂y  u(x, y)
=∂ ∂y 
Hhni (x, y) Uhni = Bhni (x, y) Uhni , (7.28)
 
γxy   ∂ ∂  v(x, y) ∂ 
∂y ∂x   ∂y ∂x 
εz 1 1
r 0 r 0

Versión 2014 Preliminar - 75 -

7.4. Simetría de revolución El Método de los Elementos Finitos

La ecuación constitutiva para un estado con simetría de revolución, es:

    
σx 1−ν ν 0 ν εx
 σy 
 = E(1 − ν)  ν
 1−ν 0 ν   εy 
 
.
τxy  (1 + ν)(1 − 2ν)  0 1−2ν (7.29)
0 2 0  γxy 
σz ν ν 0 1−ν εz

Con las ecuaciones (7.28) y (7.29) podemos obtener la variación de la energía interna de
deformación y luego aplicar el Teorema de los Desplazamientos Virtuales. La variación de
la energía interna de deformación por unidad de volumen es:
Z ZZZ
T
δU U = δεε σ dV = δεεT σ r dθ dy dx. (7.30)
V

Si definimos que:
 
1−ν ν 0 ν
E(1 − ν)  ν
 1−ν 0 ν 
C= 1−2ν
,
(1 + ν)(1 − 2ν)  0 0 2 0 
ν ν 0 1−ν

podemos expresar la variación de energía interna por unidad de radián:

ZZ ZZ
hni T T
X
δU U = T
δεε σ r dy dx = δU Bhni C Bhni r dy dx Uhni . (7.31)
m

La variación del trabajo de las fuerzas exteriores por unidad de radián está dada por:

T Rhni
δU W = δUhni , (7.32)
2π
donde Rhni son las cargas nodales. En caso de tener una carga distribuida (ver figura 7.6), hay
que obtener las cargas nodales equivalentes.

Figura 7.6: Cargas distribuidas.

Para ello, tomemos las funciones unas funciones adicionales, u(s) y v(s):
hni
 y hni hni y
u(y) = u1 1 − + u2 0 + u3
L L (7.33)
hni
 y hni hni y
v(y) = v1 1− + v2 0 + v3 ,
L L

- 76 - Preliminar Versión 2014

El Método de los Elementos Finitos 7. ELEMENTO DE ESTADO PLANO

que podemos resumir en una expresión matricial:

U(y) = Hhni (y) Uhni (7.34)

Con estas funciones desplazamiento, hallemos la variación del trabajo de las fuerza exte-
riores:
Z Z
1 1 T T
δU W = δu(y)p(y) dy = δUhni Hhni (y) p(y) dy
2π 2π
T (7.35)
δUhni
Z
hni T
δU W = H (y) p(y) dy.
2π
Para resolver el problema, basta con aplicar el Teorema de los Desplazamientos
Virtuales, igualando las expresiones (7.31) y (7.32) o (7.31) y (7.35). Así obtenemos los despla-
zamientos o incógnitas.

Versión 2014 Preliminar - 77 -

 CAPÍTULO 8

ELEMENTO DE PLACA

8.1. Introducción
En este capítulo nos ocuparemos de los elementos de placa plana. Primeramente vamos
a recordar la definición de una placa:

Una placa es el lugar geométrico de los puntos comprendidos entre dos superficies que
equidistan de un plano, llamado plano medio. Dichas superficies se encuentran muy
próximas entre sí a una distancia llamada espesor y que mucho menor que las otras
dimensiones que se pueden medir sobre el plano medio.Las cargas son principalmente
en la dirección normal al plano medio.

Se pueden considerar los siguientes tipos:

1. Placas planas delgadas: Una placa se considera delgada cuando su espesor es menor que
un décimo de la mínima dimensión en su plano medio. En este caso se parte de la hipótesis
de Kirchhoff.

2. Placas flexibles o muy delgadas: Se trata de placas cuyas deflexiones son mayores a
0,5 h, siendo h el espesor de la misma.

3. Membranas: Son placas que no tienen rigidez a flexión y sólo resisten las cargas normales
al plano medio con esfuerzos normales y tangenciales.

4. Placas gruesas: Cuando el espesor es mayor a un décimo de la mínima dimensión en su

plano medio. En este caso se utiliza la hipótesis de Mindlin.

Analizaremos primero la placa sometida a solicitaciones de flexión solamente y hallaremos

la matriz de rigidez, aplicando la hipótesis de Kirchhoff. Plantearemos como elemento uno de
forma rectangular y obtendremos como solución una función para determinar los desplazamientos
en el interior de la placa a partir de los desplazamientos nodales.
En segundo término, desarrollaremos un elemento de placa que tenga en cuenta la de-
formación por corte, que corresponde al caso 4, es decir, tomaremos en cuenta la hipótesis de
Mindlin, y nuevamente, hallaremos una función desplazamiento que tenga en cuenta dichas de-
formaciones por corte. En este caso se analizará la convergencia de los resultados y los métodos
para asegurar dicha convergencia.

79
8.2. Placa plana con Hipótesis de Kirchhoff El Método de los Elementos Finitos

8.2. Placa plana con Hipótesis de Kirchhoff

8.2.1. Obtención de las deformaciones

Empezaremos el análisis del elemento de placa partiendo de la hipótesis de Kirchhoff: las

secciones planas y normales al plano medio se mantienen planas y perpendiculares a dicho plano
una vez que la pieza se ha deformado. A partir de esta suposición hallemos las deformadas de la
placa.

Figura 8.1: Hipótesis de Kirchhoff.

De la figura 8.1 podemos ver que los desplazamientos están descriptos por las expresiones:

ū0 = u0 ex + v0 ey + w0 ez
(8.1)
ū = u ex + v ey + w ez .

Por hipótesis, los desplazamientos son pequeños y por lo tanto se puede asumir que:

∂w
α ≈ sen α ≈ tg α ≈ − (8.2)
∂x
∂w
β ≈ sen β ≈ tg β ≈ − (8.3)
∂y

- 80 - Preliminar Versión 2014

El Método de los Elementos Finitos 8. ELEMENTO DE PLACA

Con (8.2) y (8.3) podemos expresar los desplazamientos así:

∂w
u = u0 − z ,
∂x
∂w (8.4)
v = v0 − z ,
∂y
w = w0 .

Con las expresiones de los desplazamientos vamos a obtener las deformaciones de la placa,
asumiendo que u0 , v0 y w0 son constantes:

∂2w
εx = −z ,
∂x2    ∂2w 
2 εx ∂x2
∂ w  2 
εy = −z 2 , ⇒  εy  = −z  ∂∂yw2  . (8.5)
∂y ∂2w
2
γxy 2 ∂x
∂ w ∂y
γxy = −2 z .
∂x ∂y

8.2.2. Ecuaciones de equivalencia

Para introducir las características de los materiales en el desarrollo del elemento, anali-
cemos el caso de un material homogéneo, isótropo y elástico-lineal. Es válida, entonces, la Ley
de Hooke:     
σx 1 ν 0 εx
 σy  = E ν 1 0   εy  (8.6)
1 − ν2 1−ν
τxy 0 0 2 γ xy

Analicemos ahora las ecuaciones de equivalencia entre tensiones y solicitaciones. En nues-

tro caso vamos a trabajar con los momentos Mx , My y Mxy . Podemos expresar estos momentos
como:
Z h
2
Mx = σx z dz,
−h
2
Z h
2
My = σy z dz, (8.7)
−h
2
Z h
2
Mxy = τxy z dz,
−h
2

que matricialmente escribimos así:

Z h
2
M= σ z dz, (8.8)
−h
2

con 
Mx
M =  My  ,
Mxy
y
 
σx
σ =  σy  .
τxy

Versión 2014 Preliminar - 81 -

8.2. Placa plana con Hipótesis de Kirchhoff El Método de los Elementos Finitos

Ahora reemplacemos la (8.5) en la (8.6):

     ∂2w 
σx 1 ν 0 ∂x2
 σy  = E ν 1  ∂2w 
0 (−z)  ∂y2  .
 (8.9)
1 − ν2 1−ν ∂2w
τxy 0 0 2 2 ∂x ∂y

Si definimos  
1 ν 0
E 
C= ν 1 0  (8.10)
1 − ν2 1−ν
0 0 2
y
∂2w
 
∂x2
 ∂2w 
w00 =  ∂y2  = B(x; y) U, (8.11)
∂2w
2 ∂x ∂y
entonces la ecuación (8.8) queda así:
Z h
2
M= C B(x, y)(−z 2 ) dz U, (8.12)
−h
2

que al integrar nos deja:

h3
M=− C B(x, y) U, (8.13)
12

8.2.3. Planteo del Teorema de los Desplazamientos Virtuales

Partamos ahora de la ecuación diferencial de una placa con la hipótesis de Kirchhoff:
M00 + p = 0, (8.14)
Si aplicamos el Teorema de los Desplazamientos Virtuales y operamos algebraica-
mente, nos queda: Z Z
T
− δw00 M dA = δwT p dA, (8.15)
A A
que en notación matricial queda:
Z  3 Z
T h
− T
δU B(x, y) − C B(x, y) U dA = δUT
H(x, y)T p dA, (8.16)
A 12 A
con
w(x, y) = H(x, y) U,
y, por lo tanto
δw(x, y) = H(x, y) δU.
Si definimos
h3
Z
KL = B(x, y)T C B(x, y) dA (8.17)
12 A
y Z
R= H(x, y)T p dA, (8.18)
A
podemos reescribir la ecuación (8.16) así:
δUT (KL U − R) = 0, (8.19)
que termina en el sistema de ecuaciones lineales ya conocido:
KL U = R. (8.20)

- 82 - Preliminar Versión 2014

El Método de los Elementos Finitos 8. ELEMENTO DE PLACA

8.3. Placa plana con Hipótesis de Mindlin

En el punto, tomamos cómo hipótesis válidas que εz = γxz = γyz = 0. Esto nos permitía
resolver nuestra ecuación diferencial tomando cómo incógnitas nada más que los desplazamientos
normales a la superficie media del elemento, y los giros relativos. Al igual que en el caso de un
elemento de viga, también en placas planas podemos considerar las deformaciones introducidas
por el esfuerzo de corte.
Para ello vamos considerar la denominada hipótesis de Mindlin, que introducen las dis-
torsiones γxz 6= 0 y γyz 6= 0, quedando solamente que εz ∼
= 0.
De nuevo, vamos a proponer funciones para los desplazamientos y para los giros en cada
nodo. Así tendremos las siguientes funciones para cada caso:

Desplazamientos
X
w(x, y) = wi hi (x, y), (8.21)

y giros
X
βx (x, y) = βxi hi (x, y), (8.22)
X
βy (x, y) = βyi hi (x, y). (8.23)

Nuestro siguiente paso es definir las funciones de las deformaciones, tanto las debidas a
la flexión como a las debidas al corte, que tendrán la siguiente forma:

Flexión
   ∂βx 
εx ∂x
 ∂βy
ε =  εy  = −z   = −z Bf (x, y) U, (8.24)

∂y
∂βy
γxy
∂y+ ∂βx ∂x

Corte
  " ∂w(x,y) #
γyz ∂y − β y
γ= = ∂w(x,y) = Bc (x, y) U. (8.25)
γzx − βx
∂x

De acuerdo con la Ley de Hooke, tendremos las tensiones:

Flexión
 ∂βx
   
σx 1 ν 0 ∂x
E ∂βy
σ =  σy  = −z ν 1 0   = −z E Bf /x, y) U (8.26)
 
∂y
(1 − ν 2 ) 1−ν ∂βy ∂βx
τxy 0 0 2 +
∂y ∂x

Corte
" #
  ∂w(x,y)
σyz E ∂y − β y
τ = = = G Bc (x, y)U. (8.27)
τzx 2(1 + ν) ∂w(x,y) − βx .
∂x

Si planteamos el término de la variación de la energía interna de deformación:

Z Z
T
δU U = δεε σ dV + δγγ T τ dV, (8.28)
V V

Versión 2014 Preliminar - 83 -

8.4. Elementos isoparamétricos El Método de los Elementos Finitos

y reemplazamos las ecuaciones (8.24), (8.25), (8.26) y (8.27) en la (8.28), obtenemos:

Z Z 
T 2 T
δU U = δU Bc (x, y) E Bf (x, y) z dV + Bc (x, y) G Bc (x, y) dV U, (8.29)
V V

que puede ser reducida a:

δU U = δU [Kf + Kc ] U, (8.30)
con
 
1 ν 0
h3
Z
E T 
Kf = Bf (x, y) ν 1 0  Bf (x, y) dA, (8.31a)
12(1 − ν 2 ) A 1−ν
0 0 2
Z
Kc = G h κ Bc (x, y)T Bc (x, y) dA. (8.31b)
A

En esta expresión aparece también el coeficiente κ, como hemos visto en el capítulo de

elemento de viga.
En ambos casos, nuestra formulación de los elementos se ha hecho con funciones despla-
zamientos y giros en coordenadas x e y. Este procedimiento puede resultar engorroso cuando se
trabaja con elementos de diferentes tamaños o tipos. En esos casos se puede realizar el mismo
desarrollo mediante los denominados elementos isoparamétricos.

8.4. Elementos isoparamétricos

Hemos visto en que el desarrollo de los elementos de placa (y también de elementos de
estado plano) es notoria la influencia de las coordenadas, sobre todo cuando hay que obtener
la integración de las funciones interpolantes. Si bien la integración se realiza mayoritariamente
mediante métodos numéricos, esto no evita las complicaciones derivadas de las coordenadas de
los nodos. Por otro lado, aproximar una estructura plana mediante elementos rectangulares es
mucho mejor que hacerlo con elementos triangulares, pero con el inconveniente que los primeros
no se adaptan tan fácilmente a contornos irregulares, cosa que sí se puede hacer con los segundos.
Es por eso que surge la idea de definir un elemento rectangular peor con la característica
de independizarse de las coordenadas globales (o generales) del modelo. Para ello, se define un
sistema de coordenadas local, por ejemplo, los ejes r y s, que suelen ubicarse en el baricentro del
rectángulo. Todas las coordenadas globales de ese elemento se relacionan con el sistema local. De
acuerdo con la cantidad de nodos que se utilicen para definir el elemento rectangular en cuestión,
se definirán las funciones interpolantes para representar las coordenadas de los nodos y luego de
las incógnitas desplazamiento.
Básicamente, existen tres tipos de elementos derivados de este concepto:

Los elementos subparamétricos.

Son aquellos en los cuales la cantidad de incógnitas es menor a la cantidad de coordenadas
geométricas necesarias para definir el elemento.

Los elementos hiperparamétricos.

Son aquellos en los cuales la cantidad de incógnitas es mayor a la cantidad de coordenadas
geométricas necesarias para definir el elemento.

Los elementos isoparamétricos.

Son aquellos en los cuales la cantidad de incógnitas es igual a la cantidad de coordenadas
geométricas necesarias para definir el elemento.

- 84 - Preliminar Versión 2014

El Método de los Elementos Finitos 8. ELEMENTO DE PLACA

Nos ocuparemos, entonces de los elementos isoparamétricos. Una correcta definición de un

elemento isoparamétrico es la siguiente: es todo elemento en el cual la geometría y las incógnitas
que los describen están formulados con el mismo tipo de función polinómica, también llamadas
funciones de forma.
Veamos la figura 8.2:

Figura 8.2: Elemento isoparamétrico rectangular.

Las siguientes funciones de forma para este caso son:

1
h1 (s, s) = (1 − r)(1 − s)
4
1
h2 (s, s) = (1 + r)(1 − s)
4 (8.32)
1
h3 (s, s) = (1 + r)(1 + s)
4
1
h4 (s, s) = (1 − r)(1 + s)
4
Como vemos, las cuatro funciones (h1 a h4 ) están expresadas en términos de r y s, en
vez de expresarlas en función de x e y. Con estas cuatro funciones hi vamos a expresar todas
las funciones que son necesarias para describir nuestro elemento. Primero vamos a definir las
funciones que definen al elemento en forma geométrica:
X
X(r, s) = hi (r, s) xi , (8.33)
i
X
Y (r, s) = hi (r, s) yi , (8.34)
i

donde xi y yi son las coordenadas de los puntos 1 a 4.

De igual forma, definimos las funciones para los desplazamientos y giros:
X
w(r, s) = hi (r, s) wi , (8.35)
i
X
βx (r, s) = hi (r, s) βxi , (8.36)
i
X
βy (r, s) = hi (r, s) βyi . (8.37)
i

Recordemos primero el concepto de la matriz jacobiana (Jacobiano). Esta matriz nos

sirve para transformar las derivadas de una función respecto a dos variables dadas (por ejemplo

Versión 2014 Preliminar - 85 -

8.4. Elementos isoparamétricos El Método de los Elementos Finitos

x y y) en derivadas respecto a otras variables (por ejemplo r y s) que a su vez son funciones de
x y y. Es decir:
    
∂F (x,y) ∂x ∂y  ∂F (x,y)
 ∂r   ∂r ∂r   ∂x 
 =  , (8.38)
∂F (x,y) ∂x ∂y ∂F (x,y)
∂s ∂s ∂s ∂y

con
 ∂x ∂y 
∂r ∂r
J= .
∂x ∂y
∂s ∂s

La transformación inversa queda de la siguiente forma:

   
∂F (r,s) ∂F (r,s)
∂x
−1  ∂r 
=J  , (8.39)
 

∂F (r,s) ∂F (r,s)
∂y ∂s

Apliquemos estos conceptos a lo visto en el punto 8.3. En primer lugar hallemos la

derivada de los desplazamientos:

   
w1

∂w(r,s) ∂h1 (r,s) ∂h2 (r,s) ∂h3 (r,s) ∂h4 (r,s)
 ∂r   ∂r ∂r ∂r ∂r
w2 
 
=  . (8.40)
w3

∂w(r,s) ∂h1 (r,s) ∂h2 (r,s) ∂h3 (r,s) ∂h4 (r,s)
∂s ∂s ∂s ∂s ∂s w4

∂w(r,s) ∂w(r,s)
Para obtener ∂x y ∂y , debemos premultiplicar el miembro de la derecha por J−1 .
Nos queda:
   
w1

∂w(r,s) ∂h1 (r,s) ∂h2 (r,s) ∂h3 (r,s) ∂h4 (r,s)
∂x
−1  ∂r ∂r ∂r ∂r
w2 
 
=J   . (8.41)
 
w3

∂w(r,s) ∂h1 (r,s) ∂h2 (r,s) ∂h3 (r,s) ∂h4 (r,s)
∂y ∂s ∂s ∂s ∂s w4

En forma análoga tendremos:

   
β x1

∂βx (r,s) ∂h1 (r,s) ∂h2 (r,s) ∂h3 (r,s) ∂h4 (r,s)
∂x
−1  ∂r ∂r ∂r ∂r
βx2 
 
=J   , (8.42)
 
β x3

∂βx (r,s) ∂h1 (r,s) ∂h2 (r,s) ∂h3 (r,s) ∂h4 (r,s)
∂y ∂s ∂s ∂s ∂s β x4
    β 
∂βy (r,s) ∂h1 (r,s) ∂h2 (r,s) ∂h3 (r,s) ∂h4 (r,s) y1
∂x ∂r ∂r ∂r ∂r β
−1   y2 

=J  βy3  . (8.43)
 
 
∂βy (r,s) ∂h1 (r,s) ∂h2 (r,s) ∂h3 (r,s) ∂h4 (r,s)
∂y ∂s ∂s ∂s ∂s βy4

Para simplificar las expresiones definamos que:

 
∂h1 (r,s) ∂h2 (r,s) ∂h3 (r,s) ∂h4 (r,s)
 ∂r ∂r ∂r ∂r
M = J−1  ,

∂h1 (r,s) ∂h2 (r,s) ∂h3 (r,s) ∂h4 (r,s)
∂s ∂s ∂s ∂s

donde la matriz M es función de las variables r y s.

- 86 - Preliminar Versión 2014

El Método de los Elementos Finitos 8. ELEMENTO DE PLACA

Armemos en primer lugar cada uno de los términos correspondientes a las deformaciones.
Nos queda lo siguiente:
 
β x1
∂βx (r, s)   βx2 
= M11 M12 M13 M14  βx3  ,
 (8.44)
∂x
β x4
 
βy1
∂βy (r, s)   βy2 
= M21 M22 M23 M24  βy3  ,
 (8.45)
∂y
βy4
 
β x1
βx2 
 
βx3 
 
∂βx (r, s) ∂βy (r, s)   βx4 
+ = M21 M22 M23 M24 M11 M12 M13 M14  βy  .
 (8.46)
∂y ∂x  1
βy 
 2
βy 
3
βy4
Definamos ahora la siguiente matriz:
 
0 0 0 0 M11 M12 M13 M14 0 0 0 0
Bf (r, s) =  0 0 0 0 0 0 0 0 M21 M22 M23 M24  . (8.47)
0 0 0 0 M21 M22 M23 M24 M11 M12 M13 M14
Con ayuda de la (8.47), podemos escribir las ecuaciones (8.44), (8.45) y (8.46) en forma
matricial:
 
w1
 .. 
 . 
 
 w4 
∂βx (r,s)
   
βx   
∂x  1 εx
∂βy (r,s)  .. 
 = Bf (r, s)  .  ⇒ εy  = −z Bf (r, s) UL . (8.48)
  
 ∂y
∂βx (r,s) ∂βy (r,s)  
γ
+ ∂x βx4  xy
 
∂y
βy 
 1
 .. 
 . 
βy4
Con la definición del vector ε podemos definir también las tensiones:
   
σx 1 ν 0
 σy  = E ν 1 0  (−z) Bf (r, s) UL . (8.49)
1 − ν2 1−ν
τxy 0 0 2

Por último, si simplificamos la notación definiendo que:

 
1 ν 0
C = ν 1 0 ,
1−ν
0 0 2

podemos plantear la variación de la energía interna de deformación debida a la flexión:

Z Z Z h
2
δUf = δUL T Bf (r, s)T C Bf (r, s) (−z)2 dz |J| dr ds UL , (8.50)
s r −h
2

Versión 2014 Preliminar - 87 -

8.4. Elementos isoparamétricos El Método de los Elementos Finitos

donde Z Z Z h
2
T
Kf = Bf (r, s) C Bf (r, s) (−z)2 dz |J| dr ds, (8.51)
s r −h
2

es decir,
h3
Z Z
Kf = Bf (r, s)T C Bf (r, s) |J| dr ds, (8.52)
12 s r
que resulta ser la (8.31a).
Al hacer lo mismo pero para el caso de la rigidez por corte, con un desarrollo análogo
obtenemos la matriz Kc de la expresión (8.31b).

- 88 - Preliminar Versión 2014

 CAPÍTULO 9

TEOREMAS VARIACIONALES MIXTOS

9.1. Introducción
Hasta ahora nos hemos ocupado de resolver nuestros problemas mediante la aplicación de
métodos que utilizaban incógnitas de un solo tipo, esto es, en el caso de los sistemas estructurales,
proponiendo funciones que determinaran los desplazamientos (Teorema de los Desplazamientos
Virtuales o T.D.V.), o las tensiones (Teorema de las Tensiones Virtuales o T.T.V.). En los ca-
pítulos 3 a 8 nos ocupamos de ver diferentes métodos y de definir una cantidad de elementos
básicos para resolver algunos de los problemas más comunes en la ingeniería civil, fundamental-
mente en el campo de los sistemas estructurales. Así vimos que si aplicamos el T.D.V. en un
sistema estático para hallar la función desplazamiento y se pretende conocer el estado tensional
en una sección determinada, por ejemplo, el momento flector máximo en una viga, deberemos
hallar dicho estado a partir de la función desplazamiento, en nuestro caso, deberemos partir de la
relación M (x) = −E I w00 (x), lo que significa una pérdida de precisión en la solución aproximada
introducida por la necesidad de derivar nuestra función desplazamiento dos veces respecto a x.
Esta pérdida de precisión en muchos casos es significativa, toda vez que depende de la función
desplazamiento elegida, y así se evidencia en varios ejemplos prácticos vistos como aplicación del
método de Ritz. También hemos visto que mediante la utilización del método de los elementos
finitos es posible reducir el error en las aproximaciones de los resultados a través de discretizar
nuestro sistema de acuerdo con las necesidades del problema, pero esto requiere trabajar en for-
ma iterativa. Igualmente estos resultados también están vinculados al tipo de función propuesta
como solución para elemento considerado.
Hay también una cantidad de problemas en la ingeniería que no se ajustan al modelo
que hemos ido desarrollando hasta ahora, y que por lo tanto, no es posible resolver mediante
elementos finitos proponiendo funciones solución para un solo tipo de incógnita. Un caso bien
claro en este sentido son los problemas derivados de la mecánica de los fluidos.
Por ello, veamos a continuación otros métodos que partiendo también de principios va-
riacionales, nos permiten resolver algunos de nuestros problemas con menores errores de aproxi-
mación.

9.2. Teorema de Hellinger-Reissner

Para el estudio de estos teoremas partiremos primeramente de un caso sencillo y clásico
de sistemas estructurales que ya hemos resuelto mediante el método de Ritz. Tomemos el ejemplo
del punto 4.3.1. Se trata de una viga simplemente apoyada con una carga concentrada en el centro
de la luz.

89
9.2. Teorema de Hellinger-Reissner El Método de los Elementos Finitos

Figura 9.1: Viga sometida a carga concentrada.

La resolución de este problema mediante el método de Ritz arroja como función aproxi-
mada de los momentos flectores la siguiente expresión: M (x) = P8L
Según vimos en el capítulo 4, la función de momentos es una función lineal, mientras que
la obtenida por aplicación del método de Ritz es una función constante. Esto es debido a que
la función desplazamiento propuesta es una parábola de segundo grado, que al ser derivada dos
veces, se convierte en una constante.
Para resolver este problema y reducir el error para los momentos flectores, vamos a aplicar
el concepto de los teoremas variacionales mixtos con el objetivo de mejorar nuestra aproximación
de la función mencionada. Empezaremos de forma similar que en el caso del método de Ritz.
Propongamos una función para los desplazamientos, en este caso, una función lineal definida de
la siguiente forma:

x
w1 2 L
 si 0 ≤ x ≤ L2 ,
w(x) = (9.1)
 x

L
w1 2 1 − L si 2 < x ≤ L,


donde w1 es nuestra incógnita desplazamiento (incógnita cinemática) ubicada en el centro de la

luz de la viga en concordancia con la ubicación de la carga concentrada P . Esta función está
representada en la figura 9.2:

Figura 9.2: Desplazamiento propuesto.

Como podemos ver, la función que hemos propuesto cumple con las condiciones de borde
cinemáticas:

w(0) = 0, (9.2a)

w(L) = 0. (9.2b)

Lo nuevo es que ahora vamos a proponer una solución aproximada para los momentos
flectores. Para ello, consideremos por un momento que no hay ninguna relación entre los momen-
tos flectores y los desplazamientos, es decir, no se cumple que M (x) = −E J w00 (x). Entonces,

- 90 - Preliminar Versión 2014

El Método de los Elementos Finitos 9. TEOREMAS VARIACIONALES MIXTOS

definamos la siguiente función para los momentos:


x L
M1 2 L
 si 0 ≤ x ≤ 2,
M (x) = (9.3)

M1 2 1 − Lx L


si < x ≤ L,

2

donde M1 es nuestra incógnita momento (incógnita estática), que también está ubicada en el
centro de la luz de la viga en concordancia con la ubicación de la carga concentrada P . Esta
función está representada en la figura 9.3:

Figura 9.3: Momento flector propuesto.

Hemos propuesto nuestras funciones de aproximación o interpolación en función de nues-

tras incógnitas, es decir en función de w1 y M1 , es decir, hemos puesto en evidencia las incógnitas
de nuestro problema y a partir de ellas vamos a hallar nuestras funciones. Prosiguiendo con nues-
tro análisis, plantearemos ahora el Teorema de los Desplazamientos Virtuales,o sea, δw U = δw W .
Definamos cada uno de los términos involucrados.
En primer lugar, analicemos el término correspondiente a la variación de la energía in-
terna. Este término, expresado para un elemento de viga, resulta ser:
Z
δw U = − M (x) δw00 (x) dx. (9.4)
L

Ahora definamos el desplazamiento virtual a partir de una variación de la función des-

plazamiento:

x
δw1 2 L
 si 0 ≤ x ≤ L2 ,
δw(x) = (9.5)
 x

L
δw1 2 1 − L si 2 < x ≤ L,


que está representado gráficamente en la figura 9.4.

Figura 9.4: Desplazamiento virtual propuesto.

Analicemos ahora la ecuación (9.4). Como nuestra función solución para los momentos
flectores es una función interpolante cuyas incógnitas son discreta. Por otro lado, la función des-
plazamiento virtual o la variación de la función desplazamiento elegida es una función lineal,
por lo que la derivada segunda es nula, y nuestra integral no es posible de ser calculada. pero

Versión 2014 Preliminar - 91 -

9.2. Teorema de Hellinger-Reissner El Método de los Elementos Finitos

aprovechando justamente que las incógnitas momento son discretas, podemos discretizar el inte-
gral, o sea, convertirla en una sumatoria. Sabemos, además, que φ = − dw(x)
dx . Con estos cambios,
la (9.4) podemos expresarla así:
Z X dδw(x) X
δw U = −Mi δw00 (x) dx = − Mi = Mi δφi , (9.6)
L dx xi
i i

donde las Mi son nuestras incógnitas de momentos y las δφi son las variaciones de los giros
relativos vinculados con las Mi , que vamos a definir en función de los desplazamientos virtuales
δwi . En nuestro ejemplo, la (9.6) la reescribimos como:
!
δw1 δw1 δw1
δw1 U = M1 L
+ L = M1 4 . (9.7)
2 2
L

Ahora, expresemos el trabajo de las fuerzas exteriores. De nuevo, para nuestro ejemplo,
tenemos:
δw W = P δw1 . (9.8)
Finalmente, igualemos (9.7) y(6.8), y obtengamos M1 :
δw1 PL
M1 4 = P δw1 ⇒ M1 = . (9.9)
L 4
En este caso en particular, con sólo haber planteado el Teorema de los Desplazamientos
Virtuales hemos obtenido la solución para nuestra incógnita momento. Pero nuestro ejemplo está
definido por dos incógnitas, por lo que falta todavía hallar w1 .
Recordemos que para proponer una función interpolante de los momentos habíamos su-
puesto por un instante que no se cumplía M (x) = −E J w00 (x). Sabemos que esto no es rigu-
rosamente cierto, por lo que debemos buscar la forma de hacer que esta relación se cumpla.
Como estamos trabajando con momentos discretos, los Mi , vamos a hacer que dicha relación
sólo se cumpla en forma discreta, es decir, que se cumpla sólo punto a punto, o dicho de otro
modo, solamente donde se evidencien las incógnitas momento (nodo a nodo). Vamos a plantear
lo siguiente:
δMi U (M, w) = δMi U (M ). (9.10)
Esta expresión vincula a las funciones desplazamiento y momento flector, de modo que
sólo se verifique la relación M (x) = −E J w00 (x) en el punto considerado, es decir, no es continua
en el elemento. Esta condición se conoce como Teorema de Hellinger - Reissner.
Vamos a analizar esta igualdad. En primer lugar, propongamos una variación de la función
momento. Esta variación será:

x
δM1 2 L
 si 0 ≤ x ≤ L2 ,
δM (x) = (9.11)
 x

L
δM1 2 1 − L si 2 < x ≤ L,


que está representada en la figura 9.5.

Para el término a la izquierda de la igualdad, podemos usar el mismo razonamiento que
utilizamos con el Teorema de los Desplazamientos Virtuales. Entonces nos queda:
X
δMi U (M, w) = δMi φi , (9.12)
i

dado que ahora la variación es de la función momento (M (x)) propuesta y sabemos que podemos
expresar φi en función de wi . Para nuestro ejemplo, la (9.12) se reduce a:
w1
δM1 U = δM1 4 . (9.13)
L

- 92 - Preliminar Versión 2014

El Método de los Elementos Finitos 9. TEOREMAS VARIACIONALES MIXTOS

Figura 9.5: Variación de la función momento flector.

Nos falta la otra parte de la igualdad. Dado que estamos trabajando con la variación
de la función momento, y debemos relacionarlo con los desplazamientos punto a punto, basta
con que supongamos que en cada punto se cumpla que w00 (xi ) = − ME(xJi ) . Con esta condición,
si suponemos además que E y J son constantes, resulta que el diagrama de momentos flectores
propuesto multiplicado por una constante es también la función w00 (x). En consecuencia, la
segunda parte de la igualdad la podemos escribir así:
Z Z Z
00 M (x) 1
δMi U = − δM (x) w (x) dx = δM (x) dx = δM (x) M (x) dx, (9.14)
L L EJ EJ L
que para nuestro ejemplo no es otra cosa que integrar los diagramas representados en las figu-
ras 9.3 y 9.5. El resultado de esta integración es:
M1
δM1 U = L δM1 . (9.15)
3E J
Ahora igualemos (9.13) y (9.15), y así obtenemos w1 en función de M1 :

w1 M1 M1 L2
δM1 4 = L δM1 ⇒ w1 = . (9.16)
L 3E J 12 E J
Como por la (9.9) tenemos M1 , la reemplazamos en (9.16) y obtenemos w1 :

P L3
w1 = . (9.17)
48 E J
Si comparamos estos resultados con los obtenidos en el ejemplo mencionado y con la
solución de la ecuación diferencial o solución exacta, veremos que mediante la aplicación de
principios variacionales mixtos, los resultados obtenidos corresponden a las soluciones «exactas»
del problema.
Para resumir, podemos generalizar la aplicación de los principios variacionales mixtos con
ayuda del Teorema de Hellinger-Reissner de la siguiente forma:

1. Teorema de los Desplazamientos Virtuales

δu U = δu W (9.18a)

2. Teorema de Hellinger-Reissner

δσ U (σ, ε) = δσ U (σ) (9.18b)

9.2.1. Ecuaciones para vigas sometidas a flexión con teoría de primer orden
Luego de haber resuelto por el Teorema de Hellinger-Reissner una viga simplemente
apoyada, vamos a generalizar las expresiones para un elemento de viga. Partiremos de un elemento
genérico con las siguientes características, como se ve en la figura 9.6:

Versión 2014 Preliminar - 93 -

9.2. Teorema de Hellinger-Reissner El Método de los Elementos Finitos

Figura 9.6: Diagramas para un elemento de viga genérico.

Según hemos visto en los puntos anteriores, la expresión para la energía interna de defor-
mación es: X
δφ U = Mi δφi . (9.19)
i

La variación de los desplazamientos (desplazamiento virtual), y por consiguiente, la va-

riación del giro (giro virtual), los podemos obtener a partir de la función desplazamiento original.
En la figura 9.6 vemos el caso de δwi y los giros δφi−1 , δφi y δφi+1 . En consecuencia, obtenemos
que
δwi
δφi−1 = , (9.20a)
Li−1
δwi
δφi+1 = , (9.20b)
Li

y, por lo tanto, que

δwi δwi
δφi = + . (9.20c)
Li−1 Li

Si ahora reemplazamos estos valores en (9.19) y desarrollemos la sumatoria, nos queda:

 
δwi δwi δwi δwi
δw U = −Mi−1 + Mi + − Mi+1 . (9.21)
Li−1 Li−1 Li Li

Si reagrupamos los términos, nos queda la expresión así:

δwi δwi
δwi U = (Mi − Mi−1 ) + (Mi − Mi+1 ) . (9.22)
Li−1 Li

En esta formulación se ha en cuenta la siguiente convención de signos:

El sentido positivo para el giro y para el momento, es aquel en que las fibras inferiores de
la viga son traccionadas,

La energía es positiva cuando el giros y el momento coinciden en sus sentidos,

- 94 - Preliminar Versión 2014

El Método de los Elementos Finitos 9. TEOREMAS VARIACIONALES MIXTOS

Así, en nuestro ejemplo de análisis, δφi−1 < 0, δφi+1 < 0 y δφi > 0, en tanto que todos
los momentos son positivos.
La variación del trabajo de las fuerzas exteriores se puede expresar como:
Z Z
δwi W = q(x) δwi (x)dx + q(x) δwi (x)dx + P δwi . (9.23)
Li−1 Li

Si igualamos (9.22) y (9.23), obtenemos:

Z Z
δwi δwi
(Mi − Mi−1 ) +(Mi − Mi+1 ) = q(x) δwi (x)dx+ q(x) δwi (x)dx+P δwi . (9.24)
Li−1 Li Li−1 Li

Si Li−1 = Li , entonces la (9.24) queda:

Z
δwi
(−Mi−1 + 2 Mi − Mi+1 ) =2 q(x) δwi (x)dx + P δwi . (9.25)
Li Li

El planteo del teorema de Hellinger-Reissner lo hacemos en forma análoga. El primer

término es similar a la variación de la energía interna de deformación vista:
X
δMi U (M, w) = δMi φi . (9.26)
i

De la figura 9.6 podemos hallar los giros en función de wi , y, en este caso sólo necesitamos
conocer φi :
wi − wi−1 wi − wi+1
φi = + . (9.27)
Li−1 Li
Si reemplazamos la (9.27) en la (9.26) nos queda:
 
wi − wi−1 wi − wi+1
δMi U = δMi + . (9.28)
Li−1 Li
Nos falta la otra parte de la igualdad del Teorema de Hellinger-Reissner. Desarrollemos
el segundo término así:
Z Z
00 M (x)
δMi U (M ) = − w (x) δMi (x) dx = δMi (x) dx
L L EJ
Z Z (9.29)
M (x) M (x)
= δMiIzq (x) dx + δMiDer (x) dx,
Li−1 E J Li E J

donde M (x) es la función momento «real», y δMiIzq (x) y δMiDer (x) es la variación del momento,
como vemos en la figura 9.6.
Si integramos ambas funciones, obtenemos:
 
Li−1 Li
δMi U (M ) = (Mi−1 + 2 M − i) + (Mi+1 + 2 Mi ) δMi . (9.30)
EJ EJ
Ahora igualemos la (9.28) y la (9.30):
   
wi − wi−1 wi − wi+1 Li−1 Li
δMi + = (Mi−1 + 2 Mi ) + (Mi+1 + 2 Mi ) δMi . (9.31)
Li−1 Li EJ EJ
Nuevamente, si Li−1 = Li , la (9.31) queda:
−wi−1 + 2 wi − wi+1 Li
δMi = (Mi−1 + 4 Mi + Mi+1 ) δMi . (9.32)
Li 6E J
Como resumen, si tenemos un sistema estático con m + n incógnitas, donde m son las
incógnitas cinemáticas (desplazamientos) y n las incógnitas estáticas (en este caso, momentos),
podemos expresar la aplicación de los Principios Variacionales Mixtos con ayuda d el Teo-
rema de Hellinger-Reissner de esta forma:

Versión 2014 Preliminar - 95 -

9.2. Teorema de Hellinger-Reissner El Método de los Elementos Finitos

1. Teorema de los Desplazamientos Virtuales

  Z
Mi − Mi−1 Mi − Mi+1
Z
δwi + = q(x) δwi (x)dx + q(x) δwi (x)dx+
Li−1 Li Li−1 Li (9.33)
+ P δwi ,

con i ∈ (1, m), y;

2. Teorema de Hellinger-Reissner
  
wj − wj−1 wj − wj+1 Lj−1
δMj + = (Mj−1 + 2 Mj ) +
Lj−1 Lj EJ
 (9.34)
Lj
+ (Mj+1 + 2 Mj ) δMj ,
EJ

con j ∈ (1, n).

Como caso particular, cuando Li−1 = Li y Lj−1 = Lj , nos queda:

1. Teorema de los Desplazamientos Virtuales

−Mi−1 + 2 Mi − Mi−1
Z
δwi =2 q(x) δwi (x)dx + P δwi , (9.35)
Li−1 Li

con i ∈ (1, m), y;

2. Teorema de Hellinger-Reissner

−wj−1 + 2 wj − wj+1 Lj
δMj = (Mj−1 + 4 Mj + Mj+1 ) δMj , (9.36)
Lj EJ

con j ∈ (1, n).

9.2.2. Ecuaciones de la teoría de segundo orden para vigas

Para este caso nos vamos a limitar a estudiar que ocurre cuando se nos presenta el caso
de inestabilidad del equilibrio. Tomemos una viga empotrada en el extremo inferior y libre en el
superior con una carga concentrada axil de compresión.

Figura 9.7: Carga concentrada de compresión.

- 96 - Preliminar Versión 2014

El Método de los Elementos Finitos 9. TEOREMAS VARIACIONALES MIXTOS

Desplacemos a la viga de su posición de equilibrio, siempre dentro del campo de las

pequeñas deformaciones, como se en la figura 9.7. A continuación, afectemos al sistema despla-
zado con un desplazamiento virtual (o una variación del desplazamiento original) en las mismas
condiciones. Como hemos asumido que estamos en la pequeñas de deformaciones, tendremos:
w ∼ w
tg α ∼
= sen α = =α ⇒ α∼ = , (9.37)
L L
junto con

cos α ∼
= 1. (9.38)

Con estas condiciones, y analizando desde un punto geométrico podemos definir que:

δu = α δw. (9.39)

Si reemplazamos la (9.37) en la (9.39), tenemos:

w
δu = δw. (9.40)
L
La variación del trabajo de las fuerzas exteriores será, entonces:

δu W = P cos α δu ∼
= P δu, (9.41)

por (9.38); y tomando (9.40) podemos expresarlo como

w
δu W = P δw. (9.42)
L

Figura 9.8: Nodo genérico.

Si generalizamos esto para un nodo cualquiera (figura 9.8 nos queda:

 
wi − wi−1 wi − wi+1
δu W = P + δwi . (9.43)
Li−1 Li
Con esta última igualdad obtenemos una expresión más general para el caso de teoría de
segundo orden:
1. Teorema de los Desplazamientos Virtuales
   
Mi − Mi−1 Mi − Mi+1 wi − wi−1 wi − wi+1
δwi + =P + δwi , (9.44)
Li−1 Li Li−1 Li
con i ∈ (1, m), y;
2. Teorema de Hellinger-Reissner
  
wj − wj−1 wj − wj+1 Lj−1
δMj + = (Mj−1 + 2 Mj ) +
Lj−1 Lj EJ
 (9.45)
Lj
+ (Mj+1 + 2 Mj ) δMj ,
EJ
con j ∈ (1, n).

Versión 2014 Preliminar - 97 -

9.3. Teorema de Hu-Washizu El Método de los Elementos Finitos

Al igual que en el caso anterior, si Li = Li−1 y Lj = Lj−1 , las ecuaciones (9.44) y (9.45)
quedan:

1. Teorema de los Desplazamientos Virtuales

−Mi−1 + 2 Mi − Mi−1 −wi−1 + 2 wi − wi+1
δwi =P δwi , (9.46)
Li−1 Li

con i ∈ (1, m), y;

2. Teorema de Hellinger-Reissner
−wj−1 + 2 wj − wj+1 Lj
δMj = (Mj−1 + 4 Mj + Mj+1 ) δMj , (9.47)
Lj EJ

con j ∈ (1, n).

9.3. Teorema de Hu-Washizu

Existe un segundo teorema más general todavía que el de Hellinger-Reissner. Es el Teo-
rema de Hu-Washizu. Este teorema, además de considerar las condiciones de Hellinger-Reissner,
establece que tampoco se debe cumplir en el continuo la relación cinemática (ε 6= du dx en el
continuo).
Análogamente al teorema de Hellinger-Reissner, vamos a proponer en este caso tres fun-
ciones interpolantes: para los desplazamientos, para las tensiones, y para las deformaciones. Es
decir, vamos a desacoplar todas las funciones que intervienen en la resolución de los problemas y
luego las vamos a ir relacionando punto a punto. Esto nos permite definir tres funciones distintas,
pero de similares características para cada una de las incógnitas definidas.
Las expresiones generales que definen el Teorema de Hu-Washizu son:

1. δu U = δu W ;

2. δσ U (σ, ε) = δσ U (σ), y;

3. δε U (ε, u) = δε U (ε).

- 98 - Preliminar Versión 2014

 APÉNDICE A

EL USO DE LA COMPUTADORA

A.1. Introducción
Hasta este punto, lo que se hemos visto son los diferentes métodos que se pueden usar para
resolver un determinado problema estructural, como son los teoremas energéticos, el método de
Rayleigh-Ritz, los elementos finitos propiamente dichos, como así también la definición de varios
tipos de elementos según sus características estructurales.
Cada uno de estos elementos vincula una cierta cantidad de nodos, y el tipo de elemento
define los grados de libertad de cada uno de éstos (los desplazamientoas asociados), que no son
otra cosa que las incógnitas de nuestro problema.
Esto se traduce en que una vez armada la malla de discretización de nuestro modelo
estructural, tendremos un sistema de ecuaciones lineales con un número grande o muy grande de
incógnitas. Es decir, tendremos que resolver un sistema de cientos o miles de ecuaciones, según
la densidad de malla. Evidentemente, resolver esta cantidad de ecuaciones no resulta práctico
mediante métodos manuales o algebraícos tradicionales. Surge inmediatamente la necesidad de
contar con una herramienta de cálculo un poco más poderosa que una simple calculadora, que
es la computadora.
El desarrollo del método, que comenzó en los años sesenta (aunque el método en sí mismo
es de los años 50) y sigue hasta hoy, introdujo una serie de métodos y/o programas de computación
para resolver sistemas de ecuaciones, teniendo como base la optimización en el uso de la memoria
disponible y de la velocidad de cálculo o procesamiento; en muchos casos fueron especialmente
diseñados para la aplicación de los elementos finitos. Muchos de los métodos numéricos existentes
resultaron no tan eficientes cuando el número de ecuaciones a resolver se volvía muy grande, o
cuando la necesidad de mejorar o refinar los resultados requería una modificación de la malla,
así sea cambiando su distribución o aumentando la densidad de la misma.

A.2. Los métodos numéricos

Para resolver un sistema de ecuaciones mediante el uso de la computadora se dispone de
varios métodos, aplicables según el tipo de sistema que se deba resolver. En cursos de Análisis
numérico se estudian varios de esos métodos como parte de un curso normal. Algunos de los
métodos más conocidos son: el Método de Eliminación de Gauss, la factorización de matrices,
entre ellos el método de Cholesky, que componen los denominados métodos directos. Otro de los
métodos más conocidos son el método de Gauss-Seidel y el Método de los Gradientes Conjugados,
que integran el conjunto de los métodos iterativos.

99
A.3. Sistemas de ecuaciones lineales El Método de los Elementos Finitos

Cada uno de ellos puede aplicarse a casi todos los sistemas de ecuaciones, en especial
el método de eliminación de Gauss, el método directo más usado. En los demás, si se tiene
un sistema expresado de la forma A x = B, la matriz A debe cumplir con ciertas condiciones
adicionales. Si se pudiera conocer de antemano que esta matriz de coeficiente siempre cumple
con las condiciones necesarias para aplicar determinado sistema, podría utilizarse un método
determinado. Si bien esto no siempre es posible, resulta que en muchos casos esta matriz sí
cumple con ciertas condiciones cualquiera sea el sistema a resolver. En el caso particular de los
sistemas estructurales esta matriz A suele conocerse como la matriz de rigidez K, la cual siempre
es definida positiva y simétrica. Esto es una gran ayuda, ya veremos por qué.
Hoy en día casi no hay limitaciones para la aplicación del método de los elementos fi-
nitos. El avance en su desarrollo ha llevado, por un lado, a extenderlo a muchas disciplinas de
la ingeniería, que anteriormente hacían uso del método de las diferencias finitas para resolver
ecuaciones diferenciales, y por otro, a producir elementos con más cantidad de nodos (incógnitas)
con el fin de mejorar la capacidad de modelación y alcanzar en algunos casos mayor estabilidad
y convergencia en los resultados. Por lo tanto, según el tipo de problema que se tenga, la matriz
de coeficiente será distinta. Pero de nuevo hay una ventaja: la matrices suelen ser simétricas.
También debemos tener en cuenta que esos coeficientes que integran la matriz generalmen-
te salen de cálculos que deben efectuarse en el programa, usualmente por integración numérica.
Esto significa que no sólo se deben disponer de rutinas que resuelvan ecuaciones lineales, si no
también que calculen integrales en forma numérica. Y no debe olvidarse que el proceso de de-
finición del elemento lleva consigo la utilización de funciones de interpolación numérica, con lo
cual previamente debe definirse el tipo de función a emplear.
Con todo lo dicho, puede verse que el método de los elementos finitos es una notable forma
de integrar el conocimiento físico-matemático (mecánica, análisis estructural, hidráulica, física
en general y matemática) con procedimientos numéricos de resolución. Sin entender el problema
que debe ser resuelto en forma conceptual, difícilmente se podrá obtener una solución numérica
válida o aplicable.

A.3. Sistemas de ecuaciones lineales

Durante algún tiempo los métodos directos para resolver sistemas de ecuaciones lineales,
esto es, eliminación de Gauss o factorización de matrices, fueron los métodos más utilizados, sobre
todo para matrices completas o casi completas. El desarrollo mostró que esta característica de las
matrices no es el caso más común en el método de los elementos finitos, si bien es muy aplicado
el método de Cholesky, una variante de la factorización de matrices cuando dichas matrices son
simétricas y definidas positivas.
Pero, el armado de un modelo de elementos finitos, esto es, el diseño de la malla, suele
generar matrices denominadas ralas, es decir, con muchos coeficientes nulos. En muchos ocasiones,
esto se complementa con la generación de matrices banda o seudo–banda, donde el ancho de la
misma va cambiando a lo largo de la diagonal principal. Con lo cual un método directo puede
no resultar muy eficiente a la hora de resolver el sistema.
Esta situación hace atractiva la aplicación de métodos iterativos, como el de Gauss-Seidel
o el de sobrerrelajaciones, que si bien en principio son aproximados, disminuyen notablemente
las cantidad de operaciones que deben realizarse en el caso de matrices ralas. Pero para poder
aplicar estos métodos resulta imprescindible que la matriz de coeficientes A cumpla P con una
nueva condición: que sea estrictamente diagonal dominante es decir, que |aii | > j |aij | con
j 6= i. Un vez más, resulta que una buena parte de las matrices de coeficientes cumplen con esta
condición.
Consecuentemente, en la resolución de sistemas de ecuaciones lineales encontraremos que
los métodos más utilizados son el método de eliminación de Gauss (o variantes), el método de
Cholesky o los métodos iterativos. En los años setente se desarrollaron algunas variantes, como

- 100 - Preliminar Versión 2014

El Método de los Elementos Finitos A. EL USO DE LA COMPUTADORA

el método de eliminación frontal, una mejora para operar con eliminación de Gauss y el método
del gradiente conjugado, un método iterativo muy poderoso.

A.3.1. Métodos tradicionales

Hemos comentado que los método más usados para resolver los sistemas son el de elimi-
nación de Gauss, el método de Cholesky y el método de Gauss-Seidel.
Los dos primeros métodos operan directamente sobre la matriz de coeficientes y la trans-
forman para obtener un sistema triangular en el primer caso, y dos matrices triangulares, en
el segundo. Eliminación de Gauss transforma en realidad una matriz ampliada, pues no sólo
modifica la matriz de coeficientes (A o K), sino que también modifica el vector independiente
(B o R), crenado un nuevos sitema triangular, que luego se resolverá medainte el procedimiento
de sustitución inversa. En cambio, el método de Cholesky crea una matriz triangular S, tal que
ST S = A. Éste tiene la gran ventaja de calcular solamente una matriz triangular, pero la matriz
A debe ser simétrica definida positiva, pues los coeficientes sii son raíces cuadradas.
El método de Gauss-Seidel es iterativo, y necesita un vector inicial para poder hallar
la solución del sistema. Mediante iteraciones va corrigiendo el valor anterior (en el primer paso
corrige el vector inicial) y si converge, se obtendrá una solución muy precisa, que estará fijada
por la tolerancia establecida. La matriz A debe ser simétrica definida positiva para que converja.
Sin embargo, estos procedmientos o métodos mostraron no ser muy eficientes cuando las
matrices son de gran tamaño, tanto ralas como completas. Además, en muchos casos, las matrices
ralas son banda, o cuasi banda, lo que complica aun más el lograr un buen aprovechamiento de
la memoria disponible.
Existe otro problema. Cuando se trabaja con grandes matrices ralas, almacenar o reservar
memoria para coeficientes nulos no es la mejor forma de aprovecharla. Significa «gastar» recursos
que pueden utilizarse para agilizar el cálculo. Un segundo problema es que el cero (0) no es un
número fácil de manejar en las computadoras. Cualquier error en un paso que deje un coeficiente
nulo en las matrices o que deba guardarse un valor que se entienda como nulo pero que en
realidad no lo es, puede llevar auna solución distorionada o no válida.
La necesidad de reducir el uso de la memoria y optimizar los recursos disponible llevó
a que a fines de la década del sesenta y principios de los setenta, aparecieran algunos métodos
diseñados especialmente para mejorar los algoritmos de cálculo de los programas de análisis por
elementos finitos.

A.3.2. Métodos más modernos

Método de Sustición Frontal

Una de las técnicas surgidas para optimizar el uso de la computadora se publicó en el

año 1970 y fue desarrollada por Bruce Irons. Es la conocida como Método de Sustitución (o
Eliminación) Frontal, y se basa en optimizar el método de eliminación de Gauss para matrices
simétricas y definidas positivas.
Es una técnica específicamente desarrollada para la resolución de sistemas de ecuaciones
lineales derivadas de la aplicación el método de los elementos finitos. El problema a resolver es
el tradicional, x = b, donde A es una matriz simétrica definda positiva, la que normalmente
llamamos K.Además, esta matriz generalmente suele ser una matriz banda, con lo cual buena
parte de los coefcientes son nulos, esto es, que aij = 0 si |i − j| > n, siendo n el ancho de banda.
Para ver brevemente cómo funciona el procedimiento, analicemos el método de elimina-
ción de Gauss. Para transformar la matriz tenemos que:
ais ais
a∗ij = aij − asj y b∗i = bi − bs ,
ass ass

Versión 2014 Preliminar - 101 -

A.3. Sistemas de ecuaciones lineales El Método de los Elementos Finitos

donde los a∗ij son los coeficentes de la matriz transformada (A∗ ) y los b∗i son los coeficientes del
vector independiente transformado (B∗ ).
Por lo tanto, si cualquier ais o asj es nulo, implica que a∗ij = aij y que b∗i = bi . Esto puede
interpretarse como que los aij deben quedar en la matriz cuando ambos, ais y asj sean distintos
de cero, pues sólo en ese caso se modificarán; lo mismo pasa con los bi .
La técnica de sustición frontal consiste en ensamblar y eliminar coeficientes. Ensambla
la matriz global a medida que se van activando las incógnitas, y las va eliminando cuando las
incógnitas aparecen por última vez. En realidad no elimina los coefcientes sino que guarda los
coeficientes en disco, pues sólo las necesita cuando efectúe el proceso de sustitución inversa. De
este modo, se consigue una gran ventaja, porque deja de ser importante el ancho de banda de
la matriz A, o dicho de otro modo, no es necesario ordenar los nodos de manera de reducir el
ancho de banda. Lo que va a importar es el ordenamiento (o numeración) de los elementos.
Una técnica frontal es generalmente más eficiente que las técnicas tradicionales. Sobre
todo cuando el desarrollo de nuevos elementos se concentra en aumentar los grados de libertad
disponible por elemento, fundamentalmente en elementos volumétricos o de tres dimensiones.
Esta técnica es más eficiente que el mejor método de eliminación de Guass.
Como ejemplo se puede ver el siguiente caso. Si se tiene la malla de la figura A.1 y se
considera que la misma no es lo suficientemente densa, modificarla puede ser un problema. En

Figura A.1: Malla gruesa.

efecto, hacer la malla más densa en los alrededores del nodo 7 requiere volver a numerar todo el
modelo para lograr un ancho de banda reducido, si se quiere aplicar el método de eliminación de
Gauss como método de resolución de sistemas de ecuaciones. Esta nueva malla queda definida
como se ve en la figura A.2.

Figura A.2: Malla más densa.

Al usar la técnica frontal, en cambio, no es necesario modificar la numeración original, y

alcanza con agregar los nuevos nodos, como se ve en la figura A.3, pero sí es necesario volver a
numerar los elementos.

- 102 - Preliminar Versión 2014

El Método de los Elementos Finitos A. EL USO DE LA COMPUTADORA

Figura A.3: Malla más densa sin modificar la numeración original.

Método de los Gradientes Conjugados

El otro de los métodos de uso generalizado es el Método de los Gradiente Conjugados.
Se trata de un método iterativo, como el Gauss-Seidel pero cuya convergencia es mucho más
segura. Los métodos iterativos parten de resolver el problema que ya hemos planteado, es decir,
A x = B ⇒ B − A x = 0.
Podemos escribir lo anterior como:
Px = Px − Ax + B ⇒
P x = (P − A) x + B ⇒
−1 (A.1)
x=P (P − A) x + P−1 B o,
x = P−1 P − P−1 A x + P−1 B = I − P−1 A x + P−1 B.

Si definimos que T = I − P−1 A y C = P−1 B, entonces nos queda:

x = T x + C, (A.2)
que al convertirla en un método iterativo queda así:
xhi+1i = T xhii + C. (A.3)
De esta expresión surgen los distintos métodos iterativos. Si A se descompone en tres
matrices, A = L + D + U (o A = D(L + I + U)), donde L es una matriz triangular inferior, D es
una matriz diagonal y U es una matriz triangular superior (e I es la matriz identidad), el método
de Gauss-Seidel corresponde al caso en que la matriz P es igual a L + D (P = L + D). Otros
métodos son el de Jacobi y el de Sobrerrelajaciones sucesivas. Todos estos métodos se denominan
iterativos estacionarios.
El método del gradiente conjugado surge de plantear que P = α1 I y que:

xhi+1i = xhii + αi dhii . (A.4)

El algoritmo del Método de los Gradientes Conjugados más usual es el siguiente:
rh0i = dh0i = B − A xh0i
T
rhii rhii
αi = T
dhii A dhii
xhi+1i = xhii + αi dhii
(A.5)
rhi+1i = rhii − αi A dhii
T
rhi+1i rhi+1i
βi+1 = T
rhii rhii
dhi+1i = rhi+1i + βi+1 dhii .

Versión 2014 Preliminar - 103 -

A.4. Integración numérica El Método de los Elementos Finitos

Este método está incluido en los denominados iterativos no estacionarios, también cono-
cidos como métodos variacionales.

A.4. Integración numérica

Otro punto que debe resolverse en forma numérica es la integración para obtener las
matrices de rigidez. Esto también lleva a la utilización de métodos específicos.
En la aplicación del método de los elementos finitos se conocen las funciones de interpo-
lación, funciones que deberánn ser integradas. En efecto, como se ha visto en los capítulos an-
teriores, las funciones de interpolación utilizadas para representar las funciones desplazamientos
de los distintos elementos estudiados son polinómicas. Y para este tipo de funciones un método
muy apropiado y eficiente para integrar numéricamente es la cuadratura de Gauss–Legendre,
tanto en una dimensión como en dos.
En el caso de una dimensión, el procedimiento es sencillo. La expresión general para este
caso es: Z b n
f (x) dx ∼
X
= ci f (xi ), (A.6)
a i=1
pero mediante una adecauda elección de los puntos de evaluación en un intervalo [−1; 1] se llega
a la expresión de la cuadratura de Gauss–Legendre:
Z 1 n
f (t) dt ∼
X
= ci f (ti ), (A.7)
−1 i=1
Con esta expresión y para n = 1; 2; 3; . . ., se obtienen las raíces de un polinomio de grado
2n − 1 (t1 ; t2 ; t3 ; . . . ; tn ) y los coeficientes ci (c1 ; c2 ; c3 ; . . . ; cn ). Es decir, para un determinado se
tiene un expresión diferente, por ejemplo:
1. Para n = 2 tenemos Z 1
f (t) dt = c1 f (t1 ) + c2 f (t2 ) (A.8a)
−1

2. Para n = 3 tenemos
Z 1
f (t) dt = c1 f (t1 ) + c2 f (t2 ) + c3 f (t3 ) (A.8b)
−1

3. Para n = 4 tenemos
Z 1
f (t) dt = c1 f (t1 ) + c2 f (t2 ) + c3 f (t3 ) + c4 f (t4 ). (A.8c)
−1

El valor de n depende en gran medida de la función a integrar. En el caso de funciones

polinómicas, el resultado es exacto para polinomios de grado 2n − 1 o menor. Para cambiar el
intervalo de integración, basta con plantear un cambio de coordenadas, como el siguiente:
a = m (−1) + n, (A.9a)
b = m t + n, (A.9b)
cuya solución es m = b−a a+b b−a a+b
2 y n = 2 y entonces, x = 2 t + 2 . Por lo tanto, si reescribimos
todo en función de la variable t, para el caso n = 2, nos queda los siguiente:
Z b Z 1
b−a b−a
f (x) dx = f [x(t)] dt = (c1 f [x(t1 )] + c2 f [x(t2 )]) . (A.10)
a −1 2 2
Así como hemos aplicado esto a funciones en una dimensión, también se puede aplicar la
cuadratura de Gauss para integrales en dos y tres dimensiones 1 .
1
Para más detalles acerca de la integración numérico en dos más dimensiones con cuadratura de Gauss–
Legendre, consultar cualquier libro de Análisis numérico.

- 104 - Preliminar Versión 2014

El Método de los Elementos Finitos A. EL USO DE LA COMPUTADORA

A.5. Los programas de análisis por elementos finitos

A.5.1. Conceptos generales
Hemos visto hasta ahora los distintos elementos para poder armar un modelo de análisis
por elementos finitos, cómo se calculan las matrices de rigidez, las funciones interpolantes que se
utilizan, la forma de ensamblar las matrices locales y cómo resolver un sistema con un análisis
por teoremas variacionales mixtos. Todo dentro del análisis estructural.
Pero el método de los elementos finitos no sólo sirve para resolver problemas estructu-
rales. Se lo aplica también en otros campos de la ingeniería, como en problemas de difusión
(transmisión del calor), en electricidad, mecánica de los fluidos, geotecnia, hidráulica marítima,
etc.. Es evidente que dados el amplio rango de su utilización, existen cada vez más programas
específicos para cada caso a estudiar. Hay programas para análisis estructural exclusivamente,
otros que integran problemas estructurales y mecánicos, programas que analizan fenómenos ma-
rítimos (como la agitación en recintos casi cerrados), o que sirven para el análisis de redes de
filtración en problemas relacionados con presas.
Con todo, los programas no suelen ser muy diferentes entre sí. La organización es similar
y según el tipo de caso a resolver, pueden adoptar diferentes métodos numéricos de resolución,
particularmente en el caso de resolución de sistemas de ecuaciones lineales. Un caso aparte son
aquellos que resuelven sistemas no lineales.
Un punto importante es conocer qué tipos de datos son necesarios ingresar a los progra-
mas. No se requieren los mismos datos para un modelo de análisis estructural que para el cálculo
de redes de filtración, eso es obvio. Analicemos en detalle con ayuda de un modelo estructural
cuáles son los datos necesarios.
Hemos visto que el método de los elementos finitos es la transformación de un sistema de
ecuaciones diferenciales en un sistema de ecuaciones lineales de la forma A x = B (o vcK U = R).
También hemos visto una breve descripción de los métodos que pueden utilizarse para resolver
este tipo de sistemas, dependiendo de las características de la matriz A (o K). La matriz A (o
K, en este caso), va depender del tipo de modelo a analizar (unidimensional, espacial) y del tipo
de elemento (barra, viga, etc.). Sin embargo, hay algo que sí se sabe de antemano de la matriz
K: es simétrica y definida positiva. Aún sabiendo esto subsisten algunos problemas.
Los datos que deben ingresarse en los programas son generalmente los mismos, a sa-
ber: las coordenadas de los nodos (x, y, z), las características de los materiales, (E, G, ν, λ), las
características geométricas de los elementos (A, Iyy , Izz , Ip ) y el tipo estructural (viga, barra,
columna, placa), las condiciones de borde (apoyos), las cargas y/o las combinaciones de cargas,
y en algunos casos, el tipo de análisis (primer o segundo orden en el caso de estructuras).
Cada programa tiene, en consecuencia, su propia forma de ingresar los datos. Los más
sencillos suelen requerir, por ejemplo, que uno arme la malla y defina los nodos, con lo cual se
deben ingresar las coordenadas de todos los nodos. Esto es factible de hacer en modelos sencillos
(de pocos nodos) pero muy trabajoso en grandes modelos.
Del mismo modo, los programas más elementales también dejan en manos del usuario la
vinculación de los nodos, es decir, la definición de los elementos que vinculan los nodos mediante
un procedimiento manual.
Ingresar las coordenadas nodos y los elementos está asociado a la identificación o nume-
ración de unos y otros. Y esta numeración incide directamente en la matriz de rigidez. Hemos
dicho que las matrices para el caso del análisis estructural son simétricas, definidas positivas y
ralas, con prepondrrancia de matrices banda. Una mala numeración puede hacer que una matriz
banda tenga un ancho de banda excesivo (como se vió en el caso en la comparación del método
frontal y eliminación de Gauss). Muchos programas tienen en entre sus limitaciones, el ancho de
banda, el número de nodos y/o el número de elementos a ingresar.
Estas limitaciones (aparte de las limitaciones de las versiones para prueba) están relacio-
nadas con la capacidad de manejar las matrices resultantes. Una matriz con un ancho de banda

Versión 2014 Preliminar - 105 -

A.5. Los programas de análisis por elementos finitos El Método de los Elementos Finitos

muy grande es muy lenta en su procesamiento, no así una de ancho de banda chico. Por esta
razón existen programas complementarios cuyo único fin es diseñar la malla de elementos para
optimizar el proceso de cálculo. Algunos programas tienen procedimientos que generan mallas
automáticas muy eficientes, pero aún así se está volviendo muy común el usar programas externos
para definir la malla más eficiente para un determinado problema. Estos programas se conocen
como de preprocesamiento y una vez determinada la malla, guardan los resultados en archivos
que tienen el formato según el programa por el que se efectuará el análisis. Con esto se pueden
analizar distintas mallas para un mismo tipo de problema de forma de obtener los resultados
más adecuados o aproximados a las necesidaes del ingeniero.
Un segundo punto es que los resultados obtenidos luego del cálculo son una colección de
números a veces ininteligibles en forma individual. Las planillas con los resultados completos de
un modelo pueden ocupar cientos de páginas y ser de difícil interpretación a simple vista. Muchos
programas cuentan con unidades de procesamiento de estos resultados a través de interfaces
gráficas, que visualizan las tensiones en los elementos o los desplazamientos, y en el caso particular
de los programas de cálculo estructural, grafican los diagramas de características (M, N , Q).
También para esto existen programas de postprocesamiento, que suelen ser los mismos
que hacen el preprocesamiento. Estos programas leen los archivos de salida de los progrmas de
cálculo y luego generan gráficos, tablas, curvas a partir de los datos originales, que ayudan a una
mejor interpretación de los resultados obtenidos.
Estos programas son muy usados en otras disciplinas o ramas, com es el caso de la
hidráulica marítima, donde se cuenta con programas que solamente realizan el cálculo del sistema
de ecuaciones , en tanto que el postprocesamiento se hace con otro programa que permite graficar
curvas o superficies, que hacen la interpretación de los resultados obtenidos más amigable. Un
ejemplo de este tipo de programas es el FEMAP, un programa de pre y postprocesamiento de
datos provenientes de otros programas de análisis estructural.

A.5.2. Programas de análisis estructural por elementos finitos

Gracias al análisis de estructuras por el método de los elementos finitos, hoy se pue-
den resolver una gran cantidad de sistemas estructurales. Principalmente este avance se ve en
el análisisde estructuras tridimensionales, algo muy díficil de llevar adelante con los métodos
anteriores.
El avance y la investigación en distintos elementos con variadas características, ha hecho
posible disponer de una variada librería de elementos para modelar estructuras muy complejas.
Es común que los programas comerciales de gran envergadura dispongan en sus bibliotecas de
elementos de una gran cantidad de elementos planos, triangulares y cuadriláteros, con tres o seis
nodos, para el caso triangular, o de cuatro, ocho y nueve nodos, para el caso del cuadrilátero,
como así también la incorporación cada vez más de elementos tridimensionales como el tetraedro
y el cubo. En algunos casos también se dispone de rutinas para el cálculo de estructuras según
la teoría de segundo orden (pandeo o "buckling"), facilitando el dimensionado y verificación de
los elementos estructurales.
En el mercado existen varios programas comerciales para el análisis estructural como
así también que pueden aplicarse tanto a ingeniería civil como a ingeniería mecánica. Entre
estos programas se pueden mencionar el STRAP(Structural Analysis Programs), el STAAD Pro
(Structural Analysis And Design Software), el SAP 2000, GT STRUDL (Structural Design and
Analysis Software), orientados principalmente al análisis estructural, y el NISA, o el SolidWorks
que tienen aplicación también en otros campos, además del análisis estructural.

- 106 - Preliminar Versión 2014

El Método de los Elementos Finitos Bibliografía

BIBLIOGRAFÍA

[1] Baker, A.J. y Pepper, D.W. Finite Elements 1-2-3. McGraw-Hill. 1991.

[2] Bathe, K.J. Finite Element Procedures. Prentice–Hall. 1996.

[3] Chandrupatla, T.R. y Belegundu, A.D. Introducción al Elemento Finito en Ingeniería.

Prentic–Hall. Segunda edición. 1999.

[4] Cook, R.D., Malkus, D.S. y Plesha, M.E. Concepts and applications of finite element analy-
sis. John Wiley & Sons. 1989.

[5] Ferguson, J. A Brief Survey of the History of the Calculus of Variation and its Applications.
University of Victory. Canada. 2004. (URL: http://arxiv.org/abs/math/0402357)

[6] Fischer, J. Introduction to the Calculus of Variations.

[7] Fung, Y.C. Foundations of Solid Mechanics. Prentice–Hall. 1965.

[8] Hildebrand, F.B. Métodos de la matemática aplcada. Eudeba. 1973.

[9] Hughes, T.J.R. The Finite Element Method. Linear Static and Dynamic Finite Element
Analysis. Dover Publications, Inc. 2000.

[10] Hutton, D.H. Fundamentals of Finite Element Analysis. McGraw–Hill. 2004.

[11] Irons, B.M. A Frontal Solution for Finite Element Analysis.International Journal for Nu-
merical Methods in Egineering, Vol. 2, 5–32. 1970.

[12] Krasnov, M.L., Makarenko, G.I. y Kiseliov, A.I. Cálculo variacional (ejemplos y problemas).
Editorial Mir. 1976.

[13] O’Connor, J.J. y Robertson, E.F. The Brachistrocrone problem. The MacTutor History of
Mathetmatics History. University of Saint Andrews. Scotland. 2002. (URL: http://www-
history.mcs.st-andrews.ac.uk/history/HistTopics/Brachistochrone.html)

[14] Prathap, G. Finite Element Analysis as Computation.

URL: http://www.cmmacs.ernet.in/cmmacs/online_contents.html.

[15] Prathap, G. y Mukherjee, S. Analysis of delayed convergence in the three-noded Timoshenko

beam element using the function space approach. Sadhana, Vol 27, Part 5, October 2002.

[16] Puppo, A. H. Los Teoremas de Energía en la Mecánica del Sólido. Anales de la Sociedad
Científica. T. CXCIV, Entrega V-VI. 1972.

Versión 2014 Preliminar - 107 -

BIBLIOGRAFÍA El Método de los Elementos Finitos

[17] Samarski, A.A. Introducción a los Métodos Numéricos. Editorial Mir. 1986.

[18] Timoshenko, S. Resistencia de Materiales - Tomos I y II. Espasa-Calpe S.A. 1957.

- 108 - Preliminar Versión 2014