Difusión de una sustancia contaminante (Grupo 24C)

Trabajo realizado por estudiantes
Título	Difusión de una sustancia contaminante. Grupo 24-C
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2014-15
Autores	Jose Antonio Martinez Montalvo 1494 Jorge Sempere Ruíz 4 Isaac Rebollo Palos 1522 Daniel Pascual Cobos 1690 Rodrigo Bellot Rodriguez 1270
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

El trabajo tiene como objetivo el estudio de la mezcla de dos sustancias, una de ellas contaminante, en un tubo largo.

Su planteamiento se realizará de forma similar al de la transmisión de calor a lo largo de una varilla, visto en clase.

Analizaremos el modelo que define nuestro problema (ecuación de difusión), varios puntos teóricos del mismo (la conservación de la masa contaminante, la solución estacionaria) y terminar por plantear un problema concreto el cual resolveremos mediante métodos numéricos, modificando según que condiciones para poder interpretar las distintas situaciones que pueden modelizar nuestro problema.

2 Ecuación de difusión

2.1 Deducción

Sea tubo largo de longitud L, orientado en la dirección x y cuya sección es constante desde x=0 hasta x=L. Supondremos el tubo aislado tanto superficialmente como lateralmente.

Configuración inicial de las concentraciones de contaminante en el tubo.

Sean dos sustancias contenidas en el tubo, una de ellas contaminante. Denotaremos por u(x,t) la concentración de contaminante, la cual dependerá únicamente de la posición del tubo y del tiempo, manteniéndose constante en cada sección transversal del tubo. Medido en mol/m2s.

Sea F(x,t) el flujo de contaminante (análogo al flujo de calor), definiéndose como la cantidad de contaminante que atraviesa una sección transversal por unidad de tiempo y area, dependiendo del número de moles, siendo éste último proporcional a la masa de la sustancia. Debido a los aislantes laterales y superficiales, el flujo sera a lo largo del eje x.

Debido al principio de conservación de la masa la variación de la concentración de contaminante en cada posición del tubo en función del tiempo, es igual a la suma del flujo de contaminante a través de los extremos del tubo por unidad de tiempo, más la concentración de contaminante generada en el interior por unidad de tiempo. Suponemos igual a 0 las perdidas y/o ganancias en el interior del tubo (problema homogéneo). La ley de Flick (similar a Furier) determina que el flujo de difusión del contaminante es proporcional a la variación de concentración: [math]F(x,t)=-D\frac{\partial u}{\partial x}[/math]

siendo D el coeficiente de difusión (medido en m2/s), que dependerá de las propiedades químicas de los compuestos.

De esta forma y continuando con la semejanza a la transmisión de calor, la concentración de sustancia contaminante en un instante de tiempo t y en una sección transversal que dista x unidades del extremo izquierdo del tubo satisfará: [math]u_t(x,t)=-\frac{\partial }{\partial x}(-D u_x(x,t)) = D u_{xx}(x,t)[/math]

Pasando todo al termino de la izquierda obtendremos la ecuación de difusión: [math]u_t(x,t)- D u_{xx}(x,t)= 0[/math]

dónde D es el coeficiente de difusión anteriormente citado.

Para ello hemos definido [math] A [/math] como la variable de superficie, tomando una sección pequeña del tubo, designada por [math]\Delta x[/math]: La variación de la concentración de contaminante en función del tiempo es:

[math] A u(x,t) \Delta x [/math]

Derivamos respecto del tiempo:

[math] A u_t(x,t) \Delta x [/math]

A continuación, se considera nula la concentración de contaminante generada en el interior por unidad de tiempo debido a la ausencia de sumideros.

Suponemos que Δx>0 y la concentración de contaminante en un tiempo t es menor en x+ Δx que en x, entonces u(x+ Δx) – u(x,t) < 0, y como Δx es pequeño, se tiene que ux(x,t)<0 y el flujo de difusión del contaminante es positivo y va hacia la derecha en la dirección del eje x.

El flujo de calor en un intervalo [math][x, x + ∆x ] [/math] viene dado dado por:

[math]F(x,t)A-F(x+\Delta x, t)A[/math]

Igualando los términos anteriores y dividiendo por [math]A \Delta x[/math] se obtiene:

[math]A u_t(x,t) \Delta x = F(x,t)A-F(x+\Delta x, t)A[/math]:

[math]u_t(x,t)= \frac {F(x,t)-F(x+\Delta x, t)}{\Delta x}[/math]

Haciendo que [math]\Delta x[/math] tienda a 0:

[math]u_t(x,t)= -F_x(x,t)[/math]

Y aplicando la ley de Fick:

[math]u_t(x,t)=-\frac{\partial }{\partial x}(-D u_x(x,t)) = D u_xx(x,t) [/math]

Así se obtiene la ecuación de difusión de la sustancia contaminante a lo largo del tubo anteriormente citada:

[math]u_t(x,t)- D u_xx(x,t)= 0 [/math]

3 Problema Propuesto

3.1 Sistema de ecuaciones

Sea el problema:

[math] (P)\left\{\begin{matrix}\\u_t-u_{xx}=0 ,x\in\ (0,7), t\gt0\\u_x(0,t)=0, t\gt0\\u_x(5,t)=0, t\gt0\\u(x,0)=u_0, x\in\ (0,7)\end{matrix}\right. [/math]

Sean las condiciones iniciales:  :[math] u(x,0)=\left\{\begin{matrix}\\0, x≤5\\3, x\gt5\end{matrix}\right. [/math]

Pasamos a su interpretación física.

3.2 Modelización

El problema modeliza el proceso de difusión de una sustancia contaminante a lo largo de un tubo de longitud x=7, delgado, difusor, con superficie lateral aislada y con coeficiente de difusión D igual a 1.

En los puntos interiores del tubo no existentes ni fuentes ni sumideros de sustancia contaminante.

Tanto el extremo izquierdo como el extremo derecho del tubo se encuentran aislados de forma que el flujo de sustancia en ellos es igual a 0.

En el instante inicial, las secciones transversales que distan 5 o menos unidades de longitud del extremo izquierdo tendrán una concentración inicial de sustancia contaminante igual a 0. Las secciones que distan mas de 5 unidades de ese extremo tendrán una concentración de contaminante igual a 3.

3.3 Resolución

3.3.1 Método diferencias finitas

Se trata de un método numérico que nos proporciona una solución aproximada de distintos problemas de ecuaciones en derivadas parciales, aplicable a nuestro caso. Primero realizamos una discretización del espacio, tomando como longitud de paso h y desde x=0 hasta x=7, obteniendo un vector N con x+1 elementos.

Seguidamente aplicamos la ecuacion diferencial a cada nodo interior xn, aproximamos u_xx mediante: [math]u_{xx}(x,t)\simeq\frac{u(x_{n-1},t)-2u(x_n,t)+u(x_{n+1},t)}{h^2}=0[/math]

Aplicamos la notación [math] u_n(t) = u(x_n,t)[/math] resulta el sistema de N-1 ecuaciones: \begin{array}{c}u'_n(t)+\frac{-u_{n-1}(t)+2u_n(t)-u_{n+1}(t)}{h^2}=0\end{array} para n=1,2...N-1

Aplicamos las condiciones de contorno:

\begin{array}{c}u'_0(t)=0\\u'_n(t)=0\end{array} Resultando el sistema: \begin{array}{c}U'+KU=0=F\\U(0)=U^0\end{array}

Siendo [math]U^0[/math] las condiciones iniciales en el instante inicial se verifica  :[math] u(x,0)=\left\{\begin{matrix}\\0, x≤5\\3, x\gt5\end{matrix}\right. [/math]

3.3.1.1 Trapecio

Se resuelve el sistema en primer lugar utilizando el método del trapecio tomando [math]\triangle t= \triangle x/4[/math] en [math]t\in\ [0,7][/math], con [math]\triangle x=0.1[/math]

El código en MATLAB del método sería:

% Aproximar la ecuacion de difusion de una sustancia contaminante
% u_t-qu_xx=0, x en (0,L)
% u_x(0,t)=0 Condicion Neumann
% u_x(L,t)=0 Condicion Neumann
% u(x,0)=u0(x)
%%%
% Datos del problema
L=7;
T=7;
q=1;
% Datos de la discretizacion espacial
N=50;
h=L/N;
% Vector de nodos en en espacio
x=0:h:L; 
% Aproximacion de -q*u_xx
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
 K(1,2)=-2;
 K(N+1,N)=-2;
 K=q*K/h^2;
 % Calculamos u0
 u0=[zeros(1,1+N*3/5),3*ones(1,N*2/5)]';
 % Discretizacion temporal
dt=h/4; % Paso
t=0:dt:T; % Vector de tiempos
% Definimos F
F=zeros(N+1,1); 
% Guardamos la solucion
sol(1,:)=u0'; 
U=u0; % Inicializacion
% Calculamos U_n ----> U_n+1
for j=1:length(t)-1 % Voy hasta la longitud de t menos 1 porque ya conozco un valor
U=(eye(N+1)+(dt*K)/2)\((eye(N+1)-(K*dt)/2)*U); % Trapecio
sol(j+1,:)= U'; % Guardamos solucion
end
 % Dibujamos la solucion
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol);

Se obtiene la siguiente gráfica concentración / espacio-tiempo

3.3.1.2 Euler explícito

% Aproximar la ecuacion de la difusión de una sustancia contaminante
% u_t-qu_xx=0, x en (0,L)
% u_x(0,t)=0 Condicion Neumann
% u_x(L,t)=0 Condicion Neumann
% u(x,0)=u0(x)
%%%
% Datos del problema
L=7;
T=7;
q=1;
% Datos de la discretizacion espacial
N=50;
h=L/N;
% Vector de nodos en en espacio
x=0:h:L; 
% Aproximacion de -q*u_xx
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
 K(1,2)=-2;
 K(N+1,N)=-2; 
 K=q*K/h^2;
 % Calculamos u0
 u0=[zeros(1,1+N*3/5),3*ones(1,N*2/5)]';
 % Discretizacion temporal
dt=h/20; % Paso para euler explicito
t=0:dt:T; % Vector de tiempos
% Definimos F
F=zeros(N+1,1);
% Guardamos la solucion
sol(1,:)=u0'; 
U=u0; % Inicializacion
% Calculamos U_n ----> U_n+1
for j=1:length(t)-1 % Voy hasta la longitud de t menos 1 porque ya conozco un valor
U=U-dt*K*U; % metodo explicito
sol(j+1,:)= U'; % Guardamos solucion
end
 % Dibujamos la solucion
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol);

La gráfica concentración / espacio-tiempo obtenida es: