Tesis: Caracterización Mecánica de Materiales Constituídos por partículas
Autor: Alejo O. Sfriso
Director: Eduardo Núñez (Univ. Buenos Aires)
Codirector: Guillermo G. Weber (Univ. Texas en Brownsville)
Comité de seguimiento: Eduardo N. Dvorkin
Fecha de presentación: 15/12/2008
Facultad de Ingeniería, Universidad de Buenos Aires
Contenido
1 Introducción 7
1.1 Motivación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2 Objeto de esta tesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.3 Organización y presentación . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.4 Notación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.5 Agradecimientos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2 Comportamiento de las arenas 15
2.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.2 Definiciones iniciales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.2.1 Caracterización del empaquetamiento . . . . . . . . . . . . . . . . 16
2.2.2 La definición de resistencia en los materiales friccionales . . . . . 17
2.3 Resistencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.3.1 Componentes de la resistencia al corte . . . . . . . . . . . . . . . 18
2.3.2 Modelo de los dos bloques . . . . . . . . . . . . . . . . . . . . . . 18
2.3.3 Ensayo triaxial . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.3.4 Estado crítico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.3.5 Estado de deformación constante . . . . . . . . . . . . . . . . . . 23
2.3.6 Resistencia máxima y presión media . . . . . . . . . . . . . . . . 24
2.3.7 Estado triaxial general . . . . . . . . . . . . . . . . . . . . . . . . 27
2.4 Rigidez . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2
2.4.1 Elasticidad y contacto Hertziano . . . . . . . . . . . . . . . . . . 30
2.4.2 Rango elástico y disipación de energía . . . . . . . . . . . . . . . 32
2.4.3 Módulo de Young inicial . . . . . . . . . . . . . . . . . . . . . . . 34
2.4.4 Teoría de tensión — dilatancia de Rowe . . . . . . . . . . . . . . . 36
2.4.5 Cambio de volumen en descarga . . . . . . . . . . . . . . . . . . . 40
2.4.6 Reducción de la rigidez secante por deformación . . . . . . . . . . 41
2.4.7 Módulo secante E50 . . . . . . . . . . . . . . . . . . . . . . . . . . 46
2.4.8 Compresión proporcional . . . . . . . . . . . . . . . . . . . . . . . 46
2.4.9 Anisotropía y endurecimiento cinemático . . . . . . . . . . . . . . 50
3 Modelización de la compresión isotrópica 53
3.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.2 Elementos del modelo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.2.1 Cinemática . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.2.2 Variables de estado . . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.2.3 Elasticidad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.2.4 Función de fluencia y potencial plástico . . . . . . . . . . . . . . . 54
3.2.5 Ecuaciones de evolución . . . . . . . . . . . . . . . . . . . . . . . 55
3.3 Implementación numérica . . . . . . . . . . . . . . . . . . . . . . . . . . 57
3.3.1 Operador de integración . . . . . . . . . . . . . . . . . . . . . . . 57
3.3.2 Algoritmo de integración . . . . . . . . . . . . . . . . . . . . . . . 57
3.4 Resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
3.4.1 Influencia de los parámetros de entrada . . . . . . . . . . . . . . . 60
3.4.2 Comportamiento del algortimo . . . . . . . . . . . . . . . . . . . . 61
3.A Apéndice: Selección de e0 como variable de estado . . . . . . . . . . . . . 65
4 Modelización de la falla por corte 68
4.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
4.2 Estructura matemática del modelo . . . . . . . . . . . . . . . . . . . . . 69
3
4.2.1 Cinemática . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
4.2.2 Variables de estado . . . . . . . . . . . . . . . . . . . . . . . . . . 70
4.2.3 Relación tensión-deformación . . . . . . . . . . . . . . . . . . . . 70
4.2.4 Función de fluencia . . . . . . . . . . . . . . . . . . . . . . . . . . 70
4.2.5 Regla de flujo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
4.2.6 Condición de plasticidad . . . . . . . . . . . . . . . . . . . . . . . 71
4.2.7 Ecuaciones de evolución . . . . . . . . . . . . . . . . . . . . . . . 71
4.3 Funciones de estado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
4.3.1 Elasticidad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
4.3.2 Ángulo de fricción interna . . . . . . . . . . . . . . . . . . . . . . 76
4.3.3 Función de fluencia . . . . . . . . . . . . . . . . . . . . . . . . . . 79
4.3.4 Regla de flujo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
4.4 Implementación numérica . . . . . . . . . . . . . . . . . . . . . . . . . . 88
4.4.1 Operador de integración . . . . . . . . . . . . . . . . . . . . . . . 88
4.4.2 Algoritmo de integración . . . . . . . . . . . . . . . . . . . . . . . 89
4.5 Resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
4.5.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
4.5.2 Simulación de ensayos triaxiales . . . . . . . . . . . . . . . . . . . 100
4.5.3 Simulación de un ensayo de deformación plana . . . . . . . . . . . 107
4.5.4 Simulación del comportamiento de una presa CFRD sujeta a un
sismo destructivo . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
4.A Apéndice: Deducciones de algunas funciones de estado . . . . . . . . . . 115
4.A.1 Hiperelasticidad . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
4.A.2 Función de fluencia . . . . . . . . . . . . . . . . . . . . . . . . . . 119
4.A.3 Regla de flujo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
4.A.4 Determinación de pr y φc para arenas del Puelchense . . . . . . . 123
5 Modelización de la carga monotónica de corte 127
5.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
4
5.2 Estructura matemática del modelo . . . . . . . . . . . . . . . . . . . . . 128
5.2.1 Cinemática . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
5.2.2 Variables de estado . . . . . . . . . . . . . . . . . . . . . . . . . . 128
5.2.3 Relación tensión-deformación . . . . . . . . . . . . . . . . . . . . 128
5.2.4 Función de fluencia . . . . . . . . . . . . . . . . . . . . . . . . . . 128
5.2.5 Regla de flujo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
5.2.6 Condición de plasticidad . . . . . . . . . . . . . . . . . . . . . . . 129
5.2.7 Ecuaciones de evolución . . . . . . . . . . . . . . . . . . . . . . . 129
5.3 Funciones de estado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
5.3.1 Variable de estado y función de fluencia . . . . . . . . . . . . . . 130
5.3.2 Relación tensión-deformación . . . . . . . . . . . . . . . . . . . . 132
5.3.3 Función de endurecimiento para ς . . . . . . . . . . . . . . . . . . 134
5.3.4 Resistencia máxima y evolución de variables de estado . . . . . . 136
5.4 Implementación numérica . . . . . . . . . . . . . . . . . . . . . . . . . . 138
5.4.1 Operador de integración . . . . . . . . . . . . . . . . . . . . . . . 138
5.4.2 Algoritmo de integración . . . . . . . . . . . . . . . . . . . . . . . 139
5.5 Resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
5.5.1 Lista de parámetros . . . . . . . . . . . . . . . . . . . . . . . . . . 145
5.5.2 Análisis de parámetros materiales . . . . . . . . . . . . . . . . . . 145
5.5.3 Análisis de parámetros de estado . . . . . . . . . . . . . . . . . . 149
5.5.4 Calibración para la arena Puelchense . . . . . . . . . . . . . . . . 157
5.5.5 Calibración para la arena Toyoura . . . . . . . . . . . . . . . . . . 157
5.5.6 Análisis de las variables de estado . . . . . . . . . . . . . . . . . . 164
5.5.7 Análisis del tamaño del paso . . . . . . . . . . . . . . . . . . . . . 169
5.5.8 Simulación de un ensayo de compresión plana . . . . . . . . . . . 171
5.5.9 Simulación de un problema de valores de contorno complejo . . . 174
6 Consideraciones finales y oportunidades de investigación futura 179
6.1 Resumen del comportamiento de las arenas . . . . . . . . . . . . . . . . . 179
5
6.2 Resumen de los modelos constitutivos para arenas . . . . . . . . . . . . . 181
6.3 Características del modelo constitutivo desarrollado en esta tesis . . . . . 182
6.4 Caminos de investigación futura . . . . . . . . . . . . . . . . . . . . . . . 184
7 Apéndice: Algebra Tensorial 185
7.A Tensores unitarios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 185
7.B Invariantes de tensores de segundo orden . . . . . . . . . . . . . . . . . . 186
7.C Derivadas de invariantes . . . . . . . . . . . . . . . . . . . . . . . . . . . 187
7.D Medidas de tensión y deformación . . . . . . . . . . . . . . . . . . . . . . 188
7.E Derivadas de las medidas de tensión . . . . . . . . . . . . . . . . . . . . . 189
6
Capítulo 1
Introducción
1.1 Motivación
En los últimos diez años, la geomecánica computacional se ha convertido en una
herramienta para el análisis rutinario de problemas geotécnicos, hasta el punto que
hoy existen numerosos programas computacionales que aplican métodos numéricos
específicamente orientados a geomateriales, como Plaxis [18], FLAC [30] y otros. Estos
programas permiten la simulación numérica del comportamiento de estructuras típicas
de la ingeniería geotécnica como muros, fundaciones, excavaciones, túneles y presas,
construidas en suelos de todo tipo, desde arcillas blandas hasta macizos rocosos.
Sin embargo, la capacidad predictiva de esta técnica de análisis todavía es poco
confiable, aún para problemas rutinarios. Una prueba de ésto se dio en 1994, cuando la
FHWA1 construyó y ensayó cinco zapatas cimentadas sobre arenas uniformes. Antes de
efectuar los ensayos, se le pidió a 31 participantes que predijeran la carga que produciría
un determinado asentamiento. En el informe final, Briaud y Gibbens [17] informaron:
“Nadie entregó un conjunto de respuestas que cayera consistentemente dentro del 20%
de desvío respecto de los valores medidos. Dos participantes (de 31) tuvieron el 80% de sus
respuestas dentro de este rango.... Como 31 participantes predijeron el comportamiento
1Asociación de Carreteras Federales de los Estados Unidos.
7
de 5 fundaciones, existieron 155 valores para el coeficiente de seguridad. Sólo uno de los
coeficientes de seguridad fue menor que uno; el siguiente peor fue 1.6; el promedio fue
de 5.4. Por lo tanto parece que nuestra profesión conoce cómo diseñar bases aisladas
de manera muy segura... Considerando los elevados coeficientes de seguridad y los
valores reducidos de asentamiento, la carga de diseño pudo haber sido considerablemente
superior. Por lo tanto parece que nuestra profesión podría diseñar bases aisladas más
económicamente.”
Aunque la predicción global fue pobre, peor fue el resultado del sub-grupo que empleó
métodos numéricos: el promedio de las predicciones fue del orden del 50% de los valores
medidos2.
La razón básica para este pobre desempeño está en la complejidad intrínseca de
los materiales friccionales: son aquellos cuya resistencia y rigidez dependen del estado
tensional y cambian de manera muy significativa durante los procesos de deformación.
Dentro de los materiales friccionales, las arenas presentan una característica distintiva:
su rigidez y resistencia al corte son nulas cuando la presión es nula. Esta característica,
que permite por ejemplo la fabricación de relojes de arena, complica enormemente la
toma de muestras, la ejecución de ensayos experimentales repetibles y, sobre todo, su
interpretación.
Hay dos características del comportamiento de las arenas que son bien conocidas:
1. Todas las arenas tienen un ángulo de reposo, que es independiente de los procesos
de carga y deformación que pueda haber sufrido el material: cuando termina su
tiempo, el reloj de arena siempre muestra el mismo cono en su mitad inferior.
2. Independientemente del estado inicial, las arenas se deforman a volumen constante
luego de una distorsión suficientemente grande. Si se sacude un reloj de arena, el
cono se asienta y el volumen de arena disminuye; si luego se invierte el reloj (y se
2La experiencia se repitió en 2006 en la Universidad de Western Australia con resultadossorprendentemente similares.
8
espera), el volumen de arena que queda en la mitad inferior es igual al volumen que
había antes de sacudir el reloj.
Todos los modelos constitutivos sofisticados para arenas contemplan este
comportamiento básico y tienen, entre sus parámetros materiales, al menos dos que
permiten calcular la pendiente del ángulo de reposo y la densidad terminal del material.
Sin embargo, en la práctica profesional se prefiere el empleo de modelos constitutivos
simples de elasticidad lineal y plasticidad perfecta - del tipo de Mohr-Coulomb - o modelos
constitutivos de endurecimiento isotrópico y elasticidad dependiente de la presión como
el modelo hiperbólico de Duncan-Chang. Esto ocurre porque estos modelos constitutivos
son bien conocidos por los usuarios y porque existe una amplia base de datos para sus
parámetros de entrada, y ocurre a pesar que ninguno de estos modelos podría predecir
el ángulo de reposo o la densidad terminal.
Además, los modelos constitutivos mencionados tienen el severo inconveniente de
que los parámetros de entrada dependen del problema que se estudia. Por ejemplo, se
requiere un conjunto de parámetros para caracterizar la arena que rodea el fuste de un
pilote y otro conjunto para la arena que rodea la punta, aún cuando todo el pilote esté
embebido en un depósito uniforme. Esto ocurre porque algunos de los parámetros de
entrada dependen de la densidad del material y de su estado tensional, variables ambas
que alcanzan valores muy diferentes en distintos puntos de la malla y que cambian de
manera muy significativa durante el proceso de deformación que se pretende modelizar.
Recientemente, el modelo hiperbólico disponible en Plaxis fue aumentado con elasticidad
para baja deformación, al costo de dos parámetros materiales adicionales a los ocho que
ya tenía [6]. Debe reconocerse que esta contribución sugiere que el modelo hiperbólico
está alcanzando un grado de madurez excesivo, y que su reemplazo podría ser bienvenido
por la industria.
El campo de los modelos constitutivos para arenas presenta oportunidades de
desarrollo en varios aspectos, que se pueden agrupar en:
1. La búsqueda del “modelo universal” que capture el comportamiento de las arenas
9
para todas las escalas de tensión, deformación y velocidad de deformación.
2. La implementación robusta de algunos de los modelos constitutivos avanzados
disponibles.
3. El desarrollo de modelos predictivos, prácticos y simples orientados a la industria.
La mayoría de los modelos constitutivos sofisticados que la academia ha aportado en
estos años fueron diseñados como herramientas para la comprensión del comportamiento
de las arenas más que para resolver problemas industriales. Por eso, tienen algunos
parámetros que no pueden ser medidos en ensayos convencionales y nunca han sido
implementados en códigos generales.
Un modelo constitutivo para arenas que aspire a ser empleado por la industria debe
tener, como mínimo, las siguientes características:
1. Respecto de los parámetros de entrada
(a) Deben ser determinables mediante procedimientos normalizados y conocidos
por la comunidad geotécnica.
(b) Deben ser independientes del problema que se estudia: un único juego de
parámetros debe capturar el comportamiento del material en todo el rango de
tensión, densidad y deformación de interés ingenieril.
2. Respecto de la capacidad predictiva
(a) La resistencia al corte debe ser reproducida con precisión o conservadoramente
para cualquier trayectoria de carga que produzca la falla del material.
(b) La deformación continuada debe conducir a un estado terminal, denominado
estado crítico, en el que todas las variables de estado alcanzan un valor
estacionario.
10
(c) Las curvas “tensión-deformación” deben ajustarse razonablemente, al menos
para los ensayos de laboratorio y trayectorias de carga que ocurren con mayor
frecuencia en las aplicaciones prácticas.
3. Respecto de la implementación numérica
(a) Debe ser suficientemente robusta como para que se puedan resolver problemas
de valores de contorno complejos.
(b) Debe emplear variables disponibles en los códigos de elementos finitos
generales para geotecnia, que son básicamente códigos de elastoplasticidad.
1.2 Objeto de esta tesis
El primer objeto de esta tesis es el desarrollo de un modelo constitutivo para arenas con
las características descriptas en la sección anterior, y orientado al análisis rutinario de
problemas de interacción suelo — estructura. Los problemas que deben ser simulados son,
entre otros: i) cálculo de empujes y deformaciones de muros y tablestacas; ii) cálculo
de resistencia y rigidez de bases y pilotes sometidos a carga estática; iii) problemas
de interacción estática y dinámica entre suelo y estructuras; y iv) asentamiento de
terraplenes y presas de materiales sueltos sometidos a acción sísmica; v) cálculo de
coeficientes de seguridad; y vi) otros problemas en los que la deformación plástica por
corte monotónico domina el compotamiento del material.
El segundo objeto de la tesis es la implementación del modelo en un programa de
elementos finitos general, con una robustez suficiente como para que pueda ser aplicado
a problemas reales de cualquier naturaleza. Se ha elegido efectuar la implementación en
el programa de elementos finitos Plaxis, del cual la Facultad de Ingeniería UBA posee
una licencia académica.
11
1.3 Organización y presentación
El documento está organizado en seis capítulos: i) esta introducción; ii) una descripción
del comportamiento mecánico de los materiales friccionales; iii) la formulación e
implementación de un modelo constitutivo para compresión isotrópica; iv) la formulación
e implementación de un modelo constitutivo para falla por corte con deformaciones
grandes; v) la formulación e implementación de un modelo constitutivo para corte
monotónico y deformaciones moderadas; y vi) algunas consideraciones finales y
oportunidades de investigación futura.
1.4 Notación
Las relaciones funcionales se indican con corchetes f = f [x]. Los tensores de segundo
orden se escriben en negrita minúscula; los de cuarto orden en negrita mayúscula. Los
tensores unitarios de segundo y cuarto orden se escriben respectivamente como 1 y I. La
doble contracción, la contracción simple y la diada se escriben como a : b, a · b and ab.
La compresión es positiva. La presión media se define como
p =1
3σ : 1 > 0 (1.1)
donde σ es el tensor de tensiones efectivas. No se admiten tensiones de tracción. La
deformación volumétrica se define como
εv = ε : 1 (1.2)
donde ε es el tensor de deformaciones infinitesimales y εv > 0 implica reducción de
volumen.
Todas las operaciones son tensoriales excepto que se indique otra cosa en el texto.
Para la presentación de resultados, los autovalores de σ y ε se ordenan σ1 ≥ σ2 ≥ σ3 > 0
12
y ε1 ≥ ε2 ≥ ε3.La proyección deviatórica de cuarto orden
Idev = I− 1311 (1.3)
se emplea para definir el tensor de tensiones y deformaciones deviatóricas
s = Idev : σ (1.4)
εd = Idev : ε (1.5)
También se emplea el tensor de oblicuidad de tensiones
r = s/p (1.6)
y su norma r = ‖r‖ como medidas del estado tensional.
1.5 Agradecimientos
Es imposible aquí ensayar palabras de gratitud para Eduardo Núñez, mi Director de
Tesis3. Núñez es mi maestro. Agregó conmigo un eslabón a una larga cadena que sigue
hacia atrás con A. Bolognesi, primer profesor de Mecánica de Suelos de la Universidad de
Buenos Aires. Núñez define esta circunstancia de una manera impersonal cuando dice:
“usted contrae una deuda con la ingeniería”. No es cierto. Yo contraigo una deuda con
él. Mi promesa es que nuevos eslabones premiarán su esfuerzo, pues así es la relación
entre un maestro y su discípulo, y porque esa es la retribución que él merece y espera.
Otro profesor que tuvo una muy fuerte influencia en mi trabajo fue Eduardo Dvorkin,
quien me transmitió su pasión por la mecánica computacional y, con infinita paciencia,
me brindó su consejo y orientación a lo largo de los años. Fue la combinación de mi
3Aquí abandono la redacción impersonal.
13
experiencia con Núñez y Dvorkin lo que me orientó a desarrollar un trabajo en el campo
de la geomecánica computacional.
Emir Macari, profesor entonces en Georgia Tech y autoridad en el ámbito de los
modelos constitutivos para geomateriales, me invitó a una visita de estudio en 1998 que
me permitió conocer el sistema universitario estadounidense, hacer algo de experiencia
en los laboratorios de Georgia Tech y aprovechar su vasta biblioteca. En una segunda
visita efectuada en 2006, cuando era el decano de ingeniería de la Universidad de Texas
en Brownsville, Emir Macari me presentó a Guillermo Weber, quien fue desde entonces
mi co-director de tesis, guía y fuente de referencia para la implementación numérica del
modelo. El trabajo conjunto con Guillermo Weber fue continuo desde 2006 hasta la fecha
y fue esencial en el logro del objetivo del programa de tesis.
A. Verri, P. Sanz y V. Calo - primero mis estudiantes y pronto mis amigos y
compañeros de estudio - tuvieron largas discusiones conmigo sobre muchos aspectos de
geotecnia, plasticidad y métodos numéricos. Mis colaboradores J. Laiún, M. Rousset, A.
Comas, G. Quaglia, P. Serigós, A. Abad y J. Hasbani me ayudaron a preparar muchos
de los gráficos.
Por último, mi familia. Betty me ayudó a terminar algo que comenzó antes que
ella llegase a mi vida. Mis cinco hijos me prestaron buena parte del tiempo que les
correspondía y que yo utilizé para desarrollar este trabajo. Mi padre, ingeniero geotécnico
como yo, me enseño a ver y a aprender de lo que veía. Mucho más importante, me enseñó
a vivir y a disfrutar de mi trabajo.
Para todas estas personas, el reconocimiento volcado en estas líneas es relativamente
poco importante. Sin embargo, es importante para mí reconocer y agradecer su
invalorable ayuda y apoyo.
14
Capítulo 2
Comportamiento de las arenas
2.1 Introducción
Se puede presentar el problema conceptual del comportamiento mecánico de las arenas
con dos referencias históricas:
Reynolds ([117], citado en [82]) observó que: “Un empaquetamiento denso
de granos (colocado) dentro de una bolsa flexible invariablemente incrementa
su volumen cuando la envoltura es deformada; si la envoltura es inextensible,
pero no inflexible, no es posible ninguna deformación hasta que las fuerzas
aplicadas rompen la bolsa o fracturan los granos.”
Mohr publicó ([96], citado en [151]): “Las deformaciones observadas
en un cuerpo homogéneo luego de que se alcanza el límite de elasticidad
no están confinadas a los más pequeños dominios del cuerpo. Consisten
aproximadamente en el hecho de que partes del cuerpo de dimensiones finitas
se desplazan respecto a otras en dos grupos de bandas de deslizamiento...”
Reynolds mostró por primera vez una característica de los materiales granulares:
la dilatancia, que implica que un cambio en las tensiones de corte produce un cambio
de volumen. Como el empaquetamiento de Reynolds es denso, el volumen crece y la
15
dilatancia es positiva. Los empaquetamientos sueltos, por el contrario, exhiben dilatancia
negativa, lo que quiere decir que su volumen disminuye cuando se les aplica un esfuerzo
de corte. La deformación no uniforme descripta por Mohr se observa claramente en los
materiales granulares densos, fundamentalmente debido a que algunas zonas incrementan
su volumen más que otras, lo que conduce a inestabilidad y localización de deformaciones.
En este capítulo se describe el comportamiento mecánico de las arenas. Las fuentes
bibliográficas principales son [59][77][89][92][128][131][144][147][148][151].
2.2 Definiciones iniciales
2.2.1 Caracterización del empaquetamiento
Los materiales granulares son medios porosos compuestos por partículas sólidas y vacíos.
La densidad del empaquetamiento se caracteriza mediante la relación de vacíos
e = vv/vs (2.1)
donde vv es el volumen de los vacíos y vs es el volumen de las partículas sólidas. También
se utiliza la porosidad n = vv/ (vv + vs) y el volúmen específico v = (vv + vs) /vs. La
relación de vacíos tiene un valor máximo normalizado emax y uno mínimo emin que
permiten definir la densidad relativa
Dr =emax − eemax − emin
(2.2)
De acuerdo con esta definición, una arena es suelta si Dr < 1/3, medianamente densa
si 1/3 < Dr < 2/3, densa si 2/3 < Dr < 3/4, y muy densa si Dr > 3/4.
16
2.2.2 La definición de resistencia en los materiales friccionales
En ingeniería se denomina “resistencia” a la capacidad que tiene un sólido de resistir
acciones que tienden a cambiar su forma. En el caso de los materiales friccionales, ésta
“resistencia” se identifica con la resistencia al corte. La resistencia al corte de un material
friccional incoherente no es una propiedad del material, sino que depende de: i) la presión
integranular; ii) la densidad; y iii) la dirección del plano de eventual cizallamiento con
respecto a las direcciones principales de anisotropía material.
La tensión de corte τ f que produce el cizallamiento de un plano se calcula con la
ecuación de Mohr — Coulomb
τ f = σn tan [φ] (2.3)
donde σn es la tensión normal que actúa en el plano y φ es el ángulo de fricción interna.
Esta expresión se aplica en todos los cálculos prácticos de la ingeniería geotécnica; en
cada caso, se elige φ en función del problema en estudio1.
[ ]tanf nτ σ φ=
fτ
τφ
3σ 1σ
3σ2σ
1σ3σ nσ σ
1σ
Figura 2-1: La ec. (2.3) del criterio de Mohr - Coulomb en coordenadas τ − σ.
1Esta es una de las prácticas que requieren mayor experiencia profesional.
17
2.3 Resistencia
2.3.1 Componentes de la resistencia al corte
La resistencia al corte de las arenas proviene de la interacción de mecanismos físicos que
se puede reducir de manera simplista a una contribución de la fricción mineral entre
partículas y una contribución de la interferencia entre las trayectorias de las partículas
que se desplazan. La fricción se caracteriza mediante el coeficiente de fricción µ o del
ángulo de fricción mineral
φµ = atan [µ] (2.4)
φµ depende de la composición mineralógica de las partículas, de sus contaminantes
superficiales y, en menor medida, de la composición del fluido intergranular.
La interferencia es la restricción cinemática al movimiento de partículas vecinas;
depende de la forma de las partículas, de su distribución granulométrica y de la densidad
del empaquetamiento. En una muestra suelta, las partículas en movimiento encuentran
vacíos, caen en ellos, y el volumen del cuerpo disminuye. En una muestra densa las
partículas se montan unas sobre otras, con trayectorias que se apartan de la dirección
de desplazamiento macroscópico y que se traducen en dilatancia. El cambio de volumen
consume trabajo que debe ser aportado por la acción externa, por lo que la dilatancia
inducida por la interferencia intergranular constituye una componente de la resistencia
del material.
2.3.2 Modelo de los dos bloques
La resistencia al deslizamiento entre dos cuerpos en contacto plano es proporcional a
la fuerza de compresión normal al contacto. En el esquema de la figura 2-2, el bloque
superior está en equilibrio si la tensión de corte τn cumple con
τn < τ f = σn tan[φµ
](2.5)
18
P
S
nτ
φψ R
nσψ
1A a x=
a
Figura 2-2: El modelo de los dos bloques para el análisis de la resistencia friccional.
Si se definen las tensiones τ = S/A y σ = P/A referidas al área horizontal A del
bloque superior, el equilibrio ahora requiere que
τ < s = σ tan [φ] (2.6)
donde s es una medida convencional de la resistencia al corte y
φ = φµ + ψ (2.7)
Por lo tanto, la resistencia al corte s depende del material, caracterizado por φµ, de
la tensión “normal”2 σ , y de la cinemática del movimiento, caracterizada por ψ.
2.3.3 Ensayo triaxial
En un ensayo triaxial típico, la tensión principal mayor σ1 crece de manera monotónica
mientras se mantienen constantes σ2 = σ3 hasta que se produce la falla del material. De
manera idealizada puede considerarse que, sobre las superficies de deslizamiento que se
2Tensión normal al plano horizontal. Obsérvese que la tensión normal al plano de deslizamiento reales σn y no σ.
19
producen en el interior de la muestra, el conjunto de partículas presenta la cinemática
que se muestra en la figura 2-3. En esta concepción, el ángulo ψ “aparente” depende del
progreso de la deformación y de la topología de los contactos en cada instante.
P
S
ψ
x∆
y∆
[ ] y x tan ψ∆ = ∆
superficie de
deslizamiento
Figura 2-3: El modelo de los dos bloques en un conjunto de partículas.
En el momento de la falla, el cociente σ1/σ3 alcanza su valor máximo, que puede ser
calculado con la ec. (2.5) y que es igual a
(σ1/σ3)f = Nmax (2.8)
donde
Nmax = tan2 [π/4 + φ/2] (2.9)
Si se ejecutan ensayos distintos del ensayo triaxial de compresión, (σ1/σ3)f alcanza
valores distintos de Nmax, lo que prueba que la ec. (2.8) no configura un criterio de falla
exacto para arenas. Aunque en la geotecnia clásica se define comúnmente un “ángulo de
fricción interna” particular para cada uno de los ensayos que se ejecutan, en este trabajo
φ se define únicamente para el ensayo triaxial de compresión.
2.3.4 Estado crítico
Las deformaciones volumétricas debidas al reordenamiento de partículas tienden a cesar
a medida que los contactos se disponen de forma estadísticamente paralela a la dirección
20
0.80
0.85
0.90
0.95
1.00
10 100 1000 10000[ ]kPap
ce
Arena Toyoura
Figura 2-4: Relación de vacíos crítica en función de la presión para arena Toyoura [154].
macroscópica de movimiento. El estado de densidad resultante fue descripto por primera
vez por Casagrande [19][20] como estado crítico y está caracterizado por la relación de
vacíos crítica ec y por el ángulo de fricción interna crítico φc. En la analogía de los bloques,
se cumple que φc = φµ, mientras que para un material compuesto por partículas φc > φµ.
φc es aproximadamente independiente de la presión. En cambio, ec depende fuertemente
de la presión media debido a los procesos de ruptura de partículas que ocurren a altas
presiones. Como ejemplo, en la figura 2-4 se muestra la dependencia de la relación de
vacíos crítica con respecto a la presión para arena Toyoura [154]. Cada punto es un
ensayo. Aunque para cada punto se alcanzó un valor de ec determinado, en todos los
ensayos se midió φc ≈ 31◦ [154].
La contribución de la fricción y de la interferencia se comprende mejor con la ayuda
de la siguiente experiencia: Una caja llena con arena se inclina lentamente hasta que los
granos comienzan a rodar por la superficie del talud. Se reduce entonces la inclinación
de la caja hasta que las partículas se detengan. El ángulo máximo alcanzado por la
caja es φ, mientras que el ángulo que produce el cese el movimiento es φc. La diferencia
entre ambos es el efecto del empaquetamiento, mayor cuanto mayor es la densidad de
21
Figura 2-5: Nmax, Ncv, Nc y el cambio de volumen de una arena suelta y otra densa.
la arena que llena la caja. Puesto que las partículas al rodar pierden completamente
sus posiciones relativas, φc es independiente de la densidad inicial de la arena y es una
verdadera propiedad del material.
Existen dos estados en los que el material se deforma a volumen constante durante la
carga monotónica: i) uno transitorio cuando el material cambia de contractivo a dilatante,
en el que se define el ángulo de fricción a volumen constante φcv; y ii) el estado crítico
final, en el que se define φc. En la figura 2-5 se muestra la curva σ1/σ3 − ε1 y εv − ε1de ensayos triaxiales efectuados sobre arenas sueltas y densas, junto con los puntos que
definen a Nmax y a
Ncv = tan2 [π/4 + φcv/2] (2.10)
Nc = tan2 [π/4 + φc/2] (2.11)
En la figura 2-6 [123] se aprecia la relación entre φc, φcv, φµ, φ y la porosidad n. Existen
numerosas fórmulas, no todas consistentes entre sí, que permiten predecir φµ a partir de
22
φc (ej. [8]). Horne [55][56] encontró una relación funcional teórica entre φµ y φc que, para
arenas naturales con 24◦ < φc < 40◦, es numéricamente igual a φc =
89
(φµ + 8
◦).
0%rD =
cφ incremento de
debido a la d ilatancia
φ
cvφreacomodamiento
de partículas
µφfricción
intergranular
1 0 0 %rD =n
ángulo de fricción a volumen constante
ángulo de fricción crítico
ángulo de fricción intergranular
ángulo de fricción
cv
c
µ
φφφφ
φ
Figura 2-6: Influencia de la densidad relativa sobre la fricción, tomada de [123].
2.3.5 Estado de deformación constante
El estado de deformación constante es, para la compresión triaxial no drenada, el
equivalente del estado crítico para la compresión drenada [19][20]. Fue definido por
Poulos [115]:
“El estado constante de deformación ... es ese estado en el cual la
masa se deforma continuamente a volumen constante, tensión normal efectiva
constante, tensión de corte constante y velocidad constante. El estado
constante de deformación se logra solamente después que la orientación de
todas las partículas ha alcanzado estadísticamente una condición de estado
23
constante y despues que toda la rotura de partículas, si existe, se complete,
de manera que la tensión de corte necesaria para continuar la deformación y
la velocidad de deformación permanecen constantes”.
El estado constante de deformación es el elemento básico que se emplea en la
actualidad para la caracterización del comportamiento de arenas sueltas bajo carga cíclica
[21][22][28][29][58][59][115][154]. Sin embargo, los estados constante de deformación y
crítico se deben a fenómenos físicos comunes [104][105], por lo que la caracterización del
comportamiento de arenas en compresión no drenada puede hacerse mediante ensayos
drenados [104][105][106].
Durante la compresión triaxial drenada, el material pasa por un estado transitorio
en el que εv = 0 y en donde σ1/σ3 = Ncv. El mismo estado transitorio se aprecia en la
compresión triaxial no drenada - en la que siempre es εv = 0 - como un punto estacionario
en el que p = 0 kPa. En ese punto puede definirse un ángulo de transformación de fase
φpt tal que
σ1/σ3 = tan2[π/4 + φpt/2
](2.12)
Al estado caracterizado por φpt se lo denomina estado de cuasi deformación constante.
En la figura 2-7 [59] se presenta el comportamiento esquemático de una arena ensayada
a compresión triaxial no drenada, junto con las definiciones de estado de deformación
constante y estado de cuasi deformación constate.
2.3.6 Resistencia máxima y presión media
Si se toma una muestra de arena densa, se la somete a una presión hidrostática reducida
y se le aplica una tensión de corte monotónicamente creciente, los granos se deslizan
unos sobre otros y el material falla con un gran desarrollo de dilatancia (por lo que
φ≫ φc). Luego de una cierta distorsión, el arreglo de partículas se desarma por completo
y eventualmente se alcanza el estado crítico.
Si se repite la experiencia bajo una presión hidrostática muy elevada, se observa que
24
máximo local
estado de deformación constante
estado de cuasi
deformación constante
dσ dσ
1εp
Figura 2-7: Definición de estado constante y cuasi estado constante de deformación [59].
las tensiones de contacto fracturan granos, aristas y salientes, que la muestra rompe sin
dilatar y que exhibe un ángulo de fricción interna máximo φ = φc [20][89][123][125].
Existe una relación (casi) única entre densidad inicial, presión media y dilatancia, y una
presión que suprime la dilatancia para cualquier densidad [7][35][144][156]. La figura 2-8
muestra resultados de ensayos de compresión triaxial de arena de Sacramento [79]. Se
aprecia que el aumento de σ3 produce una disminución del ángulo de fricción interna
máximo, medido a través del cociente (σ1/σ3)f .
La dependencia de φ con respecto a la presión es tenida en cuenta en la mayoría de
los análisis modernos de presas de materiales sueltos [24][25][26][41] [80][101][107][126]
[133][138][141][144]. Las expresiones preferidas son de la forma
φ = φ0 −∆φ log [p/patm] (2.13)
donde φ0 y ∆φ son parámetros y patm es la presión atmosférica.
Bolton [11][12] sistematizó algunos resultados de arenas de distinta procedencia y
características físicas e introdujo una definición de densidad relativa dependiente de la
25
1ε
[ ]3 MPaσ
dilatancia positiva
dilatancia negativa
ε
vε
todas las muestras densas
100%rD =
[ ]3 MPaσ
1
3
σσ
Figura 2-8: Ensayos triaxiales de arena Sacramento densa a diferentes presiones deconfinamiento, adaptado de Lee [78] [79].
26
presión cuya expresión es
IR = Dr (Q− ln [p/ (1 kPa)])−R (2.14)
donde Q y R son parámetros de ajuste. A partir de esta definición, Bolton propuso la
siguiente relación funcional
φ = φc + 3◦ IR (2.15)
Been y Jefferies [4] definieron el parámetro de estado
Ψ = e− (ess − λ ln [p/patm]) (2.16)
donde λ y ess son parámetros materiales. Establecieron que dos muestras de arena con
un mismo valor de Ψ presentan el mismo comportamiento si son ensayadas al corte no
drenado. Si Ψ ≃ 0 el material está en estado crítico; siΨ > 0 la muestra contrae; si Ψ < 0
el material dilata [4][5][61]. El término (ess − λ ln [p/patm]) caracteriza la dependencia
de la relación de vacíos crítica con respecto a la presión.
2.3.7 Estado triaxial general
El criterio de falla de Mohr-Coulomb ignora la contribución de la tensión intermedia σ2
en la resistencia del material. Sin embargo, se sabe que (σ1/σ3)f es sistemáticamente
superior en deformación plana que en condiciones triaxiales, para un mismo estado de
densidad y de presión media [34][44][63][70][82][90][109][124][142]. Se acepta que éste
fenómeno se produce porque en estados distintos al de compresión triaxial las partículas
experimentan una mayor restricción al movimiento relativo, lo que se traduce en un mayor
trabajo de deformación y en una mayor disipación de energía por fricción intergranular.
Lade y Duncan [70][71] y Lade [72][73][74] presentaron el criterio de falla
I31 − κ I3 = 0 (2.17)
27
donde I1 e I3 son los invariantes canónicos de σ y κ es un parámetro material. Lade
y Duncan [70][71], Lade y Kim [75][76] y Kim y Lade [62] también mostraron que la
deformación plástica es aproximadamente normal a la traza deviatórica de la ec. (2.17)
. Matsuoka y Nakai [90][91][99] presentaron la expresión
I1I2 −Mf I3 = 0 (2.18)
donde el parámetro material Mf , calibrado en un ensayo de compresión triaxial, vale
Mf = 9 + 8 tan2 [φ] (2.19)
Las ecs. (2.17) y (2.18) son muy similares. La diferencia más importante es que la
ec. (2.17) produce un valor de (σ1/σ3)f mayor en extensión triaxial que en compresión
triaxial, mientras que la ec. (2.18) predice que el cociente (σ1/σ3)f en extensión triaxial
es igual al de compresión triaxial. La comparación entre los criterios de Tresca, Von
Mises, Mohr—Coulomb y Matsuoka—Nakai puede verse en la figura 2-9 [91].
De acuerdo a esta interpretación, la falla de los materiales friccionales ocurre cuando
se alcanza la resistencia al corte en el plano espacial movilizado (SMP ), punto P de la
mitad inferior de la figura 2-9, análogo al plano octaédrico de Von Mises que se define en
el punto P de la mitad superior de la figura.
Las tensiones en el SMP se calculan con las expresiones
σSMP =3I3I2
(2.20)
τSMP =
√I1I3I2
−(3I3I2
)2(2.21)
y cumplen con la relación
τSMP =
√8
9σSMP tan [φSMP ] (2.22)
28
maxτ
τ
1σ
3σ2σ
3σ σ1σ
octτ
2σ o c tσ
,oct octσ τTresca
Von Mises2P
1P3P
P
Tresca
Von Mises1σ
2σ 3σ
σ1σ3σ 2σSMPσ
1P3P
2PP
,SMP SMPσ τ
1στ
SMPτ Mohr-Coulomb
Matsuoka-Nakai
Mohr-Coulomb
Matsuoka-Nakai
Figura 2-9: Comparación entre los criterios de Tresca, Von Mises, Mohr-Coulomb yMatsuoka-Nakai [91].
lo que permite apreciar la similitud formal que existe entre el criterio de Matsuoka—Nakai
y el criterio de Mohr—Coulomb.
2.4 Rigidez
La rigidez macroscópica de un conjunto de partículas obedece a una compleja interrelación
entre granos, que incluye su compresión elástica, deslizamientos y ruptura. Para un
mismo estado tensional de mesoescala, la cantidad de contactos activos crece con la
densidad relativa mientras que las fuerzas que actúan en los contactos y las tensiones
dentro de las partículas disminuyen. Por todo ésto, las arenas densas tienen mayor
29
rigidez y un rango de comportamiento “elástico” mayor que las arenas sueltas.
2.4.1 Elasticidad y contacto Hertziano
Dos esferas elásticas de radio R,módulo de Young E y módulo de Poisson ν, comprimidas
entre sí por una fuerzaN , se deforman como se muestra en la figura 2-10. El acercamiento
u entre los centros de las esferas es
u = 2
(3
4
(1− ν2)NE√R
) 2
3
(2.23)
mientras que la rigidez volumétrica asociada es
K = CK
(p
pref
)m
pref (2.24)
donde CK es un parámetro que depende de E, ν, R y de la densidad del empaquetamiento,
p es la presión media, pref es una presión de referencia para adimensionalización y
m = 1/3. La misma teoría predice que cuando un cono indenta una superficie roma,
el exponente de la ec. (2.24) es m = 1/2.
La rigidez volumétrica de mesoescala está formada por la interacción de muchas
partículas de diferentes formas y tamaños, y puede calcularse de forma aproximada con
la ec. (2.24). En este caso, CK y m deben ser determinados experimentalmente. Se ha
indicado [163] que 1/2 ≤ m ≤ 1/2, lo que sugiere que las condiciones de contacto esfera— esfera y cono — plano constituyen casos extremos del comportamiento de los materiales
reales.
Pestana estudió la influencia de e sobre CK y propuso la expresión [109][110]
CK = cb1 + e
e(2.25)
donde cb es un parámetro material.
La rigidez al corte G se origina en los mismos fenómenos físicos que controlan a la
30
RN N
2 -R u
Figura 2-10: Compresión de dos esferas elásticas.
rigidez volumétrica K. Para materiales isótropos, está bien establecido que G puede
caracterizarse mediante expresiones del tipo
G ∝(σiσj
pref
)m
pref (2.26)
donde σi y σj son las tensiones actuantes en el plano de propagación y polarización de la
onda de corte [127][128]. Dentro de las limitaciones de la elasticidad isotrópica, son muy
empleadas las formas [49][50][66][118]
G = CG
(p
pref
)m
pref (2.27)
en la que CG es constante o función de la relación de vacíos. De las muchas expresiones
propuestas para CG, la más utilizada es [49][50][66][118]
CG = cs(ce − e)21 + e
(2.28)
con parámetros materiales cs y ce.
Santamarina y Cascante [127] revisaron teorías que permiten calcular la relación
de Poisson elástica de empaquetamientos regulares e irregulares de esferas de tamaño
uniforme. Mostraron que la relación de Poisson ν del empaquetamiento depende de
la relación de Poisson νm del material que forma las esferas y de la topología del
empaquetamiento. Los valores de ν que se obtienen con νm = 0.30, un valor típico
31
para minerales que forman arenas, van desde ν = 0.00 para un empaquetamiento cúbico
a ν = 0.05 para un empaquetamiento cúbico centrado en el cuerpo [127][128]. Las arenas
naturales pueden tener empaquetamientos más densos que el cúbico centrado en el cuerpo,
por lo que pueden esperarse valores máximos ν > 0.10.
2.4.2 Rango elástico y disipación de energía
Las arenas son materiales disipativos aún para pequeños rangos de deformación. El
fenómeno físico que justifica este comportamiento puede explicarse con el siguiente
experimento ideal.
Se toma un empaquetamiento desordenado y denso de esferas elásticas y se lo somete
a una presión hidrostática. Se observan algunos deslizamientos de contactos, que cesan
cuando las tensiones de corte locales son iguales a la resistencia al corte dada por la
fricción entre esferas. Los contactos que deslizan quedan en “falla” por corte mientras
que otros contactos quedan con fuerzas tangenciales menores a su resistencia friccional.
Analizado en mesoescala, el material está sometido a una presión hidrostática, y por lo
tanto libre de tensiones de corte; sin embargo, tiene deformaciones irreversibles dadas
por los mencionados deslizamientos de partículas.
Se analiza ahora el comportamiento de uno cualquiera de los contactos que no está
en falla por corte como el que se esquematiza en la figura 2-11 [39]. El contacto tiene
una fuerza normal N y una fuerza tangencial
T < N tan[φµ
](2.29)
La superficie del contacto, según la teoría de contacto de Hertz, es un círculo de radio
a =3
√3
4
1− ν2E
NR (2.30)
La distribución de tensiones normales dentro del contacto sigue una ley parabólica,
32
a
N N
T
T
límite friccional
solución elástica
τnσ 0δ =
0δ >
0δ >
Figura 2-11: Distribución de tensiones normales y tangenciales en el contacto de dosesferas elásticas [39].
con un máximo en el centro y un valor nulo en el perímetro. La componente T paralela al
contacto produce tensiones de corte τ , infinitas en el perímetro según la teoría de contacto
elástico. La resistencia friccional limita el valor en cada punto a τ f = σ tan[φµ
], por lo
que el contacto desarrolla una corona de deslizamiento relativo (τ = τ f) y un núcleo en
estado elástico (τ < τ f ).
Debe notarse que no se requiere que haya un movimiento de cuerpo rígido entre las
dos esferas para que se produzca la corona de plastificación en el contacto, y con ella
la disipación de energía en forma de trabajo plástico. Por lo tanto, se concluye que la
compresión hidrostática de esferas rígidas es un fenómeno disipativo para cualquier escala
de deformación.
33
2.4.3 Módulo de Young inicial
El módulo de Young fue el primer parámetro que se usó para caracterizar la rigidez de
los geomateriales. Parece ser que la primera caracterización del módulo de Young para
arenas se debe a Scheidig, y se conoce a través de la figura 2-12, incluida en el libro de
Terzaghi y Peck [148], en el que se relacionan la rigidez inicial en el ensayo triaxial, la
densidad relativa y la presión de confinamiento σ3. La rigidez inicial en el ensayo triaxial
se define como
Ei =dσd
dε1
⌋
ε1=0
(2.31)
donde
σd = σ1 − σ3 (2.32)
y ε1 es la deformación principal mayor. La figura 2-12 muestra tres características de las
arenas: i) Ei crece con Dr; ii) Ei crece con p; y iii) la dependencia entre Ei y p disminuye
a medida que crece Dr. Estas conclusiones coinciden con las de muchos otros trabajos
posteriores (ej. [35][64][104][119]). Es sorprendente la precisión de la figura 2-12, sobre
todo si se tiene en cuenta que fue escrita cuando las técnicas de medición de tensiones y
deformaciones en el laboratorio eran muy rudimentarias.
Jambu [60], en base a trabajos anteriores de Schultze y Moussa [132], desarrolló la
expresión
Ei = CE
(σ3patm
)m
patm (2.33)
donde Ei es el módulo de Young inicial, mientras que CE y m son parámetros materiales.
Núñez [104] [105] estableció las correlaciones
100 + 2000D2r ≤ CE/pref ≤ 100 + 2000Dr (2.34)
1/3 + 2/3 (Dr − 1)2 ≤ m ≤ 1− 2/3Dr (2.35)
En la definición de Ei se incluyen tanto deformaciones reversibles como irreversibles,
lo que explica porqué el exponente m puede caer por fuera de los límites dados por la
34
Figura 2-12: Relación entre el módulo de deformación, la densidad y la presión deconfinamiento (aporte de Scheidig, publicado en [148]).
teoría de contacto elástico.
La rigidez “elástica”, medida en descarga, es mayor que Ei. Fue caracterizada por
Duncan [41] mediante la expresión
Eur = Kur
(σ3patm
)n
patm (2.36)
donde Kur es un parámetro. Al respecto, Duncan [41] escribió:
“El valor deKur es siempre mayor a CE3 (para primera carga). Kur puede
ser 20% mayor que CE para suelos rígidos como (por ejemplo) arenas densas.
Para suelos compresibles como (por ejemplo) arenas sueltas, Kur puede ser
tres veces más grande queCE. El valor del exponente n es siempre muy similar
para primera carya y descarga, y en las relaciones hiperbólicas se asume que
es el mismo.”
3La notación se ajustó por consistencia con el resto del documento.
35
1dδ
3σ
1σ1 1 1
v A dδ δ=
3 3 3v A dδ δ=
3dδ
Figura 2-13: Dispositivo esquemático para ilustrar la teoría de tensión-dilatancia de Rowe[123] [124] [125].
2.4.4 Teoría de tensión — dilatancia de Rowe
La teoría de tensión — dilatancia de Rowe [123][124][125] relaciona el estado de tensiones
y la velocidad de dilatancia. Puede ser introducida de manera simplificada con el auxilio
de la figura 2-13, que muestra un sistema de dos pistones que encierran un material de
prueba.
El primer material de prueba es un fluido incompresible. Puede observarse que el
desplazamiento del pistón superior ´hacia abajo es posible únicamente si el pistón lateral
se desplaza hacia afuera. En este caso, el trabajo incremental efectuado por las tensiones
principales mayor y menor es
win = ‖σ1δv1‖ (2.37)
wout = ‖σ3δv3‖ (2.38)
donde δv1 = A1δd1 y δv3 = A3δd3. Como no se disipa trabajo en la deformación, se
verifica que
win = wout (2.39)
36
Como el fluido es incompresible, δv3 = δv1, por lo que σ1 = σ3. El ensayo se repite,
reemplazando al fluido incompresible por un fluido dilatante. Por continuidad,
δv3 = (1− ∂v/∂d1) δv1 (2.40)
donde v = v1+ v3 y ∂v/∂d1 es un coeficiente de dilatancia. La aplicación de la ec. (2.39)
arroja
σ1 = (1− ∂v/∂d1) σ3 (2.41)
lo que muestra que el flujo de un material dilatante sólo es posible si σ1 > σ3, aún si el
material es un fluido sin resistencia al corte.
Si el dispositivo se llena con partículas rígidas e irrompibles de arena en el estado
crítico, el flujo isocórico ocurre para la relación de tensiones
σ1 = Ncσ3 (2.42)
Si se multiplica el primer término de la ec. (2.42) por δv1 y el segundo término por
δv3 y se tiene en cuenta que el flujo isocórico implica δv3 = δv1, queda
win = Ncwout (2.43)
que es la expresión general de la teoría de tensión dilatancia de Rowe [123][125].
Para un material friccional dilatante, la condición de continuidad δv3 = δv1 se
reemplaza por la ec. (2.40), por lo que
σ1 = (1− ∂v/∂d1)Ncσ3 (2.44)
que es la forma reducida de la teoría tensión-dilatancia de Rowe para flujo 2D [123][125].
37
Si se inserta la ec. (2.8) en la ec. (2.44) se obtiene
Nmax = (1− ∂v/∂d1)maxNc (2.45)
lo que muestra que el ángulo de fricción interna φ que define a Nmax depende del ángulo
de fricción interna crítico φc que define a Nc y de la cinemática del flujo de partículas
caracterizada por ∂v/∂d1 [38][123][125].
Como se muestra en la figura 2-5, el punto de la curva σ1/σ3−ε1 en el que ∂v/∂d1 = 0ocurre para σ1/σ3 = Ncv por lo que la ec. (2.43) debe ser reemplazada por
win = Ncvwout (2.46)
La ec. (2.46) se puede particularizar para los ensayos de compresión triaxial, extensión
triaxial y deformación plana [125]. Quedan las expresiones
compresión triaxial σ1ε12σ3ε3
= Ncv
extensión triaxial 2σ1ε1σ3ε3
= Ncv
deformación plana σ1ε1σ3ε3
= Ncv
(2.47)
El fundamento mecanicista de la eq. (2.44) puede seguirse con la ayuda de la figura
2-14, en la que un arreglo de partículas rígidas está sometido a tensiones exteriores [123].
Como lo explica el modelo de los bloques, el desplazamiento relativo de dos partículas
requiere que
P1 = P3 tan[φµ + β
](2.48)
La dirección del deslizamiento inicial está dada por el ángulo β, por lo que el cambio
de volumen es
∂εv/∂ε1 = 1− tan [α] tan [β] (2.49)
La geometría del empaquetamiento permite calcular el cociente entre tensiones
38
y
x
α
β
1Py
R
β
µφ
3P
x
1P3P
Figura 2-14: Esquema de Rowe que indica la relación entre el plano de deslizamiento yla dirección del deslizamiento.
principales durante el deslizamiento como
σ1σ3=P1/d1P3/d3
= tan[φµ + β
]tan [α] (2.50)
Rowe [123] postuló que la razón entre el trabajo recibido P1dy y el ejercido contra el
confinamiento P3dx es un mínimo
min
[E]=P1dy
P3dx=tan[φµ + β
]
tan [β](2.51)
lo que permite calcular β = π/4 − φµ/2 que, reemplazado en las ecs. (2.50) y (2.49)
arroja
σ1/σ3 = tan2[45 + φµ/2
](1− ∂εv/∂ε1) (2.52)
donde
Nµ = tan2[45 + φµ/2
](2.53)
39
Otra vez, se encuentra que la teoría tiene mayor capacidad predictiva si se reemplaza
Nµ por Ncv. Una discusión sobre la relación entre φcv y φµ puede verse en [55][56].
La teoría tensión-dilatancia puede ser empleada en el marco de la teoría de la
plasticidad. Por ejemplo, para estado triaxial de compresión, la teoría de tensión-
dilatancia es equivalente a [36][155]
εpv = −
σ1 − σ3 − (σ1 + σ3) sin [φcv]
σ1 + σ3 − (σ1 − σ3) sin [φcv]γp
V (2.54)
donde εpv es la deformación plástica volumétrica y
γpV = 2ε
p1 − εp
v (2.55)
es una medida de distorsión plástica. Existen muchas otras expresiones que relacionan
εv con la relación σ1/σ3 (ej. [27][131][144]). Sin embargo, la teoría tensión-dilatancia
de Rowe permite una mejor comprensión del comportamiento del material, a la vez que
tiene una adecuada capacidad predictiva, al menos para los ensayos de laboratorio más
utilizados.
2.4.5 Cambio de volumen en descarga
Durante la descarga también se producen cambios de posición entre partículas que a
su vez producen cambios de volumen irreversibles. Como todo proceso irreversible, el
cambio de volumen en descarga está asociado a disipación de energía [23][61][108]. Por lo
tanto, existe una restricción termodinámica que limita el cambio de volumen irreversible
debido a un cambio de oblicuidad tensional. En efecto, el trabajo plástico instantáneo es
w = σ : εp = pεpv + s : ε
pd (2.56)
40
y tiene que ser positivo. Por lo tanto, durante la descarga tiene que cumplirse que
εpv + r : ε
pd > 0 (2.57)
donde εpd es la deformación plástica deviatórica. Esta restricción se observa en los
resultados experimentales. En la figura 2-15 se muestra un resultado típico de un ensayo
triaxial cíclico no drenado, en el que se aprecia que los diferentes ciclos producen una
acumulación de presión neutra y, por lo tanto, una disminución de la presión efectiva
[136]. Puede observarse que todas las ramas de descarga son puramente contractivas,
mientras que las ramas de carga cambian de contractivas a dilatantes para una dada
oblicuidad tensional.
-20
-10
0
10
20
30
0 10 20 30 40
[ ]kPa
dσ
[ ]kPap
Ensayo triaxial cíclio no drenadoArena NevadaD r =40%
Figura 2-15: Ensayo triaxial no drenado cíclico de arena Nevada [2], que prueba que lasarenas tienden a contraer en descarga.
2.4.6 Reducción de la rigidez secante por deformación
La evolución de los mecanismos disipativos de los materiales friccionales se estudia a
través de las curvas de reducción de rigidez al corte secante o del módulo de deformación
41
secante
Gs = τ/γ (2.58)
Es = σd/ε1 (2.59)
en ensayos de compresión triaxial, deformación plana, torsión anular y de columna
resonante. Las curvas se normalizan siempre con respecto a la rigidez máxima, que
se asume en este trabajo igual a la rigidez elástica G o E, respectivamente, porque se
aprecia que esta normalización elimina el efecto de la relación de vacíos y produce una
(casi) única curva Gs/G− γ o Es/E − ε1 [41][49][50][58][66][118].
Se aprecia que Gs, medido en ensayos de carga cíclica, depende del nivel máximo
de distorsión alcanzado [49][50][58][163]: para distorsiones menores a 10−5 − 5·10−4 elcomportamiento es (casi) reversible; en el rango 10−4−10−2 se desarrollan deformacionesirreversibles y una fuerte reducción deGs; si la distorsión es mayor. sobreviene la ruptura,
con eventual ablandamiento y localización de deformaciones. La figura 2-16 muestra la
relación Gs/G− γ a partir de ensayos triaxiales cíclicos de arena Toyoura, adaptados de
[66].
La curva llena de la figura 2-16 corresponde al modelo hiperbólico para carga cíclica
[49][50]. Su expresión es
Gs =1
1/G+ γ/τ f
= G (1− τ/τ f ) (2.60)
donde γ es la distorsión y τ f es la tensión de corte de falla. La rigidez tangente
Gt = ∂τ/∂γ correspondiente, expresada como función del estado tensional, es
Gt = G (1− τ/τ f )2 (2.61)
La evolución de Gs/G puede estudiarse también en un gráfico Gs/G −τ/τ f , como se
muestra en la figura 2-17 [142]. Se aprecia que la curva Gs/G−τ/τ f es muy distinta para
42
aγ
aγ
deformación reversible
con disipación
τsG
sG
G
e=0.640-0.649
e=0.696
e=0.742
e=0.793
0
Arena Toyoura
=100KPa; =10Nσ
falla
deformación
irreversible0.80
0.60
0.40
0.20
6
10− 5
10− 4
10− 3
10− 2
10−
Figura 2-16: Relación entre el módulo secante y el módulo inicial, para 10 ciclos de cargay para diferentes relaciones de vacíos. Arena Toyoura, de [66].
carga monotónica y cíclica. Se entiende que la causa de esta diferencia es que durante la
carga monotónica los contactos entre partículas se reorientan progresivamente hasta que
se genera una superficie de falla, mientras que en carga cíclica se producen inversiones
que regeneran el núcleo elástico de algunos contactos y aumentan la rigidez macroscópica.
El acoplamiento entre el cambio de las tensiones de corte y el cambio de volumen
se traduce en los siguientes fenómenos: i) una carga alternativa moderada produce
densificación y rigidez creciente; ii) una carga alternativa cercana a la falla produce
dilatancia y falla con ablandamiento [165][166]. La figura 2-18 muestra el efecto de la
amplitud de carga alternativa en el comportamiento del material.
Kondner y Zelasko [67] y Kondner [68] propusieron una relación tensión — deformación
hiperbólica de expresión
σd =ε1
a+ b ε1(2.62)
43
/ fτ τ
τ0
Arena Toyoura
normalmente consolidada
167 ; 0.79; 65%KPa e Drσ = = =
carga cíclica
carga
m onotón ica
sG
G
γ
0 .40
0 .20
0 .60
0 .80
0 .2 0 0 .40 0 .60 0 .80
Figura 2-17: Módulo secante en ensayos monotónicos y cíclicos de arena suelta, según[142].
en la que a y b son parámetros. Luego, Duncan y Chang [40] y Duncan [41] desarrollaron
esta expresión y la convirtieron en una ecuación constitutiva hipoelástica de enorme
difusión y utilización en ingeniería práctica. En condiciones triaxiales, la expresión de
Duncan - Chang es
σd =ε1
1/Ei +Rf ε1/σdf
(2.63)
donde
σdf = (σ1 − σ3)⌋falla (2.64)
y Rf < 1 es un parámetro. Los módulos de Young secante y tangente correspondientes
son
Es = Ei (1− Rf σd/σdf) (2.65)
Et = ∂σd/∂ε1 = Ei (1− Rf σd/σdf)2 (2.66)
44
amplitud reducida
amplitud grandedσ
monotónico
amplitud
media
1ε
Figura 2-18: Comportamiento de una arena densa bajo carga triaxial monotónica y cargascíclicas de distintas amplitudes.
En la figura 2-19 se muestra la hipérbola de Kondner con los parámetros de Duncan—
Chang, junto con una curva σd − ε1 típica. σu = σdf/Rf es la asíntota horizontal de la
hipérbola.
Núñez [104][105][106] desarrolló expresiones similares y las aplicó a suelos locales en
una variedad de problemas prácticos. Debido a esta circunstancia, el modelo hiperbólico
es bien conocido en Argentina, y existe experiencia en la calibración de sus parámetros
para muchos problemas de ingeniería geotécnica aplicada.
El modelo hiperbólico es universalmente empleado para caracterizar arenas. Entre
las muchas correlaciones disponibles para sus parámetros (ej. [41][104][130][133][149]) es
de particular interés, desde un punto de vista conceptual, la de Trautmann y Kulhawy
[69][149] que propusieron
CE = 300 + 900φ− 25◦45◦ − 25◦ (2.67)
La idea implícita en la ec. (2.67) es que cualquier cambio de presión y/o densidad
que produzca un efecto en φ producirá un efecto similar en CE.
45
resultado
experimental
hiperbola de
Kondner
iE
1
1
1
df
i df
R
E σ
εσε
=+
dσ
dfσuσ
1εfε
Figura 2-19: Modelo hiperbólico y su comparación con el resultado de un ensayo triaxialtípico [40][67][68].
2.4.7 Módulo secante E50
Otra medida de rigidez elastoplástica ampliamente utilizada en geotecnia práctica es el
módulo secante E50, definido como
E50 =∂σd
∂ε1
⌋
σd/σdf=0.5
(2.68)
E50 puede ser adimensionalizado si se lo divide por el módulo de Young elástico E.
El cociente E50/E resultante depende de la densidad relativa y, en menor medida, de la
presión de confinamiento [28][29][84][85][164].
46
2.4.8 Compresión proporcional
La compresión proporcional consiste en el incremento de las tensiones principales mientras
se mantienen constantes las relaciones σ1/σ2 y σ1/σ3. Un caso particular de compresión
proporcional es la compresión hidrostática, en la que se aplica un presión uniforme
exterior. Si el material es isótropo, una compresión hidrostática produce una deformación
puramente volumétrica. En los párrafos siguientes se trata únicamente con materiales
isótropos.
Si se toma una muestra de arena densa sometida a una presión p1 relativamente
baja y se incrementa p1 hasta un valor mayor p2, las fuerza en los contactos aumentan su
magnitud pero se mantienen su dirección. Sólo ocurren pequeñas deformaciones asociadas
a la compresión elástica de los granos, al deslizamiento de unos pocos puntos de contacto
y a la ruptura de algunos granos [64][65][119].
Una arena muy suelta (por ejemplo, con Dr < 20%), sometida al mismo ensayo, sufre
el colapso de grupos de partículas en equilibrio metaestable y el subsecuente deslizamiento
entre granos, por lo que la deformación volumétrica es mayor y más irreversible. Si la
presión p1 es elevada, es mayor la cantidad de granos que se rompe y la deformación
reológica de los contactos, por lo que se observan mayores deformaciones irreversibles y un
comportamiento menos acrónico [7][51][52][100]. A muy altas presiones, muchos granos
sufren fragmentación y la topología original de la muestra se destruye por completo. En
la figura 2-20 se presentan tres ensayos de compresión isotrópica de la arena Toyoura [94]
que permiten apreciar el comportamiento macroscópico del material. Alcanzadas estas
altas presiones, son idénticos el comportamiento de una muestra originalmente suelta y el
de una originalmente densa: ambas recorren la denomindada curva de compresión límite
cuya expresión, propuesta por Pestana, es [109][110][111]
e =
(p
pr pref
)−ρ
(2.69)
con parámetros materiales pr y ρ.
47
0.50
0.60
0.70
0.80
0.90
1.00
100 1000 10000 100000
e = 0.83e = 0.77e = 0.59
Arena Toyoura
r refp p
[ ]kPap
e
Figura 2-20: Compresión isotrópica de arena Toyoura [94] y línea de compresión límite(ec. (2.69)) [109][110]
La dureza de grano pr determina la presión a la que el fenómeno de ruptura
de partículas controla el comportamiento general; depende fundamentalmente del
tamaño medio de las partículas y de la resistencia del material que forma los granos
[88][89][109][110]. La figura 2-21 muestra la evolución típica de una arena remoldeada en
compresión isotrópica y la curva de compresión límite [109].
El primer modelo constitutivo que tuvo en cuenta la compresión de geomateriales fue
el Cam Clay y luego el Cam Clay modificado [120][121][122][131]. Este último se expresa
en términos de tensiones principales como
(σd/p)2 +M (p− pc)
2 = 0 (2.70)
donde M y pc son parámetros materiales. El modelo Cam Clay reproduce bien el
comportamiento de arcillas blandas. Murata y otros [98][167][168] presentaron una
adaptación de la ec. (2.70) que se ajusta mejor a los resultados experimentales de arenas
48
r r efp p [ ]p MPa
e
curva limite:
ruptura masiva de granoscompresión elástica de granos
y deslizamiento de partículas
arena densa
arena suelta
1
10.01.00.1
0.10
1.00
100.0
ρ
Figura 2-21: Compresión isótropa de arenas. Significado de pr.[109][110]
consolidadas isotrópicamente bajo un amplio rango de tensiones. La expresión es
(σd/p)2 + η ln [p/pc] = 0 (2.71)
en la que η es un parámetro. La figura 2-22 muestra la forma de la ec. (2.71), su
comparación con la ec. (2.70), con otros criterios de fluencia y con algunos resultados
experimentales [98].
En la figura 2-22 puede también que las diferentes funciones de fluencia tienen formas
muy distintas, producto que algunas están calibradas para arcillas [120][121][122][131]
mientras que otras están calibradas para arenas [98][93][94][114]. En el caso particular
de la función de Poorooshasb, ésta es en realidad un potencial plástico determinado
experimentalmente, que incluye tanto plasticidad por corte como por compresión[114].
La ec. (2.71) utiliza el cociente σd/p como medida de la movilización de la resistencia
al corte, por lo que su traza deviatórica es una función facetada. Sfriso [135][136][137] se
presentó una modificación del criterio de Murata-Miura en el que se reemplaza σd/p por
49
d
cp
σ
/ cp p
Murata
Miura
arena Aio
densa
experimental
Pooroosshab
0.60
0.40
0.20
0.20 0.40 0.60 0.80
Cam-clay
modificado
Cam-clay
Figura 2-22: Criterio de Murata-Miura y su comparación con otras funciones de fluenciapara arenas, adaptado de [98].
τSMP/σSMP , que es un equivalente conceptual compatible con el criterio de Matsuoka —
Nakai. La función de fluencia que queda es
(τSMP
σSMP
)2+ η ln [p/pc] = 0 (2.72)
que puede reescribirse como
Fc = 3r2 − 3J3r
1− 12r2 + J3r
+ η ln [p/pc] = 0 (2.73)
La figura 2-23 muestra una vista frontal de la superficie dada por la ec. (2.73),
mientras que la figura 2-24 muestra una vista lateral.
50
Figura 2-23: Superficie de fluencia definida por la ec. (2.73) [135][136][137]. Vista frontal.
Figura 2-24: Superficie de fluencia definida por la ec. (2.73) [135][136][137]. Vista lateral.
2.4.9 Anisotropía y endurecimiento cinemático
Los materiales granulares naturales se forman por deposición de partículas en ambientes
áridos, húmedos o bajo agua. Los agentes de transporte pueden ser el viento, el agua, la
nieve, los eventos piroclásticos y las remociones en masa. Las partículas gruesas contienen
asociaciones de muchos minerales y formas muy variables, generalmente aplanadas4.
Arrastradas por el agua, se depositan con sus caras chatas en posición subhorizontal,
y son posteriormente desplazadas, deformadas y rotas por carga, intemperismo y nuevos
4La naturaleza hizo esto con las piedras para que podamos arrojarlas al agua en vuelo rasante...
51
procesos de transporte. Luego de esta somera descripción, no sorprende que los materiales
friccionales genéticamente isótropos sean escasos en la naturaleza.
La anisotropía genética no implica, necesariamente, anisotropía mecánica, aunque
ésta es frecuente en los materiales granulares gruesos. Hay opiniones diversas (p.ej.
[57][161][163][168][169]) acerca del efecto práctico de la anisotropía genética sobre las
propiedades mecánicas, y en la mayoría de los casos prácticos se la ignora. En los menos,
se la caracteriza como ortotropía transversal, con eje principal en la dirección de la
deposición [33][113][146].
Las ecuaciones 2.70, 2.71 y 2.73 son funciones de fluencia con endurecimiento
isotrópico. Para este tipo de funciones, la compresión del material en una dada
trayectoria de carga aumenta su dominio elástico en todas las direcciones de carga.
Este comportamiento es contrario a la evidencia experimental (ej. [161]), lo que
ha motivado el desarrollo de numerosos modelos de endurecimiento cinemático (ejs.
[31][32][87][61][97][111][162]). No existe aún consenso acerca de cuales son los parámetros
de estado que controlan la compresión de arenas. Los candidatos más firmes son las
medidas de ruptura de partículas y los tensores de orientación de cadenas de contactos.
En todos los casos, estos parámetros son de casi imposible medición experimental,
por lo que los modelos de endurecimiento cinemático con frecuencia tienen numerosos
parámetros con poco significado físico que impiden su aplicación práctica.
52
Capítulo 3
Modelización de la compresión
isotrópica
3.1 Introducción
En este capítulo se formula e implementa un modelo simple para compresión isotrópica.
El modelo tiene cuatro parámetros mecánicos y se basa en el modelo de compresión de
Pestana—Whittle [109][110]. Este desarrollo no tiene aplicación práctica por sí mismo.
El objeto de su presentación es mostrar, a través de un ejemplo unidimensional, algunos
aspectos de las consecuencias teóricas y numéricas de la dependencia de la rigidez respecto
de la presión que es común a todos los materiales friccionales.
3.2 Elementos del modelo
3.2.1 Cinemática
Se adopta la descomposición aditiva de las deformaciones infinitesimales
εv = εev + ε
pv (3.1)
53
donde εv es la deformación volumétrica total, εev es la deformación volumétrica elástica
y εpv es la deformación volumétrica plástica. εv > 0 implica disminución de volumen.
3.2.2 Variables de estado
Las variables de estado son la presión p, la relación de vacíos para presión nula e0 y la
presión de preconsolidación pc. La relación de vacíos para presión nula e0 se define como
la relación de vacíos que tiene el material si se congelan todos sus mecanismos plásticos y
se reduce la presión a cero. Su evolución depende únicamente de la deformación plástica
según
e0 = − (1 + e0) εpv (3.2)
3.2.3 Elasticidad
El módulo de rigidez volumétrico se calcula con [109][110]
K = cb1 + e0e0
(p
pref
)m
pref (3.3)
Como el problema es unidimensional, se puede integrar εev = p/K y calcular la relación
analítica p− εev que puede escribirse como
εev =
p
K (1−m) (3.4)
p =
(cb1 + e0e0
(1−m) εev
) 1
1−m
pref (3.5)
3.2.4 Función de fluencia y potencial plástico
La función de fluencia es
F = p− pc = 0 (3.6)
En un problema unidimensional el potencial plástico es trivial. La ec. (3.6) es también
el potencial plástico.
54
3.2.5 Ecuaciones de evolución
e0 evoluciona según la ec. (3.2). La ecuación de evolución de pc es
pc = Kpεp
v (3.7)
donde Kp es el módulo de endurecimiento plástico.
Para la deducción de Kp se define la función de estado
χ = p/pult (3.8)
donde pult se calcula con [109]
pult = e−1/ρ0 pr pref (3.9)
que es la ec. (2.69) invertida, en la que se reemplazó e por e0. La figura 3-1 repite la
información de la figura 2-20 junto con líneas correspondientes a tres valores de χ. En
particular, la línea χ = 1 es la descripta por la ec. (3.9).
0.50
0.60
0.70
0.80
0.90
1.00
100 1000 10000 100000
e = 0.83e = 0.77e = 0.59
Arena Toyoura
0 χ←
r refp p
[ ]kPap
e 1
2 3
1 3
χχχ
===
Figura 3-1: Compresión isotrópica de arena Toyoura, líneas de igual χ [94][109][110].
Puede apreciarse que cuando χ → 0 la ruptura de partículas es despreciable y el
55
comportamiento es “elástico”, o Kp →∞. En el otro extremo, cuando χ = 1 se alcanza
la línea de compresión límite dada por las ecs. (2.69) o (3.9).
La ecuación de evolución del material sobre la curva de compresión límite se obtiene
por derivación directa de la ec. (3.9) y reemplazo de e0 por la ec. (3.2). Luego de algo
de álgebra, queda
pc⌋χ→1 =1 + e0e0
pult
ρεp
v (3.10)
que, insertada en la ec. (3.7) permite calcular
Kpult =
1 + e0e0
pult
ρ(3.11)
donde Kpult es el valor de K
p cuando χ = 1.
El comportamiento del material sobre la curva límite es repetible y está
adecuadamente caracterizado para propósitos prácticos. Sin embargo, los fenómenos
físicos asociados a la transición entre el estado inicial (χ→ 0) y el estado último (χ→ 1)
son muy complejos y dependientes de muchas características físicas y mecánicas de los
granos, incluyendo su rugosidad superficial, esfericidad y distribución granulométrica. La
caracterización precisa de todos estos mecanismos requeriría la formulación de modelos
con un gran número de parámetros, para los que no existirían ensayos de laboratorio que
permitieran su determinación de manera confiable. La calibración de un modelo así se
basaría necesariamente en retroanálisis, y con frecuencia se encontraría que más de un
juego de parámetros permitiría reproducir un mismo resultado experimental.
Ante esta circunstancia, el camino más adecuado es desarrollar una ecuación
fenomenológica que reproduzca con suficiente precisión el comportamiento del material
para el fin propuesto por el modelo, aunque no capture por completo la física del
fenómeno. En este ejercicio teórico se propone la expresión
Kp =Kp
ult
χb(3.12)
56
donde b es un parámetro material. La ec. (3.12) cumple con las restriciciones de los
estados inicial y final
χ = 0 =⇒ Kp =∞ (3.13)
χ = 1 =⇒ Kp = Kpult (3.14)
3.3 Implementación numérica
3.3.1 Operador de integración
Se adopta un esquema de integración exacta para p e implícito para pc. Se debe resolver
el conjunto de ecuaciones
n+1p =
(cb1 + n+1e0
n+1e0(1−m)
(nεe
v +∆εv − n+1∆εpv
)) 1
1−m
pref (3.15)
n+1e0 = (1 + ne0) exp [−∆εpv]− 1 (3.16)
n+1pc = npc +Kp[
n+1p, n+1e0]∆εp
v (3.17)
sujeto a la restricciónn+1p− n+1pc = 0 (3.18)
3.3.2 Algoritmo de integración
La implementación numérica del modelo se basa en un algoritmo de Newton-Raphson
convencional para rigidez dependiente de la presión.
1. El algoritmo es implícito, por lo que todas las variables se evalúan en (n+ 1) excepto
que se aclare lo contrario. Se parte de un paso n convergido con variables de estado
{np, ne0,npc} (3.19)
57
2. Se inicializa
n+1e(0)0 = ne0 (3.20)
n+1p(0)c = npc (3.21)
εpv = 0 (3.22)
3. Se calcula la deformación volumétrica elástica con la ec. (3.4)
nεev =
npnK (1−m) (3.23)
4. Se asume que el incremento de deformación volumétrica total ∆εv es elástico y se
calcula
εe(0)v = nεe
v +∆εv (3.24)
5. Se calcula la presión de prueba p mediante la ec.(3.5), que es la integración exacta
de p = Kεev
p =
(cb1 + e0e0
(1−m) εe(0)v
) 1
1−m
pref (3.25)
6. Se calcula la función de fluencia
F (0) = p− p(0)c (3.26)
7. Si F (0) < 0, el paso es elástico:
(a) Se actualiza n+1p = p como una integral exacta.
(b) Las variables de estado e0 y pc mantienen sus valores convergidos ne0 y npc.
(c) El algoritmo termina.
8. Si F (0) > 0, el paso es elastoplástico:
58
(a) Se plantea el residuo del algoritmo de actualización para la iteración (i)
ζ(i) = p(i) − p(i)c (3.27)
(b) Se anula el residuo para la iteración (i+ 1) mediante una expansión de Taylor
de primer orden
ζ(i+1) = ζ(i) +dζ
dεpvδεp
v = 0 (3.28)
lo que permite calcular
δεpv = −
ζ(i)
dζ/dεpv
(3.29)
(c) Se incrementa la deformación volumétrica plástica
∆εp(i+1)v = ∆εp(i)
v + δεpv (3.30)
(d) Se actualizan εev y e0, éste último como integración exacta en el paso
εe(i+1)v = nεe
v +(∆εv −∆εp(i+1)
v
)(3.31)
e(i+1)0 = (1 + ne0) exp
[−∆εp(i+1)
v
]− 1 (3.32)
(e) Se actualizan las variables de estado. En este caso y como ejemplo, la
presión se actualiza mediante una fórmula exacta mientras que la presión
de preconsolidación se actualiza con Kp evaluado al final del paso, con las
expresiones
p(i+1) =
(cb1 + e
(i+1)0
e(i+1)0
(1−m) εe(i+1)v
) 1
1−m
pref (3.33)
p(i+1)c = npc +Kp(i+1)∆εp
v (3.34)
(f) Se calcula el nuevo residuo y se itera hasta que ζ(i+1) < tol.
59
Para que el algoritmo sea operativo se requiere la expresión de
dζ
dεpv=dp
dεpv− dpc
dεpv
(3.35)
donde
dp
dεpv=
∂p
∂εpv
∂p
∂e0
∂e0∂εp
v= −K + K
e0(nεe
v +∆εv −∆εpv) (3.36)
dpc
dεpv=
∂pc
∂εpv+∂pc
∂e0
∂e0∂εp
v= Kp +
∂pc
∂e0(1 + e0) (3.37)
En la ec. (3.37) falta ∂pc/∂e0. Si se deriva la ec. (3.17) queda
∂pc
∂e0=∂Kp
∂e0∆εp
v = −1 + b+ e0 + be0 + ρ
χbe20ρ2
pult∆εpv
por lo que la expresión final es
dζ
dεpv=
(nεe
v +∆εv −∆εpv
e0− 1)K −Kp +
1 + b+ e0 + be0 + ρ
χbe20ρ2
pult (1 + e0)∆εpv (3.38)
3.4 Resultados
3.4.1 Influencia de los parámetros de entrada
En las figuras 3-2 a 3-6 se muestra el efecto de los diferentes parámetros sobre la curva
e0 − p. Puede verse que los parámetros elásticos cb y m no tienen influencia sobre la
curva en la amplia escala de presiones estudiada. Por otro lado, la disponibilidad de tres
parámetros - pr, ρ y b - para el rango de alta presión permite ajustar cualquier curva
experimental. Puede apreciarse que los parámetros pr y ρ tienen un claro significado
físico, mientras que b es únicamente un parámetro de ajuste del modelo.
En la figura 3-7 se muestra la comparación entre predicción y resultados
experimentales para arena Toyoura, con cb = 800, m = 0.5, pr = 55, ρ = 0.4 y b = 1.
Si el objetivo de la modelización fuera la predicción del comportamiento del material
60
0.20
0.40
0.60
0.80
1.00
100 1000 10000 100000[ ]kPap
0e
Figura 3-2: Efecto del parámetro cb para tres relaciones de vacíos.
en el rango de muy altas presiones, debería investigarse la dependencia de b respecto de la
densidad relativa para un conjunto grande de arenas y gravas, y, si es posible, determinar
una función b [e0] que tuviera una capacidad predictiva acorde con ese objetivo. Si, en
cambio, el objetivo del modelo fuera la predicción del comportamiento del material en
el rango de presiones que interesa en ingeniería geotécnica - 0− 12 MPa - el parámetro bpuede tomar un valor por defecto. En la figura 3-7 se aprecia que b = 1 es una elección
razonable, al menos para los ensayos utilizados para la calibración. En la figura 3-8 se
presenta la calibración del modelo para la arena Sacramento [79]. Los parámetros son
cb = 700, m = 0.5, pr = 35, ρ = 0.4, b = 1.
3.4.2 Comportamiento del algortimo
En la figura 3-9 se presenta la curva e0−p que se obtiene con incrementos de deformación∆εv = 1%, ∆εv = 0.5% y ∆εv = 0.1%. En la figura se indica, para una de las relaciones
de vacíos, los puntos convergidos para la discretización temporal en 50 pasos.Puede
apreciarse que las curvas se superponen y que la diferencia es indistinguible en la escala
de la figura.
61
0.20
0.40
0.60
0.80
1.00
100 1000 10000 100000[ ]kPap
0e
Figura 3-3: Efecto del parámetro m para tres relaciones de vacíos.
0.20
0.40
0.60
0.80
1.00
100 1000 10000 100000[ ]kPap
0e
Figura 3-4: Efecto del parámetro pr para tres relaciones de vacíos.
62
0.20
0.40
0.60
0.80
1.00
100 1000 10000 100000
0.3ρ =0.4ρ =
0.5ρ =
[ ]kPap
0e
Figura 3-5: Efecto del parámetro ρ para tres relaciones de vacíos.
0.20
0.40
0.60
0.80
1.00
100 1000 10000 100000[ ]kPap
0e
Figura 3-6: Efecto del parámetro b para tres relaciones de vacíos.
63
0.50
0.60
0.70
0.80
0.90
1.00
100 1000 10000 100000
e = 0.83e = 0.77e = 0.59
Arena Toyoura
r refp p
[ ]kPap
e
Figura 3-7: Comparación entre simulación numérica y resultados experimentales paraarena Toyoura [94][109][110].
0.50
0.60
0.70
0.80
0.90
1.00
100 1000 10000 100000
e=0.87
e=0.78
e=0.71
e=0.61
r refp p
[ ]kPap
eArena Sacramento
Figura 3-8: Comparación entre simulación numérica y resultados experimentales paraarena Sacramento [79].
64
3.A Apéndice: Selección de e0 como variable de
estado
Si d es el tensor velocidad de deformación que corresponde a un paso finito n → n + 1,
se verifica que ∫ n+1
n
d : 1 dt = −ln[1 + n+1e
1 + ne
]= ∆εv (3.39)
donde εv es la deformación volumétrica natural, ne es la relación de vacíos al comienzo
del paso y n+1e es la relación de vacíos al final del paso [53]. Por lo tanto, el cambio de
relación de vacíos es una medida de deformación volumétrica finita que puede calcularse
mediante la expresiónn+1e = (1 + ne) exp [−∆εv]− 1 (3.40)
La tasa de deformación volumétrica correspondiente es
e = − (1 + e) εv (3.41)
que puede escribirse como
e = − (1 + e) (εpv + p/K) (3.42)
Se aprecia que e tiene una componente elástica y otra plástica y que depende de p,
que es a su vez una variable de estado. Este acoplamiento entre mecanismos elásticos
y plásticos tiene serias desventajas, porque produce modelos cuya interpretación y
calibración es más compleja que los modelos cuyas variables de estado corresponden
a mecanismos físicos desacoplados.
La componente elástica de e puede ser removida si se define una nueva variable de
estado que se mide siempre a la misma presión. En este trabajo se eligió e0, la relación
de vacíos que se mediría si se descargara la presión a cero (y esta descarga no produjera
65
deformaciones plásticas). La relación entre ambas es
1 + e01 + e
= exp [εev] (3.43)
Para un paso finito, e0 se actualiza mediante
n+1e0 = (1 +ne0) exp [−∆εp
v]− 1 (3.44)
donde ∆εpv es la deformación volumétrica plástica del paso.
66
0.20
0.40
0.60
0.80
1.00
100 1000 10000 100000[ ]kPap
0e
Figura 3-9: Comportamiento numérico del algoritmo. Las curvas e0 − p obtenidas con50, 100 y 1000 pasos son indistinguibles.
67
Capítulo 4
Modelización de la falla por corte
4.1 Introducción
Los problemas rutinarios de ingeniería geotécnica se resuelven en dos etapas: i) se verifica
la “estabilidad” de un determinado problema con su geometría, cargas y materiales; y ii)
una vez asegurada la estabilidad, se estiman las “deformaciones” que pueden producirse
por acción de esas cargas.
Para la primera verificación se emplean fórmulas de equilibrio límite desarrolladas
a partir de la hipótesis de que el material es rígido-plástico perfecto. Para la segunda
verificación se asume que el material es elástico lineal y se aplican fórmulas de la teoría
de la elasticidad. Es tarea del analista elegir los parámetros mecánicos, de resistencia
en el primer análisis y de rigidez en el segundo, que representen razonablemente el
comportamiento del material para el problema que se estudia1.
La comparación entre los cálculos de estabilidad efectuados mediante métodos
analíticos y numéricos sólo puede hacerse si en los últimos se emplean ecuaciones
constitutivas muy simples, en las que los parámetros mecánicos también se eligen en
función del problema que se estudia. El problema es que los modelos numéricos con
frecuencia resuelven varios problemas que serían tratados por separado en un análisis
1La experiencia profesional interviene de manera decisiva en esta etapa del análisis.
68
manual. Por ejemplo, el problema de una fundación ubicada cerca de un talud
es, en realidad, un problema doble: i) el comportamiento de la fundación; y ii) el
comportamiento del talud. Este problema se resuelve analíticamente con ciertas fórmulas
específicas para el problema de la fundación y con otras fórmulas específicas para el
problema del talud. En ambos casos se emplean parámetros mecánicos diferentes, aunque
el suelo es el mismo. En un modelo numérico del problema completo, un mismo juego de
parámetros, distinto de los dos juegos anteriores, sería empleado en toda la malla para
resolver ambos problemas en conjunto. No sorprende entonces que los resultados no sean
comparables y que los errores de selección de parámetros sean muy frecuentes.
En este capítulo se formula e implementa un modelo constitutivo cuyo objetivo es
el análisis de problemas en los que interesa el comportamiento del material en estado
de falla por corte. El modelo tiene los elementos necesarios para que sus parámetros
materiales sean independientes del problema que se estudia, pero no está diseñado para
la predicción de deformaciones y desplazamientos de estructuras en su estado de servicio.
Los elementos principales del modelo son: i) hiperelasticidad isotrópica, ii) ángulo de
fricción interna dependiente de la presión y relación de vacíos; iii) capacidad de simulación
del estado crítico; iv) teoría tensión-dilatancia de Rowe.
4.2 Estructura matemática del modelo
4.2.1 Cinemática
Se adopta la descomposición aditiva de las deformaciones infinitesimales
εv = εe + εp (4.1)
donde ε es el tensor de deformación total, εe es el tensor de deformación elástica y εp es
el tensor de deformación plástica.
69
4.2.2 Variables de estado
El estado del material se caracteriza mediante el tensor de tensiones efectivas σ y la
relación de vacíos para presión nula e0.
4.2.3 Relación tensión-deformación
Se define una función de energía complementaria de deformación
Ws = Ws [σ, e0] (4.2)
con la que se calcula el tensor de deformación elástica
εe = ∂Ws/∂σ (4.3)
4.2.4 Función de fluencia
La función de fluencia es de la forma
Ff [σ, e0] = 0 (4.4)
4.2.5 Regla de flujo
Se adopta asociatividad deviatórica y no asociatividad volumétrica. El incremento de
deformación plástica se calcula con
εp = λm (4.5)
donde λ es un multiplicador plástico ym es el tensor de dirección de deformación plástica.
Su expresión es
m =md + β1 (4.6)
70
donde β es un factor de dilatancia,
md = nd/ ‖nd‖ (4.7)
es un tensor unitario deviatórico,
nd = Idev : n (4.8)
y
n = ∂Ff/∂σ (4.9)
es la normal (exterior) a la función de fluencia.
4.2.6 Condición de plasticidad
Se emplea la condición
λ =
λ si: Ff = 0 y: Ff > 0
0 en otro caso(4.10)
4.2.7 Ecuaciones de evolución
La evolución de e0 está dada por la expresión
e0 = − (1 + e0) εpv (4.11)
4.3 Funciones de estado
4.3.1 Elasticidad
Hipoelasticidad vs. hiperelasticidad
Las arenas - conjuntos de partículas en contacto - no exhiben un comportamiento
puramente elástico para ninguna trayectoria ni escala de deformación. Por lo tanto, la
71
incorporación de elementos de elasticidad para arenas se justifica principalmente porque
se los necesita para la integración de las ecuaciones constitutivas.
Existen dos marcos teóricos fundamentales: hipoelasticidad e hiperelasticidad. En el
primero se postula la relación incremental
σ = D [σ, e0] : εe (4.12)
mientras que en el segundo se postula la relación integral
σ = D [σ, e0] : εe (4.13)
En ambos, el tensor elástico de cuarto orden D [σ, e0] es función de las variables de
estado del modelo. Para que la ec. (4.13) pueda ser formulada, debe existir un potencial
elástico escalar Ws [σ, e0] tal que
εe =
∂Ws [σ, e0]
∂σ: σ (4.14)
A esta función potencial se la denomina “función de energía complementaria
de deformación”. Existen muchas funciones de este tipo para geomateriales (p.ej.
[95][103][116]), que tienen en cuenta la dependencia de los módulos elásticos respecto del
estado tensional. Algunas de estas funciones permiten modelar anisotropía y distintos
grados de acoplamiento entre los módulos volumétricos K y de corte G.
Los modelos hipoelásticos no requieren la formulación de potenciales elásticos, y por
lo tanto permiten el empleo de funciones simples para K y G. Su desventaja radica
en que no satisfacen de manera automática el requerimiento de que un ciclo cerrado de
deformación produzca un comportamiento conservativo. En cambio, las formulaciones
integrales de la forma de la ec. (4.13) aseguran un comportamiento conservativo para
cualquier ciclo cerrado de deformación.
72
Si Ws produce un tensor elástico de la forma
D [σ, e0] = K [σ, e0]11 + 2G [σ, e0] Idev (4.15)
se puede escribir
σ = K [σ, e0] (εv − εpv)1 + 2G [σ, e0] (εd − εp
d) (4.16)
o, alternativamente,
σ = σ − (K [σ, e0] εpv1 + 2G [σ, e0] ε
pd) (4.17)
donde
σ = K [σ, e0] εv1 + 2G [σ, e0] εd (4.18)
es la “tensión de prueba”.
Como Ws es una función escalar, el tensor εe hereda los autovectores de σ a través
de la ec. (4.14). La ec. (4.18) implica que σ tiene los mismos autovectores que εe y, por
lo tanto, que σ. La ec. (4.17), entonces, se puede escribir como
D [σ, e0] : εpd = σ − σ (4.19)
e implica que εp comparte autovectores con σ, σ y ε. La consecuencia inmediata de esta
deducción es que los autovectores de la solución pueden ser calculados con información
del paso de prueba. Esta característica constituye una importante ventaja para la
formulación algorítmica de los modelos, que sólo requiere la existencia de las formas
funcionales dadas por las ecs. (4.14) y (4.15).
Función de energía complementaria de deformación
Se adopta la función de energía complementaria de deformación de Molenkamp [95]
Ws =p2
G
[S
(1−m) (2−m) +r2
4
](4.20)
73
donde G es el módulo de rigidez al corte para el que se propone la expresión de Hardin
[49][50][118]
G = cs(ce − e0)21 + e0
(p
pref
)m
pref (4.21)
S es una función auxiliar a definir, cs, ce ym son parámetros materiales, y pref = 100 kPa
es una presión de referencia.
Deformación elástica
La derivación de la ec. (4.20) respecto a σ permite calcular el tensor de deformación
elástica
εe =
(1
2GIdev+
1
K11
): σ (4.22)
donde K el el módulo de rigidez volumétrica elástica, relacionado con G mediante
K =G
S/ (1−m)−mr2/4 (4.23)
La inversión de la ec. (4.22) permite calcular
σ = D : εe (4.24)
donde
D = K11+ 2GIdev (4.25)
es el tensor elástico.
Debe notarse que la ec. (4.23) arroja K →∞ para
rlim →√4S/ (m (1−m)) (4.26)
Esta oblicuidad límite, sin embargo, no tiene influencia en los casos prácticos porque
74
es mucho más alta que la que cualquier material friccional puede soportar2.
La relación de Poisson resultante de las ecs. (4.21) y (4.23) es
ν =3− 4S + r2/44S + 6− r2/4 (4.27)
que resulta dependiente de r. Para estado isotrópico, (r = 0) queda
ν =3− 4S6 + 4S
(4.28)
El parámetro S que la determina no es un dato de entrada sino la función de estado
S =3−Dr
4(4.29)
donde
Dr =emax − e0emax − emin
(4.30)
y emax y emin son parámetros del modelo. Esta expresión arroja
Dr = 0%→ S =3
4→ ν = 0 (4.31)
Dr = 100%→ S =1
2→ ν = 0.125 (4.32)
que son razonablemente coincidentes con los valores determinados en [127]. Debe notarse
que, debido a la dependencia de ν respecto de r, la relación de Poisson puede alcanzar
ν = 0.21 en el caso de arenas densas en falla.
2Se requeriría φ > 73◦ para que rlim correspondiera a un estado admisible.
75
4.3.2 Ángulo de fricción interna
Definición
El ángulo de fricción interna φ se define como una función de estado de la forma
φ = φc + ψ (4.33)
donde φc es el ángulo de fricción interna crítico y ψ [p, e0] es una función de estado que
tiene en cuenta la dependencia que φ tiene respecto de la dilatancia. Es importante
destacar que la función de estado φ se define únicamente para el estado de tensiones
caracterizado por σ1 > σ2 = σ3, correspondiente a la compresión triaxial.
Dilatancia y ángulo de fricción interna
La expresión de Bolton para el ángulo de fricción interna (ec. (2.15)) [11] es
φ = φc +∆φ Dr (Q− ln [p/pref ])− b (4.34)
donde ∆φ = 3◦, b = 3◦ para la mayoría de las arenas [11] y Q es un parámetro material
que tiene en cuenta la resistencia de los granos. El significado de Q se enfatiza si las ecs.
(4.33) y (4.34) se combinan en
ψ = −∆φ Drln [χB]− b (4.35)
donde
χB =p
exp [Q] pref
(4.36)
es una medida del nivel de tensiones.
En el modelo de compresión isotrópica se mostró que la presión que produce un
determinado nivel de ruptura de partículas depende tanto de la presión como de la
relación de vacíos y se introdujo la ec. (3.8) que define la función de estado χ. Por lo
76
tanto, en la ec. (4.35) conviene reemplazar a χB por χ con lo que queda
ψ = −∆φDrln [χ]− b (4.37)
En la definición de χ hay dos parámetros materiales, pr y ρ. Pestana [109] mostró
que ρ cae en el rango aproximado 0.33 < ρc < 0.45 para la mayoría de las arenas. Como
ρ tiene poca influencia en el comportamiento de las arenas en el rango de presiones de
interés ingenieril, en este modelo se adopta un valor constante ρc = 0.40. La calibración
de la ec. (4.37) con los datos usados por Bolton para desarrollar la ec. (4.34) arroja
∆φ = 3◦ y b = 2◦, que son parámetros por defecto de este modelo, por lo que los
parámetros materiales que intervienen en la ec. (4.37) son emin, emax, φc y pr.
La figura 4-1 muestra ψ para la arena del Río Sacramento [79] y para la arena Toyoura
[12][44], junto con las predicciones de la ec. (4.37) [137][138][139][140].
0
2
4
6
8
10
12
14
16
100 1000 10000
SS: e = 0.87SS: e = 0.61TS: e = 0.86TS: e = 0.78TS: e = 0.70TS: e = 0.66
ψ
[ ]kPap
SS: 0.61 1.03 32 35
TS: 0.61 0.98 33 55min max c r
min max c r
e e p
e e p
φφ
= = = ° == = = ° =
Figura 4-1: ψ = φ − φc en función de p y Dr para las arenas Sacramento y Toyoura,junto con la predicción de la ec. (4.37). Datos de [12][79].
Bolton [11][12] limitó la validez de la ec. (4.34) a p > 150 kPa para evitar la sobre-
estimación de la dilatancia a muy bajas presiones. En principio, la misma limitación se
aplica a la ec. (4.37).
77
Estado crítico
La función de estado ψ [p, e0] debe cumplir con las restricciones
ψ [p, e0c] = 0 (4.38)
∂ψ/∂e0 < 0 (4.39)
∂ψ/∂p < 0 (4.40)
donde e0c es la relación de vacíos para presión nula en el estado crítico. Tanto la ec.
(4.35) como la ec. (4.37) pueden usarse para predecir la relación e0c − p si se resuelve laecuación implícita
ψ [p, e0c] = 0 (4.41)
La ec. (4.35) arroja
p⌋e0c= exp
[Q− b
∆φ Dr [e0c]
]pref (4.42)
mientras que la ec. (4.37) arroja
p⌋e0c=pr
e2.50cexp
[ −b∆φ Dr [e0c]
]pref (4.43)
La figura 4-2 muestra determinaciones de pares de valores e0c−p para la arena Toyoura[154], junto con la predicción dada por las ecs. (4.42) y (4.43) con sus parámetros
por defecto. Puede observarse que la ec. (4.42) erra totalmente el comportamiento
experimental mientras que la ec. (4.43) lo ajusta razonablemente, a pesar de que en su
deducción no se tuvo en consideración la información experimental contenida en la figura
4-2. Este mejor ajuste prueba que la función de estado χ es una mejor medida del efecto
de las variables de estado {p, e0} sobre la ruptura de granos que el parámetro χB. La
línea de estado crítico de una dada arena puede ser reproducida con mayor exactitud si
se modifican los parámetros por defecto {∆φ, b} [137][138][139][140].
78
0.70
0.75
0.80
0.85
0.90
0.95
10 100 1000 10000
Experimental
[ ]kPap
0ce
Bolton
Este trabajo
Figura 4-2: Líneas de estado crítico predichas por las ecs. (4.42) y (4.43).
4.3.3 Función de fluencia
Se adopta la función de fluencia de Matsuoka-Nakai [90] que, escrita en términos de r,
tiene la expresión
Ff =(µf + 6
)J2r −
(µf + 9
)J3r − µf = 0 (4.44)
donde J2r = 12r : r y J3r = 1
3r · r : r y
µf = 8 tan2 [φ] (4.45)
es una función de estado que determina la apertura de la superficie de fluencia.
La figura 4-3 muestra una vista 3D de la superficie de fluencia definida por la ec.
(4.44). La curvatura de la directriz está controlada por la dependencia de Ff respecto de
J3r, mientras que la curvatura en la dirección de la generatriz está dada por la dependencia
de φ respecto de la presión.
µf evoluciona de acuerdo a
µf =∂µf
∂φ
∂φ
∂ψ
(∂ψ
∂e0e0 +
∂ψ
∂pp
)(4.46)
79
Figura 4-3: Superficie de fluencia, ec. (4.44) [90][91].
En una trayectoria a presión constante, la ec. (4.46) toma la forma
µf =∂µ
∂φ
∂φ
∂ψ
∂ψ
∂e0e0 (4.47)
Las derivadas parciales de la ec. (4.47) cumplen con las condiciones
∂µf/∂φ > 0 (4.48)
∂φ/∂ψ > 0 (4.49)
∂ψ/∂e0 < 0 (4.50)
Por lo tanto, el signo de µf depende del signo de e0. Pueden distinguirse dos
situaciones: i) una muestra suelta reduce su relación de vacíos, por lo que e0 < 0 y
por lo tanto µf > 0, o sea, el cono se expande y la muestra endurece; y ii) una muestra
densa dilata, por lo que e0 > 0 y por lo tanto µf < 0, o sea, el cono se contrae y la
muestra ablanda.
La solución completa de la ec. (4.44) se presenta en la figura 4-4 para tres valores
de µf . Puede observarse que Ff = 0 genera una familia de superficies, de las cuales
80
sólo una está orientada en el eje p > 0. Esta característica matemática de la ec. (4.44)
es incompatible con los algoritmos de retorno convencionales, por lo que se requiere la
formulación de algorítmos específicos, de los cuales uno se presenta más adelante.
Figura 4-4: Todas las soluciones de la ec. (4.44) evaluadas para µf = 3, 6 y 9.
4.3.4 Regla de flujo
Introducción
La forma general de la teoría de tensión-dilatancia de Rowe está dada por la ec. (2.46)
que se repite a continuación [125]
win = Ncvwout (4.51)
Rowe calculó los términos win y wout como la suma de las componentes positivas y
negativas del triplete
wI = σI εI (I no suma) (4.52)
81
o sea,
win =∑
wI>0
wI (4.53)
wout = −∑
wI<0
wI (4.54)
Dentro de la teoría de la plasticidad, el trabajo de disipación está dado por el producto
σ : εp. Rowe [125] justificó el empleo de ε en lugar de εp en el hecho de que la energía
elástica almacenada es muy pequeña comparada con la energía disipada en los procesos
de deformación por corte de arenas. En lo que sigue, se asume que en la ec. (4.52) se
emplean componentes de εp.
El cálculo del trabajo de deformación como el producto escalar de las componentes
de σ y εp sólo tiene sentido físico si ambos tensores comparten autovectores y si los
autovalores tienen el mismo ordenamiento. En cualquier otro caso, el trabajo es de la
forma σ : εp y no puede ser reducido a la suma de los elementos de wI .
Para que σ y εp compartan autovectores es necesario que el potencial plástico sea
una función isotrópica de σ. Aunque el problema no es trivial para el caso de la no
asociatividad general, puede resolverse si se postula asociatividad deviatórica (en εpd) y
no asociatividad volumétrica (en εpv). Como la proyección volumétrica no aumenta el
espacio espectral, la no asociatividad volumétrica no rompe la herencia de propiedades
espectrales de σ y εp. En el caso de no asociatividad volumétrica, la disipación plástica
unitaria puede escribirse como
wp = σ1εp1 + σ2ε
p3 + σ3ε
p3 (4.55)
y la ec. (4.52) puede ser redefinida como
wI = σI εpI (I no suma) (4.56)
que ahora contiene las componentes del trabajo plástico disipado.
82
Cálculo del término de dilatancia
Si se sustituye la ec. (4.6) en la ec. (4.5) y luego en la ec. (4.56) queda
wp = λ {σ1 (md1 + β) , σ2 (md2 + β) , σ3 (md3 + β)} (4.57)
donde {md1, md2, md3} son los autovalores de md y β es un término de dilatancia
dependiente de σ.
El deslizamiento de partículas, como mecanismo plástico, requiere que no todos los
componentes del triplete {εp1, ε
p2, ε
p3} tengan el mismo signo. Como los autovalores de σ
son positivos, debe cumplirse que σ1 (md1 + β) > 0 (aporta a win) y σ3 (md3 + β) < 0
(aporta a wout), mientras que el signo de σ2 (md2 + β) es desconocido.
Se asume primero que (md2 + β) ≤ 0. Si se inserta la ec. (4.56) en las ecs. (4.53) y
(4.54) y luego en la ec. (4.51) se obtiene, luego de algo de álgebra
β = −σ1md1 +Ncv (σ2md2 + σ3md3)
σ1 +Ncv (σ2 + σ3)(4.58)
mientras que, si (md2 + β) ≥ 0 se obtiene
β = −σ1md1 + σ2md2 +Ncvσ3md3
σ1 + σ2 +Ncvσ3(4.59)
Si (md2 + β) = 0, tanto la ec. (4.58) como la ec. (4.59) arrojan
β = −σ1md1 +Ncvσ3md3
σ1 +Ncvσ3(4.60)
La ec. (4.58) puede asociarse a estados cercanos a compresión triaxial; la ec. (4.59)
a estados cercanos a extensión triaxial; y la ec. (4.60) a deformación plana. Como las
tres obedecen a la teoría de la tensión-dilatancia de Rowe [124], un mismo cociente σ1/σ3
implica diferentes valores de β para los tres estados.
83
Compresión y extensión triaxial
Los autovalores de md para compresión triaxial son md1 =√2/3, md2 = md3 = −
√1/6,
y por lo tanto
εp1 = λ
(√2/3 + βTC
)(4.61)
εp3 = λ
(−√1/6 + βTC
)(4.62)
donde βTC es el término de dilatancia para compresión triaxial. Si se insertan las ecs.
(4.61) y (4.62) en la ec. (2.47) queda
βTC =√2/3
Ncvσ3 − σ1σ1 + 2Ncvσ3
(4.63)
Para extensión triaxial, el término de dilatancia queda
βTE =√2/3
Ncvσ3 − σ12σ1 +Ncvσ3
(4.64)
La anterior no es la primera extensión de la teoría de tensión-dilatanica a estados
tridimensionales de tensiones. Guo y Stolle [47] presentaron una extensión de la hipótesis
de la ec. (4.51) y de la ecuación tensión-dilatancia por medio del tensor auxiliar
πij = σikεkj (4.65)
cuyos autovalores coinciden, para asociatividad deviatórica, con
πI = σI εI (I no suma) (4.66)
La mayor ventaja de esta solución, comparada con las ecs. (4.58), (4.59) y (4.60)
es que el término de dilatancia que se obtiene es una forma tensorial derivable. Dos
desventajas son: i) el procedimiento de generalización de la ec. (4.51) se aparta de
84
la formulación de Rowe; y ii) se incluye la componente elástica de la deformación
incremental, en contradicción con el hecho aceptado de que la dilatancia es un mecanismo
plástico. Guo y Wan [48] investigaron arreglos granulares 2D con diferentes disposiciones
de partículas regulares y demostraron que la hipótesis de mínimo cociente de trabajos
implica una cota inferior de la dilatancia de materiales anisótropos. Desarrollaron una
nueva formulación tensión-dilatancia basada en micromecánica 2D y dedujeron una regla
de flujo consistente con esta formulación. Comparada con las ecs. (4.58), (4.59) y (4.60),
la regla de flujo de Guo y Wan tiene mejor sustento en la microescala, pero requiere
la definición y calibración de un tensor de fábrica. Además, no puede ser extendida a
estados tridimensionales de tensiones de manera sencilla. Por el contrario, las ecs. (4.58),
(4.59) y (4.60) sólo se basan en la hipótesis de la asociatividad deviatórica y no requieren
parámetros adicionales, por lo que pueden ser implementadas en códigos numéricos de
manera expeditiva.
Término de dilatancia interpolado
El conjunto de ecuaciones (4.58), (4.59) y (4.60) forman una función continua de σ con
respecto al al ángulo de Lode θ. Sin embargo, se produce una discontinuidad en la primera
derivada ∂β/∂θ en el punto común de las ecs. (4.58) y (4.59), donde vale la ec. (4.60).
Este pico es producido por el rápido cambio de dirección de md que existe cerca del
punto que representa el estado de compresión triaxial. Es más, la teoría de Rowe no
podría utilizarse junto con el criterio de Mohr-Coulomb - para el que ha sido deducida
- sin enfrentar el severo inconveniente de que la superficie de fluencia facetada de Mohr-
Coulomb produce un cambio finito en los autovalores demd para un cambio infinitesimal
de θ en los seis vértices de compresión/extensión triaxial.
Excepto en la vecindad del punto donde md2 + β = 0, β es una función suave del
ángulo de Lode que es numéricamente indistinguible de la interpolación lineal
β = βTC
1− sin [3θ]2
+ βTE
1 + sin [3θ]
2(4.67)
85
Tanto para compresión como extensión triaxial, el criterio de Matsuoka-Nakai arroja
σ1/σ3 = Ntx, donde
Ntx =1
4
(4 + µf +
õf
√8 + µf
)(4.68)
Si se combinan las ecs. (4.63), (4.64), (4.67) y (4.68) queda
β =Ncv −Ntx√
6
(1− sin [3θ]Ntx + 2Ncv
+1 + sin [3θ]
2Ntx +Ncv
)(4.69)
donde sin [3θ] = −3√6J3r/r
3. Woodward propuso una expresión similar [162].
En la figura 4-5 se muestra β como función de θ para tres valores µc = 8 tan2 [φc] ,
32µc y
12µc. El tercio superior muestra los valores calculados con las ecs. (4.58), (4.59),
y (4.60), mientras que los dos tercios inferiores muestran valores calculados mediante la
ec. (4.69).
Figura 4-5: Término de dilatancia β. Tercio superior: ecs. (4.58) a (4.60). Dos terciosinferiores, ec. (4.69).
86
Ángulo de fricción interna a volumen constante
La función de estado
Ncv = tan2 [π/4 + φcv/2] (4.70)
depende de la relación de vacíos, porque φcv depende de la relación de vacíos como se
muestra en la figura 2-6 [123][129]. Como φcv → φc a medida que la muestra se acerca al
estado crítico (figura 2-5), φcv es una función de estado que puede ponerse en la forma
φcv = φc − f (4.71)
f [p, e0] debe cumplir con las restricciones
f [p, e0c] = 0 (4.72)
∂f/∂p < 0 (4.73)
∂f/∂e0 < 0 (4.74)
para que se reproduzca el comportamiento observado. Si se comparan la definición y
propiedades de ψ [p, e0] y f [p, e0] se aprecia que f puede ser tomada como proporcional
a ψ. Por lo tanto, se propone la expresión
φcv = φc − αψ (4.75)
donde α es un parámetro de calibración que en este modelo se adopta en α = 0.5. La
principal ventaja de esta formulación es que asegura que tanto φcv como φ tienden a φc a
medida que el material se acerca a e0c, sin necesidad de recurrir a parámetros materiales
adicionales.
87
4.4 Implementación numérica
4.4.1 Operador de integración
El operador de integración es implícito. Se resuelve el conjunto de ecuaciones
n+1σ = D
[n+1σ, n+1e0
]: (nε+∆ε−∆εp) (4.76)
n+1e0 = (1 + ne0) exp[
n+1∆εp : 1]− 1 (4.77)
sujeto a la restricción
Ff
[n+1σ, n+1µf
]= 0 (4.78)
donden+1µf = µf
[n+1σ, n+1e0
](4.79)
Las ecs. (4.76) a (4.78) están expresadas en términos del tensor de deformación
plástica
∆εp = λ(md
[n+1σ, n+1e0
]+ β[
n+1σ, n+1e0
]1)
(4.80)
que tiene una única variable independiente λ, por lo que, en un paso elastoplástico, el
conjunto de ecuaciones (4.76) a (4.78) se reduce a
n+1σ = σ [λ] (4.81)
n+1e0 = e0 [λ] (4.82)
con lo que la ec. (4.80) se reduce a
Ff [λ] = 0 (4.83)
Este operador genera un algoritmo encestado que resuelve una variable escalar a
la vez. El algoritmo pertenece a una familia de métodos de integración multinivel
88
[43][81][83][159][160], diferente a la mayoría de los algoritmos de actualización de tensiones
para geomateriales, en los que se itera sobre un conjunto de variables simultáneas
[1][13][14][15][16][37]. En [43] puede verse una discusión sobre los algoritmos de nivel
simple y múltiple y el planteo de algoritmos encestados, mientras que en [54] se
presenta un estudio completo de las técnicas modernas de integración de modelos para
geomateriales.
4.4.2 Algoritmo de integración
Se presenta el algoritmo de integración del modelo entre el paso n y el paso n + 1.
El algoritmo es implícito, por lo que las variables se evalúan en n+1 (∗) , aunque este
superíndice se omite por claridad. Las iteraciones se indican como un superíndice derecho
entre paréntesis (∗)(i) .
El algoritmo de integración está organizado como una serie de bucles encestados con
la siguiente estructura:
1. Se inicializan las funciones de estado
(a) Se parte de un estado convergido
{ nσ, ne0} (4.84)
(b) Se impone un incremento de deformación ∆ε
(c) Se calculan las funciones de estado
µ(0)f = nµf (4.85)
β(0) = nβ (4.86)
89
(d) Se calculan las funciones elásticas
K(0) = nK (4.87)
G(0) = nG (4.88)
D(0) = nD (4.89)
(e) Se calcula la deformación elástica convergida
nε
e = nD−1 : nσ (4.90)
2. Se calcula la deformación elástica con congelamiento de variables plásticas
ε = nε
e +∆ε (4.91)
3. Se cálcula la tensión elástica σe con congelamiento de variables plásticas (Algoritmo
A)
4. Se controla la función de fluencia
F(0)f
[σ
e, nµf
]≤ 0 =⇒
σ = σe
e0 =ne0
y FIN (4.92)
5. Si F (0)f > 0
(a) Se calcula la tensión de prueba
σ(0) = D(0) : ε (4.93)
(b) Se estima una cota inferior del multiplicador λ(0) y de εp(0) (Algoritmo B)
(c) Se calculan los autovectores unitarios gi de s(0) = Idev : σ(0) (gi tienen módulo
90
1 y traza 0)
(d) Mientras
abs[F (i)]< tol (4.94)
i. Se actualizan
σ(i) λ(i) ε
p(i) (4.95)
con
D(i) e(i)0 µ
(i)f β(i) (4.96)
congelados (Algoritmo C)
ii. Se actualiza
e(i+1)0 K(i+1) G(i+1) D(i+1) (4.97)
iii. Se actualiza
σ(i+1) = D(i+1) : ε (4.98)
iv. Se actualiza
σ(i+1) = D(i+1) :
(ε− εp(i+1)
)(4.99)
v. Se actualizan las funciones de estado
µ(i+1)f β(i+1) (4.100)
vi. Se calcula la función de fluencia
F(i+1)f
[σ(i+1), µ
(i+1)f
]⋚ 0 (4.101)
vii. Se incrementa i e itera.
(e) Se informa σ, e0 y FIN.
91
Algoritmo A: Cálculo de σe
Se debe encontrar σe = s+ p1 tal que se cumpla la igualdad
εe = D : (s+ p1) (4.102)
con εe fijo. Como los términos volumétricos y deviatóricos están desacoplados, la ec.
(4.102) se reduce a
εed =
s
2G(4.103)
εev =
p
K(4.104)
Para eso se plantea el residuo
ζ = p−K [p, r [p]] εev (4.105)
y se lo anula para la iteración (i+ 1) (i es local del algoritmo)
ζ(i+1) = ζ(i) +
(1− εe
v
(dK
dp
)(i))δp = 0 (4.106)
con lo que se obtiene el incremento de presión
δp(i+1) =−ζ(i)
1− εev (dK/dp)
(i)= − p−Kεe
v
1− εev
mKp
(1 + (m− 1) Kr2
2G
) (4.107)
Debido a las fuertes no linealidades presentes en la formulación, el algoritmo debe ser
completado con un esquema de cotas de error como sigue:
1. Se elige una cota inferior de la presión p(0)inf y una cota superior p(0)sup suficientemente
amplias.
92
2. Se calcula
G(0) r(0) = 2G(0)εed/p
(0) K(0) (4.108)
3. Se calcula el residuo
ζ(0) = p(0) −K(0)εev (4.109)
4. Se corrigen cotas
ζ(0) < 0 =⇒ p(0)inf = max
[p(0)inf , p
(0)]
(4.110)
ζ(0) > 0 =⇒ p(0)sup = min[p(0)sup, p
(0)]
(4.111)
5. Mientras abs[ζ(i)]> tol
(a) Se calcula
p(i+1) = p(i) − ζ(i)/(1− εe
v (dK/dp)(i))
(4.112)
(b) Se controlan cotas
p(i+1) > p(i)sup ∨ p(i+1) < p(i)inf =⇒ p(i+1) =
p(i)sup + p
(i)inf
2(4.113)
(c) Se actualizan
G(i+1) r(i+1) K(i+1) (4.114)
(d) Se actualiza
ζ(i+1) = p(i+1) −K(i+1)εev (4.115)
(e) Se corrigen cotas
ζ(i+1) < 0 =⇒ p(i+1)inf = max
[p(i)inf , p
(i+1)]
(4.116)
ζ(i+1) > 0 =⇒ p(i+1)sup = min[p(i)sup, p
(i+1)]
(4.117)
93
(f) se incrementa i y se itera.
Algoritmo B: Estimación de εp
La figura 4-4 muestra que existen varios dominios en el espacio {σ} en los que Ff [σ] < 0,
de los cuales sólo el que tiene su eje coincidente con el eje p > 0 kPa tiene sentido físico.
Esta solución queda circunscripta por la traza del cono de sección circular
FDP = ‖s‖ − ρDP p (4.118)
donde
ρDP =2√6µf
3√µf + 8−
õf
(4.119)
Se puede calcular una primera estimación de λ mediante un retorno radial a la ec.
(4.118). Para ello se calculan
s(0) = s− 2G s
‖s‖λ(0) (4.120)
p(0) = p− 3Kβλ(0) (4.121)
se insertan en la ec. (4.118) y se despeja
λ(0) =‖s‖ − ρDP p
2G− 3KρDPβ(4.122)
Esta expresión permite un cálculo de λ(0) aún en el caso en que p < 0 kPa. La
deformación plástica se estima con
εp(0) = λ(0)
(s
‖s‖ + β1)
(4.123)
94
Algoritmo C: Actualización de σ con funciones de estado constantes
En este algoritmo se actualiza σ, λ y εp mientras se mantienen congelados e0, los módulos
elásticos K, G y, por lo tanto, el tensor elástico D y σ. El algoritmo recibe σ, σ, e0, K,
G, D, µf , β. Se utiliza el siguiente procedimiento (i es un contador local):
1. Se calcula el versor de deformación plástica deviatórica m(0)d (Algoritmo D)
2. Se corrige λ(0) y εp(0)mediante una proyección a la nueva dirección plástica:
λ(i) = εp(0) :m
(0)d (4.124)
εp(0) = λ(i)m
(0)d (4.125)
3. Se calcula
σ(0) = σ −D : λ(0)m
(i)d (4.126)
4. Mientras
abs
[F(0)f
]> tol (4.127)
(a) Se calcula el incremento de multiplicador
∆λ = F(i)f /(n(i) : D :m(i)
)(4.128)
(b) Se actualiza el multiplicador
λ(i+1) = λ(i) +∆λ (4.129)
(c) Se actualiza
σ(i+1) = σ −D : λ(i+1)m(i) (4.130)
(d) Se actualiza
n(i+1) m(i+1) (4.131)
95
(e) Se evalúa la función de fluencia
F(i+1)f ⋚ 0 (4.132)
Para que el código sea operativo se requiere el valor de n = ∂Ff/∂σ que es
n =1
p
((r2 +
(µf + 9
)J3r)1 +(µf + 6
)r−(µf + 9
)r · r)
(4.133)
Algoritmo D: Cálculo de md
En esta sección se propone un algoritmo que permite la determinación del tensor de
dirección plástica deviatórica md para el criterio de Matsuoka-Nakai. Como este criterio
tiene más de una superficie de fluencia asociada (ver figura 4-4), el campo tensorial
n = ∂F/∂σ no apunta siempre a la superficie objetivo, lo que constituye un severo
inconveniente para la mayoría de los algoritmos de retorno convencionales.
Para solucionar este inconveniente, el algoritmo que se presenta emplea únicamente
información contenida en la superficie de fluencia y el valor de la tensión de prueba σ.
El procedimiento es:
1. Se calculan los seis vértices de la traza de Ff
v1 =
√2
3ρout g1 v2 = −
√2
3ρin g3 v3 =
√2
3ρout g2
v4 = −√2
3ρin g1 v5 =
√2
3ρout g3 v6 = −
√2
3ρin g2
donde gi son los autovectores deviatóricos de σ y
ρout =2√6µf
3√µf + 8−
õf
ρin =2√6µf
3õf + 8 +
õf
son los valores de r en los vértices exteriores e interiores de Ff .
96
2. Se calcula
r = s/p (4.134)
3. Se calcula la distancia entre r y los seis vértices
di = r : vi (4.135)
y se eligen los dos vértices más cercanos vI y vII .
4. Se calcula el versor vAB de la recta que une vI y vII .
5. Se calcula el ángulo entre el vértice más cercano vI y r como
α = acos [vI : r/ ‖r‖] (4.136)
6. Como el ángulo entre dos vértices es 60◦, la coordenada local del punto donde r
intersecta a vII − vI , medida desde vI , es
0 ≤ x = 3α/π ≤ 1/2 (4.137)
como se muestra en la figura 4-6
7. Se eligen dos coordenadas cercanas
{xA =
x
2, xB =
x
2+1
4
}(4.138)
97
Figura 4-6: Determinación de coordenada local x0.
8. Se calcula r y md en las dos coordenadas
rA = r [xA] = vI + xA (vII − vI) (4.139)
rB = r [xB] = vI + xB (vII − vI) (4.140)
mAd = md [rA] (4.141)
mBd = md [rB] (4.142)
9. Se calcula el error de paralelismo entre (r− ri) , mAd y mB
d
ζA =((r− rA)−
((r− rA) :m
Ad
)mA
d
):vII − vI
‖vII − vI‖(4.143)
ζB =((r− rB)−
((r− rB) :m
Bd
)mB
d
):vII − vI
‖vII − vI‖(4.144)
10. Mientras
ζ > tol (4.145)
98
(a) Se calcula una nueva coordenada x por interpolación lineal
x =xA ζB − xB ζA
ζB − ζA
como se muestra en la figura 4-7
Figura 4-7: Determinación del error de paralelismo de las normales mAd y mB
d y cálculode mx
d.
(b) Se calculan
rx = r [x] = vI + x (vII − vI) (4.146)
mxd = mx [rx] (4.147)
(c) Se calcula el error de paralelismo
ζ = ((r− rx)− ((r− rx) :mxd)m
xd) :
vII − vI
‖vII − vI‖(4.148)
99
(d) Se actualizan los umbrales
x < xA → xB = xA ∧ xA = x (4.149)
x > xB → xA = xB ∧ xB = x (4.150)
xA < x < xB →sign [ζA ζ] > 0→ xA = x
sign [ζA ζ] < 0→ xB = x(4.151)
(e) Se itera.
Con rx se conoce el punto sobre Ff de retorno de σ, por lo que se puede calcular
md [rx] de manera exacta.
4.5 Resultados
4.5.1 Introducción
En esta sección se presentan algunas simulaciones en las que se resaltan los aspectos del
comportamiento de las arenas que dependen de la evolución de la relación de vacíos, y
que por lo tanto pueden ser capturados con el modelo que se presenta en este capítulo. La
evaluación del efecto de los parámetros de entrada y las variables de estado y el análisis
numérico del algoritmo se tratan en el capítulo 5, en el que este modelo se complementa
con endurecimiento isotrópico previo a la falla.
4.5.2 Simulación de ensayos triaxiales
Se ejecutaron simulaciones numéricas de ensayos convencionales de compresión triaxial,
con la malla axilsimétrica que se muestra en la figura 4-8 (66 elementos de seis nodos,
161 nodos en total). La tecnología del elemento se describe en [18]. Los desplazamientos
verticales y horizontales del borde inferior están restringidos; en el borde lateral existe
un contorno de presión constante σ3; el borde superior tiene desplazamientos verticales
100
impuestos monotónicamente crecientes y desplazamientos horizontales nulos. Estas
condiciones de borde son las típicas de un ensayo triaxial convencional.
Figura 4-8: Malla de elementos finitos axilsimétrica para la simulación de ensayos decompresión triaxial. 66 elementos de seis nodos (detalles en [18]).
Se adoptaron los siguientes parámetros materiales: emin = 0.58, emax = 0.98, cs = 840,
ce = 2.17, m = 0.5, φc = 31◦ y pr = 55.
En la figura 4-9 se presentan las curvas σ1/σ3−ε1 y εv−ε1 de una simulación efectuadapara una densidad relativa Dr = 100% y una presión de confinamiento σ3 = 100 kPa para
siete puntos de la malla. Puede apreciarse que los diferentes puntos tienen la misma curva
σ− ε hasta cerca del pico, en el que se producen las primeras deformaciones plásticas. A
partir de ese paso, los puntos cercanos a los bordes se descargan mientras que los puntos
centrales continúan su carga plástica con ablandamiento producido por el aumento de
e0. También se muestra la respuesta global que se obtendría con mediciones externas
a la muestra. Este comportamiento, en el que diferentes puntos de la probeta tienen
trayectorias de tensiones ampliamente diferentes, es típico de los materiales densos en los
que la falla es localizada [151][152][153][157].
En las figuras 4-10 y 4-11 se presentan las curvas σ1/σ3 − ε1 y εv − ε1 para el punto
101
1
2
3
4
5
6
0.00 0.02 0.04 0.06 0.08
3 100 kPaσ =
1ε
1
3
σσ C
D
EB
F
ABCDEFG
AG 3
P H
A H
δσ
−
-0.08
-0.06
-0.04
-0.02
0.00
0.02
0.00 0.02 0.04 0.06 0.08
3 100 kPaσ =
1ε
vε
0.58 0.98 840 55 31min max s r ce e c p φ= = = = = °
C D E
Figura 4-9: Curvas σ1/σ3 − ε1 y εv − ε1 para siete puntos de la malla indicada en la fig.(4-8) y respuesta global. Ensayo triaxial.
D central de la malla, para cuatro densidades relativas y dos presiones de confinamiento.
Puede observarse que, excepto por la parte inicial de la curva σ − ε que en el modelo
es muy rígida, la respuesta global reproduce el comportamiento observado de las arenas
(p.ej. figura 2-8[78] y, más abajo, 4-27) Las características salientes que son correctamente
reproducidas son: i) el valor máximo de σ1/σ3 depende de la presión y de la densidad
relativa; ii) todas las simulaciones tienden a un mismo valor final;y iii) la tendencia a
dilatar crece con la densidad relativa y disminuye con el aumento de la presión media.
En las simulaciones mostradas se alcanzó una deformación global ∆H/H = 20%, que
102
1
2
3
4
5
6
0.00 0.05 0.10 0.15 0.20
3 100 kPaσ =
1ε
0.58 0.98 840 55 31min max s r ce e c p φ= = = = = °
1
3
σσ D
25%rD =
75%rD =
100%rD =
50%rD =
-0.12
-0.10
-0.08
-0.06
-0.04
-0.02
0.00
0.02
0.00 0.05 0.10 0.15 0.20
25%rD =
75%rD =
100%rD =3 100 kPaσ =
1ε
vε
50%rD =
Figura 4-10: Curvas σ1/σ3− ε1 y εv − ε1 en el punto D para cuatro densidades relativasy una presión de confinamiento constante σ3 = 100 kPa. Ensayo triaxial.
103
1
2
3
4
5
0.00 0.05 0.10 0.15 0.201ε
0.58 0.98 840 55 31min max s r ce e c p φ= = = = = °
3 1000 kPa σ =
D1
3
σσ
25%rD =
75%rD =100%rD =
50%rD =
-0.06
-0.04
-0.02
0.00
0.02
0.00 0.05 0.10 0.15 0.201ε
vε
3 1000 kPa σ =25%rD =
75%rD =
100%rD =
50%rD =
Figura 4-11: Curvas σ1/σ3− ε1 y εv − ε1 en el punto D para cuatro densidades relativasy una presión de confinamiento constante σ3 = 1000 kPa. Ensayo triaxial.
104
correspondió a deformaciones ε1 > 70% para el punto D. En la figura 4-12 se extiende
la figura 4-10 hasta ese nivel de deformaciones locales. Puede observarse que las curvas
σ1/σ3 aún no alcanzan el estado crítico teórico
(σ1/σ3)c = tan2 [π/4 + φc/2] = 3.12 (4.152)
1
2
3
4
5
6
0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70
25%rD =50%rD =
75%rD =
100%rD = 3 100 kPaσ =
1ε
0.58 0.98 840 55 31min max s r ce e c p φ= = = = = °
1
3
σσ D
-0.20
-0.15
-0.10
-0.05
0.00
0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70
25%rD =
75%rD =
100%rD =
3 100 kPaσ =
1ε
vε
50%rD =
Figura 4-12: Curvas σ1/σ3− ε1 y εv − ε1 en el punto D para cuatro densidades relativasy una presión de confinamiento constante σ3 = 100 kPa. ε1 = 70%. Ensayo triaxial.
Esto se debe a dos razones:
1. Las funciones de estado φ [p, e0] y φcv [p, e0] tienden al estado crítico de manera
105
asintótica por lo que, teóricamente, la deformación necesaria para alcanzar el estado
crítico es infinita.
2. La malla produce una fuerte localización de deformaciones que no puede ser
adecuadamente reproducida por un modelo axilsimétrico.
En la figura 4-13 se muestran la curva local en el punto D y global para Dr = 100%.
Puede apreciarse que la curva tensión-deformación global tiene una rama descendente
mucho más abrupta que la local. Como el algoritmo no está equipado con mecanismos
para la regularización de la deformación localizada, esta rama descendente depende de
la malla empleada para la modelización.
1
2
3
4
5
6
0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80
3 100 kPaσ =
1ε
1
3
σσ
D
ABCDEFG
3
P H
A H
δσ
−
Figura 4-13: Deformaciones locales vs globales para Dr = 100% y σ3 = 100 kPa. Ensayotriaxial.
En la figura 4-14 se muestra la malla deformada, el contorno de desplazamientos
totales, el contorno de la norma de la deformación total ‖ε‖ y el contorno de la relación
de vacíos para esta muestra. La figura se presenta únicamente con carácter informativo
y como prueba del funcionamiento de los algoritmos, puesto que el tamaño y tipo de
elementos y la configuración axilsimétrica invalidan la simulación para deformaciones
grandes.
106
Figura 4-14: Simulación de un ensayo de compresión triaxial con Dr = 100% yσ3 = 100 kPa. a) malla deformada; b) desplazamientos; c) ‖ε‖; y d) e0. Escala: colorado:mayor; azul: menor. Para e0 : colorado: emax; azul: emin.
La prueba de que el modelo alcanza el estado crítico se presenta con la simulación
de un ensayo triaxial ideal, sin restricciones al desplazamiento horizontal en los bordes,
en el que se alcanza una deformación principal mayor ε1 = 90%3 En la figura 4-15 se
presentan las curvas σ1/σ3 − ε1, εv − ε1, φ− ε1 y e0 − ε1.
4.5.3 Simulación de un ensayo de deformación plana
Se ejecutó la simulación numérica de un ensayo de deformación plana con la malla que
se muestra en la figura 4-16 (568 elementos de seis nodos, 1203 nodos en total). Los
parámetros materiales son los mismos que en los ejemplos anteriores. Las condiciones de
borde son: desplazamientos del borde inferior restringidos; bordes laterales con presión
constante σ3 = 100 kPa; desplazamientos verticales impuestos en el borde superior, y
desplazamientos horizontales nulos. Los parámetros materiales son los mismos que los
3El modelo está formulado con la cinemática de las deformaciones infinitesimales.
107
1
2
3
4
5
6
7
0.00 0.20 0.40 0.60 0.80 1.00
3 100 kPaσ =
1ε
1
3
σσ 100%rD =
-0.20
-0.15
-0.10
-0.05
0.00
0.00 0.20 0.40 0.60 0.80 1.001ε
vε3 100 kPaσ =
100%rD =
0.58 0.98 840 55 31min max s r ce e c p φ= = = = = °
30
35
40
45
0.00 0.20 0.40 0.60 0.80
0.58
0.68
0.78
0.88
0.98
1ε
3 100 kPaσ =
100%rD =
φ
cφ
0e
φ
0e
Figura 4-15: Ensayo de compresión triaxial sin restricción de cabezales, Dr = 100%,σ3 = 100 kPa.
108
de las simulaciones de ensayos triaxiales descriptos en la sección anterior, y la densidad
relativa es Dr = 75%. Se activó la opción de actualización de la malla durante el análisis.
En la figura 4-16 se presentan la malla deformada, el contorno de desplazamientos,
el contorno de la norma de deformación total ‖ε‖ y el contorno de relación de vacíos e0.
Puede apreciarse que el modelo captura correctamente la localización de deformaciones
y el cambio de densidad (cambio de e0) asociado al desarrollo de bandas de localización
en materiales dilatantes. El patrón de localización puede ser comparado con mediciones
experimentales como las de la figura 4-17.
Figura 4-16: Simulación de un ensayo de compresión plana conDr = 75% y σ3 = 100 kPa.a) malla deformada; b) desplazamientos; c) ‖ε‖; y d) e0. Escala: colorado: mayor; azul:menor. Para e0 : colorado: emax; azul: emin.
En la figura 4-18 se muestra la evolución de la densidad relativa a lo largo del proceso
de deformación.
109
Figura 4-17: Curva tensión-deformación para un ensayo de compresión plana y análisistomográfico de la evolución de bandas de localización para geomateriales [157].
110
Figura 4-18: Evolución de Dr con la deformación. ∆h/h0 = 0%, 0.5%, 1.0%, 1.5%, 2.0%y 2.5%. Azul: Dr = 100%. Colorado: Dr = 0%.
4.5.4 Simulación del comportamiento de una presa CFRD
sujeta a un sismo destructivo
En esta sección se presenta un problema típico de la ingeniería geotécnica en la que el
modelo introducido en este capítulo puede ser aplicado: El diseño sísmico de presas de
enrocado con pantalla de hormigón en la cara de aguas arriba (presas CFRD).
Las presas CFRD deben ser verificadas para la acción del denominado “terremoto de
seguridad”, que es una acción sísmica destructiva de muy baja probabilidad de ocurrencia.
El efecto de un terremoto de gran intensidad es la inestabilidad temporaria de los taludes
de la presa, lo que se traduce en un asentamiento del coronamiento. El objeto de la
verificación es, por lo tanto, estimar ese asentamiento con el fin de incrementar la altura
de la presa y minimizar así el riesgo de sobrepaso [24][25][26][45][106][133][141].
Existe un amplio rango de procedimientos de cálculo para este problema, desde
métodos analíticos simples hasta complejos modelos numéricos tridimensionales. Los
métodos analíticos desacoplan el análisis de propagación de ondas en el cuerpo de la
111
presa - que produce amplificación de la aceleración basal - del análisis del efecto de
esta aceleración sobre la estabilidad de los taludes [41][86][102][107][150][134]. Entre
los métodos numéricos, algunos resuelven el problema de manera desacoplada - modelo
elástico para propagación de ondas y modelo elastoplástico para estabilidad - mientras
que otros resuelven el problema acoplado propagación-equilibrio mediante modelos
elastoplásticos.
En el caso particular de la Presa Los Caracoles, recientemente construida en San
Juan, Argentina, se efectuaron dos determinaciones independientes del asentamiento
del coronamiento para su terremoto de seguridad: una predicción analítica [107] y
una predicción numérica desacoplada [9]. Mas tarde, Sfriso [138] desarrolló un método
numérico simplificado basado en el modelo que se presenta en este capítulo y lo aplicó al
mismo problema. En la figura 4-19 se presenta una sección transversal de la presa [9], en
la figura 4-20 se presenta la malla de elementos finitos. Los parámetros adoptados son
[9][138]
emin emax cs ce m φc pr e0
3B 0.10 0.45 650 2.17 0.5 41.3◦ 5.7 0.16
3D 0.10 0.45 600 2.17 0.5 41.0◦ 2.5 0.20
3L 0.10 0.40 650 2.17 0.5 36.5◦ 3.7 0.20
Lecho 0.12 0.46 700 2.17 0.5 37.1◦ 4.0 0.20
En la figuras 4-21 y 4-22 se muestran la distribución de G y φ. Puede observarse que
la subdivisión de la malla responde exclusivamente a la zonificación de la presa, puesto
que no es necesario efectuar zonificaciones adicionales para simular la dependencia de G
y φ respecto de la presión. Se aplicó una aceleración basal compuesta por siete pulsos
sinusoidales de aceleración máxima PGA = 0.8g en dirección horizontal, con un período
T = 0.56s [138]. Luego se simuló la oscilación libre de la presa por un lapso de 7T.
El asentamiento calculado por Núñez [107] es 2.37 m − 3.39 m; el calculado por Bissio y
Tejada [9] es 2.13 m; el calculado con este procedimiento es 2.38 m [138]. En la figura 4-23
se muestra el contorno de deformaciones plásticas al fin de la simulación.
112
Figura 4-19: Sección transversal de la presa Los Caracoles, San Juan, Argentina [9].
Figura 4-20: Malla de elementos finitos para el análisis sísmico de la presa Los Caracoles[138].
Figura 4-21: Distribución de módulo de corte en la sección del modelo. Azul: G =200 MPa; Colorado: G = 1300 MPa [138].
113
Figura 4-22: Distribución de φ en la sección del modelo. Azul: φ = 41◦; Colorado:φ = 48◦ [138].
Figura 4-23: Análisis dinámico de Presa Los Caracoles. a) Malla deformada, δ = 2.38 m.b) Contorno de desplazamientos totales ‖δ‖ . c) Contorno de deformaciones totales ‖ε‖ .Azul: menor. Colorado: mayor.
114
4.A Apéndice: Deducciones de algunas funciones de
estado
4.A.1 Hiperelasticidad
Deducción de D a partir de Ws
La función de energía complementaria de deformación es
Ws =p2
G [p]
(S
(1−m) (2−m) +r2
4
)(4.153)
donde
G = Gref (p/pref)m pref (4.154)
Gref , pref , m y S no son función de σ. Ws se puede poner como
Ws =p2−mpm−1
ref
Gref
(S
(1−m) (2−m) +τ2
4p2
)(4.155)
donde
τ = ‖s‖ (4.156)
La derivada de Ws respecto de σ es el tensor de deformación
ε = ∂Ws/∂σ (4.157)
Esta derivada se expande como
∂Ws
∂σ=∂Ws
∂p
∂p
∂σ+∂Ws
∂τ
∂τ
∂σ(4.158)
115
donde
∂Ws
∂p= (2−m)
p1−mpm−1ref
Gref
(S
(1−m) (2−m) +τ2
4p2
)−p−mpm−1
ref
Gref
τ 2
2p
=1
G
(S
1−mp−mτ 2
4p
)(4.159)
∂Ws
∂τ=
τ
2G
∂p
∂σ=1
31
∂τ
∂σ=s
τ(4.160)
Simplificando
∂Ws
∂σ=1
G
(S
1−m −m τ2
4p2
)p
31+
τ
2G
s
τ
=
((S
1−m −mr2
4
)1
3GIesf +
1
2GIdev
): σ (4.161)
Esta ecuación se invierte a
σ =
(3
GS
1−m−mr2
4
Iesf + 2GIdev
): ε (4.162)
donde
De = 3G
S1−m
−mr2
4
Iesf + 2GIdev (4.163)
es el operador elástico. Si se define
K =G
S1−m
− m4r2
(4.164)
queda finalmente
De = 3KIesf + 2GIdev (4.165)
116
Deducción de σ
En plasticidad convencional con rigidez constante, el incremento de tensión es
σ = D : εe (4.166)
por lo que la tensión σe - variables plásticas congeladas - se calcula asumiendo un paso
elástico con la expresión
σe = D : (nεe +∆ε) (4.167)
Cuando el tensor elástico D depende de σ, el incremento de tensión es
σ = D : ε+ ε : ∂D/∂σ : σ (4.168)
que se puede reordenar como
σ = T−1 : D : ε (4.169)
donde
T = I− ε : ∂D/∂σ (4.170)
La determinación de T requiere el cálculo de ε : ∂D/∂σ. Como
∂D/∂σ = 311∂K/∂σ + 2Idev∂G/∂σ (4.171)
puede escribirse
ε : ∂D/∂σ = εv1∂K/∂σ + 2εd∂G/∂σ (4.172)
Pero∂G
∂σ=mGref
3
(p
pref
)m−1
1 =mG
3p1 (4.173)
117
∂K
∂σ=
∂G/∂σS
1−m− m
4r2+
mrG
2p(
S1−m
− m4r2)2(rr− r31)
=mK
p
((1− r2
2S1−m
− m2r2
)1
31 +
12S1−m
− m2r2r
)
Entonces puede escribirse
ε :∂D [σ]
∂σ= εv1
mK
p
((1 +
r2
m2r2 − 2S
1−m
)1
31+
12S1−m
− m2r2r
)+ 2εd
mG
3p1 (4.174)
Como σ = De : ε, se puede poner
p = Kεv (4.175)
r = 2Gεd/p (4.176)
por lo que queda
ε :∂D [σ]
∂σ=
(1 +
r2
m2r2 − 2S
1−m
)m
311+
m2S1−m
− m2r21r+
m
3r1 (4.177)
lo que permite escribir finalmente
T =
(1−m+ K
2Gmr2)Iesf + Idev − K
2Gm1r− m
3r1 (4.178)
Deducción de dK/dp
La derivada del módulo volumétrico respecto de la presión es
dK
dp=mK
p
(1 + (m− 1) Kr
2
2G
)(4.179)
donde las derivadas parciales son
∂K
∂p=mK
p
∂K
∂r=mK2
2Gr
∂G
∂p=mG
p(4.180)
118
∂r
∂p=∂
∂p
[2G ‖εe
d‖p
]= 2 ‖εe
d‖G(m− 1p2
)=m− 1p
r (4.181)
4.A.2 Función de fluencia
Linealización de la función de fluencia
Como muestra la figura 4-4, la ec. (4.44) define múltiples superficies en el espacio {σ},de las cuales sólo una tiene sentido físico, la orientada en el eje positivo de la presión
hidrostática. la figura 4-24 muestra Ff y ‖n‖ en función de r para compresión triaxial.
Puede observarse que Ff tiene un máximo local cuando ‖n‖ = 0 y un segundo cero en
r =√6.
Figura 4-24: Función de fluencia y norma de su normal en función de r para tres valoresde µf .
El problema del segundo cero se salva si la ec. (4.44) se redefine según
r < rlim → Ff =µf + 6
2r2 −
(µf + 9
)J3r − µf = 0 (4.182)
r > rlim → Ff =µf + 6
2rlimr −
(µf + 9
) r2limr2J3r − µf = 0 (4.183)
donde rlim es una oblicuidad umbral. En la figura 4-25 se muestra la función de fluencia
suavizada con este procedimiento.
119
Figura 4-25: Ec. (4.44) en su forma original y suavizada para el vértice de compresióntriaxial.
Derivadas de Ff
Las derivadas de Ff son
Ff,µ = −1 +1
2r2 − J3r (4.184)
n = Ff,σ =1
p
((r2 +
(µf + 9
)J3r)1 +(µf + 6
)r−(µf + 9
)r · r)
(4.185)
nd =1
p
((µf + 9
) r231+(µf + 6
)r−(µf + 9
)r · r)
(4.186)
µ,σ = −Ff,σ
Ff,µ
=n
1− 12r2 + J3r
(4.187)
∂n
∂σ=1
p2
−((µf + 12
)r2 + 4
(µf + 9
)J3r)11/3
+(µf + 6
)Idev + 2 (r1+ 1r)
+(µf + 9
)(r · r1+ 1r · r)
−(µf + 9
)(δikrlj + δjlrik) eiejekel
120
4.A.3 Regla de flujo
Demostración de que β es único
Hipótesis
1. Los autovalores de σ están ordenados σ1 ≥ σ2 ≥ σ3 > 0.
2. Los autovalores de εp están ordenados εp1 ≥ εp
2 ≥ εp3.
3. Existe coaxialidad deviatórica. md =Idev:∂F/∂σ
‖Idev:∂F/∂σ‖ con autovalores md1 ≥ md2 ≥md3, por lo que
m2d1 +m
2d2 +m
2d3 = 1 (4.188)
md1 +md2 +md3 = 0 (4.189)
md1 ≥√1/6 (4.190)
md3 ≤ −√1/6 (4.191)
4. Existe no asociatividad volumétrica m =md + β1.
5. No todos los componentes {εp1, ε
p2, ε
p3} pueden ser del mismo signo, por lo que
εp1 = λ (md1 + β) > 0 (4.192)
εp3 = λ (md3 + β) < 0 (4.193)
Tesis
∀σ ∃β y es único. β está dado por
md2 + β ≤ 0→ β = −σ1md1 +Ncvσ2md2 +Ncvσ3md3
σ1 +Ncvσ2 +Ncvσ3(4.194)
md2 + β ≥ 0→ β = −σ1md1 + σ2md2 +Ncvσ3md3
σ1 + σ2 +Ncvσ3(4.195)
121
Demostración
Se asumen dos valores β1 y β3 distintos para β
md2 + β1 < 0→ β1 = −σ1md1 +Ncvσ2md2 +Ncvσ3md3
σ1 +Ncvσ2 +Ncvσ3(4.196)
md2 + β3 ≥ 0→ β3 = −σ1md1 + σ2md2 +Ncvσ3md3
σ1 + σ2 +Ncvσ3(4.197)
El numerador de β3 es mayor que el de β1 y el denominador es menor. Por lo tanto,
β3 < β1. Entonces, hay tres posibilidades
1. md2 + β1 ≤ 0 =⇒ md2 + β3 ≤ 0 =⇒ β = β1 y elimina β3.
2. md2 + β3 ≥ 0 =⇒ md2 + β1 ≥ 0 =⇒ β = β2 y elimina β1.
3. md2 + β3 ≤ 0 ∧md2 + β1 ≥ 0.
(a) Se reemplaza β1 en md2 + β1 ≥ 0. Queda
md2 + β1 ≥ 0→ (md2 −md1) σ1 + (md2 −md3)Ncvσ3 ≥ 0 (4.198)
(b) Se reemplaza β3 en md2 + β3 ≤ 0. Queda
md2 + β3 ≤ 0→ (md2 −md1) σ1 + (md2 −md3)Ncvσ3 ≤ 0 (4.199)
(c) a∧b implican
(md2 −md1) σ1 + (md2 −md3)Ncvσ3 = 0 (4.200)
que se reordena como
md1σ1 +md3Ncvσ3 = md2Ncvσ3 +md2σ1 (4.201)
122
(d) Se reemplaza en el primer término en β1 y β3
β1 = −σ1 +Ncvσ2 +Ncvσ3σ1 +Ncvσ2 +Ncvσ3
md2 → β1 = −md2 (4.202)
β3 = −σ1 + σ2 +Ncvσ3σ1 + σ2 +Ncvσ3
md2 → β3 = −md2 (4.203)
O sea
(md2 + β3 ≤ 0) ∧ (md2 + β1 ≥ 0) =⇒ β1 = β3 = −md2 =⇒ β es único (4.204)
4.A.4 Determinación de pr y φc para arenas del Puelchense
En este apartado se demuestra la técnica de calibración del modelo para las arenas del
Puelchense. El material es una arena silícea de origen marino del terciario superior con
granos subredondeados. Se ensayó la fracción pasante por el tamiz #40 y retenida en el
tamiz # 200 de muestras recuperadas durante la ejecución de ensayos SPT en la localidad
de Campana. Las propiedades índice son: γs = 26.6 kN/m3, emin = 0.58, emax = 0.89.
Clasificación SUCS: SP.
Las muestras sueltas se prepararon mediante vuelco de material seco a través de un
embudo de 6 mm de boca desde 50−80 mmde altura, con lo que se obtuvo 32% < Dr < 42%.
Las muestras densas se prepararon mediante compactación manual en cinco capas hasta
alcanzar Dr > 95%. Todas las probetas tuvieron 38 mm±3% de diámetro y 100 mm±5%de
altura.
El ensayo consistió en la saturación de las muestras con agua desaireada; consolidación
isotrópica hasta la presión de ensayo; compresión drenada a σ3 constante y ε =
0.5 mm/min ± 10%. Se emplearon cabezales planos rígidos sin membranas deslizantes;
medición de deformación entre cabezales con flexímetro de sensibilidad 10µ/div;
medición de carga con aro dinamométrico de sensibilidad 2 N/div; medición de volumen
con bureta calibrada conectada a ambas cabezas de sensibilidad 10 mm3/div; registro
manual de resultados; cámara triaxial operada con aire comprimido de 5 dm3 y presión
123
máxima de 1.0 MPa construida por el autor.
En la figura 4-26 se muestra el ángulo de fricción interna correspondiente a cada uno
de los ensayos y el ajuste de la ec. (4.37) con φc = 30◦ y pr = 50.
-2
0
2
4
6
8
10
12
14
100 1000
Dr=35Dr=97
ψ
[ ]kPap500
0.58 0.89 30 50min max c re e pφ= = = ° =
Figura 4-26: Calibración de φ y pr para arena Puelchense.
En las figuras 4-27 y 4-28 se muestran las curvas σ1/σ3−ε1 y εv−ε1 correspondientesa todos los ensayos.
124
1.00
2.00
3.00
4.00
5.00
6.00
0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10
600 400 300 200 100 75 50 25
1
3
σσ
1ε
Presión de confinamiento
-5.0%
-4.0%
-3.0%
-2.0%
-1.0%
0.0%
1.0%
0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.101ε
vε
Figura 4-27: Ensayos triaxiales convencionales de arena Puelchense densa (95% < Dr <100%).
125
1.00
2.00
3.00
4.00
0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10
600 400 300 200 100 75 50 25
1ε
1
3
σσ
Presión de confinamiento
-2.0%
-1.0%
0.0%
1.0%
2.0%
3.0%
4.0%
0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.101ε
vε
Figura 4-28: Ensayos triaxiales convencionales de arena Puelchense suelta (Dr ≈ 35%).
126
Capítulo 5
Modelización de la carga
monotónica de corte
5.1 Introducción
En esta sección se presenta un modelo constitutivo para la simulación del comportamiento
de arenas en carga de corte monotónica. El modelo está orientado a problemas en los
que interesa el desarrollo de plasticidad previo a la falla, como por ejemplo: i) problemas
cuasi-elásticos de interacción suelo-estructura como la propagación de ondas mecánicas o
el diseño de fundaciones de máquinas; ii) la deformación de pilotes, zapatas y muros
de contención sujetos a carga estática; iii) la deformación de taludes, terraplenes y
excavaciones; y iv) otros problemas en los que las trayectorias de descarga, si existen,
son pequeñas frente a la magnitud de las deformaciones en carga.
El modelo que se presenta en este capítulo se construyó sobre el modelo presentado
en el Capítulo 4, al que se le agrega una nueva variable de estado y su función de
endurecimiento asociada.
127
5.2 Estructura matemática del modelo
En esta sección se repite la estructura matemática del modelo de falla con corte y se
agregan los elementos nuevos de este modelo.
5.2.1 Cinemática
Se adopta la descomposición aditiva de las deformaciones infinitesimales
εv = εe + εp (5.1)
5.2.2 Variables de estado
El estado del material se caracteriza mediante el tensor de tensiones efectivas σ, la
relación de vacíos para presión nula e0 y el parámetro escalar de movilización de
resistencia al corte ς.
5.2.3 Relación tensión-deformación
Se mantiene la función de energía complementaria de deformación
Ws = Ws [σ, e0] (5.2)
con la que se calcula el tensor de deformación elástica
εe = ∂Ws/∂σ (5.3)
5.2.4 Función de fluencia
La función de fluencia es de la forma
Fs [σ, e0, ς] = 0 (5.4)
128
5.2.5 Regla de flujo
El incremento de deformación plástica se calcula con
εp = λm (5.5)
donde
m =md + β1 (5.6)
β es un factor de dilatancia,
md = nd/ ‖nd‖ (5.7)
nd = Idev : n (5.8)
n = ∂Fs/∂σ (5.9)
5.2.6 Condición de plasticidad
Se emplea la condición
λ =
λ si: Fs = 0 y: Fs > 0
0 en otro caso(5.10)
5.2.7 Ecuaciones de evolución
La evolución de e0 está dada por la expresión
e0 = − (1 + e0) εpv (5.11)
La evolución de ς está dada por la expresión
ς = h ‖εpd‖ (5.12)
donde h [σ, e0, ς] es una función de endurecimiento.
129
5.3 Funciones de estado
5.3.1 Variable de estado y función de fluencia
Aunque el criterio de Mohr-Coulomb dado por la ec. (2.3) es un criterio de falla, se
lo emplea como función de fluencia mediante un cambio de variables, de φ al “ángulo
(equivalente) de fricción interna movilizado” φmob. Queda
FMC = τ/σn − tan [φmob] = 0 (5.13)
En la ec. (5.13), φmob actúa como una variable interna que endurece durante la carga
elastoplástica y alcanza el valor φ cuando el material falla.
φmob no puede ser utilizado como variable interna cuando la superficie de falla es
curva en el espacio τ − σ (cuando φ no es constante), como lo prueba un ensayo triaxial
conceptual que puede seguirse con ayuda de la figura 5-1 [140].
1. Una muestra de arena es consolidada isotrópicamente hasta una presión de cámara
σ1 = σ3 = σ3a.
2. Se aplica un esfuerzo desviador hasta alcanzar un estado cercano a la falla (círculos
a a c); la variable interna φmob se endurece para mantener FMC = 0;
3. Se disminuye σ1 levemente hasta σ1d, por lo que la muestra queda en estado elástico,
y FMC < 0 (círculo d)
4. Se aplica una trayectoria de compresión proporcional manteniendo constante la
oblicuidad σ1/σ3 (círculos d a f), siempre por debajo de la superficie de fluencia,
que queda fija con apertura φmob.
5. A pesar de que d − f es una trayectoria elástica, se encuentra la falla cuando el
círculo toca la superficie de falla curva en el punto f . Este es un absurdo: se alcanza
la falla por corte en una compresión proporcional elástica con FMC < 0.
130
σ
τ
f
Superficie de fallaDominio elástico
mobφ
e
d
cb
a
f
1aσ3aσ 1cσ1dσ
Figura 5-1: Prueba que φmob no puede ser utilizado como variable interna para superficiede falla curvas.
Puede definirse un factor de movilización de resistencia al corte ς como una variable
de estado independiente, y postularse que la superficie de fluencia es
FMC = τ/σn − ς tan [φ] = 0 (5.14)
Por simple inspección de las ec. (5.13) y (5.14), es evidente que φmob y ς están relacionados
entre sí mediante la ecuación
tan [φmob] = ς tan [φ] (5.15)
y que 0 ≤ ς ≤ 1.
Con esta nueva definición de función de fluencia, la compresión proporcional
(trayectoria de incremento de presión con σ1/σ3 constante) no es una trayectoria
elástica, ya que implica un crecimiento monotónico de la variable de estado ς (o sea,
endurecimiento plástico) para mantener φmob constante con φ decreciente.
Para el empleo de ς junto con el criterio de Mohr-Coulomb para estados
131
tridimensionales de tensiones debe definirse
FMC = σ1/σ3 − tan2[π/4 + tan−1 [ς tan [φ]] /2
]= 0 (5.16)
que es una función complicada y poco elegante. La variable de estado ς, en cambio,
puede introducirse de manera expeditiva en el criterio de Matsuoka-Nakai. La función
de fluencia se define como
Fs = (µ+ 6)J2r − (µ+ 9)J3r − µ = 0 (5.17)
que es similar a la ec. (4.44) pero en donde se reemplazó µf por [140]
µ = ς2µf (5.18)
5.3.2 Relación tensión-deformación
La relación tensión-deformación más empleada para arenas en carga monotónica es el
modelo hiperbólico dado por la ec. (2.63). Sin embargo, este modelo presenta dos
inconvenientes: i) subestima la rigidez a baja deformación; y ii) emplea el parámetro
Rf < 1, de dudoso significado físico pero necesario para evitar que la carga máxima se
alcance para deformaciones infinitas. Tatsuoka [142][143] presentó resultados de ensayos
torsionales sobre muestras cilíndricas de arena Toyoura y propuso el empleo de tres
funciones hiperbólicas [142][145] como se muestra en la figura 5-2. Aunque la propuesta
de Tatsuoka resuelve el problema, lo hace a costa de un número excesivamente grande
de parámetros que no tienen un claro significado físico.
La información experimental firme es que la rigidez inicial está dada por la ec. (4.21)
y que la rigidez tangente en falla es nula, por lo que la rigidez tangente elastoplástica
132
0
0.2
0.4
0.6
0.8
1
0.0 0.2 0.4 0.6 0.8 1.0
Experimental, NC TSExperimental, OC TSHiperbólico 1Hiperbólico 2Hiperbólico 3
sG
G
ς
Figura 5-2: Ensayos de torsión de de arena Toyoura normalmente consolidada (NC) ysobreconsolidada (OC) y tres funciones hiperbólicas para su simulación [142][143][145].
Gt = ∂τ/∂γ debe cumplir con
Gt⌋ς=0 = G (5.19)
Gt⌋ς→1 = 0 kPa (5.20)
Una ecuación simple que cumple con estas restricciones es [140]
Gt =1− ς1 + ας
G (5.21)
donde α es una función auxiliar de calibración. El módulo secante Gs = τ/γ
correspondiente es [140]
Gs =−ς
(1 + α) ln [1− ς] + αςG (5.22)
En la figura 5-3 se repite la información de la figura 5-2 con el agregado de la ec.
(5.22) con α = 10 (curva llena). Pese a que en la figura no es evidente por la escala, la
ec. (5.22) vale exactamente 0 para ς = 1.0.
Con la ec. (5.21) y la relación Gp = GGt/ (G−Gt) puede calcularse el módulo
133
0
0.2
0.4
0.6
0.8
1
0.0 0.2 0.4 0.6 0.8 1.0
Experimental, NC TSExperimental, OC TSEste modeloHiperbólico 1Hiperbólico 2Hiperbólico 3
sG
G
ς
Figura 5-3: Ensayos experimentales, ec. (5.22) y comparación con las tres funcioneshiperbólicas [140][142][143][145]
plástico
Gp = 22τ
‖εpd‖=1/ς − 11 + α
G (5.23)
5.3.3 Función de endurecimiento para ς
Para convertir la ec. (5.23) en una función de endurecimiento se restringe el análisis a
una trayectoria de compresión triaxial a presión constante, para la que vale la expresión
aproximada
µ ≈ 3r2 (5.24)
y r = τ /p. La evolución de µ se simplifica entonces a
µ = 6r
pτ (5.25)
Como τ = 2Gp ‖εpd‖ queda
µ = 12
õ
3
Gp
p‖εp
d‖ (5.26)
134
De la ec. (5.18),
ς =√µ/µf (5.27)
por lo que
ς =1
2√µµf
µ (5.28)
La combinación de las ecs. (5.26), (5.28) y (5.23) arroja
ς =
√12
µf
1/ς − 11 + α
G
p‖εp
d‖ (5.29)
que puede escribirse como
ς = h ‖εpd‖
donde
h =
√12
µf
1/ς − 11 + α
G
p(5.30)
es la función de endurecimiento buscada [140].
La simulación numérica de ensayos triaxiales drenados permite establecer la siguiente
relación aproximada entre la variable auxiliar α y el cociente E50/E
E50E
=12 +Dr
13 + 4α(5.31)
El cociente E50/E depende de la densidad relativa [85]. En este modelo se propone
la interpolación lineal [140]
E50E
= cd − 2 (1−Dr) (cd − cl) (5.32)
que tiene dos parámetros materiales
cd = (E50/E)Dr=100%(5.33)
cl = (E50/E)Dr=50%(5.34)
135
que se miden en ensayos triaxiales. Aunque existe, se desprecia la débil dependencia que
el cociente E50/E tiene con respecto a la presión media dentro del rango de presiones de
interés ingenieril [84].
5.3.4 Resistencia máxima y evolución de variables de estado
La ec. (5.17) tiene una sola función interna, µ, que ahora depende de dos variables de
estado. Por lo tanto, la condición de falla µ [p, e0, ς] = 0 ahora se expande a
µ =∂µ
∂ςς +
∂µ
∂φ
∂φ
∂ψ
(∂ψ
∂e0e0 +
∂ψ
∂pp
)= 0 (5.35)
Para presión constante queda
µ =∂µ
∂ςς +
∂µ
∂e0e0 = 0 (5.36)
Las derivadas de la ec. (5.36) cumplen con
∂µ/∂ς > 0 (5.37)
∂µ/∂φ > 0 (5.38)
Mientras que el primer término de la ec. (5.36) es no negativo, el signo del segundo
término depende de e0. Se presentan tres situaciones que se analizan a continuación.
Muestra suelta
Una muestra suelta alcanza el estado ς = 1 para una relación de vacíos e0 > e0c. Luego
de ese punto se deforma en falla mientras reduce e0, por lo que el segundo término de la
ec. (5.36) es positivo hasta que se anula en el estado crítico. Por lo tanto, la ec. (5.36)
sólo se satisface cuando e0 = 0, o sea, la resistencia máxima y el estado crítico ocurren
simultáneamente para muy grandes deformaciones.
136
Muestra con la relación de vacíos crítica
Una muestra que alcance la condición ς = 1 con una relación de vacíos e0 = e0c satisface
la ec. (5.36) sin cambios posteriores de volumen. Por lo tanto, la resistencia máxima y
el estado crítico se alcanzan simultáneamente para una deformación moderada.
Muestra densa
Una muestra densa dilata (e0 > 0) durante la etapa de endurecimiento de ς (ς < 1). En
este caso, la evaluación de la ec. (5.36) arroja
ς = −∂µ/∂e0∂µ/∂ς
e0 > 0 (5.39)
Por lo tanto, la resistencia máxima ocurre cuando la tasa de crecimiento de Fs debido
al endurecimiento de ς iguala a la tasa de ablandamiento producida por la dilatancia.
Conclusión
Se concluye entonces que la condición µ = 0 ocurre luego que la condición ς = 1 para
arenas sueltas, simultáneamente con la condición ς = 1 para arenas en el estado crítico
y antes de ς = 1 para arenas densas. Debe notarse que las muestras más sueltas que el
estado crítico tienen ψ < 0 y por lo tanto φ < φc. En otras palabras, estas muestras
sólo pueden alcanzar el estado crítico luego de una disminución de volumen, dado que
una deformación puramente isocórica implicaría la licuefacción estática con resistencia
residual nula.
137
5.4 Implementación numérica
5.4.1 Operador de integración
El operador de integración es implícito. Se resuelve el conjunto de ecuaciones
n+1σ = D
[n+1σ, n+1e0
]: (nε+∆ε−∆εp) (5.40)
n+1e0 = (1 + ne0) exp[
n+1∆εp : 1]− 1 (5.41)
n+1ς = nς + n+1h[
n+1σ, n+1e0,
n+1ς] ∥∥Idev : ∆εp
∥∥ (5.42)
sujeto a la restricción
Fs
[n+1σ, n+1µ
]= 0 (5.43)
donde n+1µ = µ [ n+1σ, n+1e0,
n+1ς] .
Las ecs. (5.40) a (5.43) están expresadas en términos del tensor de deformación
plástica
∆εp = λ(md
[n+1σ, n+1e0,
n+1ς]+ β[
n+1σ, n+1e0
]1)
(5.44)
que tiene una única variable independiente λ, por lo que, en un paso elastoplástico, el
conjunto de ecuaciones (5.40) a (5.42) se reduce a
n+1σ = σ [λ] (5.45)
n+1e0 = e0 [λ] (5.46)
n+1ς = ς [λ] (5.47)
con lo que la ec. (5.43) se reduce a
Fs [λ] = 0 (5.48)
138
5.4.2 Algoritmo de integración
El algoritmo de integración de este modelo está construido sobre la base del algoritmo
presentado en el capítulo 4, al que se agregan los siguientes elementos:
1. El estado convergido es
{ nσ, ne0,
nς} (5.49)
2. Al inicio del paso se agrega el cálculo de
nµ = nς2 nµf (5.50)
3. Cambia el Algoritmo C, que ahora actualiza las variables
σ(i) λ(i) ε
p(i) ς(i) (5.51)
con
D(i) e(i)0 µ(i)f (5.52)
congelados.
Algoritmo C: Actualización de σ con funciones de estado constantes
1. Se calcula el versor de deformación plástica deviatórica m(0)d
2. Se corrige el multiplicador plástico proyectándolo sobre la nueva dirección plástica:
λ(0) = εp(i) :m(0)d (5.53)
3. Se calcula
s(0) = s− 2G λ(0) m(0)d (5.54)
139
4. Se calcula una cota a ∆λ para que se cumpla s : s > 0.5 ‖s‖ ‖s‖ (se limitan las
oscilaciones del algoritmo)
∆λlim =‖s‖ − 0.5 ‖s‖m(0)d : s
‖s‖2G
(5.55)
5. Se calcula β(0) y p(0) mediante un algoritmo de bisección (Algoritmo E)
6. Se calcula
σ(0) = σ − λ(0)
(s(0) + p(0)1
)(5.56)
7. Mientras
abs[F (0)s
]> tol (5.57)
(a) Se calcula la función de endurecimiento h(i)
(b) Se calcula
∆λ(i+1) = min
[F(i)s
n(j) : D :m(j) − Fs,ς h(i),∆λlim
](5.58)
(c) Se actualiza
λ(i+1) = λ(i) +∆λ(i+1) (5.59)
s(i+1) = s− 2G λ(i+1) m(i)d (5.60)
(d) Se actualiza p(i+1) y β(i+1) mediante el Algoritmo E
(e) Se actualiza
σ(i+1) = s(i+1) + p(i+1)1 (5.61)
(f) Se actualiza ς(i+1) mediante un Newtor-Raphson local con λ constante
(Algoritmo F)
140
(g) Se actualiza
µ(i+1) =(ς(i))2µf (5.62)
(h) Se evalúa la función de fluencia
F (i+1)s ⋚ 0 (5.63)
(i) Se itera hasta convergencia.
Algoritmo E: Actualización de {p, β}
El objetivo del algoritmo es encontrar
p(i+1) = p− 3Kβ(i+1) [σ, e0]λ (5.64)
con λ y K fijos. Dado que β [σ, e0] es una función altamente no lineal de p, se opta por
un algoritmo de bisección como sigue:
1. Se definen cotas p(0)infy p(0)sup suficientemente amplias
2. Se inicializan dos variables auxiliares
p(0)1 = p(0) p
(0)2 = p(0) (5.65)
3. Se calcula
β(0)1 = β
[p(0)1
](5.66)
4. Se define el error
ζ(0)1 = p
(0)1 −
(p− 3Kβ(0)1 λ
)(5.67)
141
5. Se corrigen las cotas
ζ(0) < 0 =⇒ p(0)inf = max
[p(0)inf , p
(0)1
](5.68)
ζ(0) > 0 =⇒ p(0)sup = min[p(0)sup, p
(0)1
](5.69)
6. Se calcula
p(0)2 = p− 3Kβ(0)1 λ (5.70)
7. Se controlan cotas
p(0)2 > p(0)sup ∨ p
(0)2 > p(0)sup =⇒ p
(0)2 =
p(0)sup + p
(0)inf
2(5.71)
8. Se calcula
β(0)2 = β
[p(0)2
](5.72)
9. Se calcula el error
ζ(0)2 = p
(0)2 −
(p− 3Kβ(0)2 λ
)(5.73)
10. Se corrigen cotas
ζ(0)2 < 0 =⇒ p
(0)inf = max
[p(0)inf , p
(0)2
](5.74)
ζ(0)2 > 0 =⇒ p(0)sup = min
[p(0)sup, p
(0)1
](5.75)
11. Se inicializa
p(0) = p(0)2 β(0) = β
(0)2 ζ(0) = ζ
(0)2 (5.76)
12. Mientras
abs
[ζ(i)]> tol (5.77)
142
(a) Se interpola linealmente
ζ(0)2 �= ζ
(0)1 =⇒ p(i+1) =
ζ(i)1 p
(i)2 − ζ(i)2 p(i)1ζ(i)2 − ζ
(i)1
(5.78)
ζ(0)2 = ζ
(0)1 =⇒ p(i+1) = p− 3Kβ(i)λ (5.79)
(b) Si
p(i+1) > p(i)sup ∨ p(i+1) < p(i)inf =⇒ p(i+1) =p(i)sup + p
(i)inf
2(5.80)
(c) Se actualiza
β(i+1) = β[p(i+1)
](5.81)
(d) Se actualiza
ζ(i+1) = p(i+1) −(p− 3Kβ(i+1)λ
)(5.82)
(e) Se corrigen cotas
ζ(i+1) < 0 =⇒ p(i+1)inf = max
[p(i)inf , p
(i+1)]
(5.83)
ζ(i) > 0 =⇒ p(i+1)sup = min[p(i)sup, p
(i+1)]
(5.84)
(f) Se actualiza
p(i+1)1 = p
(i+1)2 ζ
(i+1)1 = ζ
(i+1)2 p
(i+1)2 = p(i+1) ζ
(i+1)2 = ζ(i+1) (5.85)
(g) Se itera hasta convergencia.
13. Se informa p(i+1), β(i+1) y FIN
143
Algoritmo F: Actualización de ς
El objetivo del algoritmo es encontrar
ς (i+1) = nς + h[σ, e0, ς
(i+1)]λ (5.86)
con λ fijo. Se emplea el algoritmo de Newton-Raphson convencional:
1. Se calcula el residuo
ζ(0) = ς (0) −(
nς + h(0) λ)
(5.87)
2. Mientras
ζ(i) > tol (5.88)
(a) Se calcula
∆ς = −ξ(i)/(1− λ (∂h/∂ς)(i)
)(5.89)
(b) Se actualiza
ς (i+1) = max[
nς, min[1, ς(i) +∆ς
]](5.90)
h(i+1) = h[p, ς(i+1)
](5.91)
ζ(i+1) = ς(i+1) −(
nς + h(i+1) λ)
(5.92)
(c) Se itera hasta convergencia.
3. Se informa ς (i+1), h[σ, e0,ς
(i+1)]y FIN.
144
5.5 Resultados
5.5.1 Lista de parámetros
Los parámetros del modelo son
Par. unidad Descripción
emin − relación de vacíos mínima
emax − relación de vacíos máxima
cs − parámetro de rigidez al corte
φc◦ ángulo de fricción interna crítico
pr − parámetro de dureza de grano
cd − parámetro de endurecimiento (E50/E)Dr=100%
cl − parámetro de endurecimiento (E50/E)Dr=50%
Excepto por emin y emax, que son determinados por procedimientos normalizados,
todos los parámetros pueden ser calibrados en una serie de ensayos de compresión triaxial
drenada con medición local de deformaciones.
Existen cuatro parámetros adicionales con valores por defecto. En la ec. (4.21),
ce = 2.17 y m = 0.5; en la ec. (4.37), ∆φ = 3◦ y b = 2◦. Si se requiere, estos parámetros
pueden ser ajustados para un cierto material mediante una serie más extensa de ensayos
triaxiales drenados y ensayos de columna resonante.
5.5.2 Análisis de parámetros materiales
Se evalúa el efecto que los cinco parámetros materiales cs, φc, pr, cd y cl tienen
sobre la respuesta del modelo, para ensayos de compresión triaxial drenada. Para los
parámetros materiales se eligieron valores previamente publicados para la arena Toyoura
[12][28][29][44][66][109][110][130][142][143][145][154][158][164] y se los hizo variar dentro
de su rango posible. En particular, para cd y cl estimaron valores en base a la información
145
publicada en [84][85][142][143][145]. La lista es
Par. min real max
cs 640 840 1040
φc 28◦ 31◦ 34◦
pr 35 55 75
cd 0.50 0.40 0.30
cl 0.35 0.25 0.15
Se empleó una presión de confinamiento constante σ3 = 100 kPa y dos densidades
relativas, Dr = 50% y Dr = 100%. Las figuras 5-4 a 5-8 muestran los resultados de
las simulaciones. Puede observarse que los diferentes parámetros producen efectos no
acoplados sobre la curva σ1/σ3− ε1, lo que facilita el proceso de selección de parámetros
y la calibración del modelo.
146
1
2
3
4
5
6
0.00 0.01 0.02 0.03 0.04 0.05
50%rD =
100%rD =
3 100 kPaσ =
1ε
1040
840
640
s
s
s
c
c
c
===
1
3
σσ
Figura 5-4: Efecto del parámetro cs que controla la rigidez al corte.
1
2
3
4
5
6
0.00 0.05 0.10 0.15
50%rD =100%rD =
3 100 kPaσ =
34
31
28
c
c
c
φφφ
= °= °= °
1ε
1
3
σσ
Figura 5-5: Efecto del parámetro φc que controla el estado crítico.
147
1
2
3
4
5
6
0.00 0.01 0.02 0.03 0.04 0.05
50%rD =
100%rD =
3 100 kPaσ =
75
55
35
r
r
r
p
p
p
===
1ε
1
3
σσ
Figura 5-6: Efecto del parámetro pr de dureza de grano que controla la diferencia φ−φc.
1
2
3
4
5
6
0.00 0.01 0.02 0.03 0.04 0.05
50%rD =
100%rD =
3 100 kPaσ =
0.50
0.40
0.30
d
d
d
c
c
c
===
1ε
1
3
σσ
Figura 5-7: Efecto del parámetro cd que controla la curva σ − ε para Dr = 100%.
148
1
2
3
4
5
6
0.00 0.01 0.02 0.03 0.04 0.05
50%rD =
100%rD =
3 100 kPaσ =
0.35
0.25
0.15
l
l
l
c
c
c
===
1ε
1
3
σσ
Figura 5-8: Efecto del parámetro cl que controla la curva σ − ε para Dr = 50%.
5.5.3 Análisis de parámetros de estado
Ensayos triaxiales drenados
Se simularon ensayos triaxiales drenados convencionales a presión de confinamiento
constante y para un punto de integración. Se emplearon cuatro densidades relativas
y tres presiones de confinamiento para mostrar la respuesta del modelo para diferentes
parámetros de estado.
La figura 5-9 muestra la influencia de la densidad relativa sobre las curvas σ1/σ3− ε1y εv − ε1. Puede observarse que el modelo reproduce el comportamiento de las
arenas en compresión drenada con suficiente precisión para aplicaciones prácticas.
Como las simulaciones se efectuaron en un punto de integración, la respuesta de la
rama descendente de las curvas σ1/σ3 − ε1 no puede ser comparada con resultados
experimentales.
La figura 5-10 muestra la influencia de la presión de confinamiento sobre la respuesta
del modelo.
149
1
2
3
4
5
6
0.00 0.05 0.10 0.15
25%rD =50%rD =
75%rD =
100%rD =
3 100 kPaσ =
1ε
1
3
σσ
-0.10
-0.08
-0.06
-0.04
-0.02
0.00
0.02
0.00 0.05 0.10 0.15
25%rD =
75%rD =
100%rD =
3 100 kPaσ =
1ε
vε
50%rD =
Figura 5-9: Efecto de la densidad relativa sobre la respuesta del modelo. Ensayostriaxiales drenados.
150
1
2
3
4
5
6
0.00 0.05 0.10 0.151ε
3 3 310 kPa 100 kPa 1000 kPaσ σ σ= = =
50%rD =100%rD =
1
3
σσ
-0.10
-0.08
-0.06
-0.04
-0.02
0.00
0.02
0.04
0.00 0.05 0.10 0.151ε
3 3 310 kPa 100 kPa 1000 kPaσ σ σ= = =vε
50%rD =100%rD =
Figura 5-10: Efecto de la presión de confinamiento sobre la respuesta del modelo. Ensayostriaxiales drenados.
151
Tnsayos triaxiales no drenados
Se simularon ensayos triaxiales no drenados convencionales para un punto de integración.
Se emplearon cinco densidades relativas y tres presiones de confinamiento.
La figura 5-11 muestra la influencia de la densidad relativa sobre la respuesta del
modelo. En la figura no se aprecia la curva correspondiente a Dr = 15% porque el
modelo informa que el suelo licúa con esa densidad relativa . En la parte superior de la
figura se muestra el diagrama q − ε1 completo. Se observa que todas las simulaciones
alcanzan un estado de deformación constante cuando la presión media es suficientemente
alta como para impedir la dilatancia. En la parte media de la figura se muestra una vista
ampliada del mismo diagrama, donde puede verse que la muestra con Dr = 15% alcanza
un pico de resistencia reducido y luego decae hasta q = 0 kPa. En la simulación con
Dr = 25% se aprecia el desarrollo del estado de cuasi-deformación constante, mientras
que las muestras más densas presentan comportamiento dilatante en toda la curva. El
efecto de la densidad relativa sobre la tendencia del material a contraer o dilatar se
aprecia más claramente en el diagrama p − q que se muestra en la parte inferior de la
figura.
La realidad física es que aún con densidades muy bajas como Dr = 15% el material
exhibe una resistencia al corte residual no nula [21][22][58][59][115][154]. El modelo no
captura este comportamiento porque se emplea una presión mínima en la evaluación de la
función ψ [p, e0], necesaria para evitar la sobre-estimación del ángulo de fricción interna
máximo.
La figura 5-12 muestra la influencia de la presión de confinamiento sobre la respuesta
del modelo. Se aprecia que la influencia de la presión de confinamiento inicial se limita
a la porción inicial de las curvas q − ε1 hasta el estado de cuasi-deformación constante
(ς < 1). A partir del desarrollo de la plasticidad en falla, todas las curvas q−ε1 coincideny predicen correctamente la evolución hasta el estado de deformación constante.
152
0
2000
4000
6000
8000
10000
12000
0.00 0.05 0.10 0.151ε
100%
75%
50%
25%
15%
r
r
r
r
r
D
D
D
D
D
=====
[ ]kPa
q
0
50
100
150
200
0.00 0.01 0.02 0.031ε
100%
75%
50%
25%
15%
r
r
r
r
r
D
D
D
D
D
=====
[ ]kPa
q
0
50
100
150
200
0 50 100 150 200 250 300
100%
75%
50%
25%
15%
r
r
r
r
r
D
D
D
D
D
=====
[ ]kPap
[ ]kPa
q
Figura 5-11: Simulación de ensayos de compresión triaxial no drenada para diferentesdensidades relativas.
153
0
50
100
150
200
0 50 100 150 200 250 300[ ]kPap
[ ]kPa
q
15%rD =
25%rD =50%rD =
0
50
100
150
200
0.00 0.01 0.02 0.031ε
[ ]kPa
q
15%rD =
25%rD =
50%rD =
0
500
1000
1500
2000
2500
3000
0.00 0.02 0.04 0.06 0.08 0.101ε
[ ]kPa
q
15%rD =25%rD =
50%rD =
Figura 5-12: Simulación de ensayos de compresión triaxial no drenada para diferentespresiones de confinamiento y densidades relativas.
154
Se concluye que el modelo también reproduce el comportamiento de las arenas
en compresión no drenada con suficiente precisión para aplicaciones prácticas. Es
interesante destacar que el modelo captura correctamente la dependencia del ángulo de
transformación de fase respecto de la densidad relativa, la licuefacción estática de arenas
muy sueltas, el estado de cuasi-deformación constante de las arenas sueltas y el estado
de deformación constante de todas las arenas [58][59][115][154].
Ensayos de deformación plana drenados
Las ecs. (4.58) a (4.59) tienen un punto común en el caso de la deformación plana, donde
se define la ec. (4.60), para el que β presenta una discontinuidad en su derivada. En la
ec. (4.69) se presentó una formulación alternativa, que emplea los factores de dilatancia
en compresión y extensión triaxial y calcula una solución interpolada para otros ángulos
de Lode. En la figura 5-13 se comparan las dos fórmulas mediante la simulación de
ensayos de deformación plana drenados para cuatro densidades relativas y una presión
de confinamiento σ3 = 100 kPa. En la figura, β es el factor de dilatancia de las ecs.
(4.58) y (4.59), mientras que β según ec. (4.69).
Ambas predicciones producen idénticas curvas σ1/σ3 − ε1, hasta el punto que sólo
puede apreciarse la existencia de dos formulaciones como un engrosamiento en la línea
de la curva σ− ε para Dr = 100%. En cambio, las dos formulaciones producen diferentes
líneas εv − ε1, como puede apreciarse en la mitad inferior de la figura 5-13. La diferencia
entre los resultados de las ecs. (4.58) y (4.59) respecto de la ec. (4.69) es reducida e
ignorable para todos los fines prácticos. Es más, no existe evidencia contundente que
permita establecer cual de las dos formulaciones reproduce mejor el comportamiento real
de las arenas, porque la diferencia entre ambas es, aproximadamente, de igual orden
que el error de las mediciones. Por lo tanto, se concluye que la ec. (4.69) es adecuada
para la aplicación práctica del modelo. Es interesante destacar que el modelo predice
correctamente que (σ1/σ3)f en deformación plana es aproximadamente un 20% superior
a lo que se obtiene en compresión triaxial (figura 5-9).
155
1
2
3
4
5
6
7
0.00 0.02 0.04 0.06 0.08 0.10
25%rD =50%rD =
75%rD =
100%rD =
3 100 kPaσ =
1ε
1
3
σσ
-0.06
-0.04
-0.02
0.00
0.02
0.00 0.02 0.04 0.06 0.08 0.10
25%rD =
75%rD =
100%rD =
3 100 kPaσ =
1ε
vε
50%rD =
ββ
Figura 5-13: Simulación de ensayos de deformación plana drenados. β según ecs. (4.58)y (4.59); β según ec. (4.69).
156
5.5.4 Calibración para la arena Puelchense
En esta sección se presenta la calibración del modelo para la serie de ensayos
experimentales presentada en las figuras 4-27 y 4-28. Para esta arena se determinó
emin = 0.58 y emax = 0.89 mediante ensayos convencionales de laboratorio. En la figura
4-26 se mostró la determinación de otros dos parámetros: φc = 30◦, pr = 50. Como no
se dispone de elementos para medirlo, se estima a partir de la bibliografía que cs = 740
[49][59][66].Por lo tanto, los parámetros a ajustar son cd y cl.
Únicamente se presentan las simulaciones de los ensayos ejecutados con presión de
confinamiento σ3 100 kPa, 300 kPa y 600 kPa para que los gráficos resulten legibles. Los
parámetros de mejor ajuste fueron cd = 0.35 y cl = 0.15. Aunque estos valores son
algo bajos, debe tenerse en cuenta que las mediciones de desplazamientos de los ensayos
experimentales fueron externas, por lo que las curvas σ1/σ3 − ε1 incluyen los errores
introducidos por los contactos cabezal-muestra. En las figuras 5-14 y 5-15 se muestran
los resultados numéricos vs. los experimentales. Puede observarse que el ajuste del
modelo es razonable para Dr = 97%, mientras que para Dr = 35% el modelo informa
más dilatancia que los ensayos.
5.5.5 Calibración para la arena Toyoura
El modelo que se presenta en este capítulo está expresado en términos de presiones
efectivas, desarrollado alrededor de la teoría de estado crítico y tiene siete parámetros
materiales. Por lo tanto, la reproducción de resultados experimentales de ensayos
triaxiales drenados es relativamente trivial y no permite apreciar la verdadera capacidad
predictiva de la formulación.
Distinto es el caso de la compresión triaxial no drenada de arenas. En este ensayo, el
material es saturado con agua, comprimido isotrópicamente hasta una cierta presión de
confinamiento, y luego sometido a un esfuerzo q = σ1 − σ3 monotónicamente crecientemientras se mantiene constante el volumen de la muestra. Por lo tanto, la tendencia del
material a contraer se traduce en un aumento de la presión del agua intersticial y en
157
1.00
2.00
3.00
4.00
5.00
6.00
0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10
600 300 100
1
3
σσ
1ε
Presión de confinamiento
-5.0%
-4.0%
-3.0%
-2.0%
-1.0%
0.0%
1.0%
0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.101ε
vε
Figura 5-14: Ensayos triaxiales drenados de arena del Puelchense con Dr ≈ 97%.Comparación entre mediciones experimentales y resultados de la simulación numérica.
158
1.00
2.00
3.00
4.00
0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10
600 300 100
1ε
1
3
σσ
Presión de confinamiento
-2.0%
-1.0%
0.0%
1.0%
2.0%
3.0%
4.0%
0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.101ε
vε
Figura 5-15: Ensayos triaxiales drenados de arena del Puelchense con Dr ≈ 35%.Comparación entre mediciones experimentales y resultados de la simulación numérica.
159
una consiguiente disminución de la presión efectiva media. En cambio, la tendencia del
material a dilatar se traduce en una disminución de la presión intersticial de agua y en el
consiguiente aumento de la presión efectiva. Si las partículas fueran rígidas, el aumento
de la presión efectiva sólo terminaría con la cavitación del agua intersticial. Como las
partículas se rompen, el aumento de la presión efectiva termina cuando se alcanza un nivel
de presión que suprime la dilatancia o, en otras palabras, cuando la muestra alcanza el
estado de deformación constante. Por lo tanto, en el ensayo de compresión no drenada se
puede probar el funcionamiento conjunto de las funciones de estado para rigidez elástica,
tensión-dilatancia, ángulo de fricción interna y función de endurecimiento y, sobre todo,
la dependencia de todas estas funciones de estado respecto de la presión y la densidad
relativa.
En las figuras 5-16 a 5-18 se presenta la simulación numérica de ensayos de compresión
no drenada de la arena Toyoura, para cuatro presiones de confinamiento y tres densidades
relativas. Los datos experimentales fueron obtenidos de Verdugo & Ishihara [154]. Tres de
los parámetros de entrada del modelo están reportados en [154]: emin = 0.60, emax = 0.98
y φc = 31◦. Los otros cuatro parámetros se seleccionaron para obtener el mejor ajuste
a los datos experimentales. Quedó cs = 350, pr = 50, cd = 0.93 y cl = 0.56. Para cada
densidad relativa, los ensayos fueron ejecutados a una misma relación de vacíos ec luego
de la consolidación hasta la presión pc. Para cada simulación se calculó el parámetro e0
correspondiente mediante la expresión
e0 = (1 + ec) exp [pc/K]− 1 (5.93)
Los ensayos triaxiales de Verdugo e Ishihara fueron ejecutados con medición de la
deformación axial [154]. Esta técnica experimental produce una fuerte subestimación
de la rigidez inicialde los. Para lograr un buen ajuste debió adoptarse cs = 350.
Consecuentemente, los parámetros cd y cl que miden el cociente E50/E debieron
aumentarse. Si se duplicara cs a cs = 700 y se redujera a la mitad cd = 0.46 y cl = 0.28
se obtendría el resultado que se muestra en la figura 5-19.
160
Figura 5-16: Compresión no drenada de Arena Toyoura, Dr = 18.5%. Comparaciónnumérica vs experimental. Data de [154].
161
Figura 5-17: Compresión no drenada de Arena Toyoura, Dr = 38%. Comparaciónnumérica vs experimental. Data de [154].
162
Figura 5-18: Compresión no drenada de Arena Toyoura, Dr = 63.5%. Comparaciónnumérica vs experimental. Data de [154].
163
Figura 5-19: Compresión no drenada de Arena Toyoura, Dr = 63.5%. Comparaciónnumérica vs experimental. Parámetros realistas. Data de [154].
Nunca se logra un ajuste perfecto de las curvas tensión-deformación. Esto es debido
a la falta de flexibilidad del modelo inducida por el reducido número de parámetros
materiales. Esta es una limitación aceptada para un modelo desarrollado para su
aplicación en la práctica profesional, donde es muy apreciado que el modelo tenga pocos
parámetros de calibración simple y donde los ensayos experimentales de alta calidad como
los usados para esta comparación [154] son muy escasos.
5.5.6 Análisis de las variables de estado
El modelo tiene tres variables de estado: σ, e0 y ς. Mientras que la primera está presente
en todos los modelos de la mecánica del continuo, las otras dos son específicas de este
modelo y requieren el siguiente análisis.
La variable de estado e0 es una medida del estado de densidad del material, y es
la única que registra los cambios de volumen asociados a la deformación plástica. Su
principal ventaja es que su determinación es muy simple; la segunda, que e0 es una
medida directa de la deformación finita del material.
164
En otros modelos para geomateriales se emplean medidas de la deformación plástica
acumulada como variables de estado. Sin embargo, la deformación plástica acumulada
no es una buena variable de estado para arenas. Aunque hay muchos argumentos
que permitirían respaldar esta afirmación, no es necesario listarlos, sino que basta con
reflexionar acerca de estas dos preguntas: ¿Cual es la deformación plástica acumulada de
la arena de un reloj de arena que fue invertido una vez? ¿y luego de 100 inversiones? e0
es el mismo en ambos casos y el estado del material también.
La variable de estado ς es indirecta. Para medirla es necesario calcular µf [p, e0] y
aplicar la expresión
ς =
√3
µf
r2 − 3J3r1− r2/2 + J3r
(5.94)
o recurrir a la mecánica de suelos clásica. Queda
ς =tan [φmob]
tan [φ](5.95)
No existe un estado físico en el material que pueda ser asociado a ς, porque ς es
una variable de estado escalar - necesariamente asociada a una eventual característica
isotrópica del material - y las arenas sometidas a trayectorias generales de carga no
presentan un comportamiento isotrópico.
En carga monotónica de corte, es plausible que haya una relación entre la proporción
de los contactos que sufren deslizamiento y ς. De todas maneras, ς es una variable de
estado que, por definición, alcanza su valor de saturación antes del estado crítico y por
lo tanto no controla los aspectos básicos del comportamiento del material para grandes
deformaciones. Otro aspecto de gran importancia es que las variables de estado estén
asociadas a procesos claramente identificables y, en la medida de lo posible, desacoplados.
El grado de acoplamiento entre variables de estado se puede estudiar de manera
gráfica en mapas de evolución relativa. En la figura 5-20 se muestran los mapas ς − e0 yσ1/σ3− e0 para los ensayos triaxiales drenados que se presentaron en la figura 5-9. En la
parte superior puede apreciarse que las trayectorias son prácticamente ortogonales entre
165
sí, que ς controla el comportamiento del material hasta la falla y que e0 lo controla desde
la falla y hacia el estado crítico. En la parte inferior de la figura puede apreciarse que
una vez que se alcanza la falla y el material comienza a dilatar, describe una única curva.
Cualquier punto en esa curva límite corresponde a un material que tiene un estado único,
independiente de la relación de vacíos inicial del ensayo.
En la figura 5-21 se muestran los mapas ς − e0 y q− e0 para los ensayos triaxiales nodrenados que se presentaron en la figura 5-11. Puede apreciarse que en este caso e0 no
cambia significativamente en la falla, como corresponde a un análisis no drenado, en el
que la deformación volumétrica está limitada por la rigidez volumétrica del agua, y es
prácticamente nula para todos los ensayos corrientes de laboratorio.
166
0.00
0.20
0.40
0.60
0.80
1.00
0.50 0.60 0.70 0.80 0.90 1.00
ς
25%rD =50%rD =75%rD =100%rD =
0e
3 100 kPaσ =
1.00
2.00
3.00
4.00
5.00
0.50 0.60 0.70 0.80 0.90 1.00
25%rD =50%rD =75%rD =100%rD =
1
3
σσ
0e
3 100 kPaσ =
Figura 5-20: Mapas de variables de estado. Ensayo triaxial de compresión a presión deconfinamiento constante.
167
0.00
0.20
0.40
0.60
0.80
1.00
0.50 0.60 0.70 0.80 0.90 1.00
25%rD =50%rD =75%rD =100%rD =
ς
0e
0
5000
10000
15000
0.50 0.60 0.70 0.80 0.90 1.00
[ ]kPa
q
0e
25%rD =50%rD =75%rD =100%rD =
Figura 5-21: Mapas de variables de estado. Ensayo triaxial de compresión no drenada.
168
5.5.7 Análisis del tamaño del paso
Tanto el operador de integración como su algoritmo numérico son implícitos en todos sus
niveles, por lo que la estabilidad numérica está asegurada. La consistencia del algoritmo
y la dependencia de la respuesta del modelo respecto del tamaño del paso se muestran
en las figuras 5-22 y 5-23 para el ensayo triaxial de compresión no drenada, en la que se
muestra la predicción del modelo para una deformación total ‖∆εd‖ = 5% en 5, 10 y 100
pasos.
0
100
200
300
400
500
0 100 200 300 400 500 600 700 800 900 1000[ ]kPap
[ ]kPa
q
0
50
100
150
200
250
300
0.00 0.01 0.02 0.03 0.041ε
[ ]kPa
q
25%rD =
Figura 5-22: Influencia del tamaño del paso sobre la respuesta del modelo. ‖εd‖ = 5%en 5 pasos, 10 pasos y 100 pasos. Dr = 25%.
169
0
1000
2000
3000
4000
5000
0 2000 4000 6000 8000[ ]kPap
[ ]kPa
q
0
1000
2000
3000
4000
5000
0.00 0.01 0.02 0.03 0.041ε
[ ]kPa
q
75%rD =
Figura 5-23: Influencia del tamaño del paso sobre la respuesta del modelo. ‖εd‖ = 5%en 5 pasos, 10 pasos y 100 pasos. Dr = 75%.
170
5.5.8 Simulación de un ensayo de compresión plana
En esta sección se repite la simulación de un ensayo de compresión plana drenado que se
presentó en la sección 4.5.3, con los mismos parámetros materiales, variables de estado
y malla. La figura 5-24 muestra la curva tensión-deformación que puede calcularse con
información del contorno de la muestra. Si se asume una deformación uniforme, el valor
terminal de P/A se puede estimar con las fórmulas corrientes de la mecánica de suelos,
y asumiendo φps = 1.15φc [69]. Queda
(P
A0
)
c
≈ 1
1−∆H/H tan2
[π
4+1.15 · φc
2
]σ3 = 413 kPa (5.96)
0
100
200
300
400
500
600
0.00 0.02 0.04 0.06 0.08 0.10H H∆
[ ]0
kPa
P A
Figura 5-24: Ensayo de deformación plana. Curva tensión-deformación observable desdelos contornos de la muestra.
En las figuras que siguen se muestra la evolución de algunas variables para
acortamientos de la muestra de 2%, 4%, 6% y 8% de la altura inicial. En la figura 5-25 se
muestra la malla deformada. En la figura 5-26 se muestra el contorno de desplazamientos
verticales como porcentual de la altura inicial. En la figura 5-27 se muestra la evolución
de ‖ε‖ . En la figura 5-28 se muestra la evolución de e0. Por último, en la figura 5-29 se
muestra la evolución de φ.
171
Figura 5-25: Ensayo de deformación plana. Malla deformada. ∆H/H = 2%, 4%, 6% y8%.
Figura 5-26: Ensayo de deformación plana. Contorno de desplazamientos. ∆H/H = 2%,4%, 6% y 8%.
172
Figura 5-27: Ensayo de deformación plana. Contorno de ‖ε‖. ∆H/H = 2%, 4%, 6% y8%. Azul: ‖ε‖ = 0. Colorado: ‖ε‖ = 100%.
Figura 5-28: Ensayo de deformación plana. Contorno de e0. ∆H/H = 2%, 4%, 6% y8%. Azul: Dr = 0%. Colorado: Dr = 100%.
173
Figura 5-29: Ensayo de deformación plana. Contorno de φ. ∆H/H = 2%, 4%, 6% y 8%.Azul: φ = φc = 31
◦. Colorado:φ = 41◦.
5.5.9 Simulación de un problema de valores de contorno
complejo
En esta sección se presenta la simulación de un problema de valores de contorno complejo,
típico de la ingeniería geotécnica práctica. Sólo se informan algunas salidas gráficas
representativas, porque el único objeto de esta sección es mostrar la robustez de la
implementación numérica modelo y su capacidad para resolver problemas prácticos. Se
simula la construcción de una excavación cercana a un túnel existente. Los parámetros
materiales son los de la arena de la Formación Puelchense. La relación de vacíos de todo
el depósito es e0 = 0.75. Se asumió un estado tensional inicial con σ3 = σ1/2 y σ1 = γz,
donde γ es el peso unitario del material y z es la profundidad de cada punto. Las etapas
del proyecto que se simula son:
1. Construcción del túnel
(a) Se excava el terreno y se permite una relajación parcial del estado tensional
174
(b) Se coloca el revestimiento del túnel
2. Ejecución de la excavación
(a) Colocación de la tablestaca
(b) Primera etapa de excavación
(c) Colocación de un puntal
(d) Segunda etapa de excavación
En la figura 5-30 se muestra la malla y las etapas del modelo. En la figura 5-
31 se muestra la trayectoria de tensiones principales, el diagrama de empujes sobre la
tablestaca, el diagrama de empujes pasivos del lado excavación, los desplazamientos de la
tablestaca y su diagrama de momentos flexores. También se muestran los desplazamientos
del túnel, su diagrama de esfuerzos normales y su diagrama de momentos flexores. Como
el gráfico se presenta con fines ilustrativos, no se indican los valores máximos de las
diferentes variables.
En la figura 5-32 se muestra el contorno de deformaciones. Puede observarse
que la deformación máxima es del orden de ‖ε‖ = 2.5% y ocurre únicamente en
el pie de la excavación, mientras que más del 80% de la malla tiene deformaciones
menores a 0.5%. Este resultado, típico de los problemas rutinarios de geotecnia, explica
porqué es importante mejorar los modelos constitutivos disponibles para predecir el
comportamiento en el rango de baja deformación.
Por último, en la figura 5-33 se muestra el contorno de ς. Puede observarse que el pie
de la excavación y el coronamiento alcanzaron localmente la falla. También se aprecia que
algunos puntos en la superficie tienen ς = 1. Esto se debe a que, durante la construcción
del túnel, esos puntos tuvieron un incremento de tensiones de compresión horizontal.
Como está en la superficie, el material tiene rigidez y resistencia prácticamente nulas.
175
Figura 5-30: Malla de elementos finitos y etapas de construcción. Las mallas estándeformadas con un factor de amplificación 50.
176
Figura 5-31: Estado final del modelo. Se muestran la trayectoria de tensiones, losempujes del terreno lado tierra y lado excavación, los desplazamientos de la tablestaca,sus momentos flexores y los desplazamientos, esfuerzos normales y momentos flexores deltúnel.
177
Figura 5-32: Contorno de deformaciones ‖ε‖ . La deformación máxima es ‖ε‖ = 2.5%.La escala está multiplicada por 105.
Figura 5-33: Movilización de la resistencia al corte ς para el fin de construcción.
178
Capítulo 6
Consideraciones finales y
oportunidades de investigación
futura
6.1 Resumen del comportamiento de las arenas
El comportamiento macroscópico de las arenas es la resultante de la interacción entre
granos aislados que, sometidos a cargas, rotan, se desplazan y se rompen. El estudio de
este comportamiento con las herramientas de la mecánica del continuo es, por lo tanto,
necesariamente complejo y de alcance limitado.
Hay algunos aspectos básicos del comportamiento de las arenas que están bien
establecidos, mientras que otros son aún objeto de investigación. Entre los primeros
están:
1. Las arenas sometidas a procesos de deformación por corte continuada alcanzan un
estado terminal, denominado estado crítico, que es independiente del estado inicial
del material, y que sólo depende de la presión media que actúa en el estado crítico.
2. En el estado crítico, tanto la relación de vacíos como el tensor de tensiones alcanzan
179
valores estacionarios.
3. Existe una interrelación bien establecida entre la resistencia máxima y la tendencia
del material a cambiar de volumen.
4. Las arenas son materiales disipativos - sin dominio elástico - para cualquier escala
de deformación.
5. La resistencia al corte depende de la tensión y la densidad.
6. La rigidez instantánea depende de la tensión, la densidad, y de la trayectoria de
carga.
Muchos otros aspectos del comportamiento de las arenas son objeto de investigación
actual. Entre los que se requieren para desarrollar nuevos modelos constitutivos para
arenas están:
1. El desarrollo de técnicas de medición experimental de nuevas variables de estado que
caractericen la topología del material y su evolución. El número de coordinación, los
cambios de granulometría, los cambios de esfericidad y redondez de las partículas
y la densidad direccional de contactos son candidatos posibles, para los que hay
desarrollos continuos en sus técnicas de medición.
2. La comprensión de la interrelación entre esas variables de estado y los efectos de
memoria que el material exhibe cuando se lo somete a trayectorias complejas de
carga. La simulación de estos mecanismos se hace hoy con modelos que tienen un
número de variables internas que no pueden medirse experimentalmente.
3. La caracterización de los efectos de velocidad de carga sobre la respuesta del
material.
180
6.2 Resumen de los modelos constitutivos para
arenas
Existen varias familias de modelos constitutivos para arenas, desde el más simple -
el modelo de elasticidad lineal y plasticidad de Mohr-Coulomb - hasta modelos muy
complejos que tienen múltiples superficies de fluencia móviles, ecuaciones de evolución
generalizadas y numerosas variables de estado tensoriales.
La industria resuelve la enorme mayoría de sus problemas con el modelo de Mohr-
Coulomb, aún a sabiendas de que tiene baja capacidad predictiva. Los problemas que
no se resuelven con el modelo de Mohr-Coulomb se resuelven con el modelo hiperbólico
o modelos similares, que tienen mejor capacidad predictiva pero que, como el primero,
no reproducen el comportamiento del material ni aún en los aspectos que están mejor
comprendidos y universalmente aceptados. En particular, en estos modelos:
1. El ángulo de fricción interna es un parámetro material, no una función de estado.
Por lo tanto, es constante a lo largo del proceso de deformación, aunque cambie la
presión y la relación de vacíos.
2. El ángulo de fricción interna crítico no es un parámetro de material. Por lo tanto,
los modelos no simulan la evolución desde el estado inicial al estado crítico.
Por estas y otras razones, estos modelos arrastran una limitación severa que es común
en geotecnia aplicada pero que es poco comprensible para analistas ajenos a este campo:
Los parámetros de entrada dependen del problema que se estudia.
En el otro extremo están los modelos constitutivos para arenas que ofrece la academia.
Todos los modelos simulan el comportamiento básico del material, y muchos tienen
notable capacidad predictiva aún para trayectorias muy complejas de carga. Sin embargo,
muy pocos están implementados en códigos disponibles para su empleo en la industria,
hasta el punto que casi no se reportan aplicaciones de modelos constitutivos avanzados
para construcciones reales.
181
6.3 Características del modelo constitutivo
desarrollado en esta tesis
El modelo constitutivo desarrollado en esta tesis es una contribución que llena
parcialmente el espacio que hay entre lo que produce la academia y lo que emplea la
industria. Está orientado al análisis de problemas de ingeniería geotécnica en los que el
material está sometido a procesos de carga por corte dominante, desde las vibraciones
producidas por fundaciones y hasta el análisis de procesos de falla. Los procesos que
pueden ser simulados son los mismos que se simulan con los modelos de Mohr-Coulomb
e hiperbólico y constituyen, por lejos, el grupo de problemas que se presenta con más
frecuencia en la geotecnia práctica.
Las principales cualidades de este modelo constitutivo son:
1. Formulación
(a) Tiene siete parámetros de entrada, independientes del problema que se estudia
y determinables mediante ensayos triaxiales convencionales. Los parámetros
no superponen su efecto sobre la respuesta del modelo para los ensayos de
calibración.
(b) Tiene tres variables de estado: la tensión, la relación de vacíos para tensión
cero y la movilización de resistencia al corte, que pueden ser determinadas
con mediciones convencionales de laboratorio. Como lo requiere la teoría de
estado crítico, las variables de estado alcanzan un valor estacionario para el
estado crítico.
(c) El ángulo de fricción interna crítico es un parámetro de entrada (queda fijo
durante el proceso de deformación).
(d) El ángulo de fricción interna máximo es una función de estado; varía durante
el proceso de deformación y alcanza el valor crítico para deformaciones
suficientemente grandes.
182
(e) El ángulo de fricción interna máximo controla la evolución de una superficie de
fluencia de tres invariantes que reproduce la resistencia al corte para cualquier
presión, densidad y trayectoria de carga, tanto drenada como no drenada.
(f) El ángulo de fricción interna a volumen constante controla el término de
dilatancia que se calcula con la teoría de tensión-dilatancia de Rowe en su
forma más general.
(g) La formulación hiperelástica reproduce la rigidez máxima del material, lo que
permite la simulación de procesos de baja deformación, como propagación
de ondas y análisis de vibraciones elásticas. El módulo de rigidez al corte
elástico utilizado tiene procedimientos normalizados de medición, tanto en
campo como en laboratorio.
(h) El modelo reproduce la degradación de la rigidez aparente con la deformación
desde la rigidez elástica hasta la falla y luego hasta el estado crítico.
2. Respecto de la implementación
(a) Está completamente implementado en un programa de elementos finitos
general y probado para los problemas rutinarios de geotecnia.
(b) Tiene un algoritmo encestado consistente que resuelve una incógnita escalar
por vez.
(c) La implementación es suficientemente robusta como para permitir el estudio
de problemas en los que ocurre falla localizada y otros en los que el material
alcanza el estado σ = 0 kPa.
(d) El algoritmo entrega una respuesta para cualquier juego de parámetros de
entrada y variables de estado y cualquier incremento de deformación.
183
6.4 Caminos de investigación futura
El modelo constitutivo presentado en esta tesis resuelve satisfactoriamente los problemas
para los que fue diseñado. Como todos los modelos que producen curvas tensión-
deformación con ablandamiento, requiere técnicas de regularización de deformaciones
localizadas. La implementación de una de estas técnicas es el primer paso en la evolución
de este trabajo. El segundo es la implementación del modelo en otros códigos generales
de elementos finitos.
Este modelo no resuelve los siguientes problemas:
1. Carga cíclica en régimen irreversible, como la que ocurre en los depósitos de arenas
sometidos a terremotos.
2. Compresión plástica para alta presión, como la que ocurre en la punta de pilotes
cargados hasta la falla y en los niveles inferiores de las presas de materiales sueltos
altas.
Sfriso [137] presentó un modelo para corte y compresión plástica que comparte
elementos con este trabajo. Sfriso [136] también presentó un modelo constitutivo
orientado a la simulación de los procesos de carga cíclica, basado en plasticidad
generalizada. Estos modelos no son compatibles entre sí, porque emplean diferentes
variables de estado.
Por lo tanto, la unificación de los elementos presentados en [136] y [137] para el
desarrollo e implementación de un modelo para compresión plástica y carga cíclica son
el tercer paso en el camino de investigación futura en el campo de la modelización
constitutiva para geomateriales.
184
Capítulo 7
Apéndice: Algebra Tensorial
Se presentan los elementos básicos de álgebra tensorial que se emplearon en las
deducciones de fórmulas y en la implementación del modelo constitutivo. Una parte
importante de las definiciones y fórmulas fueron obtenidas o derivadas de [42].
7.A Tensores unitarios
El tensor unitario de segundo orden se define como
1 = δij eiej δij = ei · ej (7.1)
El tensor unitario de cuarto orden simétrico se define como
I =1
2(δik δjl + δil δjk) eiejekel (7.2)
y puede separarse en sus proyecciones esféricas y deviatóricas I = Iesf + Idev donde
Iesf =1
311 Idev = I− Iesf (7.3)
185
7.B Invariantes de tensores de segundo orden
Para un tensor simétrico de segundo orden se definen los invariantes
I1 = T : 1 = Tii (7.4)
I2 =1
2
(I21 −T : T
)(7.5)
I3 =1
6
(2 (T ·T) : T− 3I1 (T : T) + I31
)= det [T] (7.6)
J1 = tr(Td)= 0 (7.7)
J2 =1
2Td : Td (7.8)
J3 =1
3
(Td ·Td
): Td (7.9)
Estos invariantes están relacionados mediante
J2 =1
3
(I21 − 3I2
)(7.10)
J3 =1
27
(2I31 − 9I1I2 + 27I3
)(7.11)
I2 =1
3
(I21 − 3J2
)(7.12)
I3 =1
27
(I31 − 9I1J2 + 27J3
)(7.13)
Otras operaciones frecuentes entre invariantes son
T ·T=Td ·Td +I2191+
2
3I1T
d (7.14)
(T ·T) : 1 = T : T (7.15)
186
T : T = I21 − 2I2 =1
3I21 + 2J2 (7.16)
7.C Derivadas de invariantes
Las siguientes son las derivadas primeras de los principales invariantes de tensores
simétricos de segundo orden
I1,T = 1 (7.17)
I2,T = I11−T (7.18)
I3,T = T ·T+ I21− I1T (7.19)
J2,T = Td (7.20)
J2,Td = Td (7.21)
J3,T = T ·T+(I219− 23J2
)1− 2
3I1T = T
d ·Td − 23J21 (7.22)
J3,Td = Td ·Td − 23J21 = T
d ·Td − 13
(Td : Td
)1 (7.23)
Otras derivadas de interés son
(T : T),T = 2T (7.24)
Td,T = I
dev (7.25)
187
(T ·T),T =∂ (timtmj)
∂tkl= (δiktlj + δjltik) eiejekel = (δikslj + δjlsik) eiejekel + 2pI
=
2t1 0 0 2t4 0 2t6
0 2t2 0 2t4 2t5 0
0 0 2t3 0 2t5 2t6
t4 t4 0 t1 + t2 t6 t5
0 t5 t5 t6 t2 + t3 t4
t6 0 t6 t5 t4 t1 + t3
=
2s1 0 0 2s4 0 2s6
0 2s2 0 2s4 2s5 0
0 0 2s3 0 2s5 2s6
s4 s4 0 −s3 s6 s5
0 s5 s5 s6 −s1 s4
s6 0 s6 s5 s4 −s2
+ 2pI
‖T‖,T =∂√T : T
∂T=
T√T : T
= T (7.26)
(T
‖T‖
)
,T
=∂T∂T‖T‖ −T∂‖T‖
∂T
‖T‖2=I− TT‖T‖ (7.27)
‖A‖,T =∂√A : A
∂T=A : A,T
‖A‖ (7.28)
(A
‖A‖
)
,T
=A,T
‖A‖ −(A
‖A‖A
‖A‖
):A,T
‖A‖ (7.29)
7.D Medidas de tensión y deformación
Existen diferentes medidas del estado de tensiones y deformaciones que se emplean
corrientemente en modelos para geotecnia. Algunas son
r =s
pr =
√r : r (7.30)
q =√3J2 =
√3
2‖s‖ (7.31)
γpoct =
2√3‖εp
d‖ (7.32)
188
Para el estado triaxial de compresión (σ1 > σ2 = σ3), valen
N = σ1/σ3 K = σ3/σ1 q = σ1 − σ3 (7.33)
σ =3p
N + 2{N, 1, 1} = 3p
1 + 2K{1, K,K} (7.34)
τ = ‖s‖ =√6N − 1N + 2
p (7.35)
r =1−K1 + 2K
{2,−1,−1} r =√61−K1 + 2K
(7.36)
J3r = 2
(1−K1 + 2K
)3=
1
3√6r3 (7.37)
r · r = {4, 1, 1}(1−K1 + 2K
)2(7.38)
7.E Derivadas de las medidas de tensión
(σ · σ),σ = (δikslj + δjlsik) eiejekel + 2pI (7.39)
(s · s),σ = (σ · σ),σ −2
3s1−2pI = (δikslj + δjlsik) eiejekel −
2
3s1 (7.40)
τ ,σ =(√s : s)
,σ=s
τ(7.41)
r,σ =∂ (τ/p)
∂σ=1
p
(rr− r31)
(7.42)
r,σ =∂ (s/p)
∂σ=1
p
(Idev − 1
3r1
)(7.43)
J3r,σ =∂ (J3/p
3)
∂σ=1
p
(r · r− r
2
3I− J3r1
)(7.44)
(r · r),σ =(s · sp2
)
,σ
=1
p
((δikrlj + δjlrik) eiejekel −
2
3r1− 2
3r · r1
)(7.45)
189
Referencias
[1] Ahadi A, Krenk S (2003) Implicit integration of plasticity models for granular
materials. Comp Meth Apl Mech Engrg 192:3471-3488
[2] Arulmoni K et al (1992) VELACS Laboratory testing program. ETP 90-0562
[3] Bauer E (1996) Calibration of a comprehensive hypoplastic model for granular
materials. Soil Found 36(1):13-26
[4] Been K, Jefferies G (1985) A state parameter for sands. Geotechnique 35(2):99-112
[5] Been K Jefferies G, Hachey J (1991) The critical state of sands. Geotechnique
41(3):365-381
[6] Benz T, Vermeer A, Schwab R (2008) A small-strain overlay model. Int J Numer
Anal Meth Geomech doi 10.1002/nag.701
[7] Bishop A, Webb D, Skinner A (1965) Triaxial tests on soil at elevated cell pressures.
VI ICSMFE Montreal, 1:170-174
[8] Bishop A (1954) Correspondence on shear characteristics of a saturated silt,
measured in triaxial compression. Geotechnique 36(1):43-45
[9] Bissio J,Tejada C (2006) Caracoles, análisis dinámico de la presa. IV CAP, Posadas
(Argentina), 57-74
[10] Bolognesi A (1989) Materiales granulares en condición no drenada. X CAMSIF
(1):1-31
190
[11] Bolton M (1986) The strength and dilatancy of sands. Geotechnique 36(1):65-78
[12] Bolton M (1987) The strength and dilatancy of sands: Discussion. Geotechnique
37(2):219-226
[13] Borja R, Lee S (1990). Cam-Clay plasticity I: Implicit integration of elastoplastic
constitutive relations. Comp meth appl mech eng 78:49-72
[14] Borja R, Lee S (1991) Cam-Clay plasticity II: Implicit integration of constitutive
equation based on a nonlinear elastic stress predictor. Comp meth appl mech eng
88:225-240
[15] Borja R, Lin C, Montans F (2001) Cam-Clay plasticity IV: Implicit integration of
anisotropic bounding surface model with nonlinear hyperelasticity and ellipsoidal
loading function”. Comp meth appl mech eng 190:3293-3323
[16] Borja R, Sama K, Sanz P (2003) On the numerical integration of three-invariant
elastoplastic constitutive models. Comp meth appl mech eng 192:1227-1258
[17] Briaud J, Gibbens R (1995) Predicted and measured behavior of five spread footings
on sand. ASCE Geot Sp Publ 41
[18] Brinkgreve R (2002) Plaxis users manual. Balkema, Rotterdam
[19] Casagrande A (1936) Characteristics of cohesionless soils affecting the stability of
slopes and earth fills. J. Boston Soc Civil Eng 23(1):13-32
[20] Casagrande A (1975) Liquefaction and cyclic deformation of sands - a critical
review. V PCSMFE, Buenos Aires, 1-25
[21] Castro G (1975) Liquefaction and cyclic mobility of saturated sands. ASCE J Geot
Eng 101(GT6):551-569
[22] Castro G, Poulos S (1977) Factors affecting liquefaction and cyclic mobility. ASCE
J Geot Eng 103(GT6):501-516
191
[23] Collins I, Muhunthan B (2003) On the relationship between stress-dilatancy,
anisotropy, and plastic dissipation for granular materials. Geotechnique 53(7)611-
618
[24] Cooke J (1984) Progress in rockfill dams. ASCE J Geot Eng 110(10):1381-1414
[25] Cooke J, Sherard J (1987) Concrete-face rockfill dam: II. Design. ASCE J Geot
Eng, 113(10):1113-1132
[26] Cooke J (1997) The concrete face rockfill dam. 17 USCOLD San Diego (USA),
117-132
[27] Cornforth D (1964) Some experiments on the influence of strain conditions on the
strength of sand. Geotechnique 14(2):143-167
[28] Cubrinovski M, Ishihara K (1998) Modelling of sand behaviour based on state
concept. Soils Found 38(3):115-127
[29] Cubrinovski M, Ishihara K (1998) State concept and modified elastoplasticity for
sand modelling. Soils Found 38(4):213-225
[30] Cundall P (2000) FLAC 4.0 Users Manual. Itasca Consulting Group, Minessota
[31] Dafalias Y, Popov E (1975) A model of nonlinearly hardening materials for complex
loadings. Acta Mech 21:173—192
[32] Dafalias Y, Manzari M (2004) Simple Plasticity Sand Model Accounting for Fabric
Change Effects. J Eng Mech 130(6):622-634
[33] Dakoulas P, Yu S (1995) Stress dependency of elastic moduli for cross-anisotropic
soils. Geotechnique 45(2):325-332
[34] Davoudzadeh F (1982) Response of sand to three independently controlled principal
stresses. Ph.D. Thesis U College, London
192
[35] De Beer E (1965) Influence of the mean normal stress on the shearing resistance of
sand. VI ICSMFE Montreal, 1:165-169
[36] De Borst R (1986) Nonlinear analysis of frictional materials. Ph. D. Thesis, TU
Delt
[37] De Borst R, Groen A (2000) Computational strategies for standard soil plasticity
models. En: Modeling in Geomechanics, Wiley
[38] De Josseling de Jong G (1976) Rowe’s stress—dilatancy relation based on friction.
Geotechnique 26(3):527—534
[39] Dobry R et al (1982) Prediction of pore water pressure buildup and liquefaction
of sands during earthquakes by the cyclic strain method. Build Sc 138, NBS,
Washington
[40] Duncan J, Chang C (1970) Nonlinear analysis of stress and strain in soils. ASCE
J Soil Mech Found Div 96(SM5):1629-1653
[41] Duncan J, Byrne P, Wong K, Mabry P (1980) Strength, stress-strain and bulk
modulus parameters for finite element analyses of stresses and movements in soil
masses. Report UCB/GT/80-01 Berkeley
[42] Dvorkin E (1993) Notas de cursos de mecánica del continuo y elementos finitos.
FCEFN UBA
[43] Etse G, Willam K (1996) Integration algorithms for concrete plasticity. Eng Comp
13(8):38-65
[44] Fukushima S, Tatsuoka F (1984) Strength and deformation characteristics of
saturated sand at extremely low pressures. Soil Found 24(4):30-48
[45] Gazetas G, Dakoulas P (1992) Seismic analysis and design of rockfill dams: State
of the art. Soil Dyn Earthq Eng 11(1):27-61
193
[46] Gudehus G (1996) A comprehensive constitutive equation for granular materials.
Soil Found, 36, 1, 1-12
[47] Guo P, Stolle D (2004) The extension of Rowe’s stress-dilatancy model to general
stress condition. Soil Found 44(4):1-10
[48] Guo P, Wan R (2007) Rational approach to stress-dilatancy modelling using an
explicit micromechanical formulation. In Bifurcations, Instabilities, Degradation in
Geomechanics, Springer, doi 10.1007/978-3-540-49342-6_10
[49] Hardin B, Richart F (1963) Elastic wave velocities in granular soils. ASCE J Soil
Mech Found Div 89(SM1):33-65
[50] Hardin B, Drnevich V (1972) Shear modulus and damping in soils: Measurement
and parameter effects. ASCE J Soil Mech Found Div 98(SM6):603-624
[51] Hardin B (1985) Crushing of soil particles. ASCE J Geot Eng 111(10):1117-1192
[52] Hardin B (1987) 1-D strain in normally consolidated cohesionless soils. ASCE J
Geot Eng 113(12):1449-1467
[53] Hashiguchi K (1995) On the linear relations of v-ln p and ln v - ln p for isotropic
consolidation of soils. Int J Num Anal Meth Geomech 19:367-376
[54] Heeres O (2001) Modern Strategies for the Numerical Modeling of the Cyclic and
Transient Behavior of Soils. PhD Thesis, TU Delft.
[55] Horne M (1965) The behaviour of an assembly of rotund, rigid, cohessionless
particles I, II. Proc Roy Soc London A(286):68-97
[56] Horne M (1969) The behaviour of an assembly of rotund, rigid, cohessionless
particles III. Proc Roy Soc London A(310):21-34
[57] Ishibashi I, Chen Y, Chen M (1991) Anisotropic behavior of ottawa sand in
comparison with glass spheres. Soil Found 31(1):145-155
194
[58] Ishihara K (1993) Liquefaction and flow failure during earthquakes, The 33rd
Rankine Lecture. Geotechnique 43:349-415
[59] Ishihara K (1996) Soil behaviour in earthquake geotechnics. Clarendon Press,
Oxford
[60] Jambu N (1963) Soil compresibility as determined by oedometer and triaxial tests.
Proc Eur Conf Wiesbaden I:19-25
[61] Jeferies M (1993) Nor-Sand: a simple critical state model for sand. Geotechnique
43(1):91-103
[62] Kim M, Lade P (1988) Single hardening constitutive model for frictional materials
I. Plastic potential function. Comp Geotech 5(4):307-324
[63] Kirkpatrick W (1957) The condition of failure for sands. IV Int Conf Soil Mech
Found Eng London 1:172-178
[64] Ko H, Scott R (1968) Deformation of sand at failure. ASCE J Soil Mech Found
Div 94(SM4):883-898
[65] Ko H, Scott R (1969) Deformation of sand in hydrostatic compression. ASCE J
Soil Mech Found Div 97(SM3):137-156
[66] Kokusho T (1980) Cyclic triaxial test of dynamic soil properties for wide strain
range. Soil Found, 20:45-60
[67] Kondner R, Zelasko J (1963) Void ratio effects on the hyperbolic stress-strain
response of a sand. Lab Shear Testing Soils ASTM, 250-257
[68] Kondner R (1963) Hyperbolic stress-strain response; cohesive soils. ASCE J Soil
Mech Found Div 78(SM1):115-143
[69] Kulhawi F, Mayne P (1990) Manual on estimating soil properties for foundation
design. Report EL-6800/1493-6, EPRI
195
[70] Lade P, Duncan J (1973) Cubical triaxial tests on cohesionless soils. ASCE J Soil
Mech Found Div, 99(SM10):793-811
[71] Lade P, Duncan M (1975) Elastoplastic stress strain theory for cohesionless soils.
ASCE J Geot Eng Div 101(GT10):1037-1053
[72] Lade P (1977) Elastoplastic stress strain theories for cohessionless soils with curved
yield surfaces. Int J Solids Struct 13:1014-1035
[73] Lade P(1978) Prediction of undrained behavior of sand. ASCE J Geot Eng Div
104(GT6):721-735
[74] Lade P (1988) Double hardening constitutive model for soils: parameter
determination and predictions for two sands. Proc. Cleveland Workshop Const
Eq Granular Non-Cohesive Soils 367-382
[75] Lade P, Kim M (1988) Single hardening constitutive model for frictional materials
II. Yield criterion and plastic work contours. Computers and Geotechnics 6(1):13-30
[76] Lade P, Kim M (1988) Single hardening constitutive model for frictional materials
III. Comparison with experimental data. Computers and Geotechnics 6(1):31-48
[77] Lambe W, Whitman R (1969) Soil Mechanics. Wiley, New York
[78] Lee, K(1965) Triaxial compressive strength of saturated sands under seismic loading
conditions. PhD Dissertation, Berkeley
[79] Lee K, Seed H (1967) Drained strength characteristics of sands. ASCE J Soil Mech
Found Div 93(SM6):117-141
[80] Leps T (1970) Review of Shearing Strength of Rockfill. ASCE J Soil Mech Found
Div 96(4):1159-1170
196
[81] Lush A, Weber G, Anand L (1989) An implicit time-integration procedure for a set
of internal variable constitutive equations for isotropic elasto-viscoplasticity. Int J
Plast 5:521-549
[82] Macari E (1989) Behavior of granular materials in a reduced gravity environment
and under low effective stresses. PhD Thesis, U Colorado
[83] Macari E, Weihe S, Arduino P (1997) Implicit integration of elastoplastic
constitutive models for frictional materials with highly nonlinear hardening
functions. Mech coh frict mat 2:1-29
[84] Maeda K, Miura K (1999) Confining stress dependency of mechanical properties of
sands. Soil Found, 39(1):53-67
[85] Maeda K, Miura K (1999) Relative density dependency of mechanical properties
of sands. Soil Found, 39(1):69-79
[86] Makdisi F, Seed B (1978) Simplified procedure for estimating dam and embankment
earthquakes induced deformation. ASCE J Geot Eng 104(7):849-867
[87] Manzari M, Dafalias Y (1997) A two-surface critical plasticity model for sand.
Geotechnique 47(2):255—272
[88] Marsal, R (1967) Large scale testing of rockfill materials. ASCE J Soil Mech Found
Div, 93, SM2, 27-43
[89] Marsal R, Reséndiz Núñez D (1983) Presas de tierra y enrocamiento. Limusa,
México
[90] Matsuoka H, Nakai T (1974) Stress-deformation and strength characteristics of soil
under three different principal stresses. Proc Japan Soc Civil Eng 233:59-70
[91] Matsuoka H, Nakai T (1985) Relationship among Tresca, Mises, Mohr-Coulomb
and Matsuoka-Nakai failure criteria. Soil Found, 25(4):123-128
197
[92] Mitchell J (1993) Fundamentals of soil behavior. 2◦ Ed, Wiley, New York
[93] Miura N(1979) A consideration of the stress - strain relation of a sand under high
pressures. Proc Japan Soc Civ Eng 282(2):127-130
[94] Miura N, Murata H, Yasufuku N (1984) Stress - strain characteristics of sand in a
particle crushing region. Soil Found, 24(1):77-89
[95] Molenkamp F (1988) A simple model for isotropic non-linear elasticity of frictional
materials. Int J Numer Anal Meth Geomech 12(5):467-475
[96] Mohr O (1900) Welche umstände bedingen die elastizitäsgrenze und den bruch
eines materials?. ZeitschrVereines deutschIngenieure, 44:1-12 citado en [151]
[97] Mroz Z, Norris V, Zienkiewicz O (1981) An anisotropic critical state model for soils
subjected to cyclic loading. Geotechnique, 31(4):451-469
[98] Murata H, Miura N, Hyodo M, Yasufuku N (1989) Experimental study on yielding
of sands. En: Mech Granular Mat, ISSMFE TC13, 173-178
[99] Nakai T, Matsuoka H (1986) A generalized elastoplastic constitutive model for clay
in three dimensional stresses. Soil Found, 26(3):91-98
[100] Nakata Y, Hyodo M, Hyde A, Kato Y, Murata H (2001) Microscopic particle
crushing of sand subjected to high pressure one dimensional compression. Soil
Found, 41(1):69-82
[101] Naylor D, Marahna Das Neves D, Mattar D, Veiga Pinto A (1986) Prediction of
construction perfomance of Beliche Dam. Geotechnique 36(3):359-376
[102] Newmark N (1965) Effects of earthquakes on dams and embankments.
Geotechnique 15(2):139-160
[103] Niemunis A, Cudny M (1998) On hyperelasticity for clays. Comp Geot 23:221-236
198
[104] Núñez E (1995) Propiedades mecánicas de materiales granulares incoherentes. Acad
Nac Cs Ex Fís Nat Bs As 46:71-89
[105] Núñez E (2000) Investigación aplicada en la ingeniería geotécnica. Acad Nac Cs Ex
Fís Nat Bs As 52:29-66
[106] Núñez E (2007) Uncertainties and approximations in geotechnics. XIII PCSMGE,
Margarita (Venezuela) 26-39
[107] Núñez E (2007) Behavior of coarse alluvium slopes subjected to earthquakes. XIII
PCSMGE, Margarita (Venezuela) 862-867
[108] Okada N, Nemat-Nasser S (1994) Energy dissipation in inelastic flow of saturated
cohesionless granular media. Geotechnique 44(1):1-19
[109] Pestana J (1994) A unified constitutive model for clays and sands. DScThesis MIT
[110] Pestana J, Whittle A (1995) Compression model for cohesionless soils.
Geotechnique 45(4):611-631
[111] Pestana J, Whittle A (1999) Formulation of a unified constitutive model for clays
and sands. Int J Numer Anal Meth Geomech 23(1):1215-1243
[112] Pestana J, Whittle A, Salvati L (2002) Evaluation of a constitutive model for clays
and sands: Part I — sand behaviour. Int J Numer Anal Meth Geomech 26(2):1097-
1121
[113] Pietruszczak S, Mroz Z (2000) Formulation of anisotropic failure criteria
incorporating a microstructure tensor. Comp Geot 26:105-112
[114] Poorooshasb H, Holubec I, Sherbourne A (1966) Yielding and flow of sand in triaxial
compression. Can Geot J 3:179-190
[115] Poulos S (1981) The steady state of deformation. ASCE J Geot Eng Div 107(5):553-
562
199
[116] Puzrin A, Tatsuoka F (1998) Elastic strain energy potential for uncemented
granular materials. Soil Found 38(4):267-274
[117] Reynolds O (1886) Experiments showing dilatancy properties of granular material.
Proc Royal Inst London, 2:354-363, citado en [?]
[118] Richart F, Hall J, Woods D (1970) Vibrations of Soil Found. Int Series Theo Appl
Mech, Prentice-Hall
[119] Roberts J, De Souza J (1958) The compressibility of sands. Proc ASTM, 58:1269-
1277
[120] Roscoe K, Schofield A, Wroth C (1958) On the yielding of soils. Geotechnique
8(1):22-53
[121] Roscoe K, Schofield A, Thurairajah A (1963) Yielding of clays in state wetter than
critical. Geotechnique 13(3):211-240
[122] Roscoe K, Burland J (1968) On the generalized stress strain behaviour of wet clay.
Geotechnique 18(4):535-608
[123] Rowe P (1962) The stress dilatancy relation for static equilibrium of an assembly
of particles in contact. Proc Royal Soc London 269:500-527
[124] Rowe P, Barden L, Lee K (1964) Energy components during the triaxial cell and
direct shear tests. Geotechnique 14(3):247-261
[125] Rowe P (1971) Stress - strain relatioships for particulate materials at equilibrium.
Spec Conf Perfomance Earth Earth Supported Struct 3:327-359
[126] Saboya F, Byrne P (1993) Parameters for stress and deformation analysis of rockfill
dams. Can Geotl J 30:690-701
[127] Santamarina J C, Cascante G (1996) Stress anisotropy and wave propagation — A
micromechanical review. Can Geot J 33:770-782
200
[128] Santamarina J, Klein K, Fam M (2001) Soils and Waves: Particulate Materials
Behavior, Characterization and Process Monitoring. ISBN: 978-0-471-49058-6,
Wiley
[129] Schanz T, Vermeer P (1996) Angles of friction and dilatancy of sand. Geotechnique
46(1):145-151
[130] Schanz T, Vermeer P, Bonnier P (1999) The hardening soil model: formulation and
verification. Proc Plaxis Symp Beyond 2000 Comp Geotech Balkema, 55-58
[131] Schofield A, Wroth C (1968) Critical State Soil Mechanics. McGraw-Hill, London
[132] Schultze E, Moussa A (1961) Factors affecting the compressibility of sand, 5th Int
Conf Soil Mech Found Eng, Paris, I:335-340
[133] Seed H (1979) Considerations in the earthquake-resistant design of earth and rockfill
dams. Geotechnique 29(3), 215-283
[134] Seed H et al (1984) Moduli and damping factors for dynamic analyses of
cohesionless soils. Report UCB/EERC-84/14, Berkeley
[135] Sfriso A (1996) Una ecuación constitutiva para arenas. Enc Geot Arg GT96 I:38-46
[136] Sfriso A (2006) Calibration of ARENA for Nevada sand based on VELACS project
results. VII World Conf Comp Mech Los Angeles (CD-ROM)
[137] Sfriso A (2007) A constitutive model for sands: Evaluation of predictive capability.
XIII Conf Panam Mec Suelos Ing Geot, 242-247
[138] Sfriso A (2008) Numerical assessment of the deformation of CFRD dams during
earthquakes. XII Int Assoc Comp Meth Adv Geomech 4054:4061
[139] Sfriso A (2009) The friction angle and dilatancy of sands. XVII ICSMGE,
Alexandria, Egypt, 433-435
201
[140] Sfriso A, Weber G (2010) Formulation and validation of a constitutive model for
sands in monotonic shear. Acta Geotechnica, in press
[141] Sherard J, Cooke J (1987) Concrete-face rockfill dam: I Assessment. ASCE J Geot
Eng 113(10):1096-1112
[142] Tatsuoka F, Shibuya S (1992) Deformation characteristics of soils and rocks from
field and laboratory tests. Inst Ind Sc U Tokio, 37:1
[143] Tatsuoka F, Siddiquee M, Park C, Sakamoto M, Abe F (1993) Modelling stress
strain relations of sand. Soil Found 33(2):60-81
[144] Taylor D (1948) Fundamentals of soil mechanics. Wiley, New York
[145] Teachavorasinskum S, Shibuya S, Tatsuoka F (1991) Stiffness of sands in monotonic
and cyclic loading in simple shear. ASCE Geot Eng Congress Boulder (1):7863-7878
[146] Tobita Y, Yanagisawa E (1992) Modified stress tensors for anisotropic behavior of
granular materials. Soil Found 32(1):85-99
[147] Terzaghi K (1943) Theoretical soil mechanics. Wiley, New York
[148] Terzaghi K, Peck R (1948) Soil Mechanics in Engineering Practice. Wiley, New
York (1996), 3ra Ed
[149] Trautmann C, Kulhawy F (1987) CUFAD - A Computer program for compression
and uplift foundation analysis and design. Report EPRI EL-4540-CCM
[150] Uddin N, Gazetas G (1995) Dynamic response of CFRD to strong seismic
excitation. ASCE J Geot Eng 121(2):185-197
[151] Vardoulakis I, Sulem J (1995) Bifurcation analysis in geomechanics. Blackie
Academic & Professional, UK
202
[152] Vardoulakis I (1996) Deformation of water-saturated sand: I. uniform undrained
deformation and shear banding. Geotechnique 46(3):441-456
[153] Vardoulakis I (1996) Deformation of water-saturated sand: II. effect of pore water
flow and shear banding. Geotechnique 46(3):457-472
[154] Verdugo R, Ishihara K (1996) The steady state of sandy soils. Soil Found, 36(2):81-
91
[155] Vermeer P (1978) A double hardening model for sand. Geotechnique 28(4):413-433
[156] Vesic A, Clough G (1968) Behaviour of granular materials under high stresses.
ASCE J Soil Mech Found Div, 94(SM3):661-688
[157] Viggiani C, Hall S (2008) Full-field measurements, a new tool for laboratory
experimental geomechanics. IV Intl Symp Def Charac Geomat CH001:3-26
[158] Wan R, Guo P (1999) A pressure and density dependent dilatancy model for
granular materials. Soil Found 39(6):1-11
[159] Weber G, Anand L (1990) Finite deformation constitutive equations, and a time
integration procedure for isotropic, hyperelastic viscoplastic solids. Comp Meth
App Mech Eng 79:173—202
[160] Weber G, Lush A, Zavaliangos A, Anand L (1990) An objective time-
integration procedure for isotropic rate-independent and rate-dependent elastic-
plastic constitutive equations. Int J Plast 6:701—744
[161] Wong R, Arthur J (1985) Induced and inherent anisotropy in sand. Geotechnique,
35(4):471-481
[162] Woodward P, Molenkamp F (1999) Application of an advanced multi surface
kinematic constitutive soil model. Int J Numer Anal Meth Geomech 23(15):1995-
2043
203
[163] Wroth C, Houlsby G (1985) Soil mechanics property characterization and analysis
procedures. 11 Int Conf Soil Mech Found Eng, San Francisco, 1:1-55
[164] Yamashita S, Jamiolkowski M, Lo Presti D (2000) Stiffness nonlinearity of three
sands. J Geot Geoenv Eng ASCE 126(10):929-938
[165] Yasuda N, Matsumoto N (1993) Dynamic deformation characteristics of sands and
rockfill materials. Can Geot J 30:747-757
[166] Yasuda N, Matsumoto N (1994) Comparison of deformation characteristics of
rockfill materials using monotonic and cyclic loading laboratory tests and in situ
tests. Can Geot J, 31:162-174
[167] Yasufutu N, Murata H, Hyodo M, Hyde A (1991) A stress strain relationship for
anisotropically consolidated sand over a wide stress region. Soil Found, 31(4):75-92
[168] Yasufutu N, Murata H, Hyodo M (1991) Yield characteristics of anisotropically
consolidated sand under low and high stresses. Soil Found, 31(1):95-109
[169] Zlatovic S, Ishihara K (1997) Normalized behavior of very loose non plastic soils:
effects of fabric. Soil Found 37(4):47-56
204