Filtración por elementos finitos — documentación de OGR Slip2D

En vez de recibir como datos las presiones intersticiales de un talud, OGR Slip2D puede calcularlas resolviendo, por elementos finitos y sobre una malla del modelo, las ecuaciones del flujo de agua subterránea en suelo saturado y no saturado, en régimen permanente o transitorio. Los apartados siguientes plantean las ecuaciones y las decisiones numéricas que se toman para resolverlas (elemento, malla, condiciones de contorno, iteración sobre las superficies de rezume, funciones de permeabilidad e integración en el tiempo); el paso del campo resultante al análisis de equilibrio límite se describe al final y en la Sección 3.6.

3.7.1Ecuaciones de gobierno#

El agua se desplaza por los poros por efecto de los gradientes de altura piezométrica total. Con la altura de presión $P = u/\gamma_w$ , la altura total a la cota $y$ es

H = y + P, \qquad P = \dfrac{u}{\gamma_w}

El caudal por unidad de superficie, o velocidad de Darcy, sigue la ley de Darcy (1856), extendida al suelo no saturado mediante una conductividad que depende de la succión:

\mathbf{v} = -\,k_r(\psi)\,\mathbf{K}\,\nabla H

donde

$\mathbf{v}$: velocidad de Darcy (m/s)
$\mathbf{K}$: tensor de conductividad saturada del material (m/s), Sección 3.7.2
$k_r$: permeabilidad relativa, igual a 1 en suelo saturado (—), Sección 3.7.7
$\psi = -P$: altura de succión, donde la altura de presión es negativa (m)

La conservación de la masa (el flujo neto que sale de un volumen es igual al ritmo al que este pierde agua almacenada) conduce a la ecuación de Richards (1931), escrita aquí en forma mixta, con el contenido volumétrico de agua $\theta$ como magnitud almacenada:

\dfrac{\partial \theta}{\partial t} = \nabla\cdot\left[\,k_r(\psi)\,\mathbf{K}\,\nabla H\,\right]

En régimen permanente el miembro izquierdo se anula:

\nabla\cdot\left[\,k_r(\psi)\,\mathbf{K}\,\nabla H\,\right] = 0

Con $k_r = 1$ en todo el dominio, la ecuación es lineal y describe un flujo saturado y confinado; con una función de permeabilidad, la conductividad depende de la solución y la ecuación se ha de resolver de forma iterativa (Sección 3.7.6). La misma ecuación rige por encima y por debajo de la superficie freática, por lo que la superficie libre no es un contorno del problema, sino la línea $P = 0$ de su solución.

En OGR Slip2D

Las alturas se dan en la unidad de longitud del modelo. El motor de cálculo no depende de la unidad elegida, siempre que las conductividades y los tiempos de un análisis transitorio estén en unidades coherentes entre sí (por ejemplo, m/s y segundos). Las presiones intersticiales son $u = \gamma_w P$ , con el peso específico del agua intersticial definido en el proyecto.

Código: ogr_fem2d/solvers/seepage.py

3.7.2Conductividad anisótropa#

Cada material se define mediante su permeabilidad saturada $k_s$ (la conductividad principal mayor $K_1$ ), la relación $K_2/K_1$ y el ángulo $\theta_1$ de $K_1$ , medido desde el eje $+x$ . En los ejes del modelo, el tensor es el tensor principal girado ese ángulo:

\mathbf{K} = \mathbf{R}(\theta_1)\begin{bmatrix}K_1 & 0\\ 0 & K_2\end{bmatrix}\mathbf{R}(\theta_1)^{\mathsf T}, \qquad \begin{aligned} K_{xx} &= K_1\cos^2\theta_1 + K_2\sin^2\theta_1\\ K_{yy} &= K_1\sin^2\theta_1 + K_2\cos^2\theta_1\\ K_{xy} &= (K_1 - K_2)\sin\theta_1\cos\theta_1 \end{aligned}

Es simétrico y definido positivo si $K_1, K_2 > 0$ . En la zona no saturada todo el tensor se multiplica por $k_r$ , por lo que se conservan la relación de anisotropía y su dirección.

En OGR Slip2D

En Groundwater → Define Hydraulic Properties… se fijan Saturated permeability Ks (por defecto $1\times10^{-6}$ ), K2 / K1 (por defecto 1) y K1 angle (deg from +X) (por defecto 0). Un material sin propiedades hidráulicas se calcula con estos valores por defecto, y el mensaje posterior a Compute Groundwater menciona ese material.

Código: ogr_core/hydraulic/hydraulic_properties.py · ogr_slip2d/transient_stability.py

3.7.3Triángulos lineales#

El modelo se divide en triángulos de tres nodos con funciones de forma lineales $N_i$ . La formulación de Galerkin de la Ec. (3.7.4) sobre un elemento da su matriz de conductancia. Como los gradientes de las funciones de forma lineales son constantes en el elemento, la matriz se obtiene de forma analítica y exacta, sin integración numérica (Bathe y Khoshgoftaar, 1979):

k^{e}_{ij} = A_e\,k_r^{e}\,(\nabla N_i)^{\mathsf T}\,\mathbf{K}\,\nabla N_j, \qquad \dfrac{\partial N_i}{\partial x} = \dfrac{y_j - y_k}{2A_e}, \qquad \dfrac{\partial N_i}{\partial y} = \dfrac{x_k - x_j}{2A_e}

donde

$A_e$: área del elemento (m²)
$k_r^{e}$: permeabilidad relativa del elemento (—)
$(i, j, k)$: permutación cíclica de los tres nodos, numerados en sentido antihorario

Al ensamblar se obtiene el sistema lineal $\mathbf{K}\,\mathbf{H} = \mathbf{Q}$ , donde $\mathbf{Q}$ recoge los caudales nodales de las condiciones de caudal. La altura es lineal dentro de cada elemento, por lo que la velocidad de la Ec. (3.7.2) es constante en él. La permeabilidad relativa de un elemento se evalúa con la media de las alturas de presión de sus tres nodos, que en un campo lineal es el valor en el centroide.

En OGR Slip2D

El sistema se resuelve con un método directo para matrices dispersas. Las alturas prescritas se imponen por eliminación, lo que mantiene la simetría del sistema reducido. Un problema sin ninguna altura prescrita es singular (la altura solo queda definida salvo una constante), y el programa lo señala como un error.

Código: ogr_fem2d/solvers/seepage.py · ogr_fem2d/mesh/mesh.py

3.7.4Mallado#

La malla se genera región por región sobre los materiales del modelo, de modo que ningún elemento queda a caballo de un contorno de material y el material de cada elemento es exacto:

Cada arista de cada región se divide en segmentos de longitud próxima al tamaño objetivo $h$ . Todos los nodos del contorno se dan de alta en un único registro, indexado por sus coordenadas, de modo que una arista compartida por dos regiones recibe los mismos nodos en ambas: la malla es conforme en las interfaces entre materiales, sin necesidad de unirla después.
Cada región se rellena con puntos interiores en una retícula al tresbolillo, casi equilátera, de paso $h$ , que se mantienen al menos a $0.45h$ del contorno; después se unen todos los puntos de la región con una triangulación de Delaunay. Se descartan los triángulos cuyo centroide cae fuera de la región o en un hueco de esta.
En los triángulos cuyo ángulo menor es inferior a 25° se inserta su circuncentro como punto nuevo y la región se vuelve a triangular, hasta tres veces (Chew, 1989; Ruppert, 1995). No se inserta el circuncentro si queda fuera de la región, a menos de $0.35h$ del contorno o a menos de $0.45h$ de un punto existente, lo que garantiza que el refinamiento termine.

Si no se indica un tamaño, $h$ es el lado de los triángulos equiláteros que llenarían el área $A$ del modelo con el número de elementos solicitado $n_e$ :

h = \sqrt{\dfrac{4A}{\sqrt{3}\;n_e}}

En OGR Slip2D

Groundwater → Mesh → Generate FE Mesh… pide el número aproximado de elementos, 1000 por defecto (de 20 a 200 000). El ángulo mínimo (25°) y el número de pasadas de refinamiento (3) se pueden cambiar desde las interfaces de scripting. Una lente se malla una sola vez, como un hueco en la región que la rodea. Al generar una malla nueva se descartan las condiciones de contorno y los campos de la anterior, porque están ligados a la numeración de sus nodos.

Código: ogr_fem2d/mesh/generator.py · ogr_fem2d/mesh/mesh.py

Resultado de Interpret Groundwater en la presa de tierra del ejemplo 3: isolíneas de altura total sobre la malla triangular, vectores de flujo y superficie libre — Figura 3.7.1 Presa de tierra del ejemplo 3: sección homogénea de 15 m de altura, con el embalse a 13 m y un dren de pie de 16 m de longitud. Malla de 1 562 triángulos (873 nodos), isolíneas de altura total, vectores de flujo y superficie libre, calculados en 25 iteraciones.

3.7.5Condiciones de contorno#

A los nodos del contorno se les pueden aplicar seis tipos de condición (Tabla 3.7.1). Tres prescriben la altura, dos un caudal y la sexta corresponde a una posible superficie de rezume.

Tabla 3.7.1 Condiciones de contorno del análisis de filtración. Los caudales son positivos hacia el interior del modelo.

Condición	Identificador	Tipo	Prescribe
Total Head	`total_head`	altura	$H$ igual al valor
Pressure Head	`pressure_head`	altura	$H = y + P$ , con $P$ igual al valor
Zero Pressure	`zero_pressure`	altura	$H = y$
Nodal Flow Rate	`nodal_flow`	caudal	un caudal por unidad de espesor del modelo en el nodo
Infiltration	`infiltration`	caudal	un flujo por unidad de longitud de contorno, repartido entre los dos nodos de cada arista como $qL/2$ , que es la carga consistente de un elemento lineal
Unknown (P=0 or Q=0)	`unknown`	unilateral	o bien $P = 0$ con salida de agua, o bien flujo nulo (Sección 3.7.6)

Al generar una malla, cada nodo del contorno recibe una condición por defecto. Los nodos de los bordes izquierdo, derecho e inferior del rectángulo envolvente del modelo reciben un caudal nodal nulo, con lo que esos bordes son impermeables; todos los demás, incluida la superficie del terreno con el paramento del talud, reciben la condición Unknown. Las condiciones por defecto no prescriben ninguna altura, así que antes de resolver el problema hay que prescribir una en algún punto.

Un embalse se define por su nivel y por el lado en que se encuentra. El programa prescribe una altura total igual al nivel en el perímetro mojado, que determina recorriendo el contorno desde el borde vertical de ese lado y subiendo por la superficie del terreno mientras su cota no supere el nivel. El recorrido se detiene en el primer nodo situado por encima del nivel y nunca entra en la base ni en el borde opuesto del modelo. La parte sumergida del borde vertical se incluye por defecto. Para bajar un embalse basta repetir la asignación con un nivel más bajo. La Figura 3.7.2 muestra estas condiciones sobre un terraplén y la solución que calcula el programa con ellas. Una condición de caudal nodal también se puede marcar como posible superficie de rezume, con lo que su nodo pasa a formar parte de la conmutación de la Sección 3.7.6.

Figura 3.7.2 Terraplén con un embalse a 25 m, resuelto con OGR Slip2D 0.1.235 sobre una malla de 448 elementos, con la función de van Genuchten y sus parámetros por defecto. Altura total en el perímetro mojado y condiciones por defecto en el resto del contorno: Unknown en el tramo restante de la superficie del terreno y flujo nulo en el borde derecho y en la base. La superficie libre es la línea

P = 0

de la solución, y los puntos son los nodos mantenidos en

P = 0

en la superficie de rezume, que empieza en el paramento de aguas abajo, por encima del pie.

En OGR Slip2D

Con Groundwater → Set Boundary Conditions… se asigna una condición a un lado del modelo (Left edge, Right edge, Bottom edge, Ground surface) o se coloca un embalse a la izquierda o a la derecha; Restore defaults devuelve las condiciones descritas más arriba. Las interfaces de scripting permiten además actuar sobre nodos individuales y sobre el tramo de contorno que sigue una polilínea. Cada conjunto de condiciones pertenece a una malla: si nombra nodos que la malla actual no tiene, se sustituye por las condiciones por defecto.

Código: ogr_fem2d/solvers/seepage.py · ogr_fem2d/solvers/bc_targets.py

3.7.6Iteración de Picard y superficies de rezume#

Hay dos no linealidades que se resuelven a la vez: la conductividad depende de la succión, y cada nodo Unknown ha de quedar como superficie de rezume o como impermeable. El programa parte de la solución lineal con $k_r = 1$ y con todos los nodos Unknown impermeables. Después, en cada iteración $m$ , evalúa $k_r$ en cada elemento a partir de las alturas actuales, resuelve el sistema lineal con los nodos que en ese momento están en la superficie de rezume fijados en $P = 0$ , y relaja la actualización:

\mathbf{K}\!\left(\mathbf{H}^{m}\right)\mathbf{H}^{*} = \mathbf{Q}, \qquad \mathbf{H}^{m+1} = (1 - \omega)\,\mathbf{H}^{m} + \omega\,\mathbf{H}^{*}

donde

$\mathbf{H}^{*}$: alturas de la solución lineal con las conductividades de la iteración $m$
$\omega$: factor de relajación, entre 0.05 y 1 (—)

La subrelajación evita que la iteración oscile cuando la función de permeabilidad tiene mucha pendiente, como ocurre en las arenas.

La superficie de rezume#

Donde la superficie libre llega a un paramento expuesto, el agua sale del modelo a presión atmosférica. La extensión de esa superficie de rezume no se conoce de antemano, por lo que los nodos Unknown obedecen una condición unilateral: o bien $P = 0$ con salida de agua, o bien flujo nulo con $P \le 0$ . El programa la determina por conmutación nodal (Neuman, 1973; Bathe y Khoshgoftaar, 1979): un nodo libre cuya altura de presión se hace positiva se fija en $P = 0$ , y un nodo fijado cuya reacción nodal indica que entra agua en el modelo se libera y pasa a flujo nulo. Dos bandas de tolerancia, una altura de presión del 2 % del tamaño del elemento y una reacción de $10^{-3}$ veces la mayor reacción de la primera solución, impiden que la conmutación oscile sin cesar, y un nodo que ha conmutado 25 veces se congela.

La iteración converge cuando el mayor cambio de altura es inferior a la tolerancia y ningún nodo ha conmutado. Después se comprueba el estado final con la propia condición unilateral: si algún nodo Unknown aún la incumple, el cálculo se da por no convergido, sea cual sea el motivo por el que se detuvo la iteración. La malla nunca se mueve; la superficie libre se traza a continuación como la línea donde la altura de presión cambia de signo.

En OGR Slip2D

Groundwater → Compute Groundwater resuelve el problema en régimen permanente con $\omega = 0.4$ , un máximo de 200 iteraciones y una tolerancia de $10^{-5}$ m en la altura; si el cálculo no converge, se indica el motivo. La superficie libre que dibuja Interpret Groundwater es, sobre una serie de líneas verticales, el punto más alto en que la altura de presión cambia de signo.

Código: ogr_fem2d/solvers/seepage.py · ogr_slip2d/transient_stability.py

3.7.7Funciones de permeabilidad#

La permeabilidad relativa $k_r = k/k_s$ de cada material es una de siete funciones de la succión. Dos de ellas se expresan en altura de succión $h$ , en metros de columna de agua, y las demás en succión matricial $s$ , en kPa; el motor de cálculo usa para cada función su propia variable, con $s = \gamma_w h$ . Todas valen 1 con succión nula y están acotadas inferiormente por $k_{r,\min} = 10^{-6}$ , de modo que la matriz de conductancia nunca se hace singular en zonas muy secas. La Tabla 3.7.2 recoge las funciones con sus parámetros por defecto, y la Figura 3.7.3 las representa.

Constante. $k_r = 1$ : el material conduce como si estuviera saturado cualquiera que sea la succión, y el problema es lineal si todos los materiales usan esta función.

Simple. Función propia del programa, que no procede de la bibliografía: la permeabilidad disminuye linealmente en escala logarítmica $D$ órdenes de magnitud a lo largo de un intervalo de succión $s_\text{ref}$ y se mantiene constante a partir de él,

k_r = 10^{-D\,\min\left(1,\ s/s_\text{ref}\right)}

con $(s_\text{ref},\ D)$ = (100 kPa, 1) para el tipo de suelo general, (10 kPa, 4) para la arena, (50 kPa, 3) para el limo, (1500 kPa, 2) para la arcilla y (100 kPa, 3) para el suelo franco.

Brooks-Corey. A partir de la curva de retención potencial de Brooks y Corey (1964), con el índice de tamaño de poros $\lambda$ y la presión de burbujeo $s_b$ :

k_r = \begin{cases} 1 & s \le s_b\\[4pt] \left(\dfrac{s_b}{s}\right)^{2 + 3\lambda} & s > s_b \end{cases}

Fredlund-Xing. La permeabilidad relativa toma la forma de la curva característica suelo–agua de Fredlund y Xing (1994), normalizada y sin el factor de corrección $C(\psi)$ de la ecuación original:

k_r = \dfrac{1}{\left\{\ln\!\left[e + (s/A)^{B}\right]\right\}^{C}}

Gardner. La función de Gardner (1958), en altura de succión:

k_r = \dfrac{1}{1 + a\,h^{n}}

van Genuchten. La expresión analítica de van Genuchten (1980), deducida de su curva de retención con el modelo de conductividad de Mualem (1976), expresada en altura de succión:

S_e = \left[1 + (\alpha h)^{n}\right]^{-m}, \qquad k_r = S_e^{1/2}\left[1 - \left(1 - S_e^{1/m}\right)^{m}\right]^{2}, \qquad m = 1 - \dfrac{1}{n}

La expresión analítica es válida para $m = 1 - 1/n$ . El programa permite dar a $m$ un valor distinto, y entonces aplica la misma expresión con el valor dado.

Definida por el usuario. Una tabla de puntos $(s_i,\ k_i)$ . La permeabilidad relativa es la permeabilidad tabulada dividida por la del primer punto, con interpolación lineal en la succión y logarítmica en la permeabilidad; vale 1 por debajo de la primera succión y mantiene la última relación más allá de la última succión.

Tabla 3.7.2 Funciones de permeabilidad, variable en que se expresa cada una y sus parámetros por defecto.

Función	Identificador	Variable	Parámetros (por defecto)
Constant	`constant`	—	—
Simple	`simple`	$s$ (kPa)	tipo de suelo (general)
Brooks-Corey	`brooks_corey`	$s$ (kPa)	$\lambda = 0.6$ , $s_b = 30$ kPa
Fredlund-Xing	`fredlund_xing`	$s$ (kPa)	$A = 50$ kPa, $B = 2$ , $C = 1$
Gardner	`gardner`	$h$ (m)	$a = 0.01$ m⁻ⁿ, $n = 2$
van Genuchten	`van_genuchten`	$h$ (m)	$\alpha = 3.6$ m⁻¹, $n = 1.56$ , $m = 1 - 1/n$
User Defined	`user_defined`	$s$ (kPa)	tabla de al menos dos puntos

Figura 3.7.3 Permeabilidad relativa frente a la succión matricial para todas las funciones, salvo la tabla definida por el usuario, con los parámetros por defecto y calculada con la misma llamada que usa el motor de cálculo (OGR Slip2D 0.1.235). El eje superior da la misma succión como altura de agua. La función de van Genuchten por defecto, que corresponde a un suelo franco, alcanza el límite inferior

k_{r,\min}

a unos 87 kPa (8.8 m de columna de agua).

En OGR Slip2D

La función se elige en Groundwater → Define Hydraulic Properties…; Plot la dibuja en su propia variable y Pick carga valores representativos para Brooks-Corey, Fredlund-Xing, Gardner y van Genuchten (las texturas de van Genuchten son las de Carsel y Parrish, 1988, convertidas a m⁻¹). En una tabla User Defined, la conductividad es $k_s$ multiplicada por la relación tabulada.

Código: ogr_core/hydraulic/permeability_models.py · ogr_core/hydraulic/hydraulic_properties.py

3.7.8Análisis transitorio#

Un análisis transitorio sigue la evolución del campo en el tiempo a partir de un estado inicial, por etapas cuyas condiciones de contorno pueden ser distintas (un desembalse, un período de lluvia). Integra la Ec. (3.7.3) en forma mixta, con el agua almacenada como magnitud acumulada. Esta forma conserva el balance de masa de la solución discreta, mientras que la forma en que $\partial\theta/\partial t$ se sustituye por $C\,\partial P/\partial t$ acumula errores de balance de masa (Celia, Bouloutas y Zarba, 1990).

Almacenamiento#

El agua almacenada es una única función de la altura de presión, que reúne el almacenamiento capilar de la zona no saturada y el almacenamiento elástico de la saturada:

W(P) = \begin{cases} \theta(\psi = -P) & P < 0\\ \theta_s + S_s\,P & P \ge 0 \end{cases} \qquad\qquad \dfrac{\mathrm{d}W}{\mathrm{d}P} = \begin{cases} \max\!\left[C(\psi),\ 10^{-3}S_s\right] & P < 0\\ S_s & P \ge 0 \end{cases}

donde

$\theta_s,\ \theta_r$: contenidos volumétricos de agua saturado y residual (—)
$S_s$: almacenamiento específico del suelo saturado (m⁻¹)
$C(\psi)$: capacidad hídrica específica, la pendiente de la curva de retención (m⁻¹)

Por debajo del nivel freático el contenido de agua es constante; sin el almacenamiento específico, la ecuación transitoria se reduciría allí a la de régimen permanente. Cualquiera que sea su función de permeabilidad, el contenido de agua de cada material sigue la curva de retención de van Genuchten (1980), con los valores $\alpha$ , $n$ y $m$ de ese material (Figura 3.7.4):

\theta(\psi) = \theta_r + (\theta_s - \theta_r)\left[1 + (\alpha\psi)^{n}\right]^{-m}, \qquad C(\psi) = (\theta_s - \theta_r)\,m\,n\,\alpha\,(\alpha\psi)^{n-1}\left[1 + (\alpha\psi)^{n}\right]^{-m-1}

Figura 3.7.4 Curva de retención que usa cada material en un análisis transitorio, aquí con los parámetros por defecto (

\alpha = 3.6

m⁻¹,

n = 1.56

\theta_s = 0.40

\theta_r = 0.05

), calculada con OGR Slip2D 0.1.235. Su pendiente es el coeficiente de almacenamiento de la zona no saturada.

Integración en el tiempo#

El tiempo se integra con el método de Euler implícito (hacia atrás), que es incondicionalmente estable. Cada paso se resuelve con la iteración de Picard modificada de Celia, Bouloutas y Zarba (1990), que linealiza el almacenamiento en torno a la iteración en curso pero trata de forma explícita la variación del agua almacenada:

\left[\dfrac{\mathbf{M}\,\mathbf{C}^{m}}{\Delta t} + \mathbf{K}^{m}\right]\mathbf{H}^{m+1} = \mathbf{Q} + \dfrac{\mathbf{M}\,\mathbf{C}^{m}}{\Delta t}\,\mathbf{H}^{m} - \dfrac{\mathbf{M}\left(\mathbf{W}^{m} - \mathbf{W}^{n}\right)}{\Delta t}

donde

$n,\ m$: nivel de tiempo anterior e iteración de Picard en curso en el nuevo nivel
$\mathbf{M}$: matriz de masa concentrada (diagonal), en la que cada nodo recibe un tercio del área de cada elemento al que pertenece (m²)
$\mathbf{C}^{m}$: diagonal de los coeficientes de almacenamiento nodales $\mathrm{d}W/\mathrm{d}P$ en la iteración en curso (m⁻¹)
$\mathbf{W}^{m},\ \mathbf{W}^{n}$: agua almacenada en los nodos en la iteración en curso y en el nivel de tiempo anterior (—)
$\mathbf{K}^{m}$: matriz de conductancia con las permeabilidades relativas de la iteración en curso
$\Delta t$: paso de tiempo

La actualización se relaja igual que en régimen permanente, la conmutación de la superficie de rezume se aplica en cada iteración de cada paso, y un paso converge cuando el mayor cambio de altura es inferior a la tolerancia y ningún nodo ha conmutado.

Etapas y pasos de tiempo#

Cada etapa termina en el tiempo $t_k$ y puede llevar sus propias condiciones de contorno; si no las tiene, usa las vigentes. Cada etapa se divide en pasos iguales, en número fijo o en un número elegido de modo que un frente de humedecimiento no pueda recorrer mucho más de un elemento por paso:

n_k = \min\!\left[200,\ \max\!\left(4,\ \left\lceil \dfrac{t_k - t_{k-1}}{t_e} \right\rceil\right)\right], \qquad t_e = \dfrac{C_{\min}\,h^2}{k_{s,\max}}

donde

$h$: tamaño objetivo de los elementos de la malla (m)
$k_{s,\max}$: mayor permeabilidad saturada del modelo
$C_{\min}$: menor coeficiente de almacenamiento de los materiales con altura de presión de −1 m, sin superar nunca el valor 1 (m⁻¹)

El estado inicial es la solución en régimen permanente con las condiciones de contorno iniciales, que son las guardadas con el análisis transitorio o, si no guarda ninguna, las vigentes. Una etapa de duración nula devuelve el estado que hereda.

En OGR Slip2D

Las etapas (tiempo, etiqueta, Calculate SF y condiciones propias) se definen en Groundwater → Transient Groundwater…. En Analysis → Project Settings → Transient se fijan la tolerancia ( $10^{-5}$ m), el número máximo de iteraciones de Picard por paso (30) y el número de pasos de tiempo por etapa (Auto por defecto). El motor de cálculo transitorio usa $\omega = 0.5$ , también en el estado inicial permanente, que puede necesitar hasta 200 iteraciones. Cada nodo toma los parámetros de almacenamiento de uno de los elementos que lo comparten, el primero de la lista de elementos. En cada etapa se registran el agua almacenada y su variación, para comprobar el balance de masa. El análisis transitorio, el exceso de presión intersticial y el desembalse rápido son incompatibles entre sí.

Código: ogr_fem2d/solvers/seepage.py · ogr_core/hydraulic/hydraulic_properties.py · ogr_slip2d/transient_stability.py · ogr_core/project/settings.py

3.7.9Acoplamiento con el análisis de estabilidad#

Las presiones intersticiales nodales $u_i = \gamma_w (H_i - y_i)$ pasan al análisis de equilibrio límite mediante los materiales de tipo fem, interpoladas en el punto medio de la base de cada dovela (Ec. (3.6.9)); los valores negativos se tratan según la resistencia de los suelos no saturados (Sección 3.6.9). Las condiciones de contorno definen dónde hay agua sobre el talud: hay agua embalsada allí donde la altura total o la altura de presión prescrita supera la cota del nodo (Sección 3.6.7).

En un análisis transitorio, cada etapa marcada con Calculate SF se analiza con los métodos y la búsqueda del proyecto, con el campo de esa etapa y con las condiciones que lo produjeron, de modo que una etapa de desembalse no incluye ni las presiones intersticiales ni el peso del embalse ya vaciado. El resultado es la evolución del factor de seguridad en el tiempo (Sección 3.8.7).

Figura 3.7.5 Superficie crítica del talud de aguas abajo de la misma presa (Spencer,

F = 1.799

) y presión intersticial tomada del campo de filtración en la base de cada una de sus 30 dovelas. Por encima de la superficie libre la presión es nula, y también lo es sobre el dren.

En OGR Slip2D

Los factores de seguridad de cada etapa se guardan con su resultado y se muestran en Interpret Groundwater. Si ningún material es de tipo fem, no se calculan los factores de las etapas y un aviso explica por qué.

Código: ogr_slip2d/transient_stability.py · ogr_core/hydraulic/ponded_water.py

3.7.10Verificación#

El motor de cálculo se verifica con soluciones analíticas e identidades (Sección 4.3):

el flujo confinado unidimensional, en el que la altura es lineal y el caudal sigue la ley de Darcy;
una prueba del parche con un tensor totalmente anisótropo ( $K_{xy} \ne 0$ ), en la que se reproducen exactamente un gradiente de altura uniforme y su velocidad de Darcy;
los suelos estratificados, cuya conductividad equivalente es la media armónica de las capas para el flujo perpendicular a ellas y la media aritmética para el flujo paralelo, con conservación del caudal a través de un contraste de permeabilidad de 100;
el caudal a través de una presa rectangular sobre una base impermeable, que tiene una expresión cerrada exacta y se alcanza a medida que se hace más abrupta la función de permeabilidad, y la salida de la superficie libre por el paramento de aguas abajo por encima del nivel de aguas abajo;
la respuesta de un acuífero confinado a un cambio brusco de altura, que sigue la función de error complementaria, y la solución transitoria a tiempos largos, que debe reproducir la de régimen permanente;
la curva de retención y su derivada.

Criterio de ingeniería

Las permeabilidades y las condiciones de contorno de un análisis de filtración suelen ser más inciertas que su error numérico. Compara los niveles de agua calculados con las lecturas de los piezómetros donde las haya, y comprueba hasta qué punto es sensible el factor de seguridad a los contrastes de permeabilidad que controlan la superficie libre.

3.7.11Referencias#

Bathe, K.J. & Khoshgoftaar, M.R. (1979). Finite element free surface seepage analysis without mesh iteration. International Journal for Numerical and Analytical Methods in Geomechanics, 3(1), 13–22.
Brooks, R.H. & Corey, A.T. (1964). Hydraulic Properties of Porous Media. Hydrology Paper 3, Colorado State University, Fort Collins.
Carsel, R.F. & Parrish, R.S. (1988). Developing joint probability distributions of soil water retention characteristics. Water Resources Research, 24(5), 755–769.
Celia, M.A., Bouloutas, E.T. & Zarba, R.L. (1990). A general mass-conservative numerical solution for the unsaturated flow equation. Water Resources Research, 26(7), 1483–1496.
Chew, L.P. (1989). Guaranteed-Quality Triangular Meshes. Technical Report TR 89-983, Department of Computer Science, Cornell University, Ithaca.
Darcy, H. (1856). Les fontaines publiques de la ville de Dijon. Victor Dalmont, Paris.
Fredlund, D.G. & Xing, A. (1994). Equations for the soil-water characteristic curve. Canadian Geotechnical Journal, 31(4), 521–532.
Gardner, W.R. (1958). Some steady-state solutions of the unsaturated moisture flow equation with application to evaporation from a water table. Soil Science, 85(4), 228–232.
Mualem, Y. (1976). A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resources Research, 12(3), 513–522.
Neuman, S.P. (1973). Saturated-unsaturated seepage by finite elements. Journal of the Hydraulics Division, ASCE, 99(HY12), 2233–2250.
Richards, L.A. (1931). Capillary conduction of liquids through porous mediums. Physics, 1(5), 318–333.
Ruppert, J. (1995). A Delaunay refinement algorithm for quality 2-dimensional mesh generation. Journal of Algorithms, 18(3), 548–585.
van Genuchten, M.Th. (1980). A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal, 44(5), 892–898.