Dinamica No Lineal (Sergio Oller)
Dinamica No Lineal (Sergio Oller)
Dinamica No Lineal (Sergio Oller)
S. Oller
Sergio Oller
Profesor de la Escuela Técnica Superior de Ingenieros
de Caminos, Canales y Puertos de Barcelona.
BARCELONA, España
Fractura mecánica. Un enfoque global
Sergio Oller
Primera edición, Enero 2001
Impreso por: Artes Gráficas Torres S.A., Morales 17, 08029 Barcelona, España
ISBN:
Depósito legal:
Presentación
Esta monografía trata la “dinámica no-lineal” de los sistemas estructurales. Existen diversos
enfoques para esta materia y por ello intento que este trabajo aporte un punto de vista más
al estudio dinámico no-lineal.
La motivación para escribir estas páginas se basa en la necesidad de contar con un material
de estudio para la asignatura de “Dinámica no-lineal” que dicto en doctorado del
departamento de “Resistencia de Materiales y Estructuras en la Ingeniería” de la
Universidad Politécnica de Cataluña.
Espero que estas notas ayuden a la mejor comprensión de la dinámica e incentiven a lector
a una mayor profundización.
Sergio Oller
VI DINAMICA NO-LINEAL
Contenido.
B 1 Introducción. .......................................................................................................1-1
B1.1 Presentación. .....................................................................................................................1-1
B1.1 Presentación.
La dinámica de estructuras estudia el equilibrio estructural a lo largo del tiempo entre las
acciones externas, las fuerzas elásticas, las fuerzas másicas y las fuerzas de
amortiguamiento, para un sistema estructural discreto en forma de puntos vinculados
internamente entre sí y todos ellos a un sistema de referencia fijo. Estos vínculos internos
entre los puntos que describen el sistema estructural pueden o no ser elástico, en el caso
que no lo sean, el comportamiento del sistema de puntos es no conservativo y por lo tanto
se dice que el material de la estructura tiene un comportamiento constitutivo no lineal disipativo.
Además de este comportamiento no lineal, también existe el comportamiento no lineal disipativo
por influencia de la viscosidad del material, que da lugar a fuerzas de amortiguamiento
dependientes de la velocidad del sistema. En algunos casos más simples, la no linealidad
por amortiguamiento se debe al desarrollo de fuerzas viscosas proporcionales a la
velocidad, pero en otros casos más complejos puede ocurrir que el propio término de
viscosidad es dependiente del tiempo. La no linealidad del sistema, también se manifiesta
en aquellos casos en que hay grandes movimientos y el sistema trabaja fuera de su
configuración geométrica inicial, motivando un comportamiento cinemático no lineal. Esta
última no linealidad se profundiza más en aquellos casos en que además de los grandes
movimientos, ocurren grandes deformaciones, situación que hace más compleja la solución
del problema dinámico de la estructura.
Toda esta descripción que antecede será tratada con más profundidad a lo largo de este
trabajo, cuyos conceptos se fundamentan en la dinámica lineal de las estructuras, en la mecánica
de medios continuos y en las técnicas numéricas entre las cuales está el método de los elementos finitos.
Existen diversos enfoques sobre el contenido y el desarrollo de los conceptos que debe
contener un curso de dinámica no lineal de las estructuras y todos ellos son válidos en
consecuencia con los objetivos que se proponen alcanzar. El presente trabajo se ha
decidido trata los conceptos necesarios para completar la formación básica de quien haya
sido formado en la dinámica lineal de las estructuras, en la mecánica de medios continuos y
en el método de los elementos finitos. Es por esta razón que a lo largo del desarrollo de
este trabajo, no se insiste sobre temas que están incluidos en ésta formación básica de las
estructuras y que ya se suponen sabidos por el lector.
A continuación se hace una descripción breve del contenido que se desarrolla en las
páginas de este libro. En el capítulo B2, se hace una introducción a las bases termodinámicas de
la ecuación del movimiento. Este capítulo fundamental muestra los orígenes del problema y lo
enmarca dentro de una formulación estructurada que permite abordar con coherencia
todos los temas restantes. En el capítulo B3 se muestra con detalle los métodos para
resolver la ecuación del movimiento, presentándose los procedimientos implícito y
explícito, como así también las ventajas y desventajas que presentan cada uno de ellos. El
capítulo B4 estudia el concepto de estabilidad de la solución de sistemas conservativos
para distintos métodos de resolución de la ecuación del movimiento. Una vez establecidas
1-2 DINÁMICA NO LINEAL
las bases de estudio de la estabilidad de solución para sistemas lineales, se hace una
aproximación al problema no lineal y se dan criterios para el estudio de la estabilidad
exigiendo la conservación de la energía como requerimiento fundamental. Esto mismo, da
lugar a la “formulación de métodos de solución conservativos”, que actualmente se están
utilizando en la dinámica no lineal. Establecidas las bases fundamentales de la dinámica no
lineal, el capítulo B5 aborda la formulación constitutiva independiente del tiempo, tales
como la plasticidad y el daño, y muestra la forma en que estos comportamientos afectan a
la no-linealidad de la estructura. De manera análoga, el capítulo B6 detalla el
comportamiento constitutivo de los materiales dependientes del tiempo, tales como la
elasticidad retardada y la relajación, donde se introduce en forma natural la no linealidad
por amortiguamiento. Sobre este punto se insiste pues, se considera que es uno de los sitios
de la dinámica no lineal donde se encuentra un gran vacío conceptual.
B 2 Bases Termodinámicas
de la Ecuación del
Movimiento.
B2.1 Introducción.
Este capítulo introduce las bases termodinámicas necesarias para definir el
comportamiento lineal o no lineal de un sólido durante el tiempo que dure el proceso
mecánico. Los conceptos que aquí se sintetizan ayudan a comprender el comportamiento
no lineal de los sólidos y permiten establecer claramente el equilibrio del mismo en cada
instante.
Se considera conveniente hacer una reseña sobre la cinemática del sólido deformable,
estableciendo la notación que se utilizará, como así también las definiciones de la mecánica
de medios continuos que se consideran importantes recordar. De la misma manera se hace
una muy breve descripción de termodinámica para puntualizar los aspectos importantes a
tener en cuenta en la formulación de los modelos constitutivos que representan la no-
linealidad en el comportamiento del sólido. Es muy recomendable recurrir a las referencias
de mecánica de medios continuos y termodinámica,1,2,3 en caso de que se quiera mayor
profundidad y amplitud en los conceptos que aquí se establecen.
1 Malvern, L. (1969). Introduction to the mechanics of continuous medium. Prentice Hall, Englewood Cliffs, NJ.
2 Lubliner, J. (1990). Plasticity theory. MacMillan, New York.
3 Maugin, G. A. (1992). The thermomechanics of plasticity and fracture. Cambridge University Press.
2-2 Bases termodinámicas de la ecuación del movimiento
∂x ACTUAL x ;t
F =
∂x 0
Ωt ⊂ R3
x
∂x
x 0 ;t 0 Fe =
∂x t
X
Ω0 ⊂ R3
REFERENCIAL
y X
xt ; t
P∂x
x F = t Ωt ⊂ R 3
∂x 0
INTERMEDIA
En =
n
(
1 n
U −I ) (B2.8)
para : n = 0 ⇒ E 0 = ln U , Def. natural
En =
n
(
1 n
U −I )
⇒ para : n = 1 ⇒ E1 = U − I (B2.9)
1
para : n = 2 ⇒ E = E 2 = (C − I ) , Def. de Green - St. Venant
2
La deformación euleriana medida en la configuración actualizada se expresa en la forma de
Almansi. Esto es,
e=
1
2
(
I − B −1 ) (B2.10)
ε = ∇Su =
1
2
(
∇ 0 u + ∇ T0 u ) (B2.11)
∂u ∂x
∇0u = = − I = (F − I ) = j (B2.12)
∂x 0 ∂x 0
∂x ∂
F= = (x 0 + u ) = I + ∂u = I + j
∂x 0 ∂x 0 ∂x 0
(B2.13)
1
(
ε = ∇ S u = j + jT
2
)
e
D
∂x
F = x σ ; τ = F ⋅σ
∂x 0
c
ACTUAL
x; t
E
∂x
X
&
E Fe =
∂x t
S
C E
REFERENCIAL D
X
x 0 ;t 0 σ
P ∂x
F = t C
y ∂x 0
INTERMEDIA
x xt ; t
Figura B2.2 – Relación entre variables mecánicas entre las distintas configuraciones.
CONF.
φ( A # ) A# ESPACIAL
Lv ( A # )
CONF. REFERENCIAL
d
dt
[φ( A # ) ] d
[ ]
φ φ( A # )
dt
x0
Sólo es posible hacer una derivada temporal
objetiva en la configuración de referencia
∂
[
L (A # ) = F −T ⋅ F T ⋅ (A # ) ⋅ F ⋅ F −1
v
]
∂ t
(B2.15)
L A
v
( ) # ∂
[ ( )
= F ⋅ F −1 ⋅ A # ⋅ F −T ⋅ F T]
∂ t
donde puede verse en la primera de las anteriores la derivación temporal de un tensor
covariante, en tanto en la segunda se tiene la derivación temporal de un tensor contravariante.
B2.2.5 Velocidad.
El cálculo de la velocidad necesita una derivación temporal objetiva como la definida
anteriormente o bien la clásica derivación material objetiva1,2,3. La velocidad en la configuración
referencial o formulación lagrangeana se define como x& 0 = ∂x 0 ∂t , y su correspondiente
imagen en la configuración actualizada como v = x& . Se define ahora la derivada material
objetiva de una variable cualquiera como
derivada
D ∂φ ∂φ
φ( xi ; t ) material
→ φ= + vi (B2.16)
Dt ∂t ∂xi
Tal que la velocidad del gradiente de deformación resulta entonces,
DINÁMICA NO LINEAL 2-7
∂x ∂ 2x ∂x& ∂v ∂v ∂x
F = ∇0x = ⇒ F& = = = = = L⋅F (B2.17)
∂x 0 ∂t ∂x 0 ∂x 0 ∂x 0 ∂x ∂x 0
de donde el gradiente espacial de velocidades puede también escribirse como:
( )
∂v & −1 r r
L=∇ v= = F⋅F ∈ L V;V (B2.18)
∂x
Se considera en la configuración material la posición de un punto x 0 , su cambio
infinitesimal de posición x 0 + dx 0 y la correspondencia de este movimiento en la
configuración actualizada x → x + dx . Por otro lado, si el movimiento diferencial puede
expresarse según la ecuación (B2.17) como dx = Fdx 0 , entonces la velocidad en la
configuración actualizada resulta luego de aplicar la correspondiente derivación objetiva,
Dx Dt = F& dx 0 = Ldx (ver también ecuaciones (B2.17) y (B2.18)).
El tensor velocidad de deformación se define como,
D=
1
2
( )
L + LT = {L}S (
r r
∈ L V; V ) (B2.19)
Ω=
1
( )
L − LT = {L}A
2
1
(
r r
∈ L V;V ) (B2.20)
Ω = ∇×v
2
En vista de estas dos últimas ecuaciones, también puede escribirse el tensor gradiente de
velocidades, como,
L=D+Ω (B2.21)
La velocidad de deformación lagrangeana se puede expresar a partir de la ecuación (B2.9)
como,
E
∂t ∂t n
( )
& = ∂E = ∂ 1 U n − I = φ(D ) = F T ⋅ (D) ⋅ F (B2.22)
& = 1C
E
2 2 ∂t
( 2
) (
& = 1 ∂ F T ⋅ F = 1 F& T ⋅ F + F T ⋅ F& ) (B2.23)
E
2
( )
& = 1 F T ⋅ LT ⋅ F + F T ⋅ L ⋅ F = 1 F T
2
⋅ (LT )
+ L ⋅ F = FT ⋅ (D) ⋅ F = Lv (e ) (B2.24)
r 1
τ = φS = F ⋅ S ⋅ F T ; σ= F ⋅ S ⋅ FT
J (B2.25)
(
S& = f REF E& ; C T ) σ& = f ACT (e& ; c )
En esta última puede verse dos variaciones temporales de las relaciones constitutivas
que pueden establecerse en forma equivalente en una u otra configuración, es decir en la
configuración referencial, (REF), se relaciona la tensión de Piola-Kirchoff con la
deformación de Green-Lagrange y en la configuración actualizada, (ACT), la tensión de
Cauchy con la deformación de Almansi. Las tres medidas de la tensión coinciden en
pequeñas deformaciones.
Axioma: Existen:
DINÁMICA NO LINEAL 2-9
1. Una cantidad de calor propio Q prop. (calor propio más transferencia global
de calor), regida por leyes físicas bien definidas (Newton, Fourier),
2. Y una cantidad de energía global interna W , función del estado físico del
sólido, que depende de la variable de estado ω que representa densidad de
energía interna o densidad de energía interna por unidad de masa ρ y
volumen V. Tal que se cumple la siguiente relación,
d dW
∫ ρω dV = = W& = Q prop + Pd (B2.26)
dt V dt
Relación entre:
La potencia Mecánica La potencia No-
deformativa: Mecánica o cantidad
Y de Calor Propio:
Pd
Q prop
con
El cambio de energía interna que
experimenta el cuerpo:
W&
Pint = ∫ t ⋅ v dS + ∫ ρ b ⋅ v dV = ∫ t i v i dS + ∫ ρ bi vi dV (B2.27)
S V S V
En el caso que la aceleración sea nula, o el problema sea independiente del tiempo, se
obtiene de la expresión (B2.29), la clásica ecuación de equilibrio estático de Cauchy,
∂v i ∂σ ij
=0 ⇒ = −ρbi
1∂t23 ∂x j (B2.30)
14243
problema equilibrio de Cauchy
cuasi−estático
Donde K& es la potencia cinética y Pd la potencia deformativa, tal que esta última resulta
entonces,
Pint = K& + Pd ⇒ Pd = Pint − K& (B2.32)
Observando ésta última expresión, resultan tres casos típicos de la mecánica,
∂v i
• Problemas de comportamiento de sólido rígido: = 0 ⇒ Pint ≡ K&
∂x j
• Problemas de comportamiento cuasi estático: K& = 0 ⇒ Pint = P 0
• Problemas de vibración libre: Pint = 0 ⇒ K& = − P0
siendo q el flujo de calor que sale por la frontera, q i el vector de flujo de calor, ni la
normal saliente a la superficie que envuelve la frontera y r la fuente de calor por unidad de
masa.
Sustituyendo Pd y Q prop en la ecuación (B2.26), resulta la ecuación de conservación en
problemas cuasi-estáticos
∂vi
∫ ρ ω& dV = ∫ ρ r dV − ∫ qi ni dS + ∫ σ ij ∂x j dV (B2.34)
V 14442S444
V 3 V
142
4 434
Q prop Pd
∂vi
ρω
& = ρ r − div(q i ) + σ ij (B2.35)
∂x j
∂vi
σ ij = σ ij Lij = σ ij Dij + σ ij Ω ij (B2.36)
∂x j 123
0
o bien
2-12 Bases termodinámicas de la ecuación del movimiento
& = ρ r dV − q i n dS
ℑ ∫ θ ∫θ i (B2.42)
V S
d &
ℑ ≥ ℑ
dt
{ {
INCREMENTO INCREMENTO
DE ENTROPÍA DE ENTROPÍA
EN EL SISTEMA INTRODUCIDA
σ D
Ξ m = θ η& − ω& + ij ij ≥0 → Disipación Mecánica
ρ (B2.48)
Ξ θ = − q i ∇θ ≥ 0 → Disipación Térmica
( ) { }
def
Ψ = Ψ eij ; θ; p i = ω − θη con p i = FijP ; α i ; d i (B2.49)
Ξm = ρ − Ψ ( )
& − θ& η + σ D ≥ 0
ij ij (B2.52)
Sustituyendo en esta una forma general de la variación temporal de la energía libre,
donde eij y θ representan las variables libres mecánicas y térmicas y pi las variables
internas del proceso no lineal. Considerando esto, la disipación resulta entonces,
∂Ψ ∂Ψ ∂Ψ
Ξ m = σ ij − ρ Dij − ρ + η θ& − ρ p& i ≥ 0 (B2.54)
∂e ∂θ ∂p
ij i
Dij y θ& son las variaciones temporales de las variables libres, por lo tanto, para garantizar el
cumplimiento de la desigualdad de Clausius-Duhem, sus multiplicadores deben ser
idénticamente nulos. Es decir, de aquí resultan las ecuaciones constitutivas y la disipación
mecánica,
σ ij − ρ ∂Ψ = 0 ⇒ σ ij = ρ ∂Ψ
∂eij ∂eij
Ec. Constitutiva (B2.55)
∂Ψ ∂ Ψ
+ η = 0 ⇒ η = −
∂θ ∂θ
2-14 Bases termodinámicas de la ecuación del movimiento
∂Ψ
Ξm = ρ p& i ≥ 0 Disipación (B2.56)
∂p i
∂x
σ : D = tr (σ : D ) ; D = F −T ⋅ E
& ⋅ F −1 ; F= (B2.57)
∂x 0
donde F representa el gradiente de deformación que es un tensor bipuntual que relaciona
un punto de una configuración en referencia x 0 , con el mismo punto en una configuración
actualizada x , (esto se verá más adelante y puede también consultarse en las
referencias1,2,3). Operando algebraicamente se llega a,
(
σ : D = tr σ ⋅ F −T ⋅ E ) ( & = tr 1 S : E
& ⋅ F −1 = tr F −1 ⋅ σ ⋅ F −T ⋅ E )
& (B2.58)
J
donde J , determinante de la matriz Jacobiana y S el tensor de tensiones en la
configuración de referencia o también llamado tensor de Piola-Kirchoff.
1 & dV
σ:D= S:E ; J = det(F) = (B2.59)
J dV0
Sustituyendo esta última expresión en la forma Euleriana de la disipación, se obtiene la
forma Lagrangeana, a veces mas apropiada para resolver algunos problemas estructurales.
Esto es,
Ξ m = Jρ − Ψ (
& − θ& η + S E& ≥ 0
ij ij ) (B2.60)
Donde la energía libre de Helmholtz Ψ se define ahora como,
( ) { }
def
Ψ = Ψ E ij ; θ; Pi = ω − θη con Pi = FijP ; α i ; d i (B2.61)
∂Ψ
Ξ m = S ij − ρ 0 E& ij − ρ 0 ∂Ψ + η θ& − ρ 0 ∂Ψ P&i ≥ 0 (B2.63)
∂E ij ∂θ ∂Pi
E& ij y θ& son las variaciones temporales de las variables libres, por lo tanto, para garantizar el
cumplimiento de la desigualdad de Clausius-Duhem, sus multiplicadores deben ser
idénticamente nulos, resultando de aquí las ecuaciones constitutivas y la disipación en la
configuración lagrangeana,
S − ρ ∂Ψ = 0 ⇒ S = ρ ∂Ψ
ij 0
∂E ij ij 0
∂E ij
Ec. Constitutiva (B2.64)
∂Ψ ∂Ψ
+ η = 0 ⇒η= −
∂θ ∂θ
∂Ψ &
Ξ m = ρ0 Pi ≥ 0 Disipación (B2.65)
∂Pi
∂u& i (B2.67)
∫ ρω& dV − ∫ ρ r dV − ∫ qi ni dS = ∫ σ ij Dij dV = ∫ t i u& i dS + ∫ ρbi u& i dV − ∫ ρ u& i
∂t
dV
V V4
14444 S
42444 3
444 V
14243 S 4442
1 V4443 V
144244 3
Pot. Mecánica Pot. Deformativa Pot. Introducida Pot. Cinética
u j ( x, y , z ) = N jk ( x, y, z ) U k ⇒ eij = ∇ iS u j = ∇ iS N jk U k (B2.70)
Ωe Ωe Ωe Ωe Ωe
5 Zienkiewicz, O. and Taylor, R. (1989). The finite element method. McGraw-Hill, Vol I y II.
DINÁMICA NO LINEAL 2-17
∫ σ ij ∇ iS N jk dV U& k = ∫ t i N ik dS + ∫ ρbi N ik dV − ∫ ρ N ki N ij U&& j dV U& k (B2.71)
Ωe Ωe
V e e S e Ve Ve
e
Ω Ω
Pero esta ecuación se cumple para cualquier velocidad U& k , por lo tanto la igualdad
Ωe
establecida en la ecuación (B2.71) es independiente de esta velocidad, obteniéndose de aquí
la siguiente ecuación de equilibrio dinámico para el sólido discreto
M kj e
Ω
64447444 8
siendo f kint , f kmas y f kext los conjuntos ordenados, en forma de matrices columna,
Ωe Ωe Ωe
de las fuerzas interna, másica y externa que se desarrollan en cada punto del sistema
discreto que aproxima el continuo, U&& j e la aceleración en dichos puntos, M kj e la masa
Ω Ω
elemental y Bijk Ωe
= ∇ iS N jk
Ωe
el tensor de compatibilidad de deformaciones o gradiente
simétrico de la función de forma.
La ecuación (B2.72) representa la ecuación elemental de equilibrio dinámico en la
configuración actualizada, que expresada en la configuración de referencia adquiere la siguiente
forma,
&&
∫ ij i jk 0 ∫ i ik 0 ∫ 0 i ik 0 ∫ 0 ki ij 0 U j
S
S ∇ N dV = t N dS + ρ b N dV − ρ N N dV
Ωe
Ve e Se Ve e V e e 0
0 Ω 0 0 Ω 0 Ω
0 0 0
(B2.73)
tensor constitutivo C ijkl (ver Tabla B2.1), y por ello establece una relación no lineal
entre tensiones y deformaciones (ver ecuación (B2.25)). Además, estos cambios de
configuración son producidos por grandes movimientos, traslaciones y rotaciones,
que también producen cambios en el sistema de referencia local en los puntos del
sólido, afectando por ello al tensor de compatibilidad de deformaciones B ijk .
- No linealidad por grandes desplazamientos, que a diferencia de las grandes
deformaciones, sólo afecta al tensor de compatibilidad de deformaciones B ijk ,
porque en este caso sólo ocurren cambios en el sistema de referencia local de los
puntos del sólido como consecuencia de grandes movimientos.
Estas posibles no linealidades que pueden ocurrir todas juntas o por separado, se
resumen en la forma que se muestra en el siguiente cuadro descriptivo,
Ωe
[
0 = Α f kmas + f kint − f kext ]
Ωe
= ∆f k Ω (B2.74)
DINÁMICA NO LINEAL 2-19
fuerzas másicas f kmas y las exteriores f kext . Este desequilibrio, en un cierto instante de
Ω Ω
tiempo “t” del proceso dinámico, puede eliminarse mediante la linealización de esta fuerza
residual ∆f k Ω (B2.74), en la vecindad del estado de equilibrio actual (i+1). Para ello es
necesario forzar el equilibrio en el estado actual (i+1) y expresar dicha condición mediante
una expansión en serie de Taylor truncada en su primera variación,
i ∂ (∆f ) t i +1
0 =Α
Ωe
i +1
[∆f ]
k
t
Ωe
≅ Α ∆f k
Ωe
i
[ ] t
Ωe
+ Α
Ω ∂U r e
e
k
⋅ ∆U r [ ]
Ωe
t
Ω
(B2.75)
i t
∂U&& j ∂f kint ∂f kint ∂U& j ∂f kext i +1
0=
i +1
[∆f ] ≅ [∆f ]
k Ω
t i
k Ω
t
+ Α M kj
Ωe
+ + −
∂U r ∂U r ∂U& j ∂U r ∂U r e
⋅ ∆U r [ ] t
Ωe
Ω
donde la aceleración y la velocidad deben expresarse mediante una aproximación lineal en
diferencias finitas, ver más adelante, en el apartado (B3.3.21), el método de Newmark como
un ejemplo de esta aproximación. Sustituyendo en esta ecuación las fuerzas internas y
másicas expresadas en la ecuación (B2.72) (En forma análoga puede procederse con la
ecuación (B2.73) ), se tiene,
i t
0 =Α M kjU&& j + ∫ σ ij ∇ iS N jk dV − f kext +
Ωe Ve
e
Ω
i t
∂U&& j & ext
+Α ∫ ρ N ki N ij dV
∂ σ ∇ S N dV + ∂ σ ∇ S N dV ∂U m − ∂f k ⋅
∫e ∫e
+
Ωe
Ve ∂U r ∂U r
ij i jk
∂U& m ij i jk
∂U r ∂U r
V V Ωe
⋅Α
i +1
[∆U r ]tΩe
Ωe
i t
0 =Α M kj U&& j + ∫ σ ij ∇ iS N jk dV − f kext +
Ωe Ve
e
Ω
i t
∂U&& j ∂σ ij ∂e ∂σ ij ∂D ∂U& ∂f ext
+Α ∫ ρ N ki N ij dV + ∫ st
∇ iS N jk dV + ∫ st
∇ iS N jk dV m
− k ⋅
Ωe e ∂U r e ∂e st ∂U r e ∂D st ∂U& m ∂U r ∂U r e
V V V Ω
⋅Α
i +1
[∆U r ]tΩe
Ωe
Tal que particularizando esta ecuación de equilibrio dinámico para un material cuya ley
constitutiva visco elasto-plástica es del tipo σ ij = ρ(∂Ψ (eij , pi ) / ∂ eij ) = C ijkl : ekle + ξ ijkl : Dkl
para una relación cinemática del tipo eij = ∇ iS u j = ∇ iS N jk U k , y Dij = ∇ iS u& j
= ∇ iS N jk U& k , resulta,
2-20 Bases termodinámicas de la ecuación del movimiento
i
i
t &&
0 =Α M kj U&& j + ∫ σ ij ∇ iS N jk dV − f kext + Α ρ N N dV ∂U j +
Ωe e Ωe ∫e ki ij
∂U r
Ve Ω V
t
∂U&
e
( ) e
(
+ ∫ ∇ Ss N tr C Tijst ∇ iS N jk dV + ∫ ∇ Ss N tr ξ ) ( ) T
ijst (∇ S
i N jk ) dV m ∂f ext
− k
∂U r ∂U r
⋅
(B2.76)
V V Ω e
Αe i +1 [∆U r ]tΩe
Ω
[ ]
0 = i ∆f k
t
Ω
[ ]
+ i J krT
t
Ω ⋅
i +1
[∆U r ]Ωt
donde ξ Tijst = ξ T es el tensor de viscosidad tangente y J krT = J T es el operador jacobiano
tangente. Esta ecuación puede también presentarse en la siguiente forma matricial, donde
se detallan los operadores que contribuyen a la definición del jacobiano,
i t
∂U&& & ∂f ext
T ∂U
0=
i +1
[∆f ]t
Ω≅
i
[∆f ] t
Ω + M T
+ K +D − ⋅
i +1
[∆U ]tΩ
∂U ∂U ∂U
Ω
(B2.77)
144444424444443
i Jt
Ω
Αe ∫V e ρ N : N dV es la matriz de masa,
Ω
[D ]
T
Ω
= Α∫
Ωe V
e (∇ N ): ξ : (∇ N ) dV
S T S
es la matriz de
desequilibrada en el sólido
i +1
[ ] t
∆f k Ω se elimina siguiendo una resolución por Newton-
5
Raphson hasta que este residuo resulte despreciable, situación que se conoce como
convergencia del proceso linealizado hacia la solución exacta (ver Figura B2.4).
En la Figura B2.4 se describe el equilibrio espacial, dejando el tratamiento de la
convergencia en el tiempo para ser tratado al estudiar los métodos de resolución en el
tiempo de la ecuación de equilibrio dinámico.
1 ∂u ∂u
T
ε=
1
2
(F F T
)
− I = ∇Su = +
2 ∂x ∂x
(B2.79)
i +1
[∆f k ]tΩ = λt ⋅ [ f kext ]Ω Iteración: “i=0”
[∆f ]
i t
[∆U r ]tΩ = [∆U r ]Ω + J krT ⋅
i +1 i t −1 i +1 t
Ω k Ω
Nueva iteración:
i +1
[U r ]tΩ = [U r ]tΩ− ∆t + i +1 [∆U r ]tΩ
“i+1” i +1
[x r ]tΩ = [X r ]0Ω + i +1 [U r ]tΩ
Cálculo de la deformación
i +1 i +1 t
P.G . t ∇ S N U
e = P.G .
ij Ωe i jk k Ωe
Ecuación constitutiva
i +1 P .G . t i +1 P .G . t
e → σ
ij Ωe ij Ωe
Fuerza residual
i +1
[∆f k ]tΩ = Αe i +1 [∆f k ]tΩe ≅ Αe i +1 [ f kmas + f kint − f kext ]
t
Ωe
Ω Ω
Verificación de equilibrio
y convergencia
NO
¿ i+1[∆fk ]tΩ →0 ?
Nuevo incremento:
t+∆t ⇒“n+1”
FIN
Tiempo de retardo t ξ t
=r
ε C
σ
σ0 ε0
t t
Modelo de Kelvin Modelo de Maxwell
σ e (t ) = C ε e (t )
σp = Cεp
σ
σlim σlim
(1 − d )C 0
e
σ = Cε
C C C0
εp εe ε ε ε
ε
σ = C ε e = C (ε − ε p ) σ = (1 − d )C ε = C ε − d Cε
dispositivo de fricción
Parte degradada del muelle
σlim σ σ
σ C σ
εe (1 − d )C dC
εp
ε ε
Problema de Viscoplasticidad
σ
dispositivo de fricción
σ p = C ε vp σlim
C
σ σ
ξ
σlim
σ = C εe εe ε vp
C C ε
σ = C ε e = C (ε − ε vp )
εp εe ε
ε
B3.1 Introducción.
En este capítulo se estudia la resolución en el tiempo de la ecuación del movimiento en
su forma semi-discreta (ver ecuación de equilibrio, apartado B2.5). A continuación se
considera el ensamblaje de la ecuación B2.72 (o la B2.73 si se hace el equilibrio en la
configuración referencial), que permite definir el equilibrio en todo el sólido en el instante
t + ∆t ,
i t + ∆t i t + ∆t i t + ∆t
t + ∆t
Αe ∫ σ ij ∇ i N jk dV
S
=Α ∫ t i N ik dS + ∫ ρ bi N ik dV −Α ∫ ρ N ki N ij dV U&& j
Ωe
Ω Ωe Ωe
V e e S e Ve e V e e
Ω Ω Ω
t+∆t (x, y, z)
=cte
(x, i=1,n
y, z)
Resolución
Resolución temporal
Problema:
iterativa del
Discretizado en el espacio (MEF)
del problema
y
problema dinámico.
Linealizado en sus variables
no-lineal mecánicas.
linealizado.
U(t + ∆t ) = U (t ) + ∆t U & (t + α ∆t )
& (B3.4)
& (t ) + α U
U(t + ∆t ) = (1 − α ) U & (t + ∆ t )
U(t+∆t)
U(t+α ∆t) U(t+∆t)
U(t+α ∆t)
U(t) α=1
t t+α ∆t t+∆t t
Debido a la robustez que ofrece en la solución los métodos implícitos, se orientará el resto de
esta presentación a presentar sus particularidades. En todos los casos se admitirá un problema
espacial aproximado mediante una forma discreta tal como lo formula el método de los
elementos finitos y sólo se estudiará la forma de resolver el problema temporal, dejando la
solución del problema espacial para un tratado más orientado a problemas independientes del
tiempo.
Se admite que los coeficientes k U y k V son no nulos, que las aceleraciones U && t + ∆t y
velocidades U& t + ∆t son funciones del desplazamiento U t + ∆t en el tiempo (t + ∆t ) y que la fuerza
residual en dicho instante i [∆f ]tΩ+ ∆t , linealizada para cada iteración "i" , está definida en todo el
dominio Ω (ecuación B2.77), por la expresión
0=
i +1
[∆f ]Ωt + ∆t ≅ i [∆f ]Ωt + ∆t + i J tΩ+ ∆t ⋅ i +1 [∆U ]tΩ+ ∆t
(B3.6)
0= ∆f i +1 t + ∆t
≅ ∆fi t + ∆t i
+ J t + ∆t
⋅ [ i +1
U t + ∆t i
− U t + ∆t
]
En el resto de este capítulo y para simplificar la notación, siempre que se haga referencia al
equilibrio global, se prescindirá del subíndice Ω en las formulaciones. Este equilibrio global bajo
comportamiento no lineal puede alcanzarse satisfactoriamente mediante el procedimiento
iterativo de Newton-Raphson, el cual permite aproximar la solución en la vecindad del punto
“(i+1)”, mediante la linealización descrita en la ecuación (B3.6). En esta solución linealizada, el
operador jacobiano tangente tiene la siguiente forma,
i t + ∆t
i t + ∆t
∂U & ∂f ext
&& ∂f int ∂f int ∂U
i
J t + ∆t
=J U( i t + ∆t
) ∂ ∆f
= = M + +
& ∂U
− (B3.7)
∂U ∂U ∂U ∂U ∂U
T
∂f int
Tal que se entiende por operador tangente de rigidez a la relación K = , operador
∂U
∂f int ∂f ext
tangente de amortiguación a la relación D T = y a la relación como la influencia del
∂U& ∂U
cambio de la posición de las fuerzas externas con los sucesivos cambios de configuración. Este
último término puede considerarse nulo en pequeños desplazamientos, pues en este caso los
cambios de posición de las cargas son casi despreciables frente al tamaño de la pieza.
Sustituyendo estos operadores tangentes en la ecuación (B3.7), la matriz jacobiana toma la
siguiente forma,
i t + ∆t
i t + ∆t
∂U & ∂f ext
( )
i t + ∆t i t + ∆t ∂ ∆f && T T ∂U
J =J U = = M + K +D − (B3.8)
∂U ∂U ∂U ∂U
En problemas dinámicos elásticos lineales se simplifica aun más, por que se transforma en el
siguiente operador constante,
&&
∂U &
∂U
i
J t + ∆t ≡ J 0 = M + K 0 + D0 (B3.9)
∂U ∂U
y en el caso de problemas cuasi estáticos, donde && = U
U & =0 y f ext = cte , el operador
jacobiano tiende a la clásica matriz de rigidez.
t + ∆t
[ ]
i
i
( ∂ ∆f
J t + ∆t = J i U t + ∆t = ) i
= KT
t + ∆t
(B3.10)
∂U
3-6 RESOLUCIÓN DE LA ECUACIÓN DEL MOVIMIENTO
U
Término de 1er orden 2do orden
(
&& t + ∆t , U
f3 U && t )
&t &t
U t + ∆t U
U
t t+∆t τ
Figura B3.3 – Aproximación temporal del desplazamiento.
La diferencia fundamental con el clásico método de diferencia centrada1 se basa en que la
velocidad y el desplazamiento no resultan de la aceleración instantánea (una derivada), sino de la
integral de la aceleración que hace la solución más estable. Así, la velocidad y el desplazamiento
se definen siguiendo el procedimiento que a continuación se resume,
DINÁMICA NO LINEAL 3-7
• Cálculo de la velocidad:
& =U
dU && (τ) dτ
&
&& (τ) = d U(τ)
U
dτ t + ∆t t + ∆t
τ
Área del diagrama
t dτ t+∆t
• Cálculo del desplazamiento:
& (τ) dτ
dU = U
& ( τ) = d U ( τ)
U
dτ
t + ∆t t + ∆t
( t + ∆t ) − τ ∫ dU = ∫ U& (τ) dτ
t t
t + ∆t
(B3.13)
&& (τ)
U
U t + ∆t = U t + ∫ U& (τ) dτ
t
Momento de primer orden.
Figura (B3.3)
τ
t dτ t+∆t
t + ∆t
Sustituyendo la velocidad U& t + ∆t = U& t + ∫ U&& (τ) dτ , previamente obtenida, en ésta última
t
expresión, resulta la expresión del desplazamiento en función de la aceleración,
t + ∆t
τ
U t + ∆t = U t + ∫ &
U ( τ) + ∫ && (τ) dτ dτ
U
t 0
∆t ∆t τ
(B3.14)
= U + ∫U ∫ [ ∫ U&& (τ) ] dτ
t + ∆t & (τ) dτ +
t
U dτ
0 0 0
f(τ)
&& (τ)
U &&
∆U
1
t τ t+∆t τ t t+∆t τ
&
U = &
U t
+ ∫ &
U& t
dτ + ∫ f ( τ) U (
&& t + ∆t − U
&& t dτ )
0 0
(B3.17)
[ )] dτ
∆t
t + ∆t
U = Ut + U& t ∆t + [(t + ∆t ) − τ ] ⋅ U
∫ (
&& t + ∆t − U
&& t + f (τ) U && t
0
f(τ)
1
∆t
g = ∫ f (τ) dτ = γ ∆t
0
t t+∆t τ
∫ f (τ) dτ = k U ∆t = γ ∆t
0
(B3.19)
∆t τ
∫ ∫
2 2
f ( τ) dτ dτ = k V ∆t = β ∆t
0 0
i && t + ∆t + i f int
∆f t + ∆t = M i U [ ] t + ∆t
−[f ] =
i +1 ext t + ∆t
1 i t + ∆t
= M 2
∆ U
( & ∆t ) − 1 − 1
− U i t
2β
&& t +
U i
β ∆t (B3.24)
+∫
V
i
(
σ t + ∆t ∇ S N dV − ) i +1
[f ] ext t + ∆t
Puede probarse que la solución en el tiempo que ofrece este método es incondicionalmente
estable para valores de γ ≥ 1 2 y β ≥ 1 4 (0.5 + γ )2 (ver cap-B-4).
El esquema de integración en el tiempo se realiza siguiendo los siguientes pasos, luego de
imponer la condición de fuerza residual nula ( i +1∆f t + ∆t = 0 , ecuación (B3.23)),
1 i +1 t + ∆t i +1 & t
∆f t + ∆t = M
i +1
(
∆ U
1
− U ∆t −
− 1 ) i +1 && t
U +
2β
2
β ∆t
+∫
V
i +1
(
σ t + ∆t ∇ S N dV − ) i +1
[f ]ext t + ∆t
i +1
5. Si ∆f t + ∆t > TOL , entonces “IR a 2”; caso contrario “IR a 1” e incrementar el tiempo ( t + ∆t ).
∑ δ k [i +1∆f t + ∆t ]k =
m
i +1 1
∆f ∗ =
δ0 k =1
[f ] [ ] (B3.26)
m
1
∑ δ k M =0
i +1 && ( t + ∆t ) − k ⋅∆t i +1 int ( t + ∆t ) − k ⋅∆t ( t + ∆t ) − k ⋅∆t
= U + − f ext
δ0 k =1
y de acuerdo a la ecuación de la aceleración (B3.25), puede reescribirse las fuerzas residuales sólo
en función de los desplazamientos e implicitamente de las velocidades. Esto es,
γ
[f ] [ ] = 0
m
1
U (t + ∆t ) − k ⋅∆t + δ k
i +1 int ( t + ∆t ) − k ⋅∆t ( t + ∆t ) − k ⋅∆t
i +1
∆f ∗ =
δ0
∑ ∆tk2 M i +1
− f ext
(B3.27)
k =1
& t + ∆t = α 0 U t + ∆t + 1 ∑ α U (t + ∆t ) − k∆t − ∆t β U ( )
m
U & (t + ∆t ) − k∆t = 0 ⇒
k k
β 0 ∆t β 0 ∆t k =1
(B3.29)
& t + ∆t
∂U α
⇒ = 0
∂U β 0 ∆t
i +1
∆f ∗ = − i J t + ∆t ⋅ ∆ i +1 U t + ∆t (B3.31)
El esquema de integración en el tiempo se realiza siguiendo los siguientes pasos, luego de
forzar la condición de residuo nulo ( i +1∆f t + ∆t = 0 ), en la ecuación (B3.26),
i +1
∆f ∗ =
1
δ0
m
γ
∑ ∆tk2 M
k =1
i +1
U (t + ∆t ) − k + δ k
i +1
[∫V
i +1
( ) ]
σ t + ∆t ∇ S N dV
( t + ∆t ) − k
−
i +1
[f ]
ext ( t + ∆t ) − k
i +1
7. Si ∆f ∗ > TOL , entonces “IR a 3”; caso contrario “IR a 1” e incrementar el tiempo ( t + ∆t ).
2Press W., Flannery B., Teukolsky S., Vetterling W. (1987). Numerical Recipes – The art of scientific computing.
Cambridge University Press.
DINÁMICA NO LINEAL 3-13
∂U ∂U ∂U
Ω3
1 4444442444444
i J t + ∆t
Ω
f
2
∆f t + ∆t 3
∆f t + ∆t
f t + ∆t
3
J tΩ+ ∆t t+∆t
3
2
J tΩ+ ∆t
1
∆f t + ∆t 2
1
J tΩ+ ∆t
ft t 1
1
[∆U]t + ∆t 2
[∆U]t +∆t
1
[U ]t +∆t 2
[U]t +∆t
Ut U t + ∆t U
Figura B3.6 – Método de Newton-Raphson.
f t =3 3er incremento
1
J tΩ+ ∆t
f t =2 2nd incremento
1
J tΩ+ ∆t
1er incremento
f t =1
1
J tΩ+ ∆t
U t =1 U t =2 U t =3
f t =3 3er incremento
3
J tΩ+ ∆t
f t =2
2nd incremento
2
J tΩ+ ∆t
f t =1 1er incremento
1
J tΩ+ ∆t
U t =1 U t =2 U t =3
Sin acelerador:
i +1
[U ]tΩ+∆t = i [U ]tΩ+∆t + i+1 [∆U ]tΩ+∆t (B3.37)
Con acelerador:
i +1
[U ]tΩ+ ∆t = i [U ]tΩ+ ∆t + A ⋅i +1 [∆U ]tΩ+ ∆t (B3.38)
donde A representa la “matriz de aceleración” en el instante (t). Entre las distintas formas
de definir esta matriz de aceleración, está el denominado acelerador de Aitken2, que se
presenta a continuación.
3-16 RESOLUCIÓN DE LA ECUACIÓN DEL MOVIMIENTO
[ J ] ⋅ [ [∆U]
t + ∆t ∗
Ω
i t + ∆t
Ω −
i +1
[∆U]tΩ+ ∆t ] = [ J Ωt + ∆t ]⋅i [∆U]tΩ+ ∆t
y de esta última se define la matriz de aceleración o extrapolación, como
[
A = J tΩ+ ∆t ] ⋅[J ]
∗ t + ∆t −1
Ω ⇒ A⋅ [ [∆U]
i t + ∆t
Ω −
i +1
[∆U]tΩ+ ∆t ]= i [∆U]tΩ+ ∆t (B3.42)
sustituyendo esta matriz de aceleración en la ecuación del cambio de la fuerza residual entre
dos iteraciones sucesivas (ecuación (B3.40)), se tiene la siguiente ecuación de equilibrio
linealizada
i
[∆f ]Ωt +∆t [
= − J tΩ+ ∆t ⋅ ] ∗ i +1
[∆U ]tΩ+∆t = − A ⋅ [ J Ωt +∆t ]⋅i+1 [∆U ]tΩ+ ∆t (B3.43)
De esta forma siempre se resuelve el sistema con un operador secante, cuya expresión se
resume en la siguiente ecuación,
[J ]= A ⋅ [J ]
t + ∆t
Ω
−1
n
t + ∆t ∗
Ω
(B3.44)
Las distintas formas de obtener el operador de aceleración A , da lugar a distintos
algoritmos.
aproximación inversa del operador jacobiano a partir de la igualdad (B3.40). De esta forma,
tras varios pasos de manipulación matemática, se llega a la siguiente fórmula de
actualización de la inversa del jacobiabo,
[ i
J tΩ+ ∆t ]
-1
=i AT [ i −1
J tΩ+ ∆t ]⋅A
-1 i
; i
[
A = I + i V⋅ i W T ] (B3.45)
Siendo,
i
V = 1 −
[∆U ]tΩ+∆t ⋅ [i−1 [∆f ]Ωt +∆t − i [∆f ]Ωt +∆t ]⋅i−1 [∆f ]t +∆t − i [∆f ] t +∆t
i −1
i −1
[∆U ]tΩ+∆t ⋅i−1 [∆f ]Ωt +∆t
Ω Ω
(B3.46)
i
W=
i −1
[∆U]tΩ+∆t
i −1
[∆U ]tΩ+∆t [i−1 [∆f ]Ωt +∆t − i [∆f ]Ωt +∆t ]
Entre las ventajas de este método de Newton secante se encuentra el ahorro de cálculo, pues
sólo debe invertirse la matriz de rigidez inicial una sola vez, además, la forma de obtener la matriz
jacobiana mediante la ecuación (B3.45), preserva la simetría de la matriz original. Por el contrario,
la naturaleza de esta actualización no preserva la distribución de los elementos de la matriz o
topología (forma “sparsa”). Por esta razón es conveniente recuperar la topología original de la
matriz jacobiana 1 J tΩ+ ∆t , en cada iteración y reformular la ecuación (B3.45) en la siguiente forma,
[ ] ( ) i− j
[∆U ]tΩ+ ∆t = ∏ ( i − j A T ) ⋅ i b
i
J tΩ+ ∆t ⋅ ∏ A ⋅ [∆f ]Ω
i 1 j j t + ∆t i
b= ⇒ (B3.47)
j =2 j =0
i
f t + ∆t
i
[∆f ]Ωt + ∆t
i +1
i
i −1
[∆f ]Ωt + ∆t
[J ] t + ∆t
Ω
i −1
[∆f ]Ωt + ∆t − i [∆f ]Ωt + ∆t
i −1 t + ∆t
f
i −1
i
[∆U]tΩ+∆t
i +1
[∆U]tΩ+∆t
i −1
[U]tΩ+∆t i
[U]tΩ+ ∆t i +1
[U]tΩ+∆t [U]Ω
Figura B3.9 – Cuasi-Newton – Algoritmo de Aitken.
b ⋅ bT
[ J ] =[ J ]
i t + ∆t −1 1 t + ∆t −1
+
Ω Ω
bT [ [∆f ]
i t + ∆t i −1
Ω − [∆f ]Ωt +∆t ]
(B3.48)
siendo,
b = [∆U ]Ω + [∆U ]Ω − [∆U]tΩ+∆t
t i t + ∆t i −1
(B3.49)
tal que los incrementos de desplazamiento se mencionan en esta ecuación se obtienen de la
siguiente forma,
i
[∆U]tΩ+∆t = −[ i J tΩ+∆t ]−1 ⋅i [∆f ]Ωt +∆t ,
i −1
[∆U]tΩ+∆t = −[ i J tΩ+∆t ]−1 ⋅i−1 [∆f ]Ωt +∆t (B3.50)
El método consiste en obtener i [∆U ]tΩ+∆t por cualquier procedimiento antes citado. Este
incremento de desplazamiento no será considerado como incremento en sí mismo, sino
como dirección de búsqueda. Dentro de esta nueva dirección se búsqueda el mínimo
direccional P , que no es el mínimo absoluto, se obtiene minimizando µ(s ) = Π[U(s )] , de
donde resulta la posición del mínimo s* . Al igual que en el acelerador de Aitken (ecuación
(B3.38)), el incremento de desplazamiento se escribe como,
i +1
[U]tΩ+∆t = i [U ]tΩ+∆t + i [s ∗ ]t +∆t ⋅i+1 [∆U]tΩ+∆t (B3.52)
Π ( [U] ) Ω
Π ( [U] ) Ω
P
Λ
[U ]Ω
Π[U(s )] = µ(s )
µ (s ) U [s ] dirección de la búsqueda
P s =1
s=0
s = s*
i −1 i j
U [s ]
Figura B3.11 – Plano Λ - Plano de búsqueda (ver Figura (B3.10)).
Teniendo ahora en cuenta las dos ecuaciones anteriores, resulta la siguiente ecuación en
función de la coordenada (s),
(
∂Π [U ]Ω + [s ]
i t + ∆t i t + ∆t
⋅
i +1
[∆U]tΩ+ ∆t )= i+1 ∂Π ∂U t + ∆t = i+1 [∆f ]t + ∆t ⋅i+1 [∆U ]t +∆t = G (s ) = 0 (B3.55)
∂s ∂U ∂s Ω Ω
De esta manera puede verse que se trata de solucionar una ecuación G (s ) , no lineal en ( s ),
cuya solución es s ∗ (ver Figura B3.12).
3-20 RESOLUCIÓN DE LA ECUACIÓN DEL MOVIMIENTO
G (s )
s*
s
[
i +1 ∆f (U
&& , U
Ω]
& , U, λ) t + ∆t = M i +1
[U&& ]
t + ∆t
Ω +
i +1
[f int & , U)
(U Ω
]
t + ∆t
−
i +1
[λ f ]
ext t + ∆t
Ω
=0
(B3.56)
i +1
[c(U, λ)]Ω = 0t + ∆t
( [U]+ ) ( )
i i
∆f
i i +1
[δU ] , [λ]+
i i +1
[δλ] = 0 = ∆f [U ], [λ] + ∂∆f ⋅i +1 [δU ]tΩ+ ∆t + ∂∆f ⋅i +1 [δλ] (B3.57)
i i
∂U ∂λ
Pero en esta última ecuación puede identificarse las siguientes relaciones (ver ecuación
(B3.34)),
i i
∂∆f i ∂∆f ext
= J ; = −f (B3.58)
∂U ∂λ
Que sustituidas en la ecuación (B3.57), se rescribe el sistema de ecuaciones de equilibrio
con restricciones, como
(
0 = ∆f i U , i λ + i J ⋅ ) i +1
[δU ] − f ext ⋅i +1 [δλ]
i +1 (B3.59)
0 = [c(U, λ )]
4 Crisfield M. A. (1983). An arc-length method including line searches and accelerations. Int. Journal for
numerical method in engineering. Vol. 19, pp. 1269, 1289.
5 NOTA: Para simplificar la presentación, se omitirá a continuación el superíndice (t+∆t) en todos los
desarrollos.
DINÁMICA NO LINEAL 3-21
i
J⋅ i
+1
[δU] = −∆f ( i U , i λ ) + f ext ⋅i +1 [δλ]
i
J⋅i
+1
[δU] = − M (U
i
[U&& ]+ [f
& , U) − λ f ext + i J⋅i
i
int
] [ i
] +1
[U]TOT ⋅ i +1 [δλ] (B3.60)
144444424444443
i +1
i
J ⋅ δU ˆ
Donde i +1
[U]TOT es el desplazamiento total obtenido con el último, o máximo, valor de
i +1
fuerza exterior, δÛ es la solución del sistema de ecuaciones sin corrección y δ i +1 λ es
cambio del factor de aplicación de la carga.
Una vez definidas las ecuaciones generales del método, es necesario precisar sobre la forma de
la ecuación de restricción y para ello se introduce el siguiente apartado.
C1 =
i +1
[U]TTOT ⋅i+1 [U]TOT
T
C 2 = 2 [∆U ]+ δU ˆ ⋅i +1 [U ]
i i +1
TOT
(B3.66)
T
C3 = [∆U ]+ δU ˆ ⋅ i [∆U ]+ δU ˆ − ∆l 2
i i +1 i +1
i +1
2
[δ λ ] = − C 2 ± C 2 − 4C1C 2 ( )1/ 2
⇒
i +1 [δλ ]1
i +1 (B3.67)
2C1 2C1 [δλ ]2
{ [λ]+
i i +1
[δ λ]1}f
∆l
{ [λ]+
i i +1
[δ λ]2 }f
i +1 i +1
i
[∆U]+ ˆ i +1
δU +
[U ]TOT ⋅i +1[δ λ]2 i
[∆U]+ ˆ i +1
δU +
[U]TOT ⋅i +1[δ λ]1
Una vez obtenido los dos factores de carga posible i +1 [δ λ ]1, 2 , es necesario determinar cual de
ellos es el correcto. Para esta finalidad se explora a cerca de la dirección de avance correcta,
[δU ]=
i +1
i +1
[U]TOT ⋅i+1 [δ λ ]1,2 , que se supone que es aquella cuyo ángulo con el incremento
ˆ i +1
δU +
θ max ∝ γ max
{ [λ]+
i i +1
[δ λ ]1 } f
∆l
{ [λ]+
i i +1
[δ λ]2 } f
θ min ∝ γ min
[f ] t
i +1
i
[∆U ]
i
[∆U]+ ˆ
δU
i
[∆U]+i +1 δUˆ +i +1[U]TOT ⋅i +1[δ λ] 2 i
[∆U]+ i +1 δUˆ + i +1[U ]TOT ⋅i +1[δ λ ]1
Figura B3.14 – “Arc-Length” camino esférico – Detalle de avance.
(B3.68)
γ 2 = [∆U ]2 ⋅ [∆U ]
i +1 T i
B4.1 Introducción.
En la primera parte de este capítulo se particulariza la ecuación del movimiento (B3.1) a
problemas elásticos lineales con el objetivo de estudiar la convergencia de la solución para
distintos métodos numéricos de resolución en el tiempo. El concepto de convergencia no
puede garantizarse en el sentido estricto en las ecuaciones diferenciales de segundo orden
no lineales, tal como la ecuación estudiada en el capítulo B3, debido a que la convergencia
implica estabilidad en la solución y es precisamente ésta la que no puede garantizarse. No
obstante, se estudia el concepto de “estabilidad linealizada”, que es el más utilizado y que
garantiza sólo condiciones necesarias de estabilidad, pero no suficientes.
El análisis de estabilidad es importante en el estudio de la solución dinámica en el tiempo
de un problema, porque junto a la condición de consistencia, determinan la convergencia del
algoritmo de solución. Dicho de otra manera, cuando en un sistema dinámico lineal se
cumplen las condiciones de estabilidad y consistencia, se tiene garantizada la convergencia
de la solución.
Existen varios caminos para estudiar la estabilidad de un algoritmo de solución de la
ecuación diferencial del equilibrio dinámico, pero el denominado estudio espectral de
Fourier es el mas utilizado y es también el que se presentará en este capítulo.
&& (t ) + D U
MU & (t ) + K U (t ) − f ext (t ) = 0 ∈ Ω (B4.1)
además se exige en esta ecuación que la matriz de masa M sea simétrica y definida positiva, en
tanto las matrices de amortiguamiento D y rigidez K deberán ser simétricas y semi-definidas
positiva. Al igual que en el capítulo anterior, este sistema de ecuaciones diferenciales ordinarias a
coeficientes constantes ya contiene la discretización espacial, tal que el campo de desplazamientos
& (t ) y aceleraciones U
U(t ) , velocidades U && (t ) , están definidos sólo en algunos puntos del medio
continuo (nodos), más precisamente en aquellos puntos donde se soporta el polinomio de
aproximación local o también denominado función de forma de un elemento finito (ver
apartado B2.5). Es necesario puntualizar que en la ecuación (B4.1) el campo temporal es
continuo. Dentro del espacio discreto elástico lineal, aproximado mediante funciones
polinómicas, se representan las siguientes magnitudes (ver ecuación B2.77),
[ ]
Las fuerzas externas : f ext (t + ∆t ) = A ∫ e N ⋅ t dS + ∫ e ρ N ⋅ b dV
Ωe S V
las magnitudes que intervienen en esta ecuación han sido definidas en los capítulos previos.
ϕ1 (t )
M
U (t ) = U ⋅ Φ(t ) = [U 1 L U h L U n ]⋅ ϕ h (t ) (B4.5)
M
ϕ n (t )
♣
Nota: Las propiedades de ortogonalidad de los autovectores establece,
0 ∀ i ≠ j 0 ∀ i≠ j
U Ti M U j = ; U Ti K U j = 2
1 ∀ i = j ω i = λ i ∀ i = j
siempre que estos autovectores estén normalizados respecto a la masa.
4-4 ANÁLISIS DE ESTABILIDAD EN LA SOLUCIÓN
&& (t ) + D U
MU & (t ) + K U (t ) − f ext (t ) = 0 && t + ∆t + D U
MU & t + ∆t + K U t + ∆t − [f ext ]t + ∆t = 0
U (t = 0) = U 0 & (t = 0) = V
; U & t + ∆t = U
U & t + (1 − γ ) U && t + ∆t ∆t
&& t ∆t + γ U
0
& t ∆t + 1 − β U
U t + ∆t = U t + U && t + ∆t ∆t 2
&& t ∆t 2 + β U
Descomposición MODAL
2
Discretización temporal. U(t = 0) = U 0 ; &
U(t = 0) = V0
(Newmark)
Descomposición MODAL
donde τ t representa la matriz de error local y el asterísco (*) implica que se trata de una
magnitud aproximada. Se dice que el algoritmo de solución cumple con la condición de consistencia, si
la matriz A es espectralmente estable y para un tiempo t + θ∆t , con θ = 1 , se cumple la siguiente
condición,
τ t = c ⋅ ∆t K +1 ∀ t ∈[0, T ] (B4.11)
donde c es una constante independiente del incremento de tiempo ∆t y K > 0 es la
relación de convergencia.
La ecuación que expresa el error, en la solución obtenida por diferencias finitas, surge de
establecer la diferencia entre la solución correcta y la aproximada, esto es, restando miembro a
miembro ambas soluciones se tiene,
6 T. Belytschkp and T. Hughes (1983). Computational methods for transient analysis. North-Holand.
7 M. Greenberg (1978). Foundation of applied mathematics. Prentice-Hall.
4-6 ANÁLISIS DE ESTABILIDAD EN LA SOLUCIÓN
X t + θ∆t − A θ ⋅ X t − Lt +θ∆t = 0
∗ t +θ∆t (B4.12)
X − A θ ⋅∗ X t − Lt + θ∆t − τ t = 0
de donde surge la siguiente expresión del error en el instante actual, es decir en el tiempo
( t + θ∆t ),
(X t + θ∆t
) (
− ∗ X t +θ∆t − A θ ⋅ X t − ∗ X t − τ t = 0 ) (B4.13)
e t +θ∆t − A θ ⋅ e t − τ t = 0
Siguiendo la misma idea conceptual, se obtiene el error en el paso de tiempo previo, es decir
en ( t ),
e t − A θ ⋅ e t −θ∆t − τ t −θ∆t = 0 (B4.14)
Sustituyendo esta última ecuación, que define el error en el paso previo (B4.14), en la
ecuación del error en el instante actual (B4.13), con el objetivo de eliminar e t de esta última, se
tiene,
[
e t +θ∆t − A θ ⋅ A θ ⋅ e t −θ∆t + τ t −θ∆t − τ t = 0 ]
(B4.15)
e t +θ∆t − [A ] ⋅ e
θ 2 t −θ∆t
[ ]
− A θ ⋅ τ t −θ∆t − τ t = 0
Repitiendo ahora el proceso deductivo para eliminar el error e t −θ∆t de la ecuación anterior,
resulta,
[ ] 3
[ ] 2
[ ]
e t +θ∆t − A θ ⋅ e t − 2θ∆t − A θ ⋅ τ t − 2θ∆t − A θ ⋅ τ t −θ∆t − τ t = 0 (B4.16)
Y siguiendo así sucesivamente se puede llegar a escribir el error acumulado hasta el tiempo
actual, ( t n = n ⋅ ∆t ). Es decir durante todo el proceso de solución,
[ ] [ ]
n
n +1
⋅ e t =0 − ∑ A θ
i
e t +θ∆t − A θ ⋅ τ t −i θ∆t = 0 (B4.17)
i =1
Puesto que el segundo término es mucho menor que la acumulación del error en los restantes,
[A ] θ n +1
[ ]
⋅ e t =0 << ∑i =1 A θ ⋅ τ t −i θ∆t , puede despreciarse y directamente obtener el error
n i
∑ [A θ ]
n −1
i
et = ⋅ τ t −(i −1) θ∆t ≤ t ⋅ c ⋅ ∆t K (B4.18)
i =1
Esta conclusión, conocida también como “teorema de equivalencia de Lax” para problemas
de valores iniciales7, conduce a un hecho obvio, pues cuando el incremento de tiempo tiende a
cero desaparece el error en la solución aproximada por diferencias finitas ( e t → 0, si ∆t → 0 ).
DINÁMICA NO LINEAL 4-7
Para que esta solución sea estable debe ocurrir que la matriz A θ este acotado para θ → ∞ , y
para ello el máximo autovalor max µ k de la matriz A (o radio espectral ρ(A ) ), debe ser
k
menor o igual que la norma de dicha matriz ρ(A ) ≤ A . Así, la estabilidad está relacionada con la
relación de crecimiento, o disminución de A θ en función de la potencia θ y recibe también el
nombre de condición de estabilidad espectral. Para probar la convergencia es entonces necesario
que se cumpla la siguiente condición,
A θ ≤ const. , ∀θ (B4.20)
Este requerimiento se satisface si se cumple:
1. que el radio espectral sea menor que la unidad, ρ(A ) ≤ 1 ; siendo el radio
espectral de la matriz A , la magnitud de su máximo autovalor. Esto es,
ρ(A ) = max µ k con k = 1, 2,3 (B4.21)
k
Una matriz A que cumple con estas dos condiciones se dice que es estable o espectral-
mente estable.
A= A1−1 ⋅ A2 , L = t
A1−1
∆t 2
2
[(1 − 2β) [a
+ 2β[a hext ]t + ∆t
ext t
h ]
] (B4.24)
∆t [(1 − γ ) [a ext t
h ] + γ[ a ext t + ∆t
h ] ]
Donde las matrices que componen la matriz de integración A , tienen las siguientes
expresiones,
1 + ∆t 2 βω h2 2∆t 2 βξ h ω h
A1 = 2
∆t γ ω h 1 + 2 ∆t γ ξ h ω h
∆t 2 (B4.25)
1− (1 − 2 β) ω h2 ∆t [1 − ∆t (1 − 2β )ξ h ω h ]
A2 = 2
− ∆t ( 1 − γ ) ω h 1 − 2∆t (1 − γ ) ξ h ω h
2
♣
NOTA: Un camino alternativo es no eliminar la aceleración y sustituir las ecuaciones (B4.23) en la (B4.8),
&& th+ ∆t y luego sustituir nevamente ésta en la (B4.23), de donde se obtiene
de donde se obtiene la aceleración ϕ
ϕ th+ ∆t y ϕ& th+ ∆t . Reordenando estas ecuaciones así obtenidas y sustituyendola en (B4.9), para θ = 1 , se tiene,
β
ϕ && th + ∆t
[( ) 1
− − β x − 2(1 − γ ) y
2
] [− x − 2 y ] 1
∆t
−
x
2
∆t ϕ && th
ω 2 ∆t 2
h
t + ∆t
ϕ& h
[ ( ) 1
= ∆t 1 − γ − 2 − β γ x − 2(1 − γ ) γ y ] [1 − x γ − 2 γ y ] −
x γ t
⋅ ϕ& h
∆t t
β γ
+ 2
ext t + ∆t
⋅ [a h ]
t + ∆t
ϕ h
[ ( )
∆t 2 − β − − β β x − 2(1 − γ ) β y
1
2
1
2
] ∆t [1 − x β − 2 βy ] (1 − β x) ϕ h
ω h ∆t
β x
2
ω h
−1
1 2 ξh γ ξh x
Siendo x = 2 2 + + β e y = . El inconveniente de este procedimiento es que la
ωh ∆t ωh ∆t ω h ∆t
obtención del radio espectral resulta de una ecuación cúbica, siendo más simple utilizar un procedimiento que
condense la aceleración y permita obtener una ecuación cuadrática para el radio espectral.
DINÁMICA NO LINEAL 4-9
2
A + A22 A + A22
µ1, 2 = 11 ± 11 − ( A11 A22 − A12 A21 ) = a1 ± a12 − a 2 (B4.27)
2 2
con a1 = ( A11 + A22 ) / 2 y a 2 = ( A11 A22 − A12 A21 ) . De esta última ecuación y de la condición
de estabilidad ρ(A ) ≤ 1 , puede establecerse gráficamente el dominio donde hay estabilidad de la
solución y cuya frontera queda delimitada a partir de hacer µ = 1 y adoptar el signo negativo en la
ecuación (B4.27) (ver Figura 4.2). De esto resulta la siguiente relación,
0 = 1 − 2a1 + a 2 (B4.28)
Graficando las posibles soluciones en el espacio a1 - a 2 , se observa que se tiene estabilidad en
la solución de este algoritmo cuando la ecuación (B4.27) tiene raíces dentro de la zona gris de la
Figura 4.2. Es decir, para que sea incondicionalmente estable debe cumplirse la siguiente relación,
que establece los vértices de la Figura 4.2,
a2 + 1 a +1
− ≤ a1 ≤ 2 (B4.29)
2 2
Sustituyendo en la ecuación (B4.28) 2 a1 = ( A11 + A22 ) y a 2 = ( A11 A22 − A12 A21 ) por los
coeficientes de la matriz A (ecuación (B4.24)), resultan las magnitudes γ y β que garantizan una
solución estable con el algoritmo de Newmark. Los extremos de variación de estos coeficientes
cumplen con la siguiente desigualdad,
2β ≥ γ ≥ 1
2 (B4.30)
Un ejemplo de la resolución numérica de este problema puede verse en el “Apéndice 1”,
de este capítulo, donde se describe un pequeño programa en FORTRAN sobre la
obtención del radio espectral y con el que puede verificarse la condición expresada por la
ecuación (B4.30).
ϕ th+ ∆t = ϕ th +
∆t
2
[ϕ& t
h + ϕ& th+ ∆t ] (B4.32)
a2
a2 = 1
(ρ = 1)
a1=-1 a1=1
RAÍCES COMPLEJAS
a12 = a2 ρ<1
Raíces REGIÓN
dobles REGIÓN
INESTABLE
reales ESTABLE
a1
RAÍCES REALES
ρ<1
0 = 1 − 2a1 + a 2 0 = 1 − 2a1 + a 2
(ρ = 1) (ρ = 1)
a 2 = −1
(ρ = 1)
[
aˆ ht + ∆t = 2ξ h ω h ∆t + 2ω 2h ∆t 2 ϕ ]
&& th − + 2ξ h ω h ∆t 2 ϕ& th +
4 (B4.34)
∆t
t t 2 t ext t + ∆t
+ϕ&& h + 2ξ h ω h ϕ& h + ω h ϕ h − [a h ]
1444424444 3
aht
0 = fˆ t + ∆t + J t + ∆t ⋅ ∆U t + ∆t
t + ∆t
4 2
con : J t + ∆t = M 2 + D T + KT
∆ t ∆ t
(B4.35)
4
(
y : fˆ t + ∆t = M 2 ∆U t + ∆t − U && t +
& t ∆t − U
)
∆t
[ ]
+ f int
t + ∆t
[ ]
− f ext
t + ∆t
2 − ω h2 ∆t 2 − 1 1 0
det − µk = 0
1 0 0 1
(B4.40)
( )
µ k2 − µ k 2 − ω h2 ∆t 2 + 1 = 0 ⇒ µ1, 2 =
(2 − ω h2 ∆t 2
)± (
2 − ω h2 ∆t )
2 2
−1
2 4
Imponiendo a los autovalores previamente obtenidos la condición de radio espectral menor o
igual que la unidad (apartado B4.5.1), se obtiene la medida del incremento de tiempo ∆t que
garantiza la estabilidad de la solución para el método de diferencias finitas centradas con
amortiguamiento nulo ξ h = 0 . Esto es,
2 ∆t 1
ρ(A ) = max µ1 , µ 2 ≤ 1 ⇒ ∆t ≤ o también ≤ (B4.41)
ωh T π
∆t x x
2 1 − 2 − y 1 − 2 − y
A= 2
(B4.43)
∆t
∆t
2
−1
ξ ξh x
Siendo en esta última x = 21 2 + h e y= .
ω ∆t ω h ∆t ω h ∆t
h
1 0
A 1 = ∆t 2
2 ω h 1 + ∆t ξ h ω h
∆t 2 2 ⇒
1 − ωh ∆t [1 − ∆tξ h ω h ]
A2 = 2
− ∆t 2
ωh 1 − ∆t ξ h ω h (B4.44)
2
∆t 2 2
1− ωh ∆t [1 − ∆tξ h ω h ]
2
⇒ A = A 1−1 A 2 = 2
(2 2
ω h ∆t − 2 + ω h ∆t − 2∆t ) ( )
− ω 2h ∆t 2 + ω3h ∆t 3 ξ h + 2 − 2ω h ∆t ξ h
4[1 + ∆tξ h ω h ] 2[1 + ∆tξ h ω h ]
Tal que luego de obtener el radio espectral, se concluye que sólo se consigue estabilidad en la
solución siempre que se cumpla un incremento de tiempo en función de la máxima frecuencia
angular y del término de amortiguación,
(ω h )max ∆t ≤ 2( 1 + ξ 2 − ξ) (B4.45)
Esta condición comporta un coste elevado de cálculo cuando se resuelve un problema de
elementos finitos, por ejemplo, en un problema con viscosidad nula ξ = 0 , primero habrá que
obtener la máxima frecuencia angular de la estructura para cada modo “h” y luego exigir que el
paso de tiempo cumpla con (ω h )max ∆t ≤ 2 (ver este mismo resultado en la ecuación (B4.41)).
Esto podría solucionarse a través de una forma simple, obteniendo las frecuencias para cada
elemento finito y buscando la máxima de ellas, pues esta siempre cumplirá con la condición9,
(ω h ) max
≤ max ω
∀e
( )e max
≤
2( 1 + ξ 2 − ξ )
∆t
(B4.46)
9 F. Cesari (1982). Metodi di calcolo nella dinamica delle strutture. Ed. Pitagora, Bologna.
4-14 ANÁLISIS DE ESTABILIDAD EN LA SOLUCIÓN
∂U Ω
i t + ∆t
∂U&& & ∂f ext
T ∂U
(B4.51)
=
i
[∆f ] t + ∆t
Ω + M T
+ K +D − ⋅
i +1
[∆U ]tΩ+ ∆t
∂U ∂U ∂U
Ω
1444444 424444444 3
i J t + ∆t
Ω
Recordar que en esta ecuación (i) representa el contador de iteración y (t) el tiempo. La
solución de la misma resulta por inversión del operador jacobiano i J tΩ+ ∆t . Esto es,
i +1
[∆U]tΩ+ ∆t = − [ i J tΩ+ ∆t ]−1 ⋅i [∆f ]Ωt + ∆t (B4.52)
Diversos algoritmos en el tiempo pueden utilizarse para la solución de esta última ecuación,
pero con ninguno de ellos puede garantizarse la estabilidad en la solución. Se supone en esta
formulación que M es constante; M , K y D son simétricas; M y K son definidas positivas; y
D es semi definida positiva.
La utilización del concepto de “estabilidad de la ecuación linealizada” es sin duda el método
mas difundido para probar estabilidad en la solución en problemas no lineales. En este caso, el
concepto de estabilidad para un algoritmo lineal, estudiado en apartados previos, sólo puede
garantizarse la condición necesaria, pero no suficiente, para que exista estabilidad.
∂ ∂ ∂ int ∂
f int ( U) =
σ ∫
∂U t
Pint dt =
∂U V
∫ ρ ω dV + ∫ K& dt ⇒ K T =
f (U) =
∂U σ
∂U V∫ ρω dV
(B4.53)
144 424t 44 3
&)
E (U,U
En esta última puede verse que para una fuerza exterior nula la energía se conserva y es
constante,
Sí f ext (τ) = 0 ⇒ & ) = E (U , U
E ( U, U & )
0 0 (B4.55)
En esta situación, la regla del trapecio puede escribirse como,
[
&& t + ∆t + f int ( U)
∆f = M U
σ
] t + ∆t
[ ]
− f ext
t + ∆t
=0 (B4.56)
que linealizada conduce a la forma (B4.35). En la ecuación (B4.56), según la regla del trapecio, la
velocidad y el desplazamiento se expresan como,
4-16 ANÁLISIS DE ESTABILIDAD EN LA SOLUCIÓN
t + ∆t
U = Ut +
∆t & t & t + ∆t
2
U +U [ ]
(B4.57)
U
& t + ∆t U
& t + ∆t = U
2
&& t + ∆t
&& t + U [ ]
Cuyos valores iniciales son,
U (t = 0 ) = U
0
U& (t = 0 ) = U& (B4.58)
[ ]
0
&& && −1 ext int
U(t = 0) = U 0 = M f − f σ (U 0 )
El problema algebraico no lineal en U t + ∆t se obtiene por sustitución de las ecuaciones
(B4.57) en (B4.56).
En lo que sigue conviene observar el problema algebraico no-lineal, linealizado en el tiempo
mediante el procedimiento de Newmark como un funcional de Euler-Lagrange. Esta forma
coincide también con el denominado funcional de Hamilton7 expresado en ( t + ∆t ), cuya
expresión es,
π(U t +∆t ) = E (U t +∆t , U& t + ∆t ) − [U t +∆t ]T ⋅ { fˆ t + ∆t
[ ]
− f int
t + ∆t
} (B4.59)
Donde f̂ t + ∆t es la expresión de la fuerza externa linealizada para el algoritmo de Newmark,
cuya expresión se muestra en la ecuación (B4.35). Sustituyendo en la ecuación (B4.59) la
expresión de E ( U t + ∆t , U& t + ∆t ) y f̂ t + ∆t , resulta el funcional buscado,
t + ∆t
π(U t + ∆t 2
[
) = 2 U t + ∆t
∆t
]
T
⋅M⋅U t + ∆t
+ ∫ ρω dV −
V
(B4.60)
− U[ ] [ ]
t + ∆t T
⋅ f ext
t + ∆t && t
+ M ⋅ U
4 &t
+ 2 ∆U t + ∆t U
∆
( )
t
Cuya primera variación, que hace estacionario el funcional, resulta
δπ(U t + ∆t ) = δ U[ t + ∆t T
] ⋅ ∂∂Uπ =δU [ t + ∆t T
] ∆4t 2
M ⋅ U t + ∆t + f int (U t + ∆t ) −
σ
t + ∆t
(B4.61)
[ ]
− f ext
t + ∆t && t
− M ⋅ U
4 &t
+ 2 ∆U t + ∆t U
∆t
( ) = 0
Tal que esta última es nula para cualquier variación de δ U t + ∆t ≠ 0 . [ ]
Regla del trapecio modificada.
Esta condición L (U t +∆t ) = 0 , combinada con las ecuaciones (B4.57), implica que la
1
condición de conservación de la energía, E t +∆t − E t =
2
(∆U t + ∆t )T ⋅ ∆f ext
t + ∆t
[
, queda satisfecha. ]
Considerando ahora el funcional de Euler-Lagrange π(U t + ∆t ) , que define el problema
algebraico no lineal mediante la regla del trapecio, y la ecuación de restricción L (U t +∆t ) = 0 ; se
formula a continuación un funcional de Euler-lagrange con restricción de conservación de la energía a través
del uso del multiplicador de Lagrange mλ . Esto es,
π(U t +∆t ) = π(U t +∆t ) + mλ L (U t +∆t )
1 424 3 1424 3 (B4.64)
funcional restricción
Haciendo nula la primera variación, se encuentra el campo de desplazamiento que hace
estacionario este funcional.
∂π t + ∆t t + ∆t
δπ(U t + ∆t
[
)=δU ]
t + ∆t T
⋅ + mλ
∂L t + ∆t
+ δmλ L (U ) = 0 (B4.65)
δU δU
Sustituyendo en esta los respectivos funcionales, y haciendo nula la expresión para cualquier
[ ] T
variación no nula de δ U t + ∆t y δmλ , resulta el siguiente sistema algebraico que representa la
ecuación del movimiento con la condición de conservación de la energía,
4 m
0 = (1 + mλ ) 2 M ⋅ U t + ∆t + f int (U t + ∆t ) − 1 − λ f ext
σ
[ ] t + ∆t
− f [ ]
mλ ext t
−
∆t 2 2
(B4.66)
t m & t
− M ⋅ U&& + 4 2 (1 + mλ )U t + ∆t 1 + λ U = 0
∆t 2
La solución simultánea de las ecuaciones (B4.63) y (B4.66), mediante un algoritmo iterativo,
permiten obtener el multiplicador de Lagrange mλ que controla la fuerza introducida al sistema, y
los campos de desplazamiento, velocidad y aceleración buscados.
4-18 ANÁLISIS DE ESTABILIDAD EN LA SOLUCIÓN
APÉNDICE - 1
Programa para la obtención del radio espectral de una matriz.
∆t / T
PROGRAM estabilidad
C-----------------------------------------------------------------------
C
C ESTE PROGRAMA RESUELVE EL PROBLEMA DE ESTABILIDAD
C POR EL METODO DEL RADIO ESPECTRAL
C
C ESTA APLICADO SOBRE EL METODO DE NEWMARK
C
C Entrada de datos: archivo.dts
C Salida de resultados RADIO ESPECTRAL vs DT/T: archivo.res
C Salida de resultados de Vect. y Val. propios: archivo.sal
C
C MANUAL DE ENTRADA DE DATOS
C
C Ng : Numero de grados de libertad del sistema,
C Dtini : Incremento de tiempo inicial
C Dtfin : Incremento de tiempo final
C T : Periodo
C B : Beta Newmark
C G : Gama Newmark
C Psi : Amortiguamiento
C Iread : 1 (calcula para la matriz de Newmark)
C 2 (calcula para la matriz de Dif. Centrada)
c------------------------
C si Iread = 2, entonces:
C
C COEFICIENTES PARA MONTAR LA matriz
C
C COEFmatrix(1 : 1)
C COEFmatrix(1 : 2)
C COEFmatrix(1 : 3)
C .
C .
C COEFmatrix(Ng:Ng)
C
C-----------------------------------------------------------------------
IMPLICIT INTEGER*4 (I - N)
IMPLICIT REAL*8(A-H,O-Z)
complex (8) RAIZ,vl1,vl2
PARAMETER(NP=2)
DIMENSION AUX1(NP,NP),AUX2(NP,NP),AUX3(NP,NP),
. COEFmatrix(NP,NP),A1(NP,NP),A2(NP,NP),
. A1inv(NP,NP)
DIMENSION DVAL(NP),DVEC(NP,NP),RES(NP)
DIMENSION IEIG(NP)
CHARACTER ARCH1*20,ARCH2*20,ARCH3*20,CN*20
C
PI=DACOS(-1.0D00)
C
C*** PUESTA A CERO DE ALGUNAS MATRICES
C
DO Ig=1,NP
DO Jg=1,NP
AUX1(Ig,Jg)=0.0
AUX2(Ig,Jg)=0.0
AUX3(Ig,Jg)=0.0
IF(Jg.EQ.Ig)AUX2(Ig,Jg)=1.0
ENDDO
ENDDO
4-20 ANÁLISIS DE ESTABILIDAD EN LA SOLUCIÓN
1111 CONTINUE
IF(Iread.EQ.1)THEN !NEWMARK
Dt2=Dt*Dt
W2=W*W
C matriz A1
A1(1,1)=1.0+(Dt2*B*W2)
A1(1,2)=2.0*Dt2*B*Psi*W
A1(2,1)=Dt*G*W2
A1(2,2)=1.0+(2.0*Dt*G*Psi*W)
C matriz A2
A2(1,1)=1.0-((0.5*Dt2)*(1.0-2*B)*W2)
A2(1,2)=Dt*(1.0-(Dt*(1.0-2.0*B)*Psi*W))
A2(2,1)=-Dt*(1.0-G)*W2
DINÁMICA NO LINEAL 4-21
A2(2,2)=1.0-(2.0*Dt*(1.0-G)*Psi*W)
C Inversión de la matriz A1
CALL INGAUS(A1,A1inv,NP,NG)
DO Ig=1,Ng
DO Jg=1,Ng
COEFmatrix(Ig,Jg)=0.0
DO Kg=1,Ng
COEFmatrix(Ig,Jg)=COEFmatrix(Ig,Jg)+
. A1inv(Ig,Kg)*A2(Kg,Jg)
ENDDO
ENDDO
ENDDO
ELSE ! Diferencias Centradas
Dt2=Dt*Dt
W2=W*W
C matriz A1
x=1.0/((1.0/(W2*Dt2))+(Psi/(W*Dt)))
y=Psi*x/(W*Dt)
A1(1,1)=0.5*Dt*(1.0-x*0.5-y)
A1(1,2)=(1.0-x*0.5-y)
A1(2,1)=0.5*Dt2
A1(2,2)=Dt
DO Ig=1,Ng
DO Jg=1,Ng
COEFmatrix(Ig,Jg)=A1(Ig,Jg)
ENDDO
ENDDO
ENDIF
WRITE(3,903)Dt
903 FORMAT(' Matriz Espectral a Analizar, para',1F12.5)
DO 3 Ig=1,Ng
WRITE(3,999)(COEFmatrix(Ig,Jg),Jg=1,Ng)
DO 3 Jg=1,Ng
AUX1(Ig,Jg)=COEFmatrix(Ig,Jg)
3 CONTINUE
A1coef=0.5*(AUX1(1,1)+AUX1(2,2))
A2coef=AUX1(1,1)*AUX1(2,2)-AUX1(1,2)*AUX1(2,1)
DIS=(A1coef*A1coef)-A2coef
cc
RAIZ=CDSQRT(DCMPLX(DIS))
vl1=(A1coef+RAIZ)
vl2=(A1coef-RAIZ)
c
4-22 ANÁLISIS DE ESTABILIDAD EN LA SOLUCIÓN
Pent1=DREAL(vl1)
Pima1=DIMAG(vl1)
vmod1=dsqrt(Pent1*Pent1+Pima1*Pima1)
c
Pent2=DREAL(vl2)
Pima2=DIMAG(vl2)
vmod2=dsqrt(Pent2*Pent2+Pima2*Pima2)
c
RADIO=vmod1
IF(vmod2.ge.radio)RADIO=vmod2
C
C*** IMPRESION DE LA RESPUESTA
C
c WRITE(3,850)1,Pent1,Pima1
c WRITE(3,850)2,Pent2,Pima2
WRITE(2,100)(Dt/T),radio
C
Dt=Dt+Dinc
IF(Dt.LT.Dtfin)GOTO 1111
C
STOP
100 FORMAT(3(1X,F12.5))
850 FORMAT(/,5x,I5,2x,' Auto VALOR -- w=',1x,1E12.5,'i',1E12.5)
999 FORMAT( 30(1X,1E12.5))
END
C-----------------------------------------------------------------------
SUBROUTINE INGAUS(A,B,NP,N)
C***********************************************************************
C SUBRUTINA PARA INVERSION DE MATRICES POR ELIMINACION GAUSIANA
C PARAMETROS:
C
C A ====> MATRIZ A INVERTIR
C B ====> MATRIZ INVERTIDA
C N ====> ORDEN DE LA MATRIZ
C***********************************************************************
IMPLICIT INTEGER*4 (I - N)
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION A(NP,NP),B(NP,NP)
C
DATA TOLSL /1.0E-30/
DO I=1,N
IF(DABS(A(I,I)).LT.TOLSL)A(I,I)=1.0E-15
ENDDO
DO I=1,N
DO J=1,N
B(I,J)=0.0
IF(I.EQ.J)THEN
B(I,J)=1.0D0
ENDIF
ENDDO
ENDDO
DO I=1,N
P=A(I,I)
IF(ABS(P).LT.TOLSL) THEN
PRINT '(A72)',
. ' ERROR EN LA RUTINA DE INVERSION DE MATRICES'
DINÁMICA NO LINEAL 4-23
STOP
ENDIF
TEMP=1.0D0/P
DO J=1,N
B(I,J)=B(I,J)*TEMP
A(I,J)=A(I,J)*TEMP
ENDDO
DO K=1,N
IF(K.NE.I) THEN
P=A(K,I)
DO J=1,N
A(K,J)=A(K,J)-P*A(I,J)
B(K,J)=B(K,J)-P*B(I,J)
ENDDO
ENDIF
ENDDO
ENDDO
RETURN
C
END
C-----------------------------------------------------------------------
4-24 ANÁLISIS DE ESTABILIDAD EN LA SOLUCIÓN
APÉNDICE - 2
Ecuación diferencial vinculada a un funcional – Funcional natural de Euler-Lagrange.
x =b
∂I
x =b
d ∂I ∂I
δπ = 0 = ∫ − δf dx + δf x (B4.69)
∂f
x=a
dx ∂f x ∂f x x=a
∂I d ∂I
a) La ecuación de Euler - Lagrange − = 0 , ∀ a ≤ x ≤ b,
∂f dx ∂f x
x =b (B4.70)
∂I
b) La condición de contorno δf x =0 en : x = a y x=b
∂f x x=a
B 5 Modelos
independientes del
tiempo((**)).
B 5.1 Introducción.
En este capítulo se hace un recordatorio de algunos conceptos básicos de la teoría de la
elasticidad y sus variables mecánicas, y principalmente se presenta en forma breve la teoría de
la plasticidad clásica y una modificación de la misma que la hace más general y también se
hace una breve presentación de la teoría de daño continuo. Todos estos desarrollos se harán
dentro del marco de una cinemática con pequeños movimientos, que supondremos
también que introduce pequeñas deformaciones. Es aconsejable apoyar estas lecturas con
consultas a las referencias básicas de la mecánica de medios continuos1,2,3, donde se
encuentras las respuestas a cuestiones más profundas sobre este tema. No obstante es
obligado establecer los criterios, hipótesis y notaciones, y recordar aquellos conceptos que
se consideran de mayor importancia para el tema que se desarrolla en este trabajo.
1 - Condición de equilibrio
∂σ ij
= −ρbi (B5.1)
∂x j
1 ∂u i ∂u j
ε = ∇Su =
1
2
(
∇ 0 u + ∇ T0 u ) ; ε ij = ∇ Sj u i = +
2 ∂x j ∂x i
(B5.3)
tal que en estas tres últimas ecuaciones t i es la fuerza de superficie aplicada sobre el
contorno S , σ ij es el tensor de tensiones de Cauchy, n j el vector normal a la superficie S
que envuelve el sólido1,2,3, bi fuerzas de volumen por unidad de masa; ρ = ∂ M ∂V la
densidad de masa, M la masa y V el volumen; ui es el campo de desplazamientos y ε ij el
tensor de deformación infinitesimal.
Las ecuaciones que definen los modelos constitutivos elásticos clásicos están bien
establecidas y son las siguientes:
σ ij = f ij (ε kl ) ; ε ij = f ij−1 (σ kl ) (B5.4)
Así resulta la tensión a partir de un modelo cuya variable libre es la deformación,
o se obtiene la deformación en los modelos basados en variables libres de tensión.
Estas relaciones son invertibles y también reversibles, por lo que en elasticidad no
hay disipación de energía.
“Modelo elástico de Green”, o también llamado hiperelástico, en el cual la variable
del problema depende de una densidad de potencial Ψ = Ψ (ε ij ) , o de su complemento
( )
Ψ = Ψ σ ij , que debe ser preestablecida,
( )
∂Ψ ε ij ( )
∂Ψ σ ij
σ ij = ; ε ij = (B5.5)
∂ε ij ∂σ ij
(
σ& ij = f ij ε& ij ; σ mn ) ; (
ε& ij = g ij ε ij ; σ& mn ) (B5.6)
En todos los modelos antes citados siempre debe cumplirse el concepto de reversibilidad
termodinámica e independencia entre tensiones y trayectoria.
Una relación lineal muy establecida y que cumple con las tres definiciones anteriores es la
ley de Hooke generalizada,
DINÁMICA NO LINEAL 5-3
σ ij = C ijkl ε kl (B5.7)
donde aparece un operador mecánico, denominado tensor constitutivo C ijkl . Este tensor es
de cuarto orden y tiene 81 componentes en su expresión más general, pero para el caso de
elasticidad lineal e isótropa se reduce a dos componentes independientes. Estas
componentes son el módulo de Young E y el coeficiente de Poisson ν. En este caso de
isotropía y haciendo uso de las condiciones de simetría del tensor constitutivo de cuarto
orden2, puede escribirse como un tensor de segundo orden (matriz cuadrada) de 36
componentes.
Otras formas de escribir la relación elástica lineal (B5.7), es a partir de las constantes de
Lamé 1, en cuyo caso la tensión suele descomponerse en la contribución de dos sumandos,
uno que se refiere al efecto de los cambios volumétricos y otro al efecto de las desviaciones
o distorsiones. Esto es,
1 λ
σ ij = λε kk δ ij + 2µ(ε ij − 13 ε kk δ ij ) ⇒ ε ij = σ ij − σ kk δ ij
144244 3 2µ 2µ(3λ + 2µ ) (B5.8)
eij
donde λ = E 3(1 − 2ν) y µ = E 2(1 + ν) son las constantes de Lamé. También puede
definirse esta misma ley constitutiva elástica en la siguiente forma,
3νκ 1 ν
σ ij = 2Gε ij + ε kk δ ij ⇒ ε ij = σ ij − σ kk δ ij (B5.9)
1+ ν 2G 3κ(1 − 2µ )
3λ + 2µ p
Donde κ = = es el módulo volumétrico y en el se identifica la deformación
3 εv
volumétrica ε v = tr (ε ij ) = ε kk δ ij = 3 ε oct y la presión o tensión octaédrica p = 13 tr(σ ij ) =
τ ij
= 13 σ kk = σ oct . Además, en la ecuación (B5.9) se tiene G = µ = , en la cual se identifica
γ ij
γ ij = 2ε ij ; τ ij = σ ij .
σ 22 = 0, ε 22 ≠ 0
x2 σ11
σ 11 σ11 , ε11
a) E
1
x1 ε11
x2
σ12 σ12
b) µ
γ
1
x1 ε12
x3
p ε v = 3 ε oct - deformación volumétrica
p
p = σ oct - presión
c)
p ( )
ε v = tr ε ij = ε ii
κ
1
1 p = σ kk
x2 3
p
εv
x1
Como se sabe, esta forma tensorial es sólo una representación convencional del estado
tensional y deformacional de un punto de un sólido a través de sus componentes, según
tres planos ortogonales. Así pues, puede darse la paradoja de observar como el estado de
tensión o deformación en un punto varía según el plano que se considere, pero esto no es
real pues sólo varía la forma en que se representa dicho estado tensional. No obstante hay
un modo inequívoco de independizarse de esta forma inobjetiva de medición y esto es a
DINÁMICA NO LINEAL 5-5
J1 = 0 J 1' = ε ii = 0
1 1
J 2 = s ij s ij J 2' = e ij e ij (B5.13)
2 2
1 ' 1
J 3 = s ij s jk s ki J 3 = e ij e jk e ki
3 3
donde el tensor desviador de tensiones se representa mediante la siguiente expresión,
σ kk
s ij = σ ij − φδ ij = σ ij − δ ij , y el tensor desviador de deformación se expresa como,
3
ε ε
e ij = ε ij − v δ ij = ε ij − kk δ ij . Existen también otras formas matemáticas de expresar estos
3 3
invariantes y algunas de ellas también se verán en otras secciones de este libro.
B 5.3.1 Introducción.
A continuación se muestra una forma de representar el comportamiento constitutivo
elástico no-lineal ─reversible─ de un material ideal. Para ello es conveniente partir desde
una formulación hiperelástica que depende de la densidad de energía, cuya definición puede
realizarse en función del campo de deformación (variable libre del problema), como
t
Ψ= ∫ σ ij ε& ij dt (B5.14)
t =0
También podría definirse una ley constitutiva elástica no-lineal siguiendo los conceptos
de la hipoelasticidad, es decir, generalizando las ecuaciones lineales.
5-6 MODELOS INDEPENDIENTES DEL TIEMPO
ε
Figura B5.2 – Energía complementaria y de deformación.
(
Ψ = Ψ I 1' ; I 2' ; I 3' ) ⇒ σ ij =
∂Ψ
∂ε ij
Basado en deformación
(B5.16)
(
Ψ = Ψ I1 ; I 2 ; I 3 ) ⇒ ε ij =
∂Ψ
∂σ ij
Basado en tensión
∂Ψ ∂Ψ ∂I 1 ∂Ψ ∂J 2 ∂Ψ ∂Ψ
ε ij = = + = δ ij + δ ij
∂σ ij ∂I 1 ∂σ ij ∂J 2 ∂σ ij ∂I 1 ∂J 2 (B5.17)
ε ij = bJ 2 δ ij + (a + bI 1 )δ ij
1
I1 = σ ; J 2 = σ2
3
(B5.18)
2 σ
s 11 = σ ; s 22 = s 33 =− ; s 12 = s 13 = s 23 = 0
3 3
b 2 2 2a
ε11 = ε = σ + (a + bσ ) σ = σ 2 b + σ (B5.19)
3 3 3
donde es necesario parametrizar el modelo (obtener los parámetros a y b) y para ello
debe realizarse un ensayo de laboratorio, del tipo del que se muestra en la Figura B5.3.
DINÁMICA NO LINEAL 5-7
σ
σ
σB B
2a
ε A = σ2Ab + σA
3 ⇒ a , b
σA 2a
A ε B = σ2Bb + σB
3
σ
εA εB ε
Figura B5.3 – Modelo elástico no-lineal trabajando a tracción.
Requisitos de Estabilidad:
1. El trabajo realizado por el cambio de magnitud en los agentes externos, debe
ser siempre positivo.
∂Ψ (ε )
Ψ (ε ) = ∫ σ : dε ⇒ σ (ε ) =
ε
∂ε
(B5.22)
∂Ψ (σ )
Ψ ( σ ) = ∫ ε : dσ ⇒ ε (σ ) =
σ
∂σ
Toda relación constitutiva del tipo hiperelástica o de Green, cumple los criterios de
estabilidad antes citados, siempre que los potenciales de energía sean definidos positivos.
Para probar esto, considérese una relación constitutiva del tipo σ (ε ) = ∂Ψ (ε ) ∂ε , tal que
cualquier variación en los agentes externo provoca el siguiente cambio en la tensión,
∂σ (ε ) & ∂ 2 Ψ (ε ) &
σ& (ε ) = :ε = :ε (B5.23)
∂ε ∂ε ⊗ ∂ε
La condición necesaria y suficiente para que se cumpla el criterio de estabilidad de cargas
para todo el volumen, ecuación (B5.20) y también para cargas cíclicas, ecuación (B5.21), es
que todos y cada uno de los punto de este sólido realicen un trabajo específico de segundo
orden positivo,
∂ 2 Ψ (ε )
H ijkl = H ; det(H ) = det >0 (B5.26)
∂ε ⊗ ∂ε
Condición de convexidad:
A
A C
Ψ (σ ) = cte Ψ (σ ) = cte
∆σ : ε a > 0 ∆σ : ε a ≤ 0
σ& ESTABLE
ε& σ& : ε& > 0 σ&
ε&
ε ε
σ σ
σ&
ε&
INESTABLE
σ& σ& : ε& < 0
ε&
ε
ε
B 5.4.1 Introducción.
La teoría de la plasticidad representa el comportamiento mecánico de los sólidos
cargados dentro de un rango de aplicación que va más allá del que define la teoría de la
elasticidad. Existen varias teorías que, sin ser la plasticidad, cumplen con esta idea, por lo
tanto su definición necesita mayor precisión y para ello será presentada brevemente en este
apartado. Esta teoría matemática fue inicialmente formulada en 1872 para representar el
fenómeno de distorsión que sufre la red cristalina de los metales4 cuando son sometidos a
deformaciones. Actualmente se utiliza su estructura matemática para representar el
comportamiento macroscópico de un punto sólido sin que necesariamente, a escala
microscópica, se esté desarrollando un fenómeno de plasticidad de metales en el sentido
estricto1,2,3. La plasticidad en pequeñas deformaciones se caracteriza por suponer que las
deformaciones en un punto, ε = ε e + ε p , se descomponen en una parte elástica ε e y otra
plástica ε p irreversible. Esta última es la que induce a un comportamiento energético no
conservativo dependiente del camino recorrido. La formulación de la teoría de la
plasticidad se basa en la mecánica de los sólidos continuos y representa el comportamiento
físico macroscópico de un sólido ideal, con los siguientes rasgos característicos:
El punto A’ que marca la separación entre estos dos estados mecánicos se lo conoce
como limite de fluencia para los materiales metálicos, y como limite de discontinuidad para
los materiales friccionales, quedando definido a través de una función en el espacio de
tensiones que recibe el nombre de función de fluencia plástica o función de discontinuidad,
respectivamente.
4 Tresca, H. E. (1872). Mémorie sur l’coulement des corps solides. Acad. Sci. Paris, 20, 75-135.
DINÁMICA NO LINEAL 5-11
con ablandamiento
con endurecimiento
Zona plástica
Zona plástica
elástica
σ
Zona
B
C D
A'
E
O
ε
P e
ε ε
Figura B5.6 – Comportamiento uniaxial esquemático de un material elasto-plástico ideal.
La Figura B5.6, muestra en forma esquemática el comportamiento uniaxial de un punto
correspondiente a un material elasto-plástico ideal. Al comenzar el proceso de carga
presenta una zona elástica lineal que se mantiene hasta el punto A, llamado limite de
proporcionalidad. Seguidamente inicia un proceso elástico no-lineal hasta alcanzar el punto
A' llamado limite de elasticidad, a partir del cual comienza el comportamiento elasto-
plástico caracterizado por un decrecimiento sostenido del módulo de rigidez tangente
debido a la acción de los mecanismos inelásticos irreversibles. Si durante el
comportamiento elasto-plástico del punto se inicia un proceso de descarga, se observa que
sólo se recupera la parte elástica ε e del total de la deformación ε , quedando la
deformación plástica ε p como la parte no recuperable. Dentro del período elasto-plástico
se pueden distinguir tres regiones (ver Figura B5.6):
- Una donde hay crecimiento de la tensión, tramo (A'-C), que recibe el nombre de
zona elasto-plástica con endurecimiento.
- Otra donde el punto que se analiza no experimenta cambio de tensión, tramo (C-
D), y recibe el nombre de zona elasto-plástica perfecta o de endurecimiento nulo.
- Por ultimo, una zona donde la tensión decrece bajo crecimiento sostenido de las
deformación, tramo (D-E), período elasto-plástico con ablandamiento.
evolución de las fronteras del dominio elástico dentro del espacio de tensiones
(ver Figura B5.7).
F(σ; q ) = 0 (B5.29)
donde σ es tensor de tensiones de Cauchy, q el conjunto de variables internas
agrupadas en forma de matriz columna (“back stress”). Esta función establece el límite a
partir del cual se inicia el comportamiento no-lineal. Cualquier estado tensional fuera del
recinto encerrado por esta superficie es inadmisible, por lo tanto el proceso de resolución
elasto-plástico consiste en forzar a que el estado tensional se sitúe en la frontera o en el
interior de la función elasto-plástica (ver Figura B5.7). En un proceso uniaxial esta función
está perfectamente establecida por tratarse de un escalar que se compara con otro escalar
que representa el umbral entre un comportamiento elástico y otro plástico. Para
comportamientos multiaxiales, la función de fluencia se comporta como un traductor de
estados multiaxiales a uniaxiales. Un vez obtenida la tensión uniaxial equivalente al estado
multiaxial, se compara con un escalar que representa el umbral obtenido en laboratorio
para un problema uniaxial.
σ3
σ1
σ3
PLANO: σ1 − σ 2
σ1
σT
σc
− σI
σ1 = σ 2
σc σ2
σc
− σ II
plano octaédrico
σ1 + σ 2 + σ 3 = cte meridiano de tracción
ξ
ρ min
σc
σ1 = σ 2 = σ 3
ρ
− σ II
Planos meridianos de tracción máxima: Estos planos son ortogonales a los planos
octaédricos, y quedan inequívocamente definidos por la recta que describe el espacio
diagonal ( σ1 = σ 2 = σ 3 ), y por cada recta que describe el radio octaédrico ρ, cuando
θ = − 1 π 6 , − 5 π 6 , − 9 π 6 . Dichos planos cortan a los ejes de tensiones principales σ i , en
puntos de igual valor al de la tensión de tracción uniaxial σ T . La intersección de estos
planos con la superficie de fluencia, define curvas en el espacio de tensiones, que se
denominan funciones de fluencia según planos meridianos de tracción.
Planos principales: Son aquellos que quedan definidos por la intersección de dos, de las
tres direcciones de tensión principal. La intersección de estos planos con la superficie de
fluencia, define curvas en el espacio de tensiones que se denominan funciones de fluencia
según planos de tensión principal.
σ3
− σI ξ0
− σI σ1
meridiano de
tracción − σc
− σI
σ2
− σc
plano: Π
− σ II
plano
octaédrico
ξ
σc
σ1 = σ 2 = σ 3 − σ II
ρ
meridiano de compresión
− σ II
ε& iP = λ& s i ( ) (
→ ε& iP − ε& Pj = λ& s i − s j ) {i, j ∈1,2,3} (B5.40)
elevando al cuadrado ambos miembros y sumando componente a componente se obtiene,
DINÁMICA NO LINEAL 5-17
(ε& P
1 − ε& 2P ) + (ε&
2 P
2 − ε& 3P ) + (ε&
2 P
1 − ε& 3P )
2
[
= λ& 2 (s 2 − s 1 ) + (s 2 − s 3 ) + (s 1 − s 3 )
2 2 2
] (B5.41)
pero recordando que (s i − s j ) = (σ i − σ j ) , resulta de la ecuación anterior el factor de
consistencia plástico,
P
P J&2'
6 J& 2' = λ& 6 J 2 ⇒ λ& = (B5.42)
J2
donde J& 2' = (e& P : e& P ) es el segundo invariante del incremento temporal del tensor
P 1
2
1
desviador de deformaciones plásticas e& p y J 2 = (s : s ) el segundo invariante del tensor
2
desviador de tensiones s . Sustituyendo (B5.42) en (B5.39), se obtiene el tensor de
deformación plástica para la teoría de Prandtl-Reus,
P
J& 2' P s
ε& =P
s = J& 2' (B5.43)
J2 J2
ε = C −1 : σ + ε p = ε e + ε P (B5.44)
donde la deformación plástica ε P representa la variable interna fundamental del problema
elasto-plástico, cuya definición tiene la siguiente forma,
∂G(σ , q) &
ε& P = λ& H P = λ& =λ g (B5.45)
ε ∂σ
esta expresión, que también recibe el nombre de regla de normalidad normal a la
superficie de potencial plástico G(σ , q) = cte. , es una generalización de la ecuación
(B5.39), y λ es un escalar no negativo denominado parámetro de consistencia plástica que
se determina a partir de la propia condición de consistencia de Prager que se presentará
más adelante y que da la magnitud del incremento temporal de deformación plástica ε& P . La
función de potencial plástico se determina a partir de estudios experimentales y es la que
define la dirección del incremento temporal de deformación plástica.
Un caso particular de flujo plástico se produce cuando por hipótesis se adopta la
superficie de fluencia plástica como superficie de potencial plástico G(σ , q) ≡ F(σ , q) . En
este caso particular la ecuación (B5.45) se reduce a
5-18 MODELOS INDEPENDIENTES DEL TIEMPO
∂F(σ , q) &
ε& P = λ& =λ f (B5.46)
∂σ
y se dice que la regla de flujo es asociada a la superficie de fluencia plástica. En caso
contrario se dice que el flujo es no-asociado.
Si se adopta la función de von Mises como función de potencial, se tiene el caso
particular del flujo de Prandtl-Reus (ecuación (B5.39)),
ω ( )
& = σ : ε& = σ : ε& e + ε& P = σ : ε& e + σ : ε& P = ω
&e +ω
&p (B5.48)
Se conoce esta forma de escribir la variación temporal de la energía como elasticidad
desacoplada y sólo vale para problemas elasto-plásticos cuyas deformaciones elásticas son
infinitesimales y por lo tanto se acepta la hipótesis de aditividad de las deformaciones
(ecuación (B5.44) ).
Concentrando la atención en la parte plástica de la energía, se desarrolla a continuación la
expresión que define el trabajo plástico para un material metálico ideal de Prandtl-Reus,
1
& p = λ& (s : s) = ε& P : ε& P =
ω
λ&
( ) &
J2
'P
(
ε& P : ε& P ) (B5.50)
J 2
siendo:
P 1
(
J& 2' = e& P : e& P
2
) 1
e& ijP = ε& ijP − ε& vP δ ij
3 23 (B5.51)
1
0
Tal que en metales puede escribirse el segundo invariante del tensor desviador de
deformaciones como J& 2' =
P 1 P P
2
( )
ε& ij ε& ij , con lo que el trabajo plástico, ecuación (B5.50),
resulta
DINÁMICA NO LINEAL 5-19
& p = 2J 2 ⋅
ω
(ε& P
: ε& P )= 2J 2 ε& P : ε& P (B5.52)
P P
ε& : ε&
La función potencial G que sustituida en (B5.45) da un flujo plástico equivalente a
Prandtl-Reus, es la función de von Mises, G = [G]von Mises ,
∂[G ]
von Mises
g=s= (B5.53)
∂σ
donde,
[G]von Mises = J 2 − K 2 = 0 → [G]
o también
von Mises
= 3J 2 − 3 K = 0 (B5.54)
2
&P =
ω σ ⋅ ε& P : ε& P = σ ε& P (B5.55)
3
Esta expresión permite escribir en forma general la deformación plástica efectiva, como:
( )
t t t
&P
ω 2 &P &P
ε P = ∫ ε& P dt = ∫ dt = ∫ ε : ε dt (B5.57)
0 0
σ 0
3
εP =
3
(
2 P P
ε :ε ) (B5.58)
F(σ; q ) = f (σ ) − K = 0 (B5.59)
Para ello se establece la función de tensión f (σ ) como un traductor de un estado
tensorial de tensiones a otro escalar equivalente. Este escalar se utiliza para ser comparado
con la evolución del endurecimiento plástico K , que inequívocamente se relaciona con la
evolución de la tensión uniaxial equivalente σ = K .
κ p = ε P
( )
K( κ p ) = f κ p con : p (B5.60)
κ = ω p
Definiendo la función de endurecimiento como una variable interna del proceso plástico,
resulta una formulación mucho más general que permite mayores posibilidades de
representación del comportamiento de una gran diversidad de sólidos,
♣
NOTA: f (σ ) es una función homogénea de grado n en las tensiones, siempre que si, f (α, σ ) ≡ α n f (σ ) .
DINÁMICA NO LINEAL 5-21
∂G(σ ; κ p )
κ& p = λ& H κ (σ ; q ) = λ& h κ (σ ; q ) :
∂σ (B5.61)
K = λ H K (σ ; q ) = hK (σ ; q ) κ&
& & p
donde la función tensorial h κ (σ; q ) y la función escalar hK (σ; q ) , dependen del estado de
tensiones actualizado y de las variables internas. En el caso más simple de la teoría de la
plasticidad se identifican las siguientes relaciones,
& p = σ ε& P
h κ ≡ σ ⇒ κ& p = h κ : ε& P ≡ σ : ε& P = ω
de donde resulta (B5.62)
∂K( κ p ) p
K& = hK κ& p = κ&
∂κ p
tal que en esta última K( κ p ) = f ( κ p ) es una función tal como lo expresa la ecuación
(B5.60).
ISÓTROPO CINEMÁTICO
σ2 σ2
F(σ; κ ) = f (σ ) − K( κ ) = 0
p
F(σ , κ, η) = f (σ - η) − K( κ p ) = 0
t+∆t
σ
t+∆t
σ t Sup. actual
t σ
σ
Sup. inicial σ1 σ1
Sup. inicial
Sup. actual
F(σ; q ) = f (σ − η) − K = 0 (B5.63)
donde el endurecimiento plástico puede definirse, según Prager y Melan1, como
η& = β κ& p = ck ε& P , con β = c k ε& P / ε& p . La expresión de c k depende del tipo de función de
potencial plástico que se utilice. En el caso más simple, para una función potencial de von
Mises, tiene la siguiente forma,
5-22 MODELOS INDEPENDIENTES DEL TIEMPO
2
c k = hk (B5.64)
3
donde hk representa una propiedad del material a determinar mediante ensayos.
Utilizando esta misma propiedad, puede definirse la siguiente expresión más general para
5
c k , que ajusta mejor el comportamiento de los metales con endurecimiento cinemático ,
1 σ ij ε ijp
ck = p p ⋅ ⋅ hk (B5.65)
ε& rs ε& rs f (σ kl − η kl )
El origen de la función de carga plástica puede también ajustarse para representar el efecto
Bauschinger, siguiendo una expresión que descomponga el comportamiento cinemático en
la siguiente forma2,6,
F(σ ; q ) = f (σ - η) − K = 0
∂F & ∂F & ∂F & ∂F & ∂F & &
F& = :σ + :η+ K = 0 ⇒ :σ + :η−K = 0 (B5.67)
∂σ ∂η ∂
{K ∂σ ∂η
−1
Sustituyendo en esta la ecuación (B5.61) y la ecuación η& = β κ& p = c k ε& P , resulta
∂F
∂σ
( )
: C : ε& - ε& P + c k
∂F & P
∂η
(
: ε − hK h κ : ε& P = 0 )
(B5.68)
∂F & & ∂F ∂G ∂F ∂G ∂G
∂σ : C : ε - λ ∂σ : C : ∂σ − c k ∂η : ∂σ + hK h κ : ∂σ = 0
De esta última expresión se puede obtener el factor de consistencia plástica λ , que
puede interpretarse como el factor que evalúa la distancia que hay entre un estado tensional
inadmisible, fuera del dominio elástico y la superficie de carga plástica. Esto es
5 Chaboche J. L. (1983). On the constitutive equations of materials under monotonic or cyclic loading. Rech.
Formulation and basic features for ratcheting behavior. International Journal of Plasticity. Vol. 9 pp. 375-390.
DINÁMICA NO LINEAL 5-23
∂F
: C : ε&
λ& = ∂σ donde λ& ≥ 0
∂F ∂G ∂G ∂F ∂G (B5.69)
− c k ∂η : ∂σ + hK h κ : ∂σ + ∂σ : C : ∂σ
444442444443
1
A
donde A es el parámetro de endurecimiento plástico. En un caso simple de la teoría de la
plasticidad clásica sin endurecimiento cinemático ck = 0 , este parámetro resulta ser la
pendiente de la curva tensión uniaxial equivalente σ ( ε p ) = K( ε p ) vs. ε p . Para demostrar
esto, se considera una función de endurecimiento K( ε p ) = f ( ε p ) y se define su pendiente,
dK( ε p ) dK( ε p ) dκ p
A= = (B5.70)
dε p dκ p d ε p
De esta última y la ecuación (B5.62), se verifica el denominador de la ecuación (B5.69)
∂F ∂G
A=− σ = hK h κ : (B5.71)
∂κ p
∂σ
Sustituyendo la ecuación (B5.69) en la ecuación constitutiva tangente,
(
σ& = C : ε& - ε& P ) (B5.72)
Resulta la ley elasto-plástica tangente,
∂G ∂F
C: ⊗ : C
∂σ ∂σ &
σ& = C − :ε ⇒ σ& = C T : ε& (B5.73)
∂F ∂G ∂G ∂ F ∂G
− ck : + hK h κ : + :C:
∂η ∂σ ∂σ ∂σ ∂σ
k −1 [σ ]
i t + ∆t
= C: [ε ] − ε p
t + ∆t i −1 t + ∆t
[ ]
i
k −1 [q]t + ∆t = i −1 [q]t + ∆t
2. Verificación de la condición de fluencia plástica:
i [σ ]t + ∆t = k −1i [σ ] t + ∆t
a. Si: F ( i
[σ ]t + ∆t ; k −1i [q]t + ∆t ) < 0 , entonces i t + ∆t t + ∆t
y va a la SALIDA
[q] = k −1 [q]
k −1 i
b. Si: F ( i
k −1 [σ ]t + ∆t ; k −1i [q]t + ∆t ) ≥ 0 , entonces inicia la integración de la ec. constitutiva
k −1 ∂σ
4. Actualiza las variables internas y el tensor constitutivo tangente con la nueva tensión
k [q ] = k −1 [q]
i t + ∆t i t + ∆t
i t + ∆t
∂G ∂F
C : ∂σ ⊗ ∂σ : C
i
[ ]t + ∆t
= −
k C T C
∂F ∂G ∂G ∂F ∂G
− ck : + hK h κ : + :C:
k
∂η ∂σ ∂σ ∂σ ∂σ
5. Hace k = k + 1 y regresa al punto 2
( )
σ : ε& p ≥ σ ∗ : ε& p → σ − σ * : ε& P ≥ 0 (B5.74)
donde puede verse que necesariamente se exige que el estado tensional posterior σ , sea
siempre mayor que el anterior σ ∗ .
Haciendo ahora la siguiente aproximación,
y sustituyendo esta última en la ecuación (B5.76), resulta la siguiente forma particular del
2ndo postulado de Drucker
∂Ξ
ε& : e ≥ 0 (B5.77)
∂ε
donde la disipación plástica Ξ , para problemas sin degradación de rigidez, se escribe como,
& ≥0
Ξ = σ : ε& P − Ψ (B5.78)
Sustituyendo ésta última en la expresión de la M.D.P., se tiene
∂G
ε& : (C ) : ε& p = ε& : (C ) : λ& ≥0 (B5.80)
∂σ
Pero el factor de consistencia plástica λ es un escalar no negativo. Por ello la desigualdad
anterior puede también escribirse, como
∂G
λ& ≥ 0 , ε& : (C ) : ≥0 (B5.81)
∂σ
Además, analizando el propio factor de consistencia plástica, se deduce que,
∂F
: C : ε&
&
λ= ∂σ ≥0 (B5.82)
∂G ∂F
A+ :C:
∂σ ∂σ
si el endurecimiento es positivo A > 0 , y G y F son funciones convexas, para que el factor
de consistencia plástico sea negativo, debe cumplirse necesariamente que
∂G
ε& : C : ≥0 (B5.83)
∂σ
Para garantizar que las desigualdades (B5.81) y (B5.83) se cumplan a la vez, debe ocurrir
que el flujo plástico sea asociado. Es decir,
∂G ∂F
∝ , Flujo asociado (B5.84)
∂σ ∂σ
1 2 1
Π = Π * + δΠ + δ Π + L ⇒ ∆Π = Π − Π * ≅ δ{
Π + δ2Π + L (B5.85)
2! 0 2!
DINÁMICA NO LINEAL 5-27
siempre es posible conocer la expresión del funcional de energía potencial total7, si que es
posible conocer sus respectivas variaciones, por ejemplo mediante el principio de los
trabajos virtuales que es lo que se ha establecido en las ecuaciones anteriores. Sustituyendo
en (B5.85) la condición de equilibrio, se tiene que el incremento total de trabajo virtual es
igual a la segunda variación del funcional y que es a la vez la condición de estabilidad
(concavidad o convexidad del funcional),
1 1 1
∆Π ≅ ∫ δσ : δε dV = ∫ δσ : δε dV + ∫ δσ : δε dV ≅ ∆Π + ∆Π >0 (B5.88)
2 2 2 v0 vp
v v v
0 p
De todo esto se puede concluir que, en un punto del sólido es posible que se viole la
condición de estabilidad de Drucker (ecuación (B5.83)), sin que necesariamente la respuesta
global sea inestable.
λ& ≥ 0
F(σ ; q ) ≤ 0 (B5.90)
&
λ F(σ ; q ) = 0
de estas tres condiciones se deduce brevemente lo siguiente,
F < 0 ⇒ λ& = 0 comportamiento elástico o decarga,
λ& > 0 comportamiento plástico o carga,
F = 0 ⇒ & (B5.91)
λ = 0 carga plástica neutra,
F > 0 ⇒ estado incompati ble.
DINÁMICA NO LINEAL 5-29
( ) π
F I 1 ; J 2 ; θ; σ Tmax = 2 3 J 2 cos θ + + I 1 − 3 σ Tmax ( κ) = 0
6
(B5.93)
( ) π
F ρ; ξ; θ; σ Tmax = 2 ρ cos θ + + ξ − 3 σ Tmax ( κ) = 0
6
(B5.94)
2 J2
donde ξ = 3 σ oct = 3 ( I 1 3) = I 1 3 , el radio octaédrico ρ = 3 τ oct = 3 = 2 J2 ,
3
[( )
y el ángulo de similaridad de Lode θ = arcsen 3 3 J 3 (2 J 23 / 2 ) . ]
DINÁMICA NO LINEAL 5-31
− σ1
a) plano octaédrico
I 1 = 0 (plano Π )
− σ2
− σI
− σ II
σ I = σ II = σ III π
−
6
corte puro
π
+
6 − σ3
− σ III
b)
meridiano de
tracción − π 3 max
6 ρ t0 = σt
2
φT ξ 0 = 3σ tmax
ξ
φC
ρc0 = + 6σtmax
meridiano de
compresión + π
6
− σ II +
5π
− σI
6
ρ 0c ρ t0
c)
π
−
6
− σ III θ=0
( ) 1
F σ ; τ max = max σ i −σ j
max
− τ ( κ) = 0 (B5.95)
2
(
F J 2 ; θ; τ max =) J2
3
cos θ − τ max ( κ) = 0
(B5.96)
o multiplicando por 2 2 , resulta en función de la tensión uniaxial efectiva,
F(J 2 ; θ; σ ) = 2 J 2 cos θ − σ ( κ) = 0
[( ) ]
siendo ρ = 3 τoct = 2 J 2 , y el ángulo de similaridad de Lode θ = arcsen 3 3 J 3 (2 J 23 / 2 ) .
La insensibilidad a la presión hace que el plano octaédrico se mantenga constante e igual
al plano π. Este plano octaédrico representa un hexágono regular. De la intersección del
plano meridiano de tracción (θ = −π / 6) , con la superficie de fluencia, surge una recta de
pendiente nula, paralela a la que resulta de la intersección del plano meridiano de
compresión (θ = + π / 6) , con la superficie de fluencia. Ambas rectas meridianas cortan al eje
de tensión de corte octaédrico en ρ C0 = ρ T0 = ± 2 / 3 σ ( κ) . En el plano principal
σ1 , σ 3 , σ 2 = 0 , representa un hexágono deformado según el eje de tensiones σ1 = σ 3 .
DINÁMICA NO LINEAL 5-33
− σ1
plano octaédrico
I 1 = 0 (plano Π )
a) Espacio de Westergaard
− σ2
− σI
− σ II
ξ
σ I = σ II = σ III ρ
π
−
6
corte puro
π
+ − σ III − σ3
6
b) Planos Meridianos
π Meridiano de tracción 2
θ=− ρ t0 = σ(x )
6 3
ξ
1
ρ t0 = − σ( x )
θ = 0 Meridiano de corte puro 2
π Meridiano de compressión ρ t0 =
2
σ(x )
θ=− 3
6
c) Plano Octaédrico
π
θ=5
6
− σ II
− σI
ρ t0
ρ 0c
π
θ=−
6
π
θ=+
− σ III θ = 0 corte puro
( max
F σ ; τ oct) 16 [(σ −σ
= 1 2 )2 + (σ 2 −σ 3 )2 + (σ 3 −σ1 )2 ] − [τ oct
max
( κ) ]
2
=0 (5.98)
siendo ρ = 3 τoct = 2 J 2 .
La insensibilidad a la presión hace que el plano octaédrico se mantenga constante e igual
al plano Π . Este plano octaédrico representa un círculo. De la intersección del plano
meridiano de tracción (θ = −π / 6) , con la superficie de fluencia, surge una recta de
pendiente nula, paralela a la que resulta de la intersección del plano meridiano de
compresión (θ = + π / 6) , con la superficie de fluencia. Ambas rectas meridianas cortan al eje
de tensión de corte octaédrico en ρ C0 = ρ T0 = ± 2 / 3 σ ( κ) al igual que el criterio de Tresca.
En el plano principal σ1 , σ 3 , σ 2 = 0 , representa una elipse cuyo eje mayor coincide con el
eje de tensiones σ1 = σ 3 .
DINÁMICA NO LINEAL 5-35
− σ1
a)
− σI
− σ2
O ξ
π
−
6
σ1 = σ2 = σ3
− σ3
π
+
6
− σ III
ρ
b)
meridiano de tracción
π 2
θ=− ρ t0 = σ(κ ) ≡ ρ 0τ
6 3
ξ
π 2
θ=+ ρ 0c = σ(κ ) ≡ ρ 0τ
6 3
meridiano de compresión
− σI
− σ II
c)
O
2 π
ρ 0c = ρt0 = ρ 0τ = σ θ=−
3 ρ 6
π
θ=+
6
− σ III
τ
σ1 > σ 3 > σ 3
τ = +c
φ
σ1 − σ3 φ
τ
2
− σn φ σ
− σ3 − σ1
A O σ on = c.cotgφ
τ = −c
σ + σ3
− 1
2
Figura B5.14 – Representación plana del estado tensional en un punto, según el criterio
de Mohr-Coulomb.
Observando la Figura B5.14, se puede rescribir la ecuación (B5.101) en función de las
tensiones principales,
σ − σ 3 σ + σ 3 σ1 − σ 3
F(σ ; c; φ) = 1 cos φ − c − − 1 − sen φ = 0
2 2 2
(B5.102)
σ − σ 3 σ1 + σ 3
⇒ F(σ ; c; φ) = 1 + sen φ − c φ cos = 0
2 2
donde σ1 y σ 3 representan la tensión principal mayor y menor respectivamente. De esta se
deduce que el criterio de Mohr-Coulomb ignora el efecto de la tensión principal intermedia
σ 2 , lo que es una gran limitación. No obstante, este problema se soluciona si se expresa su
formulación en función de los invariantes, del tensor de tensiones y sus desviadores,
DINÁMICA NO LINEAL 5-37
I1 sen θ sen φ
F(I 1 ; J 2 ; θ; c; φ) = sen φ + J 2 cos θ − − 6 c( κ ) cos φ = 0 (B5.103)
3 3
sen θ sen φ
F(ρ; ξ; θ; c; φ) = 2 ξ sen φ + 3 ρ cos θ − − 6 c( κ ) cos φ = 0 (B5.104)
3
[(
de similaridad de Lode θ = arcsen 3 3 J 3 (2 J 23 / 2 ) . ) ]
Estas funciones describen en el espacio de tensiones principales una pirámide de base
hexagonal distorsionada, cuyo eje coincide con el de presiones isostáticas σ1 = σ 2 = σ 3 . El
aumento de presión hace que el plano octaédrico crezca. Este plano octaédrico representa
un hexágono deformado. De la intersección del plano meridiano de tracción (θ = −π / 6) ,
con la superficie de fluencia, surge una recta de pendiente ( 2 2 sen φ ) /( 3 + sen φ ) ,
que corta el eje de tensión tangencial octaédrica en ρ T0 = ( 2 c 6 cosφ) /(3 + senφ) y el eje de
presiones en ξ 0 = 3 c cotgφ . De la intersección del plano meridiano de compresión
(θ = + π / 6) , con la superficie de fluencia, resulta una recta de pendiente
( 2 2 senφ) /(3 − senφ) , mayor que la del meridiano de tracción y que corta el eje de tensión
tangencial octaédrica en ρ C0 = ( 2 c 6 cosφ) /(3 − senφ) . En el plano principal σ1 , σ 3 , σ 2 = 0 ,
representa una hexágono deformado cuyo eje mayor coincide con el eje de tensiones
σ1 = σ 3 .
De las funciones que describen el criterio de fluencia de Mohr-Coulomb, resulta claro
que su principal característica es la capacidad para distinguir el comportamiento a tracción
del de compresión. De aquí resulta que el criterio admite implícitamente que la relación de
resistencia a tracción y compresión, cumplen con la siguiente expresión (ver Figura B5.15),
σ C0 π φ
R Mohr = = tan 2 + (B5.105)
σ 0
T 4 2
Mohr-Coulomb Mohr-Coulomb
Modificado Standard
σ oc
Ro =
σ To π φ
R = tan 2 +
2 π φ 4 2
R = α tan +
4 2
Equivalente a α = 1,0
14 α = 3,61 (60 o
, ≈ 14 )
α = 2,16 σ II α = α1
Rango del hormigón
10 α = α2
σI
8
α 2 > α1
6 (45 , ≈ 6)
o
3
(30 ,3)
o
Dominio del
Hormigón
− σ1
plano octaédrico
I 1 = 0 (plano Π )
− σ2
− σI 3c cot gφ
a)
− σ II
ξ
σ I = σ II = σ III ρ
π
−
6
corte puro
π − σ3
+ − σ III
6
ρ+ξ
ρ−ξ
ρ = 3τ oct
TRACCIÓN
meridiano de TOTAL
tracción − π 2c 6 cos φ
6 ρt0 =
(3 − senφ )
COMPRESIÓN
TOTAL
45 0 φT
ξ 0 = 3c cot g φ b)
ξ = 3σ oct
φC
dε
φC
ρc0 3 + senφ
R0 = =
ρto 3 − senφ
2c 4 cos φ componente que
ρ0c =
(3 − senφ) provoca dilatancia
meridiano de
compresión + π
6 σ3
I1 = 0 φ=0
− σ II − σI
σ1 ≤ 0 ≤ σ 3 0
2c R Mohr
0
2c R Mohr
σ1
ρ t0
σ1 ≤ σ 3 ≤ 0 2c
ρ 0c
σ 3 ≤ 0 ≤ σ1
θ=0 σ 3 ≤ σ1 ≤ 0
c) − σ III d) 0 π φ σ
o
RMohr = tan 2 + = 0c
4 2 σT
Figura B5.16 – Superficie de fluencia de Mohr-Coulomb: a) En el espacio de tensiones
principales, b) Según los meridianos de tracción y compresión Máxima, c) Según el plano
octaédrico I1=0 o plano Π .
5-40 MODELOS INDEPENDIENTES DEL TIEMPO
− σ1
a) plano octaédrico
I 1 = 0 (plano Π )
− σ2
− σI
K
3α
− σ II
ξ
ρ π
−
6
b) ρ+ξ = 0 ρ−ξ = 0
ρ = 3τ oct
ξ = 3σ oct
ρc0
R0 = 1 =
ρto
ρ t0 = 2K
componente que
provoca dilatancia
meridiano de
compresión + π
6
− σ II − σI
O
2 π
ρ 0c = ρt0 = ρ 0τ = σ θ=−
3 ρ 6
c) θ=+
π
6
− σ III
F
F
Material
Friccional
F F
PLASTICIDAD
MODELO DE DAÑO
DAÑO
MATERIALES PLÁSTICO
DÚCTILES MATERIALES
FRÁGILES
Figura B5.19 – Representación simple de las teorías que contribuyen a la definición del
“modelo de daño plástico”.
Luccioni, B. (1993). Formulación de un Modelo Constitutivo para Materiales Ortótropos, Tesis Doctoral.
11
12 Lubliner, J.; Oliver J.; Oller S. and Oñate E. (1989). A Plastic-Damage Model for Concrete. International
Journal of Solids and Structures, Vol.25, No.3, pp. 299,326.
13 Oller S. (2001). Fracura mecánica – Un enfoque global. CIMNE – Ediciones UPC. Barcelona
♣
Nota: Se dice que el comportamiento de un punto del sólido es “radial” o “proporcional”, si durante todo
0 0
el proceso de carga se cumple la relación σ11 σ11 = σ12 σ12 = L = σ 33 σ 033 , entre las componentes del tensor
de tensiones en el estado actual y el estado inicial, respectivamente.
5-44 MODELOS INDEPENDIENTES DEL TIEMPO
Puntos
dañados Lugar geométrico
del daño
MECÁNICA DEL
SÓLIDO
TEORÍA
DE LA TÉCNICAS NUMÉRICAS
PLASTICIDAD
Y DAÑO
• MEF
• Resolución del Sistema de Ecuaciones
• Control de Desplazamiento
• Control de Plastificación
Formulación • Integración de la Ecuación Constitutiva
MODELO Implementación
DE DAÑO
PLÁSTICO
14 Malvern, L. (1969). Introduction to the mechanics of continuous medium. Prentice Hall, Englewood Cliffs, NJ.
DINÁMICA NO LINEAL 5-45
Como puede observarse, todos los ítem aquí nominados no son simples, pero si posibles,
de conseguir dentro de una formulación mecánica.
masa del sólido en el punto de análisis (discontinuidad física). A continuación se hace una
breve presentación, que para mayor detalle pueden consultarse las fuentes1,15,16.
El criterio de fluencia plástica presentado en apartados anteriores, es ahora tratado en
este modelo a través de una expresión matemática que puede ser escrita en la siguiente
forma general,
F(σ , c ) = f (σ ) − c = 0 (B5.108)
donde f (σ ) es una función escalar homogénea de primer grado en las componentes del
tensor de tensiones, que permite definir la cohesión c , o una tensión uniaxial escalada,
como una función de endurecimiento plástico o como una variable interna dependiente de
la evolución del proceso mecánico.
Los criterios de plasticidad de Mohr-Coulomb y Drucker-Prager, también pueden
representarse por una expresión similar a la ecuación (B5.108), pero no aproximan
adecuadamente el comportamiento real de muchos de los materiales friccionales.
Numerosos criterios de fluencia plástica han tratado de mejorar esta aproximación para los
materiales cohesivos-friccionales. Hay algunas funciones f (σ ) no-homogéneas y de grado
superior en las componentes del tensor de tensiones, característica que impide definir una
función de endurecimiento plástico c con una interpretación física directa.
El modelo constitutivo de daño plástico puede utilizarse cualquier criterio de fluencia,
pero para mejorar la aproximación del comportamiento particular del hormigón se pueden
definir criterios apropiados a cada caso1.
La cohesión c es tratada como una magnitud escalada con la resistencia inicial a
compresión uniaxial del hormigón σ C0 (umbral de discontinuidad tensional), que es el nivel
de tensiones para el cual la deformación volumétrica εV es máxima. Consecuentemente se
define la cohesión inicial, o cohesión del material virgen, como c 0 ∝ σ C0 para κ p = 0 ,
situación que establece la posición inicial del criterio de fluencia, y la cohesión final, o
cohesión del material totalmente deteriorado, c u = 0 para κ p = 1 , situación que define la
posición final del criterio de fluencia.
A diferencia de la plasticidad clásica con endurecimiento isótropo, la cohesión no es una
simple función de la variable de endurecimiento plástico c( κ p ) , sino una variable interna que
depende de la evolución del proceso elasto-plástico, gobernada por una ecuación de
evolución (ecuación diferencial).
El ángulo de rozamiento interno φ también podría definirse como una variable
interna a partir de una ley de evolución que dependa del proceso elasto-plástico, pero dada
la evidencia del comportamiento físico de este fenómeno en hormigones17, sólo se plantea
una simple función explícita de la variable de daño plástico φ( κ p ) . Con esta hipótesis se
obtiene una fricción inicial nula φ0 = 0 cuando aún la cohesión c 0 no permite la
movilización de la fricción, y máxima al final del proceso elasto-plástico
15 Lubliner, J.; Oliver J.; Oller S. and Oñate E. (1989). A Plastic-Damage Model for Concrete. International
Structures Using a ``Plastic-Damage Model''. Engineering Fracture Mechanics, Vol 35; pp 219-231.
17 Borst, R. De and Vermeer, P. (1984). Non associated plasticity for soils, concrete and rock. Heron. Vol. 29, Delf.
Netherlands.
DINÁMICA NO LINEAL 5-47
El modelo constitutivo que resulta de estas definiciones básicas, consigue una muy
buena respuesta durante el proceso de fractura. A modo de resumen puede decirse que el
modelo reúne las siguientes características,
• Definición de una ley constitutiva que depende de las variables internas de cohesión
y daño plástico, permitiendo así representar situaciones de carga complejas no-
radiales.
• Trata en forma unificada los estados complejos de tensión multiaxial.
• Admite que los materiales tiene distintos limites de resistencia máxima y de
deformación última, dependiendo del proceso mecánica que este desarrollándose.
• Admite la posibilidad de considerar distintos criterios de fluencia plástica, no siendo
esta una característica del modelo, sino una variable más del mismo que necesita ser
preestablecida.
• Considera un flujo plástico no-asociado, que permite el control del fenómeno de
dilatancia.
• Permite obtener toda la información relacionada con el deterioro de un punto a
través de un post-proceso de la información mecánica del punto.
plástico que contemple la posibilidad de tener una deformación última diferenciada para
cada proceso en particular. Esto se consigue estableciendo una variable interna que define
una disipación normalizada a la unidad, es decir que se trata de la relación que hay entre la
densidad de energía disipada en un determinado instante del proceso respecto del máximo
que podría disipar el punto del sólido. Así, puede decirse que la variable de daño plástico es una
“medida relativa” de la energía disipada durante el proceso plástico. Por simplicidad y
orden en la presentación, es conveniente definir esta variable para procesos uniaxiales
simples y luego generalizarla a procesos multiaxiales.
1 t
0 ≤ κ p = p ∫t =0 σT ε& T dt ≤ 1
p
(B5.115)
gf
Resultando así una variable normalizada con respecto a la energía específica máxima a
tracción, con valores comprendido entre 0 ≤ κ p ≤ 1 para el inicio y fin del proceso plástico,
respectivamente. Con κ p como variable puede ahora transformarse las curvas de respuesta
uniaxial de tensión-deformación plástica en otras curvas que dependan de esta variable de
daño plástico, en la cual que se cumplen los siguientes extremos, σ T ( κ p = 0) = σ T0 y
σ T ( κ p = 1) = σUltima
T en el intervalo [0,1].
En el problema de compresión uniaxial, al igual que en el problema de tracción, la
variable de daño plástico resulta
1 t
0 ≤ κ p = p ∫t =0 σC ε& C dt ≤ 1
p
(B5.116)
gC
∞
siendo g Cp = ∫t =0 σ C ε& Cp dt la energía disipada en el proceso de compresión. Al igual que en
el proceso de tracción pura, puede ahora definirse la respuesta uniaxial de tensión-
deformación plástica mediante otra curva que dependa del daño plástico en lugar de la
deformación plástica, en la cual que se cumplen los siguientes extremos, σ C ( κ p = 0) = σ C0 y
σ C ( κ p = 1) = σUltima
C en el intervalo [0,1] (ver Figura B5.22 ).
5-50 MODELOS INDEPENDIENTES DEL TIEMPO
σ σ
pic
σC σC( ε p ) σC pic
c
c Cpic
c 0 = f (σ ) 0
cC(κp)
curva de compresión
cT(κp)
curva de tracción
1 κp
c) Cohesión en función de la variable de daño plástica
r (σ ) 1 − r (σ )
κ& p = h κ : ε& p = p + ⋅ Ξm (B5.118)
g f g Cp
escalar que define los estado de comportamiento de un punto en función del estado
tensional, siendo x = 0,5 [x + x ] la función de McAully. Obsérvese los siguientes casos
particulares, r (σ ) = 1 para problemas de tracción pura, r (σ ) = 0 para compresión pura y
r (σ ) = 0,5 para estado de corte puro. De esta manera, la disipación plástica siempre estará
normalizada respecto de la máxima energía correspondiente al proceso que está realizando
en cada momento.
donde hc ( σ , κ p , c ) es una función escalar del estado actual de la variable libre de tensión σ
y de las variables internas κ p y c . La expresión adoptada para la ley de evolución de la
variable interna de cohesión resulta de elegir la siguiente expresión para hc ,
r ( σ ) dc T 1 − r ( σ ) dc C
hc = c ⋅ p
+ (B5.120)
cT d κ cC d κp
1
c C (κ p ) = σC (κ p ) (B5.121)
ℵ
Tal que ℵ es un coeficiente que depende del tipo de criterio de umbral de fluencia y
representa el factor de escala entre cohesión y tensión uniaxial de compresión1. A modo de
ejemplo para Tresca y von-Mises vale ℵ = 1 , para Mohr Coulomb ℵ = 2 RMohr donde
R Mohr es la relación de resistencias (ver ecuación 6.111), para Drucker-Prager inscrita en la
superficie de Mohr-Coulomb ℵ = 6 cos φ /(sin φ − 3) y para Drucker-Prager circunscrita en la
superficie de Mohr-Coulomb ℵ = 6 cos φ /(3 sin φ − 3) . Así, para cualquier criterio umbral de
fluencia, debe definirse este coeficiente.
La función c T ( κ p ) (ver Figura B5.22) puede obtenerse en forma explícita y representa la
evolución de la cohesión durante un ensayo uniaxial de tracción simple. La relación entre
cohesión y tensión uniaxial de tracción viene dada por la siguiente expresión
R0
c T (κ p ) = σT ( κ p ) (B5.122)
ℵ
[ ] [ ]
donde R 0 = f C0 f T0 = σ C ( κ p = 0) σ T ( κ p = 0) es la relación entre resistencias uniaxiales.
A modo de ejemplo para Tresca y von-Mises vale R 0 / ℵ = 1 , para Mohr Coulomb
R 0 / ℵ = RMohr / 2 donde R Mohr es la relación de resistencias (ver ecuación 6.111), para
Drucker-Prager inscrita en la superficie de Mohr-Coulomb R 0 / ℵ = (3 + 3 sin φ) / 6 cos φ y
para Drucker-Prager circunscrita en la superficie de Mohr-Coulomb
R / ℵ = (3 + sin φ) / 6 cos φ .
Algunos investigadores sostienen que las curvas de resistencia a tracción y compresión
simple del hormigón, que resultan de ensayos experimentales uniaxiales, tienen formas
análogas18,19, lo que equivale a afirmar que la relación de escalas entre ellas es una constante
durante todo el proceso de carga cuasi-estático, y viene dado por
σC (κ p )
R( κ p ) = = cte. ⇒ R ( κ p ) = R 0 (B5.123)
σT ( κ p )
y en tal las funciones explícitas de la cohesión a tracción uniaxial y a compresión uniaxial
coinciden.
18 State of the art report on : Finite Element Analysis of Reinforced Concrete. ASCE (1982).
19 Tasuji, E.; Slate F. and Nilson A. (1978). Stress strain response and fracture of concrete. Journal of the
structural division. ASCE – Vol 75, No 7, pp. 306-312.
DINÁMICA NO LINEAL 5-53
senφ
senφ max
κp
L
c κ Zona donde la fricción está
completamente movilizada y
c pic la cohesión es nula
κp
κL
Figura B5.23 – Función que define la evolución del ángulo de rozamiento interno en
función de la variable de daño plástico y su relación con la cohesión.
En forma general se podría formular una variable interna de fricción, con una ley de
evolución del tipo φ& = λ& H φ (σ ; q ) = hφ (σ ; q )⋅ κ& p , que represente el mecanismo de
incremento de la fricción arriba mencionado. Sin embargo, hay evidencias
experimentales1,20 que muestran que es suficiente proponer una simple función de la
variable de daño plástico para representar la evolución del ángulo de fricción interna,
κ pκL
p 2 p sen φ max ; ∀ κ p ≤ κL
sen φ( κ ) = κ + κ L (B5.124)
max p L
sen φ ;∀ κ >κ
Borst, R. De and Vermeer, P. (1984). Non associated plasticity for soils, concrete and rock. Heron. Vol. 29, Delf.
20
Netherlands.
5-54 MODELOS INDEPENDIENTES DEL TIEMPO
π
meridiano de tracción −
6
c ( κ 2p ) ρ = 3 τ OCT = 2 J 2 ρ−ξ
c( κ3p ) c( κ1p )
TRACCIÓN
TOTAL
2 c( κ p ) 6 cos φ
ρt0 =
(3 − sen φ)
COMPRESIÓN
TOTAL ξ 0 = 3 c ( κ p ) cotg φ
a)Evolución de la cohesión 450 φT
c( κ1p ) > c( κ 2p ) > c( κ 3p ) φC
3
ξ = 3 σ OCT = I1
3
ρ 0c 3 + sen φ
R0 = =
ρ to 3 − sen φ 2 c( κ p ) 4 cos φ
ρ C0 =
(3 − sen φ)
π
meridiano de compresión +
6
π
meridiano de tracción −
6 ρ+ξ ρ = 3 τ OCT = 2 J 2
ρ−ξ
φ(κp)Max
TRACCIÓN
TOTAL
2c 6 cos φ( κ p )
ρ t0 =
COMPRESIÓN
TOTAL (3 − sen φ(κ ))
p
π 2c 4 cos φ( κ p )
meridiano de compresión + ρ C0 =
6 (3 − sen φ(κ ))
p
Figura B5.24 – Movilidad del criterio de fluencia de Mohr-Coulomb según sus planos
meridianos. a) Movilidad isótropa por evolución de la cohesión. b) Movilidad anisótropa
por cambio en el rozamiento interno.
tipo de la de Mohr-Coulomb, en la cual puede verse una movilidad isótropa por efecto de
la evolución de la cohesión y otra movilidad anisótropa provocada por el incremento del
rozamiento interno. Este último fenómeno provoca un movimiento de acercamiento de la
cúspide de la superficie hacia el origen, acompañada de un crecimiento de su base (ver
Figura B5.24 y Figura B5.25). Este efecto se presenta como un fenómeno de
endurecimiento para procesos de compresión y como un ablandamiento para procesos de
tracción.
− σ1
− σI
− σ2
3 c cotgφ( κ p )
− σ II
σ I = σ II = σ III ξ
π
−
6
corte puro
π − σ3
+
6
− σ III
σ3
I1 = 0 φ=0
σ1 ≤ 0 ≤ σ 3 0
2c RMohr
0
2c R Mohr σ1
2c
σ 3 ≤ 0 ≤ σ1
o
0 π φ σ
σ 3 ≤ σ1 ≤ 0 R Mohr = tan 2 + = 0c
4 2 σT
En materiales frágiles con alta cohesión inicial, como los hormigones, cerámicos, etc., es
posible utilizar durante todo el proceso plástico un ángulo de rozamiento interno constante
y máximo sin que esto induzca a resultados insatisfactorios en la resolución de problemas
multiaxiales.
5-56 MODELOS INDEPENDIENTES DEL TIEMPO
A τ
p
(ε )ρ A τ
(εp)ξ A
εp
ψ
A
τ
(ε ) p
ξ : Parte volumétrica de la def. plástica (dilatancia)
(ε ) p
ρ : Parte desviadora de la def. plástica
ρ+ξ (ε p ) ρ εp ρ−ξ
meridiano de ψ ρ = 3τ oct = 2 J 2
π
tracción −
6 (ε p ) ξ
TRACCIÓN
TOTAL
COMPRESIÓN 2 c 6 cos φ
ρ ρ t0 =
TOTAL (3 − sen φ)
ξ 0 = 3 c cotgφ
ξ 45 0
φT
φC I1
ρ 0c 3 + sen φ ξ = 3σ oct =
R0 = = 3
ρ to 3 − sen φ
G(σ ) = cte
2 c 4 cos φ
ρ 0c =
(3 − sen φ)
meridiano de
π
compresión +
6
Figura B5.26 – Representación de la dilatancia en forma esquemática y su relación con la
deformación plástica en el plano meridiano.
Una forma apropiada de evaluar este fenómeno, es mediante el ángulo de dilatancia ψ , que
fue inicialmente introducido por B. Hansen21 y que representa la relación que hay entre el
incremento de volumen plástico y la distorsión plástica.
21Hansen, B. (1958). Line ruptures regarded as narrow ruptures zone – Basic equations based on kinematics
considerations. Proc. Brussels Conf. 58 on earth pressures problems.
DINÁMICA NO LINEAL 5-57
siendo φ max la máxima fricción y ψ max es la máxima dilatancia. Como valores orientativos,
para el hormigón estas magnitudes valen φ max = 35 0 y ψ max = 13 0 .
B 5.15.1 Introducción.
Los resultados experimentales muestran que los materiales cohesivos-friccionales tienen
una pérdida de rigidez aún en el campo de comportamiento elástico. Este efecto se
incrementa mucho más cuando se produce la descohesión entre partículas y por lo tanto
cuando comienza el período plástico. Esta evidencia induce a pensar que hay dos
fenómenos de degradación de rigidez que actúan sobre el material. Uno que depende de la
energía acumulada, que se denomina degradación elástica, y otro que depende de la
movilización de la fricción y que recibe el nombre de degradación plástica. Basado en estos
dos fenómenos, se modificará en este apartado el modelo de daño plástico antes presentado.
Rowe, P. (1972). Theoretical meaning and observed values of deformation parameter for soil. Proc. Rascoe
22
d& ie = Φ i k ie : ε& e
(B5.127)
d& ip = λ& H i = k ip : ε& p
σ0
E0
E0
ε pic εu ε
Figura B5.27 – Evolución de la resistencia uniaxial en un punto por efecto de la plasticidad
y la degradación de rigidez.
donde “i” representa el índice del mecanismo de daño que se está desarrollando y el
problema puede estar compuesto de un número finito de mecanismos distintos. El tensor
k ie representa la dirección de degradación elástica del mecanismo i esimo , Φ i un escalar
positivo a definir y H i = k ip : (∂G / ∂σ ) una función escalar elasto-plástica de argumentos
tensoriales que tiene en cuenta la dirección de degradación k ip inducida por la plasticidad.
23
Kachanov,L.M. (1958).Time Rupture Process under Creep Conditions, (in Russian). Izv.ARad.SSSR
Teckh.Nauk.,8 ,26-31.
24 Kachanov,L.M. (1986 ). Introduction to Continuum Damage Mechanics. Martinus Nijho Publishers, Dordrecht.
DINÁMICA NO LINEAL 5-59
Ψ( ε e , q α , q β ) = Ψe (ε e , d p , q β ) + Ψp (q α ) (B5.128)
& ≥0
Ξ = σ : ε& − Ψ (B5.129)
Esta desigualdad expresa el balance de entropía para el continuo de Cauchy y es válida
para cualquier proceso de carga admisible. Sustituyendo la derivada temporal de la energía
libre (B5.128) en la ecuación (B5.129) resulta la siguiente expresión para la disipación,
∂Ψ & ∂Ψ & p ∂Ψ ∂Ψ
Ξ = σ - e
:ε+ :ε − q& α − q& β ≥ 0 (B5.130)
∂ε ∂ε e
∂q α ∂q β
∂Ψ
σ − ≥0 ∀ ε& (B5.131)
∂ε e
de donde se deduce que la tensión vale,
∂Ψ
σ= (B5.132)
∂ε e
tal que la energía libre para un sólido elasto-plástico con degradación de rigidez puede ser
escrita en pequeñas deformaciones, como
1
Ψ(ε e , q α , q β ) = (ε − ε p ) : C(d ie , d ip ) : (ε − ε p ) + Ψ p (q α ) (B5.133)
2
Sustituyendo esta última en la ecuación (B5.132) queda expresada la ecuación constitutiva
secante para un sólido elasto-plástico con degradación de rigidez, como
σ = C( d ie , d ip ) : (ε − ε p ) (B5.134)
A partir de esta ecuación se obtiene la variación temporal de la tensión,
5-60 MODELOS INDEPENDIENTES DEL TIEMPO
∂C ∂C p
σ& = ∑ e d ie + p
d i : (ε - ε p ) + C : (ε& - ε& p ) (B5.135)
i ∂
i
d ∂d i
Se considera ahora un caso particular de daño simple, con un único mecanismo de daño
(i=1), entonces las reglas de evolución (B5.127) pueden simplificarse de la siguiente forma,
e e p Φ
C T = C(d , d ) − 1 − d e (σ ⊗ σ )
σ& = C Te (d e , d p ) : ε& - C Tp (d e , d p ) : ε& p (B5.137)
p e p hc
C T = C(d , d ) − c (σ ⊗ h κ )
d e = Φ ωe (B5.138)
La magnitud de Φ resulta de un proceso de carga en el que se congela la variable de
degradación plástica y se deja que sólo evolucione la degradación elástica. Esto es,
σ = C( d e ) : ε e = (1 − Φ ω e ) ⋅ C 0 : ε e ⇒ ω [
& e = ε e : (1 − Φ ω e ) C 0 : ε& e ] (B5.139)
de donde puede obtenerse,
[ ]
t t
&e
ω e 0 e
- (Φ ε : C : ε ) / 2
∫ (1 − Φ ωe ) ∫ ε : C : ε& dt ⇒
e 0 e
dt = d e = Φ ωe = 1 − e (B5.140)
t =0 t =0
1 - (Φ ε :C :ε )/2
e 0 e
1 e
ω d = ∫ Ξ e dt = 1− e 0
− 2 ε :C :ε
e
(B5.142)
t
Φ
DINÁMICA NO LINEAL 5-61
Esta degradación elástica puede también ajustarse más a la realidad de algunos procesos
mecánicos, estableciendo una cuota de degradación para comportamientos volumétricos y
otra distinta para el comportamiento desviador. Ejemplos de una formulación de este tipo
pueden verse en la referencia1.
Energía disipada ωd
E0 ε ε
F(σ, c ) = f (σ ) − c = 0
∂F & ∂F & ∂F & &
F& = :σ + c=0 ⇒ :σ − c = 0 (B5.143)
∂σ ∂c ∂σ
{
−1
∂F
∂σ
[ ] ( )
: C Te (d e , d p ) : ε& − C Tp (d e , d p ) : ε& p − hc h κ : ε& P = 0
(B5.144)
∂F e & & ∂F p ∂G ∂G
∂σ : C T : ε − λ ⋅ ∂σ : C T : ∂σ + hc h κ : ∂σ = 0
Resultando de esta última expresión el factor de consistencia plástica λ, que representa
una medida de la distancia que hay entre un estado tensional inadmisible, fuera del dominio
elástico, y la superficie de carga plástica. Esto es
5-62 MODELOS INDEPENDIENTES DEL TIEMPO
∂F
: C Te : ε&
λ& = ∂σ con λ& ≥ 0
∂G ∂F p ∂G (B5.145)
hc h κ : ∂σ + ∂σ : C T : ∂σ
42
1 4 44 3
A
donde A es el parámetro de endurecimiento plástico. Sustituyendo la ecuación (B5.69) en la
ecuación constitutiva tangente (B5.137), se tiene la siguiente ley elasto-plástica tangente con
daño,
p ∂G ∂F
CT : ⊗ : C Te
∂σ ∂σ : ε&
σ& = C Te − ⇒ σ& = C Tep : ε&
∂G ∂F p ∂G (B5.146)
hc h κ : + : CT :
∂σ ∂σ ∂σ
donde C Tep es el tensor constitutivo tangente continuo. Como puede observarse, se tendrá
rigidez tangente simétrica si se cumple la siguiente proporcionalidad,
p ∂G ∂F e
C T : ∂σ ∝ ∂σ : C T . Con lo que puede verse que no es suficiente que se cumpla la
clásica regla de flujo asociada para garantizar dicha simetría.
1 I1 sin(θ)sin(φ)
f (σ ; φ) = K 3 + J 2 K1 cos(θ) − K 2 (B5.148)
cos φ 3 3
Siendo los invariantes ya definidos en apartados previos y los factores K i , para la
función de Mohr-Coulomb clásica, los que a continuación se presentan
K 1 = f 1 (α R ; φ ) K 1 = 1
si α R =1
K 2 = f 2 (α R ; φ ) →K 2 = (B5.149)
K 3 = f 3 (α R ; φ) K = sin (φ)
3
donde α R = R' / RMohr representa el cociente entre la relación de resistencia requerida R ' y
relación de resistencia propia de la función clásica de Mohr-Coulomb RMohr .
Esta nueva función de Mohr-Coulomb Modificada, permite entonces establecer
cualquier relación de resistencias requerida por los distintos materiales, con la sola
modificación de los K i , sin que esto experimente un incremento de la dilatancia (ver
aspecto de la función modificada en la Figura B5.31). A continuación se puede ver las
expresiones que resultan de esta modificación y la forma que adquiere la función de Mohr
al modificar la relación α R = R' / RMohr .
σ II α Reducción
π φ
J1 = 0 Nφ = tan 2 +
4 2
(σ )C
c=
(σ )T 2 Nφ
(σ )T Nφ
2c 2c σI c=α
α Nφ 2
Nφ
(σ )C
R= =α
(σ )T
(σ )C 1≤ α ≤ R
R
α=
Nφ
2c N φ
(σ )T depende de la
σ I ≡ σ II
(σ )C fricción interna
I sin(θ)sin(φ)
F(σ; c; φ) = 1 K 3 + J 2 K 1 cos(θ) − K 2 − c cos(φ) = 0 (B5.150)
3 3
donde
5-64 MODELOS INDEPENDIENTES DEL TIEMPO
(1 + α ) (1 − α )
K1 = − sin(φ)
2 2
(1 + α ) (1 − α ) 1
K2 = − (B5.151)
2 2 sin(φ)
(1 + α ) (1 − α)
K3 = sin (φ) −
2 2
Para mayores detalles en la obtención de estos coeficientes, se recomienda consultar la
referencia de origen1.
Mohr-Coulomb Mohr-Coulomb
Modificado Standard
σ oc
Ro =
σ To π φ
R = tan 2 +
2 π φ 4 2
R = α tan +
4 2
Equivalente a α = 1,0
14 α = 3,61 (60 o
, ≈ 14 )
α = 2,16 σ II α = α1
Rango del hormigón
10 α = α2
8
α 2 > α1
6 (45 , ≈ 6)
o
3
(30 ,3)
o
Dominio del
Hormigón
Mohr-Coulomb
0 α 2 > α1
α2
Mohr-Coulomb α1
modificada
σ II
σI
Línea de puro corte
( θ ≡ 0) ( I 1 ≡ 0)
DINÁMICA NO LINEAL 5-65
σ III
c cot φ
σ II
σI
σ II
J1 = 0
π φ
σ I = σ II Nφ = tan 2 +
4 2
2c
(σ )C − 1 + sin φ
Nφ c=
3 cos φ
2 cos φ
2c
α 2 > α1
3 + sin φ (σ )T 3 + sin φ
c=
2 3 cos φ
cos φ (σ )C 3 + sin φ
2c ≡ 2c Nφ R1 =
− 1 + sin φ (σ )T 3 cos φ
Drucker-Prager Mohr-Coulomb
Standard
( σ ) COM
R=
( σ ) TEN π φ
R = tan 2 +
π φ 4 2
R = α tan 2 +
4 2
Equivalente a α = 1,0
14 α = 3,61 (60 o
, ≈ 14 )
α = 2,16 σ II α = α1
Rango del hormigón
10 α = α2
σI
8
α 2 > α1
6 (45 , ≈ 6)
o
3
(30 ,3)
o
Dominio del
Hormigón
F(I 1 ; J 2 ; c ; φ) = α( φ) I 1 + J 2 − K ( κ; φ) = 0 (B5.152)
donde las funciones de endurecimiento α(φ) y K ( κ; φ) , luego de ser ajustadas con el
criterio de Mohr-Coulomb, resultan
Parámetros
α K
Aplicación
Para compresión triaxial 2 sin(φ) tan (φ)
convencional 3[3 − 6 sin(φ)] [9 + 12 tan (φ)]2 1/ 2
En los últimos años los modelos constitutivos conocidos como de daño continuo han
sido ampliamente aceptados para simular el complejo comportamiento constitutivo de
muchos materiales que se utilizan en ingeniería4,5,6,7. Estos modelo se caracterizan por su
simplicidad en la implementación, versatilidad y consistencia, ya que están basados en la
mecánica de daño continuo.
(∗)
Este apartado ha sido escrito con la colaboración del Dr. Eduardo Car, investigador de la Universidad
Politécnica de Cataluña.
25 Maugin, G. A. (1992). The termodinamics of plasticity and fracture. Cambridge University Press.
26 Kachanov, L. M. (1958). Time of rupture process under creep conditions. Izvestia Akaademii Nauk; Otd Tech
Nauk, 8 26-31.
27 Lemaitre, J and Chaboche, J. L. (1978). Aspects phénoménologiques de la rupture par endommagement. J.
Appl., 2, 317-365.
28 Chaboche, J. (1988). Continuum damage mechanics part I. General Concepts. Journal of Applied Mechanics
55, 59-64.
29 Chaboche, J. (1988). Continuum damage mechanics part II. Damage Growth. Journal of Applied Mechanics 55,
65-72.
30 Simo, J. and Ju, J. (1987). Strain and stress based continuum damage models – I Formulation. Int. J. Solids
analysis of concrete. Second international conference on Computer Aided Analysis and Design of Concrete Structures.
5-68 MODELOS INDEPENDIENTES DEL TIEMPO
σ 0 = M −1 : σ (B5.153)
donde M es el tensor de cuarto orden del modelo de daño anisótropo. Para el caso del
modelo de daño isótropo, la degradación del material se desarrolla en todas las direcciones
por igual y sólo depende de una variable escalar de daño d, con lo que el tensor M se
reduce a M = (1 − d )I y la ecuación de daño anisótropo (B5.153) queda:
σ
σ0 = (B5.154)
(1 − d )
ε ε
σ σ σ0 σ0
CS C0
0 ≤ d ≤1 (B5.155)
33 Rabotnov, I. (1963). On the equation of state for creep. Progress in Applied Mechanics. The Prager Anniversary
Volume. Pp 307-315.
34 Salomón, O.; Oller S.; Car E.; Oñate E. (1999). Thermomechanical fatigue analysis based on continuum mechanics.
donde d=1 representa un estado del material completamente degradado y define la rotura
local completa y d=0 representa un material no dañado. El concepto de tensión efectiva se
formuló por primera vez en conexión con la hipótesis de equivalencia de deformaciones
por Lemaitre-Chaboche3 en 1978 en la siguiente forma, “...la deformación asociada a un
estado dañado bajo una tensión aplicada σ es equivalente a la deformación asociada con el
estado no dañado sometido a una tensión efectiva σ 0 ”. En la Figura B5.34 se observa una
representación esquemática de la hipótesis de tensión efectiva.
Ψ = Ψ (ε; p i ) con p i = {d }
(B5.156)
Ψ = Ψ (ε; d ) = (1 − d ) Ψ0 (ε )
1
Ψ0 (ε ) = ε : C0 : ε (B5.157)
2
donde C 0 es el tensor constitutivo elástico del material en estado no dañado. Para
problemas térmicamente estables es válida la siguiente forma de la desigualdad de Clausius-
Plank,
∂Ψ & ∂Ψ &
Ξ = σ − :ε− d ≥0 (B5.158)
∂ε ∂d
Esta expresión de la potencia disipativa permite hacer las siguientes consideraciones:
a. La inecuación (B5.158) debe cumplirse para cualquier variación temporal de la variable
libre ε , con lo que el multiplicador ε& tiene que ser nulo (método de Coleman, ver
Maugin1). Esta condición proporciona la ley constitutiva hiperelástica para el problema
de daño escalar,
∂Ψ ∂Ψ
σ= , = − Ψ0 ≤ 0 ⇒ − Ψ0 conjugada de d (B5.159)
∂ε ∂d
b. Considerando la ley constitutiva, el valor de la disipación del modelo de degradación
resulta,
Ξ = Ψ0 d& ≥ 0 (B5.160)
Teniendo en cuenta la ecuación (B5.159) se obtiene la siguiente forma de la ecuación
constitutiva
∂Ψ ∂Ψ
σ= = (1 − d ) 0 = (1 − d ) C 0 : ε (B5.161)
∂ε ∂ε
5-70 MODELOS INDEPENDIENTES DEL TIEMPO
σd = + σd
∫t Ξ dt ∫t Ξ dt
d C0
σ σ σ
-1 -1 -1
Cs=(1-d) C0 Cd=[ Cs – C0 ]
C0 C0
ε εu ε εe ε ε εd ε εu
está relacionado con su resistencia a compresión según sea la función umbral de daño que
se elija.
La ecuación (B5.33) representa una superficie límite en el espacio de las deformaciones o
de las tensiones no dañadas. El daño en el material se verifica cuando el valor de f (σ 0 ) es
igual o mayor que c max = σ max por primera vez. Una expresión equivalente a la (B5.162)
está dada por la siguiente expresión,
F (σ 0 ; q ) = G[ f (σ 0 )] − G [c (d )] ≤ 0 , con q ≡ {d } (B5.164)
donde G[χ] es una función escalar, invertible, positiva y de derivada positiva y monótona
creciente.
∂F (σ 0 ; q ) ∂G[ f (σ 0 )]
d& = µ& ≡ µ& (B5.165)
∂[ f (σ 0 )] ∂[ f (σ 0 )]
donde µ es un escalar no negativo denominado parámetro de consistencia de daño,
análogo al factor de consistencia plástico λ , y se utiliza para definir las condiciones de
carga, descarga y recarga a través de las condiciones de Kuhn-Tucker,
∂G[ f (σ 0 )] ∂G[c(d )]
F (σ 0 ; q ) = 0 ⇒ G[ f (σ 0 )] = G[c(d )] ⇒ f (σ 0 ) = c(d ) ⇒ = (B5.167)
∂f (σ 0 ) ∂c(d )
De la condición de permanencia sobre la superficie umbral de daño se deduce que,
∂G[ f (σ 0 )] &
G& [ f (σ 0 )] = f (σ 0 )
∂f (σ 0 )
⇒ d& ≡ G& [ f (σ 0 )] ⇒ µ& ≡ f& (σ 0 ) (B5.169)
d& = µ&
∂G [ f (σ 0 )]
∂[ f (σ 0 )]
∂G[ f (σ 0 )] &
σ& = (1 − d ) C 0 : ε& − f (σ 0 )⋅ [C 0 : ε ] (B5.175)
∂[ f (σ 0 )]
Teniendo en cuenta que la variación temporal de la función umbral se escribe como
∂f (σ 0 ) ∂f (C 0 : ε )
f& (σ 0 ) = : σ& 0 = : ε& (B5.176)
∂σ 0 ∂ε
Reemplazando en la ecuación (B5.175) se tiene:
∂G[ f (σ 0 )] ∂f (C 0 : ε )
σ& = (1 − d ) C 0 : ε& − : ε& ⋅ [C 0 : ε ] (B5.177)
∂[ f (σ 0 )] ∂ε
DINÁMICA NO LINEAL 5-73
r (σ ) 1 − r (σ )
κ& = K (σ 0 ) ⋅ Ξ m = 0
+ 0
⋅ Ξm (B5.179)
g f gC
que define los estado de comportamiento de un punto en función del estado tensional,
siendo x = 0,5 [x + x ] la función de McAully. Las magnitudes g f y g c representan la
máxima disipación de un punto sometido a tracción y a compresión respectivamente. De
esta manera, la disipación de daño siempre estará normalizada respecto de la máxima
energía correspondiente al proceso mecánico que esté realizando en cada momento.
Utilizando la variable κ como variable auxiliar, puede ahora definirse G[χ] en la siguiente
forma general,
c ( κ)
G[c(κ )] = 1 − (B5.180)
f (σ 0 )
b. Si:
i
k −1
[F(σ 0;q )]t + ∆t = k −1i [G[ f (σ 0 )] − G[c(κ )]]t + ∆t
>0 ,
entonces inicia la integración de la ec. constitutiva
t + ∆t ∂G[ f (σ 0 )]
t + ∆t
∂F (σ 0 ; q )
i
i
k [κ]t + ∆t = k −1i [κ]t + ∆t + ki [κ]t + ∆t
i
[σ ]t + ∆t = (1− ki [d ]t + ∆t ) [σ 0 ]t + ∆t
∂G[ f (σ 0 )] ∂f (C 0 : ε )
i t + ∆t
i
[C ] = (1 − d ) C 0 - ∂[ f (σ )] [C 0 : ε] ⊗ ∂ε
T t + ∆t
0
5. Hace k = k + 1 y regresa al punto 2
G [c (d )]= 1 −
c max A 1− c
e
( c(d )
max ) con 0 ≤ c max ≤ c( d ) (B5.181)
c (d )
DINÁMICA NO LINEAL 5-75
G[ f (σ 0 )] = 1 − e (σ 0 ) = c max (B5.182)
0
0
0
con f
f (σ 0 )
siendo A un parámetro que depende de la energía de fractura del material8. El valor de
f 0 (σ 0 ) = c max se obtiene a partir del cumplimiento del criterio de daño para el primer
umbral de degradación, situación en que se cumple G f 0 (σ 0 ) − G c max [ ] [ ] = 0 y también que
[
G f 0
(σ 0 )] = G[c max ] ≡ 0.
En la Tabla 5.3 se presenta la algoritmia de este modelo ya integrado, siendo por lo tanto
más simple su utilización, aunque menos general pues sólo utiliza un ablandamiento
exponencial.
τ
4. Actualiza la tensión y el tensor constitutivo tangente
i
[σ ]t + ∆t = (1− i [d ]t + ∆t ) [σ 0 ]t + ∆t
∂G[ f (σ 0 )]
[C 0 : ε ] ⊗ ∂f (C 0 : ε )
i t + ∆t
i
[C ]
T t + ∆t
= (1 − d ) C 0 -
∂[ f (σ 0 )] ∂ε
c max
1−
c (d ) (B5.183)
G [c (d )]= con 0 ≤ c max ≤ c( d )
1+ A
o también puede expresarse como,
f 0 (σ 0 )
1−
f (σ 0 ) (B5.184)
G [ f (σ 0 )] = con f 0
(σ 0 ) = c max
1+ A
siendo A un parámetro que depende de la energía de fractura del material. El valor inicial
de f 0 (σ 0 ) se obtiene del criterio de daño expresado en la ecuación (B5.33) ó (B5.164) para
el primer umbral de degradación, situación en que se cumple G f 0 (σ 0 ) − G c max [ ] [ ] =0 y
[ ] [ ]
también que G f 0 (σ 0 ) = G c max ≡0.
τ = f (σ 0 ) = 2 Ψ0 (ε ) = ε : C 0 : ε (B5.185)
En este caso el tensor constitutivo tangente resulta teniendo en cuenta la ecuación
(B5.178), esto es:
∂G[τ] 1
C T = (1 − d ) C 0 − [C 0 : ε ] ⊗ [C 0 : ε ] (B5.186)
∂τ τ
σ tmax
τ = f (σ 0 ) = 2 Ψ0 ( ε) = ε C 0 ε = ⇒ σ tmax = τ C 0 (B5.187)
C0
1 1
Ψ0 = ε C0 ε = σ tmax ε =
1 σ tmax ( ) = (τ
2
C0 )
2
=
1 2
τ (B5.188)
2 2 2 C0 C0 2
La disipación total se obtiene integrado la expresión de la disipación y teniendo en cuenta
la ecuación (B5.171) y (B5.172),
1 2 ∂G [τ]
∞ ∞ ∞
∫ Ξ dt = ∫ Ψ0 d& dt = ∫ 2
τ
∂τ
dτ (B5.189)
t =0 t =0 τ0
1
∞
= τ 2 G[τ] − ∫
∞
( τ−τ e 0 ( ) ) dτ
A 1−
τ
τ0
2 τ τ 0 0
1
= ( τ2 − τ τ0 e
A 1−
τ
( )) τ
0
∞ τ2
−
∞
1
( ) ( )
2 A 1−
+ τ0 e τ
τ
0
∞
(B5.190)
2 τ0 2 τ0
A τ0
= τ2 − −
( ) ( )
τ2 τ0 2 τ0 2
− = τ 0 2 1 + 1 ( )
2 2 A 2 A
donde τ 0 = f 0 (σ 0 ) es el valor inicial del criterio de daño (ecuación (B5.33) ó (B5.164)) para
el primer umbral de degradación.
La máxima energía disipada por cada punto será g f y de aquí se tiene
(τ ) 12 + 1A = g
0 2
f ⇒ A=
gf
1
1
(B5.191)
−
(τ ) 0 2 2
Para el caso en que G (τ) sea una función lineal (ver ecuación (B5.183)), el extremo
superior de la misma se obtiene teniendo en cuenta el valor máximo de la función (B5.171),
de donde resulta
τ0
1 −
τ (B5.192)
τ0
G[τ → ∞ ] = 1 = ⇒ τ=−
(1 + A) A
La disipación total en este caso resulta entonces,
5-78 MODELOS INDEPENDIENTES DEL TIEMPO
τ0
−
1 2 ∂G [τ]
∞ A
∫ Ξ dt = ∫ 2
τ
∂τ
dτ
t =0 τ0
τ0
τ0 −
−
1 A
τ − τ0
= τ 2 G[τ]
A
2 τ
− ∫ 1+ A
dτ
0
τ0 (B5.193)
τ 0 τ0
− −
1 A 1 τ2 A
= τ 2 G [τ] − − τ 0 τ
2 τ 0 1+ A 2 τ0
=
1 τ ( )
−
1 0 2
(1 + A ) τ
= −
1 τ0 ( ) 0 2
( ) 2
2 A2 2 A2 2 A
De la misma forma que se hizo anteriormente, se admite que la máxima energía a disipar
por un punto del sólido será g f , por lo tanto igualando con la máxima disipación
(ecuación (B5.193)), se tiene
−
1 τ0 ( )
= gf ⇒
2
A=−
1 τ0 ( ) 2
(B5.194)
2 A 2 gf
T
C = (1 − d ) C 0 − e
A 1−
τ
( )
τ
τ0 + A τ 1
0
[C 0 : ε ] ⊗ [C 0 : ε ] , Abland. Exponencial
τ 2 τ
0
(B5.195)
τ 1
C T = (1 − d ) C −
0 [C 0 : ε ] ⊗ [C 0 : ε ] , Abland. Lineal
τ (1 + A) τ
2
Ambas formas del tensor constitutivo tangente son simétricas. La simetría de este tensor
depende de la norma τ = f (σ 0 ) .
τ = f (σ 0 ) = ε : ε (B5.196)
Esta norma conduce a un tensor constitutivo tangente no simétrico, que para el caso de
un comportamiento post-pico (para τ ≥ τ 0 ) exponencial (ecuación (B5.181)) o lineal
(ecuación (B5.183)), resulta
DINÁMICA NO LINEAL 5-79
T
C = (1 − d ) C 0 − e
A 1−( )
τ
τ
τ0 + A τ 1
0
[C 0 : ε ] ⊗ [ε ] , Abland. Exponencial
τ 2 τ
0
(B5.197)
τ 1
C T = (1 − d ) C −
0 [C 0 : ε ] ⊗ [ε ] , Abland. Lineal
τ (1 + A) τ
2
τ = f (σ 0 ) = σ tmax (B5.198)
1 1
Ψ0 = ε C 0 ε = σ tmax ε =
1 σ tmax ( )
2
=
1 (τ )
2
(B5.199)
2 2 2 C0 2 C0
La disipación total se obtiene integrado la expresión de la disipación en el tiempo, esto es
teniendo en cuenta la ecuación (B5.171) y (B5.172),
1 (τ ) ∂G[τ]
∞ ∞ ∞ 2
∫ Ξ dt = ∫ Ψ0 d& dt = ∫ 2 C 0 ∂τ
dτ (B5.200)
t =0 t =0 τ0
36 Luccioni, B. (1993). Formulación de un modelo constitutivo para materiales ortótropos - Tesis Doctoral,
Universidad Nacional de Tucumán. Argentina.
37 Car, E. (2000). Modelo Constitutivo para el Estudio del Comportamiento Mecánico de los Materiales
∞
1 (τ )2 1 (τ )2 (τ)
∞ ∞
2
1 (τ )2
=
∞
G [τ] − ∫
1
∞
( τ−τ e 0 ( ) ) dτ
A 1−
τ
τ0
2 C0 τ τ C0 0 0
(B5.201)
=
1
( τ2 − τ τ0 e
A 1−( ))τ
τ0
∞
−
1 τ2
∞
1
( ) ( )
2 A 1−
+ τ0 e τ
τ
0
∞
2 C0 τ0
C0 2 τ0
A τ0
=
1 2 1
τ − − =( ) ( )
τ2 τ0 2 τ0 2 τ0 2 1 1
− +
( )
C0 C0 2 2 A C0 2 A
donde τ 0 = f 0 (σ 0 ) es el valor inicial del criterio de daño (ecuación (B5.33) o (B5.164)) para
el primer umbral de degradación.
Al igual que en los otros modelos, la máxima energía disipada por cada punto será g f y
de aquí se tiene la expresión del parámetro A que garantiza una disipación controlada,
(τ )
0 2
1 1
+ = gf ⇒ A=
1
C0 2 A C0 g f 1 (B5.202)
−
(τ ) 0 2 2
Para el caso en que G (τ) sea una función lineal (ver ecuación (B5.183)), el extremo
superior de la misma se obtiene teniendo en cuenta el valor máximo de la función (B5.171),
de donde resulta al igual que en la ecuación (B5.192),
τ0
1 −
τ (B5.203)
τ0
G[τ → ∞ ] = 1 = ⇒ τ=−
(1 + A) A
La disipación total en este caso resulta entonces,
τ0
−
1 (τ ) ∂G [τ]
∞ A 2
∫ Ξ dt = ∫ 2 C0 ∂τ
dτ
t =0 τ0
τ0 τ0
− −
1 (τ )
2 A A
1 τ − τ0
= G [τ] − ∫ dτ
2 C0 τ C0 1 + A
0 τ0 (B5.204)
τ0 τ0
− −
1 (τ )2 A 1 τ2 A
= G [τ] − − τ 0 τ
2 C0 τ 0
C 0 (1 + A) 2 τ 0
=
1 τ ( )
1 0 2
τ
− (1 + A) 2 =−
1 τ ( ) 0 2
( )
0 2
2
2 A C0 2 A C0 2 A C0
De la misma forma que se hizo anteriormente, se admite que la máxima energía a disipar
por un punto del sólido será g f , por lo tanto igualando con la máxima disipación
(ecuación (B5.193)), se tiene
DINÁMICA NO LINEAL 5-81
−
( )
1 τ0
2
= gf ⇒ A=−
( )
1 τ0
2
(B5.205)
2 A C0 2 g f C0
∂τ ∂τ ∂σ 0 ∂σ 0
= : = C 0 : (B5.206)
∂ε ∂σ 0 ∂ε ∂ε
De donde se obtiene el tensor constitutivo tangente considerando un comportamiento
post-pico (para τ ≥ τ 0 ) exponencial (ecuación (B5.181)) o lineal (ecuación (B5.183)). Esto
es,
T
C = (1 − d ) C 0 - e
( )
A 1−
τ
τ
0 τ0 + A τ
[ : ε ] ⊗
:
∂σ 0
, Abland. Exponencial
C 0 C 0
τ2 ∂ε
T τ0 ∂σ 0 (B5.207)
C = (1 − d ) C 0 - [C 0 : ε ] ⊗ C 0 : , Abland. Lineal
τ 2 (1 + A) ∂ε
B6.1 Introducción.
Se ha visto en los capítulos 3 y 4 que la no-linealidad en dinámica está motivada por
cambios de orientación en las fuerzas exteriores debido a grandes movimientos y por no
linealidades en las fuerzas interiores, causada por fenómenos independientes del tiempo,
estudiados en el capítulo 5, y también por fenómenos sensibles al tiempo, que será
presentado en éste capítulo.
Uno de los comportamientos que provoca no-linealidad en la respuesta en el tiempo de
los materiales se debe a la viscoelasticidad, que estudia el comportamiento reológico de los
materiales, es decir, aquellos comportamientos afectados por el transcurso del tiempo.
Existe un gran un trabajo extenso sobre este tema, que puede consultarse en libros
especialmente dedicados a estudiar la influencia del tiempo en los materiales1,2.
En una primera parte de este capítulo, la presentación se limitará a estados que pueden
ser descrito por una sola componente de tensión y deformación, esta forma simplificada
permitirá introducir el concepto y luego se extenderá la formulación al comportamiento
multiaxial. En estos modelos simplificados se hace un símil, tal que la fuerza representa la
tensión y el desplazamiento representa la deformación, y de esta forma se puede explicar en
forma simple el concepto fenomenológico que describe su ecuación de comportamiento.
1 G. Creus (1986). Viscoelasticity – Basic theory and applications to concrete structures. Ed. By C. Brebbia
and S. Orszag. Springer-Verlag. Berlin.
2 R. M. Christensen (1982). Theory of viscoelasticity. An introduction. Academic Press, Inc. N. York.
6-2 MODELOS DEPENDIENTES DEL TIEMPO
ε(t ) = 0 ∀ τ < τ0
t
1 −(t − s ) / r (B6.4)
ε ( t ) = ∫ ξe ⋅ σ( s ) ds ∀ τ ≥ τ0
−∞
DINÁMICA NO LINEAL 6-3
σ e (t ) = C ε e (t )
σ C σ
ε ξ
ξ =r
σ C
ε0
vis vis
σ (t ) = ξ ε& (t )
σ0
τ0 t τ0 Tiempo de retardo t
σ(t ) = 0 ∀ t < τ0
t
− (t − s ) / τ (B6.7)
σ (t ) = ∫C e ε( s ) ds ∀ t ≥ τ0
−∞
Siendo esta, al igual que la ecuación (B6.4), una convolución en el tiempo, cuya solución
exige un alto coste computacional.
En este modelo puede observarse que su respuesta se relaja en el tiempo si se mantiene
constante la imposición de la deformación (ver Figura B6.2). En este sentido aparece un
tiempo denominado de relajación r , cuya expresión es la misma que el tiempo de retardo
introducido en el modelo de Kelvin. Aunque no es evidente en la formulación, como ya se
ha dicho en la introducción, el modelo de Maxwell establece la forma mecánica inversa del
modelo de Kelvin, pudiéndose probarse esto numéricamente.
6-4 MODELOS DEPENDIENTES DEL TIEMPO
σ e (t ) = C ε e (t ) = σ vis (t ) = ξ ε& (t )
σ σ
ξ
ε C σ Tiempo de relajación
σ0
εe ε vis = 1
∫
ξ t
σ( s) ds
ε0
τ0 τ0
t ξ t
=r
C
C 1 ε i (t )
C 0 ε e (t )
σ C1 σ
ξ
C0
ξ ε& i (t )
εe εi
σ(t ) − C1 εi (t ) C0 εe (t ) − C1 εi (t )
σ(t ) = C0 ε e (t ) = C1 εi (t ) + ξ ε& i (t ) ⇒ ε& i (t ) = = (B6.9)
ξ ξ
DINÁMICA NO LINEAL 6-5
Donde la deformación inelástica εi puede considerarse como una variable interna del
modelo. Esta deformación resulta de la solución de la ecuación diferencial (B6.9), para
una imposición tensional σ(t ) a partir de un tiempo t ≥ τ 0 . Esto es,
ε i (t ) = 0 ∀ τ <τ0
i t
1 − (t − s ) / r1 (B6.10)
ε (t ) = ∫ e ⋅ σ( s ) ds ∀ τ ≥ τ0
−∞ ξ
Que comparada con la ecuación (B6.4), puede verse que el modelo clásico de Kelvin
representa ahora, convencionalmente, la parte inelástica de la deformación de este modelo
más general. Al igual que en aquella expresión, r1 = ξ / C1 representa el tiempo de retardo.
La deformación total resulta entonces,
ε(t ) = 0 ∀ τ < τ0
σ(t ) t 1 −(t − s ) / r1 (B6.11)
C 0 −∫∞ ξ
e i
ε ( t ) = ε (t ) + ε ( t ) = + e ⋅ σ( s ) ds ∀ τ ≥ τ0
Definiendo ahora la función uniaxial de fluencia J (t ) como,
1 1 −t / r
J (t ) = + 1− e 1 (B6.12)
C 0 C1
Y realizando una integración por partes en la ecuación (B6.11), suponiendo además que
σ(−∞) = 0 , la deformación total puede escribirse como,
ε(t ) = 0 ∀ τ < τ0
t
dσ( s ) (B6.13)
ε ( t ) = ∫ J (t − s) ⋅ ds ds ∀ τ ≥ τ0
−∞
En el caso particular que la tensión impuesta σ(t ) cumpla con la siguiente simpli-ficación
(ver Figura B6.4),
σ(t ) = 0 ∀ t < τ0 = 0
(B6.14)
σ(t ) = σ 0 ∀ t ≥ τ0 = 0
ε(t ) = 0 ∀ τ < τ0
ε(t ) = ε e (t ) + ε i (t ) = σ 0 + 1 (1-e −t / r1 ) ⋅ σ ∀ τ ≥ τ0
0 (B6.15)
C 0 C1
1442443
ε i (t )
6-6 MODELOS DEPENDIENTES DEL TIEMPO
ε ξ
= r1
σ C1
σ0 σ0
ε= +
C 0 C1
σ0
σ0
ε0 =
C0
τ0 = 0 t τ 0 = 0 Tiempo de retardo t
Figura B6.4 – Respuesta del modelo de Kelvin Generalizado para una tensión impuesta
constante.
ε i (t ) = 0 ∀ τ < τ0
i n − Kelvin (B6.16)
ε (t ) = ∑ ε iα (t ) ∀ τ ≥ τ0
α =1
Por analogía con las expresiones obtenidas previamente, en particular con la (B6.9), se
tiene la variación temporal de la deformación inelástica para un sub-modelo α cualquiera.
Esto es,
σ(t ) − C α ε iα (t )
σ(t ) = C 0 ε e (t ) = C α ε iα (t ) + ξ α ε& iα (t ) ⇒ ε& iα (t ) = (B6.17)
ξα
C1 ε1i (t ) C 2 ε i2 (t ) C n ε in (t )
C 0 ε e (t )
σ C1 C2 Cn σ
ξ1 ξ2 ξn
C0
ξ1 ε& 1i (t ) ξ 2 ε& i2 (t ) ξ n ε& in (t )
εe ε1i ε i2 ε in
ε iα (t ) = 0 ∀ τ < τ0
t
i 1 −(t − s ) / rα (B6.18)
ε
α ( t ) = ∫ ξα e ⋅ σ( s ) ds ∀ τ≥τ
−∞ 0
ε(t ) = 0 ∀ τ < τ0
n − Kelvin
σ(t ) t n − Kelvin 1 −(t − s ) / rα (B6.19)
ε ( t ) = ε e
(t ) + ∑ α ε i
( t ) = + ∫ ∑
ξ
e ⋅ σ( s ) ds ∀ τ ≥ τ0
α =1 C 0 −∞ α =1 α
Definiendo ahora la función uniaxial de fluencia J (t ) para este problema general, como
1 n − Kelvin 1 −t / r
J (t ) = + ∑ 1− e α (B6.20)
C0 α =1 C α
y realizando una integración por partes en la ecuación (B6.19), suponiendo que σ(−∞) = 0 ,
la deformación total puede escribirse como,
ε(t ) = 0 ∀ τ < τ0
t
dσ ( s ) (B6.21)
ε ( t ) = ∫ J (t − s) ⋅ ds ds ∀ τ ≥ τ0
−∞
Puede verse que esta ecuación coincide con la forma general de la ecuación (B6.13), pero
que en este caso la función de fluencia J (t ) es distinta a aquella. En el caso particular que la
tensión impuesta σ(t ) cumpla con la siguiente simplificación,
σ(t ) = 0 ∀ t < τ0 = 0
(B6.22)
σ(t ) = σ 0 ∀ t ≥ τ0 = 0
ε(t ) = 0 ∀ τ < τ0
(B6.23)
ε(t ) = J (t ) σ 0 ∀ τ ≥ τ0
C ∞ ε(t )
σ C∞ σ
C1 ξ1
C1 ε e (t ) ξ1 ε& i (t )
e i
ε =ε−ε ε i
σ ∞ (t ) = C ∞ ε(t )
i (B6.24)
(
σ (t ) = C1 ε(t ) − ε i (t ) = ξ1 ε& i (t ) )
La condición de equilibrio requiere el cumplimiento de la siguiente relación,
( )
σ(t ) = σ i (t ) + σ ∞ (t ) = C1 ε(t ) − ε i (t ) + C ∞ ε(t ) = ξ1 ε& i (t ) + C ∞ ε(t ) (B6.25)
Denominando C 0 = C ∞ + C1 y operando algebraicamente en la última ecuación, resulta
la ecuación de la tensión, cuya expresión es
ε(t ) ε i (t )
C1 ε(t ) = C1 ε i (t ) + ξ1 ε& i (t ) ⇒ = + ε& i (t ) (B6.27)
r1 r1
ε i (t ) = 0 ∀ τ < τ0
i t
1 −(t − s ) / r1 (B6.28)
ε (t ) = ∫ r1 e ⋅ ε( s ) ds ∀ τ ≥ τ0
−∞
σ(t ) = 0 ∀ τ < τ0
t
C1 − (t − s ) / r (B6.29)
σ(t ) = C 0 ε(t ) − r ∫e 1
⋅ ε( s ) ds ∀ τ ≥ τ0
1 −∞
−t / r
G (t ) = [J (t )] = C ∞ + C1 e
−1 1 (B6.30)
Considerando la inversión de la función de relajación para este modelo particular, se
obtiene la siguiente función de fluencia uniaxial,
- ∞ t
C
1 C1 r1C 0
J (t ) = 1 − e (B6.31)
C∞ C0
Y realizando una integración por partes en la ecuación (B6.29), la tensión puede
escribirse en la siguiente forma compacta,
σ(t ) = 0 ∀ τ < τ0
t
dε( s ) (B6.32)
σ (t ) = ∫ G (t − s) ⋅ ds ds ∀ τ ≥ τ0
−∞
t t
dε ( s ) dσ ( s )
ds ⇒ ε(t ) = [σ(t )] = ∫ J (t − s ) ⋅
−1
σ(t ) = ∫ G (t − s ) ⋅ ds (B6.33)
−∞
ds −∞
ds
ε(t ) = 0 ∀ t < τ0 = 0
(B6.34)
ε(t ) = ε 0 ∀ t ≥ τ0 = 0
σ
ε
σ(t = 0) = (C ∞ + C 1 ) ε 0
ε0
C∞ε0
ξ1
= r1
C1
τ0 = 0 t τ 0 = 0 Tiempo de relajación t
Figura B6.7 – Respuesta del modelo de Maxwell Generalizado para una deformación
impuesta constante.
la expresión de la tensión (B6.11) se reduce a la siguiente forma simple,
σ(t ) = 0 ∀ τ < τ0
−t / r (B6.35)
σ(t ) = C ∞ + C1e ε0 ∀ τ ≥ τ0
1
σ(t ) = 0 ∀ τ < τ0
[ ]
n − Maxw
σ(t ) = C ∞ ε(t ) + ∑ C α ε(t ) − ε α (t ) =
i
(B6.36)
α =1
n − Maxw
= C 0 ε(t ) − ∑ C α ε iα (t ) ∀ τ ≥ τ0
α =1
Siendo C 0 = C ∞ + ∑α =1
n − Maxw
C α . Por otro lado la deformación ε es única para todas las
cadenas de Maxwell y la deformación inelástica en la cadena α esima es ε iα , cuya expresión
resulta de la resolución de la siguiente ecuación diferencial (ver ecuación (B6.27)) que
deriva del equilibrio en cada cadena de Maxwell,
ε(t ) ε iα (t )
C α ε(t ) = C α ε (t ) + ξ α ε& iα (t ) ⇒
i
= + ε& iα (t ) (B6.37)
rα rα
ε iα (t ) = 0 ∀ τ < τ0
i t
1 −(t − s ) / rα (B6.38)
ε α (t ) = ∫ e ⋅ ε( s ) ds ∀ τ ≥ τ0
r
−∞ α
σ(t ) = 0 ∀ τ < τ0
n − Maxw
σ(t ) = C 0 ε(t ) − ∑ C α ε α (t ) =
i
(B6.39)
α =1
n − Maxw
C α t −(t − s ) / rα
= C 0 ε(t ) − ∑
r1 −∫∞
e ⋅ ε( s ) ds ∀ τ ≥ τ0
α =1
n − Maxw
−t / r
G (t ) = [J (t )] = C ∞ + ∑
−1 α
Cαe (B6.40)
α =1
t
dε ( s )
σ(t ) = ∫ G (t − s ) ⋅ ds (B6.41)
−∞
ds
C ∞ ε(t )
C∞
C1 ξ1
C1 ε e (t ) ξ1 ε& i (t )
C2 ξ2
σ σ
Cn ξn
εe = ε − εi εi
ε(t ) = 0 ∀ t < τ0 = 0
(B6.42)
ε(t ) = ε 0 ∀ t ≥ τ0 = 0
σ(t ) = 0 ∀ τ < τ0
(B6.43)
σ(t ) = G (t ) ε 0 ∀ τ ≥ τ0
∂Ψ (ε, ε iα , t )
[ ]
n − Maxw n − Maxw
σ(t ) = = C ∞ ε(t ) + ∑ C α ε(t ) − ε iα (t ) ≡ ∑ ξ α ε& iα (t ) (B6.45)
∂ε α =1 α =1
[ ]
n − Maxw n − Maxw
∑ σ α (t ) ⋅ ε iα (t ) ≡ ∑
2
Ξ m (t ) = ξ α ε& iα (t ) ≥ 0 (B6.46)
α =1 α =1
Esta misma expresión puede obtenerse por otro camino, también a partir de la definición
de la propia disipación (Capítulo 2, ecuación B2.51, para θ& = cte ). Esto es,
n − Maxw
& (t ) = σ(t ) ε& (t ) − ∂Ψ ε& (t ) −
Ξ m = σ(t ) ε& (t ) − Ψ ∑
∂Ψ i
ε& α (t )
∂ε α =1 ∂ε iα
n − Maxw
∂Ψ ∂Ψ i
= σ(t ) − ε& (t ) − ∑ ε& α (t ) ≥ 0
14243 ∂ε α =1 ∂ε α
i
0 (B6.47)
n − Maxw n − Maxw
∂Ψ i
Ξm = − ∑ i
ε
& α ( t ) = − ∑ -C α (ε(t ) − ε iα (t )) ε& iα (t ) =
α =1 ∂ε α α =1
[ ]
n − Maxw n − Maxw
∑ ∑
2
Ξm = σ α (t ) ε& iα (t ) = ξ α ε& iα (t ) ≥ 0
α =1 α =1
( ) (
J ijkl (t ) = J V (t ) δ ij δ kl + J D (t ) δ ik δ jl + δ il δ jk )
(B6.48)
V
( ) D
(
Gijkl (t ) = G (t ) δ ij δ kl + G (t ) δ ik δ jl + δ il δ jk )
4 Malvern, L. (1969). Introduction to the mechanics of continuous medium. Prentice Hall, Englewood Cliffs, NJ.
5 Lubliner, J. (1990). Plasticity theory. MacMillan, New York.
DINÁMICA NO LINEAL 6-13
t
dε ( s )
σ
ij
( t ) = ∫ Gijkl (t − s ) ⋅ kl
ds
ds : Modelo viscoelástico de relajación
−∞
t
(B6.49)
dσ kl ( s )
ε
ij (t ) = ∫ J ijkl (t − s) ⋅ ds ds : Modelo viscoelástico de fluencia
−∞
Como puede verse, esta definición es muy general, pero es poco práctica, puesto que
determinar las funciones de fluencia y relajación para un material afectado por fenómenos
viscoelásticos no es una tarea fácil. No obstante esto que se acaba de decir, existe un
camino alternativo basado en la experimentación para la definición de las funciones de
fluencia y relajación. A continuación y a modo de ejemplo se introduce la forma que
establece el “Comité Europeo del Hormigón” (CEB-1978) y la “Federación Internacional
del Pretensado” (FIP-1978), para función de relajación,
[ ]h
J ijkl (t ) = 0.4 β d (t − t c ) C ijkl
−1
(B6.50)
h
Siendo t c el tiempo inicial de referencia, C ijkl el tensor constitutivo del hormigón y
β d una función que tiene en cuenta los efectos reológicos del material. La forma de esta
última es como se muestra en la Figura B6.9.
βd
1
log(tc )
log(t )
Figura B6.9 – Función de respuesta en el tiempo para el hormigón, según (CEB-FIP).
Desafortunadamente esta no es una función exponencial que haga fácil la resolución de
la integral (B6.49), siendo inevitable resolverlas a través de una integral de convolución.
De la ecuación (B6.50) se deduce que siempre es posible aproximar el problema
multiaxial mediante una función tensorial de fluencia, o de relajación, que dependa de una
función escalar. En este sentido se puede practicar otra aproximación basada en la
suposición de un material isótropo, cuyo coeficiente de Poisson ν → 0 . En este caso es
posibles definir una función de deformación de fluencia única para todas las direcciones.
6-14 MODELOS DEPENDIENTES DEL TIEMPO
{ } ∫ J (t − s) ⋅ dσds(s) ds
t
-1
ε ij (t ) ≅ C ijkl ⋅C kl
(B6.51)
−∞
f ( a + b ) = f ( a ) ⋅ f (b ) (B6.53)
se puede evitar la integral de convolución, haciendo una integral en el tiempo de la
siguiente forma,
t t
I (t ) = ∫ f (t − s ) ⋅ g ( s ) ds = I (t − ∆t ) + ∫ f (t − s) ⋅ g (s)ds (B6.54)
−∞ t − ∆t
En cada paso se utiliza la integral del paso de tiempo anterior I (t − ∆t ) y sólo se integra
durante el intervalo de tiempo actual ∆t .
Las funciones exponenciales que aparecen en los modelos viscoelásticos, basados en
analogías muelle-amortiguador, dan como resultado funciones de relajación y fluencia
exponenciales que cumplen con la propiedad del semigrupo. No es este el caso de los
materiales reales, como el hormigón cuya función es como aquella expresada en la ecuación
(B6.50), tal que si se quiere utilizar esta propiedad para evitar la convolución, es necesario
admitir una simplificación.
deformación,
t + ∆t
−1 C − (t + ∆t − s ) / r
∫
−1
ε ij (t + ∆t ) = C ijkl σ kl (t + ∆t ) + C ijkl e 1
⋅ σ kl ( s ) ds
1442443 ξ −∞ (B6.55)
εije (t + ∆t ) 144444 42444444 3
εiij (t + ∆t )
6 J. C. Simo, and T.J.R Hughes, (1998). Elastoplasticity and Viscoplasticity. Computational Aspects. Springer Verlag.
DINÁMICA NO LINEAL 6-15
Tal que la integral de convolución que define la deformación inelástica, puede ahora
resolverse como,
t + ∆t
− (∆t ) / r C − (t + ∆t − s ) / r
∫
−1
ε iij (t + ∆t ) = ε iij (t ) e 1
+ C ijkl e 1
⋅ σ kl ( s ) ds (B6.56)
ξ t
1. Inicialización
0
[ε ]
i t + ∆t
ij [ ]
≡ ε iij
t
0
[σ ]
i t + ∆t
ij ≡ [σ ] i t
ij
C t + ∆t
−(t + ∆t − s ) / r
σ ij (t + ∆t ) = C ijkl ε ij (t + ∆t ) − 1
C0ξ ∫ e 1
⋅ ε kl ( s ) ds
(B6.58)
−∞
En la cual las variables que participan son las ya definidas para el problema uniaxial,
apartado B6.2.5. La solución de esta integral puede llevarse a cabo sin necesidad de hacer
una convolución, para ello se la rescribe como,
6-16 MODELOS DEPENDIENTES DEL TIEMPO
C t
− (t − s ) / r
σ ij (t + ∆t ) = C ijkl ε ij (t + ∆t ) − C ijkl 1
C0ξ ∫e 1
⋅ ε kl ( s ) ds −
−∞
t + ∆t
(B6.59)
C1 − (t + ∆t − s ) / r
− C ijkl
C0ξ ∫ e 1
⋅ ε kl ( s ) ds
t
Luego de integrar por la regla del trapecio el tercer término y reordenar las expresiones,
se escribe la tensión como,
− (∆t ) / r − ( ∆t ) / r
σ ij (t + ∆t ) = C ijkl ε ij (t + ∆t ) − C ijkl ⋅ ε kl (t ) ⋅ e 1
− σ ij (t ) ⋅ e 1
−
C1 −(∆t ) / r1 ∆t
− C ijkl e ⋅ ε kl (t ) − ε kl (t + ∆t ) ⋅
C0ξ
2
(B6.60)
− (∆t ) / r C 1 ∆t C1 ∆t
σ ij (t + ∆t ) = −C ijkl ε kl (t ) ⋅ e 1
1 + + C ijkl ε kl (t + ∆t ) ⋅ 1 − +
C0ξ 2 C0ξ 2
− (∆t ) / r
+ σ(t ) ⋅ e 1
Un posible algoritmo para resolver este problema es el que se muestra en la Tabla 6.2.
1. Obtención de la deformación
[ε ] ij
t + ∆t
2. Integración de la tensión,
[σ ]ij
t + ∆t
[ ]
= σ ij
t
⋅e
− (∆t ) / r
1
− C ijkl [ε kl ] ⋅ e
t − (∆t ) / r
1
1 +
C 1 ∆t
+ C ijkl [ε kl ]
t + ∆t C ∆t
⋅ 1 − 1
C0ξ 2 C0ξ 2
Si se impone una fuerza que de lugar a una tensión constante, este modelo representa
correctamente el problema de fluencia en el tiempo, y si se impone un desplazamiento que
de lugar a una deformación constante se recupera el comportamiento por relajación.
En lo que sigue se presenta la influencia del modelo multiaxial de Kelvin (B6.2.2), cuya
formulación es sencilla y apropiada para representar el amortiguamiento estructural y a la
vez permite establecer una comparación con el ya mencionado concepto de
amortiguamiento de Rayleigh.
Siguiendo la formulación del modelo presentado en el apartado (B6.2.1), puede
sintetizarse su ecuación constitutiva mediante las siguientes expresiones,
ξ
σ(t ) = σ e (t ) + σ vis (t ) = C ε(t ) + ξ ε& (t ) = C ε(t ) + ε& (t )
C
t (B6.61)
σ(t ) ξ 1 −(t − s ) / r1
e vis
ε(t ) = ε (t ) = ε (t ) = − ε& (t ) = ∫ e ⋅ σ( s ) ds
C C −∞
ξ
ξ
σ ij (t ) = σ ije (t ) + σ ijvis (t ) = C ijkl ε kl (t ) + ε& kl (t ) (B6.62)
C
ξ
Siendo r = el tiempo de retardo ya definido en apartados previos.
C
La resolución de este problema en dinámica, exige que se obtenga el campo de
desplazamientos y velocidades a partir de la solución de la ecuación del movimiento, y para
ello puede utilizarse el método de Newmark (ver capítulo B 3). Es decir, la algoritmia es
muy simple y se resume en el siguiente cuadro,
[ε& ]
ij
t + ∆t
= ∇ [u& ]
s
i j
t + ∆t
Tabla 6.3 – Algoritmo para obtener la tensión en problemas dinámicos a través del
modelo de Kelvin .
1
Ψ= ε ij C ijkl ε kl (B6.63)
2
6-18 MODELOS DEPENDIENTES DEL TIEMPO
0 = ∆f k = M kj U&& j (t + ∆t ) + f kint (U
& , U, t + ∆t ) − f ext (t + ∆t )
k ∈ Ω
&& (t + ∆t ) + f
0 = ∆f = M U int & , U , t + ∆t ) − f
(U ext
(t + ∆t ) ∈ Ω
i t + ∆t i t + ∆t i t + ∆t
t + ∆t
0 =Α ∫ ρ N ki N ij dV U&& j +Α ∫ σ in ∇ i N nk dV
S
−Α ∫ t i N ik dS + ∫ ρ bi N ik dV
Ωe Ωe Ωe Ωe
V e e V e e S e Ve e
Ω Ω Ω
y luego operando sobre el término de las fuerzas internas, se tiene la siguiente forma
conocida para la ecuación de equilibrio
DINÁMICA NO LINEAL 6-19
i t + ∆t
ξ S S
∆f k = M kj U&& j (t + ∆t ) + Α ∫
S
C inrs ∇ r N sjU j + C inrs
&
∇ r N sjU j ∇ i N nk dV −
Ωe
V e C e (B6.67)
Ω
− f kext (t + ∆t ) = 0
i t + ∆t
0 = ∆f k = M kj U&& j (t + ∆t ) + Α
Ωe
[ ]
∫ C inrs ∇ r N sj ∇ i N nk dV
S S
U j
Ωe
+
V e
Ωe
i t + ∆t
(B6.68)
ξ S S
+Α ∫ C inrs ∇ r N sj ∇ i N nk dV U& j − f kext (t + ∆t ) = 0
Ω
e
Ωe
V e C
Ωe
t + ∆t
(B6.70)
& (t + ∆t ) = Α B : (r C ) : B dV
f vis (t + ∆t ) = D ⋅ U ∫
Ωe e
&
⋅U [] t + ∆t
Ωe
V e
Ω
En ésta última puede verse que la matriz de amortiguamiento para el modelo de Kelvin es
igual a la de rigidez inicial por el tiempo de retardo.
D=rK (B6.71)
Así, en el caso particular del amortiguamiento de Rayleigh con coeficientes
α = 0 y β = r = ξ/C , conduce al amortiguamiento viscoso de Kelvin,
α =0 , β= r
D = αM +βK → D=rK (B6.72)
ξ && (t + ∆t ) + f int ( U
& , U, t + ∆t ) − f ext (t + ∆t )
Kelvin con r = en : 0 = ∆f = M U
C
≡
Rayleigh con α = 0 y β = r : && (t + ∆t ) + K U(t + ∆t ) + D U
0 = ∆f = M U & (t + ∆t ) − f ext (t + ∆t )
6-20 MODELOS DEPENDIENTES DEL TIEMPO
& , U)
f int (U
644474448
[&& ⋅ U
0= M U
1424
]
3 1444
[
& + f elast (U ) + f vis ( U
424444
&) ⋅U
3 1
]
424
[ ]
& − f ext ⋅ U
3
& (B6.74)
K& Pd P int
{
Pd = f elast
} & = Α B : C : B dV
& = { K ⋅ U} ⋅ U
(U) ⋅ U
e ∫
&
⋅ [U ] e ⋅ U +
Ω V e e Ω
Ω
(B6.75)
{vis &
}
& = D⋅ U
+ f (U) ⋅ U & = Α B : (r C ) : B dV
& ⋅U { e ∫
} []
&
⋅U
Ωe
&
⋅U
Ω V e e
Ω
Resultando naturalmente de esta ecuación el término correspondiente a la potencia de
amortiguamiento.
DINÁMICA NO LINEAL 6-21
& , U)
f int (U
64748
[&& ⋅ U
0= M U
1424
]
3 14
[
& + f elast (U ) ⋅ U
4244
]
3 1 424
[ ]
& − f ext ⋅ U
3
& (B6.76)
K& Pd P int
Pd = f { elast
}
(U ) ⋅ U = {K ⋅ U} ⋅ U = Α ∫ B : C : B dV
& & &
⋅ [U ] e ⋅ U (B6.77)
Ω V e e
e Ω
Ω
Donde no se tiene el término de la potencia de amortiguamiento y por lo tanto no hay
potencia disipativa en el sistema.. Para completar esta formulación debe introducirse un
término disipativo al balance de potencia, pero debe quedar muy claro que es externo al
sistema y no proviene del propio material. Así, la ecuación (B6.76) se escribe,
Pd = f { elast
} & R & &{ & }
(U ) ⋅ U + Pdisip = {K ⋅ U} ⋅ U + D ⋅ U ⋅ U = Α ∫ B : C : B dV
&
⋅ [U ] e ⋅ U +
Ω V e e
e Ω
Ω
(B6.78)
+ [αM + β K ] ⋅ U
& [] Ωe
⋅U
&
La diferencia fundamental entre la utilización de un modelo viscoelástico y la utilización
de un modelo elástico más amortiguamiento estructural (ver ecuación (B6.73)), se presenta
al evaluar la disipación. En el primer caso la disipación proviene del modelo constitutivo.
En el segundo caso de un término añadido. No obstante esto, debe cumplirse que para
toda la estructura, la energía disipada en ambos casos sea la misma,
t
t
E dis (t ) = ∫ Α ∫
C inrs
ξ S
{ }
∇ r N sj ∇ iS N nk dV ⋅ U& j ⋅ Α U& k
Ωe
dt ≡ ( )
∫ D : U& ⋅ U& dt (B6.79)
e
0 Ω V e C Ωe Ω
e
0
En los ejemplos que a continuación se muestran, puede notarse las consideraciones que
se han resaltado sobre los modelos de Kelvin y Rayleigh.
Figura B6.13 – Pórtico. Tensiones vs. Tiempo para: 1) La segunda y tercera respuesta
con amortiguamiento estructural con β = 0.14 y modelo viscoelástico con r = 0.14 ,
respectivamente. 2) La cuarta y quinta respuesta con amortiguamiento estructural con
β = 0.25 y modelo viscoelástico con r = 0.25 , respectivamente.
B6.5 Viscoplasticidad.
La teoría de la viscoplasticidad es un caso general de la teoría de la elasticidad y de la
plasticidad a la vez (ver Figura B6.14).
La viscoplasticidad, a diferencia de la elasticidad y la plasticidad, incorpora el parámetro
de viscosidad ξ como variable del modelo. Esto hace que el modelo sea sensible al
tiempo. Como en los casos anteriores, el estudio que aquí se presenta se centra en la
solución viscoplástica para pequeñas deformaciones.
Las características básicas del modelo viscoplástico en pequeñas deformaciones son:
ε = ε e + ε vp (B6.80)
Φ[F(σ , q)]
λvp = (B6.82)
ξ
Donde Φ[F(σ , q)] es la “función de sobretensión” introducida por Perzyna8 (1963) y
F(σ , q) = [ f (σ ) / K ] − 1 es la función de fluencia viscoplástica, análoga a la función umbral
de plasticidad (ecuación B5.59). Los paréntesis de Mc Aully de definen como
x = 0,5 [x + x ] , y K es la variable interna de endurecimiento que al igual que en
plasticidad puede relacionarse con una resistencia uniaxial equivalente σ = K . La
evolución de esta variable depende de una regla definida a partir de su variación temporal y
puede ser lineal o cuadrática o seguir una forma asignada según la experimentación. La
función de sobretensión se define como,
0 ∀ F ≤ 0
Φ[F(σ , q)] = (B6.83)
Φ[F(σ , q)] ∀ F > 0
dispositivo de fricción
σlim
C ε e (t ) µ
σ σ
ξ
C
εe ε i = ε vp
8Perzyna P. (1963). The constitutive equation for rate sensitive plastic materials. Quart. Appl. Math., 20, 321-
332.
6-26 MODELOS DEPENDIENTES DEL TIEMPO
Φ[F(σ , q)]
λvp = = µ& (B6.84)
ξ
estableciendo que durante el comportamiento viscoplástico ocurra siempre que Φ > 0 , y
permitiendo que se cumpla en cada instante la siguiente relación, cuyo cumplimiento
garantiza un estado de equilibrio viscoplástico:
σ Tiempo de relajación
σ Respuesta viscoplástica
K
Respuesta plástica t→∞
ξ t
=r
C
Φ[F(σ , K )] Φ[F(σ , K ) = 0]
λvp = lim = lim = λ& p ⇒ ε& vp ≡ ε& p (B6.87)
ξ →0 ξ ξ→0 ξ
En esta última expresión puede observarse la necesidad de imponer la
condición de fluencia plástica F(σ , K ) = 0 para que exista el límite buscado.
Siendo éste el único camino para garantizar la existencia del parámetro
viscoplástico, ahora con la forma de factor de consistencia plástica
λvp dt = dλp , en el extremo de viscosidad nula. Desde el punto de vista de la
optimización, puede entenderse la viscoplasticidad como un problema
elastoplástico regularizado. Es decir, se relaja el problema elastoplástico
F(σ , K ) ≥ 0 , con lo que se permite la existencia de una solución fuera del
espacio convexo elástico, pero penalizada por el parámetro de viscosidad ξ .
Estas situaciones límites hacen del modelo viscoplástico una formulación versátil para
abordar problemas comprendidos entre la elasticidad, plasticidad y la propia
viscoplasticidad.
∂G(σ , q) F(σ , K )
ε& vp = λvp ⋅ = ⋅g (B6.90)
∂σ ξ
Por lo tanto la tensión queda definida como,
F(σ , K )
( )
σ& = C:ε& e = C: ε& − ε& vp = C: ε& e −
ξ
⋅ g (B6.91)
9 Maugin, G. A. (1992). The termodinamics of plasticity and fracture. Cambridge University Press.
6-28 MODELOS DEPENDIENTES DEL TIEMPO
(σ )t + ∆t = C: ε t + ∆t − (ε )
i −1 i −1 vp t + ∆t
b. Condición de fluencia viscoplástica
i
[F(σ , K )] t + ∆t f
= i −1
i −1
[(σ)
t + ∆t
] ; donde ε vp = γ (ε vp : ε vp )
−1
K ε
vp
( )
t + ∆t
[F(σ , K )]t + ∆t < 0 , comportamiento elástico, i (σ )t + ∆t ≡ (σ * )
i i t + ∆t
Sí: ⇒ , Ir
al punto “ f. ”
Sí: i
[F(σ , K )]t + ∆t ≥ 0 Entonces se tiene comportamiento viscoplástico
Hacer i
[r ]t + ∆t = i [ F(σ , K ) ]t + ∆t
ξ
i
[ε ] vp t + ∆t
=
i −1
[ε ]
vp t + ∆t
+ λvp
i
[ ] t + ∆t
⋅ [g]
i t + ∆t
∆t = [∆µ ]
i
⋅ vp t + ∆t i
[g] t + ∆t
d. Corrección de la tensión
i
(σ )t + ∆t = i −1 (σ )t + ∆t − C : λvp
i
1444
t + ∆t i
[ ]
⋅ [g] ∆t
t + ∆t
424444 3
t + ∆t i
⋅ [g]
i t + ∆t
vp
∆µ
f (σ ) J2
Φ[F(σ , q )] ≡ F(σ , q) = vp − 1 ≡ vp
−1 (B6.92)
K (ε ) K ( ε )
DINÁMICA NO LINEAL 6-29
1 3 2
J 2 = s : s = τ oct o también J2 =
s
=
(s2
1 + s 22 + s 32 ) 1
2
(B6.93)
2 2 2 2
donde s = σ - 13 tr (σ ) I , o s ij = σ ij − 13 σ kk δ ij , es la tensión desviadora.
Considerando la ecuación (B6.90) y una función de potencial viscoplástica asociada a la
función de fluencia, ecuación (B6.92), se obtiene el siguiente tensor de flujo viscoplástico,
∂G ∂G ∂s
G(σ , q) = J 2 − Kˆ ( ε vp ) ⇒ g= = =s (B6.94)
∂σ ∂s ∂σ
sustituyendo esta última en la ecuación (B6.90), resulta la siguiente deformación
viscoplástica,
J2
vp
−1
∂G(σ , q) F (σ , K ) K ( ε ) (B6.95)
ε& vp = λvp ⋅ = ⋅g = ⋅s
∂σ ξ ξ
Debido a que se propone un potencial viscoplástico del tipo von Mises, se escribe una
ley constitutiva apropiada para representar un comportamiento con desviación dominante.
De esta forma, se tendrá un comportamiento volumétrico elástico y sólo la parte
desviadora sufrirá los efectos de la viscoplasticidad11. Por ésta razón sólo se describirá e
integrará la ley constitutiva viscoplástica desviadora, pero es importante recordar que para
obtener la tensión completa hay que añadir el comportamiento elástico volumétrico σ ijvol .
La tensión desviadora se escribe entonces,
s = 2 G (e − e vp ) ⇒ ∆s = 2 G (∆e − ∆e vp ) (B6.96)
Además, la tensión desviadora en el instante t + ∆t valdrá,
s t + ∆t = s t + ∆s t + ∆t = s t + 2 G ∆e t + ∆t − 2 G (∆e vp ) t + ∆t
1442443 (B6.97)
(s ∗ ) t + ∆t
10 S.Oller (2001). Fractura mecánica – Un enfoque global. CIMNE – Ediciones UPC. Barcelona.
11 La ley constitutiva descompuesta en su parte volumétrica y desviadora se escribe como,
σ ij = λLame ε ekk δ ij + 2 µ Lame (ε ije − 13 ε ekk δ ij ) = λLame ε ekk δ ij + 2 µ Lame eije
siendo las constantes de Lame λLame = E 3(1 − 2 ν ) y µ Lame = G = E 2(1 + ν) , donde E y ν son
respectivamente el módulo de Young y el coeficiente de Poisson. Debido a que el flujo de von Mises es
vp vp vp
[ ]
dominantemente desviador ( ε& kk = 0 ⇒ ε& ij ≡ e&ij = Φ ξ s ), la ley constitutiva anterior queda
particularizada en la siguiente forma,
σ ij = λLame ε kk δ ij + 2 G (eij − eijvp ) = σ ijvol − s ij
6-30 MODELOS DEPENDIENTES DEL TIEMPO
[
Φ F t + ∆t ]
s t + ∆t
= (s ) ∗ t + ∆t
− 2G
ξ
s t + ∆t ∆ t , con [ ]
Φ F t + ∆t ≡ F t + ∆t = F(s ∗ ) t + ∆t (B6.98)
( s ∗ ) t + ∆t
s t + ∆t =
∆t (B6.99)
1+ 2 G F ( s ∗ ) t + ∆t
ξ
En esta última ecuación puede verse las siguientes características de comportamiento
viscoplástico,
1. Sí ∆t → ∞ o ξ → 0 , para que haya solución debe alcanzarse la condición de
ΦF [ t + ∆t
]
= 0 y para ello hay que aplicar la condición de consistencia plástica y por lo
tanto transformarlo en un problema elastoplástico (ver ecuación (B6.84)).
2. Sí ∆t → 0 o ξ → ∞ , ocurre que la tensión buscada coincide con la predictora
t + ∆t ∗ t + ∆t
s = (s ) , lo que implica que el comportamiento durante el intervalo ∆t ha
sido elástico lineal.