Agua subterránea y comprobaciones analíticas

El motor de agua subterránea y otras herramientas de OGR Slip2D se contrastan con resultados que no admiten aproximación: una solución analítica, una identidad que debe cumplirse exactamente o un requisito de coherencia.

Cada comprobación tiene un valor exacto o calculable de forma independiente (Sección 4.1.2). Lo exigido es la tolerancia escrita en el test; lo obtenido es el valor de la magnitud que el test comprueba, medido al ejecutar el modelo del test con la versión 0.1.235. El archivo de cada test se enlaza en su sección, y cada test puede ejecutarse por separado (Sección 4.1.7).

4.3.1Las comprobaciones de un vistazo#

4.3.2Flujo confinado: la ley de Darcy en una dimensión#

En una capa confinada y saturada, limitada por dos bordes con altura impuesta y sin flujo a través de los demás contornos, la altura es lineal y el caudal obedece la ley de Darcy (1856):

donde

$H_0,\ H_1$

alturas totales impuestas en $x = 0$ y en $x = L$ (m)

$K$

conductividad hidráulica de la capa (m/s)

$h$

espesor de la capa, transversal al flujo (m)

$Q$

caudal por unidad de longitud de la sección (m²/s)

$u,\ \gamma_w,\ y$

presión intersticial (kPa), peso específico del agua (kN/m³) y cota (m)

Una altura lineal pertenece al espacio de funciones de los triángulos lineales, por lo que la solución por elementos finitos de este problema es exacta salvo errores de redondeo, sea cual sea la malla; cualquier error en el ensamblaje, en las condiciones de contorno o en el cálculo del caudal aparece como una diferencia.

4.3.3Medios estratificados: medias armónica y aritmética#

Dos capas de conductividades $K_1$ y $K_2$ dan lugar a dos casos límite, ambos consecuencia de la ley de Darcy y de la continuidad. Cuando el flujo atraviesa las capas (en serie), el mismo caudal pasa por todas ellas y las pérdidas de altura se suman, de lo que resulta la media armónica de las conductividades; cuando discurre a lo largo de las capas (en paralelo), todas están sometidas a la misma altura y los caudales se suman, de lo que resulta la media aritmética:

donde

$L_i$

longitud de la capa $i$ en la dirección del flujo, con $L = \sum L_i$ (m)

$b_i$

espesor de la capa $i$ transversal al flujo, con $h = \sum b_i$ (m)

$r_i$

resistencia hidráulica de la capa $i$ (s)

$H_{\mathrm{int}}$

altura en la interfaz de dos capas en serie (m)

La comprobación emplea un contraste de 100, que exige mucho al método numérico: casi toda la pérdida de altura tiene que producirse en la capa de baja conductividad, y el caudal debe conservarse a través de la interfaz.

4.3.4Anisotropía: la prueba del parche con un tensor completo#

La conductividad de un material es un tensor cuyos valores principales $K_1$ y $K_2 = (K_2/K_1)\,K_1$ están girados un ángulo $\theta_1$, que es el que forma $K_1$ con el eje $+x$, lo que da $K_{xy} = (K_1 - K_2)\sin\theta_1\cos\theta_1$ (Ec. (3.7.5)). Un campo de alturas lineal satisface la ecuación en régimen permanente para cualquier tensor constante, puesto que su gradiente es constante:

Si ese campo se impone en todo el contorno, la solución por elementos finitos debe reproducirlo en cada nodo, y la velocidad de Darcy debe ser en cada elemento el vector constante de la Ec. (4.3.4). Es la prueba del parche de un elemento. Con $K_{xy} \ne 0$ pone a prueba el giro del tensor y los términos cruzados de la matriz del elemento, algo que los casos isótropos no pueden hacer, y por eso comprueba tanto el giro como el ensamblaje de Galerkin.

Exigido. Una diferencia de altura inferior a 10⁻⁹ m en cada nodo, para los cinco tensores; una diferencia entre la velocidad de cada elemento y la exacta inferior a 10⁻⁶, relativa a su componente mayor; y, en el caso isótropo, un módulo del gradiente igual a $\sqrt{A^2 + B^2}$ con una diferencia inferior a 10⁻⁹. Obtenido. Diferencias de altura de hasta 2.5×10⁻¹⁴ m y de velocidad de hasta 3.6×10⁻¹⁴, en los cinco.

4.3.5Flujo transitorio: la respuesta a un escalón de altura#

En un acuífero confinado y totalmente saturado, la ecuación transitoria se reduce a la ecuación de difusión, con la difusividad hidráulica $D = K/S_s$, donde $S_s$ es el almacenamiento específico. Para un acuífero semiinfinito con altura uniforme $H_0$, cuyo contorno en $x = 0$ se eleva a $H_1$ en $t = 0$, la solución es el perfil clásico de la función de error complementaria (Carslaw y Jaeger, 1959):

El frente avanza como $2\sqrt{D\,t}$. La solución solo vale mientras el frente está lejos del otro extremo del acuífero, y es la solución analítica, no el motor de cálculo, la que deja de ser aplicable más allá: el test compara solo donde se cumple la hipótesis de acuífero semiinfinito.

El mismo archivo comprueba que el frente avanza con el tiempo y que un almacenamiento diez veces mayor lo retrasa. El mismo acuífero se construye también mediante la capa de operaciones, con el mismo límite del 1 % en los mismos dos tiempos, y para un tercer tiempo de 10⁶ s debe reproducirse el perfil lineal del régimen permanente con una diferencia inferior a 10⁻³ del escalón de altura.

4.3.6Tiempo largo: el régimen permanente#

Con el paso del tiempo, el agua almacenada deja de cambiar, y la solución transitoria tiene que aproximarse a la solución en régimen permanente de las mismas condiciones de contorno, que es la de la Ec. (3.7.4):

El problema no saturado en régimen permanente no tiene solución cerrada, por lo que esta es una comprobación de coherencia y no de una solución cerrada: relaciona el esquema transitorio con el motor de cálculo en régimen permanente, que a su vez se comprueba con el caudal exacto de una presa rectangular (Sección 4.3.11).

Modelo. Un rectángulo de 20 m por 12 m, con una malla de 0.9 m. La función de permeabilidad es la de Gardner, con $K_s = 1\times10^{-5}$ m/s, $a = 1$ y $n = 3$; la curva de retención es la de van Genuchten, con $\alpha = 0.1$ m⁻¹ y $n = 2$, $\theta_s = 0.4$ y $\theta_r = 0.05$; el almacenamiento específico es $1\times10^{-4}$ m⁻¹. Al comienzo, la altura es de 10 m en ambos lados; las condiciones finales son 10 m a la izquierda y 3 m a la derecha, impuestas hasta esos niveles y libres por encima de ellos. La solución tras una única etapa de 10⁹ s se compara con la solución en régimen permanente de las condiciones finales, calculada con el motor de cálculo no saturado en régimen permanente (relajación 0.4, como máximo 200 iteraciones, tolerancia 10⁻⁶).

Exigido. Que la mayor diferencia de altura total entre las dos soluciones sea inferior a 0.05 m. Obtenido. 2.8×10⁻⁷ m. Mediante la capa de operaciones, el acuífero confinado de la sección anterior debe quedar, en $t = 10^6$ s, a menos de 10⁻³ del escalón de altura, 0.01 m, de su perfil lineal en régimen permanente. Otro test exige que el drenaje disminuya el agua almacenada de forma monótona.

Test: tests/test_transient_v130.py, TestAsymptoticConsistency.

4.3.7La malla de elementos finitos conserva el área#

Una malla conforme de triángulos lineales que cubre exactamente una región debe sumar el área de esa región:

La identidad no se cumple si la malla tiene huecos, solapes o triángulos que salen del dominio, así que ninguno de estos defectos puede pasar la comprobación. Se comprueba en tres modelos:

• un cuadrado de 10 m por 10 m, con un tamaño objetivo de elemento de 1.5 m: 100 m²;
• dos cuadrados apilados que comparten una arista, con lo que se comprueba la conformidad en la interfaz entre dos regiones: 200 m²;
• un talud de la batería de tests con tres regiones de material, con un objetivo de 200, 800 y 2 000 elementos: el área total debe ser la misma en cada refinamiento e igual al área de las regiones tal como la calcula el programa.

Exigido. Una diferencia inferior a 10⁻⁶ m² en cada caso, y una dispersión de las áreas inferior a 10⁻⁶ m² entre los tres refinamientos. Obtenido. Ninguna diferencia, a la precisión de la coma flotante, en ninguno de ellos. El mismo archivo comprueba que ninguna arista pertenece a más de dos elementos, que ningún elemento cruza la interfaz entre regiones y que los nodos de una arista compartida son comunes a ambos lados, que todo elemento tiene área positiva y está numerado en sentido antihorario, que el ángulo mínimo no baja de 20° en el cuadrado ni de 15° en el talud, y que las funciones de forma reproducen exactamente un campo lineal: el gradiente de $3x - 2y + 7$ resulta correcto con una diferencia inferior a 10⁻⁶.

4.3.8Retroanálisis de soportes: las fuerzas activa y pasiva coinciden con un factor de 1#

El retroanálisis fija el factor de seguridad en un valor objetivo $F_t$ y despeja la fuerza horizontal que hay que añadir para alcanzarlo (Sección 3.10.6). Si $R$ y $D$ son las sumas resistente y desestabilizadora del método, evaluadas en $F_t$, y $a$ el factor de brazo de esa sección,

las dos hipótesis dan la misma fuerza con $F_t = 1$, sean cuales sean el talud, la superficie y el método, y la fuerza pasiva es la mayor de las dos para cualquier objetivo superior a 1 cuando hace falta alguna fuerza. La identidad dejaría de cumplirse ante casi cualquier error algebraico, porque $D - R/1$ y $1\cdot D - R$ son la misma expresión. La fuerza es nula cuando $F_t$ es el factor de seguridad de la superficie sin soporte.

Modelos. En el test del motor de cálculo, un talud de la batería de tests con tres regiones de material, el círculo de centro (88.0, 70.5) m y radio 47.212 m, 25 dovelas, el método de Bishop simplificado y la fuerza aplicada a una cota de 25 m. En el test mediante la capa de operaciones, la operación back_analysis_run sobre un talud de 10 m de altura con 1V:2H ($c' = 6$ kPa, $\phi' = 22^\circ$, $\gamma = 19$ kN/m³) con una Grid Search de 3 × 3 intervalos, un objetivo de 1.0 y la fuerza a una cota de 5 m.

Exigido. Con $F_t = 1.0$, una diferencia entre la fuerza activa y la pasiva inferior a 10⁻⁶ (en unidades de fuerza por unidad de longitud), en ambos tests. En el test del motor de cálculo, además: una fuerza pasiva mayor que la activa para objetivos de 1.2, 1.5 y 2.0; una fuerza que crece con el objetivo (1.0, 1.2, 1.4, 1.6); y una fuerza nula, y no negativa, para un objetivo que la superficie ya supera, con menos de 1 kN/m justo en el propio factor sin soporte. Obtenido. Una diferencia de 0 con $F_t = 1$ y una relación $T_{\mathrm{pas}}/T_{\mathrm{act}}$ de 1.2, 1.5 y 2.0 para esos objetivos, con un error de hasta 4×10⁻¹⁶.

4.3.9Constantes de Hoek-Brown para roca intacta#

La Parameter Calculator obtiene las constantes del criterio de Hoek-Brown generalizado a partir del índice de resistencia geológica (GSI), la constante de la roca intacta $m_i$ y el factor de alteración $D$ (Hoek, Carranza-Torres y Corkum, 2002; Ec. (3.5.14) y Ec. (3.5.15)):

Para $\mathrm{GSI} = 100$ y $D = 0$ se anulan los exponentes de $m_b$ y $s$, y las dos exponenciales de $a$ son el mismo número, por lo que el macizo rocoso se reduce a la roca intacta: $m_b = m_i$, $s = 1$ y $a = 1/2$, de forma exacta.

Exigido. Para $\mathrm{GSI} = 100$, $m_i = 10$ y $D = 0$: $m_b$ a menos de 10⁻⁹ de 10, $s$ a menos de 10⁻⁹ de 1 y $a$ a menos de 10⁻³ de 0.5. Obtenido. Exactamente 10, 1 y 0.5. Un segundo test usa un valor de Hoek et al. (2002): para $\mathrm{GSI} = 50$, $m_i = 10$ y $D = 0$, las constantes son $m_b = 1.6768$ con una tolerancia de 10⁻³, $s = 0.003866$ con 10⁻⁵ y $a = 0.5057$ con 10⁻³, y el programa da 1.67677, 0.0038659 y 0.50573. Otro test pasa por el botón del diálogo de materiales y comprueba, con la solución cerrada $\exp[(\mathrm{GSI} - 100)/(9 - 3D)]$, el valor de $s$ que llega al material, y no solo el que muestra la calculadora; con un GSI bajo, cuya $s$ es muy pequeña, el valor debe llegar íntegro.

4.3.10Ida y vuelta de DXF#

Un modelo exportado a DXF e importado de nuevo debe recuperar la misma geometría. En rigor, el dibujo no es idéntico vértice a vértice, y no tiene por qué serlo: el importador divide los contornos en sus cruces, con lo que el contorno externo gana nodos donde lo corta un contorno de material o una línea de agua. Lo que tiene que cumplirse es que la forma no cambie:

• el área encerrada es la misma, calculada con $A = \tfrac12\left|\sum_i (x_i\,y_{i+1} - x_{i+1}\,y_i)\right|$;
• todos los vértices originales siguen presentes;
• todo vértice añadido está sobre un segmento original.

Es una condición más exigente que pedir listas idénticas, porque sigue detectando una coordenada que se desplace, se pierda o se redondee.

Modelo. Un talud de la batería de tests cuyo contorno externo tiene ocho vértices, con dos contornos de material y un nivel freático de (0, 32) a (120, 22), exportado en metros e importado con una tolerancia de soldadura del 0.05 % y sin simplificación; el área se comprueba también tras una ida y vuelta en milímetros. Exigido. Los mismos tipos y números de contornos; una diferencia relativa del área encerrada inferior a 10⁻⁹; todos los vértices originales presentes, con una coincidencia de 10⁻⁶; todo vértice añadido a menos de 10⁻⁶ m de un segmento original; los contornos de material sin cambios hasta el noveno decimal; el carácter cerrado o abierto de cada contorno conservado; y las regiones aún construibles. Obtenido. Una diferencia de área de 0, en metros y en milímetros, y un contorno externo que pasa de 8 a 13 vértices, con los cinco vértices añadidos a una distancia de hasta 4×10⁻¹⁵ m del perfil original.

4.3.11Una presa rectangular: el caudal de Dupuit y Charnyi#

En el flujo a través de una presa rectangular de anchura $L$ sobre una base impermeable, con un nivel $H_1$ aguas arriba y un nivel $H_2$ aguas abajo, medidos sobre la base, el caudal por unidad de longitud es el que da la expresión siguiente (resultado de Charnyi, que coincide con el caudal de la aproximación de Dupuit), siempre que la zona situada por encima de la superficie libre no conduzca:

La formulación no saturada de OGR Slip2D permite que esa zona conduzca, por lo que cabe esperar que el caudal calculado resulte mayor. Al hacer más abrupta la función de permeabilidad, el caudal debe tender de forma monótona a la solución cerrada; el test comprueba esa convergencia, que valida el motor de cálculo y explica el exceso.

Modelo. Una presa de 20 m de anchura y 12 m de altura sobre una base impermeable, $K = 1\times10^{-5}$ m/s, $H_1 = 10$ m en la cara izquierda y $H_2 = 2$ m en la derecha; en la coronación y en las caras por encima de los niveles de agua se impone la condición de posible superficie de rezume; malla de 0.6 m; el caudal por la sección vertical del centro de la presa. La solución cerrada es $q = 10^{-5}(10^2 - 2^2)/(2\cdot 20) = 2.4\times10^{-5}$ m²/s. La función de permeabilidad es la de Gardner (1958), con un límite inferior $k_r = 10^{-8}$, en cuatro variantes cada vez más abruptas, que figuran en la Tabla 4.3.6.

Exigido. La diferencia disminuye a medida que la función se vuelve más abrupta (se tolera un aumento de hasta 10⁻³ entre casos consecutivos) y en el caso más abrupto queda por debajo del 5 %. El mismo archivo exige que la superficie libre salga por la cara de aguas abajo por encima del nivel de aguas abajo y por debajo del de aguas arriba, que parta cerca de la altura de aguas arriba, que descienda de forma monótona, con tolerancia para el ruido de la traza discreta, y que se den los nodos de la superficie de rezume. Mediante la capa de operaciones, la misma presa, con la función más abrupta, debe dar la solución cerrada con una diferencia inferior al 5 %.

4.3.12La curva de retención y su derivada#

La curva de retención de van Genuchten (1980) da el contenido de agua en función de la altura de succión $\psi$, y su derivada es la capacidad hídrica específica que aparece como término de almacenamiento de la zona no saturada (Ec. (3.7.15)):

Ambas expresiones valen para $\psi \gt 0$; donde el suelo está saturado, $\theta = \theta_s$ y $C = 0$. Modelo. $\theta_s = 0.45$, $\theta_r = 0.08$, $\alpha = 0.05$ m⁻¹ y $n = 2$. Exigido. $\theta = 0.45$ con una diferencia inferior a 10⁻¹² para succión nula y para presión intersticial positiva, y un valor inferior a 0.09 para una succión de 10⁶ m; y la capacidad analítica a menos de 10⁻³, en términos relativos, de la diferencia central de la curva (con un paso relativo de 10⁻⁴) en $\psi = 2,\ 10,\ 40$ y 150 m. Obtenido. Diferencias de 1.6×10⁻¹⁰, 2.0×10⁻⁹, 4.0×10⁻⁹ y 9.4×10⁻⁹, y $\theta = 0.08001$ a 10⁶ m. En el mismo archivo se comprueba que el agua almacenada es la integral del almacenamiento en ambas zonas, con una diferencia de hasta el 1 %, condición que mantiene la validez del esquema de Picard modificado a través del nivel freático.

4.3.13Desembalse rápido: el factor de seguridad en función del nivel y del tiempo#

Se comprueban dos propiedades del desembalse rápido, una frente a un resultado publicado y otra como requisito físico.

El desembalse total de un talud homogéneo. Para un talud homogéneo del tipo que recogen los ábacos de Morgenstern (1963), los valores de referencia que usa la batería de tests sitúan el factor de seguridad más bajo en el desembalse total, 1.20, frente a 1.41 con medio desembalse. El barrido de niveles de desembalse repite toda la búsqueda en cada nivel (Sección 3.8.6). Modelo. Un talud de 3 en horizontal por 1 en vertical, con $c'/\gamma H = 0.025$, $\phi' = 30^\circ$ y $\gamma/\gamma_w = 2$, comportamiento sin drenaje con $\bar B = 1$, y un embalse lleno hasta la coronación; método de Bishop, una rejilla de 3 × 4 intervalos, 18 dovelas y cuatro niveles desde la coronación hasta el pie. Exigido. El nivel crítico es el desembalse total; su factor de seguridad difiere menos del 6 % de 1.20; el factor en el nivel intermedio más próximo a medio desembalse es mayor; y el barrido de niveles no revela ningún margen inseguro superior al 1 %. Obtenido. 1.2103 en el desembalse total, un 0.86 % por encima de 1.20, y 1.304, 1.673 y 2.588 a un tercio, a dos tercios y a la altura total del embalse.

El factor de seguridad en el tiempo. Con el análisis transitorio, el factor de seguridad de las etapas marcadas con Calculate SF es el más bajo justo después del descenso del nivel, cuando las presiones aún no se han disipado, y se recupera a medida que lo hacen (Sección 3.8.7 y Figura 3.8.4). El test construye un desembalse sobre un talud de la batería con tres materiales, cuyas alturas piezométricas pasan de 45 m y 30 m a cada lado a 28 m y 26 m, con dos etapas marcadas a 10⁴ y 10⁶ unidades de tiempo, y exige que el factor de seguridad de la segunda etapa sea mayor que el de la primera. Es un requisito de orden y no una solución cerrada, y se cumple en la versión 0.1.235.

4.3.14Referencias#

• Darcy, H. (1856). Les fontaines publiques de la ville de Dijon. Victor Dalmont, Paris.
• 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.
• Hoek, E., Carranza-Torres, C. & Corkum, B. (2002). Hoek–Brown failure criterion – 2002 edition. Proceedings of the 5th North American Rock Mechanics Symposium and 17th Tunnelling Association of Canada Conference (NARMS-TAC 2002), Toronto, Vol. 1, 267–273.
• Morgenstern, N. (1963). Stability charts for earth slopes during rapid drawdown. Géotechnique, 13(2), 121–131.
• 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.
• Carslaw, H.S. & Jaeger, J.C. (1959). Conduction of Heat in Solids, 2nd ed. Oxford University Press, Oxford.