INSTITUTO TECONOLÓGICO Y DE ESTUDIOS SUPERIORES DE ...

90
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

Transcript of INSTITUTO TECONOLÓGICO Y DE ESTUDIOS SUPERIORES DE ...

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 ....................................................................................................................................65

5.- 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 analizables aplicando

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

tdta

∂∂= (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

tdtV

∂∂ )(

)( = (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)

nmedt∆ : 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)

tmedn : Tiempo medio del intervalo n, definido por:

2

1−−=nn

nmed

ttt (2.1.2)

Despejando la aceleración de la ecuación (2.1) obtenemos:

)(1 ndamp

nnet

n ffa −= −M (2.1.3)

donde:

an: Aceleración del intervalo n.

M: Matriz de masa.

nnetf : Fuerza neta del intervalo n, definida por:

nnext

nnet fff int−= (2.1.4)

nextf : Fuerza externa del intervalo n.

nf int : Fuerza interna del intervalo n.

ndampf : 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:

nmed

nnn

med t

ddV

∆−=

−1

(2.1.5)

donde:

nmedV : 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 media discreta

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

nmed

nmedn

t

VVa

∆−

=+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

11

−−

−−

=nn

med

nnmedn

tt

VVa (2.1.7)

nmed

n

nmed

nn

tt

VVa

−−

= (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 −−− −+= nnmed

nnnmed ttaVV (2.1.7b)

( )nmed

nnnmed

n 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)

nnext

nnet 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

endamp 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) ( )ndamp

nnet

n 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 −−− −+= nnmed

nnnmed 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

nnet 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

endamp 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) ( )ndamp

nnet

n ffa −= −1M

13. Se resuelve la velocidad al tiempo tn con la ecuación (2.1.8b) ( )nmed

nnnmed

n 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:

( )eee

ecritt ξξ

ω−+≤∆ 1

2min 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

oe

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

eintK : 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 tipo barra, 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)

ρσ xE

c−

= (2.3.4b)

l

c2max =ω (2.3.5)

l

c32max =ω (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

υρ += E

c (2.3.6)

t

c2max =ω (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

next

Tnnnext

next ffddWW (2.3.9)

Método de Integración Explícito Capítulo 2

16

)()(2

1 111 −−− −−+= ndamp

ndamp

Tnnndamp

ndamp ffddWW (2.3.10)

)()(2

1 nTnnkin VMVW = (2.3.11)

Donde nWint , nextW , 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 tn 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

iii 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 TTδδ (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 TTTδδ 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:

xu

∂∂ε = (2.4.2.1)

tomando la definición de los desplazamientos de (2.4.1.1b) y sustituyendo en (2.4.2.1) obtenemos:

BuuN ==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 funciones de 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 =001

2334

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

endamp 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

oi

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:

dSudquWS

TT ++ +Ω=Ω

τδδδ (2.4.4.1)

Sustituyendo (2.4.1.1b):

dSdqWS

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:

dSdqfS

TText ++ +Ω=

Ω

τ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 GRABAR CADA [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:

dxNNAML

Te += ρ (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

NNN (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( ξξ iiN

+= (3.2.3)

Cambiando el diferencial de longitud al espacio natural isoparamétrico tenemos:

Elemento Barra Capítulo 3

34

ξξ

dL

dxLXX

d

dx

22212 =∴=

−= (3.2.4)

Integrando la ecuación (3.2.1) obtenemos:

/#

.-"

,21

12

6ALρ

(3.2.5)

////

#

.

----

"

,

2010

0201

1020

0102

6ALρ

(3.2.5b)

////////

#

.

--------

"

,

200100

020010

002001

100200

010020

001002

6ALρ

(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

2ALρ

(3.2.6)

////

#

.

----

"

,

1000

0100

0010

0001

2ALρ

(3.2.6b)

////////

#

.

--------

"

,

100000

010000

001000

000100

000010

000001

2ALρ

(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 matrices fueron 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:

dxBAfL

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

NB

∂∂

∂∂ 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

121

1

2int AdxL

Af

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 MASDimensiones 2Numero de elementos 5Numero de nodos 6Materiales 1NOD_Elemento 2GDL Nodal 2

fin!Malla

Incidencias1 1 22 2 33 3 44 4 55 5 6

fin!

Coordenadas1 0.00000 0.000002 0.10000 0.000003 0.20000 0.000004 0.30000 0.000005 0.40000 0.000006 0.50000 0.00000

fin!

Propiedades!! Propiedades para la barra elástica de dosnodos

Numero 1 Tipo 1 Propiedades E2E+11 A 1.5E-4 Ro 7850.0 Masa Concentrada SalidaGauss!!

Elemento 1 Propiedad 1Elemento 2 Propiedad 1Elemento 3 Propiedad 1Elemento 4 Propiedad 1Elemento 5 Propiedad 1

fin!! Inicia definicion de los juegos de restricciones!

Restricciones ApoyoPuntuales

1 U = 01 V = 01 V = 02 V = 03 V = 04 V = 05 V = 06 V = 0

finfin

!! Definición de condiciones iniciales!

Condicion Inicial DespUnitDesplazamientos

2 U = 0.0023 U = 0.0044 U = 0.0065 U = 0.0086 U = 0.01

finfin

!! Definición de juegos de cargas!

Carga MasaConc!

!! Masas concentradas!

Nodal6 U = Masa 1000

fin

!! Termina definicion de cargas

fin! Termina comandos de mallafin! -----------------------------------------------------------Comandos

! Pasos previos para empezar analisis dinamicoControlAnalisis DinamicoMetodo ExplicitoTolerancia 0.01A0 0.0A1 0.0Grabar Cada 50 IteracionesIgnorar EsfuerzosIgnorar Deformacionesfin

Selecciona Restricciones ApoyoSelecciona Carga MasaConcSelecciona Condicion Inicial DespUnitMasa TodosRigidez TodosEnsambla

Inicializa ProcesoITERA 6501incrementa pasovelocidad mediadesplazamientosFuerzasAmortiguamientoaceleracionesvelocidadbalanceTERMINACierra ProcesoExporta GID

finfin

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:

( )tPdt

d

t

d =++ 400005.85.22

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 −++−= −− t

tCosetSened 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 −−++−= −− t

tCosetSened 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:

ttSen

d 5.122.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 BARRADimensiones 2Numero de elementos 1Numero de nodos 2Materiales 1NOD_Elemento 2GDL Nodal 2

fin! ---------------------------------------------------------------Malla

Incidencias1 1 22 2 33 3 44 4 55 5 6

fin!

Coordenadas1 0.00000 0.000002 500.00000 0.00000

fin!

! Definicion de FuncionesFuncion Fun1Archivo "c:\disco_local\P.txt"fin

Propiedades!! Propiedades para la barra elástica de dos! nodos

Numero 1 Tipo 1 Propiedades E200000 A 10.0 Ro 0.001 Masa Concentrada Salida Gauss

!!

Elemento 1 Propiedad 1fin

!! Inicia definicion de los juegos de restricciones!

Restricciones ApoyoPuntuales

1 U = 01 V = 02 V = 0

finfin

!! Definición de juegos de cargas!

Carga P!!! Funciones nodales!

Nodal2 U = funcion Fun1 Amplificacion 1.0

Desfazamiento 0.0fin

!! Termina definicion de cargas

fin! Termina comandos de mallafin! -----------------------------------------------------------Comandos

! Pasos previos para empezar analisis dinamicoControlAnalisis DinamicoMetodo ExplicitoTolerancia 0.01A0 0.02A1 0.002Grabar Cada 0.005 SegundosIgnorar EsfuerzosIgnorar Deformacionesfin

Selecciona Restricciones ApoyoSelecciona Carga PMasa TodosRigidez TodosEnsambla

Inicializa ProcesoITERA 301incrementa pasovelocidad mediadesplazamientosFuerzasAmortiguamientoaceleracionesvelocidadbalanceTERMINACierra ProcesoExporta GID

finfin

30.0 0.00.1 5000.00.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)

00000

1

2

33333

4

5+

012

345 +

+−= !

= L

tE

n

CosL

xnSen

nn

n

2

)12(

2

)12(

)12(

)1(8

022

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 Barra2Dimensiones 2Numero de elementos 10Numero de nodos 11Materiales 1NOD_Elemento 2GDL Nodal 2

fin! ---------------------------------------------------------------Malla

Incidencias1 1 22 2 33 3 44 4 55 5 66 6 77 7 88 8 99 9 10

10 10 11fin

!Coordenadas

1 0.00000 0.000002 0.15000 0.000003 0.30000 0.000004 0.45000 0.000005 0.60000 0.000006 0.75000 0.000007 0.90000 0.000008 1.05000 0.000009 1.20000 0.00000

10 1.35000 0.0000011 1.50000 0.00000

fin!

Propiedades!! Propiedades para la barra elástica de dosnodos

Numero 1 Tipo 1 Propiedades E2E+11 A 0.00015 Ro 7850 Masa Concentrada Salida Gauss

!!

Elemento 1 Propiedad 1Elemento 2 Propiedad 1Elemento 3 Propiedad 1Elemento 4 Propiedad 1Elemento 5 Propiedad 1Elemento 6 Propiedad 1Elemento 7 Propiedad 1Elemento 8 Propiedad 1Elemento 9 Propiedad 1Elemento 10 Propiedad 1

fin!! Inicia definicion de los juegos de restricciones!

Restricciones ApoyoPuntuales

1 U = 01 V = 01 V = 02 V = 03 V = 04 V = 05 V = 06 V = 06 V = 07 V = 08 V = 09 V = 0

10 V = 011 V = 0

finfin

!! Definición de condiciones iniciales!

Condicion Inicial dtoDesplazamientos

2 U = 0.00013 U = 0.00024 U = 0.00035 U = 0.00046 U = 0.00057 U = 0.00068 U = 0.00079 U = 0.0008

10 U = 0.000911 U = 0.001

finfin

!! Definición de juegos de cargas!

Carga Fuerza!

!! Termina definicion de cargas

fin! Termina comandos de mallafin! -----------------------------------------------------------Comandos

! Pasos previos para empezar analisis dinamicoControlAnalisis DinamicoMetodo ExplicitoTolerancia 0.01A0 0.0A1 0.0Grabar Cada 0.000020 SegundosIgnorar EsfuerzosIgnorar Deformacionesfin

Selecciona Restricciones ApoyoSelecciona Carga FuerzaSelecciona Condicion Inicial dtoMasa TodosRigidez TodosEnsambla

Inicializa ProcesoITERA 501incrementa pasovelocidad mediadesplazamientosFuerzasAmortiguamientoaceleracionesvelocidadbalanceTERMINACierra ProcesoExporta GID

finfin

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 de 1,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

)()( 2112

4141

1 YYNYYNN uuu z−+−= φφφ

(4.2.4)

)()( 3223

1212

2 YYNYYNN uuu z−+−= φφφ

(4.2.5)

)()( 4334

2323

3 YYNYYNN uuu z−+−= φφφ

(4.2.6)

)()( 1441

3434

4 YYNYYNN uuu z−+−= φφφ

(4.2.7)

)()( 1212

1441

1 XXNXXNN uuv z−+−= φφφ

(4.2.8)

Elemento Cascarón Capítulo 4

55

)()( 2323

2112

2 XXNXXNN uuv z−+−= φφφ

(4.2.9)

)()( 3434

3223

3 XXNXXNN uuv z−+−= φφφ

(4.2.10)

)()( 4141

3434

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, ==! iMmi

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

AAxy

y

xT ττ++ +

$%

$&

'

$(

$)

*

=int (4.3.1.1)

donde:

[ ] [ ] [ ]" #φuu NNSB ][= (4.3.1.2)

[ ] [ ] [ ]" # [ ]" #φφτ NNNRB uuT −+= 0 (4.3.1.3)

[ ] [ ] /#

.-"

,=

4321

4321

0000

0000,

uuuu

uuuuu NNNN

NNNNNN φ (4.3.1.4)

Nui =1

4(1 +ξ iξ )(1+ ηiη) (4.3.1.5)

[ ] /#

.-"

,=

4321

4321

φφφφ

φφφφφ

vvvv

uuuuu NNNN

NNNNN (4.3.1.6)

)()( 2112

4141

1 YYNYYNN uuu −+−= φφφ (4.3.1.7)

)()( 3223

1212

2 YYNYYNN uuu −+−= φφφ (4.3.1.8)

)()( 4334

2323

3 YYNYYNN uuu −+−= φφφ (4.3.1.9)

)()( 1441

3434

4 YYNYYNN uuu −+−= φφφ (4.3.1.10)

)()( 1212

1441

1 XXNXXNN uuv −+−= φφφ (4.3.1.11)

Elemento Cascarón Capítulo 4

59

)()( 2323

2112

2 XXNXXNN uuv −+−= φφφ (4.3.1.12)

)()( 3434

3223

3 XXNXXNN uuv −+−= φφφ (4.3.1.13)

)()( 4141

3434

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

xS

∂∂

∂∂

∂∂

∂∂

0

0

(4.3.1.19)

$%

$&'

$(

$)*−

=x

yR

∂∂

∂∂

21

(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 de

la rutina a medida que se cambiaba la regla de integración a 3x3 y 4x4 puntos donde se optó por usar la

misma regla de integración de 5 puntos que con menos iteraciones logra una mejor precisión.

Figura 4.3.1-1 Regla de integración con 5 puntos de muestreo

4

1 0WW −=α (4.3.1.21)

Elemento Cascarón Capítulo 4

61

α

αW3

1= (4.3.1.22)

con W0=0.01

i ξξξξ ηηηη W 1 0.00000 00000 00000 0.00000 00000 00000 0.99750

2 -0.57807 33130 16080 -0.57807 33130 16080 0.99750

3 +0.57807 33130 16080 -0.57807 33130 16080 0.99750

4 -0.57807 33130 16080 +0.57807 33130 16080 0.99750

5 +0.57807 33130 16080 +0.57807 33130 16080 0.99750

Tabla 4.3.1-1 Posición de los puntos de integración y valor de factor de peso para regla de 5

puntos

La integración de la parte de la membrana reporta las siguientes fuerzas:

$$$$$$$$

%

$$$$$$$$

&

'

$$$$$$$$

(

$$$$$$$$

)

*

=

4

4

4

3

3

3

2

2

2

1

1

1

int

z

y

x

z

y

x

z

y

x

z

y

x

M

f

f

M

f

f

M

f

f

M

f

f

f (4.3.1.23)

Elemento Cascarón Capítulo 4

62

4.3.2.- Fuerzas Internas de Placa

La ecuación (2.4.2.4) para la parte de placa usando la formulación del elemento cascarón de la referencia

[6] se obtiene la siguiente expresión:

[ ] [ ] dA

M

M

M

BdAQ

QBf

xy

y

xT

AA y

xT

$%

$&

'

$(

$)

*

+%&'

()*

= ++ φγint (4.3.2.1)

donde:

[ ] [ ][ ] [ ]/#

.-"

,−∇

∇=

φφγ NN

NB

w

w (4.3.2.2)

[ ] [ ][ ]/#.

-"

,−

φ NLB

0 (4.3.2.3)

[ ] [ ] /#

.-"

,=

4321

4321

0000

0000,

uuuu

uuuuw NNNN

NNNNNN φ (4.3.2.4)

)1)(1(4

1 ηηξξ iiuiN ++= (4.3.2.5)

[ ] " # " # " # " #[ ]1111 φφφφφ wwwww NNNNN = (4.3.2.6)

" # " # " #)()()()( 212112

414141

1 YYXXNYYXXNN www −−+−−= φφφ (4.3.2.7)

" # " # " #)()()()( 323223

121212

2 YYXXNYYXXNN www −−+−−= φφφ (4.3.2.8)

" # " # " #)()()()( 434334

232323

3 YYXXNYYXXNN www −−+−−= φφφ (4.3.2.9)

" # " # " #)()()()( 141441

343434

4 YYXXNYYXXNN www −−+−−= φφφ (4.3.2.10)

)1)(1(16

1 212 ηξφ −−−=uN (4.3.2.11)

)1)(1(16

1 223 ξηφ +−−=uN (4.3.2.12)

Elemento Cascarón Capítulo 4

63

)1)(1(16

1 234 ηξφ +−−=uN (4.3.2.13)

)1)(1(16

1 241 ξηφ −−−=uN (4.3.2.14)

$%

$&'

$(

$)*

=∇y

x

∂∂

∂∂

(4.3.2.15)

[ ]

/////

#

.

-----

"

,

=

xy

y

xL

∂∂

∂∂

∂∂

∂∂

0

0

(4.3.2.16)

Donde los esfuerzos usados son los siguientes:

Qx: Esfuerzo cortante Sx multiplicado por el espesor del elemento

Qy: Esfuerzo cortante Sy multiplicado por el espesor del elemento

Mx: Esfuerzo de flexión en la cara x multiplicado por el espesor del elemento

My: Esfuerzo de flexión en la cara y multiplicado por el espesor del elemento

Mxy: Esfuerzo de alabeo multiplicado por el espesor del elemento

Para la integración numérica se usó al igual que la referencia una integración completa con 3x3 puntos

de gauss mostrados en la tabla 4.3.2-1.

i ξξξξ ηηηη Wξξξξ Wηηηη

1 -0.77459 66692 41483 -0.77459 66692 41483 0.55555 55555 55556 0.55555 55555 55556

2 0.00000 00000 00000 0.55555 55555 55556 0.88888 88888 88889

3 +0.77459 66692 41483 0.55555 55555 55556 0.55555 55555 55556

4 0.00000 00000 00000 -0.77459 66692 41483 0.88888 88888 88889 0.55555 55555 55556

5 0.00000 00000 00000 0.88888 88888 88889 0.88888 88888 88889

6 +0.77459 66692 41483 0.88888 88888 88889 0.55555 55555 55556

7 +0.77459 66692 41483 -0.77459 66692 41483 0.55555 55555 55556 0.55555 55555 55556

8 0.00000 00000 00000 0.55555 55555 55556 0.88888 88888 88889

9 +0.77459 66692 41483 0.55555 55555 55556 0.55555 55555 55556

Tabla 4.3.2-1 Posición de los puntos de muestreo y su factor de peso para regla de integración 3x3

Elemento Cascarón Capítulo 4

64

La parte de placa reporta las siguientes fuerzas en el sentido local del plano del cascarón

$$$$$$$$

%

$$$$$$$$

&

'

$$$$$$$$

(

$$$$$$$$

)

*

=

4

4

3

3

2

2

1

1

4

3

2

1

int,

x

y

x

y

x

y

x

y

z

z

z

z

ne

M

M

M

M

M

M

M

M

f

f

f

f

f (4.3.2.17)

La Figura 4.3.2-1 muestra el sentido de los momentos de placa, los cuales deben de ser convertidos al

sentido cartesiano local del elemento placa como se muestra en el vector (4.3.2.17)

Figura 4.3.2-1 Fuerzas de placa en el sentido local del elemento

Elemento Cascarón Capítulo 4

65

4.4.-Validación

Ahora que el procedimiento de integración explícita fue probado con las pruebas de la barra, se

procederá a validar el elemento con algunos ejemplos comparados con su solución analítica en caso de

haberla y con la solución numérica tomada de ejemplos de las referencias.

Ejemplo 4.4.1

Se analiza una viga de sección transversal rectangular doblemente empotrada modelada con 10

elementos finitos tipo cascarón sujeta a una carga uniformemente distribuida en su longitud. Con la

forma en que es modelada la viga se puede validar principalmente el esfuerzo de membrana del

elemento cascarón. Cuando un elemento se le excita con una acción constante empieza a oscilar en su

primer modo de vibración, el cual corresponde a la masa propia del elemento, lo anterior es válido ya

que la fuerza que se le imprime para la simulación no forma parte de la masa del elemento.

Figura 4.4-1 Viga doblemente empotrada sujeta a carga uniformemente distribuida

De la referencia [4] se toma el período de vibración de este elemento, que se determina mediante:

4AL

EIann ρ

ω = (4.4.1)

donde n representa el número de período, para una viga con estas características los dos primeros

períodos de vibración le corresponde los factores a1 = 22.0 y a2=61.7.

Elemento Cascarón Capítulo 4

66

Datos del problema:

b : 0.20 m Ancho de la viga h : 0.40 m Peralte de la viga L : 4.00 m Longitud de la viga

E : 2.10E+09 kg/m2 Módulo de Young

ρ : 800 kg s2 / m4 Masa de la viga

A : 0.08 m2 Área transversal I : 1.07E-03 m4 Momento de Inercia ρ A : 64 kg s2 / m2

a1 : 22.0 ω1 : 257.24 rad/seg

f : 0.0244 seg

a2 : 61.7 ω2 : 721.44 rad/seg f : 0.0087 seg

Se plantea la solución estática del sistema con el fin de verificar hacia donde converge la solución con

amortiguamiento. Para la solución de equilibrio estático se tienen los siguientes parámetros:

P : -10.00 kg/m2 Carga sobre la viga W : -2.00 kg/m Carga distribuída sobre la viga δmax : -5.95E-07 m Deflexión estática

Para la solución numérica, se presenta los datos con oscilación libre y con oscilación amortiguada, se

propuso un amortiguamiento del 2% relativos a los 2 primeros períodos de oscilación

ξ1 : 2.0% ξ2 : 2.0%

Coeficientes de Amortiguamiento a0 : 7.59 rad/seg a1 : 0.00004 seg3

Elemento Cascarón Capítulo 4

67

Donde los coeficientes de amortiguamiento se determinan resolviendo el par de ecuaciones que resultan

de aplicar la fórmula (2.3.2)22

10 n

nn

aa ωω

ξ +=

Figura 4.4-2 Resultados numéricos del Ejemplo 4.4.1

De los resultados numéricos se obtiene lo siguiente:

%error

ω : 251.61 rad/seg 2.19% f : 0.0250 seg

En la Figura 4.4-2 se puede apreciar que el resultado sin amortiguamiento oscila alrededor de la

deflexión estática mientras que el resultado con amortiguamiento tiende a dejar de oscilar y quedar

deformado en la solución estática de equilibrio

Elemento Cascarón Capítulo 4

68

Ejemplo 4.4.2

Tomemos ahora el caso de una losa cuadrada de espesor constante simplemente apoyada en sus orillas

sometida a una carga uniformemente distribuida en su superficie, el modo fundamental de vibración

tomado de la referencia [4] está dado por:

t

D

l ρπω 0

12

345=

22 2

(4.4.2)

donde:

( )2

3

112 υ−= Et

D (4.4.3)

ρ: Densidad del material

t: Espesor de la losa

E: Módulo de Young

υ: Coeficiente de Poisson

l: Lado de la losa

Datos del Problema:

L : 4.00 m t : 0.25 m E : 2.53E+09 kg/m2 ρ : 245 kg s2 / m4 ν : 0.30

D : 3.62E+06 kg-m

ω : 300 rad/seg T : 48 cps f : 2.09E-02 s

Elemento Cascarón Capítulo 4

69

Se modeló solamente un cuarto del problema con una malla de 16 elementos(4 por 4) aprovechando las

condiciones de simetría. Con esta modelación las fuerzas principales de los elementos finitos son fuerzas

de placa. Al igual que el ejemplo 4.4-1 se compara la oscilación libre con la amortiguada y con la

deflexión estática al centro de la losa. De un análisis estático y con carga distribuida se obtuvo el

siguiente valor de deflexión al centro:

P : -300.00 kg/m2 Carga uniformemente distribuida δ : -9.19E-05 m Deflexión al centro de la losa Se hicieron dos análisis dinámicos, uno con oscilación libre y el otro con amortiguamiento con los

siguientes datos:

a0 : 5.00 rad/seg a1 : 5.50E-06 seg3

Estos valores disminuyen el incremento de tiempo crítico en un 83%

La Figura 4.4-3 muestra la comparativa de los resultados:

Figura 4.4-3 Resultados numéricos del Ejemplo 4.4.2

Elemento Cascarón Capítulo 4

70

De la solución del problema podemos obtener el siguiente valor de frecuencia y compararlo con el valor

analítico de la solución:

T : 2.13E-02 s error : 1.71%

En la Figura 4.4-3 se puede apreciar lo mismo que la Figura 4.4-2, que el resultado sin amortiguamiento

oscila alrededor de la deflexión estática mientras que el resultado con amortiguamiento tiende a dejar de

oscilar y quedar deformado en la solución estática de equilibrio

Ejemplo 4.4.3

El siguiente ejemplo se tomó de la referencia [6], este ejemplo es muy utilizado en la literatura para

comparar diferentes elementos y formulaciones en régimen elastoplástico, la forma más sencilla de

analizarlo es con sólidos de revolución debido a que es un problema axisimétrico, tal y como fue hecho

originalmente. En este capítulo se analizó con elementos cascarón. La figura 4.4-4 muestra los datos del

problema y la figura 4.4-5 muestra el modelo de elementos finitos utilizado. Se modeló una cuarta parte

del domo con 48 elementos tipo cascarón y sometiéndolo a restricciones de simetría en lo planos X y Y

Figura 4.4-4 Domo bajo carga uniforme

Elemento Cascarón Capítulo 4

71

Figura 4.4-5 Modelo de elementos finitos del Ejemplo 4.4.3

No existe una solución analítica para este problema y sólo se disponen de algunos resultados numéricos

obtenidos de soluciones con diferentes elementos finitos. En la figura 4.4-5 se grafica el historial del

desplazamiento del ápice del domo en función del tiempo. Para este caso elástico, el comportamiento del

desplazamiento es similar al analizado en la referencia [7] con diferentes elementos finitos.

Elemento Cascarón Capítulo 4

72

Figura 4.4-6 Desplazamiento al centro del domo del Ejemplo 4.4.3

Ejemplo 4.4.4

Se analiza el caso de una viga en cantiliver tomada de la referencia [7], con una carga puntual en su

extremo libre produciendo flexión compuesta y esfuerzo axial, los datos de problema y el modelo de

elementos finitos utilizado se muestra a continuación:

Figura 4.4-7 Viga en cantiliver con carga puntual en el extremo libre

Elemento Cascarón Capítulo 4

73

Datos del problema:

b : 0.05 h : 0.10 L : 1.00 E : 2.00E+06 ρ : 7800 Componentes de la carga puntual en el extremo libre

P : -1.0 en X 0.1 en Y 1.0 en Z

La viga se modeló con 10 elementos cascarón, en este ejemplo el elemento finito trabaja tanto con

fuerzas de membrana como con fuerzas de placa. En la Figura 4.4-8 se muestra el desplazamiento del

extremo libre logrando una excelente similitud a los resultados de la referencia:

Figura 4.4-8 Resultados numéricos del Ejemplo 4.4.4

Elemento Cascarón Capítulo 4

74

Figura 4.4-9 Resultados del Ejemplo 4.4.4 en la Referencia [7]

En la Figura 4.4-9 se muestra el resultado en la Referencia [7], U1, U2 y U3 significan Desplazamiento

X, Y y Z respectivamente

Conclusiones Capítulo 5

75

5.- Conclusiones

En dinámica estructural la mayoría de los problemas en los que es posible encontrar una solución

analítica son ejemplos muy elementales; pero cuando se trata de encontrar una solución a un problema

real el planteamiento analítico puede resultar sumamente complicado y hasta imposible de resolver, por

lo tanto se tiene que recurrir a soluciones numéricas donde la técnica de los elementos finitos es una

excelente herramienta. Aunque en la mayoría de la literatura se omita mencionar el uso de elementos

finitos para la solución de problemas de dinámica estructural, por lo general se plantea la solución

suponiendo que se cuenta con parámetros como matrices de masa, matrices de amortiguamiento, campo

de esfuerzos, etc. parámetros que son fácilmente determinables segmentando espacialmente un dominio

en elementos finitos, es por ello que la dinámica estructural está muy ligada con los elementos finitos.

En el presente trabajo se pudo observar que para implementar una solución numérica a un problema

dinámico en general valiéndose de herramientas de cómputo se debe de partir de dos planteamientos: por

un lado se debe de desarrollar la rutina de integración, tal como se plantea en la Sección 2.2, y traducirla

a algún lenguaje de programación, por otro lado se debe de preparar los elementos finitos con rutinas

propias para proveer de resultados necesarios al algoritmo de solución, como se muestra en la Sección

2.4. El algoritmo funcionando satisfactoriamente debe de ser capaz de incluir cualquier elemento finito

que tenga rutinas propias de análisis dinámico.

Se implementó un procedimiento de solución explícito, cuyo algoritmo es sencillo de entender por la

simplicidad que presentan sus cálculos. A medida que se fueron haciendo pruebas del algoritmo fueron

surgiendo nuevas limitantes para hacer el procedimiento estable numéricamente, como se detalla en la

Sección 2.3, se comprobó que el procedimiento es muy sensible a los cambios de incrementos de tiempo

aún cuando éstos estén por debajo de los incrementos de tiempo críticos, se analizó un caso particular

donde se alternaba el salto de tiempo crítico del problema: Para una iteración se utilizó el salto de

Conclusiones Capítulo 5

76

tiempo crítico y en la siguiente iteración se usó la mitad del mismo, siguiendo este mismo patrón la

solución convergía numéricamente al inicio de la simulación pero después de un determinado número de

iteraciones el algoritmo se volvió inestable reportando deformaciones excesivas en la estructura. Cuando

se continúe el presente trabajo adicionando no-linealidad a las relaciones esfuerzo-deformación los

saltos de tiempo se verán afectados ya que éstos son función del modo de vibración de los elementos

finitos que a su vez son función de la rigidez la cual cambia dependiendo del estado de esfuerzos del

elemento. Se propone implementar un método de solución implícito para tomar en cuenta efectos de no-

linealidad del material, el cual es incondicionalmente estable a los incrementos de tiempo.

Una tarea interesante es observar qué modo de vibración predecirá el algoritmo, por lo general, bajo una

condición de carga estable la estructura puede empezar a oscilar en su modo más elevado, sin embargo,

depende mucho de la forma de cargar la estructura, la distribución de la masa, lo periódico de la carga y

otros factores para que la estructura empiece a oscilar en otro de sus modos o combinaciones de éstos.

En la solución numérica, las rotaciones nodales son numéricamente despreciables comparados con los

desplazamientos, sin embargo, es importante tenerlas en cuenta para poder tener convergencia numérica,

especialmente en aquellos elementos finitos cuyos desplazamientos y rotaciones van ligados.

Cuando se programe en algún lenguaje computacional las rutinas propias de los elementos finitos, es

importante generar pruebas, especialmente con el vector de fuerzas internas y de ser posible compararlas

con la solución estática, ya sea el caso lineal o no-lineal. En el presente trabajo se realizaron algunas

pruebas de membrana utilizando una regla de integración de fuerzas de 2x2 puntos de muestreo, aunque

el elemento estaba en equilibrio estático en cada iteración reportaba pequeños incrementos en las fuerzas

comparando la solución estática para el mismo estado de deformación, estas fuerzas al principio

imperceptible a medida que la simulación avanzaba provocaban deformaciones fuera de control, en este

caso fue un problema de inestabilidad numérica debido al elemento finito (ver Sección 4.3.1 y 4.3.2)

Conclusiones Capítulo 5

77

Al igual que un análisis estático con elementos finitos, se puede obtener mejor precisión para ciertas

zonas si se refina el malleo de elementos, sin embargo el tiempo de cómputo se eleva

considerablemente, ya que por un lado son más elementos por calcular, por otro lado el salto de tiempo

crítico es función también de las dimensiones de los elementos. Obtener una simulación de unos cuantos

segundos puede resultar de varios miles de iteraciones. Es recomendable entonces contar también con

equipos de cómputo veloces en sus procesos, con suficiente memoria para el cálculo de las simulaciones

y con buen espacio en disco para almacenar resultados. Aunque para el cálculo de las variables para un

determinado tiempo solamente es necesario conocer los valores de las mismas variables un lapso de

tiempo anterior, es importante guardar los resultados del historial para apreciar la simulación.

Contribuciones y Trabajo Futuro Capítulo 6

78

6.- Contribuciones y trabajo futuro

Al programa de análisis de elementos finitos PAEF se le ha adicionado el código necesario para poder

generar un análisis dinámico lineal, pudiendo simular problemas con desplazamientos nodales iniciales,

cargas fijas y cargas nodales dependientes del tiempo, a su vez se complementó dos elementos finitos

existentes en el programa, el elemento barra y el elemento cascarón, con las rutinas necesarias para

alimentar los parámetros requeridos en un análisis dinámico.

Con las rutinas ahora existentes se pueden modelar diversos problemas dinámicos, tales como:

1. Propagaciones de ondas en medios continuos, como chasis de carros, sistema de pisos, etc.

2. Simulación de cargas sísmicas en estructuras de tipo cascarón, tales como domos de concreto,

tanques de almacenamiento de agua, losas de entrepisos, etc.

3. Revisión de vibraciones en sistemas de pisos a base de Joist y Joist Girders con losa de concreto,

como lo muestra la Referencia [14] y [17].

Aunque se implementó un análisis dinámico lineal con resultados comprobados aceptables, el programa

puede ser complementado con otras rutinas, como podrían ser no-linealidad en material y no-linealidad

geométrica. Para implementar la no-linealidad en el material se debe de trabajar en las relaciones

constitutivas de esfuerzo deformación en la rutina de fuerzas internas y ayudándose con algún catálogo

de materiales, además de poder incluir un historial de deformaciones con respecto al tiempo. Ahora, para

implementar la no-linealidad geométrica debe incluirse en las rutinas de determinación de

desplazamientos un procedimiento iterativo que contemple los pasos necesarios para tal efecto. Un

método de integración implícito sería lo más conveniente si se pretende implementar no-linealidad en el

material, de esa manera el salto de tiempo en las soluciones puede ser más amplio con resultados

estables numéricamente, por las razones expuestas en las conclusiones y en la sección 2.3

Contribuciones y Trabajo Futuro Capítulo 6

79

A su vez, es conveniente encontrar una forma de evaluar el incremento de tiempo crítico en cascarones

incluyendo esfuerzos cuando se pretenda tomar en cuenta no-linealidad.

El programa PAEF en su inicio fue pensado para ser didáctico, sin embargo, para problemas con

considerables números de grados de libertad el análisis estático y sobre todo el dinámico se torna lento,

es conveniente entonces implementar algunas herramientas numéricas para hacer más eficiente el tiempo

de cómputo, la memoria y el espacio de almacenamiento, lo anterior es independiente del tipo de análisis

a generar, pudiendo ser las siguientes:

1. Rutinas de estructura de datos: Útiles en la búsqueda de datos en listas, disminuye los tiempos de

búsqueda resultando de soluciones más rápidas.

2. Almacenamiento de datos importantes: Útiles para hacer más eficiente el espacio de a

lmacenamiento de datos, por ejemplo, las matrices pueden guardar sólo su ancho de banda o en el

caso de matrices simétricas solamente su diagonal superior o inferior.

3. Programación orientada a objetos con memoria dinámica: Permite mantener los datos importantes en

memoria en lugar de estar llamándolos desde memoria de almacenamiento, resultando tiempos de

cómputo más rápidos.

Referencias Capítulo 7

80

7.-Referencias

1. Ted Belytschko,Wing Kam Liu y Brian Moran, Nonlinear finite elements for continua and structures, John Wiley & Sons, (2000).

2. Miguel Xicotencatl Rodriguez Paz, Análisis no-lineal y dinámico por elementos finitos: Desarrollo e implementación de los módulos de análisis para el programa “PAEF”, Tesis de maestría en ciencias en la especialidad en ingeniería civil, Instituto Tecnológico y de Estudios Superiores de Monterrey, (1996).

3. Klaus-Jürgen Bathe, Numerical methods in finite element analysis, Prentice Hall,(1996).

4. J.P Den Hartog, Mechanical Vibrations, 4th ed., McGraw Hill, (1956).

5. Ray W. Clough y Joseph Penzien, Dynamics of structures, 2nd ed, McGraw Hill, (1993).

6. Sergio Gallegos Cázares y Mariano Morán Guillaumin, “Un elemento cascarón cuadrilareral estable”, Memorias del segundo congreso internacional en métodos numéricos en ingeniería y ciencias aplicadas, V. 1, p.589, (2002).

7. F. Flores y E. Oñate, “Análisis dinámico de estructuras de láminas y vigas”, Publicación del Centro Internacional de Métodos Numéricos en Ingeniería No. 39, Octubre (1993).

8. A. Tabiei, Romil Tanov, “Sandwich shell finite element for dynamic explicit analyisis”, Int. J. Numer. Meth. Engng, 54,763-787 (2002).

9. Sergio Gallegos Cázares, “Apuntes de la clase de Elementos Finitos”, Instituto Tecnológico y de Estudios Superiores de Monterrey, (2002).

10. Tomas J.R. Hughes, The finite element method, linear static and dynamic finite element analysis, Prentice-Hall, (1987).

11. Robert Davis Cook, Concepts and applications of finite element analysis, 2nd ed, Wiley, (1981)

12. Singiresu S. Rao, Mechanical Vibrations 2nd Ed., Addison-Wesley, (1990).

13. Kenneth H. Lenzen, “Vibration of Steel Joist-Concrete Slab Floors”, American Institute of Steel Construction Inc, Volume 3, No. 3, 133-137 (1966).

14. Richard N. Wright y William H. Walker, “Vibration and Deflection of Steel Bridges”, American Institute of Steel Construction Inc, Volume 9, No. 1, (1972).

15. Thomas M. Murray, William E. Hendrick, “Floor Vibrations and Cantilivered Construction”, American Institute of Steel Construction Inc, Volume 14, No. 3, (1977).

16. Thomas M. Murray, “Acceptability Criterion for Occupant-Induced Floor Vibrations”, American Institute of Steel Construction Inc, Volume 18, No. 2, (1981).

Referencias Capítulo 7

81

17. Raed A Tolaymat, “A New Approach to Floor Vibration Analysis”, American Institute of Steel Construction Inc, Volume 25, No. 4, (1988).

18. Jeffrey A. Laman, “Design Aids for Walking Vibrations in Steel Framed Floors”, American Institute of Steel Construction Inc, Volume 36, No. 2, (1999).

19. M. Fouad Ahmad y Harry H. Hilton, “Finite Element Algorithms for Dynamic Simulation of Viscoelastic Composite Shell Structures Using Conjugated Gradient Method on Coarse Grained and Massively Parallel Machines”, Int. J. Numer. Meth. Engng, 40, 1857-1875, (1997).

20. Per-Anders Hansson y Göran Sandberg, “Mass Matrices by Minimization of Modal Errors”, Int. J. Numer. Meth. Engng, 40, 4259-4271, (1997).

21. T.C. Fung, “A Precise Time-Step Integration Method by Step-Response and Impulsive-Response Matrices for Dynamic Problems”, Int. J. Numer. Meth. Engng, 40, 4501-4527, (1997).

22. Adnan Ibrahimbegovic y Mazen Al Mikdad, “Finite Rotations in Dynamics of Beams and Implicit Time-Stepping Schemes”, Int. J. Numer. Meth. Engng, 41, 781-814, (1998)

23. Bostjan Brank, Lamberto Briseghella, Nicola Tonello y Frano B. Damjanic, “On Non-Linear Dynamics of Shells: Implementation of Energy-Momentum Conserving Algorithm For a Finite Rotation Sheel Model”, Int. J. Numer. Meth. Engng, 42, 409-442, (1998).

24. B. W. Golley, “A Weighted Residual Development of a Time-Stepping Algorithm for Structural Dynamics Using Two General Weight Functions”, Int. J. Numer. Meth. Engng, 42, 93-103, (1998).

25. Eugenio Oñate y Francisco Zárate, “Rotation-free triangular plate and shell elements”, Int. J. Numer. Meth. Engng, 47, 557-603, (2000).

26. S.Gopalakrishnan, “A deep rod finite element for structural dynamics and wave propagation problems', Int. J. Numer. Meth. Engng, 48, 731-744, (2000).

27. Antjony Gravouil y Alain combescure, “Multi-time-step explicit-implicit method for non-linear structural dynamics”, Int. J. Numer. Meth. Engng, 50, 199-225, (2001).

28. L. Liu, G. R. Liu y V. B. C. Tan, “Element free method for static and free vibration analysis of spatial thin shell structures”, Comp. Meth. in Applied Mech. and Engng, Volume 191, Issues 51-52, 5923-5942, (2002).

29. V. Balamurugan, S. Narayanan, “Shell finite element for smart piezoelectric composite plate/shell structures and its application to the study of active vibration control”, Finite Elements in Analysis and Design, Volume 47, Issue 9, 713-738, (2001).