Difusión de una sustancia contaminante - Grupo 11B
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Difusión de sustancias contaminantes. Grupo 11B |
| Asignatura | Ecuaciones Diferenciales |
| Curso | Curso 2016-17 |
| Autores | Alejandro Barnuevo Javier Carrillo Pablo Escudero Alejandro Gaitón Alejandro González Pablo Vidal |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Contenido
1 Deducción de la ecuación.
En primer lugar, sabemos que para que se produzca difusión, es necesario que exista un gradiente de concentración o diferencia entre dos puntos del medio. Definimos la variable u como la concentración de contaminante en cada posición del tubo (que consideramos unidimensional y con sección despreciable, de longitud L=8m), la cual varía a lo largo del eje de abscisas (eje x), por lo que el intervalo en el que trabajaremos será [0,L]. Teniendo en cuenta la ley de Flick, sabemos que la densidad de corriente de partículas es proporcional al gradiente de concentración F = −D x(∂u /∂x). Añadiendo la variable del tiempo, la concentración de contaminante en cada posición y en cada instante la expresaremos por u(x,t). La acumulación de partículas en cada instante que se produce por cada unidad de volumen S.dx (S es la sección) 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 por unidad de tiempo es S.dx.(du/dt) Si igualamos ambas expresiones y aplicamos la Ley de Flick nos queda d(D.(du/dx))/dx=du/dt. Obtenemos una ecuación diferencial en derivadas parciales que describe el fenómeno de la difusión. Todo esto es válido si el coeficiente de difusión D no depende de la concentración.
2 Planteamiento del sistema.
3 Método del trapecio.
Resolvemos por el método de las líneas la ecuación diferencial con dos condiciones Newman que vienen dadas por la ausencia de flujo en los extremos del tubo. Empleamos también la condición inicial: u(x,0) = 2, si x ≤ 4, 4, si x > 4 y resolvemos matricialmente el sistema de N+1 incógnitas.
% Datos
L=8;
T=7;
q=1;
% Discretizacion espacial
h=0.15;
N=round(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; % Condicion Neumann
K(N+1,N)=-2; % Condicion Neumann
K=q*K/h^2;
% Definimos F
F=zeros(N+1,1);
% Calculamos u0
u0=[2*ones(1,2+N*4/8),4*ones(1,N*4/8)]';
% 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 % 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);
4 Conservación de la masa.
Integramos en todos los instantes la cantidad de contaminante que hay en cada punto de la barra, lo que nos da la masa total de contaminante en cada instante. A partir de la ecuación u(x,t)se deduce que la masa total del contaminante se conserva a lo largo del tiempo. Dicha masa es 23.550, hallada según el programa:
% Datos
L=8;
T=12;
q=1;
% Discretizacion espacial
h=0.15;
N=round(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; % Condicion Neumann
K(N+1,N)=-2; % Condicion Neumann
K=q*K/h^2;
% Definimos F
F=zeros(N+1,1);
% Calculamos u0
u0=[2*ones(1,2+N*4/8),4*ones(1,N*4/8)]';
% 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))/2;
figure(2)
plot(t,sol(:,punto_medio))4.1 Evidencias
5 Aplicación de los métodos
Resolvemos por los métodos de Euler explícito, ímplicito y de Heun con las condiciones inicilales del ejercicio anterior.
5.1 Euler explícito
% Euler Explicito
% Datos del problema
L=8;
T=7;
q=1;
% Discretizacion espacial
h=0.15;
N=round(L/h);
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);%tienes N+1 incognitas
K(1,2)=-2;
K(N+1,N)=-2;
K=q*K/h^2;
% Calculamos u0
u0=[2*ones(1,2+N*4/8),4*ones(1,N*4/8)]';
% 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');
5.2 Euler implícito
% Euler Implícito
clear all
clf
% Datos del problema
L=8;
T=7;
q=1;
% Discretizacion espacial
h=0.15;
N=round(L/h);
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=[2*ones(1,2+N*4/8),4*ones(1,N*4/8)]';
% 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');
5.3 Heun
% Heun
clear all
% Datos del problema
L=8;
T=7;
q=1;
% Discretizacion espacial
h=0.15;
N=round(L/h);
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=[2*ones(1,2+N*4/8),3*ones(1,N*4/8)]';
% 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');
5.4 Comparación de resultados
Observando las gráficas se aprecia que apenas hay diferencia.
6 Concentración casi constante. Consecuencias.
En la imagen apreciamos como a tiempos grandes (50 segundos) el contaminante se distribuye de forma homogénea estabilizándose en 2.9438, la cantidad de contaminante total dividida la longitud del tubo, 8.
Por lo que podemos decir que el flujo es 0(muy próximo) y por la tanto la tanto la solución es homogénea siendo en cada instante y posición 2.9438. La función de la solución sería u(x,t)= 2.9438 cuando el tiempo es grande. Este estado estacionario se da cuando el flujo se aproxime a 0.
Para demostrar esta solución en régimen estacionario hemos usado los tiempos 0,1,2,10 y hemos calculado la media cuadrada entre el valor de estacionamiento y los diferentes valores de la concentración en cada posición:
Distancia para 0 segundos= 0.999499
Distancia para 1 segundos = 0.776147
Distancia para 2 segundos = 0.657469
Distancia para 10 segundos = 0.192049
% 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
L=8;
T=10;
q=1;
% Discretizacion espacial
h=0.15;
N=round(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; % Condicion Neumann
K(N+1,N)=-2; % Condicion Neumann
K=q*K/h^2;
% Definimos F
F=zeros(N+1,1);
% Calculamos u0
u0=[2*ones(1,2+N*4/8),4*ones(1,N*4/8)]';
% 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 % 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);
dist0=(sqrt(sum((sol(1,:)-(23.550/8)).^2)/(N+1)))
dist1=(sqrt(sum((sol(28,:)-(23.550/8)).^2)/(N+1)))
dist2=(sqrt(sum((sol(56,:)-(23.550/8)).^2)/(N+1)))
dist10=(sqrt(sum((sol(267,:)-(23.550/8)).^2)/(N+1)))
7 Estado estacionario. Conclusiones.
Hemos calculado el tiempo en el cual podemos considerar que se ha alcanzado el estado estacionario (es decir cuando la distancia que hemos calculado anteriormente es menor de un valor establecido). El valor establecido es el 5% de la concentración media, es decir 0.147187 .
El tiempo que nos ha dado es 11.7375 segundos.
Para ello utilizamos el siguiente programa con tiempo 50.
% 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
L=8;
T=50;
q=1;
% Discretizacion espacial
h=0.15;
N=round(L/h);
x=0:h:L;
% Discretizacion temporal
dt=0.15/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; % Condicion Neumann
K(N+1,N)=-2; % Condicion Neumann
K=q*K/h^2;
% Definimos F
F=zeros(N+1,1);
% Calculamos u0
u0=[2*ones(1,2+N*4/8),4*ones(1,N*4/8)]';
% 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 % 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
% Calculamos cual es el error del 5%
limite=0.05*(23.550/8);
% Buscamos para que tiempo se alcanza ese valor
for i=1:length(t)
dist=(sqrt(sum((sol(i,:)-(23.550/8)).^2)/(N+1)));
if dist<limite
tiempo=(i-1)*dt
break
end
endHemos querido probar si el ∆x influía y para ello lo hemos multiplicado por 10. Hemos utilizado el anterior programa y solo hemos cambiado el ∆x .
% 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
L=8;
T=400;
q=1;
% Discretizacion espacial
h=1.5;
N=round(L/h);
x=[0,1.5,3,4.5,6,7.5,8];
% Discretizacion temporal
dt=0.15/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; % Condicion Neumann
K(N+1,N)=-2; % Condicion Neumann
K=q*K/h^2;
% Definimos F
F=zeros(N+1,1);
% Calculamos u0
u0=[2*ones(1,2+N*4/8),4*ones(1,N*4/8)]';
% 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 % 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
% Calculamos cual es el error del 5%
limite=0.05*(23.550/8);
% Buscamos para que tiempo se alcanza ese valor
for i=1:length(t)
dist=(sqrt(sum((sol(i,:)-(23.550/8)).^2)/(N+1)));
if dist<limite
tiempo=(i-1)*dt
break
end
end
Y cambia el resultado ya que la distancia en ningún momento es menor del 5%, es decir de 0.147187. Por ello, aunque se podría suponer que no se alcanza el estado estacionario en realidad sí lo hace, pues lo que ocurre es que el error es mayor ya que la distancia que nos da es siempre constante e igual a 0.34375.
7.1 Aplicación de un limpiador
En el segundo 1 se sitúa en el extremo 0 un filtro que limpia dicha posición. Hasta ese instante el problema es idéntico a los anteriores, pero a partir de entonces solo la concentración será nula en el extremo inicial. El nuevo valor estacionario es 0 ya que se va eliminando contaminante.
Para obtener ambos gráficos se ha utilizado el siguiente programa:
% u_t-qu_xx=0, x en (0,L)
% u_x(0,t)=0 Condicion Neumann para t<=1
% u(0,t)=0 Condicion Dirichlet para t>1
% u_x(L,t)=0 Condicion Neumann
% u(x,0)=u0(x)
% Datos
L=8;
T=10;
q=1;
% Discretizacion espacial
h=0.15;
N=round(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=[2*ones(1,2+N*4/8),4*ones(1,N*4/8)]';
% 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:round(((length(t)-1)/T))
U=(eye(N+1)+dt*K)\(U+dt*F); % Euler Implicito
sol(j+1,:)= U'; % Guardamos la solucion
end
U(1)=[];
for j=round(((length(t)-1)/T)):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);
Además este estado estacionario (que es 0) se alcanzará con menos del 5% de error en 78.1875 segundos.
Para ello utilizamos este programa:
% u_t-qu_xx=0, x en (0,L)
% u_x(0,t)=0 Condicion Neumann para t<=1
% u(0,t)=0 Condicion Dirichlet para t>1
% u_x(L,t)=0 Condicion Neumann
% u(x,0)=u0(x)
% Datos
L=8;
T=100;
q=1;
% Discretizacion espacial
h=0.15;
N=round(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=[2*ones(1,2+N*4/8),4*ones(1,N*4/8)]';
% 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:round(((length(t)-1)/T))
U=(eye(N+1)+dt*K)\(U+dt*F); % Euler Implicito
sol(j+1,:)= U'; % Guardamos la solucion
end
U(1)=[];
for j=round(((length(t)-1)/T)):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);
% Calculamos cual es el error del 5%
limite=0.05*(23.550/8);
% 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
end
