TAREA Nº 4: ADVECCIÓN, DIFUSIÓN Y REACCIÓN Nombre: Danky Hernández Urbano UNIVERSIDAD DE VALPARAÍSO Av. Brasil 1786, Valparaíso, Chile Teléfono: (56) 32-2995630 www.ingenieriaoceanica.cl Asignatura: Modelación Hidráulica. Profesor: Patricio Winckler. Ayudante: Rodrigo Campos. Fecha: 6 de marzo de 2017.  CONTENIDOS INTRODUCCIÓN 3 1 Objetivos 4 1.1 Objetivo General 4 1.2 Objetivos específicos 4 2 PROBLEMA DE LA DIFUSIÓN 5 2.1 Solución Numérica del problema de Difusión 8 2.2 Análisis de sensibilidad para distintos coeficientes de difusión 8 2.3 Límites de Estabilidad del Algoritmo 10 2.4 Resultados cuando no se satisface la condición de Courant 11 3 PROBLEMA DE REACCIÓN 12 3.1 Solución Númerica del problema de reacción: 15 3.2 Análisis de Sensibilidad del coeficiente de decaimiento: 15 4 Problema difusión-reacción 18 4.1 Solución numérica del problema de Reacción-Difusión 20 CONCLUSIÓN 21 REFERENCIAS 22 INTRODUCCIÓN El estudio y comprensión de los procesos que ocurren cuando un contaminante es vertido dentro de masas oceánicas y hacia donde estos se propagan o disuelven, resulta ser relevante para la ingeniería oceánica y la ingeniería ambiental. Los procesos de Advección, Difusión y Reacción pueden coexistir en un proceso de mezcla, dependiendo de la velocidad media del flujo asociadas a la Advección y de las propiedades difusivas y reactivas del proceso, algunos de estos mecanismos pueden dominar sobre otro, para lo cual se definen números adimensionales los cuales cuantifican la relación entre ellos. La Advección es un proceso de transporte de una sustancia o de una propiedad ya sea el calor, la humedad, la salinidad, el momentum o la energía, por efecto del flujo. Y mientras que la Difusión se centra en la distribución de concentración generada por gradientes de densidades en el medio, la Reacción se centra en fenómenos químicos o biológicos que provocan que la concentración de un contaminante dentro de una mezcla, varíe. En el presente informe se busca entender el comportamiento de estos fenómenos mediante el uso de un modelo de diferencias finitas y así resolver las ecuaciones que gobiernan a estas variaciones en la concentración, además de encontrar sus soluciones numéricas bajo parámetros y condiciones de borde previamente establecidas. Estos estudios se realizarán individualmente y combinados a cada uno de estos procesos ya que en la naturaleza se evidencia que ocurren en conjunto. Objetivos Objetivo General Implementar los métodos de resolución numérica aplicados a ecuaciones diferenciales para representar diferentes fenómenos presentes en los procesos de mezcla; Particularmente para los fenómenos de difusión y reacción, a través de la obtención de esquemas de diferencias finitas y solución numérica, para luego determinar la relación entre el espaciamiento y tiempo para así obtener un modelo estable. Objetivos específicos Plantear la expresión en diferencias finitas para cada fenómeno. Encontrar la expresión que permita obtener la solución numérica. Definir la relación entre el salto de tiempo dt y la resolución espacial dx para garantizar la estabilidad del algoritmo. PROBLEMA DE LA DIFUSIÓN La difusión se asocia al movimiento térmico de las partículas de un fluido líquido o gas a temperaturas por encima del cero absoluto. De forma empírica se sabe que el flujo neto de las moléculas ocurre de una región de mayor concentración a una de menor concentración. Considerando un problema de difusión unidimensional, en un dominio –𝐿<𝑥<𝐿, definido por la ecuación: Ecuación SEQ1: Ec. de Gobierno de la Difusión Fuente: “Winckler, P. (2016)” Donde: es la concentración. D es el coeficiente de difusión. Además, se consideran las condiciones de borde de tipo Dirichlet: Y una condición inicial de tipo Gaussiana: Ecuación SEQ2: Condición inicial de tipo Gaussiana para la Difusión Fuente: “Winckler, P. (2016)” Donde Utilizando un esquema de diferencias finitas explícito, de tipo avanzado en el tiempo y centrado en el espacio (FTCS: forward in time, centered in space), lo que a partir de la ecuación del problema de difusión, se obtienen la primera derivada temporal del tipo avanzada y la derivada espacial mediante diferencias centradas usando series de Taylor: Ecuación SEQ3: Expresión discreta para una primera derivada en el tiempo. Fuente: “Winckler, P. (2016)” Ecuación SEQ4: Expresión discreta de una segunda derivada en el espacio. Fuente: “Winckler, P. (2016)” Luego reemplazando estas ecuaciones en la ecuación de difusión, se obtiene la solución numérica: Ecuación SEQ5: Expresión de la ecuación de difusión mediante diferencias finitas. Fuente: “Winckler, P. (2016)” mediante la siguiente figura, se explica que la ecuación anterior necesita 3 puntos denominados STENCIL asumiendo las condiciones iniciales avanzando de izquierda a derecha y pasando por cada fila: Figura SEQ Figura \* ARABIC 1: Esquema explícito del tipo FTCS con derivada espacial evaluada en el tiempo Fuente: “Winckler, P. (2016)” A continuación se muestra el algoritmo programado para el software MATLAB con el cual se calcuraron las soluciones númericas al problema de difusión: Figura SEQ Figura \* ARABIC 2: Algoritmo utilizado por Matlab para encontrar las soluciones al problema de Difusión Fuente: “Elaboración Propia” El algoritmo que se muestra en la figura anterior se compone de la definición de parámetros importantes definiendose el coeficiente de difusión igual a 0,1 [m2/s] y donde dx y dt definen los límites y los intervalos para para establecer el espacio y tiempo de desarrollo del problema respectivamente. El coeficiente k corresponde a un coeficiente arbitrario el cual garantiza la estabilidad del modelo, donde además se muestra la condición inicial gaussiana la cual funciona con L = xmax para sigma según la ecuación 2. Solución Numérica del problema de Difusión A continuación se muestra el comportamiento de la Difusión a través del tiempo utilizando los parámetros especificados en el algoritmo de la figura 2: Gráfico SEQ Gráfico \* ARABIC 1: Solución numérica Ecuación de Difusión Fuente: “Elaboración Propia” Se observa que a medida que transcurre el tiempo la concentración de una masa M y ubicada en la posición x=0, disminuye desde un tiempo igual a 0 [s] con C= 2,5 [ppm] a un tiempo igual 10 [s] con C= 1,5 [ppm], para un coeficiente de difusión de 0,1 [m2/s]. Cabe destacar que el proceso de difusión es un proceso conservativo y por lo cual la masa total del sistema no cambia (al igual que la Advección). Análisis de sensibilidad para distintos coeficientes de difusión A continuación, se muestra el comportamiento de la masa M en función de la concentración y las variaciones que sufre al cambiar el coeficiente de Difusión en cada caso: Gráfico SEQ Gráfico \* ARABIC 2: Análisis de Sensibilidad del Coeficiente de Difusión Fuente: “Elaboración Propia” Se observa la solución numérica de la ecuación de difusión para distintos valores del coeficiente de difusión, donde se concluye que este coeficiente es directamente proporcional a la difusión y cuanto mayor es el valor de este, la difusión también será mayor representada por los valores de concentración de la masa M a lo largo del eje x. Otro punto a destacar es que las mayores concentraciones se observan en el centro del eje x (x=0) mientras que disminuye hacia los bordes. Los coeficientes de difusión utilizados para este análisis fueron 0,001; 0,01; 0,02; 0,04; 0,05; 0,1 [m2/s], mostrados desde arriba hacia abajo respectivamente y bajo el supuesto de que corresponden a coeficientes de difusión de contaminantes. Límites de Estabilidad del Algoritmo Un punto importante a tratar dentro de este estudio es que el orden de aproximación del esquema numñerico denominado FTCS (Foward Time Centered Space numerical scheme) es bastante bajo, lo que implica que la exactitud de la solución aporximada es también baja. Para garantizar la estabilidad del algoritmo se debe cumplir con la condición de Courant, la cual limita la longitud del paso temporal que es necesario elegir una vez se ha elegido un paso espacial. Ecuación SEQ Ecuación \* ARABIC 6: Condición de Courant Fuente: “Winckler, P. (2016)” Al sobrepasar estas limitantes resulta en la inestabilidad del algoritmo al provocar problemas en los cálculos de las ecuaciones por consiguiente un aumento en el error a medida que se realizan.al despejar la ecuación tanto respecto del tiempo y del espacio se obtiene: Ecuación SEQ Ecuación \* ARABIC 7: Valores del espacio temporal cuando el espacio varía Fuente: “Winckler, P. (2016)” Ecuación SEQ Ecuación \* ARABIC 8: Valores del espacio cuando el espacio temporal varía Fuente: “Winckler, P. (2016)” Finalmente se establece que para asegurar la estabilidad del algoritmo en la ecuación de difusión, t debe ser menor a y debe ser mayor . Resultados cuando no se satisface la condición de Courant Utilizando los valores iniciales que se proponen en la figura 2 con un coeficiente de difusión D= 0,1 [m2/s] y un x= 0,1 [m], se obtiene que t= 0,05 [s]. Por lo tanto para encontrar la inestabilidad del algoritmo basta con superar el límite impuesto por t, asi que ultilizando un mayor espacio temporal se obtiene: Gráfico SEQ Gráfico \* ARABIC 3: Inestabilidad del algoritmo al tomar valores mayores del t Asimismo cuando se estudia la inestabilidad con los valores de x, se obtiene que el límite de x según la ecuación 8 es igual a 0,2236 [m] aproximadamente, por lo tanto si tomamos un valor de x menor se conseguira la inestabilidad como se muestra en el siguiente gráfico: Gráfico SEQ Gráfico \* ARABIC 4: Inestabilidad del Algoritmo al tomar valores menores de x  PROBLEMA DE REACCIÓN Existen reacciones de tipo químico, físico y biológico. Una reacción esta asociado a un proceso en el cual una o más sustancias reactantes se transforman en otras sustancias denominadas productos (en el caso de una reacción química). La ecuación de reacción es: Ecuación SEQ Ecuación \* ARABIC 9: Ecuación de Reacción Fuente: “Winckler, P. (2016)” Para un problema de reacción en una dimensión, en un dominio –𝐿<𝑥<𝐿, esta definido por la ecuación: Ecuación SEQ Ecuación \* ARABIC 10: Ecuación de Reacción en una Dimensión Fuente: “Winckler, P. (2016)” Donde k es el coeficiente de decaimiento de unidades [1/s], tomando un valor positivo y considerando la condición inicial evaluada en x=0, se procede a resolver mediante separación de variables: Ecuación SEQ Ecuación \* ARABIC 11: Método de separación de Variables para el problema de Reacción Fuente: “Winckler, P. (2016)” Obteniendo el resultado: Considerando que corresponde a la concentración en el instante inicial t0 y aplicando exponencial se obtiene que: Ecuación SEQ Ecuación \* ARABIC 12: Solución Numérica del problema de Reacción Fuente: “Winckler, P. (2016)” que indica que para k>0, la concentrtación disminuye exponencialmente en el tiempo y para k<0 existe un aumento exponecial de la sustancia, lo cual asume un proceso no conservativo. Si k=0 no existe reacción y la concentración permanece constante y conservativo. Gráfico SEQ Gráfico \* ARABIC 5: Comportamiento Teórico de la Concentración para diferentes valores de k en un problema de Reacción Fuente: “Winckler, P. (2016)” Igualmente este problema expresado en diferencias finitas queda como: Ecuación SEQ Ecuación \* ARABIC 13: Expresión del problema de Reacción en diferencias finitas Fuente: “Winckler, P. (2016)” Otra cosa importante a destacar es que el problema de reacción no necesita condiones de borde ya que depende solo del tiempo. A continuación se muestra el algoritmo utilizado en el software MATLAB con el cual se calcularon las soluciones numéricas al problema de Reacción. También se destacan los parámetros ultilizados tales como la constante de decaimiento k=5, el tiempo de cómputo de 1,5 [s] y el incremento temporal dt= 0,001 [s]. Figura SEQ Figura \* ARABIC 3:Algoritmo utilizado por el Programa Matlab para encontrar la solución al Problema de Reacción Fuente: “Elaboración Propia” Solución Númerica del problema de reacción: A continuación se muestra el comportamiento de la reacción a través de una ventana de tiempo de 1,5 [s] y un valor inicial de la constante de decaimiento de k=5: Gráfico SEQ Gráfico \* ARABIC 6: Solución Numérica para el problema de Reacción Fuente: “Elaboración Propia” El gráfico anterior muestra el decaimiento de la concentración de la sustancia acercandose asintoticamente a cero, el cual se asemeja al comportamiento de una función exponencial, solución de la ecuacción de reacción, validando el resultado del metodo de de diferencias finitas de manera teórica. Además se obtiene un error asociado de 0,19871 para un incremento temporal de 0,001 [s]. Análisis de Sensibilidad del coeficiente de decaimiento: Para el análisis de sensibilidad del coeficiente de decaimiento se hizo variar el parámetro k dentro del algoritmo de la figura 3, y del cual se espera obtener resultados del tipo estipulado en el gráfico 5. Para esto se evaluaron 8 casos en donde las constantes de decaimiento k, tomaron los valores de -10; 0; 5; 10; 20; 30; 40; 50 respectivamente para cada caso mostrado desde arriba hacia abajo. Gráfico SEQ Gráfico \* ARABIC 7: Análisis de Sensibilidad para el problema de Reacción Fuente: “Elaboración Propia” Se observa en el gráfico 7 como la variabilidad de la constante de decaimiento produce cambios sustanciales desde que toma un valor negativo igual a -10 para el primer caso, como para cuando toma valor igual a cero en el segundo caso y para los valores siguientes los cuales toman valores positivos en ambos casos. Lo que deja las siguientes conclusiones tales como cuando la constante de decaimiento adopta un valor positivo se dice que el proceso es conservativo ya que el sistema no pierde ni gana masa, solo se transforma. Además, a medida que aumenta la constante de decaimiento k menor es el tiempo con el cual se produce la reacción ver el gráfico 8 y la concentración se acerca asintóticamente a cero. Cuando la constante de decaimiento adopta el valor cero se dice que no hay reacción y se conserva la masa y, por último, cuando la constante de decaimiento adopta valores negativos se dice que el proceso no es conservativo ya que proceso indica que aumenta la concentración a través del tiempo y eso es prácticamente imposible a menos que se agregue masa al sistema, por lo tanto, es un proceso físicamente irreal para nuestro caso. Otro punto a destacar es el error generado por el algoritmo, ya que mientras que para un k negativo adopta un valor del orden de 1011, con k positivos y a medida que aumenta k, también aumenta el error con valores que van desde 0,19871 a 0,2013 en un mismo intervalo de tiempo dt= 0,001 [s]. Gráfico SEQ Gráfico \* ARABIC 8: Análisis de Sensibilidad para el problema de Reacción Fuente: “Elaboración Propia” Problema difusión-reacción Para el problema de reacción en una dimensión, en un dominio –𝐿<𝑥<𝐿, definido por la ecuación: Ecuación SEQ Ecuación \* ARABIC 14: Problema de Reacción-Difusión Fuente: “Winckler, P. (2016)” Donde k es el coeficiente de decaimiento y C el coeficiente de difusión, además, considerando la condición inicial de tipo gaussiana ecuación 2, se reemplaza y se obtiene: Ecuación SEQ Ecuación \* ARABIC 15: Expresión de la Ecuación 14 reemplazada con la condición inicial, la primera derivada en el tiempo y la segunda derivada en el espacio Fuente: “Winckler, P. (2016)” Finalmente por simplicidad se considera el factor = y = , de esta forma se obtiene: Ecuación SEQ Ecuación \* ARABIC 16: Expresión en diferencias finitas del problema de Reacción-Difusión Fuente: “Winckler, P. (2016)” A continuación se muestra el algoritmo que fue utilizado por el software MATLAB para calcular las soluciones numéricas del problema de Reacción-Difusión:  Figura SEQ Figura \* ARABIC 4:Algoritmo utilizado por el Programa Matlab para encontrar la solución al Problema de Reacción-Difusión Fuente: “Elaboración Propia” Solución numérica del problema de Reacción-Difusión Para la obtención de los resultados se ocuparon los siguientes parámetros, un coeficiente de Difusión D= 0,04 [m2/s], una constante de decaimiento de k=5 y un coeficiente de estabilidad de n=0,1 en un tiempo máximo de 4 [s]. Gráfico SEQ Gráfico \* ARABIC 9: Solución numérica al problema de Reacción-Difusión Fuente: “Elaboración Propia” A través de los puntos anteriores sobre las ecuaciones de Reacción y Difusión se puede obtener que los efectos de reacción suceden en menor tiempo ya que el cambio que existe entre la concentración inicial de 2,5 [ppm] a 1 [ppm] ocurre en un intervalo muy pequeño de tiempo. Otro punto importante a destacar es que ocurre mucho mas Difusión que en el caso donde no hay interacción con la reacción, debido a que el rango disminuye hasta 0,1 [ppm] mientras que en el caso de la difusión sola se alcanza una concentración mínima de 1,5 [ppm]. CONCLUSIÓN Tras haber concluido con el estudio de Difusión y Reacción se puede concluir lo siguiente: Hay una gran importancia con los rangos en donde se estudia el fenómeno adeás de la definición de los parámetros que definen el problema ya sea la constante de decaimiento o los coeficientes de difusión. La variación de los parámetros que definen el problema podrían alterar la importancia de cada fenómeno por sobre los otros, cuando estos funcionan en conjunto. Mientras mayor sean los coeficientes de Difusión mayor será la efectividad y velocidad de la Difusión del contaminante en la masa de agua. En cambio para un coeficiente de decaimiento menor implica un mayor tiempo en el que demoran en reaccionar y disminuir su concentración.  REFERENCIAS Winckler, P. (2010). CEE 6550. Transport, Mixing and Transformation in the Environment, Assignment 2: Molecular Diffusivity. Cornell University. Winckler, P. (2016). Modelado de procesos costeros. Apuntes del curso de Modelación Hidráulica. Ingeniería Civil Oceánica. Universidad de Valparaíso. Advección, Difusión y Reacción Facultad de Ingeniería Escuela de Ingeniería Civil O c e á n i c a 15