Un método de elementos finitos usando una variable auxiliar acotada para resolver la ecuación de Richards

Información Básica

• ID del Artículo: 2510.13012
• Título: Un método de elementos finitos usando una variable auxiliar acotada para resolver la ecuación de Richards
• Autores: Abderrahmane Benfanich (Universidad de Ottawa), Yves Bourgault (Universidad de Ottawa), Abdelaziz Beljadid (Universidad de Ottawa y Universidad Politécnica Mohammed VI)
• Clasificación: math.NA cs.NA physics.comp-ph
• Fecha de Publicación: 14 de octubre de 2025 (preimpresión arXiv)
• Enlace del Artículo: https://arxiv.org/abs/2510.13012

Resumen

La ecuación de Richards es una ecuación no lineal elíptico-parabólica ampliamente utilizada en la modelación de percolación en medios porosos. En este artículo se desarrolla un método de elementos finitos para resolver la ecuación de Richards mediante la introducción de una nueva variable auxiliar acotada que elimina términos no acotados en la forma débil. El método utiliza un esquema semi-implícito para la discretización temporal y el método de Newton para resolver sistemas no lineales. El método elimina la necesidad de técnicas de regularización y presenta ventajas en el tratamiento de regiones completamente secas y saturadas. En la técnica propuesta, se utiliza el método de descomposición de dominios de Schwarz sin solapamiento para modelar la percolación en suelos estratificados. Los experimentos numéricos validan la efectividad del método propuesto, incluyendo la modelación de flujo en fibras completamente secas, casos con regiones completamente saturadas y secas, así como problemas de percolación en suelos estratificados.

Antecedentes de Investigación y Motivación

Definición del Problema

La ecuación de Richards describe el flujo no saturado en medios porosos y es un modelo fundamental en hidrogeología, ingeniería ambiental y ciencia del suelo. La ecuación existe en múltiples formas:

1. Forma de carga de presión (Ψ): Aplicable a flujo completamente saturado y parcialmente no saturado, pero numéricamente inestable cuando S→0
2. Forma de saturación (S): Mejor desempeño en suelos secos, pero dificultades numéricas en regiones saturadas (S=1)
3. Forma mixta (S,Ψ): Mejora las propiedades de conservación de masa, pero aumenta la complejidad

Limitaciones de Métodos Existentes

1. Métodos de cambio de variable principal: Pueden producir soluciones no físicas en la interfaz saturada/no saturada
2. Técnicas de regularización: Requieren selección adicional de parámetros de modelado
3. Métodos de elementos finitos tradicionales: Presentan dificultades numéricas al tratar regiones completamente secas (S=0) y completamente saturadas (S=1)

Motivación de la Investigación

Desarrollar un método numérico unificado que pueda:

• Evitar cambios de variables y regularización
• Manejar establemente regiones secas y saturadas
• Mantener buenas propiedades de conservación de masa
• Ser aplicable a medios porosos heterogéneos

Contribuciones Principales

1. Propuesta de nueva forma u: Introducción de una variable auxiliar acotada u que elimina términos no acotados en la forma débil
2. Desarrollo de esquema numérico estable: Combinación de discretización temporal semi-implícita e iteración de Newton
3. Implementación de método de descomposición de dominios: Tratamiento de problemas en medios heterogéneos como suelos estratificados
4. Validación de efectividad del método: Demostración de estabilidad y precisión a través de múltiples experimentos numéricos

Detalle de la Metodología

Definición de la Tarea

Resolver la ecuación de Richards: $\frac{\partial\theta}{\partial t} + \nabla \cdot q = 0$

donde θ es el contenido de agua y q es el flujo de agua según la ley de Darcy-Buckingham: $q = -K_s(x)K_r(x,S)\nabla(\Psi(x,S) + z)$

Construcción de la Nueva Variable u

Para todos los modelos, la derivada de la función de Leverett puede escribirse como: $J'(S) = C(x)S^{-a}(1-S^c)^{-b}$

Se define la variable auxiliar: $u = \int_0^S (1-s^c)^{-b} ds$

Expresada mediante la función Beta incompleta: $u = \frac{1}{c}B(S^c, \frac{1}{c}, 1-b)$

Ecuación de Richards en Forma u

En medios homogéneos, la forma u es: $\phi \frac{\partial S(u)}{\partial t} - \nabla \cdot \left(K_sK_r(S(u))\left(h_{cap}CS(u)^{-a}\nabla u + e_z\right)\right) = 0$

Forma Variacional

Se busca u ∈ L²(I;H¹(Ω)) tal que: $\int_\Omega \phi \frac{\partial S(u)}{\partial t} v dx + \int_\Omega K_sK_r(S(u))\left(h_{cap}S(u)^{-a}\nabla u + e_z\right) \cdot \nabla v dx + \text{términos de frontera} = 0$

Discretización Temporal

Se utiliza esquema de Euler semi-implícito: $\int_\Omega \phi \frac{S(u_h^{n+1}) - S(u_h^n)}{\Delta t} v_h dx + \int_\Omega K_sK_r(S(u_h^n))\left(h_{cap}S(u_h^n)^{-a}\nabla u_h^{n+1} + e_z\right) \cdot \nabla v_h dx = 0$

Iteración de Newton

Para términos no lineales S(u^{n+1,m+1}), se utiliza expansión de Taylor de primer orden: $S(u^{n+1,m+1}) = S(u^{n+1,m}) + \frac{\partial S}{\partial u}(u^{n+1,m})(u^{n+1,m+1} - u^{n+1,m}) + O((u^{n+1,m+1} - u^{n+1,m})^2)$

Método de Descomposición de Dominios

Para suelos estratificados, se utiliza el método de Schwarz sin solapamiento:

• Resolución independiente de la ecuación de Richards en cada subdominio
• Acoplamiento mediante condiciones de interfaz: continuidad de carga de presión y conservación de flujo
• Resolución iterativa mediante condiciones de transmisión tipo Robin

Configuración Experimental

Implementación Numérica

• Plataforma de Software: FreeFEM++ y FENICSx
• Solucionadores Lineales: UMFPACK, PETSc, MUMPS
• Malla: Malla triangular uniforme, elementos finitos lineales (k=1)

Indicadores de Evaluación

Error L²: $E_{h,\Delta t}(T)^2 = \|u(\cdot,T) - u_{h,\Delta t}(\cdot,T)\|_{L^2(\Omega)}^2$

Orden de convergencia:

• Orden de convergencia espacial: $p = \log_2(E_{h,\Delta t}(T)/E_{h/2,\Delta t}(T))$
• Orden de convergencia temporal: $q = \log_2(E_{h,\Delta t}(T)/E_{h,\Delta t/2}(T))$

Casos de Prueba

1. Fibra 1D completamente no saturada: Verificación de estabilidad en condiciones secas
2. Región mixta 1D saturada-no saturada: Prueba de capacidad para manejar valores iniciales discontinuos
3. Verificación de solución fabricada: Cálculo de órdenes de convergencia
4. Flujo 2D completamente saturado: Comparación con software Hydrus
5. Problema de suelo estratificado: Validación del método de descomposición de dominios
6. Inclusión en medio heterogéneo: Prueba de manejo de geometrías complejas

Resultados Experimentales

Resultados Principales

Prueba de Fibra 1D

• Malla: 500 puntos de malla uniforme
• Paso de tiempo: Δt = 1s
• Resultados: Forma u y forma S completamente consistentes, buen acuerdo con datos experimentales

Región Mixta Saturada-No Saturada

• Malla: 5000 puntos de malla uniforme
• Paso de tiempo: Δt = 10⁻⁵ horas
• Convergencia de Newton: Promedio 1-5 iteraciones
• Estabilidad: Solución mantiene monotonicidad, sin oscilaciones numéricas

Análisis de Convergencia de Solución Fabricada

Convergencia Espacial (Δt fijo = 5×10⁻⁵):

Convergencia Temporal (h fijo = 1×10⁻³):

Comparación de Flujo 2D Saturado

• Malla: 95,142 elementos triangulares, 47,972 vértices
• Paso de tiempo: Δt = 0.001 horas
• Resultados: Altamente consistentes con resultados del software Hydrus

Prueba de Suelo Estratificado

• Malla: 89,306 elementos triangulares, 45,064 vértices
• Paso de tiempo: Δt = 0.05 horas
• Convergencia: Tolerancia 10⁻³, máximo 10 iteraciones, λ = 25
• Resultados: Consistentes con resultados en la literatura 20

Comparación de Métodos

En la prueba de inclusión en medio heterogéneo:

• Forma u: Convergencia estable incluso cuando ε=0, promedio 4 iteraciones
• Forma Ψ: Requiere >1000 iteraciones cuando ε=10⁻⁵, tolerancia residual aún 0.04

Trabajo Relacionado

Direcciones Principales de Investigación

1. Métodos de cambio de variables: Diersch & Perrochet (1999), Forsyth et al. (1995)
2. Técnicas de regularización: Schweizer (2007), Pop & Schweizer (2011)
3. Formas mixtas: Celia et al. (1990), Maina & Ackerer (2017)
4. Descomposición de dominios: Lions (1990), Manzini & Ferraris (2004)

Ventajas del Presente Trabajo

1. Tratamiento unificado: Sin necesidad de cambio de variables o regularización
2. Estabilidad numérica: Estable en todo el rango S∈0,1
3. Consistencia física: Mantiene continuidad de carga de presión
4. Eficiencia computacional: Convergencia rápida de Newton

Conclusiones y Discusión

Conclusiones Principales

1. Validez de forma u: Eliminación exitosa de términos no acotados en formas tradicionales
2. Estabilidad numérica: Mantiene estabilidad en regiones completamente secas y saturadas
3. Desempeño de convergencia: Convergencia de segundo orden espacial, primer orden temporal
4. Practicidad: Aplicable a problemas complejos en medios heterogéneos

Limitaciones

1. Dependencia del modelo: Requiere satisfacer la condición $\lim_{S\to 0} K(S)S^{-a} < \infty$
2. Resolución no lineal: Requiere iteración de Newton en regiones saturadas
3. Tratamiento de interfaz: Requiere técnicas de descomposición de dominios en medios heterogéneos
4. Sensibilidad de parámetros: El parámetro Robin λ en descomposición de dominios requiere ajuste

Direcciones Futuras

1. Análisis teórico: Análisis numérico riguroso y pruebas de convergencia
2. Acoplamiento multifísico: Extensión a transporte de solutos y transferencia de calor
3. Técnicas adaptativas: Desarrollo de malla adaptativa y pasos de tiempo adaptativos
4. Computación paralela: Implementación paralela a gran escala

Evaluación Profunda

Fortalezas

1. Fuerte innovación: Construcción ingeniosa de variable u con base teórica sólida
2. Estabilidad numérica: Evita dificultades numéricas de métodos tradicionales
3. Experimentos exhaustivos: Cubre casos desde 1D a 2D, de medios homogéneos a heterogéneos
4. Practicidad ingenieril: Comparación con software comercial Hydrus aumenta credibilidad

Insuficiencias

1. Análisis teórico limitado: Carencia de pruebas rigurosas de convergencia y estabilidad
2. Análisis de complejidad computacional: Análisis insuficiente del costo computacional de iteración de Newton
3. Selección de parámetros: Falta de orientación teórica para selección del parámetro Robin λ
4. Validación 3D: Carencia de validación en problemas tridimensionales

Impacto

1. Valor académico: Proporciona nuevas perspectivas para resolución numérica de ecuación de Richards
2. Valor práctico: Amplias perspectivas de aplicación en hidrogeología e ingeniería de suelos
3. Reproducibilidad: Proporciona detalles de implementación detallados, facilitando reproducción

Escenarios de Aplicación

1. Flujo en medios porosos: Flujo de agua subterránea, movimiento de agua en suelo
2. Aplicaciones ingenieriles: Percolación en presas, análisis de estabilidad de taludes
3. Problemas ambientales: Migración de contaminantes, remediación de agua subterránea
4. Ciencia de materiales: Percolación de fluidos en materiales fibrosos

Referencias

Referencias Clave

1. Richards, L.A. (1931). Capillary conduction of liquids through porous mediums. Physics 1(5), 318-333.
2. Celia, M.A., et al. (1990). A general mass-conservative numerical solution for the unsaturated flow equation. Water Resour. Res. 26(7), 1483-1496.
3. van Genuchten, M.T. (1980). A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44(5), 892-898.
4. Manzini, G., Ferraris, S. (2004). Mass-conservative finite volume methods on 2-d unstructured grids for the Richards' equation. Adv. Water Resour. 27, 1199-1215.

Evaluación General: Este es un artículo con importante significado innovador en el campo de la resolución numérica de la ecuación de Richards. Mediante la introducción de una variable auxiliar acotada u, resuelve ingeniosamente las dificultades numéricas de métodos tradicionales al tratar regiones secas y saturadas. El método posee base teórica sólida y valor práctico significativo, con experimentos numéricos que validan completamente la efectividad del método. Aunque el análisis teórico requiere fortalecimiento, en general se trata de un trabajo de investigación de alta calidad.