Revisión del 22:45 19 may 2014

Tubo con sustancia contaminante

La difusión implica un movimiento espontáneo y desordenado de moléculas individuales y, especialmente, nos interesa la difusión de moléculas de una sustancia disuelta en otra.

La difusión de una sustancia puede considerarse análoga al flujo de calor, y la ley de Fick establece que el ritmo de difusión por unidad de superficie, en dirección perpendicular a ésta, es proporcional al gradiente de la concentración de esa sustancia en esa dirección. La concentración es la masa de la sustancia por unidad de volumen, y el gradiente de concentración es la variación de concentración por unidad de distancia.

Para determinar como se distribuye la concentración de una sustancia contaminante tomaremos como ejemplo un tubo largo compuesto por dos sustancias de las cuáles una de ellas es contaminante.

1 Simplificación del problema

Para que este problema resulte más sencillo aplicaremos una serie de condiciones, facilitando el cálculo del mismo.

En primer lugar denotaremos su concentración [math] u [/math], concentración de contaminante en cada posición del tubo y constante en cada sección, medida en [math] \frac{mol}{m^2*s} [/math].

Supondremos que el tubo mide L=5 m, ocupando el intervalo [math] x \in (0,L) [/math], y que la concentración es la misma en cualquier punto de la sección transversal del tubo. Por tanto, la concentración dependerá unicamente de dos variables [math] u = u(x,t) [/math]. En los extremos se ha colocado un aislante que hace que el flujo de contaminante al exterior del tubo sea nulo. 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 Ley de conservación de masa

Tras una serie de experimentos realizados por Lavoisier, se pusieron de manifiesto que si tenemos en cuenta todas las sustancias que forman parte en una reacción química y todos los productos formados, nunca varía la masa. Esta ley exige que no haya una generación espontánea de materia. Es decir, 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 es la ley de la conservación de la masa, que podemos enunciarla, pues, de la siguiente manera:
"En toda reacción química la masa se conserva, esto es, la masa total de los reactivos es igual a la masa total de los productos"
La ley de conservación de masas exige que no haya una generación espontánea de materia. Es decir, 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 masa en el tiempo = flujo de masa a través de la frontera por unidad de tiempo + flujo de masa generada en el interior por unidad de tiempo.

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.

La ley de Fick establece 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]

donde D es el coeficiente de difusión (medido en [math] \frac{m^2}{s} [/math]), que depende de las propiedades físicas de los compuestos, y que se supone constante D=1.

2 Modelización de la ecuación diferencial

Aplicando el principio de conservación de la masa y la ley de Fourier, definimos 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]

Con lo que queda [math]A u_t(x,t) \Delta x = F(x,t)A-F(x+\Delta x, t)A[/math], y al dividir [math]A \Delta x[/math], nos queda [math]u_t(x,t)= \frac {F(x,t)-F(x+\Delta x, t)}{\Delta x}[/math].
Como [math]\Delta x \rightarrow 0 [/math], [math]u_t(x,t)= -F_x(x,t)[/math], y aplicando ley de Fick, [math]u_t(x,t)=-\frac{\partial }{\partial x}(-D u_x(x,t)) = D u_xx(x,t) [/math].
Con lo que resulta una ecuación en derivadas parciales que determina la concentración de sustancia: [math]u_t(x,t)- D u_xx(x,t)= 0 [/math]

2.1 Planteamiento del problema bien propuesto

Suponemos que el flujo será nula al haber una tapa en los extremos del tubo, con lo que el problema bien propuesto quedaría::

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

3 Resolución numérica

Para la resolución de ésta ecuación necesitaremos mínimo tres condiciones, dos de las cuales serán de contorno, y la otra será la condición inicial.Si suponemos [math]u_0(x,0)[/math], entonces:

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

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

En el primer caso particular que estudiaremos, trataremos de obtener una resolución numérica del problema anteriormente planteado, dónde tenemos las fronteras aisladas y las concentraciones iniciales conocidas:

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

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.

3.1.1 Resolución por método del trapecio

Suponiendo que en el instante inicial se verifica  :[math] u(x,0)=\left\{\begin{matrix}\\0, x≤3\\3, x\gt3\end{matrix}\right. [/math]

Si tomamos [math]\triangle t= \triangle x/4[/math] en [math]t\in\ [0,5][/math], con [math]\triangle x=0.1[/math]

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

%Datos
L=5; %Longitud de varilla
T=5; %tiempo final
q=1;
%Discretización temporal y espacial
dx=0.1; % Paso en espacio
N=L/dx; %Número de subintervalos
x=0:dx:L; % Vector de espacio
dt=dx/4; % Paso en tiempo
t=0:dt:T; % Vector de tiempos
% Matriz de coeficientes. 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/dx^2;
% Calculamos u0 condición inicial
u0=[zeros(1,1+N*3/5),3*ones(1,N*2/5)]';
% F=0 (no la consideramos)
sol(1,:)=u0'; % Guardamos la solucion
U=u0; % Inicializacion
% Calculamos U_n -> U_n+1
%Aplicamos método trapecio método del trapecio
for j=1:length(t)-1
U=(eye(N+1)+(dt*K)/2)\((eye(N+1)-(K*dt)/2)*U);
sol(j+1,:)= U'; % Guardamos solucion
end
%Dibujamos la solucion
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol);

Gráfica concentración / espacio-tiempo. Resolución por el método del trapecio.

3.1.2 Resolución por método de Euler explícito

%Datos
L=5; %Longitud de varilla
T=5; %tiempo final
q=1;
%Discretización temporal y espacial
dx=0.1; % Paso en espacio
N=L/dx; %Número de subintervalos
x=0:dx:L; % Vector de espacio
dt=dx/20; % Paso en tiempo
t=0:dt:T; % Vector de tiempos
% Matriz de coeficientes. 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/dx^2;
% Calculamos u0 condición inicial
u0=[zeros(1,1+N*3/5),3*ones(1,N*2/5)]';
% F=0 (no la consideramos)
sol(1,:)=u0'; % Guardamos la solucion
U=u0; % Inicializacion
% Calculamos U_n -> U_n+1
%Aplicamos método Euler explícito
for j=1:length(t)-1 
k1=-K*U; % Calculamos k1=h*f(tn,yn)
k2=-K*(U+dt*k1); % Calculamos k2=h*f(tn+h/2,yn+1/2*k1)
% Meto los datos en el vector
U=U+(dt/2)*(k1+k2);
sol(j+1,:)= U'; % Guardamos solucion
end
%Dibujamos la solucion
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol);

Gráfica concentración / espacio-tiempo. Resolución por el método de Euler explícito.

3.1.3 Resolución por método de Euler implícito

%Datos
L=5; %Longitud de varilla
T=5;
q=1;
%Discretización temporal y espacial
dx=0.1; % Paso en espacio
N=L/dx; %Número de subintervalos
h=L/N;
x=0:h:L; % Vector de espacio
dt=h/20; % Paso en tiempo
t=0:dt:T; % Vector de tiempos
% Matriz de coeficientes. 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 condición inicial
u0=[zeros(1,1+N*3/5),3*ones(1,N*2/5)]';
% F=0 (no la consideramos)
sol(1,:)=u0'; % Guardamos la solucion
U=u0; % Inicializacion
% Calculamos U_n -> U_n+1
% Aplicamos método implícito
for j=1:length(t)-1 
U=U-dt*K*U; 
sol(j+1,:)= U'; % Guardamos solucion
end
%Dibujamos la solucion
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol);

Gráfica concentración / espacio-tiempo. Resolución por el método de Euler implícito.

3.1.4 Resolución por método de Euler modificado

%Datos
L=5; %Longitud de varilla
T=5;
q=1;
%Discretización temporal y espacial
dx=0.1; % Paso en espacio
N=L/dx; %Número de subintervalos
h=L/N;
x=0:h:L; % Vector de espacio
dt=h/20; % Paso en tiempo
t=0:dt:T; % Vector de tiempos
% Matriz de coeficientes. 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 condición inicial
u0=[zeros(1,1+N*3/5),3*ones(1,N*2/5)]';
% F=0 (no la consideramos)
sol(1,:)=u0'; % Guardamos la solucion
U=u0; % Inicializacion
% Calculamos U_n -> U_n+1
for j=1:length(t)-1 
k1=-K*U; % Calculamos k1=h*f(tn,yn)
k2=-K*(U+dt*k1); % Calculamos k2=h*f(tn+h/2,yn+1/2*k1)
% Meto los datos en el vector
U=U+(dt/2)*(k1+k2);
sol(j+1,:)= U'; % Guardamos solucion
end
%Dibujamos la solucion
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol);

Gráfica concentración / espacio-tiempo. Resolución por el método de Euler modificado.

3.2 Conservación de la masa total de contaminante

A continuación demostraremos que la masa total de contaminante en el tubo se conserva a lo largo del tiempo, que se puede deducir a partir de la ecuación sin tener que resolverla. Integramos la ecuación y debido a las condiciones de frontera: [math]\frac d {dt} \int_0^5 u(x,t)\,dx= \int_0^5 u_{xx}(x,t)\,dx= u_x(x,t)|_0^5=0[/math]

%Datos
L=5; %Longitud de varilla
T=5;
q=1;
%Discretización temporal y espacial
N=500; %Número de subintervalos
h=L/N;
x=0:h:L; % Vector de espacio
dt=h/4; % Paso en tiempo
t=0:dt:T; % Vector de tiempos
% Matriz de coeficientes. 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 condición inicial
u0=[zeros(1,1+N*3/5),3*ones(1,N*2/5)]';
% F=0 (no la consideramos)
sol(1,:)=u0'; % Guardamos la solucion
U=u0; % Inicializacion
% Calculamos U_n -> U_n+1
for j=1:length(t)-1 
U=(eye(N+1)+dt*K)\(U+dt*F); 
sol(j+1,:)= U'; % Guardamos solucion
end
% Calculamos el area de la funcion en cada tiempo (area = cantidad de contaminante)
areas=zeros(1,length(t));
for i=1:length(t)
f=sol(i,:);
w=trapz(x,f);
area(i)=w;
end
%Dibujamos
plot(t,area)
xlabel('Tiempo')
ylabel('Cantidad de contaminante')

Conservación masa total de contaminante

Conservación masa total de contaminante

3.3 Concentración a lo largo del tiempo

La concentración de la sustancia contaminante en la varilla varía con el tiempo, pero en un cierto momento esta concentración se va homogeneizando, y deja de variar, entrando en un estado estacionario:: [math] u_t(x,t)- D*u_{xx}(x,t)=0 [/math] : [math] u_t(x,t)=0 \rightarrow D*u_{xx}(x,t)=0 [/math] : [math] u=c_1*x+c_2 [/math] Utizando Matlab comprobaremos que [math] c_1=0 [/math] [math] c_2=1.17 [/math] Para tiempos grandes la concentración se estabiliza en torno a 1.17 en toda la varilla, de forma que la solución se asemejará a una función lineal constante.

%Datos
L=5; %Longitud de varilla
T=50;
q=1;
%Discretización temporal y espacial
dx=0.1; % Paso en espacio
N=L/dx; %Número de subintervalos
h=L/N;
x=0:h:L; % Vector de espacio
dt=h/4; % Paso en tiempo
t=0:dt:T; % Vector de tiempos
% Matriz de coeficientes. 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 condición inicial
u0=[zeros(1,1+N*3/5),3*ones(1,N*2/5)]';
% F=0 (no la consideramos)
sol(1,:)=u0'; % Guardamos la solucion
U=u0; % Inicializacion
% Calculamos U_n -> U_n+1
%Aplicamos método trapecio método del trapecio
for j=1:length(t)-1
U=(eye(N+1)+(dt*K)/2)\((eye(N+1)-(K*dt)/2)*U);
sol(j+1,:)= U'; % Guardamos solucion
end
%Solución estacionaria
a=sol(length(t),:);
b=sol(1,:);
c=sol((length(t)-1)/50+1,:);
d=sol((length(t)-1)*2/50+1,:);
e=sol((length(t)-1)*10/50+1,:);
hold on
plot(x,a,'xk','linewidth',2);
plot(x,b,'-m','linewidth',1);
plot(x,c,'-g','linewidth',1);
plot(x,d,'-y','linewidth',1);
plot(x,e,'-b','linewidth',1);
legend('u(50s)','u(0s)','u(1s)','u(2s)','u(10s)');
xlabel('Posición en la varilla')
ylabel('Concentracion')
hold off

Soluciones a lo largo del tiempo

La distancia entre la solución estacionaria y las soluciones para los siguientes instantes: 0s, 1s, 2s y 10s; se puede apreciar en la imagen anterior.
Una vez analizado los datos, se demuestra que la sustancia contaminante se estabiliza en t=9.98s, con un error del 5%.

En el caso de que consideremos Δx'=Δx/10,el tiempo en el que se estabiliza la cantidad de sustancia contaminante será t=10s, dato que será más preciso debido a que el paso realizado en el programa es menor.

3.4 Nueva condición de frontera

Si colocamos un limpiador, que elimina la sustancia contaminante, en el extremo izquierdo de la varilla, varían las condiciones de contorno del problema. En este caso, en el extremo izquierdo, la condición inicial será de tipo Dirichlet ya que se esta variación esta provocada por un mecanismo físico externo que actúa sobre el extremo.
[math]\begin{cases} u_t(x,t)- D u_{xx}(x,t)= 0 ; x\in (0,5) \\ u(x,0)=u_0 , x\in (0,5) \\ u(0,t)=0 , u_x(5,t)=0 ; t\geq 0 \end{cases} [/math]

Dadas las nuevas condiciones, resolvemos el problema por el método de diferencias finitas con la ayuda de Matlab.

%Datos
L=5; %Longitud de varilla
T=5; %tiempo final
q=1;
%Discretización temporal y espacial
dx=0.1; % Paso en espacio
N=L/dx; %Número de subintervalos
x=0:dx:L;% Vector de espacio
dt=dx/4; % Paso en tiempo
t=0:dt:T; % Vector de tiempos
% Matriz de coeficientes. Aproximacion de -q*u_xx
K=2*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
K(N,N-1)=-2;
K=q*K/dx^2;
% Matriz F
F=zeros(N,1);
% Calculamos u0 condición inicial
u0=[zeros(1,N*3/5),3*ones(1,N*2/5)]';
sol(1,:)=[0,u0']; 
U=u0; % Inicializacion
% Calculamos U_n -> U_n+1
%Aplicamos método trapecio
for j=1:length(t)-1
U=(eye(N)+(dt*K))\(dt*F+U);
sol(j+1,:)= [0,U'];
end
%Dibujamos la solucion
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol);

Solución al variar las condiciones de contorno del problema.

Observando la gráfica, y por la presencia del limpiador en el extremo de la varilla, para tiempos grandes la cantidad de sustancia contaminante será nula en toda la varilla. En realidad, esto no es del todo cierto puesto que la sustancia no llega a eliminarse completamente, ya que el limpiador solo elimina la sustancia contaminante en el extremo izquierdo. Aunque la concentración en todos los puntos disminuye de manera exponencial, permanece asintótica a 0.

Momento en el que se estabiliza con la presencia del limpiador.

3.5 Resolución por el método de Fourier

Las condiciones de frontera del problema son homogéneas, por lo que podemos aplicar este método de cálculo, sin realizar ningún cambio. El primer paso es resolver el problema de autovalores asociado:: [math]\begin{cases} y''+λy=0 x \in (0,5) \\ y'(0)=0 \\ y'(5)=0 \end{cases} [/math] Realizando analíticamente el cálculo, obtenemos los siguientes autovalores y autofunciones:: [math] λ_k=\frac{k^2π^2}{25} → φ_k=sin\frac{kπ}{5}x , k=0,1,2,3... [/math]

Las condiciones de frontera e iniciales serían las siguientes en términos de la serie de Fourier::

[math] \sum_{k=0}^∞c_k(t)φ_k(x)=0 , \sum_{k=0}^∞a_k(t)φ_k(x)\begin{cases} 0, x≤3 \\ 3, x\gt3\end{cases} [/math]

Y ensayamos soluciones de la forma::

[math] u(x,t)=\sum_{k=0}^∞T_k(t)φ_k(x)=\sum_{k=0}^∞T_k(t)sin\frac{kπ}{5}x [/math]

Relacionando los resultados obtenidos en el problema de autovalores y el problema de difusión de la sustancia contaminante obtenemos lo siguiente:: [math]\begin{cases} T'_k(t)+μ_kT_k(t)=c_k=0 \\ T_k(0)=a_k \end{cases} [/math]

La solución de este último problema de valor inicial es::

[math]T_k(t)=a_ke^{-μ_kt}[/math]

De esta manera obtenemos la solución final:: [math] u(x,t)=\sum_{k=0}^∞a_ke^{-μ_kt}sin\frac{kπ}{5}x [/math]

Resolveremos a través de Matlab, con una serie de Fourier de 1,3,10 y 20 términos.

%Datos
L=5; %Longitud de varilla
T=0.5; %tiempo final
q=1;
%Discretización temporal y espacial
dx=0.1; % Paso en espacio
N=L/dx; %Número de subintervalos
x=0:dx:L;% Vector de espacio
dt=dx/4; % Paso en tiempo
t=0:dt:T; % Vector de tiempos
M=1; %Número de términos
sol=zeros(length(t),N+1);
%Aplicamos Fourier
for j=1:M
switch phi=sin((j*pi*x))
case x<3
a=(trapz(x,0*phi))/(trapz(x,phi.^2));
otherwise, 
a=(trapz(x,3*phi))/(trapz(x,phi.^2));
end
TT=a*exp(-(j^2)*(pi^2)*t/25);
sol=sol+TT'*phi;
end
%Dibujamos la solucion
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol);

Con un término en la serie de Fourier obtenemos la siguiente gráfica::

Método de Fourier con un término en la serie.

En la anterior gráfica apreciamos que la precisión es baja con un sólo término en la serie. En cambio si tomamos 20 términos en la serie, la aproximación es mejor, como veremos en la siguiente imagen.

Método de Fourier con 20 términos en la serie.

Si variamos las condiciones de frontera en el problema el problema de autovalores cambia igualmente:: [math]\begin{cases} y''+λy=0 x \in (0,5) \\ y'(0)=3 \\ y(5)=10sen(t) \end{cases} [/math] De este problema obtenemos, operando analíticamente, nuevos autovalores y autofunciones:: [math] λ_k=\frac{π^2(2k-1)^2}{100} → φ_k=sin\frac{(2k-1)π}{10}x , k=0,1,2,3... [/math] Siguiendo el mismo procedimiento que en el caso anterior, obtendremos la solución::

[math] u(x,t)=\sum_{k=0}^∞a_ke^{-μ_kt}sin\frac{(2k-1)π}{10}x+10sint+3x-15 [/math]