Difusión de una sustancia contaminante en un tubo largo. A14
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Difusión de una sustancia contaminante en un tubo largo. Grupo A 14 |
| Asignatura | Ecuaciones Diferenciales |
| Curso | Curso 2015-16 |
| Autores | Guillermo Alfaya Franklin 1546
Carlos Nieto Egido 1321 Jesús Pedroche Serrano 1309 Rubén Peláez Moreno 1584 |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Contenido
1 Planteamiento del probema
Según Flick para que el fenómeno de la difusión tenga lugar, la distribución espacial de moléculas no debe ser homogénea, debe existir una diferencia, o gradiente de concentración entre dos puntos del medio. Para ello suponemos que la concentración varia con la posición de la barra en el eje x (abscisas) denominaremos a u como la concentración en este caso de contaminante por cada posición del tubo. La ley de Flick afirma que la densidad de corriente de partículas es proporcional al gradiente de concentración F = −D x(∂u /∂x) Consideraremos que cada sección es muy pequeña respecto del tubo que será un objeto unidimensional de longitud en el eje x[0,L] siendo L=7m. La concentración en cada punto de la barra a lo largo del tiempo será u(x,t) La acumulación de partículas en la unidad de tiempo que se produce en el elemento de volumen S•dx es igual a la diferencia entre el flujo entrante uS, menos el flujo saliente u’S, es decir uS-u’S=(du/dx).S.dx La acumulación de partículas en la unidad de tiempo es S.dx.(du/dt) Igualando ambas expresiones y utilizando la Ley de Flick se obtiene d(D.(du/dx))/dx=du/dt Ecuación diferencial en derivadas parciales que describe el fenómeno de la difusión . Si el coeficiente de difusión D no depende de la concentración
(1/D).(du/dt)=(d^2u)/(dx^2)== (1/D).(du/dt)- (d^2u)/(dx^2)=0;
(1/D).(du/dt)=usubt==ut; (d^2u)/(dx^2)==Dusubxx==Duxx conclusión: ut − Duxx = 0.
2 Resolución del problema
Resolveremos el problema con el método de diferencias finitas. Este proporciona una aproximación de la solución numérica de un sistema continuo. Partimos el intervalo [0,7] definiendo h como anchura de paso. En segundo lugar, se aplica la ecuación diferencial a un nodo interior Xn, sustituyendo la derivada uxx por esta aproximación de segundo orden: [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] Teniendo en cuenta esta aproximación, y la notación [math] u_n(t) = u(x_n,t)[/math] , resulta el siguiente sistema de N+1 ecuaciones, obtenido al aplicar la ecuación diferencial para cada n:
\begin{array}{c}u'_n(t)+\frac{-u_{n-1}(t)+2u_n(t)-u_{n+1}(t)}{h^2}=0\end{array} n=1,2,3...,N-1
Teniendo en cuenta:
[math]u'_n(t)\simeq\frac{u_{n+1}(t)-u_{n-1}(t)}{h2}[/math]
Aplicamos las condiciones de contorno:
\begin{array}{c}u'_0(t)=0\\u'_n(t)=0\end{array} Resultando el sistema:
Sistema= \left\{\begin{matrix}\U'=-KU+F=0\\U(0)=U^0\end{matrix}\right. </math>
La matriz [math]U^0[/math] incluye las condiciones iniciales.
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]
2.1 Resolución por el método del 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 es
% Trapecio
clear all
clf
% Datos del problema
L=7;
T=7;
q=1;
% Discretizacion espacial
h=1/10;
N=70;
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+1 a partir de U_n con un bucle for
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);
xlabel('Coordenada x');
ylabel('Tiempo');
zlabel('Concentración de contaminante');
La gráfica representa la variación de la concentración en función del tiempo y de la posición.
2.2 Resolución por el método de Euler explícito
Se resuelve el sistema utilizando el Método de Euler Explícito 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 es
% Euler Explicito
clear all
clf
% Datos del problema
L=7;
T=7;
q=1;
% Discretizacion espacial
h=1/10;
N=70;
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
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+1 a partir de U_n con un bucle for
for j=1:length(t)-1 % Voy hasta la longitud de t menos 1 porque ya conozco un valor
U=U-dt*K*U; % EULER EXPLICITO
sol(j+1,:)= U'; % Guardamos solucion
end
% Dibujamos la solucion
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol);
xlabel('Coordenada x');
ylabel('Tiempo');
zlabel('Concentración de contaminante');2.3 Resolución por el método de Euler implícito
Se resuelve el sistema utilizando el Método de Euler Explícito tomando [math]\triangle t= \triangle x/4[/math] en [math]t\in\ [0,7][/math], con [math]\triangle x=0.1[/math]
% Euler Implícito
clear all
clf
% Datos del problema
L=7;
T=7;
q=1;
% Discretizacion espacial
h=1/10;
N=70;
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+1 a partir de U_n con un bucle for
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)\(U+dt*F); % EULER IMPLICITO
sol(j+1,:)= U'; % Guardamos solucion
end
% Dibujamos la solucion
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol);
xlabel('Coordenada x');
ylabel('Tiempo');
zlabel('Concentración de contaminante');2.4 Resolución por el método de Heun
Se resuelve el sistema utilizando el Método de Heun tomando [math]\triangle t= \triangle x/20[/math] en [math]t\in\ [0,7][/math], con [math]\triangle x=0.1[/math]
El código en MATLAB es
% Heun
clear all
% Datos del problema
L=7;
T=7;
q=1;
% Discretizacion espacial
h=1/10;
N=70;
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
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+1 a partir de U_n con un bucle for
for j=1:length(t)-1 % Voy hasta la longitud de t menos 1 porque ya conozco un valor
K1=K*U; % HEUN
K2=K*(U-dt*K*U);
U=U-dt*0.5*(K1+K2);
sol(j+1,:)= U'; % Guardamos solucion
end
% Dibujamos la solucion
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol);
xlabel('Coordenada x');
ylabel('Tiempo');
zlabel('Concentración de contaminante');3 Conservación de la masa
Al no dejar pasar el tubo masa por sus paredes, ni ser generada masa en el interior del mismo, la cantidad de masa de contaminante se conserva. Podemos deducir esto a partir de la ecuación u(x,t) integrando:
[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)|_0^L=0[/math]
Con el siguiente código de MATLAB podemos calcular la cantidad de contaminante integrando la concentración en la variable x mediante el método de Euler implícito.
%Cantidad de contaminante con Euler Implícito
clear all
clf
% Datos
L=7;
T=10;
q=1;
% Discretizacion espacial
h=0.1;
N=L/h;
x=0:h:L;
% Discretizacion temporal
dt=h/4; % Tamaño de paso
t=0:dt:T; % Vector de tiempos
% Matriz K
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;
% Definimos F
F=zeros(N+1,1);
% Calculamos u0
u0=[zeros(1,1+N*5/7),3*ones(1,N*2/7)]';
% Guardamos la solucion para t=0 en la matriz solucion
sol(1,:)=u0';
U=u0; % Inicializacion
% Calculamos U_n+1 a partir de U_n con un bucle for
for j=1:length(t)-1 % Hasta length(t)-1 porque ya tenemos el valor para t=0
U=(eye(N+1)+dt*K)\(U+dt*F); % EULER IMPLICITO
sol(j+1,:)= U'; % Guardamos la solucion
end
% Calculamos el area de la funcion para cada t (el area será la cantidad de contaminante)
contaminante=zeros(1,length(t));
for i=1:length(t)
f=sol(i,:);
area=trapz(x,f);
contaminante(i)=area;
end
% Obtenemos en la figura 1 la cantidad de contaminante con el tiempo,y en la
% figura 2 la concentracion en el punto medio de la barra a lo largo del tiempo
figure(1)
plot(t,contaminante)
xlabel('Tiempo')
ylabel('Cantidad de contaminante')
punto_medio=(length(x)-1)/2;
figure(2)
plot(t,sol(:,punto_medio),'g')En la gráfica se aprecia que la masa total de contaminante dentro del tubo no cambia. Vale 5.85
La segunda gráfica muestra cómo cambia la concentración de contaminante en el punto medio del tubo.
4 Concentración a lo largo del tiempo
Se puede comprobar en la siguiente gráfica que para tiempos grandes la cantidad de contaminante se irá acercando al valor estacionario de 1.17. Es decir, nuestra solución del problema debe parecerse a la función u(x,t)=1.17, es constante, luego no varía.
El código usado:
% Datos
L=7;
T=10;
q=1;
% Discretizacion espacial
h=0.1;
N=L/h;
x=0:h:L;
% Discretizacion temporal
dt=h/4; % Tamaño de paso
t=0:dt:T; % Vector de tiempos
% Matriz K
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;
% Definimos F
F=zeros(N+1,1);
% Calculamos u0
u0=[zeros(1,1+N*5/7),3*ones(1,N*2/7)]';
% Guardamos la solucion para t=0 en la matriz solucion
sol(1,:)=u0';
U=u0; % Inicializacion
% Calculamos U_n+1 a partir de U_n con un bucle for
for j=1:length(t)-1 % Hasta length(t)-1 porque ya tenemos el valor para t=0
U=(eye(N+1)+dt*K)\(U+dt*F); % EULER IMPLICITO
sol(j+1,:)= U'; % Guardamos la solucion
end
% Calculamos el area de la funcion para cada t (el area será la cantidad de contaminante)
contaminante=zeros(1,length(t));
for i=1:length(t)
f=sol(i,:);
area=trapz(x,f);
contaminante(i)=area;
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);
plot(x,b);
plot(x,c);
plot(x,d);
plot(x,e);
xlabel('Tubo')
ylabel('Concentracion')
hold off
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. Además las distancias entre la solución estacionaria y la solución u(x,t) para t=0,1,2,10 son:
dist10=0.3624
dist1=0.9650
dist2=0.7853
dist0=1.3880
Y el código Matlab usado es:
%distancias
dist0=(sqrt(sum((sol(1,:)-1.17).^2)/(N+1)));
dist1=(sqrt(sum((sol(41,:)-1.17).^2)/(N+1)));
dist2=(sqrt(sum((sol(81,:)-1.17).^2)/(N+1)));
dist10=(sqrt(sum((sol(401,:)-1.17).^2)/(N+1)));
Para averiguar el tiempo que tarda la concentración en alcanzar el estado estacionario con un error del 5% calculamos la media cuadrática de las distancias desde cada punto de la función al estado estacionario.
% El error del 5%
limite=0.05*8.25;
% Buscamos para que tiempo se alcanza ese valor
for i=1:length(t)
dist=(sqrt(sum((sol(i,:)-1.17).^2)/(N+1)));
if dist<limite
tiempo=(i-1)*dt;
break
end
end
Se obtiene un tiempo de 7.300 segundos.
Al dividir \triangle x por 10 sale un tiempo de 8.353 segundos.
5 Nuevas condiciones de frontera
5.1 Aplicación de un limpiador
Pasado un tiempo T = 1,se coloca en el extremo izquierdo un limpiador que elimina el contaminante en ese extremo,modificando las condiciones de frontera del problema generando que una condición de frontera tipo Neumann pase a ser una condición tipo Dirichlet, obteniendo el siguiente sistema:
- [math] (P)\left\{\begin{matrix}u_t-u_{xx}=0 ,x\in\ (0,7), t\gt0\\u_x(0,t)=0, t\lt1\\u(0,t)=0, t\gt1\\u_x(7,t)=0, t\gt0\\u(x,0)=u_0, x\in\ (0,7)\end{matrix}\right. [/math]
Suponiendo que en el instante inicial se verifica :[math] u(x,0)=\left\{\begin{matrix}0, x≤5\\3, x\gt5\end{matrix}\right. [/math]
Se comprueba que está bien propuesto analizando la existencia de solución, la unicidad de ésta y su estabilidad,
usaremos el 'método de diferencias finitas para llevar a cabo una aproximación del problema, siendo el código de Matlab obtenido:
clear all
% Datos
L=7;
T=1;
q=1;
% Discretizacion espacial
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L; % Vector de nodos interiores
% Discretizacion temporal
dt=h/4; % Tamaño de paso
t=0:dt:T; % Vector de tiempos
% Matriz K para t<=1
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2)=-2; % Condicion Neumann
K(N+1,N)=-2; % Condicion Neumann
K=q*K/h^2;
% Matriz K para t>1
K2=2*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
K2(N,N-1)=-2;
K2=K2/h^2;
% Definimos F
F=zeros(N+1,1);
F2=zeros(N,1);
% Calculamos u0
u0=[zeros(1,1+N*5/7),3*ones(1,N*2/7)]';
% Guardamos la solucion para t=0 en la matriz solucion
sol(1,:)=u0';
U=u0; % Inicializacion
% Calculamos U_n+1 a partir de U_n con un bucle for
for j=1:((length(t)-1)/7)
U=(eye(N+1)+dt*K)\(U+dt*F); % Euler Implicito
sol(j+1,:)= U'; % Guardamos la solucion
end
U(1)=[];
for j=((length(t)-1)/7):length(t)-1
U=(eye(N)+dt*K2)\(U+dt*F2); % Euler Implicito
sol(j+1,:)= [0,U']; % Guardamos la solucion
end
% Dibujamos la solucion
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol);
La gráfica obtenida sería:
El valor estacionario de la concentración en el tubo será de cero, debido a que se ha colocado un limpiador en uno de sus extremos que irá eliminando el contaminante.
Para averiguar el tiempo que tarda la concentración en alcanzar este otro estado estacionario con un error del 5% vamos a calcular la media cuadrática de las concentraciones y ver cual coincide con el 5% de 5.85/7.
clear all
% Datos
L=7;
T=70;
q=1;
% Discretizacion espacial
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L; % Vector de nodos interiores
% Discretizacion temporal
dt=h/4; % Tamaño de paso
t=0:dt:T; % Vector de tiempos
% Matriz K para t<=1
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2)=-2; % Condicion Neumann
K(N+1,N)=-2; % Condicion Neumann
K=q*K/h^2;
% Matriz K para t>1
K2=2*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
K2(N,N-1)=-2;
K2=K2/h^2;
% Definimos F
F=zeros(N+1,1);
F2=zeros(N,1);
% Calculamos u0
u0=[zeros(1,1+N*5/7),3*ones(1,N*2/7)]';
% Guardamos la solucion para t=0 en la matriz solucion
sol(1,:)=u0';
U=u0; % Inicializacion
% Calculamos U_n+1 a partir de U_n con un bucle for
for j=1:(1/dt)
U=(eye(N+1)+dt*K)\(U+dt*F); % Euler Implicito
sol(j+1,:)= U'; % Guardamos la solucion
end
U(1)=[];
for j=(1/dt):length(t)-1
U=(eye(N)+dt*K2)\(U+dt*F2); % Euler Implicito
sol(j+1,:)= [0,U']; % Guardamos la solucion
end
% Calculamos cual es el error del 5%
limite=0.05*(5.85/7);
% Buscamos para que tiempo se alcanza ese valor
for i=1:length(t)
dist=(sqrt(sum(sol(i,:).^2)/(N+1)));
if dist<limite
tiempo=(i-1)*dt
break
end
endObteniendo un tiempo de 65.775.
6 Variación de las condiciones de frontera tipo Neumann
Las nuevas condiciones neumann no son homogéneas. tendremos que homogeneizarlas con el siguiente cambio de variable :
- [math] u(x,t)=w(x,t)+2a(t)x^2+b(t)x [/math]
el nuevo problema quedaría de la siguiente manera :
- [math] (P)\left\{\begin{matrix}\\w_t-w_{xx}=-3/7 ,x\in\ (0,7), t\gt0\\w_x(0,t)=0, t\gt0\\w_x(7,t)=-3, t\gt0\\w(x,0)=u_0, x\in\ (0,7)\end{matrix}\right. [/math]
Suponiendo que en el instante inicial se verifica :[math] w(x,0)=\left\{\begin{matrix}\\0, x≤5\\3, x\gt5\end{matrix}\right. [/math]
Podemos observar que es un problema bien propuesto analizando su existencia de solución, unicidad y estabilidad con respecto al dato inicial.
ahora ya se trata de un sistema con condiciones de frontera homogéneas, el problema de autovalores asociado es:
- [math] (PA)\left\{\begin{matrix}\\x''+μx=0 ,x\in\ (0,7)\\x'(0)=0\\x'(7)=0\end{matrix}\right. [/math]
Y los autovalores y autofunciones salvo una constante son:
[math]μ_k=\frac{\pi^2k^2}{49}[/math]: [math]φ_k=cos(kπ/7)x[/math] k=0,1,2,3...
Se aproxima el dato inicial y el segundo miembro, se elige Q y se toma:
[math]\displaystyle\sum_{k=0}^Q c_kφ_k(x)=-3/7[/math]:
[math]\displaystyle\sum_{k=0}^Q a_kφ_k(x)=\left\{\begin{matrix}\\0, x≤5\\3, x\gt5\end{matrix}\right.[/math]
donde se deduce que [math]c_k=\frac{\displaystyle\int_{0}^{7} -3/7cos(kπ/7)x\, dx}{\displaystyle\int_{0}^{7} cos^2(kπ/7)x\, dx}[/math]
[math]a_k=\frac{\displaystyle\int_{0}^{7} u_0cos(kπ/7)x\, dx}{\displaystyle\int_{0}^{7} cos^2(kπ/7)x\, dx}[/math]
El problema aproximado queda de la forma:
- [math] (P)\left\{\begin{matrix}\\w_t-w_{xx}=\displaystyle\sum_{k=0}^Q c_kφ_k(x) ,x\in\ (0,7), t\gt0\\w_x(0,t)=0, t\gt0\\w_x(7,t)=0, t\gt0\\w(x,0)=w_0=\displaystyle\sum_{k=0}^Q a_kφ_k(x), x\in\ (0,7)\end{matrix}\right. [/math]
Luego la solución aproximada que obtenemos será:
[math]w(x,t)=\displaystyle\sum_{k=0}^Q T_k(t)φ_k(x)[/math]
Sustituyendo en la ecuación de difusión de una sustancia contaminante e igualando los coeficicientes de Fourier se obtiene el sistema:
[math]\left\{\begin{matrix}\\T'_k(t)+μ_kT_k(t)=c_k\\T_k(0)=a_k\end{matrix}\right.[/math]
Resolviendo el problema de valor inicial, se deduce que [math] T_k(t)=Ce^{-μ_k t}[/math]
La solución aproximada queda de la forma [math]w(x,t)=\displaystyle\sum_{k=0}^Q Ce^{-μ_k t} cos(kπ/7)x[/math]
El valor de la constante C depende de la condición inicial,que es la que condiciona el valor de ak.
por ultimo solo nos quedaria deshacer el cambio de variable.






