DocsTec-1086

Vista previa del material en texto

INSTITUTO TECONOLÓGICO Y DE ESTUDIOS SUPERIORES DE 
MONTERREY 
 
CAMPUS MONTERREY 
DIVISIÓN DE INGENIERÍA Y ARQUITECTURA 
PROGRAMA DE GRADUADOS DE INGENIERÍA 
 
 
 
ANÁLISIS DINÁMICO EXPLÍCITO DE CASCARONES 
 
 
TESIS 
 
 
PRESENTADA COMO REQUISITO PARCIAL PARA OBTENER EL GRADO 
ACADÉMICO DE: 
 
MAESTRO EN CIENCIAS 
ESPECIALIDAD EN INGENIERÍA Y ADMINISTRACIÓN DE LA 
CONSTRUCCIÓN CON ACENTUACIÓN EN INGENIERÍA ESTRUCTURAL 
 
 
POR: 
 
 
JESÚS FERNANDO GRACIA MENDÍVIL 
 
 
MONTERREY, N.L. MAYO DE 2003 
 
 i
INSTITUTO TECONOLÓGICO Y DE ESTUDIOS SUPERIORES DE MONTERREY 
 
CAMPUS MONTERREY 
DIVISIÓN DE INGENIERÍA Y ARQUITECTURA 
PROGRAMA DE GRADUADOS DE INGENIERÍA 
 
Los miembros del comité de tesis recomendamos que el presente proyecto de tesis presentado por el Ing. 
Jesús Fernando Gracia Mendívil sea aceptado como requisito parcial para obtener el grado académico de 
Maestro en Ciencias con especialidad en 
 
INGENIERÍA Y ADMINISTRACIÓN DE LA CONSTRUCCIÓN CON 
ACENTUACIÓN EN INGENIERÍA ESTRUCTURAL 
 
 Comité de Tesis 
 
 
 
 
 
 __________________________ 
 Sergio Gallegos Cázares, Ph.D. 
 Asesor 
 
 
 
 
 
 __________________________ __________________________ 
 Carlos Nungaray Pérez, M.Sc. Alex Elias Zúñiga,Ph.D. 
 Sinodal Sinodal 
 
 
 Aprobado: 
 
 
 
 
 
 __________________________ 
 Federico Viramontes Brown, Ph.D. 
 Director del Programa de Graduados en Ingeniería 
Mayo 2003 
 
 
 ii
ANÁLISIS DINÁMICO EXPLÍCITO DE CASCARONES 
 
Jesús Fernando Gracia Mendívil 
 
(RESUMEN) 
 
El objetivo de este trabajo de investigación es desarrollar una herramienta computacional capaz de 
modelar el comportamiento dinámico de estructuras compuestas de cascarones. Se trabaja en el código 
del programa PAEF para dicho propósito. 
El programa de análisis de elementos finitos PAEF hasta antes de la presente investigación era capaz de 
desarrollar análisis estático lineal para elementos tipo barra de dos nodos, elementos cuadrilaterales de 
cuatro nodos para esfuerzo plano y deformación plana y elementos bilineales de cuatro nodos para 
modelar el comportamiento de cascarones mediante el acoplamiento de placa y membrana. 
Se implementa las rutinas de análisis dinámico utilizando un método de integración explícito para 
resolver la ecuación de momentum lineal, enseguida, al elemento cascaron bilineal de cuatro nodos se le 
adicionan las funciones necesarias para su análisis dinámico como son la generación de la matriz de 
masa y su vector de fuerzas internas. 
Como resultado se obtiene un programa de elementos finitos capaz de generar análisis dinámico lineal 
de cascarones el cual se validó con soluciones conocidas como algunos ejemplos analíticos y otros 
problemas resueltos numéricamente de diferentes referencias. 
El análisis dinámico del programa PAEF puede complementarse con otros elementos finitos adicionando 
las rutinas de generación de matriz de masa y vector de fuerzas internas, además al programa PAEF se le 
puede adicionar comandos capaces de tomar en cuenta efectos de no-linealidad en materiales y 
geométricos. 
 
 
 iii
AGRADECIMIENTOS 
 
A Dios, por no dejarme solo y llevarme en tus brazos en los momentos más difíciles. 
 
A mis padres, por todo su apoyo incondicional, por quererse mucho durante todos estos años y 
querernos a todos sus hijos y sobre todo por creer en mi. Los quiero mucho y los llevo en mi corazón 
siempre. 
 
A mi esposa, por estar hombro a hombro conmigo, por todo el amor que me ha demostrado, por su 
paciencia y apoyo a pesar de todo este tiempo que la desatendí durante los estudios de posgrado y sobre 
todo en la realización de esta investigación. Te amo Pucharra. 
 
A mi asesor, por su enorme paciencia, sus extraordinarios consejos y por demostrar ser un gran amigo y 
consejero. Muchísimas gracias Doc. 
 
A mis sinodales, por brindarme su tiempo y sus conocimientos incondicionalmente los cuales 
complementaron la presente investigación. 
 
A Geométrica de México S.A. de C.V. que me alentaron a cursar el posgrado, por su confianza, consejos 
y apoyo, además de permitirme tomar parte de tiempo laboral y dedicarlo al estudio de mis clases. 
 
A mis hermanos, cuñados y amigos que con su ayuda directa o indirecta y con sus palabras de aliento 
me dieron ánimo antes y durante los estudios de posgrado. 
 
 iv
CONTENIDO 
 Lista de figuras ......................................................................................................................................vi 
 Lista de Tablas.....................................................................................................................................viii 
1.- Introducción .........................................................................................................................................1 
2.- Método de integración explícito..........................................................................................................3 
 2.1.- Descripción del método..................................................................................................................6 
 2.2.- Algoritmo .....................................................................................................................................11 
 2.3.- Estabilidad numérica ....................................................................................................................13 
 2.4.- Rutinas de elementos....................................................................................................................17 
 2.4.1.- Matriz de masa .....................................................................................................................19 
 2.4.2.- Vector de fuerzas internas ....................................................................................................22 
 2.4.3.- Matriz de amortiguamiento y vector de fuerzas de disipación.............................................24 
 2.4.4.- Fuerzas externas equivalentes por cargas aplicadas sobre el elemento................................26 
 2.5.- Comandos en PAEF .....................................................................................................................27 
3.- Elemento barra...................................................................................................................................32 
 3.1.- Introducción .................................................................................................................................32 
 3.2.- Matriz de masa del elemento barra ..............................................................................................33 
 3.3.- Vector de fuerzas internas del elemento barra .............................................................................36 
 3.4.- Validación ....................................................................................................................................37 
4.- Elemento cascarón .............................................................................................................................50 
 4.1.- Introducción .................................................................................................................................50 
 4.2.- Matriz de masa del elemento cascarón.........................................................................................53 
 4.3.- Vector de fuerzas internas del elemento cascarón .......................................................................58 
 
 v
 4.3.1.- Fuerzas internas de membrana .............................................................................................58 
 4.3.1.- Fuerzas internas de placa......................................................................................................62 
 4.4.- Validación ....................................................................................................................................655.- Conclusiones .......................................................................................................................................75 
6.- Contribuciones y trabajo futuro.......................................................................................................78 
7.- Referencias..........................................................................................................................................80 
 
 vi
Lista de Figuras 
Figura 2-1 Relación Esfuerzo – Deformación.............................................................................................5 
Figura 2.1-1 Incrementos de tiempo ...........................................................................................................6 
Figura 2.1-2 Desplazamientos discretos......................................................................................................7 
Figura 2.1-3 Velocidad media discreta .......................................................................................................8 
Figura 2.1-4 Velocidad discreta ..................................................................................................................9 
Figura 2.4-1 Algoritmo de integración numérica y rutinas de elementos .................................................18 
Figura 2.4.1-1 Funciones de interpolación................................................................................................19 
Figura 2.4.1-2 Aceleración del elemento ..................................................................................................20 
Figura 2.5-1 Ejemplo de función de tiempo, carga tipo pulso triangular..................................................30 
Figura 3.1-1 Elemento barra......................................................................................................................32 
Figura 3.2-1 Elemento lineal de dos nodos isoparamétrico ......................................................................33 
Figura 3.2-2 Masa uniformemente distribuida y concentrada del elemento barra....................................34 
Figura 3.4-1 Movimiento armónico simple...............................................................................................37 
Figura 3.4-2 Solución exacta y numérica para el ejemplo de movimiento armónico simple ...................38 
Figura 3.4-3 Comandos de PAEF para el ejemplo de movimiento armónico simple...............................39 
Figura 3.4-4 Masa sujeta a una función de carga con resorte y amortiguamiento....................................40 
Figura 3.4-5 Solución analítica y numérica con oscilación libre del Ejemplo 3.4.2.................................42 
Figura 3.4-6 Solución analítica y numérica con amortiguamiento del Ejemplo 3.4.2..............................43 
Figura 3.4-7 Comandos de PAEF del Ejemplo 3.4.2 ................................................................................44 
Figura 3.4-8 Función de Carga del Ejemplo 3.4.2 ....................................................................................44 
Figura 3.4-9 Barra en cantiliver con masa uniforme sujeta a un desplazamiento axial inicial en su 
extremo libre .............................................................................................................................................46 
 
 vii
Figura 3.4-10 Solución analítica y numérica del Punto a del Ejemplo 3.4.3 ............................................47 
Figura 3.4-11 Solución analítica y numérica del Punto b del Ejemplo 3.4.3............................................48 
Figura 3.4-12 Acercamiento de Figura 3.4-11 en el intervalo de tiempo 0.0004 < t < 0.0008.................48 
Figura 3.4-13 Comandos de PAEF del Ejemplo 3.4.3 ..............................................................................49 
Figura 4.1-1 Elemento cascarón de 4 nodos .............................................................................................50 
Figura 4.1-2 Desplazamientos de membrana ligados a rotaciones fuera del plano ..................................51 
Figura 4.1-3 Cinco puntos de muestreo ....................................................................................................51 
Figura 4.1-4 Elemento cascarón transformado .........................................................................................52 
Figura 4.2-1 Elemento cuadrilátero de 4 nodos isoparamétrico................................................................53 
Figura 4.2-2 Forma de las funciones de interpolación bilineal y cuadrática.............................................54 
Figura 4.2-3 Puntos de muestreo para regla de integración de 2x2 ..........................................................56 
Figura 4.3.1-1 Regla de integración con 5 puntos de muestreo ................................................................60 
Figura 4.3.2-1 Fuerzas de placa en el sentido local del elemento.............................................................64 
Figura 4.4-1 Viga doblemente empotrada sujeta a carga uniformemente distribuida ..............................65 
Figura 4.4-2 Resultados numéricos del Ejemplo 4.4.1 .............................................................................67 
Figura 4.4-3 Resultados numéricos del Ejemplo 4.4.2 .............................................................................69 
Figura 4.4-4 Domo bajo carga uniforme...................................................................................................70 
Figura 4.4-5 Modelo de elementos finitos del Ejemplo 4.4.3 ...................................................................71 
Figura 4.4-6 Dezplazamiento al centro del domo del Ejemplo 4.4.3........................................................72 
Figura 4.4-7 Viga en cantiliver con carga puntual en el extremo libre.....................................................72 
Figura 4.4-8 Resultados numéricos del Ejemplo 4.4.4 .............................................................................73 
Figura 4.4-9 Resultados del Ejemplo 4.4.4 en la Referencia [7] ..............................................................74 
 
 viii
Lista de Tablas 
Tabla 4.2-1 Valor de puntos de muestreo y factores de peso....................................................................56 
Tabla 4.3.1-1 Posición de los puntos de integración y valor de factor de peso para regla de 5 puntos ...61 
Tabla 4.3.2-1 Posición de los puntos de muestreo y su factor de peso para regla de integración 3x3 ....63 
Introducción Capítulo 1 
 1
1.- Introducción 
Toda aquella estructura laminar o elemento cuyo espesor se considera delgado con respecto a sus otras 
dimensiones, ya sea radio o longitud, se le conoce como cascarón, a este tipo de estructuras pertenecen 
los recipientes de pared delgada sujetos a presión, techos curvos, cúpulas, domos, fuselajes de aviones, 
tubos estructurales, etcétera. 
El estudio del comportamiento de los cascarones siempre se ha considerado un reto por la complejidad 
que presenta determinar de manera general todas sus propiedades. Se han hecho varias aproximaciones y 
simplificaciones de casos simples de cascarones para poder describir su comportamiento, tal es el caso 
de tubos de pared delgada sujetos a presión interna, domos esféricos, arcos, en fin, diversos tipos de 
estructuras o elementos estructurales los cuales, por distribuirse simétrica o uniformemente sus variables 
de análisis, pueden ser simplificadas y reducidas, sin embargo describir analíticamente un cascarón en 
forma general puede tornarse complicado e impráctico para representar cualquier tipo de cuerpo en la 
realidad. 
Con la técnica de los elementos finitos se ha logrado avanzar a grandes pasos en la forma que se puede 
determinar cuantitativamente el comportamiento de los cascarones. Cuando una estructura o elemento 
pertenece a este grupo se puede segmentar en pequeños dominios fácilmente analizablesaplicando 
herramientas de cómputo. 
Cuando un sistema de fuerzas externas se aplica sobre una estructura de manera inconstante o actúa 
cíclicamente puede lograr excitarla produciendo movimientos oscilatorios, si la frecuencia con que 
incurre en la estructura dicho conjunto de fuerzas externas es una tercera parte o menor de la frecuencia 
natural oscilatoria más baja de ésta, el problema puede resolverse con un análisis estático. Una gran 
parte de las estructuras de tipo cascarón pueden ser fácilmente excitadas por efectos accidentales como 
Introducción Capítulo 1 
 2
viento o sismo sobre domos o cubiertas de arco, también por cargas cíclicas como la vibración de un 
motor sobre un chasís. 
El análisis dinámico de estructuras nos permite estudiar el comportamiento de los cuerpos cuando es 
importante tomar en cuenta el tiempo como variable fundamental en el estudio. 
El presente trabajo tiene como propósito presentar el desarrollo de los planteamientos utilizados para 
implementar un análisis dinámico en un programa computacional de elementos finitos (PAEF), y en 
especial, el desarrollo de las rutinas propias de un cascarón necesarias para el análisis dinámico. 
Antes de implementar este tipo de análisis, PAEF incluía rutinas para desarrollo de análisis estático con 
relaciones constitutivas elásticas para elementos tipo barra, elementos tipo membrana para esfuerzo 
plano y deformación plana y elementos tipo cascarón, ahora es posible analizar dinámicamente 
estructuras compuestas de elementos tipo barra y cascarones 
 
Método de Integración Explícito Capítulo 2 
 3
2.- Método de integración explícito 
El método de solución para problemas dinámicos por integración explícita pertenece a los métodos 
directos de integración. Es un método directo porque la ecuación de movimiento, también conocida 
como Conservación de Momentum Lineal, permanece sin transformación a otra expresión equivalente y 
es integrada en forma directa mediante un procedimiento numérico. 
El método de integración explícito, también conocido como método de diferencia central, es el más 
popular en mecánica computacional y física, por la sencillez que presenta el algoritmo y la facilidad para 
implementarlo con la ayuda de herramientas de cómputo ya que se evita resolver ecuaciones 
simultáneas. 
La forma de integrar la ecuación de Conservación de Momentum Lineal es dividir el dominio de tiempo 
del problema en pequeños intervalos de tiempo y establecer las condiciones iniciales en cada intervalo 
para determinar las condiciones finales del mismo (ver sección 2.1) 
La Conservación de Momentum Lineal o Cantidad de Movimiento se describe mediante la siguiente 
ecuación: 
 )()()()( tFtFtVta ext=++ intCM (2.1) 
El primer término de la ecuación corresponde a fuerzas inerciales, donde: 
M: Matriz de masa 
a(t): Vector de aceleración, también descrito por: 
 
2
2 )(
)(
t
td
ta
∂
∂= (2.2) 
El segundo término de la ecuación se le conoce como fuerzas de amortiguamiento o fuerzas de 
disipación, donde: 
C: Matriz de amortiguamiento 
Método de Integración Explícito Capítulo 2 
 4
V(t): Vector de velocidad, determinado por la relación: 
 
t
td
tV
∂
∂ )(
)( = (2.3) 
El tercer término corresponde a las fuerzas internas, producto de las deformaciones que a su vez son 
funciones de los desplazamientos, esto es: 
 )()(int dftF = , siendo d = d(t) (2.4) 
Para (2.1), (2.2), (2.3) y (2.4): 
d(t): Vector de desplazamientos 
El término Fext(t) corresponde al vector de fuerzas externas las cuales pueden variar con respecto al 
tiempo. 
El método de integración directa se basa en dos ideas principales: la primera, se busca satisfacer la 
ecuación (2.1) para un intervalo de tiempo discreto ∆t en lugar de resolverse para cualquier tiempo t, 
básicamente significa un equilibrio estático el cual incluye los efectos de fuerzas inerciales y de 
disipación; la segunda idea en la cual se basa el método de integración directa es que se asume una 
variación de desplazamientos, velocidades y aceleraciones dentro de cada intervalo de tiempo ∆t, la 
forma de tales variaciones dentro del intervalo es lo que determina la precisión, estabilidad y el tiempo 
del procedimiento de solución. 
Observando la ecuación (2.1) se puede notar que si se mantienen las fuerzas externas y los 
desplazamientos constantes con respecto al tiempo, los dos primeros términos de la ecuación 
desaparecen, quedando únicamente: 
 extFdF =)(int (2.1b) 
la cual representa la solución para un problema estático, en el caso lineal, las fuerzas internas son el 
producto de la matriz de rigidez por los desplazamientos, esto es: 
Método de Integración Explícito Capítulo 2 
 5
 ddF K=)(int (2.1c) 
donde la matriz de rigidez se determina mediante las relaciones constitutivas elásticas entre esfuerzo y 
deformación. Sin embargo, se trata de dejar abierta la posibilidad de incluir efectos no-lineales en las 
propiedades del material a analizar para la implementación del algoritmo, es por eso que se optó por 
cambiar la representación del término Kd a una forma más general fint. 
 
Figura 2-1 Relación Esfuerzo – Deformación 
La Figura 2-1 muestra la comparativa entre esfuerzos calculados de manera lineal contra los calculados 
de manera no-lineal para el mismo grado de deformación, esta variación de esfuerzos podrá ser tomada 
en cuenta para calcular el estado de fuerzas internas como se describe más a detalle en las secciones 3.3 
y 4.3 
Método de Integración Explícito Capítulo 2 
 6
2.1.- Descripción del método 
El método de integración explícito consiste en dividir el dominio del problema en pequeños intervalos 
de tiempo ∆t desde n=0 hasta NPT (número de pasos de tiempo) y establecer las condiciones finales de 
cada intervalo en función de las condiciones iniciales del mismo. También se le conoce como método de 
diferencia central, este nombre se debe a la aproximación de diferencias finitas usadas para expresar la 
aceleración el tiempo tn. 
Así mismo, una estructura se segmenta espacialmente con elementos finitos interconectados por nudos, 
donde a cada nudo se le puede dar hasta 6 componentes de desplazamiento (3 traslacionales y 3 
rotacionales), de tal manera que los desplazamientos, velocidades y aceleraciones son calculadas para 
cada una de estas componentes las cuales se denominan grados de libertad. 
A continuación se describen los parámetros usados en el método de integración explícito y la relación 
que existen entre ellos: 
 
Figura 2.1-1 Incrementos de tiempo 
 
n: Contador, desde 0 hasta NPT (número de pasos de tiempo) 
n
medt∆ : Incremento de tiempo del paso n. 
nt∆ : Diferencia de tiempo entre tiempo medio del paso n+1 y tiempo medio del paso n 
 1−−=∆ nnnmed ttt (2.1.1) 
Método de Integración Explícito Capítulo 2 
 7
tn: Tiempo del intervalo n. 
 nmed
nn ttt ∆+= −1 (2.1.1b) 
tmed
n : Tiempo medio del intervalo n, definido por: 
 
2
1−−=
nn
n
med
tt
t (2.1.2) 
Despejando la aceleración de la ecuación (2.1) obtenemos: 
 )(1 ndamp
n
net
n ffa −= −M (2.1.3) 
donde: 
an: Aceleración del intervalo n. 
M: Matriz de masa. 
n
netf : Fuerza neta del intervalo n, definida por: 
 nnext
n
net fff int−= (2.1.4) 
n
extf : Fuerza externa del intervalo n. 
nf int : Fuerza interna del intervalo n. 
n
dampf : Fuerza de disipación del intervalo n. 
 
Figura 2.1-2 Desplazamientos discretos 
Método de Integración Explícito Capítulo 2 
 8
En la Figura 2.1-2 se puede observar como se suponen los desplazamientos discretos en cada intervalo 
de tiempo, por diferencias finitas podemos obtener la siguiente relación: 
 
n
med
nn
n
med t
dd
V
∆
−=
−1
 (2.1.5) 
donde: 
n
medV : Velocidad media del intervalo n. 
nd : Desplazamiento del intervalo n. 
Expresando los desplazamientos en función de la velocidad media definida por la ecuación (2.1.5) se 
obtiene: 
 nmedmed
nn Vtdd ∆+= −1 (2.1.5b) 
 
Figura 2.1-3 Velocidad mediadiscreta 
Una vez determinada la velocidad media a partir de los desplazamientos discretos, en la Figura 2.1-3 se 
aprecia cómo se determina la aceleración por diferencias finitas mediante: 
 
n
n
med
n
medn
t
VV
a
∆
−
=
+1
 (2.1.6) 
 
Método de Integración Explícito Capítulo 2 
 9
 
Figura 2.1-4 Velocidad discreta 
La ecuación (2.1.6) supone una aceleración constante para el intervalo n que va desde nmedt hasta 
1+n
medt , lo 
que hace posible determinar la velocidad para un intervalo n interpolando las velocidades medias, o lo 
que es lo mismo, es posible también definir la aceleración en función de la velocidad discreta del 
intervalo mediante: 
 
1
1
1
−
−
−
−
−
=
nn
med
nn
medn
tt
VV
a (2.1.7) 
 
n
med
n
n
med
n
n
tt
VV
a
−
−
= (2.1.8) 
donde: 
an: Aceleración del intervalo n. 
Con las ecuaciones (2.1.7) y (2.1.8) se puede expresar la velocidad media en función de los valores del 
intervalo inmediato anterior y la velocidad al final del intervalo en función de los valores del mismo: 
 ( )111 −−− −+= nnmednnnmed ttaVV (2.1.7b) 
 ( )nmednnnmedn ttaVV −+= (2.1.8b) 
Método de Integración Explícito Capítulo 2 
 10
Como puede observarse la función de los desplazamientos, velocidades y aceleraciones aún permanece 
como incógnita, pero es posible determinar valores discretos para cualquier intervalo de tiempo tn. 
Entre las ventajas del procedimiento podemos enumerar: 
1. Las ecuaciones resultantes son fácilmente programables para la implementación en un equipo de 
cómputo. 
2. Si se calcula la matriz de masa concentrada (ver secciones 3.2 y 4.2), entonces la matriz de masa M 
sólo tendrá elementos en su diagonal. 
3. Las fuerzas internas y las fuerzas de disipación pueden calcularse a nivel de cada elemento finito 
para cada grado de libertad en sus nudos, ensamblando solamente los vectores resultantes de fuerzas. 
4. El procedimiento puede ser utilizado contemplando linealidad o no-linealidad en los materiales. 
También el procedimiento tiene sus desventajas: 
1. Un parámetro importante para la estabilidad del procedimiento y para la convergencia de los 
resultados es escoger un incremento de tiempo ∆tmed adecuado y suficientemente pequeño para que 
una onda pueda viajar a través de todos los elementos finitos, en la sección 2.3 se explicará más a 
detalle este parámetro el cual resulta de intervalos de tiempo del orden de 0.00001s<∆t<0.0001s, lo 
cual para poder simular algún fenómeno de la naturaleza que dure 5 segundos puede llegar a 
consumir hasta 500,000 iteraciones. 
2. Haciendo referencia a la desventaja anterior, si se trata de avanzar con pasos de tiempo más largo el 
algoritmo se vuelve inestable y la solución diverge, resultando deformaciones completamente 
desproporcionadas. 
Método de Integración Explícito Capítulo 2 
 11
2.2.- Algoritmo 
En la sección anterior se describieron las ecuaciones necesarias para poder llevar a cabo una integración 
numérica explícita, ahora se detallarán los pasos para la implementación. 
Es necesario hacer hincapié que el procedimiento necesita algunas rutinas propias de los elementos 
finitos usados para la segmentación de la estructura, tales como matriz de masa, matriz de 
amortiguamiento, vector de fuerzas internas y vector de fuerzas de disipación, dichas rutinas son 
independientes de la rutina general de solución pero necesarias para el procedimiento, el cual se detalla a 
continuación: 
1. Para un paso inicial se determina condiciones iniciales 0t , 0medV , 
0d . 
2. Se calculan las matrices de masa de los elementos y se ensambla en una matriz de masa global M, se 
aprovecha que para sólidos la masa es constante. 
3. Se determina el vector de fuerzas netas a nivel elemento, mediante la ecuación (2.1.4) 
nn
ext
n
net fff int−= para el paso 
0t y se ensambla en un vector de fuerzas netas global. En esta misma 
rutina se determina medt∆ crítico. Las fuerzas netas deben de calcularse a nivel elemento para 
cuando se detallen cargas al elemento finito así como a nivel global cuando se apliquen directamente 
sobre los nudos de la estructura. 
4. Se calculan las fuerzas de disipación dampf a nivel elemento para después ensamblarse en un vector 
de fuerzas de disipación global, mediante 
 nmed
en
damp Vf C= (2.2.1) 
donde: 
eC : Matriz de amortiguamiento del elemento e. 
5. Se determina la aceleración mediante la ecuación (2.1.3) ( )ndampnnetn ffa −= −1M 
Método de Integración Explícito Capítulo 2 
 12
6. Se actualiza las variables de tiempo con la ecuación (2.1.1c) nmed
nn ttt ∆+= −1 , y la ecuación (2.1.2) 
( ) 21−−= nnnmed ttt 
7. Se resuelve para Velocidad media con la ecuación (2.1.7b) ( )111 −−− −+= nnmednnnmed ttaVV 
8. Se obliga a la velocidad recién calculada a cumplir las condiciones de frontera marcadas por el 
problema, esto es para apoyos fijos o móviles. 
9. Se resuelve para desplazamientos mediante la ecuación (2.1.5b) nmedmed
nn Vtdd ∆+= −1 
10. Se determina las fuerzas netas a nivel elemento, mediante la ecuación (2.1.4) nnext
n
net fff int−= para el 
tiempo nt , se ensambla en un vector de fuerzas netas global y se actualiza el incremento de tiempo 
crítico nmedt∆ de ser necesario. Tomar en cuenta la misma observación que en el paso 3 para las 
fuerzas externas 
11. De ser necesario se actualiza a nivel elemento la matriz de amortiguamiento C, se calcula el vector 
de disipación mediante la ecuación (2.2.1) nmed
en
damp Vf C= y se procede a ensamblar en un vector de 
fuerzas de disipación global. 
12. Se actualiza la aceleración con la ecuación (2.1.3) ( )ndampnnetn ffa −= −1M 
13. Se resuelve la velocidad al tiempo tn con la ecuación (2.1.8b) ( )nmednnnmedn ttaVV −+= 
14. Se hace el balance energético de comprobación 
15. Se incrementa el paso n ← n +1 
16. Continuar en el paso 6 hasta que se haya completado la simulación 
Método de Integración Explícito Capítulo 2 
 13
2.3.- Estabilidad numérica 
Un problema resuelto numéricamente es estable cuando pequeños cambios en los datos iniciales dan 
como resultado pequeños cambios en la solución numérica. 
Un importante parámetro para garantizar la estabilidad de la solución en cualquier tiempo tn del 
procedimiento es la determinación del incremento de tiempo crítico ∆tcrit. Se trata de que el incremento 
de tiempo sea tal que permita avanzar en la simulación con el menor costo de tiempo de cómputo y a la 
vez que la solución sea estable numéricamente. 
El tiempo crítico del sistema será aquel que sea el más crítico de todos los elementos, de tal manera que 
se garantice que una onda viaje por el elemento excitándolo en el más alto modo de vibración, como lo 
marca la Referencia [1], esto es: 
 ( )ee
e
e
critt ξξω
−+≤∆ 12min 2 (2.3.1) 
para todos los elementos e de la segmentación, donde: 
ωe : Frecuencia del más alto modo de vibración del elemento e. 
ξe : Porcentaje de amortiguamiento relativo a la frecuencia ωe del elemento e, expresado 
en decimales (0.02 = 2%) por la relación: 
 
22
1 e
e
o
e
aa ω
ω
ξ += (2.3.2) 
a0, a1: Coeficientes de amortiguamiento relativos a la masa y rigidez del elemento respectivamente. 
Para obtener la frecuencia natural del elemento es necesario hacer un análisis modal y encontrar los 
modos característicos de vibración mediante la relación: 
 02 =− eint MK ω
e (2.3.3) 
Para: 
Método de Integración Explícito Capítulo 2 
 14
e
intK : Rigidez interna del elemento. 
eM : Matriz de masa del elemento. 
Para rapidez en los cálculos, los valores característicos son usualmente obtenidos por fórmulas simples. 
Es de notar que la rigidez interna del elemento para el caso lineal elástico es la matriz de rigidez, sin 
embargo para el caso no-lineal dicha matriz de rigidez puede cambiar sus valores dependiendo de la 
magnitud de esfuerzos a los cuales esté sujeto, para este caso será necesario actualizar el intervalo de 
tiempo crítico para cada iteración 
Para el elemento tipobarra, definido en Capítulo 3, se toman las siguientes relaciones para obtener la 
más alta frecuencia de vibración: 
 
ρ
E
c = (2.3.4) 
 
ρ
σ xEc −= (2.3.4b) 
 
l
c2
max =ω (2.3.5) 
 
l
c32
max =ω (2.3.5b) 
 
 
donde: 
c: Velocidad de onda, en la ecuación (2.3.4) es para el caso lineal y la ecuación (2.3.4b) es 
 para el caso no-lineal. 
E: Módulo de Young del elemento 
σx: Esfuerzo axial 
Método de Integración Explícito Capítulo 2 
 15
ρ: Densidad del elemento 
l: Longitud del Elemento 
La ecuación (2.3.5) corresponde a la frecuencia máxima calculada con la matriz de masa concentrada, y 
la ecuación (2.3.5b) es calculada con la matriz de masa consistente, para la definición de matriz de masa 
concentrada y consistente ver sección 2.4, 3.2 y 4.2. Adicional a este cálculo, se le aplicó un factor de 
reducción de 0.80 como lo propone la referencia [1] 
Para el caso del elemento cascarón descrito en la sección 4, se tiene las siguientes relaciones tomadas de 
la referencia [8] 
 
)1(3412.0
1
υρ +
= Ec (2.3.6) 
 
t
c2
max =ω (2.3.7) 
Para: 
υ: Módulo de Poisson del elemento. 
t: Espesor promedio del elemento. 
Para el caso no-lineal, es necesario implementar un método numérico capaz de obtener los modos 
característicos de vibración de la ecuación (2.3.3) 02 =− eint MK ω
e aplicada a un cascarón con la 
rigidez interna modificada por los esfuerzos. 
Una manera de ir verificando la estabilidad del sistema en cada intervalo de tiempo tn, es generando un 
balance energético entre energía externa y energía interna: 
 )()(
2
1 111 −−− −−+= nnTnnnn ffddWW intintintint (2.3.8) 
 )()(
2
1 111 −−− −−+= next
n
ext
Tnnn
ext
n
ext ffddWW (2.3.9) 
Método de Integración Explícito Capítulo 2 
 16
 )()(
2
1 111 −−− −−+= ndamp
n
damp
Tnnn
damp
n
damp ffddWW (2.3.10) 
 )()(
2
1 nTnn
kin VMVW = (2.3.11) 
Donde nWint ,
n
extW ,
n
dampW son el trabajo de fuerzas internas, trabajo de fuerzas externas, trabajo de fuerzas 
de disipación y nkinW es la energía cinética para el tiempo t
n respectivamente. 
La conservación de energía requiere que: 
 ),,,max( intint dampextkindampextkin WWWWWWWW ε≤−−+ (2.3.12) 
Donde ε es una pequeña tolerancia, la referencia [1] la propone del orden de 10-2 
Método de Integración Explícito Capítulo 2 
 17
2.4.- Rutinas de elementos 
Como se analizó en la sección 2.2, cuando se segmenta espacialmente un dominio en elementos finitos, 
la rutina de integración numérica explícita requiere que los elementos alimenten de parámetros al 
algoritmo para su solución, dichos parámetros son responsabilidad única y exclusiva de cada elemento 
finito. En la Figura 2.4-1 se puede apreciar las rutinas propias del algoritmo y la de los elementos finitos 
y a que nivel del algoritmo se deben de calcular. 
Las rutinas de los elementos finitos se enumeran a continuación: 
1. Matriz de masa eM . 
2. Vector de Fuerzas internas ef int . 
3. Matriz de amortiguamiento eC y Vector de fuerzas de disipación ef damp . 
4. Fuerzas externas equivalentes por cargas aplicadas sobre el elemento. 
Método de Integración Explícito Capítulo 2 
 18
 
Figura 2.4-1 Algoritmo de integración numérica y rutinas de elementos 
Método de Integración Explícito Capítulo 2 
 19
2.4.1.- Matriz de masa 
La matriz de masa de un elemento finito es la relación que existe entre el vector de aceleraciones 
nodales y las fuerzas en el mismo sentido que la aceleración, se parte de la interpolación de valores 
nodales en cualquier punto del dominio del elemento 
 
 
Figura 2.4.1-1 Funciones de interpolación 
En la Figura 2.4.1-1 podemos apreciar como se interpolan valores en el dominio de un elemento tipo 
barra (ver Capítulo 3) mediante las funciones de interpolación de los valores nodales, donde: 
iN : Función de interpolación para el valor nodal i. 
iu : Valor nodal del desplazamiento en el nudo i. 
dando como resultado 
 !
=
=
n
i
ii uNu
1
 (2.4.1.1) 
Para u calculada en cualquier punto del dominio Ω del elemento. En formato matricial la ecuación 
(2.1.1.1) puede expresarse como: 
 " # Nu=
$
$
%
$
$
&
'
$
$
(
$
$
)
*
=
i
i
u
u
u
NNNu
...
... 1
1
21 (2.4.1.1b) 
Método de Integración Explícito Capítulo 2 
 20
 
Figura 2.4.1-2 Aceleración del elemento 
En la Figura 2.4.1-2 al elemento se le imprime un sistema de fuerzas nodales, el cual provoca un estado 
de aceleración en el dominio Ω del elemento, utilizando las mismas funciones de interpolación que se 
utilizaron para los desplazamientos podemos expresar la aceleración en cualquier punto del elemento 
mediante: 
 Na=a (2.4.1.2) 
de tal manera que la función de fuerza inercial se puede representar con la siguiente ecuación: 
 Nammaf == (2.4.1.3) 
Aplicando el principio de trabajo virtual: 
 Ω= +
Ω
dufuf T
T
δδ (2.4.1.4) 
sustituyendo en (2.4.1.4) las ecuaciones (2.4.1.3) y (2.4.1.1b) en el integrando obtenemos: 
 Ω= +
Ω
dumuf TT
T
δδ NNa (2.4.1.4) 
Sacando los valores nodales de la aceleración y el vector de desplazamientos virtuales del integrando 
obtenemos la definición de matriz de masa del elemento: 
 Ω= +
Ω
dmTe NNM (2.4.1.4) 
Método de Integración Explícito Capítulo 2 
 21
A esta matriz se le conoce como la matriz de masa consistente, ya que se formuló a partir de las mismas 
ecuaciones con la que se plantea la matriz de rigidez del elemento, quedando una matriz cuadrada, 
simétrica y de las mismas dimensiones que la matriz de rigidez. 
Es de notar que al ensamblar las matrices de masa consistente de todos los elementos obtenemos una 
matriz general de masa con términos diferentes de 0 por fuera de su diagonal principal, en la sección 2.2 
se menciona que para aumentar la rapidez de cómputo del procedimiento se propone el uso de una 
matriz de masa que contenga solamente términos en su diagonal, a dicha matriz se le conoce 
comúnmente como matriz de masa concentrada. En las secciones 3.2 y 4.2 se explicará más a detalle 
cómo obtener la matriz de masa concentrada del elemento a partir de la matriz de masa consistente. 
m11 m12 ... m1 j
m21 m22 ... m2 j
... ... ... ...
mi1 mi 2 ... mij
, 
" 
- 
- 
- 
- 
- 
. 
# 
/ 
/ 
/ 
/ 
/ 
→
mc11 0 0 0
0 mc22 ... 0
... ... ... ...
0 0 ... mcij
, 
" 
- 
- 
- 
- 
- 
. 
# 
/ 
/ 
/ 
/ 
/ 
 
Método de Integración Explícito Capítulo 2 
 22
2.4.2.- Vector de fuerzas internas 
Las fuerzas internas de un elemento finito fint es el conjunto de fuerzas nodales capaz de imprimir un 
estado de esfuerzos σ en su dominio. 
El estado de esfuerzo σ de un elemento finito puede ser calculado mediante relaciones constitutivas 
elásticas o mediante un análisis no-lineal del material, tal como se muestra en la Figura 2-1, de tal 
manera que la formulación presentada puede representar fuerzas internas producto de relaciones lineales 
de esfuerzo y deformación o también de relaciones no-lineales, de cualquier manera, la evaluación de 
dichos esfuerzos será rutina propia del elemento finito y en la presente sección se tomarán como datos 
de entrada para el cálculo de fuerzas internas. 
Definamos la deformación ε del elemento de la siguiente manera: 
 { }x
u
∂
∂ε = (2.4.2.1) 
tomando la definición de los desplazamientos de (2.4.1.1b) y sustituyendo en (2.4.2.1) obtenemos: 
 { } Buu
N ==
x∂
∂ε (2.4.2.1b) 
donde: 
 { }x∂
∂N
B = (2.4.2.2) 
 
Aplicando el principio del trabajo virtual de los esfuerzos sobre todo el dominio del elemento 
obtenemos: 
 Ω= +
Ω
dW Tσδεint (2.4.2.3) 
sustituyendo la definición de la deformación de (2.4.2.1b) obtenemos: 
Método de Integración Explícito Capítulo 2 
 23
 Ω= +
Ω
dW TT σδ Buint (2.4.2.3b) 
extrayendo del integrando el vector de desplazamiento virtual obtenemos la forma del vector de fuerzas 
internas: 
 Ω= +
Ω
df TσBint (2.4.2.4) 
Esta definición del Vector de fuerzas internas es la que se desarrollará para los elementos estudiados en 
el presente trabajo. La forma que tienen las funcionesde interpolación para desarrollar la matriz B varía 
para cada elemento y para cada tipo de esfuerzo que se estudie. 
Se puede notar que para el cálculo de las fuerzas internas los esfuerzos pueden ser calculados con 
relaciones constitutivas elásticas o mediante relaciones no-lineales. Supongamos que los esfuerzos 
provengan de un análisis lineal elástico, esto es: 
 εσ D= (2.4.2.5) 
donde D es la matriz de relaciones constitutivas del material, sustituyendo la definición de la 
deformación de la ecuación 2.4.1.1b se obtiene 
 uKDBB intint =00
1
2
33
4
5
Ω= +
Ω
udf T (2.4.2.4b) 
Donde el término entre paréntesis de la ecuación (2.4.2.4b) lo conocemos como la matriz de rigidez 
interna del elemento. 
Método de Integración Explícito Capítulo 2 
 24
2.4.3.- Matriz de amortiguamiento y vector de fuerzas de disipación 
La matriz de amortiguamiento C puede calcularse a nivel elemento sin necesidad de ensamblarla en una 
matriz de amortiguamiento general, lo único que tendrá que ensamblarse será el vector de fuerzas de 
disipación o fuerzas de amortiguamiento. 
Para el cálculo de la matriz de amortiguamiento de Rayleigh C como lo muestra la referencia [5] se 
toma las rutinas de elementos ya definidas: la matriz de masa M y la matriz de rigidez interna K para 
generar la matriz de amortiguamiento con la siguiente ecuación: 
 KMC 10 aa += (2.4.3.1) 
donde: 
a0, a1: Coeficientes de amortiguamiento proporcional a la masa y a la rigidez del elemento, 
respectivamente. 
El vector de fuerzas de disipación se obtiene con la ecuación (2.2.1) nmed
en
damp Vf C= 
Como puede observarse, al tratarse de un análisis no-lineal, la matriz de rigidez interna puede cambiar 
sus valores en cada iteración en función de los esfuerzos, así que será necesario recalcular dicha matriz 
en cada intervalo de tiempo de ser necesario. 
Los coeficientes a0, a1 se calculan mediante la relación (2.3.2) 
22
1 i
i
o
i
aa ω
ω
ξ += , proponiendo algún 
porcentaje de amortiguamiento relativo a 2 modos de vibración ωi global de la estructura, generalmente 
se opta por escoger los dos primeros modos de vibración. Sin embargo, para obtener dichos modos de 
vibración globales será necesario implementar un análisis modal, el cual implica solucionar la ecuación 
(2.3.3) 02 =− eint MK ω
e a nivel global para los primeros valores característicos, además si se está 
usando un análisis no-lineal la rigidez interna cambiaría en cada iteración generando excesiva carga de 
Método de Integración Explícito Capítulo 2 
 25
cómputo. Por estas razones los coeficientes a0, a1 deben de ser estimados como dato de entrada para la 
integración numérica explícita ya que el propósito del presente trabajo difiere de un análisis modal. 
Método de Integración Explícito Capítulo 2 
 26
2.4.4.- Fuerzas externas equivalentes por cargas aplicadas sobre el elemento 
En la formulación del algoritmo para integración explícito es necesario distinguir del ensamble de las 
fuerzas externas aplicadas sobre los nodos y las fuerzas externas aplicadas a los elementos finitos, ya 
que las primeras pueden ensamblarse directamente en el vector de fuerzas externas globales, sin 
embargo, para las cargas sobre los elementos finitos es necesario que éste evalúe la contribución nodal 
equivalente de dichas fuerzas. 
Supongamos una fuerza de cuerpo q y un sistema de tracciones τ aplicado sobre un elemento finito, 
aplicando nuevamente el principio de trabajo virtual tenemos: 
 dSudquW
S
TT ++ +Ω=
Ω
τδδδ (2.4.4.1) 
Sustituyendo (2.4.1.1b): 
 dSdqW
S
TTTT ++ +Ω=
Ω
τδδδ NuNu (2.4.4.1b) 
Extrayendo el vector de desplazamientos virtuales del integrando obtenemos la forma de evaluar las 
fuerzas externas nodales equivalentes de las cargas al elemento: 
 dSdqf
S
TT
ext ++ +Ω=
Ω
τNN (2.4.4.2) 
Método de Integración Explícito Capítulo 2 
 27
2.5.- Comandos en PAEF 
La rutina de integración explícita se desarrolló dentro de los módulos de análisis del programa PAEF 
(Programa de Análisis de Elementos Finitos), a continuación se detallan los pasos necesarios para correr 
un análisis dinámico explícito, los comandos implementados se muestran sombreados: 
COMANDO COMENTARIOS 
INICIA Se inicia la definición del problema 
NOMBRE EJEMPLO Se da nombre del problema 
DIMENSIONES 3 
NUMERO DE ELEMENTOS 3 
NUMERO DE NODOS 6 
MATERIALES 2 
NOD_ELEMENTO 4 Máximo número de nodos por elemento 
GDL NODAL 6 
FIN Se cierra proceso de datos iniciales 
MALLA Inicia entrada de datos de malla 
COORDENADAS Se asignan coordenadas nodales 
1 0.000 0.000 [NUDO] [CORD X] [CORD Y] 
2 1.000 0.000 
3 0.000 0.500 
4 1.000 0.500 
5 0.000 1.000 
6 1.000 1.000 
FIN Se cierra proceso de entrada de coordenadas 
INCIDENCIAS Se abre proceso para entrada de incidencias 
1 1 3 [ELEMENTO] [NODO 1] [NODO 2] …[NODO N] 
2 2 4 
3 3 4 6 5 
FIN Se cierra proceso de entrada de incidencias 
FUNCION PULSO FUNCION [NOMBRE] 
Se define funciones de tiempo 
ARCHIVO "C:\PULSO.TXT" ARCHIVO [PATH] 
FIN Fin de entrada de función, se puede definir tantas funciones 
como sea requerido, todas referenciadas a un archivo de 
texto 
PROPIEDADES Se abre proceso para entrada de propiedad de materiales 
NUMERO 1 TIPO 1, NUMERO [NUM PROPIEDAD] TIPO [TIPO ELEM] 
PROPIEDADES E 2E+11, PROPIEDADES E [MOD. YOUNG], 
A 1.5E-4, A [AREA], 
RO 7850.0, RO [DENSIDAD] 
MASA CONCENTRADA, MASA [CONCENTRADA, CONSISTENTE] 
SALIDA NODOS SALIDA [NODOS, GAUSS] 
En los tipos de elementos 1 corresponde a elemento tipo 
barra y 6 corresponde a elemento tipo cascarón, para ver 
definición de los tipos de elementos ver capítulo 3 y 4 
Método de Integración Explícito Capítulo 2 
 28
NUMERO 2 TIPO 6, NUMERO [NUM PROPIEDAD] TIPO [TIPO ELEM] 
PROPIEDADES EMOD 2E+11, PROPIEDADES EMOD [MOD. YOUNG] 
NU 0.3 ESPESOR 0.20, NU [COEF POISSON] ESPESOR [ESPESOR] 
CORTE 0.8333 CORTE [COEF. DE CORTE] 
RO 800.0 RO [DENSIDAD] 
MASA CONCENTRADA, MASA [CONCENTRADA, CONSISTENTE] 
SALIDA NODOS SALIDA [NODOS, GAUSS] 
FIN Se cierra proceso de entrada de materiales 
 
RESTRICCIONES APOYO RESTRICCION [NOMBRE] 
PUNTUALES Restricción de tipo puntual 
1 U=0.0 [NODO] [GDL] [VALOR DE RESTRICCION] 
1 V=0.0 
1 W=0.0 
1 ROTX=0.0 
1 ROTY=0.0 
1 ROTZ=0.0 
2 U=0.0 
2 V=0.0 
2 W=0.0 
2 ROTX=0.0 
2 ROTY=0.0 
2 ROTZ=0.0 
FIN Se cierra este tipo de restricción 
FIN Se cierra entrada de restricciones, pueden meterse tantas 
restricciones como sean requeridas 
CONDICION INICIAL DEFORMADA CONDICION INICIAL [NOMBRE] 
DESPLAZAMIENTOS DESPLAZAMIENTOS 
3 U=0.01 [NODO] [GDL] [VALOR DE DESPLAZAMIENTO] 
4 U=0.01 
FIN Termina condición inicial de desplazamientos 
FIN Termina condiciones iniciales, pueden alimentarse 
cuantas condiciones iniciales sean requeridas 
CARGA MUERTA CARGA [NOMBRE] 
NODAL 
4 U=100.0 [NODO] [GDL] [VALOR DE FUERZA NODAL] 
5 V=MASA 1000.0 [NODO] [GDL] MASA [VALOR MASA CONCENTRADA] 
Se le asigna a algún nodo un valor de masa concentrada 
en determinado grado de libertad 
6 W=FUNCION PULSO, [NODO] [GDL] FUNCION [NOMBRE], 
AMPLIFICACION 1.0, AMPLIFICACION [NUM], 
DESFASAMIENTO 0.0, DESFASAMIENTO [NUIM] 
Se le asigna al nodo una carga de tiempo de acuerdo a la 
función, la cual estará amplificada y desfasada si se requiere 
FIN Termina definición de carga nodales 
FIN Termina tipo de carga, se puede alimentar con tantas 
cargas como sea necesaria 
FIN Se cierra proceso de entrada de malla 
Método de Integración Explícito Capítulo 2 
 29
 
COMANDOS COMANDOS 
CONTROL Entrada de datos para análisis dinámico explícito 
ANALISIS DINAMICO ANALISIS [ESTÁTICO, DINÁMICO] 
Se especifica que tipo de análisis se generará 
METOD EXPLICITO METODO [EXPLICITO, IMPLICITO] 
Se especifica que se hará integración numérica explícita 
TOLERANCIA 0.01 TOLERANCIA [NUM] 
Tolerancia para el balance energético 
A0 0.02 A0 [NUM] 
A1 0.002 A1 [NUM] 
Se especifica los coeficientes de amortiguamiento 
relativos a la masa y a la rigidez respectivamente 
GRABAR CADA 10 ITERACIONES GRABARCADA [NUM] [SEGUNDOS, ITERACIONES] 
Se especifica que el reporte de datos se hará una vez que 
concluyan cierto numero de iteraciones o que se cumpla a 
cada múltiplo del tiempo alimentado 
IGNORAR ESFUERZOS [GRABAR, IGNORAR] [ESFUERZOS, DEFORMACIONES] 
IGNORAR DEFORMACIONES Se especifica si se desea reportar esfuerzos o deformaciones 
en la salida de los datos 
FIN Fin de datos de control 
 
SELECCIONA RESTRICCION APOYO De los juegos de restricciones, cargas o condiciones 
SELECCIONA CARGA MUERTA iniciales se indica cual se selecciona para hacer el análisis 
SELECCIONA CONDICION INICIAL DEFORMADA al menos una restricción y un tipo de carga deben de ser 
seleccionados 
MASA TODOS MASA [LISTA DE ELEMENTOS] 
Manda a calcular la matriz de masa de los elementos 
de la lista 
RIGIDEZ TODOS Manda a calcular la matriz de rigidez de los elementos 
de la lista 
ENSAMBLA Rutina que se encarga de ensamblar ya sea la matriz de 
rigidez o de masa según el tipo de análisis a efectuar 
INICIALIZA PROCESO Rutina para establecer las condiciones de inicio del 
proceso dinámico, ya sea condición inicial o los datos de 
los últimos cálculos hechos, también declara en memoria 
RAM los vectores de fuerzas, aceleraciones, velocidades y 
desplazamientos 
ITERA 5000 ITERA [NUM] 
Se especifica cuantas iteraciones se hará en el proceso 
INCREMENTA PASO Los últimos datos calculados los pasa al vector del paso 
tn-1, también actualiza las variables de tiempo, 
prepara los vectores del paso tn para ser calculados 
y determina si este paso será exportado 
VELOCIDAD MEDIA Calcula la velocidad media del intervalo 
DESPLAZAMIENTOS Cálculo de desplazamientos 
FUERZAS Genera el vector de fuerzas netas, pide a cada elemento 
que calcule su fuerza de externa de cuerpo y su fuerza interna 
para ensamblarlos en el vector de fuerzas netas 
Método de Integración Explícito Capítulo 2 
 30
AMORTIGUAMIENTO Pide a cada elemento su matriz de masa y rigidez (modificada 
de ser necesario) para el cálculo de fuerzas de disipación 
y ensambla el vector de fuerzas de amortiguamiento global 
ACELERACIONES Calcula las aceleraciones nodales 
VELOCIDAD Misma rutina que VELOCIDAD MEDIA, solo que ahora 
calcula la velocidad al final del intervalo 
BALANCE Comprueba el balance energético 
TERMINA Tolas los comando comprendidos entre ITERA y TERMINA 
se ejecutarán tantas veces como lo marque el comando 
ITERA [NUM] 
CIERRA PROCESO Los vectores declarados en INICIALIZA PROCESO son 
liberados de memoria 
EXPORTA GID Salida del archivo de datos para el post-procesamiento 
FIN Fin de comandos 
FIN Salida del programa 
En los comandos recién descritos hay que hacer hincapié en ciertos datos que no aparecen aquí: 
1. En la declaración de funciones se hace referencia a un archivo de texto guardado en disco, dicho 
archivo tiene los datos discretos de una función de tiempo con tiempo inicial y tiempo final, para 
cualquier tiempo no incluido dentro de los límites de la función la salida será 0. 
 
 
Figura 2.5-1 Ejemplo de función de tiempo, carga tipo pulso triangular 
La Figura 2.5-1 muestra una función de carga con respecto al tiempo, la forma que debe de ser el 
archivo para describir la ecuación es la siguiente: 
 
 
Método de Integración Explícito Capítulo 2 
 31
Archivo de texto Comentario 
3 [NUM] Número de datos que contiene el archivo 
0.0 0.0 [TIEMPO] [VALOR DE FUNCION] 
0.1 5000.0 
0.2 0.0 
 
Donde para cualquier valor de tiempo incluido en el archivo el valor de la función se interpola, y 
para cualquier valor de tiempo no incluido en la misma la función será 0. 
2. En un mismo tipo de carga también se pueden incluir diferentes cargas al elemento finito 
dependiendo del elemento, aquí únicamente se detallaron cargas y masas nodales y funciones de 
tiempo para cargas aplicadas en los nodos. 
 
Elemento Barra Capítulo 3 
 32
3.- Elemento barra 
3.1.- Introducción 
Aunque el propósito principal del trabajo es el desarrollo de un análisis dinámico y el enfoque a 
cascarones, se incluye la formulación del elemento tipo barra, ya que la simplicidad que presenta 
determinar su matriz de masa y su vector de fuerzas internas permite centrar la atención en la 
verificación del algoritmo de integración explícito. 
El elemento tipo barra es un elemento finito unidimensional definido por 2 nudos, solamente cuenta con 
rigidez en su eje axial principal, presenta una sección prismática constante de área A en toda su longitud 
L definida entre sus 2 nudos. 
 
Figura 3.1-1 Elemento barra 
En la Figura 3.1-1 se muestra el elemento tipo barra que puede estar dispuesto en 1,2 y 3 dimensiones 
aunque sus propiedades sean definidas unidimensionales. Para cada nudo le puede corresponder 1,2 o 3 
grados de libertad según sea analizado en 1,2 o 3 dimensiones. La misma Figura 3.1-1 muestra la 
numeración de los grados de libertad. 
Elemento Barra Capítulo 3 
 33
3.2.- Matriz de masa del elemento barra 
La ecuación (2.4.1.4) que define la matriz de masa consistente de un elemento puede ser expresada de la 
siguiente manera para el elemento tipo barra: 
 dxNNAM
L
T
e += ρ (3.2.1) 
donde: 
ρ: Densidad del material 
A: Área de la sección transversal de la barra 
Mapeando el elemento barra a un elemento isoparamétrico unidimensional mostrado en la figura 3.2-1 
 
Figura 3.2-1 Elemento lineal de dos nodos isoparamétrico 
la matriz de funciones de interpolación para los desplazamientos nodales es la siguiente: 
 N = N1 N2" # (3.2.2) 
 /
#
.
-
"
,
=
21
21
00
00
NN
NN
N (3.2.2b) 
 
/
/
/
#
.
-
-
-
"
,
=
21
21
21
0000
0000
0000
NN
NN
NN
N (3.2.2c) 
Las ecuaciones (3.2.2), (3.2.2b) y (3.2.2c) corresponden para el caso de 1,2 y 3 dimensiones 
respectivamente, donde: 
 
2
)1( ξξ i
iN
+
= (3.2.3) 
Cambiando el diferencial de longitud al espacio natural isoparamétrico tenemos: 
Elemento Barra Capítulo 3 
 34
 ξ
ξ
d
L
dx
LXX
d
dx
222
12 =∴=
−
= (3.2.4) 
Integrando la ecuación (3.2.1) obtenemos: 
 /
#
.
-
"
,
21
12
6
ALρ
 (3.2.5) 
 
/
/
/
/
#
.
-
-
-
-
"
,
2010
0201
1020
0102
6
ALρ
 (3.2.5b) 
 
/
/
/
/
/
/
/
/
#
.
-
-
-
-
-
-
-
-
"
,
200100
020010
002001
100200
010020
001002
6
ALρ
 (3.2.5c) 
donde (3.2.5), (3.2.5b) y (3.2.5c) corresponde a la matriz de masa consistente del elemento barra para el 
caso de 1,2 y 3 dimensiones respectivamente 
Para obtener la matriz de masa concentrada a partir de la matriz de masa consistente se pueden sumar 
todos los términos de un renglón y asignando este valor a la diagonal principal, con este planteamiento 
se supone que la masa total del elemento barra se concentra en sus nodos como lo muestra la Figura 3.2-
2 
 
Figura 3.2-2 Masa uniformemente distribuida y concentrada del elemento barra 
 
Elemento Barra Capítulo 3 
 35
Con el anterior razonamiento obtenemos: 
 /
#
.
-
"
,
10
01
2
ALρ
 (3.2.6) 
 
/
/
/
/
#
.
-
-
-
-
"
,
1000
0100
0010
0001
2
ALρ
 (3.2.6b) 
 
/
/
/
/
/
/
/
/
#
.
-
-
-
-
-
-
-
-
"
,
100000
010000
001000
000100
000010
000001
2
ALρ
 (3.2.6c) 
donde (3.2.6), (3.2.6b) y (3.2.6c) corresponde a la matriz de masa concentrada del elemento barra para el 
caso de 1,2 y 3 dimensiones respectivamente. 
Es bueno notar que tanto para el caso de la matriz de masa consistente y la matriz de masa concentrada 
si se suma la contribución de un reglón de la matriz que corresponde para un mismo grado de libertad en 
ambos nodos obtenemos la masa total del elemento, esto se debe a que la masa permanece constante 
para cada dirección de la aceleración, en este caso se está contemplando 1,2 o 3 direcciones de 
aceleración ortogonales dependiendo si es un problema de 1,2 o 3 dimensiones respectivamente, por lo 
tanto, la masa total que engloba la matriz de masa del elemento, ya sea concentrada o consistente es de 
1m, 2m o 3m para 1,2 o 3 dimensiones respectivamente. 
Es necesario recalcar que estas matricesfueron calculadas para el sistema coordenado local, por lo tanto, 
hay que aplicar una transformación al sistema coordenado global para ensamblarla a la matriz de masa 
global. 
Elemento Barra Capítulo 3 
 36
3.3.- Vector de fuerzas internas del elemento barra 
La ecuación (2.4.2.4) usada para obtener el vector de fuerzas internas puede ser expresada de la 
siguiente forma para el caso de la barra: 
 dxBAf
L
T+= σint (3.3.1) 
ya que el esfuerzo axial de una barra se considera constante en toda su longitud. A diferencia de la 
matriz de masa, solo se analiza para el caso unidimensional, ya que solo se considera esfuerzo en el 
sentido axial de la barra. 
Tomando la definición de la matriz B de la ecuación (2.4.2.2), tenemos la siguiente expresión: 
 /#
/
-"
-=
x
N
x
N
B
∂
∂
∂
∂ 21 (3.3.2) 
sin embargo, la funciones de interpolación Ni, son funciones de la variable ξ, así que hay que usar la 
regla de la cadena para obtener las derivadas con respecto al sistema local: 
 
x
N
x
N
∂
∂ξ
∂ξ
∂
∂
∂ = (3.3.3) 
donde: 
 /#
/
-"
- +−=/
#
/
-
"
-
=
2
1
2
121
∂ξ
∂
∂ξ
∂
∂ξ
∂ NNN
 (3.3.4) 
utilizando la misma relación que en la ecuación (3.2.4) el vector de fuerzas internas se obtiene: 
 
%
&
'
(
)
*
+
−
=
%
&
'
(
)
*
+
−
= + 1
12
1
1
2int
Adx
L
A
f
L
σσ (3.3.5) 
y al igual que la matriz de masa, este vector hay que expresarlo en el sistema global antes de 
ensamblarse. 
Elemento Barra Capítulo 3 
 37
3.4.- Validación 
Los ejemplos que se ilustran sirven para demostrar el algoritmo de integración explícita ayudándose de 
la facilidad de implementación de las rutinas propias del elemento finito tipo barra, que como se pudo 
ver en las secciones 3.2 y 3.3 la matriz de masa y el vector de fuerzas internas pueden obtenerse 
analíticamente. 
Ejemplo 3.4.1 
Se toma la solución de un movimiento armónico simple de un grado de libertad, donde una masa es 
sujeta a un desplazamiento inicial d y se pone a oscilar por la rigidez que le imprime un resorte de 
rigidez K como se muestra en la figura 3.4-1 
 
Figura 3.4-1 Movimiento armónico simple 
La rigidez del resorte se representa con las propiedades de una barra tal que K=AE/L, a la misma barra 
se le restringe uno de sus extremos y en el otro se le asigna un valor a la masa m y un desplazamiento d. 
Datos de entrada. 
A : 1.50E-04 m2 
E : 2.00E+11 Pa 
L : 0.50 m 
d : 0.01 m 
 
K : 6.00E+07 N/m 
m : 1000.0 kg 
 
La frecuencia natural y el período del sistema se define por: 
 
m
K=ω (3.4.1) 
Elemento Barra Capítulo 3 
 38
 
π
ω
2
=T (3.4.2) 
ω : 244.9 rad/seg 
T : 2.57E-02 cps 
 
Los resultados del problema se muestran en la Figura 3.4-2: 
 
Figura 3.4-2 Solución exacta y numérica para el ejemplo de movimiento armónico simple 
La Figura 3.4-2 muestra los resultados analíticos y numéricos del problema con un procedimiento de 
integración numérico explícito, con 5 elemento finito, 6500 iteraciones grabado a cada 50 iteraciones, en 
la figura 3.4-3 se muestra el archivo de comandos para el problema 
Elemento Barra Capítulo 3 
 39
 
Figura 3.4-3 Comandos de PAEF para el ejemplo de movimiento armónico simple 
 
INICIA
! Datos de definicion del problema
!
Inicia
Nombre MAS
Dimensiones 2
Numero de elementos 5
Numero de nodos 6
Materiales 1
NOD_Elemento 2
GDL Nodal 2
fin
!
Malla
Incidencias
1 1 2
2 2 3
3 3 4
4 4 5
5 5 6
fin
!
Coordenadas
1 0.00000 0.00000
2 0.10000 0.00000
3 0.20000 0.00000
4 0.30000 0.00000
5 0.40000 0.00000
6 0.50000 0.00000
fin
!
Propiedades
!
! Propiedades para la barra elástica de dos
nodos
Numero 1 Tipo 1 Propiedades E
2E+11 A 1.5E-4 Ro 7850.0 Masa Concentrada Salida
Gauss
!
!
Elemento 1 Propiedad 1
Elemento 2 Propiedad 1
Elemento 3 Propiedad 1
Elemento 4 Propiedad 1
Elemento 5 Propiedad 1
fin
!
! Inicia definicion de los juegos de restricciones
!
Restricciones Apoyo
Puntuales
1 U = 0
1 V = 0
1 V = 0
2 V = 0
3 V = 0
4 V = 0
5 V = 0
6 V = 0
fin
fin
 
!
! Definición de condiciones iniciales
!
Condicion Inicial DespUnit
Desplazamientos
2 U = 0.002
3 U = 0.004
4 U = 0.006
5 U = 0.008
6 U = 0.01
fin
fin
!
! Definición de juegos de cargas
!
Carga MasaConc
!
!
! Masas concentradas
!
Nodal
6 U = Masa 1000
fin
!
! Termina definicion de cargas
fin
! Termina comandos de malla
fin
! --------------------------------------------------------
---
Comandos
! Pasos previos para empezar analisis dinamico
Control
Analisis Dinamico
Metodo Explicito
Tolerancia 0.01
A0 0.0
A1 0.0
Grabar Cada 50 Iteraciones
Ignorar Esfuerzos
Ignorar Deformaciones
fin
Selecciona Restricciones Apoyo
Selecciona Carga MasaConc
Selecciona Condicion Inicial DespUnit
Masa Todos
Rigidez Todos
Ensambla
Inicializa Proceso
ITERA 6501
incrementa paso
velocidad media
desplazamientos
Fuerzas
Amortiguamiento
aceleraciones
velocidad
balance
TERMINA
Cierra Proceso
Exporta GID
fin
fin
Elemento Barra Capítulo 3 
 40
 
Ejemplo 3.4.2 
Se analiza el siguiente ejemplo, donde se somete a una masa concentrada sujeta a un resorte de rigidez K 
amortiguado por C a oscilar mediante una carga variable con respecto al tiempo, como se muestra en la 
Figura 3.3-4 
 
Figura 3.4-4 Masa sujeta a una función de carga con resorte y amortiguamiento 
En este ejemplo se determina la masa concentrada como la mitad de la masa de una barra de longitud L, 
área A y densidad ρ, con una matriz de masa concentrada. Se analizan los resultados numéricos con la 
solución analítica del problema 
Datos del problema: 
L : 500 
E : 200000 
A : 10 
ρ : 0.001 
 
ao : 0.02 
a1 : 0.002 
 
m : 2.5 
c : 8.05 
K : 4000 
La ecuación gobernante del sistema considerando amortiguamiento es la siguiente: 
 ( )tPd
t
d
t
d =++ 400005.85.2
2
2
∂
∂
∂
∂
 (3.4.3) 
Elemento Barra Capítulo 3 
 41
y la solución analítica se compone de 3 ecuaciones: 
Para 0 ≤ t < 0.1: 
( )
320
05.84000
)9676.39(0252.0)9676.39(3138.0 61.161.1
−++−= −− ttCosetSened tt (3.4.3b) 
Para 0.1 ≤ t < 0.2: 
( )( )
320
05.82.04000
)9676.39(6167.0)9676.39(7497.0 61.161.1
−−++−= −− ttCosetSened tt (3.4.3c) 
Y para 0.2 ≤ t : 
)9676.39(0378.1)9676.39(6555.0 61.161.1 tCosetSened tt −− +−= (3.4.3d) 
si al sistema se considera con oscilación libre, entonces la solución analítica es: 
Para 0 ≤ t < 0.1: 
 t
tSen
d 5.12
2.3
)40( +−= (3.4.3e) 
Para 0.1 ≤ t < 0.2: 
 )2.0(5.12)40(7210.0)40(4730.0 −+−= ttSentCosd (3.4.3f) 
Y para 0.2 ≤ t : 
 )40(6756.0)40(7822.0 tSentCosd −= (3.4.3g) 
 
 
Elemento Barra Capítulo 3 
 42
 
Figura 3.4-5 Solución analítica y numérica con oscilación libre del Ejemplo 3.4.2 
 
Elemento Barra Capítulo 3 
 43
 
Figura 3.4-6 Solución analítica y numérica con amortiguamiento del Ejemplo 3.4.2 
En las Figura 3.4-5 y 3.4-6 se muestra la comparativa de la solución analítica contra los resultados 
numéricos. 
Elemento Barra Capítulo 3 
 44
 
Figura 3.4-7 Comandos de PAEF del Ejemplo 3.4.2 
 
Figura 3.4-8 Función de Carga del Ejemplo 3.4.2 
INICIA
! Datos de definicion del problema
!
Inicia
Nombre BARRA
Dimensiones 2
Numero de elementos 1
Numero de nodos 2
Materiales 1
NOD_Elemento 2
GDL Nodal 2
fin
! ---------------------------------------------------
------------
Malla
Incidencias
1 1 2
2 2 3
3 3 4
4 4 5
5 5 6
fin
!
Coordenadas
1 0.00000 0.00000
2 500.00000 0.00000
fin
!
! Definicion de Funciones
Funcion Fun1
Archivo "c:\disco_local\P.txt"
fin
Propiedades
!
! Propiedades para la barra elástica de dos
! nodos
Numero 1 Tipo 1 Propiedades E
200000 A 10.0 Ro 0.001 Masa Concentrada Salida Gauss
!
!
Elemento 1 Propiedad 1
fin
!
! Inicia definicion de los juegos de restricciones
!
Restricciones Apoyo
Puntuales
1 U = 0
1 V = 0
2 V = 0
fin
fin 
!
! Definición de juegos de cargas
!
Carga P
!
!
! Funciones nodales
!
Nodal
2 U = funcion Fun1 Amplificacion 1.0
Desfazamiento 0.0
fin
!
! Termina definicion de cargas
fin
! Termina comandos de malla
fin
! --------------------------------------------------------
---
Comandos
! Pasosprevios para empezar analisis dinamico
Control
Analisis Dinamico
Metodo Explicito
Tolerancia 0.01
A0 0.02
A1 0.002
Grabar Cada 0.005 Segundos
Ignorar Esfuerzos
Ignorar Deformaciones
fin
Selecciona Restricciones Apoyo
Selecciona Carga P
Masa Todos
Rigidez Todos
Ensambla
Inicializa Proceso
ITERA 301
incrementa paso
velocidad media
desplazamientos
Fuerzas
Amortiguamiento
aceleraciones
velocidad
balance
TERMINA
Cierra Proceso
Exporta GID
fin
fin 
3
0.0 0.0
0.1 5000.0
0.2 0.0 
Elemento Barra Capítulo 3 
 45
La Figura 3.4-7 muestra el archivo de comandos para el ejemplo 3.4.2 y la Figura 3.4-8 es el archivo de 
datos de las cargas, se generaron 300 iteraciones a cada 0.005 segundos. 
Elemento Barra Capítulo 3 
 46
Ejemplo 3.4.3 
Se ilustra el problema 8.5 de la Referencia [12], se trata también de una barra empotrada en un extremo 
y libre en el otro con masa uniformemente distribuida en toda su longitud a la cual se le imprime un 
desplazamiento axial inicial δ0 en su extremo libre, la simulación comienza en el momento que se libera 
la fuerza necesaria para imprimir dicho desplazamiento y la barra comienza a oscilar. El desplazamiento 
δ en cualquier punto de la barra a cualquier tiempo a partir del inicio se describe mediante la ecuación 
(3.4.4) 
 
0
0
0
0
0
1
2
3
3
3
3
3
4
5
+
0
1
2
3
4
5 +
+
−= !
∞
= L
t
E
n
Cos
L
xn
Sen
nn
n
2
)12(
2
)12(
)12(
)1(8
0
22
0 ρ
π
π
π
δδ (3.4.4) 
Datos del problema 
L : 1.50 m 
E : 2.00E+11 Pa 
A : 1.50E-04 m2 
ρ : 7850.0 kg/m3 
δ0 : 0.001 m 
 
Para la solución numérica se modeló con 10 elementos tipo barra, con 500 iteraciones de 0.00002 
segundos cada una 
 
 
Figura 3.4-9 Barra en cantiliver con masa uniforme sujeta a un desplazamiento axial inicial en su 
extremo libre 
Elemento Barra Capítulo 3 
 47
La Figura 3.4-9 muestra los puntos en los cuales se reportan resultados analíticos y numéricos 
 
Figura 3.4-10 Solución analítica y numérica del Punto a del Ejemplo 3.4.3 
 
 
 
 
 
Elemento Barra Capítulo 3 
 48
 
Figura 3.4-11 Solución analítica y numérica del Punto b del Ejemplo 3.4.3 
En las Figura 3.4-10 y 3.4-11 se muestran 4 tipos diferentes de resultados correspondientes a los puntos 
a y b y a la solución numérica y analítica de cada punto. La Figura 3.4-12 muestra un acercamiento de la 
Figura 3.4-11 en la zona de 0.0004 < t < 0.0008 para apreciar la aproximación de la solución numérica. 
 
Figura 3.4-12 Acercamiento de Figura 3.4-11 en el intervalo de tiempo 0.0004 < t < 0.0008
Elemento Barra Capítulo 3 
 49
 
Figura 3.4-13 Comandos de PAEF del Ejemplo 3.4.3
INICIA
! Datos de definicion del problema
!
Inicia
Nombre Barra2
Dimensiones 2
Numero de elementos 10
Numero de nodos 11
Materiales 1
NOD_Elemento 2
GDL Nodal 2
fin
! ----------------------------------------------------
-----------
Malla
Incidencias
1 1 2
2 2 3
3 3 4
4 4 5
5 5 6
6 6 7
7 7 8
8 8 9
9 9 10
10 10 11
fin
!
Coordenadas
1 0.00000 0.00000
2 0.15000 0.00000
3 0.30000 0.00000
4 0.45000 0.00000
5 0.60000 0.00000
6 0.75000 0.00000
7 0.90000 0.00000
8 1.05000 0.00000
9 1.20000 0.00000
10 1.35000 0.00000
11 1.50000 0.00000
fin
!
Propiedades
!
! Propiedades para la barra elástica de dos
nodos
Numero 1 Tipo 1 Propiedades E
2E+11 A 0.00015 Ro 7850 Masa Concentrada Salida Gauss
!
!
Elemento 1 Propiedad 1
Elemento 2 Propiedad 1
Elemento 3 Propiedad 1
Elemento 4 Propiedad 1
Elemento 5 Propiedad 1
Elemento 6 Propiedad 1
Elemento 7 Propiedad 1
Elemento 8 Propiedad 1
Elemento 9 Propiedad 1
Elemento 10 Propiedad 1
fin
!
! Inicia definicion de los juegos de restricciones
!
Restricciones Apoyo
Puntuales
1 U = 0
1 V = 0
1 V = 0
2 V = 0
3 V = 0
4 V = 0
5 V = 0
6 V = 0
6 V = 0
7 V = 0
8 V = 0
9 V = 0
10 V = 0
11 V = 0
fin
fin 
!
! Definición de condiciones iniciales
!
Condicion Inicial dto
Desplazamientos
2 U = 0.0001
3 U = 0.0002
4 U = 0.0003
5 U = 0.0004
6 U = 0.0005
7 U = 0.0006
8 U = 0.0007
9 U = 0.0008
10 U = 0.0009
11 U = 0.001
fin
fin
!
! Definición de juegos de cargas
!
Carga Fuerza
!
!
! Termina definicion de cargas
fin
! Termina comandos de malla
fin
! --------------------------------------------------------
---
Comandos
! Pasos previos para empezar analisis dinamico
Control
Analisis Dinamico
Metodo Explicito
Tolerancia 0.01
A0 0.0
A1 0.0
Grabar Cada 0.000020 Segundos
Ignorar Esfuerzos
Ignorar Deformaciones
fin
Selecciona Restricciones Apoyo
Selecciona Carga Fuerza
Selecciona Condicion Inicial dto
Masa Todos
Rigidez Todos
Ensambla
Inicializa Proceso
ITERA 501
incrementa paso
velocidad media
desplazamientos
Fuerzas
Amortiguamiento
aceleraciones
velocidad
balance
TERMINA
Cierra Proceso
Exporta GID
fin
fin 
Elemento Cascarón Capítulo 4 
 50
4.- Elemento cascarón 
4.1.- Introducción 
En este capítulo se presenta la formulación para incluir las rutinas de cálculo de parámetros dinámicos 
para elementos tipo cascarón. 
Un elemento cascarón se compone de la superposición de un elemento membrana con esfuerzo plano y 
un elemento placa. Se adopta un elemento bilineal de cuatro nodos desarrollado en la referencia [6] para 
3 dimensiones con 6 grados de libertad por nudo, 3 traslacionales y 3 rotacionales como se muestra en la 
Figura 4.1-1. 
 
Figura 4.1-1 Elemento cascarón de 4 nodos 
 
La modelación de cascarones con elementos finitos ha sido un problema muy interesante por resolver 
por las diversas dificultades numéricas que presenta su planteamiento, el elemento desarrollado en la 
referencia [6] contiene formulaciones innovadoras que solucionan las dificultades relacionadas con su 
planteamiento: 
Elemento Cascarón Capítulo 4 
 51
1. La formulación de membrana para esfuerzo plano incluye desplazamientos nodales asociados con 
giros fuera del plano, esto permite eliminar un modo de cero energía al incluirle una rigidez 
rotacional a ese grado de libertad. 
 
Figura 4.1-2 Desplazamientos de membrana ligados a rotaciones fuera del plano 
2. Una regla de integración de 5 puntos de muestreo en lugar de una regla gaussiana de integración 
completa mediante 3x 3 puntos de integración ayuda a eliminar el bloqueo de membrana que se 
produce porque la rigidez de membrana es más grande que la rigidez de flexión. 
 
Figura 4.1-3 Cinco puntos de muestreo 
3. La combinación de las funciones de interpolación adecuadas y los criterios de integración aumentan 
la precisión, estabilidad y convergencia de los resultados y la sensibilidad a la distorsión geométrica 
de los elementos. 
Elemento Cascarón Capítulo 4 
 52
4. Se implementa una transformación para el elemento a un plano equivalente requerido en la 
formulación de membrana y placa en lo que se incluye rotaciones de los nudos, esto es debido a que 
4 nodos en el espacio forman un paraboloide hiperbólico y difícilmente están alineados formando un 
plano. Se adopta esta transformación para la salida de la matriz de masa y el vector de fuerzas 
internas. 
 
Figura 4.1-4 Elemento cascarón transformado 
El elemento se implementó en el programa de análisis por elementos finitos PAEF obteniendo buenos 
resultados en las pruebas con análisis estático, entre las rutinas con que cuenta podemos enumerar: 
1. Generación de matriz de rigidez elástica 
2. Salida de esfuerzos y deformaciones en puntos de gauss o nodos. 
3. Especificación del vector de referencia para salida de resultados. 
4. Vector de fuerzas externas equivalente para cargas aplicadas sobre el cuerpo 
En este trabajo se detalla el procedimiento utilizado para implementar las rutinas añadidas al elemento 
para su análisis dinámico, tales como matriz de masa y vector de fuerzas internas. 
Elemento Cascarón Capítulo 4 
 53
4.2.- Matriz de masa del elemento cascarón 
La formulación de la matriz de masa de la ecuación (2.4.1.4) Ω= +
Ω
dmNN TeM se desarrolla en este 
capítulo para el elemento cascarón. 
Una vez que el elemento se ha transformado a su plano de referencia se mapea a un elemento 
isoparamétrico de1,1 +≤≤− ηξ mostrado en la Figura 4.4-2 
 
Figura 4.2-1 Elemento cuadrilátero de 4 nodos isoparamétrico 
Se toma de la referencia [6] las funciones de interpolación para los desplazamientos traslacionales en 
función de las traslaciones y rotaciones de los nudos 
Elemento Cascarón Capítulo 4 
 54
 
Figura 4.2-2 Forma de las funciones de interpolación bilineal y cuadrática 
 [ ] [ ] [ ] [ ][ ]4321 NNNNN = (4.2.1) 
para: 
 Ni =
Nui 0 0 0 0 Nuφz i
0 Nvi 0 0 0 Nvφzi
0 0 Nwi Nwϕxi Nwϕy i 0
, 
" 
- 
- 
- 
- 
. 
# 
/ 
/ 
/ 
/ 
 (4.2.2) 
donde los desplazamientos en función de los desplazamientos nodales se interpolan bilinealmente 
mediante: 
 )1)(1(
4
1
,, ηηξξ iiwiviui NNN ++= (4.2.3) 
y los desplazamientos en función de las rotaciones nodales se interpolan cuadráticamente mediante 
 )()( 21
12
41
41
1 YYNYYNN uuu z −+−= φφφ (4.2.4) 
 )()( 32
23
12
12
2 YYNYYNN uuu z −+−= φφφ (4.2.5) 
 )()( 43
34
23
23
3 YYNYYNN uuu z −+−= φφφ (4.2.6) 
 )()( 14
41
34
34
4 YYNYYNN uuu z −+−= φφφ (4.2.7) 
 )()( 12
12
14
41
1 XXNXXNN uuv z −+−= φφφ (4.2.8) 
Elemento Cascarón Capítulo 4 
 55
 )()( 23
23
21
12
2 XXNXXNN uuv z −+−= φφφ (4.2.9) 
 )()( 34
34
32
23
3 XXNXXNN uuv z −+−= φφφ (4.2.10) 
 )()( 41
41
34
34
4 XXNXXNN uuv z −+−= φφφ (4.2.11) 
En la Figura 4.2-2 se muestra la forma de las funciones de interpolación cuadráticas siguientes: 
 )1)(1(
8
1 212 ηξφ −−=uN (4.2.12) 
 )1)(1(
8
1 234 ηξφ +−=uN (4.2.13) 
 )1)(1(
8
1 223 ξηφ +−=uN (4.2.14) 
 )1)(1(
8
1 241 ξηφ −−=uN (4.2.15) 
los grados de libertad nodal se enumeran en el siguiente orden: ui,vi,wi,φxi,φyi,φzi para i=1 a 4 (número de 
nodos) 
La integral para la matriz de masa consistente ahora se expresa en coordenadas naturales: 
 ηξρ
η ξ
ddNtjN Te + +=M (4.2.16) 
Donde: 
ρ: Densidad del material 
t: Espesor del cascarón 
j: Determinante del jacobiano de la transformación, el cual se define por: 
 
∂η
∂
∂η
∂
∂ξ
∂
∂ξ
∂
yx
yx
j = (4.2.17) 
Elemento Cascarón Capítulo 4 
 56
 
$
$
$
$
$
%
$
$
$
$
$
&
'
$
$
$
$
$
(
$
$
$
$
$
)
*
/
#
.
-
"
,
=
%
&
'
(
)
*
4
4
3
3
2
2
1
1
4321
4321
0000
0000
Y
X
Y
X
Y
X
Y
X
NNNN
NNNN
y
x
uuuu
uuuu (4.2.18) 
 
Para la integración de la matriz de masa se hicieron pruebas con diferentes grados de integración 
numérica y se optó por usar integración numérica mediante una regla de gauss de integración completa 
de 2 x 2 puntos de muestreo mostrados en la Figura 4.2-3 y en la tabla 4.2-1: 
 
Figura 4.2-3 Puntos de muestreo para regla de integración de 2x2 
i ξξξξ ηηηη Wξξξξ Wηηηη 
1 -0.57735 02691 89626 -0.57735 02691 89626 1.00000 00000 00000 1.00000 00000 00000
2 +0.57735 02691 89626 1.00000 00000 00000 1.00000 00000 00000
3 +0.57735 02691 89626 -0.57735 02691 89626 1.00000 00000 00000 1.00000 00000 00000
4 +0.57735 02691 89626 1.00000 00000 00000 1.00000 00000 00000
 
Tabla 4.2-1 Valor de puntos de muestreo y factores de peso 
Elemento Cascarón Capítulo 4 
 57
Si se usa el mismo procedimiento de la barra para obtener la matriz de masa concentrada, en el elemento 
cascarón la masa total de los grados de libertad asociados a las traslaciones es tres veces la masa total del 
elemento y produce ceros en la diagonal en los grados de libertad asociado a las rotaciones, esto es 
lógico ya que en el primer caso la masa del elemento se debe conservar para cada dirección considerada 
en la aceleración y en este caso se están considerando 3 direcciones ortogonales, y en el segundo caso 
una masa concentrada en un nodo no produciría ninguna fuerza inercial rotacional, además de que las 
rotaciones en un análisis dinámico no son tan importantes como lo son desplazamientos, sin embargo, la 
solución del sistema produciría singularidad en dichos grados de libertad provocando inestabilidad y 
resultados aleatorios de desplazamientos. El siguiente procedimiento propuesto en la referencia [11] 
genera una matriz de masa diagonal concentrada más sofisticada y aplica para aquellos elementos cuyos 
grados de libertad traslacionales son mutuamente paralelos, tales como vigas y cascarones. 
1. A partir de la matriz de masa consistente sumar en la diagonal principal todos aquellos grados de 
libertad asociados con los desplazamientos: 
 21,20,19,15,14,13,9,8,7,3,2,1, ==! iMm
i
ii (4.2.19) 
2. CalcularMte , la matriz total del elemento, utilizando integración numérica. 
 ηξρ
η ξ
ddtjt e + +=M (4.2.20) 
3. Escalar los coeficientes de la diagonal multiplicándolos por la razón 3Mt/m, y convirtiendo en ceros 
todo elemento que no sea parte de la diagonal principal, de tal manera que se conserva la masa total 
del elemento en cada sistema ortogonal de traslación y se conserva proporcionalmente la masa 
rotacional. 
Elemento Cascarón Capítulo 4 
 58
4.3.-Vector de fuerzas internas del elemento cascarón 
Para obtener el vector de fuerzas internas total debe de ser calculado para ambas partes que 
superimpuestas conforman un elemento cascarón: Membrana y Placa 
4.3.1.-Fuerzas internas de membrana 
La ecuación (2.4.2.4) para la parte de membrana se desarrolló a partir de la formulación del elemento 
cascarón de la referencia [6] y se llegó a la siguiente expresión: 
 [ ] [ ] dABdA
N
N
N
Bf o
T
AA
xy
y
x
T ττ++ +
$
%
$
&
'
$
(
$
)
*
=int (4.3.1.1) 
donde: 
 [ ] [ ] [ ]" #φuu NNSB ][= (4.3.1.2) 
 [ ] { } [ ] [ ]" # [ ]" #φφτ NNNRB uuT −+= 0 (4.3.1.3) 
 [ ] [ ] /
#
.
-
"
,
=
4321
4321
0000
0000
,
uuuu
uuuu
u NNNN
NNNN
NN φ (4.3.1.4) 
 Nui =
1
4
(1 +ξ iξ )(1+ ηiη) (4.3.1.5) 
 [ ] /
#
.
-
"
,
=
4321
4321
φφφφ
φφφφ
φ
vvvv
uuuu
u NNNN
NNNN
N (4.3.1.6) 
 )()( 21
12
41
41
1 YYNYYNN uuu −+−= φφφ (4.3.1.7) 
 )()( 32
23
12
12
2 YYNYYNN uuu −+−= φφφ (4.3.1.8) 
 )()( 43
34
23
23
3 YYNYYNN uuu −+−= φφφ (4.3.1.9) 
 )()( 14
41
34
34
4 YYNYYNN uuu −+−= φφφ (4.3.1.10) 
 )()( 12
12
14
41
1 XXNXXNN uuv −+−= φφφ (4.3.1.11) 
Elemento Cascarón Capítulo 4 
 59
 )()( 23
23
21
12
2 XXNXXNN uuv −+−= φφφ (4.3.1.12) 
 )()( 34
34
32
23
3 XXNXXNN uuv −+−= φφφ (4.3.1.13) 
 )()( 41
41
34
34
4 XXNXXNN uuv −+−= φφφ (4.3.1.14) 
 )1)(1(
8
1 212 ηξφ −−=uN (4.3.1.15) 
 )1)(1(
8
1 234 ηξφ +−=uN (4.3.1.16) 
 )1)(1(
8
1 223 ξηφ +−=uN (4.3.1.17) 
 )1)(1(
8
1 241 ξηφ −−=uN (4.3.1.18) 
 [ ]
/
/
/
/
/
#
.
-
-
-
-
-
"
,
=
xy
y
x
S
∂
∂
∂
∂
∂
∂
∂
∂
0
0
 (4.3.1.19) 
 { }
$%
$
&
'
$(
$
)
*−
=
x
yR
∂
∂
∂
∂
2
1
 (4.3.1.20) 
Para el vector de esfuerzos, los términos son los siguientes: 
Nx: Esfuerzo σx multiplicado por el espesor del elemento 
Ny: Esfuerzo σy multiplicado por el espesor del elemento 
Nxy: Esfuerzo γxy multiplicado por el espesor del elemento 
τ0: Esfuerzo residual de corte, este parámetro es arbitrario utilizado en la formulación de la matriz 
de rigidez del elemento cascarón de la referencia [6] para la parte de membrana, se utiliza en los cálculos 
Elemento Cascarón Capítulo 4 
 60
de las fuerzas internas, sin embargo en los ejemplos que se comprobaron su aportación fue insignificante 
con respecto a la aportación de los esfuerzos de membrana. 
De la misma referencia se tomó la regla de integración numérica de 5 puntos, como se muestra en la 
Figura 4.3.1-1 y en la tabla 4.3.1-1, ya que la matriz de rigidez de membrana del elemento es evaluada 
usando la misma regla de integración y si se quiere conservar la misma precisión es necesario evaluar 
los esfuerzos σσσσ y la matriz B en los mismos puntos. En las pruebas hechas en el presente trabajo, se 
evaluaron σσσσ y B para una regla de integración de 2x2, 3x3 y 4x4 puntos de muestreo y comparando los 
resultados con una multiplicación de la matriz de rigidez por los desplazamientos Kd; la regla de 2x2 
generaba un estado de equilibrio, sin embargo producía un estado de fuerzas que conforme avanzaba la 
simulación dinámica hacía el algoritmo inestable, se aumentaba la precisión numérica y la estabilidad