MODELAMIENTO NUMÉRICO

EN MECÁNICA DE ROCAS (MI6061)

MÉTODO DE LAS
DIFERENCIAS FINITAS
FLAC/FLAC3D

Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 1

Aspectos básicos
• Se resuelven las ecuaciones dinámicas de movimiento aún en
problemas cuasi-estáticos, dados valores iniciales y/o condiciones
de borde

• Para modelar la respuesta estática de un sistema se utiliza un

esquema de relajación en el cual se emplea amortiguación que
permite absorber energía cinética

• Esta aproximación tiene ventajas para modelar problemas que

involucran inestabilidades físicas (ej: colapsos) de una manera más
eficiente y realista que otros esquemas, por ej., métodos de
solución matriciales

Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 2

Grilla de diferencias finitas
Gridpoint: Almacena cantidades
vectoriales, por ejemplo:
fuerzas, velocidades,
desplazamientos, etc.

Zone: Almacena cantidades

escalares y tensoriales,
por ejemplo: propiedades
del material, tensor de
esfuerzos, tensor de
deformaciones, etc.

Las ecuaciones se plantean en puntos discretos del espacio en

términos de variables de campo (esfuerzo o desplazamientos)

Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 3

 Ecuaciones
Movimiento

Ley de Newton

𝑑𝑢
𝐹 𝑡 = 𝑚𝑢 = 𝑚
𝑑𝑡

Para un cuerpo continuo, ésta se puede generalizar como:

𝜌 : densidad de masa
𝜕𝜎𝑖𝑗 𝜕𝑢𝑖 𝑥𝑗 : componentes del vector de coordenadas cartesianas
+ 𝜌𝑔𝑖 = 𝜌
𝜕𝑥𝑗 𝜕𝑡 𝜎𝑖𝑗 : componentes del tensor esfuerzos
𝑔𝑖 : componentes de aceleración de gravedad (fuerzas de cuerpo)

Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 4

Ecuaciones
Relación esfuerzo-deformación

1 𝜕𝑢𝑖 𝜕𝑢𝑗
𝑒𝑖𝑗 = +
2 𝜕𝑥𝑗 𝜕𝑥𝑖

𝑒𝑖𝑗 tasa de deformación

𝑢𝑖 velocidades

𝜎𝑖𝑗 ≔ 𝑀 𝜎𝑖𝑗 , 𝑒𝑖𝑗 , 𝑘

𝑀() : forma funcional de la ley constitutiva

𝑘 : parámetro(s) de historial que pueden o no existir dependiendo de la ley en particular

Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 5

 Ecuaciones
Relación esfuerzo-deformación

Ejemplo: elasticidad isótropa

2
𝜎𝑖𝑗 ≔ 𝜎𝑖𝑗 + 𝛿𝑖𝑗 𝐾 − 𝐺 𝑒𝑘𝑘 + 2𝐺 𝑒𝑖𝑗 Δ𝑡
3

𝛿𝑖𝑗 Delta de Kronecker

Δ𝑡 Paso de tiempo

𝐺, 𝐾 Modulo de corte y compresión volumétrica, respectivamente

𝐸 𝐸
𝐺= 𝐾=
2 1+𝜈 3 1 − 2𝜈

Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 6

Ecuaciones
Condiciones de borde

Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 7

Marcha explicita en el tiempo
Cada elemento aparece físicamente aislado de sus vecinos durante un
incremento/paso de tiempo, Dt

Fuerzas constantes
durante este calculo

Dt

Tasa de deformación constantes

durante este calculo

Se escoge un paso de tiempo pequeño de manera que la información

no pueda físicamente traspasarse de un elemento a otro en ese
intervalo de tiempo (Todo material tienen una velocidad máxima a la
cual la información puede propagarse).
Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 8
Marcha explicita en el tiempo
El concepto central es que la velocidad de onda de calculo siempre es
mayor a la velocidad de onda física

Así, las ecuaciones siempre operan en valores conocidos que son

constantes durante el cálculo

Ventajas:
• No se requiere de un proceso de iterativo para calcular los esfuerzos
a partir de deformaciones en el elemento
• Se puede implementar cualquier modelo constitutivo sin ajustar el
algoritmo de solución

Desventaja:
• Debido a lo pequeño de Dt un gran número de pasos deben ser
ejecutados para conseguir una solución

Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 9

Formulación numérica
• Las ecuaciones de movimiento del medio continuo son
transformadas a ecuaciones discretas en los nodos, dando como
resultado un sistema de ecuaciones diferenciales, que son resueltas
numéricamente mediante una integración explicita en el tiempo

• El medio continuo es discretizado en elementos cuadrilaterales de

tasa de deformación constante, cuyos vértices son los nodos

• Internamente los elementos son subdivididos en dos sets de

elementos triangulares de deformación constante superpuestos

a, b, c, d: sub-elementos triangulares

Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 10

Formulación numérica – Grilla

Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 11

Formulación numérica – Ecuaciones
Las ecuaciones diferenciales para un triangulo son obtenidas del
teorema de divergencia generalizado de Gauss

Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 12

Formulación numérica – Ecuaciones

El teorema de Gauss permite definir de manera discreta el valor

promedio del gradiente de f sobre el área A para un sub-elemento
triangular

• Ds largo del lado del triangulo

• La sumatoria ocurre sobre los tres lados del triangulo
• es el promedio sobre el lado

Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 13

Formulación numérica – tasa de deformación
La ecuación anterior permite escribir 𝑒𝑖𝑗 en términos de velocidades
nodales para el sub-elemento triangular

1 𝜕𝑢𝑖 𝜕𝑢𝑗
𝑒𝑖𝑗 = +
2 𝜕𝑥𝑗 𝜕𝑥𝑖

Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 14

Formulación numérica – tasa de deformación

Se remplaza el vector de velocidad promedio de cada lado por f

𝑢𝑝𝑟𝑜𝑚 =

• La sumatoria es sobre los tres lados de la zona triangular

• (a) y (b) son dos nodos consecutivos en un lado
• La tasa de deformación de la zona es el promedio de los valores de los sub-elementos
triangulares

Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 15

Formulación numérica – esfuerzos
Con 𝑒𝑖𝑗 se obtiene el nuevo estado de esfuerzos de la ley constitutiva

𝜎𝑖𝑗 ≔ 𝑀 𝜎𝑖𝑗 , 𝑒𝑖𝑗 , 𝑘

Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 16

Formulación numérica – fuerzas nodales
Una vez que los esfuerzos son calculados en los lados se determinan
las fuerzas equivalentes aplicadas en cada nodo

• El esfuerzo en cada lado del triangulo es

equivalente a dos fuerzas iguales que actúan
en los nodos de cada lado
• Cada vértice del triangulo recibe dos
contribuciones de fuerzas (una de cada lado)

Las fuerzas de los triángulos que se encuentran en cada nodo son

promediadas para obtener la fuerza nodal del cuadrilátero

Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 17

Formulación numérica – ecuaciones de movimiento

• En cada nodo las fuerzas de todos los cuadriláteros que lo rodean

son sumadas para dar el vector de fuerzas nodales resultante

• Este vector incluye las fuerzas gravitacionales , donde

mg = 1/3 S masas de los triángulos conectados al nodo

• Si el cuerpo está en equilibrio o en flujo plástico estacionario

en el nodo será cero. Sino, el nodo será acelerado de acuerdo a:

Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 18

Formulación numérica – amortiguamiento mecánico

El proceso de solución dinámica es amortiguado para alcanzar el

equilibrio estático

Analogía mecánica

• Se comienza con una fuerza no

balanceada y se resuelve al
estado de equilibrio
• El amortiguamiento es utilizado
para alcanzar la solución

Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 19

Formulación numérica – amortiguamiento mecánico

• La fuerza de amortiguamiento en un nodo es proporcional a la

magnitud de la fuerza no balanceada

• La dirección de la fuerza de amortiguamiento es tal que

 1, si y  0

sgn  y    1, si y  0
 0, si y  0


• Fd es la fuerza de amortiguamiento, a es una constante (0.8 en

FLAC) y mn es una masa nodal ficticia (por definir)

Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 20

Formulación numérica – estabilidad numérica

• El procedimiento explicito no es incondicionalmente estable

• Requiere que la velocidad del frente de cálculo sea mayor a la
máxima velocidad en que la información se propaga
• Un paso de tiempo menor a un valor crítico debe ser seleccionado
• La condición de estabilidad para un solido elástico discretizado en
elementos de tamaño Dx es

• C es la velocidad máxima a la cual la información puede

propagarse, igual a la velocidad de las ondas-P, Cp

Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 21

Formulación numérica – estabilidad numérica

• Para un sistema individual masa-resorte, la condición de estabilidad

es

• En un sistema general, el paso de tiempo crítico esta relacionado

con el periodo natural menor de sistema (Tmin)

• Sin embargo, resulta impráctico calcular los valores propios del

sistema completo, por lo que el paso de tiempo crítico debe ser
estimado

Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 22

Formulación numérica – estabilidad numérica

• Para encontrar una solución estática las masas nodales mn pueden

ser vistas como factores de relajación en las ecuaciones de
movimiento

• Es decir, mn es ajustado para encontrar la velocidad optima de

convergencia

• Las fuerzas gravitacionales no son afectadas por este escalamiento

• La convergencia optima se obtiene cuando los valores locales del

paso de tiempo crítico son iguales, es decir, cuando los periodos
naturales son iguales en todas las partes del sistema

Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 23

Formulación numérica – estabilidad numérica

• Por conveniencia se asume un paso de tiempo unitario y se ajustan

las masas nodales para obtener este valor

• Se asume un Factor de Seguridad de 0.5 en el paso de tiempo

crítico

• Para una zona triangular de área A se obtiene que la mínima

distancia de propagación de la zona queda dada por .
Luego, el paso de tiempo crítico:

Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 24

Formulación numérica – estabilidad numérica

• Sustituyendo Dt =1 y

• La masa de una zona es

• Utilizando que la masa del nodo (mgp) es 1/3 de la masa de la zona

• La masa nodal queda dada por

Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 25

Ciclo de calculo para cada Dt
• Condiciones iniciales (esfuerzos in-situ)
• Condiciones de borde

A Fuerzas nodales
Ecuaciones de movimiento

esfuerzos aceleraciones

Modelo
consitutivo
  H (D , t , Dt ) Dt  dt

deformaciones velocidades

DL
L
desplazamientos
 dt

• El ciclo de cálculo es repetido para cada nodo

• Un ciclo dura Dt

Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 26

Criterios de convergencia – análisis estático
• El método de solución dinámica proporciona una solución estática
(se minimizan fuerzas inerciales) cuando la razón de las fuerzas no
balanceadas alcanza un valor pequeño (~0.001 o menor)

• Esto es comparable al nivel de error residual o criterio de

convergencia definido por los métodos de solución matricial usado
en muchos programas de elementos finitos

• El nivel de error se cuantifica por la razón de las fuerzas no

balanceadas. El proceso de solución estática termina cuando el
error esta por debajo de cierto valor.

Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 27

Criterios de convergencia – análisis estático
• Se monitorea la razón de la fuerza no balanceada, Fi, y la magnitud
de la fuerza aplicada, Fm

• Cuando Fi/Fm < 0.001, el modelo se considera en estado de

equilibrio

Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 28

 Tipos de diferencias finitas
Diferencia en la formulación de desplazamiento del gridpoint

Métodos Diferencias Finitas

Euleriano Lagrangiano
Max Vector = 1.833E+00
Max Vector = 2.076E+00

No actualiza las coordenadas Actualiza las coordenadas

por deformación por deformación

Problemas con Esencial para problemas de

deformación pequeña, gran deformación o colapso
tiempo de cálculo menor

Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 29

Referencias
• Potts, D., Axelsson, K., Grande, L., Schweiger, H. and Long, M. (2002). Guidelines for the use of advanced numerical
analysis. Thomas Telford.

Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 30

 Diferencias Finitas v/s Elementos Finitos
Diferencia en la formulación de diferencia finitas y elementos finitos

Diferencias Finitas Elementos Finitos

Solución Explicita Solución Implícita

Iterativo hasta solución final Método de inversión de Matriz

Forma Incremental Esencial para problemas de

Relaciones constitutivas grandes deformación o colapso

Resuelve ecuaciones de Resuelve matriz completa,

movimiento hasta alcanzar el equilibrio global
equilibrio

Ideal para problemas No- Ideal para problemas

elásticos y grandes elásticos y pequeñas
deformaciones deformaciones
Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 31
Analogía

Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 32

Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 33
Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 34
 FLAC 3D
Interpretación de Resultados
Criterios de Detención y Comprobación de estabilidad numérica del modelo

FLAC3D en un programa que “converge” a una solución mecánica de los

problemas en estudio. Dicha convergencia debe ser evaluada a partir de criterios
prácticos que permitan tiempos de operación computacional razonables.

Los principales criterios de convergencia empleados en análisis estático son:

• Control de la fuerza no balanceada del modelo o error de cierre

(UNBALANCE FORCE y MECH RATIO)

• Control de deformaciones y velocidad de deformación.

Javier Vallejos, Ph.D. – U. de Chile / Ing. de Minas MI6061 (Otoño, 2016) 35