TEMA PROBLEMA TEMA: Modelación y simulación de procesos de infiltración de agua en los suelos en zonas no saturadas Trabajo de graduación Manuel Jesús de los Ángeles Guamán Jachero PROBLEMÁTICA: La contaminación del suelo debido a la infiltración de líquidos contaminados en medios porosos es un problema grave y si es que no se toman medidas adecuadas, el suelo seguirá contaminando la producción de alimentos, afectando al ser humano, a las plantas y a los animales. De aquí se origina la necesidad de desarrollar un modelo matemático para prevenir o disminuir la contaminación. La infiltración del agua en zonas porosas – la ingeniería, arquitectura. El modelo matemático es la Ecuación de Richards unidimensional de flujo de fluidos en medios porosos. Otro problema es que, la Ecuación de Richards no tiene solución general analítica, por lo que es necesario usar aproximaciones numéricas. También es necesario que se desarrolle la discretización de la ecuación de Richards usando el método de volúmenes finitos. OBJETIVOS OBJETIVO GENERAL: Desarrollar un modelo matemático y construir un esquema numérico convergente y estable, para la simulación de procesos de infiltración de agua en los suelos en zonas no saturadas, mediante el uso de un esquema implícito en el tiempo, el método de Volúmenes Finitos en el espacio y dos esquemas de linealización diferentes. OBJETIVOS ESPECÍFICOS Elaborar el modelo matemático para el proceso de infiltración del agua en el suelo saturado y no saturado. Discretizar la ecuación de Richards unidimensional- - el método de volúmenes finitos - para aproximar su solución mediante sistemas de ecuaciones no lineales. Linealizar la ecuación de Richards - método de Picard y L esquemas. Construir un algoritmo numérico - método de Picard y L esquemas - elaborar un programa computacional para la simulación numérica. Mostrar la estabilidad y convergencia del esquema numérico - el método de volúmenes finitos para la ecuación de Richard unidimensional – el uso de experimentos numéricos. 4 ECUACION DE RICHARDS Ley de Darcy. es el gasto A es el área transversal el gradiente de presión K conductividad hidráulica L distancia En términos de potencial mátrico con coordenada vertical que crece hacia arriba, se expresa unidimensionalmente: Ley de Darcy-Buckingham. Buckingham se dio cuenta de que la ley de Darcy podía ser extendida para describir el flujo en un medio poroso no saturado. Donde la conductividad hidráulica K es una función del contenido de humedad, , y de la succión, Ecuación de continuidad: Ecuación de Richards. Dónde: C (h) es la capacidad hidráulica específica del suelo q es el caudal por unidad de área, h es el potencial mátrico debido a la succión, t es el tiempo, z es la profundidad, K(h) es la conductividad hidráulica, es la humedad volumétrica del suelo, Para demostrar la existencia y unicidad de la ecuación de Richards, se utilizo el articulo “quasilinear heliptic - parabolic differential ecuations” 5 MÉTODO DE VOLUMENES FINITOS 6 METODO DE VOLÚMENES FINITOS El contenido fue tomado de [Cendón, 2008],[Le Veque, 2000]. Este método tiene sus raíces en la dinámica de fluidos. Después, se convirtió en el método de elección para el análisis de problemas relacionados con el flujo de fluidos y otros fenómenos de transporte. El método de volúmenes finitos y el MDF tienen sus raíces iniciales en la dinámica de fluidos computacional (CFD) y los fenómenos de transporte. Los métodos de volúmenes finitos provienen inicialmente de la forma integral de la ley de conservación La forma general de un sistema de leyes de conservación es Este método es también utilizado para la discretización de la ecuación de Richards en forma temporal y espacial 7 Descripción matemática del método de volúmenes finitos Para el desarrollo de un algoritmo de volúmenes finitos para resolver un problema definido mediante ecuaciones diferenciales y condiciones de contorno se tienen los siguientes pasos: El problema debe reformularse en forma variacional. El dominio de variables independientes (dominio espacial para problemas dependientes del tiempo) debe dividirse mediante una partición en subdominios, llamados celdas o volúmenes finitos. Asociada a la partición anterior se construye un espacio vectorial de dimensión finita, llamado espacio de volúmenes finitos. Se plantea los dos tipos de esquemas a resolver luego de haber discretizado la ecuación de Richard para la forma mixta Unidimensional. El esquema explicito calcula el estado del sistema en un momento posterior al estado actual del sistema. El esquema implícito son la solución por la solución de una ecuación que ambos estados de sistema actuales y posteriores. 4. Esto da lugar a un sistema con un número finito de ecuaciones, con un número elevado de incógnitas. El número de incógnitas indica la dimensión del espacio vectorial de volúmenes finitos obtenido, cuanto mayor sea dicha dimensión mejor será la aproximación numérica obtenida. 5. El último paso es el cálculo numérico de la solución del sistema de ecuaciones. Los pasos anteriores permiten construir un problema de cálculo diferencial en un problema de álgebra lineal. Para la solución se pueden utilizar métodos de algebra lineal o métodos iterativos. 8 DISCRETIZACION DE LA ECUACION DE RICHARDS 9 DISCRETIZACION DE LA ECUACION DE RICHARDS DISCRETIZACION ESPACIAL DE LA ECUACION DE RICHARDS MEDIANTE VOLUMNES FINITOS Se construye un volumen de control entre en la ecuación de Richard (5,3) MATRIZ OBTENIDA CON EL METODO DE VOLUMENES FINITOS 10 Discretización del dominio Sea la dimensión del de . Dado un dominio con una frontera continua en el sentido de Lipschitz una partición de en, es una colección de subdominios que satisface: 1. 2. un conjunto compacto con frontera de Lipschitz continúa 3. Por sencillez del análisis, todos los tienen la misma forma, es decir su dominio de referencia son rectas subdivididas por igual longitud. Además, sobre cada elemento de la celda se considerarán algunos puntos especiales, llamados nodos y que generalmente incluirán los vértices del volumen finito (celda) y adicional se requerirá la condición de que elementos adyacentes comparten los nodos o aristas. Dominio N volúmenes finitos Subdominios Nodos Celdas Volúmenes finitos 11 Esquema de discretización numérica. Para Discretizar la ecuación de Richards se pueden realizar en dos esquemas: Como se puede ver en la figura, para el esquema explicito utiliza solo los valores en los pasos espaciales , en el paso temporal anterior, , para encontrar el valor de la función en el paso espacial y temporal actual, . En el esquema implícito, los valores espaciales en el paso temporal actual también se tienen en cuenta. Por lo tanto, la función en depende de , además de los valores en el paso temporal anterior. La función que se va a calcular en el esquema implícito depende, por lo tanto, de sí misma y de otras variables desconocidas t t t t 12 Discretización espacial y temporal de la ecuación de Richards unidimensional La ecuación unidimensional del flujo cuando es vertical el modelo es el siguiente: Discretización del tiempo Primero discretizamos la variable temporal para lo cual: Fijamos que es el número de divisiones temporales y definimos con lo cual se tiene que Luego la discretización temporal es Así la discretización en el tiempo se muestra en el siguiente gráfico T T T 0 13 ESQUEMA EXPLICITO Usando ésta malla temporal en hacemos una aproximación mediante diferencias finitas centrales la derivada temporal donde Así, nuestro problema discreto en la variable temporal nos queda multiplicando por , simplificando y por transposición términos se tiene: Discretización en la variable temporal 14 LINEALIZACION APLICAMOS EL MÉTODO DE PICARD MODIFICADO 15 Hallar tal que, sobre Para la solución de la no linealidad en nuestro problema, aplicamos el método de Picard Modificado, (5,1) El contenido de humedad se aplicará con el método de Taylor Reemplazando, en la ecuación anterior, (5,2) Reemplazando la ecuación (5,2) anterior en la ecuación (5,1) se tiene, (copiada la (5,1), para reemplazar mas fácil) Organizando términos se tiene: (5,3) 16 Discretización del espacial (i) En este apartado desarrollamos la discretización espacial, en donde: Fijamos la variable N y definimos El mallado de los volúmenes finitos es puntualizado en los siguientes conjuntos de nodos, donde, : N (5,3) 17 Integramos sobre el volumen de control entre en (5,3) Aplicando el teorema fundamental del cálculo y el teorema de la divergencia de Gauss. (5.5) Notación de variables: ; ; ; ; Se ha aplicado el método de diferencias finitas centrales, Reemplazamos en (5.5) 18 Aplicando derivadas parciales dentro del corchete y otros procesos, se tiene Multiplicado en primer corchete tenemos: Organizando los términos obtenemos: Se cambió con El corchete se hizo 4 términos Factor común Factor común 19 Condiciones de frontera Condiciones de frontera donde, ; Resolvemos las operaciones en los exponentes Para utilizar en el algoritmo computacional tomamos de la primera y segunda fila de la ecuación de arriba: 20 () Resolvemos las operaciones en los exponentes, 21 Media aritmética, geométrica y armónica Para evaluar en la ecuación anterior se puede utilizar cualquiera de las siguientes formas: Media Aritmética   Media Geométrica   Media Armónica   22 Algoritmo Computacional 23 Linealización de la ecuación de Richards mediante el L-esquema 24 Linealización de la ecuación de Richards mediante el L-esquema El método fue propuesto para la ecuación de Richards por Pop, I.S., Radu, F.A.,Knabner, P. Es un método iterativo, robusto y linealmente convergente Es el único método que explota la monotonía del contenido de humedad, definido como θ(·) Además, la tasa de convergencia no depende del tamaño de malla. Los sistemas lineales que surgen después de usar el L-esquema están mucho mejor condicionados que el sistema para Newton o métodos modificados de Picard. Sea , ser dado. Busquemos , de modo que tenemos la siguiente ecuación, Se mantiene para todos los valores de . Para garantizar la convergencia del esquema, la constante debe satisfacer . Explicaremos mas detalladamente enseguida. El elemento clave en el esquema es la adición de un termino de estabilización , que junto con la monotonía del contenido de humedad (.), asegurará la convergencia del L-esquema. 25 Contenido de humedad y su derivada La curva característica de humedad del suelo es la relación entre la humedad volumétrica θ y la cabeza de presión. Esta curva describe el estado de energía relativo al volumen de agua almacenada en un material poroso bajo condiciones variables de saturación, ver en (Stephens, D. B. Vadose zone hydrology, 1995). Conductividad Hidráulica (K) o permeabilidad Representa la mayor o menor facilidad con que un medio poroso permite el paso de un fluido a través de el por unidad de áreas transversal a la dirección del flujo. Puede expresarse bajo condiciones saturadas y no saturadas. En la zona no saturada este valor no es único ya que depende del contenido de humedad del suelo y también de la presión existente en el sitio de contacto del agua con el suelo. Curva del contenido de humedad θ(Ψ) Esta curva permite describir la capacidad que tiene el suelo de almacenar o libera agua, corresponde a la relación entre el volumen volumétrico de agua presente en el suelo y la succión presente en la matriz del suelo ( Rawls,1993). La siguiente expresión, muestra el modelo propuesto por Van Genuchten para θ(Ψ), Cuando el medio esta saturado. = Contenido volumétrico saturado de agua = Contenido volumétrico residual de agua m = 1 − 1/n, donde n es un parámetro empírico α = Factor de escala. 26 OBSERVACIONES: Obtenemos que el L-esquema es ligeramente diferente del considerado por el autor Slodicka, M.: “Un esquema de linealización robusto y eficiente para problemas parabólicos doblemente no lineales y degenerados que surgen en el flujo en medios porosos ”. Además, los L-esquemas están considerando con la transformación de Kirchhoff, que no es el curso en el presente documento. En el caso de una constante ”K”, los métodos de Newton y el modificado Picard coinciden. Además, si L es reemplazado por la matriz Jacobiana se vuelve a obtener el esquema de Picard modificado. Cualquier de los métodos de linealización presentados anteriormente conduce a un sistema de ecuaciones lineales para (más precisamente, lo desconocido será el vector con los componentes de en una base de Vh ). Los derivados del contenido de agua y conductividad hidráulica en el caso del esquema de Picard modificado y el método de Newton se pueden calcular analíticamente o mediante un enfoque de perturbación como lo sugiere, ver en (Forsyth, P.A., Wu, Y.S., Pruess, K., 1995,p´ ag: 25-38) y las integrales que se producen se van aproximando mediante una fórmula en cuadratura. Para detener las iteraciones, adoptamos un criterio general de convergencia dado por lo siguiente vector Norma euclidiana >0 constante constante para resultados mas óptimos 27 Convergencia de L-esquema Esta información fue obtenido de (Florian List, 2016). Denotamos por (5.12) El error en la iteración j. Un esquema es convergente si en, , cuando . Las siguientes proposiciones son afirmaciones sobre las funciones de coeficientes y la solución discreta están definiendo en el marco que podemos probar la convergencia del esquema L. A1 La función θ(·) es monotonía, estrictamente creciente y Lipschitz continua. A2 La función K(·) es Lipschitz continua y dos constantes tal que A3 La solución del problema (4) satisface denotando con norma Ahora podemos afirmar el resultado teórico central de este tema. Teorema 5.1. Sea n ∈ 1, ..., N y supongamos que (A1)y (A3) son verdaderos. Si la constante L y el paso de tiempo se eligen de modo que (P16) a continuación se cumpla, el esquema L (3.6) converge linealmente, con una tasa de convergencia dada por verdaderos, La demostración la encontramos en (Florian List, 2016). 28 Hallar tal que, sobre Para la solución de la no linealidad en nuestro problema, aplicamos el método L-esquemas, (5,1) El contenido de humedad se aplicará con el método de Taylor Reemplazando, en la ecuación anterior, (5,2) Reemplazando la ecuación (5,2) anterior en la ecuación (5,1) se tiene, (copiada la (5,1), para reemplazar mas fácil) Organizando términos se tiene: (5,3) 29 Discretización espacial y temporal de la ecuación de Richards unidimensional por L esquemas Discretización del tiempo y su grafica Las condiciones de frontera Media aritmetica, geométrica y armónica, El contenido de los cuatro ítems anteriores, es casi parecido al contenido que en método de Picard, el ingreso del elemento clave en el esquema es la adición de un termino de estabilización , que junto con la monotonía del contenido de humedad (.), asegurará la convergencia del L-esquema. 30 Discretización del espacial (i) En este apartado desarrollamos la discretización espacial, en donde: Fijamos la variable N y definimos El mallado de los volúmenes finitos es puntualizado en los siguientes conjuntos de nodos, donde, : N (5,3) 31 Integramos sobre el volumen de control entre en (5,3) Aplicando el teorema fundamental del cálculo y el teorema de la divergencia de Gauss. Notación de variables: ; ; ; Se ha aplicado el método de diferencias finitas centrales, 32 Aplicando derivadas parciales dentro del corchete y otros procesos, se tiene Multiplicado en primer corchete tenemos: Organizando los términos obtenemos: Se cambió con El corchete se hizo 4 términos 33 Condiciones de frontera Condiciones de frontera donde, ; Resolvemos las operaciones en los exponentes Para utilizar en el algoritmo computacional tomamos de la primera y segunda fila de la ecuación de arriba: 34 () Resolvemos las operaciones en los exponentes, 35 Media aritmética, geométrica y armónica Para evaluar en se puede utilizar cualquiera de las siguientes formas: Media Aritmética   Media Geométrica   Media Armónica   36 Algoritmo Computacional ; 37 Programa computacional y resultados numéricos Validación del algoritmo numérico Cuadros de resultados RE es error relativo RSE es error raíz cuadrada El error es muy pequeño, la función aproximada es valida El algoritmo es eficiente para el programa computacional 38 RE es error relativo Las graficas casi coinciden El error es muy pequeño, la función aproximada es valida El algoritmo es eficiente para el programa computacional 39 Validación Del Programa Computacional ANALISIS COMPARATIVO: El método de Picard converge mas rápidamente, en menor tiempo que el de L esquemas El método de Picard al crecer la malla el numero de iteraciones baja, mientras que en el de L esquemas, el numero de iteraciones de constante 40 Convergencia con el método de Picard: Convergencia con el método de L-esquemas: GRAFICA DE LA SOLUCION. Los dos métodos llegan a la misma solución El método de L esquemas converge con un mallado mas bajo, de forma mas rápida con 30 de mallado, mientras que el de Picard recién converge con un mallado de 100 41 CONCLUSIONES: Para resolver numéricamente la ecuación de Richards unidimensional usando el método de volúmenes finitos, Se debe usar primeramente en la variable temporal una discretización en el tiempo, Luego se linealiza la ecuación resultante usando los métodos de Picard y de L-esquemas para finalmente discretizar la variable espacial y de esta manera se consigue el esquema numérico que posteriormente se transforma en un programa computacional. Se consiguió construir un algoritmo numérico para el método de Picard y L esquemas y se ha podido elaborar un programa computacional para la simulación numérica. Se puede demostrar la estabilidad y convergencia del esquema numérico mediante el método de volúmenes finitos para la ecuación de Richard unidimensional El método de Picard converge mas rápidamente, en menor tiempo que el de L esquemas. El método de Picard al crecer la malla el numero de iteraciones baja, mientras que en el de L esquemas, el numero de iteraciones es constante En las graficas, el método de L esquemas converge con un mallado mas bajo, mientras que el de Picard recién converge con un mallado mas alto El trabajo presentado no puede interpretarse como si fuese el diseño completo. Más bien, el objetivo es demostrar la forma en que los modelo matemáticos pueden ser puestos en práctica. 42 RECOMENDACIONES: Realizar un estudio más profundo en proceso de errores de computación de ver la posibilidad de debilitar aun más el rendimiento CPU de manera que la ecuación de Richards tenga solución única- con una mínimo error en optimar el resultado. Investigar la posibilidad reformular la ley de Darcy obtener nuevas formas de procesar el algoritmo de la ecuación de Richards. Diseñar otros algoritmos de programación paralelos se ejecuten al mismo tiempo muchas instrucciones, con lo cual se conseguirá minimizar el error en un tiempo razonable de maquina. Averiguar cual seria el número de elementos optimo una solución razonable menor tiempo, los números de iteraciones con respecto a la ecuación de Richards. Extender los alcances del modelo matemático y ajustarlo de ser necesario para que tenga en cuenta los efectos asociados al flujo de agua en medios porosos, por ejemplo, la anisotropía de los suelos en términos de resistencia, rigidez y la permeabilidad. A fin de acortar el tiempo de ejecución de la simulación se recomienda investigar sobre nuevos métodos numéricos de resolución de grandes sistemas de ecuaciones lineales y así como mejorar los existentes 43 BIBLIOGRAFIA: Abuja, G, (2013), Ecuación de Richards bidimensional, Modelación y simulación numérica, Trabajo de graduación, Universidad Central del Ecuador, Quito-Ecuador. Castañeda M., y Reyes A., (2004), Solución numérica de la Ecuación de Richards, Trabajo de grado, Universidad Industrial Santander, Bucaramanga. Gálvez, P., (2014), Modelo Matemático, numérico y computacional para la contaminación de aguas subterráneas, Trabajo de graduación, Universidad Central del Ecuador, Quito-Ecuador. Guanotoa E., (2015), Modelo Matemático, numérico y computacional para la predecir el contenido de agua en suelos agrícolas del territorio ecuatoriano, Trabajo de graduación, Universidad Central del Ecuador, Quito-Ecuador. Pineda R., (2018), Resolución de la ecuación de Richards unidimensional por el método de volúmenes finitos, Universidad Central del Ecuador, Quito – Ecuador. Florian List, A study on iterative methods for solving Richards’ equation, Published online: 15 March 2016. Entre otros documentos. 44 GRACIAS image7.png image8.png image9.png image1.png image2.png image3.png image4.png image5.png image6.png image10.png image4.jpg image11.png image12.png image13.png image14.png image15.png image16.png image310.png image111.png image120.png image130.png image140.png image150.png image160.png image170.png image180.png image19.png image20.png image21.png image22.png image17.png image55.png image18.png image70.png image80.png image90.png image100.png image171.png image23.png image24.png image25.png image26.png image27.png image28.png image35.png image36.png image37.png image38.png image39.png image40.png image29.png image30.png image260.png image31.png image32.png image33.png image34.png image41.png image42.png image43.png image380.png image44.png image45.png image44.jpg image46.png image47.png image49.png image460.png image48.png image50.png image51.png image52.png image53.png image510.png image54.png image56.png image57.png image58.png image59.png image60.png image61.png image62.png image570.png image621.png image590.png image600.png image610.png image620.png image490.png image63.png image64.png image65.png image66.png image67.png image68.png image69.png image71.png image72.png