Difusión de una sustancia contaminante en un tubo (grupo B17)

Trabajo realizado por estudiantes
Título	Difusión de una sustancia contaminante en un tubo (grupo B17)
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2013-14
Autores	Eduardo Pereiras Navas Álvaro Ramírez Fernández de la Puente Jesús Infestas Robles Miguel Fandiño Álvarez Rachel L’Hostis
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

El estudio de la difusión de una sustancia en un medio determinado, se considera un problema bastante común dentro de diversas áreas de la ingeniería civil, como puede ser el tratamiento de aguas.

Conocer cómo se distribuye la concentración en un medio cualquiera implica la introducción de muchas variables que pueden afectar a esta. Para obtener soluciones más sencillas, consideraremos el caso particular de un tubo hueco largo (por ejemplo, una tubería), considerando un tramo de este, en el cual se aislarán las superficies laterales.

1 Planteamiento y modelización del problema

Para una modelización más sencilla, tendremos en cuenta varias consideraciones para simplificar el cálculo.

Suponemos que solamente se trata de una sustancia en suspensión, definimos su concentración con la variable [math] u (\frac{mol}{m^2})[/math], por unidad de superficie, constante en cada sección determinada del tubo. De esta forma, [math]u[/math] será una variable que exclusivamente dependerá de la sección tomada en el tubo , definida por [math]x \in(0,L)[/math], con [math]L [/math] como la longitud del tubo; y el tiempo [math]t \geq 0[/math], en que se encuentre. En este caso, consideraremos la longitud del tubo [math] L=5 [/math].

Además, consideraremos la simplificación de que se trata de un tubo aislado por sus laterales, pero abierto por las bases del cilindro.

1.1 Base física

Para poder conocer la concentración en cada instante y punto de la tubería, haremos uso de dos principios físicos: la ley de conservación de la masa, y la Ley de Fick.

1.1.1 Ley de conservación de masa

La ley de conservación de masas exige que no haya una generación espontánea de materia. Es decir, si en la región del espacio estudiada, no existe entrada o salida por un punto cualquiera (frontera o interior), la masa en esta región se mantiene constante. Esta suposición, en términos de variación de masa, puede quedar determinada como:

Variación de la de masa en el tiempo = flujo de masa a través de la frontera por unidad de tiempo + flujo de masa en los puntos interiores por unidad de tiempo.

1.1.2 Ley de Fick

Se trata de una ley cuantitativa en forma diferencial que describe la difusión de materia o energía en un medio en el que inicialmente existe equilibrio químico o físico. Esta ley se trata de una generalización de los fenómenos descritos por la Ley de Fourier, ya que esta es una consideración concreta del calor.

La ley de Fick anuncia que en el caso de existir diferencias de concentración en cualquier especie, el paso de estas especies se llevará a cabo desde las regiones con mayor concentración a las de menor. Formulado matemáticamente:

[math]F(x,t)=-D\frac{\partial u}{\partial x}[/math]

Donde [math]F[/math] es el flujo de la sustancia y [math]D[/math] el coeficiente de difusión, que representa la facilidad en que un soluto particular se mueve por un disolvente determinado.

1.2 Modelización de la ecuación diferencial

En base a las dos leyes anteriores, se puede definir a partir de la variable [math] u(x,t)[/math] en unidades de sustancia por área, es decir, concentración por superficie del tubo.

Consideramos [math] A [/math] como variable superficie, y tomamos un pequeño trozo, definido por [math]\Delta x[/math]. Entonces, para obtener la variación de la masa en el tiempo:

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

Y derivando en el tiempo:

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

obtenemos el primer término de la igualdad anteriormente enunciada.

En el segundo, consideramos nulo el flujo de masa en los puntos interiores por unidad de tiempo, pues suponemos la tubería perfecta y sin sumideros.

Como la superficie lateral del tubo está aislada, solamente habrá flujo de sustancia en la dirección del eje [math]x[/math]. Si [math]F(x,t)\gt 0[/math] pensamos que el flujo va hacia la derecha, si [math]F(x,t)\lt 0[/math], este va a hacia la izquierda. 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 resultados anteriores:

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

y seguidamente pasamos dividiendo [math]A \Delta x[/math]:

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

Si hacemos [math]\Delta x \rightarrow 0 [/math]:

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

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]

Y de esta forma podemos obtener la ecuación en derivadas parciales que nos define la concentración de sustancia:

[math]u_t(x,t)- D u_{xx}(x,t)= 0 [/math]

1.3 Planteamiento del problema bien propuesto

Esta ecuación en derivadas parciales exige al menos de tres condiciones para su resolución. Dos de ellas serán de contorno y la restante, condición inicial.

Como ejemplo, suponemos que existe una situación inicial, definida por [math]u_0(x,0)[/math], tal que:

[math] u_0(x,0)=\begin{cases} 0, x\leq 3\\ 3,x\gt3 \end{cases} [/math]

Situación inicial de las concentraciones

En la situación de contorno que utilizaremos, supondremos que existe una tapa en los extremos del tubo, por lo que el flujo de sustancia a través de ello es nulo (y en consecuencia, la variación de este, su derivada, es cero).

Por lo que al final podemos escribir el problema mediante:

[math]\begin{cases} u_t(x,t)- D u_{xx}(x,t)= 0 \\ u_0(x,0), x\in (0,L) \\ u_x(0,t)=0, u_x(L,t)=0, t\geq 0 \end{cases} [/math]

2 Resolución numérica

2.1 Caso particular 1: fronteras aisladas y concentraciones iniciales conocidas

En el primer caso particular que estudiaremos, trataremos de obtener una resolución numérica del problema anteriormente planteado:

[math]\begin{cases} u_t(x,t)- D u_{xx}(x,t)= 0 \\ u_0(x,0), x\in (0,L) \\ u_x(0,t)=0, u_x(L,t)=0, t\geq 0 \end{cases} [/math]

2.1.1 Resolución por el método de las diferencias finitas

La resolución por el método de las diferencias finitas se basa en expresar el problema continuo a un sistema discreto de ecuaciones diferenciales de primer orden.

Particularizando para nuestro caso con las condiciones dadas, el problema planteado sería:

[math]\begin{cases} u'_n(t)+\frac{-u_{n-1}(t)+2u_n(t)+u_{n+1}(t)}{h^2} \\ \frac{u_1(t)-u_-1(t)}{2h} \rightarrow u_1(t)=u_-1(t) \\ \frac{u_{N+1}(t)-u_{N-1}(t)}{2h} \rightarrow u_{N+1}(t)=u_{N-1}\\ u_n(0)=u^0(x) \end{cases} [/math]

Por lo que tenemos un sistema de la forma:

[math]\begin{cases} U'(t)+KU(t)=F(t)=0\\ U(0)=U^0 \end{cases} [/math]

Con [math]N+1[/math] ecuaciones. Los datos de la discretización que utilizaremos serán: [math] \Delta x=0.1[/math] y [math]\Delta t= \Delta x/4[/math].

2.1.1.1 Resolución por método del trapecio

clear all
%%Problema de difusión de barrera
%Datos
L=5;
T=6;
%Discretización temporal y espacial
dx=0.1;
N=L/dx;
dt=dx^2/2;
x=0:dx:L;
t=0:dt:T;
%Matriz de coeficientes
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(N+1,N)=-2; 
K(1,2)=-2;
K=K/dx^2;
%Condición inicial
U0=(3.*(x>3))';
sol(1,:)=[U0'];
U=U0;
%Bucle iterativo aplicando método de diferencias finitas con método del trapecio
for j=length(t)-1
U=(eye(N+1)+(dt*K)/2)\((eye(N+1)-(K*dt)/2)*U);
sol(j+1,:)=U';
end
[xx,tt]=meshgrid(x,t);
figure(1)
mesh(xx,tt,sol);
xlabel('x');
ylabel('t');
grid

Archivo:Grafica1b

2.1.1.2 Resolución por Euler explícito

clear all
%%Problema de difusión de barrera
%Datos
L=5;
T=6;
%Discretización temporal y espacial
dx=0.1;
N=L/dx;
dt=dx^2/2;
x=0:dx:L;
t=0:dt:T;
%Matriz de coeficientes
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(N+1,N)=-2; 
K(1,2)=-2;
K=K/dx^2;
%Condición inicial
U0=(3.*(x>3))';
sol(1,:)=[U0'];
U=U0;
%Bucle iterativo aplicando método de diferencias finitas con método del trapecio
for j=length(t)-1
U=U+dt*K*U;
sol(j+1,:)=U';
end
[xx,tt]=meshgrid(x,t);
figure(1)
mesh(xx,tt,sol);
xlabel('x');
ylabel('t');
grid

2.1.1.3 Resolución por Euler implícito

clear all
%%Problema de difusión de barrera
%Datos
L=5;
T=6;
%Discretización temporal y espacial
dx=0.1;
N=L/dx;
dt=dx^2/2;
x=0:dx:L;
t=0:dt:T;
%Matriz de coeficientes
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(N+1,N)=-2; 
K(1,2)=-2;
K=K/dx^2;
%Condición inicial
U0=(3.*(x>3))';
sol(1,:)=U0';
U=U0;
%Bucle iterativo aplicando método de diferencias finitas con método del trapecio
for j=length(t)-1
U=(eye(N+1)+dt*K)\(U);
sol(j+1,:)=U';
end
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)

2.1.1.4 Método de Euler Modificado

clear all
%%Problema de difusión de barrera
%Datos
L=5;
T=6;
%Discretización temporal y espacial
dx=0.1;
N=L/dx;
dt=dx^2/2;
x=0:dx:L;
t=0:dt:T;
%Matriz de coeficientes
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(N+1,N)=-2; 
K(1,2)=-2;
K=K/dx^2;
%Condición inicial
U0=(3.*(x>3))';
sol(1,:)=[U0'];
U=U0;
%Bucle iterativo aplicando método de diferencias finitas con método del trapecio
for j=length(t)-1
k1=-K*U;
k2=-K*(U+dt*k1);
U=U+(dt/2)*(k1+k2);
sol(j+1,:)=U';
end
[xx,tt]=meshgrid(x,t);
figure(1)
mesh(xx,tt,sol);
xlabel('x');
ylabel('t');
grid

2.1.2 Resolución por el método de Fourier

[math](P)=\begin{cases} u_t(x,t)- D u_xx(x,t)= 0 \\ u_0(x,0), x\in (0,L) \\ u(0,t)=0, u_x(L,t)=0, t\geq 0 \end{cases}[/math]

El método de Fourier se basa en buscar soluciones del tipo:

[math] u(x,t)= \phi (x) T(t)[/math]

Si consideramos [math] \phi(x) [/math] solución del problema de autovalores [math](PA)[/math] asociado a [math](P)[/math]:

[math](PA)=\begin{cases} \phi''(x)- \lambda \phi(x)= 0 \\ \phi (0)=0, \phi '(L)=0 \end{cases}[/math]

Cuyas autofunciones y autovalores, respectivamente [math]\phi_k(x)=cos(\frac{\pi x}{L})[/math] y [math] \lambda_k=\frac{\pi^2k^2}{L^2}[/math]

Entonces la solución de [math](P)[/math] será: [math] u(x,t)=\phi(x)e^{-\lambda t}[/math]

Por lo que la solución del problema del tubo, para las condiciones de contorno [math]u(0,t)=0, u_x(L,t)=0[/math], será:

[math] u(x,t)=cos(\frac{\pi x}{L})e^{-\frac{\pi^2k^2}{L^2} t}[/math]

2.1.3 Comprobación de la conservación de la masa

La propia ley de conservación de materia, en la que se expone la no creación espontánea (en el caso en que las condiciones de fronteras exigen al tubo cerrado), se puede comprobar analíticamente como la suma de todas las variaciones de sustancia en la varilla y debe ser nula. Expresado analíticamente:

[math]\frac{d}{dt}\displaystyle\int_{0}^{L}u(x,t)\,dx=\displaystyle\int_{0}^{L}u_{xx}(x,t)\,dx= u_x(x,t)]^L_0=0[/math]

2.2 Caso particular 2: Limpiador en un extremo

Consideramos ahora en el problema que, con las mismas condiciones que antes, se modifica la condición de contorno en el extremo [math] x=0 [/math] donde se coloca un limpiador. Esto matemáticamente se traduce como que en ese punto la concentración es en todo instante nula. Podemos reescribir el problema como

[math]\begin{cases} u_t(x,t)- D u_xx(x,t)= 0 \\ u_0(x,0), x\in (0,L) \\ u(0,t)=0, u_x(L,t)=0, t\geq 0 \end{cases} [/math]

2.2.1 Resolución por el método de las diferencias finitas

Esta vez, podemos reescribir el problema de la forma:

[math]\begin{cases} u'_n(t)+\frac{-u_{n-1}(t)+2u_n(t)+u_{n+1}(t)}{h^2} \\ u_0(t)=0\\ \frac{u_{N+1}(t)-u_{N-1}(t)}{2h} \rightarrow u_{N+1}(t)=u_{N-1}\\ u_n(0)=u^0(x) \end{cases} [/math]

Por lo que tenemos un sistema de la forma:

[math]\begin{cases} U'(t)+KU(t)=F(t)=0\\ U(0)=U^0 \end{cases} [/math]

Con [math]N[/math] ecuaciones.

2.2.1.1 Resolución numérica

clear all
%%Problema de difusión de barrera
%Datos
L=5;
T=6;
%Discretización temporal y espacial
dx=0.1;
N=L/dx;
dt=dx^2/2;
x=0:dx:L;
xint=dx:dx:L;
t=0:dt:T;
%Matriz de coeficientes
K=2*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
K(N,N-1)=-2; ;
K=K/dx^2;
%Condición inicial
U0=(3.*(xint>3))';
sol(1,:)=[0,U0'];
U=U0;
%Bucle iterativo aplicando método de diferencias finitas con método del trapecio
for j=length(t)-1
U=(eye(N)+dt*K)\(U);
sol(j+1,:)=U';
end
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)

2.2.2 Modificación de las condiciones de fronteras

Supondremos una ligera modificación en el problema respecto a las condiciones de contorno, en las que en el extremo izquierdo consideraremos el flujo de salida de la sustancia de forma constante e igual a tres ([math] u_x(0,t)=3[/math]).

En el extremo derecho, consideraremos el caso en que se conecta un aparato que extraiga y a la vez introduzca sustancia de forma armónica, en un intervalo que tiene como máximo de entrada [math]u(5,t)=10[/math] y como máximo de salida [math]u(5,t)=-10[/math], por lo que seguirá una función de la forma [math]u(5,t)=10sin(t)[/math].

2.2.2.1 Resolución numérica

clear all
%%Problema de difusión de barrera
%Datos
L=5;
T=10;
%Discretización temporal y espacial
dx=0.1;
N=L/dx;
dt=dx^2/2;
x=0:dx:L;
xint=dx:dx:L;
t=0:dt:T;
%Matriz de coeficientes
K=2*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
K(N,N-1)=-2; ;
K=(1/dx^2)*K;
%Condición inicial
U0=(3.*(xint>3))';
sol(1,:)=[U0',10*sin(0)];
U=U0;
%Bucle iterativo aplicando método de diferencias finitas con método del trapecio
for j=length(t)-1
U=(eye(N)+(dt/2)*K)\(U-(dt/2)*K*U);
sol(j+1,:)=[U', 10*sin(t(j+1))];
end
[xx,tt]=meshgrid(x,t);
figure(1)
mesh(xx,tt,sol);
xlabel('x');
ylabel('t');
grid