Analisis-del-Creep-en-Placas-por-Simulacion-Numerica

Vista previa del material en texto

Abril 2016 
 
Instituto Politécnico Nacional 
ESCUELA SUPERIOR É INGENIERÍA MECÁNICA Y ELÉCTRICA 
SECCIÓN DE ESTUDIOS DE POSGRADOS E INVESTIGACIÓN 
“Análisis del Creep en Placas por Simulación Numérica” 
 
 
 
T E S I S 
 
QUE PARA OBTENER EL GRADO DE 
 
MAESTRO EN CIENCIAS EN INGENIERÍA MECÁNICA 
 
PRESENTA 
 
Jonathan Saucedo Jardón 
 
ASESORES 
 
Dr. Didier Samayoa Ochoa 
 
Dr. Ernesto Pineda León 
 
 
 
 
 
 
 
Resumen
Actualmente hay en existencia un gran número de programas compu-
tacionales, que nos reproducen con exactitud el comportamiento de las de-
formaciones en fenómenos no linealeas como lo son la elastoplasticidad, vis-
coelasticidad y el creep. Destaca por su uso comercial el método del elemento
finito, técnica en la que se han hecho grandes esfuerzos en el desasarrollo de
nuevos algoritmos que alcanzen con mayor rapidez la solución a los proble-
mas de naturaleza no lineal. Aunque el proceso iterativo para la solución de
un problema no lineal se alcanze rapidamente. El usuario tendrá que correr
una nueva simulación por cada valor de tensión o condiciones de frontera que
desee o aquellas que el problema requiera, consumiendo mucho tiempo en
esta tarea que pudiera llegar a ser frustrante conciderando todas las posibi-
lidades de carga.
Es por eso que en esta tesis se intentó hallar una relación lineal entre la
tensión y el esfuerzo aplicado en un punto de la superficie de la placa con el
propósito de predecir la velocidad de deformación del creep usando los cri-
terios de time− hardening y strain− hardening. Una placa barrenada fue
probada conciderando tensión uniaxial constante durante 1000 horas usando
el programa ADINA para llevar a cabo el proceso no lineal del elemento
finito.
De los datos obtenidos en tres pruebas de tensión, se encontró una relación
lineal entre la tensión y el esfuerzo aplicado en cualquier punto x, y del do-
minio de la placa.
Para verificar nuestro argumento. Se recopiló información concerniente al
esfuerzo aplicado sobre la placa debido a la tensión agrupando nodos de
nuestro dominio discretizado por medio de rutas. De cada nodo de las rutas
contempladas se registró el de esfuerzo aplicado por cada valor de tensión,
obteniendo de esto una relación lineal. Esta nos permite calcular el esfuerzo
aplicado σyy a cualquier tensión y como consecuencia directa podemos cono-
cer las deformaciones por creep sin necesidad de llevar a cabo una prueba
numérica o experimental. La curva del creep ontenida por simulación numéri-
ca aplicando elemento finito fue comparada con nuestros resultados.
La expresión σyy ∼ mσ0 puede substituirse directamente en el criterio del
creep conocido como endurecimiento por deformación con el objetivo de co-
nocer la deformación y la velocidad de deformación por tensión uniaxial cons-
tante . Por medio de la comparación de curvas se observó que los resultados
obtenidos son acordes con los numéricos.
i
Abstract
There are lots software packages nowadays that closely depict nonlinear
strain behaviour such as elastoplastic, viscoelastic and creep phenomena.
Using numerical techniques like finite difference, finite element and boun-
dary element method among the most popular.
The finite element method in which great improvements are being made
with the goal to develope algorithms that yield faster a solution for nonli-
near inelastic problems leads in commercial use , even though the numerical
process seems to be the best choice for assess creep phenomenon in a wide
range of structural and mechanical components, FEM user has to perform
a new simulation each time the stress , boundary conditions and material
selection changes due to component design requirements or something else
making numerical simulation process a long task.
Therefore this thesis intended found by a data base the kept pattern between
remote stress σ0 and normal stress σyy in various points of the plate domain,
thereby stress term in Norton−Bailey creep equation could be replaced by
the σ0 vs σyy ratio.
Hence creep strain and creep strain rates in yy direction can be computed
applying time and strain hardening criterions at any point included in data
base without FEM or any other numerical tecnique.
The geometry used for reach the goals mentioned in last paragraph was a
plate with a hole , in which uniaxial stress magnitudes of 150 , 175 and 200
was applied and a total time of 1000 hours was input in ADINA for bring
out finite element procedure. Curve comparison of creep strain and creep
strain rate results was done in order to restate our assumption , from this
comparison has been observed that our creep behaviour is in agreement with
numerical.
ii
Introducción
En el presente trabajo se aborda el tema de las deformaciones por creep
en su etapa primaria a causa de esfuerzos de tensión aplicados en una placa
barrenada al centro y bajo la concideración de esfuerzo plano. El modelo
constitutivo que se utilizó para obtener las deformaciones fue el de endureci-
miento por tiempo(time hardening creep strain), mientras que para conocer
la velocidad de deformación recurrimos a la ley de endurecimiento por defor-
mación (strain hardening creep strain rate). También se hace mención a
las modificaciones que se realizan a ambos criterios para predecir las curvas
de deformación contra tiempo a esfuerzo variable.
Dentro del primer caṕıtulo encontraremos información concerniente al com-
portamiento mecánico de metales y aleaciones cargados a tensión y expuestos
a elevadas temperaturas, condiciones que favorecen la aparición del fenómeno
del creep y su representación deacuerdo a los modelos constitutivos mencio-
nados con anterioridad.
En el caṕıtulo 2 exponemos los constituyentes de los métodos numéricos ta-
les como diferencias finitas, elementos de frontera y elementos finitos.
En los que en primer lugar se pretende solucionar la ecuación diferencial del
Laplace la cual gobierna lo relativo a elasticidad, esfuerzos y deformaciones
por medio de condiciones de frontera bien definidas y haciendo uso de una
discretización. Habiendo asentado estas bases procedemos a la deducción de
las funciones forma de un elemento cuadrático de 8 nodos deacuerdo al plan-
teamiento isoparámetrico. De igual manera se obtuvo la funcional de enerǵıa
de un sistema elástico por medio del Principio de Trabajo V irtual (PTV ),
hasta llegar a una expresión lineal con la que se obtiene la matriz de rigidez
de un cuerpo.
En el caṕıtulo 3 se presenta una breve explicación del comportamiento no li-
neal y la ley de Hook generalizada para representar tales fenómenos, como el
creep es un tipo de deformación que depende del tiempo, se hace un análisis
del movimiento para conocer los desplazamientos en diferentes instantes. Se
menciona también el sistema no lineal de ecuaciones que gobierna al creep y
su solución por medio del método del elemento finito el
cual nos lleva a desarrollar un proceso iterativo que puede ser solucionado
por el método de Newton−Raphson y Newton−Raphson modificado.
El contenido del caṕıtulo 4 lleva a cabo la simulación del problema del creep
en su primera etapa tomando en cuenta tres diferentes magnitudes de esfuerzo
uniaxial medida en [MPa]. Se analizan los datos de las pruebas realizadas en
iii
diferentes puntos del dominio percatandonos que existe una relación lineal
entre el esfuerzo de ensayo y el esfuerzo en cualquier punto del dominio.
Entonces usando la Ley Potencial del Creep la cual contempla los esfuerzos
para el calculo de las deformaciones y por consecuencia las velocidades de
deformación según los criterios de time− hardening y strain− hardening.
Finalmente comparamos las curvas hechas por medio de la relación lineal
encontrada y planteamiento el M.E.F. Para culminar con este trabajo en el
caṕıtulo 5 se reflexiona sobre el alcance de los resultados obtenidos, aśı como
propuestas de programación para trabajos futuros.
iv
Justificación
Debido al avance industriala principios del siglo XX los componentes
estructurales como placas y cascarónes empezaron a ser recurrentes en su uso
dentro del sector petroqúımico, de generación de enerǵıa, en la elaboración
de recipientes a presión y tubos transportadores de vapor. En tiempos más
recientes su estudio se ha extendido a la turbomaquinaria ya que los álabes de
turbinas se enfrentan a corrientes de gases a elevadas temperaturas e intensos
ciclos de carga[40] .
Figura 1: Grieta en la ráız de un alabe de turbina debido al fenómeno del
creep
v
Figura 2: Tuberia sujeta a presión por la cual circula vapor a altas tempera-
turas favoreciendo la aparición del fenómeno del creep
Paulatinamente fueron desarrollándose los primeros estudios sobre el fenómeno
del creep, Andrade propuso un modelo logaŕıtmico en [16]. Los estudios pos-
teriores se enfocaron en las etapas primaria y secundaria a través de las
ecuaciones constitutivas propuestas por Norton-Bailey la cual describe un
modelo de deformación de tipo potencial el cual vaŕıa según el esfuerzo apli-
cado.
La teoŕıa del creep se ha ido nutriendo en su mayor parte por pruebas ex-
perimentales en las cuales se han recopilado datos como lo son tiempo a la
fractura del espećımen ya sea variando el esfuerzo o la temperatura, el com-
portamiento de la deformación y velocidad de deformación respecto al tiempo
y el comportamiento de la velocidad de deformación respecto al esfuerzo y
temperatura a la que se lleva a cabo el análisis experimental construyen-
do aśı familias de curvas que ayudan a comprender el comportamiento de
las deformaciones por creep en diferentes materiales. La implementación del
método del elemento finito en la teoŕıa de no linealidad material ha sido un
gran respaldo, capaz de crear un escenario del campo de deformaciones de-
pendientes del tiempo bastante aproximado al que se desarrollaŕıa a través
del método experimental, haciendo más rápido el proceso de obtención de
datos con la intención de saber el comportamiento del creep a diferentes
rangos de esfuerzo y temperatura en diversos elementos estructurales ma-
nufacturados con diferentes materiales. Aun cuando el método del elemento
finito no lineal constituye una gran ayuda al tedioso y largo procedimiento
vi
experimental, en cuanto al ahorro de un periodo conciderable de tiempo (las
pruebas experimentales pueden tomar incluso años para arrojar informaćıon
suficiente) la cantidad de simulaciones que se tienen que correr siguen siendo
conciderable.
Es por ello que recopilando información de un número prudente de simula-
ciones , se pueden generar patrones que nos permitan adquirir información
confiable del creep sobre otras posibles condiciones de esfuerzo que no se
llevaron a cabo por elemento finito pero que nos pueden arrojar curvas de
εcr vs tiempo bastante aproximadas a las que se pudieran obtener por expe-
rimentación o simulación numérica.
vii
Objetivos
Con los datos arrojados de pruebas del creep realizadas a través del méto-
do del elemento finito. Encontrar la relación que nos describe el esfuerzo que
se desarrolla en un punto espećıfico una placa deacuerdo a la tensión aplicada
a la misma , para aśı poder ingresar de manera directa valores de esfuerzo a
las las ecuaciones de Norton − Bailey con el propósito de conocer la curva
de deformación contra tiempo a esfuerzo constante y variable en la prime-
ra etapa del creep, sin la necesidad de una nueva simulación por medio de
software o prueba experimental.
viii
Índice
1. Conceptos Básicos del Creep 1
1.1. Ensayo de Tensión . . . . . . . . . . . . . . . . . . . . . . . . 3
1.1.1. Curvas del Creep . . . . . . . . . . . . . . . . . . . . . 5
1.2. Dependencia de la velocidad estacionaria del Creep con el es-
fuerzo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2. Elementos del Análisis Numérico 15
2.1. Clasificación de las Ecuaciones Diferenciales Parciales . . . . . 15
2.2. Condiciones Iniciales y de Frontera . . . . . . . . . . . . . . . 15
2.3. Técnicas de Discretización del Medio Continuo . . . . . . . . 16
2.4. Introducción al Método del elemento Finito . . . . . . . . . . 18
2.5. Formulación Variacional del M.E.F . . . . . . . . . . . . . . . 19
2.6. Elementos Isoparámetricos . . . . . . . . . . . . . . . . . . . . 22
2.7. Requisitos de Convergencia . . . . . . . . . . . . . . . . . . . 23
2.8. Elementos 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.9. Cuadrilátero Isoparamétrico cuadrático . . . . . . . . . . . . . 26
2.10. Transformación Isoparamétrica . . . . . . . . . . . . . . . . . 27
2.11. Ecuaciones Constitutivas de la Mecánica de sólidos . . . . . . 29
2.11.1. Deformaciones unitarias . . . . . . . . . . . . . . . . . 31
2.11.2. Campo de esfuerzos . . . . . . . . . . . . . . . . . . . . 32
2.12. Enerǵıa Potencial Total del Sistema Elástico . . . . . . . . . . 33
3. Comportamiento no lineal 37
3.1. Análisis del Movimiento . . . . . . . . . . . . . . . . . . . . . 39
3.2. El Gradiente de Deformación . . . . . . . . . . . . . . . . . . 40
3.3. Tensor Gradiente de Desplazamientos . . . . . . . . . . . . . . 42
3.4. Descomposición Polar del Tensor Gradiente de Deformación . 43
3.5. Variación de los tensores X t0 y H0 . . . . . . . . . . . . . . . . 43
3.6. Tensor infinitesimal de deformaciones unitarias . . . . . . . . . 44
3.7. Tensor de Deformaciones Unitarias de Green− Lagrange . . . 45
3.7.1. Variación del tensor de Green− Lagrange . . . . . . . 47
3.8. Modelo matemático del Creep . . . . . . . . . . . . . . . . . . 48
3.9. Métodos de integración del tiempo . . . . . . . . . . . . . . . 50
3.10. Implementación numérica . . . . . . . . . . . . . . . . . . . . 53
3.11. Métodos Iterativos . . . . . . . . . . . . . . . . . . . . . . . . 54
3.11.1. Método de Newton Raphson . . . . . . . . . . . . . . . 54
ix
3.11.2. Método de Newton − Raphson para un sistema de n
ecuaciones no lineales . . . . . . . . . . . . . . . . . . . 55
3.11.3. Método de Newton-Raphson Modificado . . . . . . . . 59
4. Aplicación 60
4.1. Resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5. Conclusiones 95
5.1. Propuestas para trabajos futuros . . . . . . . . . . . . . . . . 97
x
Índice de figuras
1. Grieta en la ráız de un alabe de turbina debido al fenómeno
del creep . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v
2. Tuberia sujeta a presión por la cual circula vapor a altas tem-
peraturas favoreciendo la aparición del fenómeno del creep . . vi
3. Variación de la longitud con el tiempo de una probeta sometida
a ceep . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
4. Curvas de fluencia a alta y baja temperatura Tf se refiere a la
temperatura de fusión . . . . . . . . . . . . . . . . . . . . . . 3
5. Curva de tensión t́ıpica de un material dúctil . . . . . . . . . . 4
6. Curvas convencionales del creep con temperaturas menores a
0.4Tf . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
7. Ensayo del creep en los primeros instantes de tiempo . . . . . 6
8. Curva de ensayo a tensión . . . . . . . . . . . . . . . . . . . . 6
9. Velocidades de deformación por creep a temperaturas menores
que 0.4Tf . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
10. Etapas del creep representadas en la curva εcr vs tiempo . . . 9
11. Creep a T o constante y diferentes magnitudes de σ . . . . . . 10
12. La pendiente de ln(ε̇cr) vs ln(σ) representa el coeficiente de
endurecimiento por deformación n . . . . . . . . . . . . . . . . 11
13. Descripción gráfica erronea del creep a esfuerzo variable . . . . 12
14. Deformación obtenida en curvas distintas en un tiempo t = tr 12
15. Misma deformación,diferentes valores de tiempo en curvas dis-
tintas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
16. Comparación time hardening vs strain hardening a σ variable 14
17. Armadura un tipo de estructura discrteta . . . . . .. . . . . . 17
18. Elemento placa ejemplo de estructura continua . . . . . . . . . 18
19. Métodos de discretización en sistemas continuos . . . . . . . . 18
20. Membrana en flexión . . . . . . . . . . . . . . . . . . . . . . . 19
21. Sistema de coordenadas naturales (ξ, η) isoparemétricas . . . 23
22. Representación gráfica de la transformación de coordenadas . 23
23. Elemento 2D con 8 nodos . . . . . . . . . . . . . . . . . . . . 25
24. Configuraciones a diferentes instantes de tiempo inicial, actual
e incrementada . . . . . . . . . . . . . . . . . . . . . . . . . . 39
25. Longitud y rotación de una fibra material en t, obtenidas a
partir de X t0 y posición t = 0 . . . . . . . . . . . . . . . . . . . 40
xi
26. Estado del cuerpo en el instante t = 0 y t = t , se muestran
las coordenadas deformadas xt1, x
t
2 . . . . . . . . . . . . . . . . 41
27. Paraleṕıpedo de lados diferenciales . . . . . . . . . . . . . . . 41
28. Rotación y alargamiento componentes de la deformación total 43
29. Distancia al cuadrado entre puntos vecinos en t = t y t+ ∆t . 48
30. Método exṕıcito de Euler . . . . . . . . . . . . . . . . . . . . 51
31. Método impĺıcito de Euler . . . . . . . . . . . . . . . . . . . . 52
32. Convergencia en el método simple de Newton−Raphson . . . 56
33. Esquema gráfico del método completo de Newton−Raphson 58
34. Comparación de Newton−Raphson completo y modificado . 59
35. Geometŕıa de la placa sometida a diferentes σ0[MPa] . . . . . 60
36. Dominio Ωplaca discretizado . . . . . . . . . . . . . . . . . . . 61
37. Distribución de σyy en la placa . . . . . . . . . . . . . . . . . 62
38. σ0i vs σyy en el nodo 11 con m = 1.7496 . . . . . . . . . . . . 63
39. Nodos agrupados a través de rutas . . . . . . . . . . . . . . . 64
40. σ0i vs σyy en los nodos de la ruta 1 . . . . . . . . . . . . . . . 72
41. Curva εcr,yy vs t en el punto (100[mm], 100[mm]) a σ0 =
150[MPa] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
42. Curva εcr,yy vs t en el punto (100[mm], 100[mm]) a σ0 =
175[MPa] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
43. Curva εcr,yy vs t en el punto (100[mm], 100[mm]) a σ0 =
200[MPa] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
44. ε̇cr,yy aplicando el criterio de time − hardening en el punto
(100[mm], 100[mm]) . . . . . . . . . . . . . . . . . . . . . . . 79
45. ε̇cr,yy aplicando el criterio de time − hardening en el punto
(100[mm], 100[mm]) . . . . . . . . . . . . . . . . . . . . . . . 80
46. ε̇cr,yy aplicando el criterio de time − hardening en el punto
(100[mm], 100[mm]) . . . . . . . . . . . . . . . . . . . . . . . 81
47. ε̇cr,yy aplicando el criterio de strain − hardening en el punto
(100[mm], 100[mm]) . . . . . . . . . . . . . . . . . . . . . . . 83
48. ε̇cr,yy aplicando el criterio de strain − hardening en el punto
(100[mm], 100[mm]) . . . . . . . . . . . . . . . . . . . . . . . 84
49. ε̇cr,yy aplicando el citerio de strain − hardening en el punto
(100[mm], 100[mm]) . . . . . . . . . . . . . . . . . . . . . . . 85
50. Curva εcr,yy vs t en el punto (100[mm], 0) a σ0 = 150[MPa] . . 86
51. Curva εcr,yy vs t en el punto (100[mm], 0) a σ0 = 175[MPa] . . 87
52. Curva εcr,yy vs t en el punto (100[mm], 0) a σ0 = 200[MPa] . . 88
xii
53. ε̇cr,yy aplicando el criterio de time − hardening en el punto
(20[mm], 0) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
54. ε̇cr,yy aplicando el criterio de time − hardening en el punto
(20[mm], 0) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
55. ε̇cr,yy aplicando el criterio de time − hardening en el punto
(20[mm], 0) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
56. ε̇cr,yy aplicando el criterio de strain − hardening en el punto
(20[mm], 0) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
57. ε̇cr,yy aplicando el criterio de strain − hardening en el punto
(20[mm], 0) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
58. ε̇cr,yy aplicando el criterio de strain − hardening en el punto
(20[mm], 0) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
59. Diagrama de flujo para utilizar los datos obtenidos de las
Tablas 1 a 11 y conocer el comportamiento del creep en cual-
quier punto (x, y) de la placa , siempre y cuando dicho punto
esté incluido en la base de datos . . . . . . . . . . . . . . . . . 96
xiii
Lista de Śımblos
εn Deformación nominal
εp Deformación plástica
εe Deformación elástica
εcr Deformación por creep
εc Deformación real
ε̇cr Velocidad de deformación por creep
εtotal Deformación total
l Longitud final
∆l Incremento de longitud
∆l0 Incremento de longitud instantánea
∆lc Incremento de longitud debido al creep
t Tiempo
∆t Incremento de tiempo
T 0 Temperatura
Tf Temperatura de fusión
F Fuerza
σn Esfuerzo normal
σing Esfuerzo ingenieŕıl
σY Esfuerzo de cedencia
Ac Área inicial
α1, α2, Constantes del material
A Coeficiente del material
n Coeficiente de endurecimiento por deformación
m Coeficiente de sensibilidad a la deformación
εAcr Deformación por creep en un punto en la curva A
εBcr Deformación por creep en un punto en la curva B
∆εcr Incremento en la deformación del creep
Lu Operador Laplaciano
Ω Dominio de un medio continuo
∂Ω Frontera de un medio continuo
R2 Espacio en dos dimensiones en el plano real
n̂ Vector normal
∇ Operador gradiente
W Trabajo
δ Operador variacional
xiv
I Funcional de enerǵıa
ξ, η Sistema natural de coordenadas isoparamétricas
δij Delta de Kronecker
Γ Contorno de la membrana en flexión
Φ Función forma
ν Relación de Poisson
E Módulo de elasticidad
σij Esfuerzo normal y cortante
J Matriz jacobiana
û Vector de desplazamientos
u Desplazamientos
εij Deformaciones normales y cortantes
B Matriz de derivadas de funciones forma en el análisis lineal
D Matriz de elasticidad
λ, µ Coeficientes de Lamé
Wc Trabajo debido a fuerzas concentradas
Wsp Trabajo debido a fuerzas de superficie
Wcp Trabajo debido a fuerzas de cuerpo
K Matriz de rigidez
f Fuerzas del sistema discretizado
R Cargas externas del sistema
F Estado de esfuerzos en el cuerpo
τ tij Tensor de Cauchy en un tiempo t
v Volumen de un cuerpo
ρ Densidad
xi Coordenada de un punto previo a la deformación
xti Coordenada de un punto posterior a la deformación en un tiempo t
X t0 Tensor de rotación y alargamiento en las fibras del material
H0 Tensor gradiente de desplazamientos
Rt0 Matriz de rotación en las fibras del material
U t0 Matriz de alargamiento en las fibras del material
Ĥ Vector gradiente de desplazamientos
St0 2
0 Tensor de esfuerzos de Piola Kirchhoff
G0 Matriz de derivadas de funciones forma en el análisis no lineal
Wint Trabajo interno
Wext Trabajo externo
KD Matriz de rigidez tangente
xv
Hs Matriz hessiana
σ0i Esfuerzo remoto a iésimos ensayos
m Pendiente de una ĺınea recta
σyy Esfuerzo aplicado en un punto con dirección y
xvi
Caṕıtulo 1
1. Conceptos Básicos del Creep
En la selección de materiales que se ven sometidos a altas temperaturas de
servicio, son muchos los factores que han de ternerse presentes. Entre dichos
factores se encuentra el tema de los costos, la facilidad para trabajarse, su
baja densidad (para usos dentro de la industria aeroespacial) y su resistencia
a la expocisión a condiciones ambientales.
Cuando evaluamos la resistencia a la deformación y fractura de algun ma-
terial, que se somete a carga durante lapsos prolongados de tiempo a alta
temperatura debemos poner especial atención en la contribución al análisis
mecánico que aporta el fenómeno del creep[2].
Ya que existen múltiples situaciones prácticas en las que los materiales se
ven sometidos a cargas y expuestos a temperaturas en periodos distintos de
tiempo. Un ejemplo claro ocurre en los álabes de de turbina de gas, cuya
operación se lleva a cabo soportando flujos de gas procedentes de la combus-
tión y grandes cargas de tensión.
Podemos ilustrar el fenómeno tomando una probeta de material metálico,
a la cual aplicamos un esfuerzo de tensión , dentro deun horno o cámara
que mantenga una temperatura constante y determinando estas condiciones
durante un periodo de tiempo, usualmente el tiempo de duración del ensayo
es de 1000hrs.
La direfencia con un ensayo de tensión a temperatura ambiente, es que en
este podemos apreciar la progresión de la deformación elástica εe conforme
aumenta la magnitud de los esfuerzos normales σn y que al descargarse la
probeta u objeto de estudio retoma su configuración original, mientras que
a mayores magnitudes de esfuerzo la deformación total será la suma de la
deformación elástica más la plástica εe + εp, recuperando parcialmente su
configuración inicial al efectuarse la descarga[11]. Existe además un ĺımite
de carga, que al superarse produce la fractura de la probeta, para estos ca-
sos la deformación producida solamente depende del esfuerzo a la que se ha
sometido. A más alta temperatura, aplicando las condiciones que carga y
temperatura elevada durante un lapso de tiempo[37]. Se podrá notar que la
longitud de la probeta crece gradualmente con el transcurrir del tiempo, es
1
decir la deformación producida depende del esfuerzo aplicado y el tiempo[13].
∆l = l − li = ∆lo + ∆lc (1.1)
Donde:
1. li Longitud inicial
2. l Longitud final
3. ∆lo Incremento de longitud instantánea que se produce al cargar la
probeta
4. ∆lc Incremento de longitud debido al creep
Tanto como ∆l como ∆lc vaŕıan con el tiempo y ∆lo es independiente de
este.
Figura 3: Variación de la longitud con el tiempo de una probeta sometida a
ceep
Hay que indicar que la forma de la curva anterior cambia su forma dea-
cuerdo a la temperatura, si es alta o baja. Naturalmente este concepto no es
el mismo para la gran gama de materiales disponibles en la industria. Este
fenómeno se puede ilustrar en la Figura 4
donde a bajas temperaturas, los cambios dimensionales debidos al creep
son extremadamente pequeños y la fractura rara vez ocurrirá. Por lo tanto
trabajaremos a altas temperaturas para que el creep resulte significativo y
exista el riesgo latente de fractura. En el caso de metales[22], la alta tempe-
ratura se concidera aquella magnitud cuyo valor parta de 0.4Tf (Tf se refiere
a la temperatura de fusión del material).
2
Figura 4: Curvas de fluencia a alta y baja temperatura Tf se refiere a la
temperatura de fusión
1.1. Ensayo de Tensión
Antes de describir la forma en la que el creep y la fractura debida a este
fenómeno son mesurables, discutiremos en que consiste y que se obtiene del
ensayo. Ya que esta prueba proporciona información concerniente entre el
esfuerzo aplicado y la deformación producida. Aśı, para un material ensaya-
do bajo las mismas condiciones, la carga necesaria para producir un cambio
depende de las condiciones de la pieza. Para definir las propiedades del ma-
terial, independientemente de su tamaño de la probeta, la fuerza frente al
alargamiento ∆l se convierte en gráficas de tensión ingenieŕıl σing frente a la
deformación siendo
σing =
F
A0
, ε =
∆l0
l0
(1.2)
El aspecto de una curva de tensión ingenieŕıl frente a la deformación para un
metal se muestra en la Figura 5.
Cuando las temperaturas de servicio son altas, las propiedades mecánicas
de los materiales debe determinarse cuidadosamente para que pueda garan-
tizarse que los componentes no se deformarán excesivamente[35]. Un requeri-
miento usual de diseño es lograr que el componente no sufra de deformación
plástica para las cargas que se planean aplicar durante el servicio. En este
caso, la sección transversal del elemento debe ser lo suficientemente grande
para que los esfuerzos aplicados no excedan a σY a la temperatura de traba-
jo.
Los ensayos de tensión convencionales proporcionan información necesaria
para que desde la fase de diseño pueda evitarse la deformación plástica o
3
Figura 5: Curva de tensión t́ıpica de un material dúctil
rotura. Sin embargo para largas condiciones de operación debe tomarse en
cuenta
1. Que la deformación debida al creep no provoque excesiva distorsión
durante la vida de servicio deseada
2. Que la fractura por creep no ocurra durante el servicio
El ensayo del creep conciste en representar la evolución temporal de la de-
formación de una probeta sometida a una carga y temperatura constante.
Hay que mencionar que, aunque la carga aplicada permanezca constante, el
esfuerzo real en la probeta aumenta a medida que esta se deforma y disminu-
ye su sección transversal. Conciderémos una probeta de longitud inicial l0 ,
sección transversal A0, que fluye bajo una carga total aplicada F , hasta una
nueva longitud l(l ≥ l0), y una nueva sección transversal A(A ≤ A0), después
de un tiempo t.
Admitiendo que el volumen de la probeta permanece constante y la defor-
mación es uniforme, es decir
l0A0 = lA (1.3)
entonces, el esfuerzo real que actúa sobre la probeta es
σ0 =
F
A
=
F
A0
l
l0
=
F
A0
(1 + εn) = σn(1 + εn) (1.4)
donde σn es el esfuerzo nominal en la probeta al comienzo del ensayo, y εn
es la deformación dada por (l−l0)
l0
.
Si la longitud de la probeta se incrementa desde l0 a l, durante un tiempo t, en
4
intervalos de tiempo dt resulta lógico pensar que también existirán cambios
diferenciales de longitud dl. Aśı el cambio de longitud de la probeta será l+dl
y no dl
l0
, ya que la longitud de la probeta en el instante t es l y no l0.
Entonces la deformación real, ε0 en la probeta se puede expresar como
ε0 =
∫ ε
0
dε =
∫ l
l0
= Ln
l
l0
(1.5)
Aunque se han desarrollado máquinas capaces de realizar ensayos a tensión
constante, la mayoŕıa de estos se llevan a cabo bajo carga constante, mucho
más sencillos de implementar. Esto es debido a la que si la ductilidad en creep
del material es baja (εcr ≤ 1 %) , los resultados obtenidos con carga constante
son indistinguibles. En este caso, además, no existe diferencia significatia
entre conciderar deformaciones nominales o reales.
1.1.1. Curvas del Creep
Al propiciar las condiciones que favorecen la operación de deformación
por creep, se observa en general que las formas de las curvas de εcr vs t que
en la Figura 6 se muestra la estructura t́ıpica de estas curvas de deformación
en rangos de T 0 ≤ 0.4Tf .
Figura 6: Curvas convencionales del creep con temperaturas menores a 0.4Tf
Observando las curvas vemos que inmediatamente después de la defor-
mación inicial debida a la carga aplicada e influencia de la temperatura, la
velocidad del creep decae continuamente. Bajo condicones de baja tempera-
tura , la deformación total por creep es normalmente muy baja, menor del
5
1 %. La deformación total será entonces
εtotal = εe + εcr (1.6)
donde εe es la deformación elástica instantánea que se produce al cargar el
espécimen y εcr es aquella deformación debida al creep.
Figura 7: Ensayo del creep en los primeros instantes de tiempo
Figura 8: Curva de ensayo a tensión
La curva σ vs ε de (7) para un determinado material vaŕıa con la tempe-
ratura, por lo que el valor de εe no sólo depende de la tensión aplicada sino
también de la temperatura. Esta dependencia se describe matemáticamente
como
ε = f1(σ, T
0) (1.7)
Lo cual indica que ε es función del esfuerzo y la temperatura. Sin embargo, la
deformción por creep es una función competente al esfuerzo y temperartura,
6
involucrará a la variable tiempo, por lo que modificando (1.7) conciderando
el tiempo
ε = f2(σ, T
0, t) (1.8)
Por ejemplo, la deformación del creep vaŕıa con el tiempo y la curva de ε vs t
cambia deacuerdo a al valor del la temperatura y a la magnitud del esfuerzo
aplicado como se puede ver en la Figura 6. La forma exacta de la expresión
matemática que satisface a (1.8), se calcula analizando las curvas del creep
obtenidas a diferentes esfuerzos, en ensayos realizados a baja temperatura.
Se suelen asumir diversas formas de (1.8) que separan la dependencia de la
deformación producida durante el creep con el esfuerzo, la temperatura y el
tiempo.
ε = fa(σ)fb(T
0)fc(t) (1.9)
En condicionesde baja temperatura, para muchos materiales cristalinos se ha
encontrado que la deformación del creep aumenta linealmente con el tiempo
de tal forma que
ε = α1Ln(α2 + 1) (1.10)
α1 y α2 son constantes para cada material que dependen del esfuerzo y de
la temperatura. Derivando (1.10) se obtiene la velocidad del creep en un
instante
ε̇cr =
dε
dt
=
α1α2
α2(t+ 1)
(1.11)
Expresión que se basa en el hecho que en ensayos a baja temperatura, la
velocidad del creep decrece continuamente al transcurrir el tiempo (9).
Entonces la ecuación (1.11) proporciona una descripción cuantitativa de
la dependencia de la deformación con el tiempo, por medio de α1 y α2 con el
esfuerzo y la temperatura.
Las ecuaciones que describen la dependencia de la deformación al creep con el
tiempo, el esfuerzo y la temperatura se denominan relaciones constitutivas[44].
Una vez que se conocen estas ecuaciones, se pueden conocer parámetros de
diseño como el esfuerzo máximo que el material puede experimentar a la T 0
de servicio sin alcanzar la deformación ĺımite del creep. Sin embargo, las le-
yes constitutivas exactas se requieren también para propósitos teoŕıcos. Un
7
Figura 9: Velocidades de deformación por creep a temperaturas menores que
0.4Tf
ejemplo del creep a bajas temperaturas se produce por el movimiento de dis-
locaciones dentro de la red cristalina. Un planteamiento sobre dislocaciones
debe explicar el comportamiento del creep a baja temperatura y tendrá que
predecir la dependencia logaŕıtmica que parte de εcr con t (1.10).
La ecuación (1.10) proporciona entonces una buena representación de los
cambios en la deformación con el tiempo para baja temperatura. Pero cuan-
do la temperatura incremente por encima de 0.4Tf , la curva del creep se
desv́ıa crecientemente de la curva logaŕıtmica. A altas temperaturas la de-
formación puede ser grande y el proceso generalmente culmina en ruptura.
Aśı que es prioritario desarrollar relaciones constitutivas correctas que des-
criban plenamente el comportamiento del material bajo condiciones de creep
a altas temperaturas.
Es importante mencionar que no se ha alcanzado un acuerdo general sobre
las ecuaciones que se deben utilizar para describir la forma de las curvas
del creep a altas temperaturas. Se pueden encontrar numerosas fuentes bi-
bliográficas que proponen diferentes modelos de comportamiento, por esta
razón, a alta temperatura, las curvas de creep se describen frecuentemente
haciendo referencia a tres etapas distintivas.
Como puede observarse, en (10) la deformación inicial que se produce a
cargar el espécimen provoca una deformación que va decreciendo al transcu-
rrir el tiempo a esto se le conoce como creep transitorio el cual conforma
la primera etapa. Después se llega a una etapa de estabilidad en la que se
igualan los procesos de endurecimeinto con procesos de recuperación, por lo
que la velocidad de deformación permanece constante. A esta epata se le
conoce como creep estacionario la cual corresponde a la segunda etapa de
8
Figura 10: Etapas del creep representadas en la curva εcr vs tiempo
la curva. Posteriormente vuelve a producirse un incremento en la velocidad
de deformación, consecuencia de la disminución de la sección transversal y
deterioro del material instantes antes de la fractura. Es aśı como se describe
la tercer etapa de la curva y en general la evolución de la curva con el tiempo
en términos de cambios de velocidad ya sea que esta se exprese como
ε̇cr = f3(σ, T
0, t) (1.12)
o
ε̇cr = f3(σ, T
0, εcr) (1.13)
En la mayoŕıa de los ensayos de creep, en especial en los de mayor dura-
ción, la segunda etapa es la predominante. Además, la velocidad durante
esta etapa puede utilizarse en relaciones emṕıricas , para determinar el tiem-
po a la fractura en función del material ensayado, por eso se acepta que lo
más importante es conocer la velocidad de la etapa estacionaria. Esta con-
cideración simplifica el problema de obtener ecuaciones que representen las
propiedades del creep en materiales a alta temperatura. Como la velocidad
de deformación permanece constante en esta etapa se dice que
ε̇ = f5(σ, T
0) (1.14)
es decir la variable tiempo y εcr son despreciables. Habitualmente suelen
separarse las dependencias con el tiempo y la deformación al contrario de lo
que expresan (1.12) y (1.13) por lo que (1.14) se reescribe
ε̇cr = u(σ)v(T
0) (1.15)
9
donde la función u(σ) describe la variación de ε̇cr con el esfuerzo y la función
v(T 0) lo hace con la temperatura.
1.2. Dependencia de la velocidad estacionaria del Creep
con el esfuerzo
Una vez que se ha encontrado una expresión satisfactoria para la función
v(T 0) seŕıa posible definir una forma conveniente para
ε̇crα(σ) (1.16)
Para ello se puede realizar ensayos a temperatura constante con diferentes
valores de esfuerzo. El tipo de curvas que se obtendrán se representan a través
de la Figura 11.
Figura 11: Creep a T o constante y diferentes magnitudes de σ
Desafortunadamente, no existe una única expresión para u(σ). Suele acep-
tarse una dependencia potencial del tipo α̇σn, esta se conoce como Ley de
Norton la cual porporciona la base para una ecuación que ha sido muy uti-
lizada para describir el comportamiento del creep a altas temperaturas.
La determinación de la pendiente Ln(ε̇cr) frente a Ln(σ) da el valor del
exponente n el cual es mejor conocido como coeficiente de endurecimiento
por deformación. Cuando la variación de ε̇cr con el esfuerzo se describe usan-
do la Ley potencial la ecuación queda de la siguiente manera.
εcr = σ
n (1.17)
10
Figura 12: La pendiente de ln(ε̇cr) vs ln(σ) representa el coeficiente de en-
durecimiento por deformación n
En el diseño de integridad estructural uno de los detalles más importantes
a tratar es el desempeño de la primera y segunda etapa como hemos venido
mencionando con bastante énfasis. Cuando m = 1 la Ley potencial describe
el comportamiento de las deformaciones por creep en la etapa estacionaria
justo como lo representa (1.17). Por lo que valores de m menores a uno re-
presentan las deformaciones en la etapa transitoria .
A (1.18) se les conoce como endurecimiento por tiempo, dondeA[N/mm2]nh−1
representa el coeficiente de difusión y m es el coeficinete de sensibilidad a la
deformación constantes del material que dependen de la temperatura y el
esfuerzo.
εcr = Aσ
ntm (1.18)
Derivando (1.18) con respecto a t obtenemos la velocidad de deformación
ε̇cr = mAσ
nt(m−1) (1.19)
Esto nos es útil para describir el comportamiento del creep en la primera eta-
pa, conciderando un ensayo uniaxial a esfuerzo y temperatura constante[10].
Si quiseramos conocer el comportamiento del creep a esfuerzo variable, man-
teniendo la condición de carga uniaxial, no podŕıamos aplicar (1.18) ya que no
contempla cambios en el esfuerzo, de hacerlo se obtendŕıa la gráfica descrita
en la Figura 13
Para obtener las expresiones que nos muestren correctamente los patro-
nes del comportamiento del creep a esfuerzo variable tenemos que despejar el
tiempo de (1.18) , teniendo otra aproximación de la velocidad de deformación
por creep. (1.20)representa el endurecimiento por deformación que incluye a
σ y εcr como variables.
11
Figura 13: Descripción gráfica erronea del creep a esfuerzo variable
ε̇cr = mA
1
mσ
n
m ε
(m−1)
n
cr (1.20)
Teniedo dos curvas εcr vs t cada una generada por un esfuerzo σ1 y σ2,
donde σ2 ≥ σ1, seleccionamos sobre el eje t un valor de tiempo que llamare-
mos tr[36]. Posteriormente observamos que en la curva de σ2 a tr se obtiene
una deformación εAcr, lo mismo con la curva de σ1 a un tr nos da una defor-
mación εBcr.
Figura 14: Deformación obtenida en curvas distintas en un tiempo t = tr
Después de dicha observación, ahora nos enfocaremos a encontrar el tiem-
po que tiene que transcurrir en la curva de σ2 para obtener una deformación
12
igual a εAcr, pudiéndose conseguir este valor de tiempo t1 por inspección co-
mo se puede ver en la Figura 15.Si resulta dif́ıcil encontrar t1 , basta con
despejar el tiempo de (1.18) con un σ = σ2 y εcr = ε
B
cr.
t1 = e
ln(εAcr)−ln(A)−nln(σ2)
m (1.21)
Figura 15: Misma deformación,diferentes valores de tiempo en curvas distin-
tas
De los diferentes valores de deformación εAcr y ε
B
cr en un tiempo t = tr
calculamos un ∆εcr
∆εcr = ε
A
cr − εBcr (1.22)
Aumentando ∆εcr a (1.18) obtenemos la expresión que nos describe el com-
portamiento de en las de formaciones según el criterio de endurecimiento por
tiempo (time− hardening) a esfuerzo variable.
εcr = Aσ
n
2 t
m −∆εcr (1.23)
Por medio de la Figura 15 y la ecuación (1.21) podemos conocer el tiempo
requerido para tener la misma deformación en la curva σ2 y en la curva de
σ1. Entonces podemos calcular el valor de ∆t.
∆t = tr − t1 (1.24)
Lo que nos conduce a predecir el comportamiento en las deformaciones por
endurecimiento por deformación (strain − hardening) a esfuerzo variable,
13
llegando a la siguiente expresión.
εcr = Aσ
n
2 (t−∆t)m (1.25)
Ambos criterios son utilizados en la práctica para predecir las curvas del
creep en condiciones de esfuerzo variable utilizando datos obtenidos en ensa-
yos a esfuerzo constante. La Figura 16 muestran el aumento en εcr debido a
un cambio en la magnitud del esfuerzo aplicado en un espacio de tiempo de-
terminado de manera correcta. Entonces podemos percatarnos que el criterio
de endurecimiento por deformación representa de una mejor aproximación
del comportamiento en las deformaciones a esfuerzo variable como se puede
ver en[10].
Figura 16: Comparación time hardening vs strain hardening a σ variable
14
Caṕıtulo 2
2. Elementos del Análisis Numérico
2.1. Clasificación de las Ecuaciones Diferenciales Par-
ciales
Los modelos de sistemas continuos están constituidos por sistemas de
ecuaciones diferenciales cuyo número es igual al número de propiedades in-
tensivas que intervienen en la formulación del modelo básico.
Los sistemas de ecuaciones diferenciales se clasifican en eĺıpticas, hiperbólicas
y parabólicas. Es necesario aclarar que esta clasificación no es exhaustiva; es
decir, existen sistemas de ecuaciones diferenciales que no pertenecen a nin-
guna de estas categoŕıas. Sin embargo, casi todos los modelos de sistemas
continuos, sobre todo aquellos que han recibido mayor atención, si están in-
cluidos dentro de alguna de estas categoŕıas.
2.2. Condiciones Iniciales y de Frontera
Dado un problema concreto de ecuaciones en derivadas parciales sobre
un dominio Ω, si la solución existe, esta no es única ya que generalmente
tiene un número infinito de soluciones. Para que el problema tenga una sola
solución es necesario imponer condiciones auxiliares apropiadas y estas son
las condiciones iniciales y de frontera[25].
Mencionaremos de manera general las C.I y C.F esenciales para definir un
problema de ecuaciones diferenciales: Las condiciones iniciales expresan el
valor de la función en el tiempo inicial t = 0
u(x, y, 0) = γ(x, y) (2.1)
Las condiciones de frontera especifican los valores que la función u(x, y, t) o
∇u(x, y, 0) tomarán en la frontera ∂Ω, siendo de tres tipos posibles:
15
1. Condiciones tipo Dirichlet: Especifica los valores que la función u(x, y, t)
toma en la frontera ∂Ω
u(x, y, t) = γ(x, y) (2.2)
2. Condiciones tipo Neumann: Aqúı se conoce el valor de la derivada de
la función u(x, y, t) con respecto a la normal n̂ a lo largo de la frontera
∂Ω
∇u(x, y, t)n̂ = γ(x, y) (2.3)
3. Condiciones tipo Robin: Esta condición es una combinación de las dos
anteriores
α(x, y)u(x, y, t) + β(x, y)∇u(x, y, t)n̂ = g(x, y)∀x, y ∈ ∂Ω. (2.4)
Una vez que se han planteado las ecuaciones que gobiernan el problema, las
condiciones iniciales , de frontera y mencionando los procesos que intervienen
de manera directa en el fenómeno estudiado, necesitamos que nuestro modelo
sea completo. Decimos que el modelo de un sistema es completo si define un
problema bien planteado. Un problema de valores iniciales y condiciones de
frontera es bien planteado si cumple que:
1. Existe una y sólo una solución y
2. La solución depende de manera continua de las condiciones iniciales y
de frontera del problema
Es decir, un modelo completo es aquél en el cual se incorporan condiciones
iniciales y de frontera que definen conjuntamente con las ecuaciones diferen-
ciales un problema bien planteado.
2.3. Técnicas de Discretización del Medio Continuo
Al efectuar una clasificación de estructuras, suelen dividirse en discretas,
reticulares y continuas. Las discretas son aquellas que están formadas por un
ensamblaje de elementos claramente diferenciados unos de otros y unidos en
una serie de puntos concretos, de tal manera que el sistema tiene forma de
malla o ret́ıcula. La caracteŕıstica fundamental de las estructuras discretas es
16
que su deformación puede definirse de manera exacta a través de un número
finito de parámetros, como lo son las deformaciones εij de los puntos de unión
de los elementos que conforman la estructura. De esta forma el equilibrio de
toda la estructura puede representarse a partir de las ecuaciones de equilibrio
que la conforman en las direcciones de dichas deformaciones.
Figura 17: Armadura un tipo de estructura discrteta
Comparativamente, en los sistemas continuos no es posible separar, a
priori, el sistema en un número finito de elementos estructurales discretos
si se toma una parte cualquiera del sistema, el número de puntos de unión
entre dicha parte y el resto del dominio Ω es infinito, por lo tanto es impo-
sible utilizar el mismo método que en los sistemas discretos, pues los puntos
de unión entre distintos elementos que ah́ı aparećıan de manera natural, no
existen en un sistema continuo. Las estructuras continuas son muy frecuentes
en ingenieŕıa, como lo son bastidores de máquinas, carroceŕıas de veh́ıculos,
losas de cimentación de edificios, bielas, ejes, vigas, placas y cascarones. Para
su análisis es necesario disponer de un método que tenga en cuenta su na-
turaleza continua. Los problemas del continuo pueden ser resueltos a través
de distintas técnicas de simulación numérica que se presentan a continuación
en el diagrama. En problemas lineales, el método del elemento finito domina
actualmente la escena en lo relativo a discretización espacial. El método del
elemento frontera constituyen una segunda alternativa en áreas de aplicacin
espećıficas. Para problemas no lineales, el dominio de los métodos por ele-
mento finito es enorme[12].
Los métodos por diferencias finitas aplicados en sólidos y mecánica de es-
tructuras han desaparecido virtualmente debido a que son poco prácticos.
17
Figura 18: Elemento placa ejemplo de estructura continua
Figura 19: Métodos de discretización en sistemas continuos
Sin embargo, para dinámica de fluidos estas técnicas son aún importantes.
Los métodos por volúmenes finitos, que se relacionan directamente con la
discretización de las leyes de conservación, son importantes en fluidos con un
elevado valor en su número de Reynolds.
Los métodos espectrales se basan en correspondencias que transforman di-
mensiones espaciales y/o temporales a espacios (por ejemplo, el dominio de
frecuencias), donde el problema es mucho más sencillo de resolver.
De reciente aplicación se tienen los métodos por mallado automático. Estos
combinan técnicas y herramientas del elemento finito, por ejemplo la formu-
lación variacional e interpolación, con caracteŕısticas de diferencias finitas.
2.4. Introducción al Método del elemento Finito
El método de los elementos finitos como ya se ha explicado es una técnica
numérica de aproximación de problemas gobernados por EDP en el espacio
f́ısico y temporal, de tal forma que: El medio continuo se divide en un número
finito de elementos, cuyo comportamiento se especif́ıca a través de un número
finito de parámetros asociados a ciertos puntos caracteŕısticos llamados no-
dos. Estos nodos son puntos de unión decada elemento con sus adyacentes.
La solución del sistema completo sigue las reglas de los problemas discretos.
18
El sistema completo se forma por ensamblaje de elementos. Las incógnitas
del problema dejan de ser funciones matemáticas y pasan a ser el valor de
estas funciones en los nodos. El comportamiento en el interior de cada ele-
mento queda definido a partir del comportamiento de los nodos a través de
funciones de interpolación mejor conocidas como funciones forma[3].
El M.E.F (método del elemento finito) entonces se basa en transformar un
dominio de naturaleza continua en un modelo discreto aproximado. El com-
portamiento de lo que suceda al interior del modelo aproximado se obtiene a
a partir de la interpolación de los valores conocidos en los nodos. Es entonces
una aproximación de los valores de una función a partir del conocimiento
de un número determinado y finito de puntos. Para fines de este trabajo y
visto desde un punto de vista estructural podemos entender al M.E.F como
una generalización del cálculo matricial de estructuras al análisis de sistemas
continuos. De hecho el método se gestó por evolución de aplicación a sistemas
estructurales[15].
2.5. Formulación Variacional del M.E.F
Se deducirá la funcional Enerǵıa para problemas de equilibrio o eĺıpti-
cos bajo condiciones de frontera homogéneas. Expĺıcitamente se calculará la
enerǵıa de un sistema f́ısico bajo condiciones de frontera homogéneas de la
forma que J.AOrtega lo hace en[18]. Para esto vamos a suponer que se tiene
una membrana o región Ω la cual se deflexiona a causa de una carga continua
distribuida en todo el dominio, con su frontera ∂Ω fija. Si u(x, y), representa
el desplazamiento del punto (x, y) a causa de F (x, y)
Figura 20: Membrana en flexión
19
entonces u(x(s), y(s)) a lo largo de todo (x(s), y(s)) ∈ ∂Ω. De esta manera
la ecuación diferencial parcial lineal que gobernará este fenómeno es
Lu = FenΩ (2.5)
Bu = gen∂Ω (2.6)
donde L = −∇2 el operador de Poisson y u(x, y) satisface las condiciones de
frontera de Dirichlet, Neumann y Robin en ∂Ω. Si F (x, y) es la fuerza loca-
lizada en el punto (x, y) y u(x, y) el desplazamiento en el mismo punto en el
cual se encuentra aplicada F , entonces W (x, y) = F (x, y)u(x, y) será el tra-
bajo realizado por la fuerza, de esta manera el trabajo total W desarrollado
por la flexión de la membrana es
W =
∫
Ω
F (x, y)u(x, y)dΩ (2.7)
W es en realidad una funcional, ya que para cada u, se tiene
W [u] =
∫
Ω
FudΩ (2.8)
De esta manera si tomamos la variación de W , tenemos
δW = δ
∫
Ω
FudΩ =
∫
Ω
F∆udΩ (2.9)
donde ∫
Ω
F∆udΩ =
∫
Ω
Lu∆udΩ = −
∫
Ω
∆u∇2udΩ (2.10)
y de la identidad vectorial
∇T (∆u∇u) = ∇T (∆u)∇u+ ∆u∇2u (2.11)
y de
∇T (∆u)∇u = ∆(1
2
‖ ∇u ‖2) (2.12)
se tiene ∫
Ω
F∆udΩ =
∫
Ω
∆(
1
2
‖ ∇u ‖2)dΩ−
∫
Ω
∇T (∆u∇u)dΩ (2.13)
20
Ahora el teorema de divergencia de Gauss se tendrá∫
Ω
∇T (∆u∇u)dΩ = −
∫
∂Ω
∆u(
∂u
∂n̂
)dS (2.14)
por lo que la integral queda como∫
Ω
F∆udΩ =
∫
Ω
∆(
1
2
‖ ∇u ‖2)dΩ−
∫
∂Ω
∆u(
∂u
∂n̂
)dS (2.15)
Donde ∂Ω = Γ1Γ2Γ3. En relación a la última integral se tiene
−
∫
∂Ω
(
∂u
∂n̂
)dS = −
∫
Γ1
∆u(
∂u
∂n̂
)dS −
∫
Γ2
∆u(
∂u
∂n̂
)dS −
∫
Γ3
∆u(
∂u
∂n̂
)dS (2.16)
y de las condiciones de frontera dadas por el problema, se tiene:
−
∫
Γ1
∆u(
∂u
∂n̂
)dS = 0Dirichlet
−
∫
Γ2
∆u(
∂u
∂n̂
)dS = 0Neumann
para el caso de la condición de frontera del tipo mixta se tiene.
−
∫
Γ3
∆u(
∂u
∂n̂
)dS =
∫
Γ3
∆u(σ(s)u)dS
=
∫
Γ3
σ(s)u(s)∆u(s)dS
=
∫
Γ3
σ(s)∆(
u2
2
)dS
= ∆
∫
Γ3
σ(s)u2
2
dS
Aśı la expresión de la variación W queda como
δW = ∆
∫
Ω
FudΩ =
∫
Ω
F∆udΩ (2.17)
donde ∫
Ω
F∆udΩ = ∆
∫
Ω
(
1
2
‖ ∇u ‖2)dΩ + ∆
∫
Γ3
σ
2
u2dS (2.18)
21
Entonces obtendremos
∆
∫
Ω
FudΩ = ∆
∫
Ω
‖ ∇u ‖2
2
dΩ + ∆
∫
Γ3
σ
2
u2dS (2.19)
y finalmente
∆[
∫
Ω
(‖ ∇u ‖2 −2Fu)dΩ +
∫
Γ3
σ
2
u2dS] = 0 (2.20)
Definiendo entonces a la funcional I(u) como
I(u) =
∫
Ω
(‖ ∇u ‖2 −2Fu)dΩ +
∫
Γ3
σ
2
u2dS (2.21)
De (2.21) podemos afirmar que la primera variación es nula
∆I(u) = 0
Esto permite asociar a las ecuaciones (2.5) y (2.6) con la funcional de enerǵıa
I(u). Veremos que para el caso en que el operador L sea eĺıptico, autoadjun-
to y definido positivo, existe una correspondencia biuńıvoca entre los puntos
cŕıticos de I(u) y las soluciones de (2.5) bajo (2.6). Siendo esta correspon-
dencia la que permite resolver a (2.5) y (2.6) en forma única mediante la
alternativa de determinar los puntos cŕıticos de I(u); es decir resolver la
ecuación funcional.
dI(u)
du
= 0 (2.22)
2.6. Elementos Isoparámetricos
Las coordenadas naturales isoparamétricas (ξ, η) de elementos cuadriláte-
ros que adoptaremos se muestran en la Figura 21
Las coordenadas (ξ, η) vaŕıan entre −1, 1 de un lado a otro del cuadriláte-
ro, tomando como valor nulo las ĺıneas isoparamétricas que unen los puntos
medios de lados opuestos.
Es bastante habitual representar elementos cuadriláteros de un sistema car-
tesiano de coordenadas (ξ, η), denominado espacio de referencia , en el que
todos los elementos se representan como cuadrados de lado 2 como lo mues-
tra la Figura 22.
22
Figura 21: Sistema de coordenadas naturales (ξ, η) isoparemétricas
Figura 22: Representación gráfica de la transformación de coordenadas
La transformación entre las coordenadas geométricas (x, y) y las naturales
viene dada por la transformación isoparamétrica.
x(ξ, η) =
n∑
i=1
xiΦi(ξ, η) (2.23)
2.7. Requisitos de Convergencia
Desde el punto de vista práctico la convergencia del MEF implica que.
Según como se vaya refinando la malla, la solución obtenida se aproxima
a la solución exacta del modelo matemático que se desea resolver[45]. Los
requisitos de convergencia se puede dividir en tres categoŕıas:
1. Complitud: Este requisito supone que las funciones de interpolación
utilizados tienen riqueza suficiente para capturar la solución exacta en
el ĺımite del refinamiento en el que el tamaño del elemento se aproxima
a 0.
23
2. Compatibilidad: Este requisito exige que las funciones de forma permi-
tan obtener una solución aproximada que tenga la continuidad exigida
entre los elementos de la malla que son adyacentes.
3. Estabilidad: Con este requisito se garantiza que las ecuaciones de ele-
mentos finitos están bien planteadas, teniendo las mismas propiedades
de unicidad que la solución exacta del problema se modela.
La idea básica de la formulación isoparamétrica es utilizar las mismas fun-
ciones de forma para interpolar la geometŕıa del elemento como el campo de
desplazamientos. La representación del concepto se muestra[30].
La generalización de la representación isoparamétrica de un elemento tridi-
mensional arbitrario es directa.
La ecuación (2.23) expresa la prioridad de partición de la unidad, asociada
al requisito de consistencia, las tres ecuaciones dentro del arreglo matricial
corresponden a la interpolación de la geometŕıa, y las tres últimas se refieren
a la interpolación del campo de desplazamientos.
Las funciones de forma Φ deben satisfacer los siguientes requisitos:
1. La función forma elementalΦi debe tomar el valor 1 en el nodo i y 0
en los demás nodos del elemento. Por lo que Φi evaluada en los nodos
j,k,l...n del elemento nos debe dar la matriz identidad δij
2. Φi debe ser nula en los contornos del elemento que no contenga al nodo
i
3. Las funciones de forma Φi están definidas en los elementos que com-
parten el nodo i deben ser continuas de clase C(m−1)
4. Las funciones forma Φi de un elemento deben expresar de forma exacta
cualquier polinomio de grado m en las coordenadas del elemento.
Como se verá en los siguientes apartados, para diseñar las funciones de forma
éstas se expresan como el producto de n factores Lj función de las coorde-
nadas naturales.
Φi = CiL1L2....Ln (2.24)
2.8. Elementos 2D
Las funciones forma de los elementos isoparamétricos en dos dimensiones
satisfacen los requisitos descritos anteriormente si se diseñande acuerdo con
24
los cinco subrequisitos siguientes aplicados al elemento que se presenta a
través de la Figura 23 .Ya que para este se deduciran sus funciones de
forma[23].
Figura 23: Elemento 2D con 8 nodos
1. Para la función forma Φi seleccionar Lj con el criterio que el nodo
contenga el mı́nimo número de ĺıneas isoparamétricas asociadas con la
expresión lineal en las coordenadas naturales que contenga a todos los
nodos del elemento exepto al nodo i
2. Calcular el valor del coeficiente Ci para la función de forma Φi valga 1
en i
3. Comprobar que la función forma Φi vale cero en todos los lados del
elemento que no coincidan con el nodo i
4. Comprobar el grado del polinomio correspondiente a particularizar Φ
en los lados del elemento que contengan al nodo i. Śı el número de
nodos en un lado del elemento es p , para que se valide el requisito de
compatibilidad el grado del polinomio en dicho lado debe ser p− 1
Una vez verificados los requisitos anteriores, comprobar finalmente que las
funciones de forma definidas en el elemento suman la unidad.
25
2.9. Cuadrilátero Isoparamétrico cuadrático
La deducción de la expresión de las funciones de forma para nodos situa-
dos en los vértices sugún la Figura 23 se detalla en el nodo 1 deacuerdo con
R− 1
Φ1(ξ, η) = C1L58L237 = C1(1 + ξ + η)(1− ξ)(1− η) (2.25)
Imponiendo el requisito de normalización C1 =
1
4
Φ1 = −
1
4
(1 + ξ + η)(1− ξ)(1− η) (2.26)
En los lados 374 y 263, definidos respectivamente por ξ = 1 , η = 1, es
evidentemente que Φi = 0 por lo que el requisito de compatibilidad se ve-
rificará śı al particularizar Φi en los lados 152 y 184 resultan polinomios de
segundo grado. Que corresponde a restar 1 al número de nodos en cada lado.
Para L152 con η = −1
Φ1(ξ) = −
1
2
ξ(1− ξ) (2.27)
y L184 con ξ = −1
Φ1(η) = −
1
2
ξ(1− ξ) (2.28)
Con éste mismo procedimiento aplicado en los demás vértices se obtiene:
Φ2 = −
1
4
(1− ξ + η)(1 + ξ)(1− η) (2.29)
Φ3 = −
1
4
(1− ξ − η)(1 + ξ)(1 + η) (2.30)
Φ4 = −
1
4
(1 + ξ − η)(1− ξ)(1 + η) (2.31)
Para las funciones forma correspondientes a los nodos situados en el punto.
Se detallan los cálculos para obtener Φ5.
Φ5 = C5L184L374L263 = C5(1 + ξ)(1− η)(1− ξ) (2.32)
Afirmando que Φ5 = 1 en (ξ = 0, η = −1) y C5 = 12
Φ5 =
1
2
(1− η)(1− ξ2) (2.33)
26
En L184, L374 y L263 es inmediato verificar que Φ5 = 0.
Para verificar la condición de compatibilidad la expresión de Φ5 en L152 ha
de ser un polinomio de grado 20. Las restantes funciones forma son:
Φ6 =
1
2
(1 + ξ)(1− η2) (2.34)
Φ7 =
1
2
(1− ξ2)(1 + η) (2.35)
Φ8 =
1
2
(1− ξ)(1− η2) (2.36)
La compatibilidad del elemento se verifica con el requisito de partición de la
unidad:
2.10. Transformación Isoparamétrica
Teniendo un punto dentro de nuestro elemento deacuerdo al sistema de
coordenadas (ξ, η) la expresión que nos relaciona éste punto con el sistema
global (x, y) es.
x =
8∑
i=1
xkΦk(ξ, η), y =
8∑
i=1
ykΦk(ξ, η) (2.37)
Donde cada (xk, yk) representa las coordenadas de cada nodo del k ésimo
elemento en el sistema global (x, y). Con los ocho nodos por elemento que
tenemos existen dos grados de libertad definidos por xk, yk lo que resulta un
total de dieciséis grados de libertad representando (2.37) en forma matricial
tenemos.
{
x
y
}
=
[
Φ1 0 · · · Φ8 0
0 Φ1 · · · 0 Φ8
]

x1
y1
·
·
·
x8
y8

(2.38)
Gracias a que las funciones Φi son continuas en C
0 también lo son entre
elementos adyacentes de la malla. Para calcular el jacobiano de ésta trans-
formación necesitamos las derivadas de las funciones forma.
∂Φ1(ξ, η)
∂ξ
=
1
4
(1− η)(2ξ + η) (2.39)
27
∂Φ1(ξ, η)
∂η
=
1
4
(1− ξ)(2η + ξ) (2.40)
∂Φ2(ξ, η)
∂ξ
=
1
4
(1− η)(2ξ − η) (2.41)
∂Φ2(ξ, η)
∂η
=
1
4
(1 + ξ)(2η − ξ) (2.42)
∂Φ3(ξ, η)
∂ξ
=
1
4
(1 + η)(2ξ + η) (2.43)
∂Φ3(ξ, η)
∂η
=
1
4
(1 + ξ)(2η + ξ) (2.44)
∂Φ4(ξ, η)
∂ξ
=
1
4
(1 + η)(2ξ − η) (2.45)
∂Φ4(ξ, η)
∂η
=
1
4
(1− ξ)(2η − ξ) (2.46)
∂Φ5(ξ, η)
∂ξ
= −ξ(1− η) (2.47)
∂Φ5(ξ, η)
∂η
= −1
2
(1− ξ2) (2.48)
∂Φ6(ξ, η)
∂ξ
=
1
2
(1− η2) (2.49)
∂Φ6(ξ, η)
∂η
= −η(1 + ξ) (2.50)
∂Φ7(ξ, η)
∂ξ
= −ξ(1 + η) (2.51)
∂Φ7(ξ, η)
∂η
=
1
2
(1− ξ2) (2.52)
∂Φ8(ξ, η)
∂ξ
= −1
2
(1− η2) (2.53)
∂Φ8(ξ, η)
∂η
= −η(1− ξ) (2.54)
Para el calculo de las derivadas cartesianas y de las integrales de volumen que
intervienen en la formulación de elementos finitos, es útil conocer la matriz
28
jacobiana de la transformación del espacio cartesiano al isoparamétrico que
es nuestro marco de referencia. Dicha transformación se expresa a través de
la regla de la cadena del cálculo diferencial[18].
∂Φi
∂ξ
=
∂Φi
∂x
∂x
∂ξ
+
∂Φi
∂y
∂y
∂ξ
(2.55)
∂Φi
∂η
=
∂Φi
∂x
∂x
∂η
+
∂Φi
∂y
∂y
∂η
(2.56)
Ensamblando el sistema de ecuaciones formado por (2.55) y (2.56) tendremos[
∂Φi
∂ξ
∂Φi
∂η
]
=
[
∂x
∂ξ
∂y
∂ξ
∂x
∂η
∂y
∂η
][∂Φi
∂x
∂Φi
∂y
]
(2.57)
Los elementos kl de la matriz jacobiana se calculan por medio de las siguien-
tes ecuaciones.
J11 =
∂x
∂ξ
=
8∑
i=1
xi
∂Φi
∂ξ
(2.58)
J12 =
∂y
∂ξ
=
8∑
i=1
yi
∂Φi
∂ξ
(2.59)
J21 =
∂x
∂η
=
8∑
i=1
xi
∂Φi
∂η
(2.60)
J22 =
∂y
∂η
=
8∑
i=1
yi
∂Φi
∂η
(2.61)
2.11. Ecuaciones Constitutivas de la Mecánica de sóli-
dos
Orientando la formulación isoparamétrica en problemas de elasticidad,
donde la posición de un punto está definida por las coordenadas (x, y) y su de-
formación por dos componentes u(x, y) y v(x, y). El campo de deformaciones[41]
es entonces un vector de la forma.
û = [u(x, y), v(x, y)] (2.62)
29
Recordando la expresión (2.37) que define a las coordenadas (x, y) y estás
definen a las componentes (u, v) de las deformaciones. Podemos expresar a
u, v como.
u =
n∑
k=1
Φkuk, v =
n∑
k=1
Φkvk (2.63)
De manera matricial û se puede expresar como el producto de la matriz Φ ,
la cual forma parte de la ecuación (2.38). Cuyo número de columnas se define
por el doble del número de nodos que contenga el elemento.
û = Φδe (2.64)
δe representa al vector de las deformaciones nodales de cada elemento que
conforma el sistema discreto. Su representación en forma vectorial es. δe =
[u1v1 · · · unvn]T
Ecuación fundamental de la cual se desprende todo un desarrollo para
calcular los campos incógnita que plantea el problema elástico.
Dentro de la teoŕıa elástica en dos dimensiones existen las condiciones de:
1. Esfuerzo Plano: Cuando el esfuerzo en dirección σzz localizado en sen-
tido perpendicular al plano xy es 0. Ya que el sólido puede desplazarse
libremente en el sentido de su espesor. Por lo tanto existe una defor-
mación unitaria εzz no nula en dicha dirección.
2. Deformación Plana: Sucede cuando no hay posibilidad de deformación
dirección al plano del espesor del sólido, es decir εzz = 0 por lo que se
genera un esfuerzo σzz diferente de 0.
En ambos casos el esfuerzo y la deformación en la dirección z del plano carte-
siano (x, y, z) no contribuye a la enerǵıa elástica del sistema. Las ecuaciones
diferenciales que rigen el problema son de m = 2 en las deformaciones, pues
contienen la derivada primera de los esfuerzos que a su vez son derivados
de las deformaciones. En la expresión del potencial total del sistema apare-
cerán las deformaciones unitarias, que son las derivadas de grado 10 de los
desplazamientos.
30
2.11.1. Deformaciones unitarias
Las deformaciones unitarias en un punto cualquiera del elemento, supo-
niendo una magnitud pequeña se expresan como.
ε =
εxxεyy
γxy
 =
 ∂u∂x∂u
∂y
∂u
∂y
+ ∂u
∂x
 (2.65)
De forma matricial εxxεyy
γxy
 =
 ∂∂x 00 ∂
∂y
∂
∂y
∂
∂x
[u
v
]
(2.66)
En ésta expresíıon se identifica el operador matricial ∂ el cual permite pasar
de las deformaciones de un punto aleatorio (u, v)a las deformaciones unita-
rias ε. Este operador tiene tantas filas como deformaciones unitarias haya en
el problema y tantas columnas como componentes tenga el campo de despla-
zamientos.
Substituyendo las deformaciones u, v en función de las deformaciones nodales
, a través de las funciones formaobtenemos.
ε = ∂u = ∂Φδe (2.67)
Esta expresión relaciona las deformaciones de los nodos que conforman un
elemento en particular (δe) con las deformaciones unitarias en un punto inte-
rior cualquiera dentro del dominio del elemento. Se le conoce como matriz B
la cual representa el campo de deformaciones unitarias que se supone existen
en el interior de cada elemento finito. Conociendo la estructura de la matriz
Φ, la submatriz Bi se puede expresar de la forma.
Bi =
∂Φi∂x 00 ∂Φi
∂y
∂Φi
∂y
∂Φi
∂x
 (2.68)
Para nuestro elemento elegido cuyo número de funciones forma ésta repre-
sentado por el sub́ındice i que corre de 1a8. Entonces nuestra matriz B
estará conformada por 8 submatrices del tipo Bi. El ensamble matricial
tendrá un tamaño de 3X16.
B =
∂Φ1∂x 0 · · · ∂Φ8∂x 00 ∂Φ1
∂y
· · · 0 ∂Φ8
∂y
∂Φ1
∂y
∂Φ1
∂x
· · · ∂Φ8
∂y
∂Φ8
∂x
 (2.69)
31
2.11.2. Campo de esfuerzos
Las direcciones del estado de esfuerzos en el plano se representa de la
siguiente manera.
σ =
σxxσyy
τxy
 (2.70)
Donde la ecuación constitutiva en ausencia de temperatura se representa por
la LeydeHook.
σ = Cε (2.71)
Para un material elástico lineal e isótropo el tensor de elasticidad C es cons-
tante. 
σxx
σyy
σzz
τxy
τyz
τzx
 =

λ+ 2µ λ λ 0 0 0
λ λ+ 2µ λ 0 0 0
λ λ λ+ 2µ 0 0 0
0 0 0 µ 0 0
0 0 0 0 µ 0
0 0 0 0 0 µ


εxx
εyy
εzz
γxy
γyz
γzx
 (2.72)
λ = Eν
(1+ν)(1−2ν) y µ =
E
2(1+ν)
son los coeficientes de Lamé. Suponiendo que
σzz = 0 se obtiene
λεxx + λεyy + (λ+ 2µ)εzz = 0 (2.73)
Se bebe calcular el valor de la deformación unitaria transversal al material.
εzz =
−λ
(λ+ 2µ)
(εxx + εyy) (2.74)
Substituyendo (2.74) en (2.72) que representa un estado de esfuerzos tridi-
mensional, obtenemos el tensor de elasticidad del estado plano.
C =
E
(1− ν2)
1 ν 0ν 1 0
0 0 (1−ν)
2
 (2.75)
Podemos obtener C para el caso de la deformación plana.
C =
E(1− ν)
(1 + ν)(1− 2ν)
 1 ν1−ν 0ν
1−ν 1 0
0 0 1−2ν
2−2ν
 (2.76)
32
2.12. Enerǵıa Potencial Total del Sistema Elástico
Concideremos un cuerpo plano que puede representarse a través de un
dominio bidimensional Ω discretizado por elementos finitos. La enerǵıa po-
tencial total I de un cuerpo elástico lineal se define como la suma de la enerǵıa
potencial de deformación A y la enerǵıa potencial asociada al trabajo de las
fuerzas externas W [18].
I = Wp + A (2.77)
Conciderando que el trabajo W realizado por las cargas es el negativo de la
enerǵıa potencial Wp.
I = A−W (2.78)
Partiendo el dominio Ω en elementos finitos
I =
nel∑
e=1
(Ae −We) =
nel∑
e=1
Ie (2.79)
W tiene la forma
W =
1
2
[
εxx εyy εzz γxy γyz γzx
]

σxx
σyy
σzz
τxy
τyz
τzx
 (2.80)
La expresión (2.79) está asociada a fenómenos de elasticidad lineal, para la
cual podemos enunciar el siguiente teorema de la enerǵıa potencial.
Teorema 1 De todos los desplazamientos que satisfacen a las condiciones
de frontera dadas, que satisfacen las ecuaciones de equilibrio, son aquellos
que hacen que la enerǵıa potencial tenga un valor estacionario
La enerǵıa de deformación para un elemento diferencial de volumen, viene
dado por
dA =
1
2
εTσ − 1
2
ε0σ (2.81)
Donde ε0 se refiere a la deformación unitaria en el estado inicial. Integrando
dA sobre todo el volumen obtenemos la energíıa de deformación total
A =
1
2
∫
Ω
[εTσ − 1
2
ε0σdΩ] (2.82)
33
Las expresiones que nos ayudarán a deducir las relaciones que dan forma a
la integral de la enerǵıa de deformación para cada elemento de la red son
(2.65),(2.80) y (2.68). Entonces para cada elemento de la red tenemos.
Ae =
1
2
∫
Ωe
uTBTe DBudΩ−
∫
Ωe
uTBTe Dε0dΩ (2.83)
Definiendo las componentes que contribuyen al trabajo W son las siguientes.
1. Wc Cargas concentradas
2. Wsp Cargas de superficie
3. Wcp Cargas de cuerpo
El trabajo realizado por Wc es el más sencillo de manejar, ya que en cada
punto donde se realice trabajo debido a las fuerzas concentradas, se locali-
zan nodos de la red. Entonces las cargas concentradas multiplicadas por el
desplazamiento expresan a Wc
Wc = u
T q (2.84)
Suponiendo que las fuerzas tienen componentes paralelas a los desplazamien-
tos.
Esta componente de trabajo no está incluida en la sumatoria de los elemen-
tos que conforman la región Ω, ya que las fuerzas han sido localizadas en los
nodos.
En lo que se refiere al trabajo realizado por las cargas distribuidas sobre la
superficie de cada elemento se presenta la ecuación. Donde u, v, w representan
las componentes del vector desplazamiento.
Wsp =
∫
∂Ω
(uP ex + vP
e
y + wP
e
z )dS (2.85)
De mejor forma (2.85) se puede representar por
Wsp =
∫
∂Ω
uT [Φe]T
PxPy
Pz
 dS (2.86)
El trabajo que se lleva a cabo por la acción de las fuerzas de cuerpo es
Wcp =
∫
Ωe
uT [Φe]T
xy
z
 dΩ (2.87)
34
Habiendo definido el trabajo que generan Wc, Wsp y Wcp definimos el trabajo
local mediante
We =
∫
Ωe
uT [Φe]T
xy
z
 dΩ + ∫
∂Ω
uT [Φe]T
PxPy
Pz
 dS + uT q (2.88)
Teniendo W y A podemos calcular I de manera local y global
Ie = Ae −We (2.89)
Ie =
nel∑
e=1
(Ae −We) = A−W (2.90)
La enerǵıa potencial local está representada por (2.89) y la global a través
de (2.90) . Donde su expresión completa queda como
Ie =
nel∑
e=1
[
1
2
∫
Ωe
uTBTe DBudΩ−
∫
Ωe
uTBTe Dε0dΩ (2.91)
−
∫
Ωe
uT [Φe]T
xy
z
 dΩ− ∫
∂Ω
uT [Φe]T
PxPy
Pz
 dS − uT qe
Derivando (2.91) con respecto a u e igualando a 0 obtenemos
∂I
∂u
=
∂
∂u
[
nel∑
e=1
Ie] =
nel∑
e=1
∂Ie
∂u
= 0 (2.92)
Desglosando (2.92) se puede observar que las derivadas quedan como
∂
∂u
(
1
2
∫
Ωe
uTBTe DeBeudΩ) =
∫
Ωe
BTe DeBeudΩ (2.93)
∂
∂u
(−
∫
Ωe
uTBTe DeBeε
e
0dΩ) =
∫
Ωe
BTe Deε
e
0dΩ (2.94)
∂
∂u
(−
∫
Ωe
uTΦTe
xeye
ze
 dΩ) = −∫
Ωe
ΦTe
xeye
ze
 dΩ (2.95)
35
∂
∂u
(−
∫
∂Ωe
uTΦTe
P exP ey
P ez
 dS (2.96)
.
Gracias a estas derivadas podemos obtener la matriz de rigidez de cada ele-
mento que conforma la malla.
Keu =
1
2
∫
Ωe
BTe DeBeudΩ (2.97)
Aśı como también podemos obtener las fuerzas locales con
−fe =
∫
Ωe
ΦTe
xeye
ze
 dΩ + ∫
Ωe
ΦTe
P exP ey
P ez
 dΩ + ∫
Ωe
BTe Deε
e
0dΩ (2.98)
La matriz de rigidez y el vector de fuerzas globales lo calculamos por medio
de
K =
nel∑
e=1
Ke (2.99)
f =
nel∑
e=1
fe (2.100)
Finalmente concluimos con la expresión de equilibrio entre las fuerzas exter-
nas e internas del elemento. Esta expresión será de gran utilidad en caṕıtulos
posteriores donde analizaremos métodos de convergencia de fenómenos no
lineales como es al que a esta tesis compete. Ya que dichas aproximaciones
por medio de iteraciones se linealizan para obtener
Ku = f (2.101)
36
Caṕıtulo 3
3. Comportamiento no lineal
A diferencia del comportamiento lineal en el sentido causa − efecto , al
existir un amumento x en los esfuerzos las deformaciones también variarán
conforme al x aumento en los esfuerzos nuestro modelo se representa como
xσ−xε. Los problemas fuera de éste dominio son clasificados como no lineales.
El comportamiento no lineal de un sólido o estructura puede ser caracterizado
a través de cuatro modos como lo hace K.J Bathe en[5].
1. No linealidad material
2. No linealidad geométrica
3. No linealidad de condiciones naturales de frontera
4. No linealidad en condiciones iniciales de frontera
En años recientes el método del elemento finito conciderando no linealidad
aplicado a problemas estáticos ha tenido un gran crecimiento en el área de
ingenieŕıa [12]. Numerosos paquetes computacionales son empleados en la
solución de éste tipo de problemas que en ocasiones se vuelven bastante
complejos.
Comparando la diferencia entre el comportamiento elástico y no-elástico de
los materiales radica en que el estado de esfuerzos σij puede ser evaluado
directamente de las deformaciones elásticas. Como lo indica la Ley de Hook.
σij = Eijkl(εe)kl
Mientras que en el comportamiento no-elástico, el estado de esfuerzos depen-
de del comportamiento de las deformaciones según el tipo de fenómenoque
se analice. T́ıpicamente se clasifican en:
1. Elasto-plasticidad
2. Viscoelasticidad
3. Viscoplasticidad
37
4. Creep
Cuya ecuación constitutiva en términos incrementales es la Ley de Hook ge-
neralizada
dσij = Eijkl(dεkl − (dεkl)n.e)
Donde (dεkl)n.e es la deformación no elástica según la clasificación anterior.
Desafortunadamente existen desventajas en un análisis de naturaleza no li-
neal , una de ellas puede ser el costo computacional en base al algoritmo
seleccionado para la solucionar algún problema. Es por ello que el usuario
debe tener amplios conocimientos y buen jucio al momento de ingresar datos
a un programa. Ya que sus concideraciones particulares deben aterrizar en
datos de salida confiables. De esta situación surge la necesidad de mejorar
los algoritmos ya existientes los cuales contengan una mayor exactitud y es-
tabilidad numérica[7].
Generalmente la obtención de una solución a un problema no lineal es más efi-
ciente definiendo incrementos en la carga o tiempo de modo que las variables
involucradas se irán actualizando en función de los incrementos definidos con
anterioridad como lo porpone K.J Bathe y CimentoArthur en[8]. Como
consecuencia obtendremos un patrón gráfico. Resulta de gran importancia
que las ecuaciones se satisfagan a cada incremento ya sea de carga o tiempo
con suficiente exactitud aśı se evita la acumulación de errores, lo que gene-
raŕıa inestabilidad numérica.
Las formulaciones Lagrangianas han sido ampliamente utilizadas en no li-
nealidad, autores como Moya J.F , Monleón S y Léger Sophie presentan en
[34] y [26] un tratamiento generalizado de estas descripciones adecuado para
analizar placas como elementos bidimensionales.
Estudiando la deformación de un sólido, tenemos que tomar en cuenta tres
configuraciones del cuerpo las cuales son:
1. Configuración Inicial: Cuerpo que se supone libre de esfuerzos y defor-
maciones
2. Configuración Actual t = t
3. Configuración aumentada t = t+ ∆t
38
Figura 24: Configuraciones a diferentes instantes de tiempo inicial, actual e
incrementada
3.1. Análisis del Movimiento
La descripción Lagrangiana requiere una configuración de referencia para
deducir ecuaciones de equilibrio. Esto significa que observaremos los sucesi-
vos estados de deformación que el cuerpo sufrirá desde una configuración
determinada que ocupa un espacio de tiempo t. Esta hipótesis persigue el
movimiento de una misma part́ıcula material cuya posición final xf se des-
conoce, y se trata de deducir por medio de la posición inicial xi. El plantea-
miento Lagrangiano es adecuado al estudio de la mecánica de sólidos, en el
que es necesario incluir alguna ecuación constitutiva del comportamiento de
las part́ıculas del material como se muestra en la Figura 24.
El cambio en los desplazamientos se calculan de la siguiente forma
ut+∆t = ut + ∆u (3.1)
39
3.2. El Gradiente de Deformación
Este tensor contiene información acerca del alargamiento y rotaciones que
se generan en las fibras del material. Se define a través de
X t0 =

∂tx1
∂0x1
∂tx1
∂0x2
∂tx1
∂0x3
∂tx2
∂0x1
∂tx2
∂0x2
∂tx2
∂0x3
∂tx3
∂0x1
∂tx3
∂0x2
∂tx3
∂0x3
 (3.2)
Por lo tanto el tensor gradiente de deformación se define como el gradiente de
las coordenadas deformadas respecto a las iniciales. Usando notación indical
(3.65) se escribe
(X t0)ij =
∂txi
∂0xj
= Xi,j (3.3)
Cabe destacar que es un tensor definido positivo, representado por el opera-
dor gradiente ∇0.
X t0 = (∇0(xti)T )T (3.4)
∇0 = ∂∂0xi , (x
t
i)
T = [xt1, x
t
2, x
t
3]
Figura 25: Longitud y rotación de una fibra material en t, obtenidas a partir
de X t0 y posición t = 0
Los vectores d0x y dtx de la Figura 25 representan la orientación y lon-
gitud de las fibras del material conectadas a través de
dtx = X t0d
0x (3.5)
Concidérese el diferencial de volumen d0v del estado inicial formado por
un paraleṕıpedo cuyos lados están definidos por tres vectores diferenciales
40
Figura 26: Estado del cuerpo en el instante t = 0 y t = t , se muestran las
coordenadas deformadas xt1, x
t
2
Figura 27: Paraleṕıpedo de lados diferenciales
d0xi i = 1, 2, 3. Para estos volumenes concideraremos la conservación de la
masa
ρtdtv = ρ0d0v (3.6)
Entonces asumimos que
dtv = detX t0d
0v (3.7)
Por (3.7) la densidad en la configuración de referencia es
ρ0 = ρtdetX t0 (3.8)
d0x1 =
10
0
 ds1 (3.9)
41
d0x2 =
01
0
 ds2 (3.10)
d0x3 =
00
1
 ds3 (3.11)
d0v = ds1ds2ds3 (3.12)
3.3. Tensor Gradiente de Desplazamientos
Es aquel gradiente de los desplazamientos que toma como marco de refe-
rencia las coordenadas iniciales x0 del cuerpo de estudio
H0 = ∇0u =
∂u
∂x
(3.13)
En su representación matricial, los términos del tensor son
H0 = (∇0uT )T =
∂u1∂x1 ∂u1∂x2 ∂u1∂x3∂u2
∂x1
∂u2
∂x2
∂u2
∂x3
∂u3
∂x1
∂u3
∂x2
∂u3
∂x3
 (3.14)
(H0)ij =
∂ui
∂xj
(3.15)
Substituyendo las coordenadas xti en los elementos del tensor X
t
0 en función
de los desplazamientos ui se obtiene
X t0 = (∇0(xti)T )T = (∇0(xt0 + ui)T )T = I + (∇0uT )T (3.16)
X t0 = I +
∂u
∂x0
(3.17)
X t0 = I +H0 (3.18)
Se puede definir asimismo el tensor gradiente de los desplazamientos respecto
a las coordenadas deformadas xti en el instante t (demominado gradiente de
desplazamientos temporal)
Ht = ∇u =
∂u
∂xt
= (Ht)ij =
∂ui
∂xtj
(3.19)
42
3.4. Descomposición Polar del Tensor Gradiente de
Deformación
Como se ha mencionado el gradiente de deformación X t0 es definido posi-
tivo, luego se puede descompponer de la siguiente manera
X t0 = R
t
0U
t
0 (3.20)
Rt0 es una matriz ortogonal de rotación y U
t
0 es el alargamiento , U
t
0 es simétri-
ca. la Figura 28 nos indica que la deformación total se compone de una elon-
gación y de una rotación. Este tensor nos da herramientas suficientes para
Figura 28: Rotación y alargamiento componentes de la deformación total
definir el tensor derecho de Cauchy − Green, el cual establece una relación
entre la longitud inicial de una fibra y su estado deformado.
Ct0 = (X
t
0)
TX t0 (3.21)
Expresión que solo depende del alargamiento U t0
Ct0 = (U
t
0)
T (Rt0)
T )Rt0U
t
0 = (U
t
0)
2 (3.22)
3.5. Variación de los tensores X t0 y H0
Imaginemos un cuerpo deformado en el instante t, sometido a los des-
plazamientos u = xt − x0. Śı aplicásemos una variación virtual δu a dichos
desplazamientos , producirá a su vez variación de los gradientes de deforma-
ción. Este evento sencillament es
δ(H0)ij =
∂δui
∂x0j
(3.23)
43
Mientras que la variación del gradiente de deformación es directa
δ(X t0)ij = δ(H0)ij (3.24)
La ecuación (3.24) puede operar en función de las coordenadas deformadas
xti . Efectuando la regla de la cadena
δ(X t0)ij =
∂δxti
∂x0j
=
∂δui
∂xj
=
∂δui
∂xtj
∂xtj
∂xj
(3.25)
El propósito de aplicar estas variaciones es para comprobar que la variación
del gradiente de deformación es
δ(X t0) = δHtX
t
0 (3.26)
3.6. Tensor infinitesimal de deformaciones unitarias
Sea el campo de deformaciones existente en el sólido en un tiempo cual-
quiera, el tensor de deformaciones infinitesimales se define como
εij =
1
2
(
∂ui
∂xtj
+
∂uj
∂xti
) (3.27)
El tensor representado por (3.27) es lineal e igual a la parte simétrica del
tensor gradiente de desplazamientos
ε =
1
2
(H0 +H
T
t ) (3.28)
En forma vectorial (3.27) toma la siguiente forma
ε =

ε11
ε22
ε33
2ε12
2ε13
2ε23
 (3.29)
Otra alternativa es presentar el tensor de deformaciones unitarias en forma
compacta como
ε = AHt (3.30)
44
A es una matriz constante. En dos dimensiones (3.30) queda expresada
ε =
[ ∂u1
∂x1
∂u1
∂x2
+ ∂u2
∂x2
]
=
1 0 0 00 0 0 1
0 1 1 0


∂u1
∂x1
∂u1
∂x2
∂u2
∂x1
∂u2
∂x2
 (3.31)
Aplicando un operador variacional al tensor infinitesimal de deformaciones
obtenemos
δε =
1
2
(δHt + δH
T
t ) (3.32)
Cuyos términos equivalentes son
δεij0 =
1
2
(
∂δui
∂xtj
+
∂δuj
∂xti
) (3.33)
3.7. Tensor de Deformaciones Unitarias de Green −
Lagrange
El tensor de Green − Lagrange mide las