Revista Internacional de Métodos Numéricos para Calculo y Diseño e.n Ingeniería. Vol.

7, 1, 55-66( 1991)

ESTRATEGIAS DE SOLUCION PARA

ANALISIS DE CONSOLIDACION
ACOPLADOS EN MEDIOS POROSOS

A. SAETTA
B.A. SCHREFLER
Y
R. VITALIANI

Istituto di Scienza e Tecnica delle Costruzioni,

Vía Marzolo 9,
35131 Padova, Italia.

RESUMEN

Se revisan brevemente las estrategias empleadas usualmente para resolver las ecuaciones
acopladas del problema de la consolidación de suelos saturados. Luego se introduce un
nuevo esquema basado en la solucion particionada del conjunto de ecuaciones acopladas y
se desarrollan los criterios de convergencia y estabilidad necesarios.

SUMMARY

Commonly used solution strategies of the coupled set of equations governing consolidation
of saturated soil are briefly reviewed. Then a new procedure based on the partitioned solution
of the coupled set of equations is introduced and the necessary stability and convergence criteria
are derived.

INTRODUCCION

Las primeras aplicaciones del método de elementos finitos al análisis de

consolidación de suelos elásticos isotérmicos aparecieron hacia el final de la década de
los sesenta'12. Desde entonces el tema ha suscitado mucho interés y se han conseguido
considerables progresos en los procedimientos de solución y en la complejidad de los
problemas resueltos. Descripciones de la literatura del tema más recientes pueden
encontrarse en las referencias [3] y [4]. Enfocamos aquí nuestro interés en los
procedimientos de solución de las ecuaciones acopladas que manejan el problema. Este
interés está justificado por el hecho de que los análisis de consolidación que siguen las
líneas de la teoría de Biot5, implican un número de grados de libertad por nodo que
generalmente es la suma de los grados de libertad de los dos campos interactuantes
(campos de desplazamiento y flujo), y más aún, a menudo conducen a largas historias
en el tiempo. Así, es importante disponer de esquemas de solución eficientes para el

Recibido: Noviembre 1989

OUniversitat Politecnica de Catalunya (España) ISSN 0213-1315

 A. SAETTA, B. SCHREFLER Y R. VITALIANI

análisis de problemas de ingeniería real que a menudo son no lineales y exigen una
extensión a problemas tridimensionales .
Recientemente Park y Felippa6 han propuesto un esquema de solución particionada
para el análisis de interacción entre dos o más campos. Este procedimiento se ha
adaptado con éxito para el análisis de consolidación no isotérmica4J. Mediante l a
introducción de una partición se descompone el problema en una parte de consolidación
(campos de desplazamiento y flujo) y otra térmica (campo de temperaturas), puede
conseguirse una destacable eficiencia computaciona18. Aquí investigamos la posibilidad
de extender más aún la partición, es decir, descomponer el problema de consolidación
en dos campos separados. La posible existencia del campo de temperaturas no complica
la situación, ya que el acoplamiento entre este campo y los de presión y desplazamientos
es débil8. Podemos así enfocar nuestra atención solamente en la consolidación
isotérmica. El proceso resultante también resulta atractivo para el análisis de
consolidación isotérmica por si misma, especialmente si se emplean rnicrocomputadores.
Es importante recordar que descomponer no significa desacoplar las ecuaciones, sino
que consiste en procesar los vectores de estado del problema acoplado mediante módulos
separados, llamados analizadores de campo. Estos analizadores de campo se sincronizan
para operar de manera secuencia1 o paralela6. En lo que sigue se emplea el esquema
secuencial ya que es el más apropiado paia ordenadores pequeños
Antes de tratar con la solución particionada se establecen simplemente las
ecuaciones del problema y se revisan los esquemas de solución más comunes para
destacar sus diferencias.

ECUACIONES DISCRETAS DEL PROBLEMA

Las ecuaciones que describen el problema de la consolidación isotérmica pueden

obtenerse ya sea imponiendo la nulidad de las derivadas de Gateaux de un
principio variacional conveniente113,bien mediante el empleo del principio de trabajos
virtuales4J0, o mediante. una formulación de residuos ponderados4~". La teoría
subyacente se basa o bien en la teoría del continuo de mezclas o en la teoría del continuo
de materiales con m i c r o e s t r u c t ~ r a ' ~ ~ ' ~ .
Las ecuaciones que describen el problema se establecen simplemente aquí, como se
han obtenido en la referencia [14], luego de discretizar el espacio mediante elementos
finitos.

donde u y p son, los valores nodales de las variables básicas, ,es decir el vector
desplazamientos del esqueleto sólido y el vector de presiones de poros del fluido
saturado, H y M son las matrices de rigidez y masa usuales de un problema
de equilibrio desacoplado en mecánica de sólidos y S y H son las matrices de
capacit ancia y conductividad que corresponden a un problema de filtración desacoplado.
Generalmente estas matrices son simétricas, definidas y positivas, aún si K resulta no
simétrica en el caso de plasticidad no asociada.
 ANALISIS D E CONSOLIDACION ACOPLADOS EN MEDIOS POROSOS 57

Finalmente Q es la matriz de interacción entre los dos campos. La expresión

detallada de,estos dos campos puede encontrarse en la literatura". Para los desarrollos
posteriores la ecuación (1) se escribe en una forma más concisa como

El significado de las matrices y vectores anteriores surge inmediatamente de la

inspección de la ecuación (1).
En este trabajo la matriz de masa se desprecia y el sistema de ecuaciones (2) resulta
simplemente:

Para resolver el problema de valores iniciales, deben transformarse las ecuaciones

diferenciales (3) en un sistema algebraico discreto mediante la integración en el tiempo.
Los esquemas temporales adoptados por diversos autores se revisan en la referencia (31.
Aquí utilizamos un esquema temporal de dos niveles (e-standard) y el sistema (3)
resulta:

donde At indica el incremento de tiempo, u,+l, p,+l, y u, y p, son los desplazamientos

y presiones calculados en los pasos de tiempo n + 1 y n respectivamente, 0 es un
parámetro libre que controla la estabilidad y la precisión del esquema. Generalmente
se supone O < 0 < 1.
El lector puede verificar fácilmente que en nuestro caso el acoplamiento entre los
dos campos ocurre solamente en el transitivo y desaparece cuando se alcanza el estado
estacionario.

PROCEDIMIENTOS DE SOLUCION

Diversos autores han abordado la solución directa de la ecuación (3), o de

ecuaciones ~imilares'.~.ll.Esta técnica proporciona ambas variables al mismo tiempo y
no requiere operaciones matriciales adicionales, una vez definidas las matrices.
Smith15 y Krauselo, reducen el número de grados de libertad en el caso lineal,
mediante la eliminación o bien de la presión de fluido p , o del desplazamiento del suelo
u . La primera técnica ha sido denominada método de los desplazamientos por Krauselo
y la segunda método de las presiones. En cada instante sólo se conoce la historia de
una variable. Si se ha determinado el campo de desplazamientos, el campo de presiones
se calcula mediante diferenciación numérica, con la consiguiente pérdida de exactitud.
58 A. SAETTA, B. SCHREFLER Y R. VITALIANI

En el método de las presiones, donde el número de grados de libertad se reduce a

un mínimo, los desplazamientos se obtienen mediante la solución de las ecuaciones de
equilibrio. En el método de las presiones las matrices no son dispersas.
El esquema de solución particionada, que se tratará a continuación, mantiene la
ventaja de tratar con matrices más pequeñas sin perder su característica de banda.
Más aún, provee la historia temporal para ambas variables de estado.
+
Comenzamos con la ecuación (4). La matriz D = (B B A t C ) del primer miembro
de la ecuación (4) se descompone según

donde

La ecuación (6) se sustituye ahora cn la ecuación (4) y D ~ x " + 'se transfiera al

segundo miembro, donde xn+l se sustituye mediante un predictor xi+' adecuado.

PREDICTORES

Hay una gran variedad de predictores de la solución en el instante n + 1. Dos de

los más empleados son:

con O < 7 < 1. A menudo se hace 7 = 0 o y = O. En el Último caso el predictor viene

dado por la Última solución conocida xn. La expresión general de los predictores puede
escribirse según:

donde el predictor se expresa como una función de la solución anterior. Formas típicas
de la ecuación (10) pueden obtenerse de las fórmulas de integración numérica según:

Otras formas del predictor incluyen también la primera derivada temporal de la

variable de estado x y son de la forma geneial
 ANALISIS DE CONSOLIDACION ACOPLADOS EN MEDIOS POROSOS 59

Estas son comunes en problemas con derivadas temporales de segundo orden6. Los
valores de las derivadas pueden obtenerse de la ecuación (3). Esto requiere sin embargo
que se conozca B-', lo que no es tan atractivo en nuestro caso. En lo que sigue
empleamos predictores dados por la ecuación (9).

SOLUCION ESCALONADA

Mediante la adopción de la ecuación (9) y descomponiendo la matriz D como en

la ecuación (7),podemos escribir el sistema siguiente

que pueden reescribirse como

Este esquema escalonado comienza resolviendo la primera ecuación del sistema (3)
con y = O (última solución), luego se resuelve la segunda ecuación con h+lya calculado
y con el vector actual del parámetro y.
A continuación el esquema es iterativo y continúa hasta que la solución converge.
Dado que el acoplamiento de los dos campos es normalmente fuerte, al menos
en materiales como suelos, se necesitan iteraciones dentro de cada paso de tiempo
para lograr convergencia. Los requisitos de convergencia dan un límite inferior para el
tamaño del paso de tiempo como se verá en el próximo apartado.
Las iteraciones no son una desventaja primordial, ya que en muchos problemas
prácticos algún grado de no linealidad está presente. Tales no linealidades requieren
generalmente iteraciones en dos esquemas temporales. No se requieren operaciones
matriciales adicionales excepto las necesarias para la descomposición del' primer
miembro de la ecuación (12).
 A . SAETTA, B . SCHREFLER Y R. VITALIANI

ANALISIS DE LA CONVERGENCIA

Se investiga ahora la convergencia de las iteraciones dentro de cada caso de tiempo

para el caso lineal.

Las ecuaciones (13) resultan

"+'= bu, + hp, + dp:+, + F,
ep:$ = fu, + gp, - fui+, + F p (14)

donde s indica la iteración.

Remarcando los términos implicados en la iteración, la ecuación (14) puede
escribirse simplemente como

y sustituyendo u:+, = u 2 de la primera ecuación en'la segunda

Considerando la segunda ecuación, podemos escribir

d+l
ePn+i -
-
+G
donde

donde J = (-e-' fa-ld) es la matriz iterativa del sistema.

Puede mostrarsels que el método iterativo converge si

donde p(J) es el el mayor autovalor de J .

Despreciando la compresibilidad de los gramos del sólido y la del fluido (S = O),
que usualmente resulta aceptable para suelos, la matriz de iteración resulta
 ANALISIS D E CONSOLIDACION ACOPLADOS EN MEDIOS POROSOS

y la ecuación (20) se reduce a

Debido al hecho que en nuestro caso J es real y simétrico, para los valores de 6 y
7 elegidos, puede conseguirse la convergencia del proceso iterativo sólo si se cumple

lo que representa un límite inferior para el tamaño del incremento de tiempo. Cuan
restrictiva resulta esta condición depende de los autovalores de la matriz H-'QTK-'Q
resultando así una función de las propiedades del material, permeabilidades y datos
geométricos.

ANALISIS D E ESTABILIDAD

La estabilidad numérica de las ecuaciones particionadas no sólo depende de la

fórmula de integración sino también de la descomposición y de la expresión del
predictor6. Como ejemplo se examina la estabilidad numérica de las ecuaciones (3),
reescritas según (13). Para ello se busca una solución no trivial del tipo

para el caso de respuesta libre, es decir cuando todos los términos de carga son nulos
en la ecuación (1). El parámetro y es el coeficiente de amplificación de la solución. La
sustitución de la ecuación (24) en la (13) proporciona la ecuación característica

P(X)xn = o
donde

+
P(X) = [DT y D E + (1- 7)DE - B + (1- B)AtC] (26)
Una solución no trivial de la (25) requiere

detP(X) = O (27)
Para la investigación de la ecuación (25) es conveniente introducir la transformación
 A. SAETTA, B. SCHREFLER Y R. VITALIANI

que transforma el círculo de estabilidad ) A ) 5 1 sobre el semiplano deñnido por

Re(z) 5 O". La ecuación (27) puede así reescribirse como

Desarrollando el determinante la ecuación (29) queda

Si las matrices a, b y c son definidas positivas, el polinomio (30) no tiene raices

reales positivas (Teorema de Routh).
La condición de matriz definida positiva para a , b y c, puede asegurarse imponiendo

EJEMPLOS NUMERICOS

La efectividad del procedimiento de solución particionado propuesto se demuestra

mediante dos ejemplos simples tomados de la literatura

Ejemplo lZ0
En la Figura 1 se muestra una columna de suelo sujeta a cargas superficiales q,
que varían con el tiempo (Figura 2). Las condiciones de contorno y las propiedades del
material se muestran en la misma figura.

m
0.1,O O. 4,O

Permeabilidad (dirección x) K, : = 1 m/seg

Permeabilidad (dirección y) Ky: = 1 m/seg
Módulo de Elasticidad E: = 1 KN/m2
Coeficiente de Poisson u: = O
Módulo voluiiiétrico del fluido Kf : = .1.106 KN/m2
Porosidad por. : = 0.3

21 23
-1 t=- 1,o
Figura 1. Ejeiiiplo 1: Malla y datos.
 - --

ANALISIS D E CONSOLIDACION ACOPLADOS EN MEDIOS PO&OSOS

'~oad K N / ~ ~

1.0

Time t (sec )
--

10

Figura 2. Función de carga.

Para garantizar la convergencia (p(J) 5 l), el tamaño del paso de tiempo se a

elegido igual a 90 segundos. Nótese que éste es un incremento de tiempo aceptable en
problemas de consolidación.
Los resultados numéricos se ilustran en las Figuras 3 y 4 en función de los
desplazamientos y presión nodal respectivamente.

Figura 3. Desplazamientos Nodales. Figura 4. Presiones nodales.

Ejemplo 2'l
En la Figura 5 se muestra un plano de suelo semi-infinito, sujeto a la misma carga
superficial que en el ejemplo 1. Las características del suelo se indican en dicha Figura
y también en este caso el tamaño del paso de tiempo debe ser 90 seg.
Igual que en el ejemplo anterior los desplazamientos y las presiones nodales se
muestran en las Figuras 6 y 7 respectivamente.
Los mismos resultados numéricos de los dos casos previos se han obtenido con un
procedimiento global ensayado en varios casos.
 A. SAETTA, B . SCHREFLER Y R. VITALIANI

Permeabilidad (dirección x) K, : = 0.001 m/seg

Permeabilidad (dirección y) Ky : = 0.001 m/seg
Módulo de Elasticidad E : = 300 KN/m2
Coeficiente de Poisson v : = 0.2
Módulo volumétrico del fluido Kf : = . l e l o 6 KN/m2
Porosidad por. : = 0.3

Figura 5. Ejeiiiplo 2: Malla y datos.

Figura 6. Desplazamientos nodales. Figura 7. Presiones nodales.

CONCLUSIONES

Las conclusiones d e los análisis de convergencia y estabilidad n o deberían cambiar

e n el caso d e descomponer las ecuaciones del problema de la consolidación n o isotérmica
e n tres campos individuales. E s t o es así porque el acoplamiento entre la temperatura
y los otros dos campos es débil y así n o resulta necesario iterar dentro del intervalo de
 ANALISIS DE CONSOLIDACION ACOPLADOS EN MEDIOS POROSOS 65

El producto entre la temperatura y las matrices relevantes del acoplamiento

aparece en el término independiente, que no se actualiza en la iteración y no afecta a
la matriz de iteraciónlg.
Se ha mostrado que el esquema particionado permite operar sobre matrices más
pequeñas que en el cálculo directo, y además mantiene la estructura en banda de las
matrices. Este procedimiento también proporciona simultáneamente las variables de
estado. El límite inferior al tamaño del incremento temporal que surge de los requisitos
de convergencia del procedimiento iterativo constituye sin embargo un inconveniente.
Este límite inferior es especialmente alto para el caso de suelos rígidos, como surge de
inspeccionar su estructura.
Una vez que se ha calculado el mayor autovalor p(J) de la matriz iterativa el
tamaño del incremento de tiempo puede elegirse para satisfacer la ecuación (23). Si el
incremento de tiempo resulta demasiado grande deben aplicarse otros procedimientos.
Sin embargo en muchos problemas de ingeniería la solución no es necesaria para períodos
de tiempo pequeños, y así el procedimiento particionado presentado en este trabajo
puede aplicarse eficientemente para resolver el sistema de ecuaciones acopladas que
describe el problema.

REFERENCIAS

1. R.S. Sandhu y E.L. Wilson, "Finite Element Analysis of seepage in elastic media", Int. J.
Engng. Mech. Div., ASCE, Vol. 95, pp. 641-652, (1968).
2. I.T. Christian y I.W. Boehnier,"Plane Strain Consolidation by Finite Elementsn, Int. J.
Soil Mech. Found. Div., ASCE, Vol. 96, pp. 1435-1475, (1970).
3. R.S. Sandhu, "Finite Element Analysis of Coupled Deformation and Fluid Flow in Porous
Media", "Numerical Methods in Geomechanics", (ed. I.B. Martins), D. Reidel Pubiishing
Company, pp. 203-227, (1987).
4. R.S. Lewis y B.A. Schrefler, "The Finite Element Method in the Deformation and
Consolidation of Porous Media", J. Wiley, (1987).
5. M.A. Biot, "Consolidation Settlement under a Rectangular Load Distributionn, Int. J. of
Appl. Physics, Vol. 12, pp. 426-430, (1941).
6. K.C. Park y C.A. Felippa, "Partitioned Analysis of Coupled Systems", "Computational
Methods for Transient Analysis)), (T. Belytschko y T.J. Hughes, eds.), Elsevier, pp. 157-
219, London, (1983).
7. B.A. Schrefler, "A Partitioned Solution Procedure for Geothermai Reservoir Analysis,
CANM, Vol. 1,pp. 47-59, (1985).
8. B.A. Schrefler, R.W. Lewis y C.E. Majorana, "A Comparison Between Different Simulation
for Non-isothermal Consolidation", Proc. Europ. Conf. Num. Methods in Geom., Vol. 2,
Univ. de Stuttgart, (1986).
9. B.A. Schrefler, "Partitioned Solution Procedure and Infinite Elements in the Mechanics
of Porous Media", Proc. VIII Conf. Comp. Methods in Structural Mechanics, Iadwisin,
Polonia, (1987).
10. G. Krause, "Finite Element Schemes for Porous Elastic Media", Int. J. Engng. Mech.
Div., ASCE, Vol. 104, pp. 606-620, (1978).
11. I.H. Prevost, "Consoiidation ofinelastic Porous Media", Int. J. Engng. Mech. Div., ASCE,
Vol. 1 0 7 , pp. 169-186, (1981).
 A. SAETTA, B. SCHREFLER Y R. VITALIANI

12. R. De Boer y W. Ehlers, "A historical review of the formulation of Porous Media Theories",
Acta Mecánica, Vol. 74, pp. 1-8, (1988).
12. R. De Boer y W. Ehlers, "A Historical Review of the Formulation of Porous media
Theories", Acta Mecánica, Vol. 74, pp. 1-8, (1988).
13. A. Bedford y D.S. Drumheller, "Theories of Inmiscible and Structured Mixtures", Int. J.
Engng. Sci., Vol. 21, pp. 863-960, (1983).
14. O.C. Zienkiewicz, K.M. Leung y C.T. Chang, "Liquefaction and Permanent Deformation
under Dynamic Conditions - Numerical Solution and Constitutive Relations", "Soil
Mechanics - Transient and Cyclic Loads", (G.N. Pande y O.C. Zienkiewicz eds.), J. Wiley,
(1982).
15. I.M. Smith, "Transient Phenomena of Offshore Foundations", "Numerical Methods in
Offshore Engineering", (O.C. Zienkiewicz, R.W. Lewis y K.G. Ctagg, eds.), J. Wiley, pp.
483-513, (1978).
16. G.H. Golub y C. Van Loan, "Matriz Computations", Hopkins University Press, (1983).
17. C.K. Park, "A Partitioned Transient Analysis Procedure for Coupled Filed Problems:
Stability Analysis", Int. J. Appl. Mech., Vol. 47, pp. 370-376, (1980).
18. B.A. Schefler y R. Vitaliani, "Partitioned Solution Procedure for Educational Software on
Microcomputers: The Case of a Coupled Consolidation Analysis", (EPMESC), Macao,
(1985).
19. K.C. Park, C.A. Felippa y J.A. DeRuntz, "Stabilization of Staggered Solution Procedure
for Fluid Structures Interaction Analysis", "Computational Methods for Fluid Structures
Interaction Analysis", (T. Belytschko y T.L. Geers eds.), ASME, AMD, Vol. 26, pp.
94-124, (1977).
20. I.M. Smith y D.V. Griffiths, "Programming the Finite Element Method", J. Wiley, (1982).
21. L. Simoni y B.A. Schrefler, "Land Subsidence Analysis Using Mapped Infinite Elements",
"Microcomputers in Engineering, (B.A. Schrefler y R.W. Lewis eds.), Pineridge Press, pp.
189-202, Swansea, (1986).