## Articulo

• Similares en SciELO

## versión On-line ISSN 1390-6542versión impresa ISSN 1390-9363

### Enfoque UTE vol.7 no.1 Quito ene./mar. 2016

#### https://doi.org/10.29019/enfoqueute.v7n1.87

Ingeniería General

Modelación y simulación numérica de la Ecuación de Richards para problemas de infiltración

Modeling and numerical simulation of the Richards equation for infiltration problems

Resumen:

Palabras clave: Infiltración; ecuación de Richards; modelo matemático; elementos finitos.

Abstract:

One of the most important natural resources we have is the soil and it is of great interest to the society to take care of them and not to pollute it. In the study of this issue, we are going to consider one of the most common forms of soil contamination due to an infiltration process. It is therefore that it is essential to address study and clearly understand this process by developing a mathematician model, which will be a representation of this physicist phenomenon. Then design and implement a computer program that simulates the infiltration of liquid pollutants in a given area. In this paper we will develop a mathematical model for two-dimensional infiltration in the saturated zone of porous media, based on the equation in nonlinear partial differential Richards Also, It will present a numerical solution through finite element method and first order This paper shows the computational implementation using a simulator that presents graphically the process of pollution afflicting the ground, exposed to certain pollutants, such as the oil spill in regions of eastern Ecuador, wastewater near industrial complexes, among others, over a certain period of time. Finally, this paper will allow for remedial studies in the case are already contaminated soils or preventive areas established as hazardous.

Keywords: infiltration; Richards equation; mathematical model; finite elements.

1. Introducción

A lo largo del tiempo, el Planeta se ha visto atacado por una serie de amenazas, como los contaminantes del suelo, que lo han llevado a su acelerado deterioro. Por ejemplo, la contaminación de suelos agrícolas debido a infiltraciones de agua contaminada, entre otros. Ante ello, surgen interrogantes como: ¿De qué manera se contaminan los suelos?, ¿Cuáles son los niveles de contaminación de los suelos a 10, 15 y 20 años?, ¿Qué medidas deberían adoptarse para evitar consecuencias negativas futuras? ¿Cuáles son los niveles de efectividad al aplicar medidas de prevención?

La búsqueda de respuestas conlleva al análisis de problemas producidos por infiltraciones de agentes contaminantes en el suelo. (Mulligan et al, 2004). La problemática planteada exige tomar medidas, una de ellas es el desarrollo de modelos matemáticos, modelos numéricos y programas computacionales que caractericen y simulen el fenómeno, analicen posibles escenarios y determinen cómo actúa y se desenvuelve el contaminante en el medio poroso. Todo ello en favor de garantizar la preservación de los recursos naturales y la vida del planeta (Gonzales de Vallejo, 2002).

El objetivo de este trabajo es desarrollar un modelo matemático que represente el fenómeno de infiltración de contaminante en el suelo basado en la ecuación de Richards no lineal, su modelo numérico del método de los elementos finitos y su implementación en un programa computacional que permita realizar simulaciones en diferentes escenarios y facilite una mejor comprensión del proceso de infiltración. Para que sirva como una herramienta para la toma de decisiones en beneficio de nuestro mundo.

La Ley de Darcy

Darcy (1856), demostró que el caudal que atraviesa el permeámetro era proporcional a la sección A y a la diferencia de alturas dividido por el diferencial de longitud es decir:

(1)

Donde k es la constante de proporcionalidad. Además determinó que al utilizar otro tipo de arena (más gruesa o más fina, o a su vez la mezcla de estas) y combinando todas las variables, se volvía a cumplir la ecuación anterior, pero la constante de proporcionalidad lineal era distinta. Por tanto esta, constante era una característica propia del material y la llamo permeabilidad. Dado que el caudal esta en , la sección en , y y son longitudes, se tiene que las unidades de permeabilidad son las de la velocidad . Por otra parte tomando incrementos infinitesimales se deduce que, la Ley de Darcy se expresa de esta forma:

(2)

En donde:

es el cambio de altura respecto a la horizontal

caudal que circula por de sección

El signo menos se debe a que el caudal es una magnitud vectorial, cuya dirección es hacia los decrecientes; es decir, que o es negativo y, por tanto, el caudal será positivo.

Ley de Darcy-Buckingham

La Ley de Darcy se ha obtenido para un solo fluido, sin embargo los suelos no saturados en sus poros, presentan dos fluidos (agua y aire), lo cual conlleva a la no aplicabilidad de este principio para estos casos. El efecto de taponamiento de los poros en los suelos, son causados por las burbujas de aire, que retienen el líquido impidiendo su natural permeabilidad, esta es la razón principal para que un suelo parcialmente saturado tenga menor permeabilidad que otro totalmente saturado.

Buckingham (1907), propuso que el grado de saturación va en aumento en la misma medida en que las burbujas de aire disminuyen en el tiempo ya que son arrastradas por las corrientes de agua, lo cual permite concluir que la permeabilidad de un suelo parcialmente saturado aumenta significativamente con el paso del tiempo.

El incremento de la constante de proporcionalidad (permeabilidad) en suelos semi - saturados es mayor mientras la presión del líquido aumenta, ya que esto provoca una considerable disminución del volumen que ocupan las burbujas de aire. Luego en este tipo de suelo la permeabilidad K depende de h y la ley de Darcy se corrige de la siguiente forma:

(3)

La densidad de un flujo (líquido o gas) homogéneo puede depender de mucho factores tales como la temperatura y la presión a la que está sometido. Para lo líquidos, la densidad varía muy poco dentro de amplios límites de presión y temperatura.

Con fines prácticos en la contaminación ambiental, los líquidos que no están sujetos a una fuente contaminante serán considerados de densidad constante. Por el contrario, si el líquido está sujeto a una fuente contaminante cambiará de densidad. La densidad de los gases es muy sensible a los cambios de temperatura y presión, la densidad se define como: o siendo sus unidades de denominación o .

Sea un volumen de control que se mueve con una velocidad de . La densidad del fluido contenido en es , de donde e integrando sobre , tenemos:

(4)

Derivando con respecto al tiempo, nos da:

(5)

y por el teorema de transporte de Reynolds, obtenemos, (Richards,1931)

(6)

Sistema conservativo

Un sistema se dice conservativo si la masa del sistema no se crea, ni se destruye; es decir que la masa permanece constante durante todo el proceso, esto es para todo es decir,

(7)

Como la integral precedente es válida para cualquier volumen de control, se demuestra que

(8)

Ecuación de Richards

La ecuación de Richards se obtiene de reemplazar la ley de Darcy-Buckingham (3) en la ecuación de conservación de masa o ecuación de continuidad para un sistema conservativo sobre , (Darcy, 1856),

Así la ecuación quedaría de la forma:

(9)

En el presente trabajo a fin de darle mayor generalidad a nuestro estudio se considerará la siguiente ecuación del tipo parabólico no lineal:

Sean abiertos, acotado, de frontera lipschisiana a trozos, . Ponemos , , .

Consideraremos el problema siguiente:

(10)

Formulación débil

Según Sayas (2007), para la obtención de una formulación débil, procedemos del modo siguiente, multipliquemos a la ecuación diferencial del problema (10) por una función tal que , e integremos sobre , dando como resultado:

Por el teorema de derivación bajo el signo de integración, (Richards, 1931) se tiene:

Ahora por el teorema de Green -Gauss, que se puede encontrar en (Brezis, 1983);

Y debido a que , se tiene:

Reemplazando estos resultados se obtiene la siguiente ecuación, para todo tal que .

Con estas notaciones la formulación débil del problema parabólico se escribe:

(11)

Modelo numérico y computacional

Sea un subespacio de dimensión finita una base de

La formulación aproximada se establece del modo siguiente: para cada , hallar solución de:

(12)

Sean y

La derivada en se aproxima por diferencias finitas progresivas del modo siguiente:

Para el resto de la ecuación usamos el método theta con . Luego el problema consiste en hallar solución de

.

Para (el método se llama de Crank Nicholson), se obtiene

, de donde

Definimos:

, ,

Se tiene el siguiente problema:

Hallar solución de :

Aplicando el método de Newton se sigue:

y tal que

El problema discreto se expresa del modo siguiente:

(13)

De forma más general el algoritmo a implementar seria:

(14)

Algoritmo numérico

y N (Topología de la malla)

(Solución inicial)

(Solución inicial arbitraria)

(Tiempo de simulación)

(Número de divisiones para la discretización temporal)

(Número máximo de iteraciones)

(Precisión con la que deseamos la solución)

PASOS:

PARA

REPETIR

Calcular la matriz con información de y

Calcular el vector con información de y

Resolver el sistema

HASTA QUE o

FIN DEL PARA

SALIDA: Solución del problema

Pruebas numéricas

En esta sección se va a tomar un ejemplo particular de la ecuación de Richards. Es decir, se fijan todas las funciones involucradas en el problema. De tal manera que la función cumpla con las condiciones de frontera y condiciones iniciales y se reemplazan en la ecuación para obtener la función ,con lo cual se ha construido un ejemplo particular de ecuación de Richards cuya solución es conocida, lo que nos permite validar el programa computacional y se obtiene algunos resultados numéricos.

Adicionalmente, se deja constante el tiempo y variamos el número de elementos (triángulos) de la discretización espacial, se calcula el error con la norma de para lo cual se integran en valor absoluto la diferencia de la solución exacta y aproximada sobre cada triángulo y luego se suman., Con lo cual se obtiene una muy buena aproximación de esta norma. Así mismo se registra el tiempo que se demora la computadora en realizar los cálculos, y finalmente se registran los resultados, como lo muestra la Tabla 1.

Como se puede apreciar, mientras aumentan el número de elementos en el mallado, el error disminuye de manera considerable, aunque el tiempo crece de manera exponencial, Sin embargo, este inconveniente se puede solucionar según las características del computador, en este caso la simulación se realizó en una laptop core i3.

Como segunda prueba, se analiza la convergencia del método de Newton en la solución de la ecuación de Richards y se evalúa la diferencia entre la solución numérica en el k-ésimo paso y en el k-1 ésimo paso, los resultados se verifican en la Tabla 2.

Tabla 2 Convergencia del método de Newton

En esta sección se encontrarán algunos resultados importantes sobre todo del simulador, corridas del programa con una comparación entre las soluciones exactas versus las aproximadas así como también graficas de en 2D y 3D del proceso de infiltración.

Se presentan los resultados del programa computacional que resuelve la ecuación de Richards bidimensional, se exponen las soluciones aproximadas y exactas en los nodos de la malla para lo cual se ha usado un problema cuya solución se conoce.

En la Figura 1 se presentan los resultados de la solución exacta y la aproximada; donde se nota que los resultados son muy similares, así validamos tanto el modelo numérico y su implementación.

El simulador tiene dos maneras de visualizar los resultados, la primera, en la Figura 2, muestra el resultado en tres dimensiones es decir, usa la malla para graficar la función resultante. Lógicamente como es un problema evolutivo dicha función va cambiando en el tiempo, lo cual nos una buena idea de cómo evoluciona el fenómeno.

La segunda forma de visualización es en dos dimensiones es decir a cada posición (x, y) del plano le corresponde un código de colores que varía de amarillo que significa 0 contaminación, hasta rojo intenso, que es donde la función toma los valores más altos. Lo que implica una mayor contaminación. En el caso del fenómeno de infiltración de agua en el suelo se puede observar en un corte transversal como el fluido va poco a poco ingresando en el suelo. En la Figura 3, se aprecia tres gráficas del simulador en tres tiempos diferentes; después de 1 hora la infiltración de agua en el medio poroso es leve, después de 5 horas aumenta, mientras que al cabo de 10 horas los niveles de infiltración son muy considerables.

3. Conclusiones y recomendaciones

Es posible demostrar la existencia y unicidad de soluciones de la ecuación de Richards tanto para el problema tipo elíptico coma para el problema tipo parabólico cuando las funciones involucradas cumplen con determinadas características.

Para resolver numéricamente la ecuación de Richards bidimensional usando el método de elementos finitos, se debe usar primeramente el método de Crak Nicolson en la variable temporal para lo cual se realiza una discretización en el tiempo. Luego se linealiza la ecuación resultante usando el método Newton. Finalmente se debe discretizar la variable espacial y de esta manera se consigue el esquema numérico que posteriormente se transforma en un programa computacional.

Para realizar pruebas con el simulador se usó un problema cuya solución se conocía de antemano, aquí se pudo determinar que dejando fijo el tiempo y variando el número de elementos de la malla, el error en norma disminuye considerablemente conforme aumenta el número de triángulos de la discretización espacial. Pero, así mismo el tiempo empleado por el computador aumenta considerablemente al aumentar el número de elementos.

Como era de esperarse la discretización de la variable tiempo también influye en la calidad de la solución obtenida, es decir si el paso temporal es demasiado fino, el número de operaciones elementales crece notablemente, para lo cual se necesitaría un computador con mayor capacidad de procesamiento, caso contrario el tiempo de obtención de resultados es demasiado grande.

La infiltración en un medio poroso de cualquier líquido contaminante depende básicamente de dos variables, el tipo de suelo y el tiempo de exposición de éste al contaminante.

A partir de los resultados del simulador, se recomienda realizar pruebas con datos experimentales en diferentes tipos de suelo y con diferentes tipos de contaminantes con el objetivo de contrastar resultados y de ser necesario realizar una calibración del simulador.

En base a los resultados que se obtiene en y se recomienda formular un modelo tridemensional y realizar el estudio de la existencia y unicidad de la solución de la ecuación de Richards.

Bibliografía

Brezis, H. (1983). Analyse Fonctionnelle. Masson, Paris. [ Links ]

Darcy, H. (1856). Les fontaines publiques de la ville de Dijon. Dalmont, Paris. [ Links ]

Buckingham, E. (1907). Studies on the movement of soil moisture. Washington: USDA [ Links ]

Evans, L.C. (1998). Partial differential Equations. American Mathematical Society. Providence: Graduate Studies in Mathematics. [ Links ]

Blytth, F., Freitas, M (2001). Geología Para Ingenieros. México: Compañía Editorial Continental. [ Links ]

Ehlers, W. Deformation and localization analysis of partially saturated soil, Recuperado 23 de septiembre de 2015. http://www.sciencedirect.com/science/article/pii/S0045782504001148Links ]

Gómez, J., Martin, D. (2010). Análisis Funcional y Optimización. Chile: Universidad de Chile. [ Links ]

Gonzales de Vallejo, L. I., & Ferrer, M. (2002). Ingeniería Geológica. Madrid: Pearson Educacion. [ Links ]

Pino, E., Mejía, J., Abel, E. (Enero, 2012). Modelamiento numérico espacio-temporal 1d de la infiltración basado en la ecuación de Richards y otras simplificadas. Eciperu, 9, 31-36. [ Links ]

Richards, L. (1931). Capillary conduction of liquids in porous media. USA: Physics, v. 1. [ Links ]

Richards, L.A., Gardner, W.R., Ogata, G. (1956). Physical processes determining water loss from soil. Soil Science Society of American Proceeding 20, 310-314. [ Links ]

Sayas, F. (2007). Modelos Matemáticos en Mecánica. España: Departamento de Matemática Aplicada, Universidad de Zaragoza. [ Links ]

Recibido: 30 de Octubre de 2015; Aprobado: 29 de Marzo de 2016

Este es un artículo publicado en acceso abierto bajo una licencia Creative Commons