Difusión de una sustancia contaminante - Grupo 11B

De MateWiki
Saltar a: navegación, buscar
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


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.


Planteamiento.png

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.

 1 % Datos
 2   L=8;
 3 T=7;
 4  q=1;
 5  
 6  % Discretizacion espacial
 7  h=0.15;
 8  N=round(L/h);
 9  x=0:h:L; 
10  
11  % Discretizacion temporal
12  dt=h/4; % Tamaño de paso
13  t=0:dt:T; % Vector de tiempos
14  
15  % Matriz K
16  K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
17  K(1,2)=-2; % Condicion Neumann
18  K(N+1,N)=-2; % Condicion Neumann
19  K=q*K/h^2;
20  
21  % Definimos F
22  F=zeros(N+1,1);
23   % Calculamos u0
24   u0=[2*ones(1,2+N*4/8),4*ones(1,N*4/8)]';
25  % Guardamos la solucion para t=0 en la matriz solucion
26  sol(1,:)=u0'; 
27  U=u0; % Inicializacion
28   % Calculamos U_n+1 a partir de U_n con un bucle for
29 for j=1:length(t)-1 % Voy hasta la longitud de t menos 1 porque ya conozco un valor
30  U=(eye(N+1)+(dt*K)/2)\((eye(N+1)-(K*dt)/2)*U); % Trapecio
31  sol(j+1,:)= U'; % Guardamos solucion
32  end
33   % Dibujamos la solucion
34  [xx,tt]=meshgrid(x,t);
35  surf(xx,tt,sol);


Concentración/espacio-tiempo

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:

 1 % Datos
 2 L=8;
 3 T=12;
 4 q=1;
 5 % Discretizacion espacial
 6 h=0.15;
 7 N=round(L/h);
 8 x=0:h:L; 
 9 % Discretizacion temporal
10 dt=h/4; % Tamaño de paso
11 t=0:dt:T; % Vector de tiempos
12 % Matriz K
13 K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
14 K(1,2)=-2; % Condicion Neumann
15 K(N+1,N)=-2; % Condicion Neumann
16 K=q*K/h^2;
17 % Definimos F
18 F=zeros(N+1,1);
19 % Calculamos u0
20   u0=[2*ones(1,2+N*4/8),4*ones(1,N*4/8)]';
21 % Guardamos la solucion para t=0 en la matriz solucion
22 sol(1,:)=u0'; 
23 U=u0; % Inicializacion
24 	
25 % Calculamos U_n+1 a partir de U_n con un bucle for
26 for j=1:length(t)-1 % Hasta length(t)-1 porque ya tenemos el valor para t=0
27 U=(eye(N+1)+dt*K)\(U+dt*F); % Euler Implicito
28 sol(j+1,:)= U'; % Guardamos la solucion
29 end
30 
31 % Calculamos el area de la funcion para cada t (el area será la cantidad de contaminante)
32 contaminante=zeros(1,length(t));
33 for i=1:length(t)
34 f=sol(i,:);
35 area=trapz(x,f);
36 contaminante(i)=area;
37 end
38 % Obtenemos en la figura 1 la cantidad de contaminante con el tiempo,y en la
39 % figura 2 la concentracion en el punto medio de la barra a lo largo del tiempo
40 figure(1)
41 plot(t,contaminante)
42 xlabel('Tiempo')
43 ylabel('Cantidad de contaminante')
44 punto_medio=(length(x))/2;
45 figure(2)
46 plot(t,sol(:,punto_medio))

4.1 Evidencias

"La masa es constante en el tiempo"
"Variación de la masa en el punto medio"


















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

 1 % Euler Explicito
 2   % Datos del problema
 3   L=8;
 4   T=7;
 5   q=1;
 6   % Discretizacion espacial
 7   h=0.15;
 8  N=round(L/h);
 9  x=0:h:L;
10  % Aproximacion de -q*u_xx
11  K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);%tienes N+1 incognitas
12  K(1,2)=-2; 
13  K(N+1,N)=-2;
14  K=q*K/h^2;
15  % Calculamos u0
16  u0=[2*ones(1,2+N*4/8),4*ones(1,N*4/8)]';
17  % Discretizacion temporal
18  dt=h/20; % Paso
19  t=0:dt:T; % Vector de tiempos
20  % Definimos F
21  F=zeros(N+1,1); 
22  % Guardamos la solucion
23  sol(1,:)=u0'; 
24  U=u0; % Inicializacion
25  % Calculamos U_n+1 a partir de U_n con un bucle for
26  for j=1:length(t)-1 % Voy hasta la longitud de t menos 1 porque ya conozco un valor
27  U=U-dt*K*U; % EULER EXPLICITO
28  sol(j+1,:)= U'; % Guardamos solucion
29  end
30  % Dibujamos la solucion
31  [xx,tt]=meshgrid(x,t);
32  surf(xx,tt,sol);
33  xlabel('Coordenada x');
34  ylabel('Tiempo');
35  zlabel('Concentración de contaminante');


5.2 Euler implícito

 1 % Euler Implícito
 2   clear all
 3   clf
 4   % Datos del problema
 5   L=8;
 6   T=7;
 7   q=1;
 8   % Discretizacion espacial
 9   h=0.15;
10  N=round(L/h);
11  x=0:h:L; 
12  % Aproximacion de -q*u_xx
13  K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
14  K(1,2)=-2;
15  K(N+1,N)=-2;
16  K=q*K/h^2;
17  % Calculamos u0
18  u0=[2*ones(1,2+N*4/8),4*ones(1,N*4/8)]';
19  % Discretizacion temporal
20  dt=h/4; % Paso
21  t=0:dt:T; % Vector de tiempos
22  % Definimos F
23  F=zeros(N+1,1); 
24  % Guardamos la solucion
25  sol(1,:)=u0'; 
26  U=u0; % Inicializacion
27  % Calculamos U_n+1 a partir de U_n con un bucle for
28  for j=1:length(t)-1 % Voy hasta la longitud de t menos 1 porque ya conozco un valor
29  U=(eye(N+1)+dt*K)\(U+dt*F); % EULER IMPLICITO
30  sol(j+1,:)= U'; % Guardamos solucion
31  end
32   % Dibujamos la solucion
33  [xx,tt]=meshgrid(x,t);
34  surf(xx,tt,sol);
35  xlabel('Coordenada x');
36  ylabel('Tiempo');
37  zlabel('Concentración de contaminante');


5.3 Heun

 1 % Heun
 2   clear all
 3   
 4   % Datos del problema
 5   L=8;
 6   T=7;
 7   q=1;
 8   % Discretizacion espacial
 9   h=0.15;
10   N=round(L/h);
11  x=0:h:L;
12  % Aproximacion de -q*u_xx
13  K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
14  K(1,2)=-2;
15  K(N+1,N)=-2;
16  K=q*K/h^2;
17  % Calculamos u0
18  u0=[2*ones(1,2+N*4/8),3*ones(1,N*4/8)]';
19  % Discretizacion temporal
20  dt=h/20; % Paso
21  t=0:dt:T; % Vector de tiempos
22  % Definimos F
23  F=zeros(N+1,1); 
24  % Guardamos la solucion
25  sol(1,:)=u0'; 
26  U=u0; % Inicializacion
27  % Calculamos U_n+1 a partir de U_n con un bucle for
28  for j=1:length(t)-1 % Voy hasta la longitud de t menos 1 porque ya conozco un valor
29  K1=K*U; % HEUN
30  K2=K*(U-dt*K*U);
31  U=U-dt*0.5*(K1+K2);
32  sol(j+1,:)= U'; % Guardamos solucion
33  end
34  % Dibujamos la solucion
35  [xx,tt]=meshgrid(x,t);
36  surf(xx,tt,sol);
37  xlabel('Coordenada x');
38  ylabel('Tiempo');
39  zlabel('Concentración de contaminante');


5.4 Comparación de resultados

Heun
Euler explícito
Euler implícito


                                       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.

Régimen estacionario

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



 1 % u_t-qu_xx=0, x en (0,L)
 2 % u_x(0,t)=0 Condicion Neumann
 3 % u_x(L,t)=0 Condicion Neumann
 4 % u(x,0)=u0(x)
 5 % Datos
 6 L=8;
 7 T=10;
 8 q=1;
 9 % Discretizacion espacial
10 h=0.15;
11 N=round(L/h);
12 x=0:h:L; 
13 
14 % Discretizacion temporal
15 dt=h/4; % Tamaño de paso
16 t=0:dt:T; % Vector de tiempos
17 
18 % Matriz K
19 K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
20 K(1,2)=-2; % Condicion Neumann
21 K(N+1,N)=-2; % Condicion Neumann
22 K=q*K/h^2;
23 
24 % Definimos F
25 F=zeros(N+1,1);
26 
27   % Calculamos u0
28   u0=[2*ones(1,2+N*4/8),4*ones(1,N*4/8)]';
29 % Guardamos la solucion para t=0 en la matriz solucion
30 sol(1,:)=u0'; 
31 U=u0; % Inicializacion
32 
33 % Calculamos U_n+1 a partir de U_n con un bucle for
34 for j=1:length(t)-1 % Voy hasta la longitud de t menos 1 porque ya conozco un valor
35 U=(eye(N+1)+(dt*K)/2)\((eye(N+1)-(K*dt)/2)*U); % Trapecio
36 sol(j+1,:)= U'; % Guardamos solucion
37 end
38  % Dibujamos la solucion
39 [xx,tt]=meshgrid(x,t);
40 surf(xx,tt,sol);
41 
42 dist0=(sqrt(sum((sol(1,:)-(23.550/8)).^2)/(N+1)))
43 dist1=(sqrt(sum((sol(28,:)-(23.550/8)).^2)/(N+1)))
44 dist2=(sqrt(sum((sol(56,:)-(23.550/8)).^2)/(N+1)))
45 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.

 1 % u_t-qu_xx=0, x en (0,L)
 2 % u_x(0,t)=0 Condicion Neumann
 3 % u_x(L,t)=0 Condicion Neumann
 4 % u(x,0)=u0(x)
 5 
 6 % Datos
 7 L=8;
 8 T=50;
 9 q=1;
10 
11 % Discretizacion espacial
12  h=0.15;
13  N=round(L/h);
14  x=0:h:L; 
15 % Discretizacion temporal
16 dt=0.15/4; % Tamaño de paso
17 t=0:dt:T; % Vector de tiempos
18 
19 % Matriz K
20 K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
21 K(1,2)=-2; % Condicion Neumann
22 K(N+1,N)=-2; % Condicion Neumann
23 K=q*K/h^2;
24 
25 % Definimos F
26 F=zeros(N+1,1);
27 
28 % Calculamos u0
29 
30 u0=[2*ones(1,2+N*4/8),4*ones(1,N*4/8)]';
31 % Guardamos la solucion para t=0 en la matriz solucion
32 sol(1,:)=u0'; 
33 U=u0; % Inicializacion
34 
35 % Calculamos U_n+1 a partir de U_n con un bucle for
36 for j=1:length(t)-1 % Voy hasta la longitud de t menos 1 porque ya conozco un valor
37 U=(eye(N+1)+(dt*K)/2)\((eye(N+1)-(K*dt)/2)*U); % Trapecio
38 sol(j+1,:)= U'; % Guardamos solucion
39 end
40 % Calculamos cual es el error del 5%
41 limite=0.05*(23.550/8);
42 % Buscamos para que tiempo se alcanza ese valor
43 for i=1:length(t)
44 dist=(sqrt(sum((sol(i,:)-(23.550/8)).^2)/(N+1)));
45 if dist<limite
46 tiempo=(i-1)*dt
47 break
48 end
49 end

Hemos 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 .

 1 % u_t-qu_xx=0, x en (0,L)
 2 % u_x(0,t)=0 Condicion Neumann
 3 % u_x(L,t)=0 Condicion Neumann
 4 % u(x,0)=u0(x)
 5 
 6 % Datos
 7 L=8;
 8 T=400;
 9 q=1;
10 
11 % Discretizacion espacial
12  h=1.5;
13  N=round(L/h);
14 x=[0,1.5,3,4.5,6,7.5,8]; 
15 % Discretizacion temporal
16 dt=0.15/4; % Tamaño de paso
17 t=0:dt:T; % Vector de tiempos
18 
19 % Matriz K
20 K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
21 K(1,2)=-2; % Condicion Neumann
22 K(N+1,N)=-2; % Condicion Neumann
23 K=q*K/h^2;
24 
25 % Definimos F
26 F=zeros(N+1,1);
27 
28 % Calculamos u0
29 
30 u0=[2*ones(1,2+N*4/8),4*ones(1,N*4/8)]';
31 % Guardamos la solucion para t=0 en la matriz solucion
32 sol(1,:)=u0'; 
33 U=u0; % Inicializacion
34 
35 % Calculamos U_n+1 a partir de U_n con un bucle for
36 for j=1:length(t)-1 % Voy hasta la longitud de t menos 1 porque ya conozco un valor
37 U=(eye(N+1)+(dt*K)/2)\((eye(N+1)-(K*dt)/2)*U); % Trapecio
38 sol(j+1,:)= U'; % Guardamos solucion
39 end
40 % Calculamos cual es el error del 5%
41 limite=0.05*(23.550/8);
42 % Buscamos para que tiempo se alcanza ese valor
43 for i=1:length(t)
44 dist=(sqrt(sum((sol(i,:)-(23.550/8)).^2)/(N+1)));
45 if dist<limite
46 tiempo=(i-1)*dt
47 break
48 end
49 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.

Se aprecia que en la izquierda la cantidad de contaminante es 0
Cuando el tiempo es muy grande la concentración total desciende a 0

Para obtener ambos gráficos se ha utilizado el siguiente programa:

 1 % u_t-qu_xx=0, x en (0,L)
 2 % u_x(0,t)=0 Condicion Neumann para t<=1
 3 % u(0,t)=0 Condicion Dirichlet para t>1
 4 % u_x(L,t)=0 Condicion Neumann
 5 % u(x,0)=u0(x)
 6 
 7 % Datos
 8 L=8;
 9 T=10;
10 q=1;
11 
12 % Discretizacion espacial
13 h=0.15;
14 N=round(L/h);
15 x=0:h:L; 
16 xint=h:h:L; % Vector de nodos interiores
17 
18 % Discretizacion temporal
19 dt=h/4; % Tamaño de paso
20 t=0:dt:T; % Vector de tiempos
21 
22 % Matriz K para t<=1
23 K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
24 K(1,2)=-2; % Condicion Neumann
25 K(N+1,N)=-2; % Condicion Neumann
26 K=q*K/h^2;
27 % Matriz K para t>1
28 K2=2*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
29 K2(N,N-1)=-2;
30 K2=K2/h^2;
31 % Definimos F
32 F=zeros(N+1,1);
33 F2=zeros(N,1);
34 % Calculamos u0
35 u0=[2*ones(1,2+N*4/8),4*ones(1,N*4/8)]';
36 % Guardamos la solucion para t=0 en la matriz solucion
37 sol(1,:)=u0'; 
38 U=u0; % Inicializacion
39 
40 % Calculamos U_n+1 a partir de U_n con un bucle for
41 for j=1:round(((length(t)-1)/T))
42 U=(eye(N+1)+dt*K)\(U+dt*F); % Euler Implicito
43 sol(j+1,:)= U'; % Guardamos la solucion
44 end
45 U(1)=[];
46 for j=round(((length(t)-1)/T)):length(t)-1 
47 U=(eye(N)+dt*K2)\(U+dt*F2); % Euler Implicito
48 sol(j+1,:)= [0,U']; % Guardamos la solucion
49 end
50 % Dibujamos la solucion
51 [xx,tt]=meshgrid(x,t);
52 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:

 1 % u_t-qu_xx=0, x en (0,L)
 2 % u_x(0,t)=0 Condicion Neumann para t<=1
 3 % u(0,t)=0 Condicion Dirichlet para t>1
 4 % u_x(L,t)=0 Condicion Neumann
 5 % u(x,0)=u0(x)
 6 
 7 % Datos
 8 L=8;
 9 T=100;
10 q=1;
11 
12 % Discretizacion espacial
13 h=0.15;
14 N=round(L/h);
15 x=0:h:L; 
16 xint=h:h:L; % Vector de nodos interiores
17 
18 % Discretizacion temporal
19 dt=h/4; % Tamaño de paso
20 t=0:dt:T; % Vector de tiempos
21 
22 % Matriz K para t<=1
23 K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
24 K(1,2)=-2; % Condicion Neumann
25 K(N+1,N)=-2; % Condicion Neumann
26 K=q*K/h^2;
27 % Matriz K para t>1
28 K2=2*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
29 K2(N,N-1)=-2;
30 K2=K2/h^2;
31 % Definimos F
32 F=zeros(N+1,1);
33 F2=zeros(N,1);
34 % Calculamos u0
35 u0=[2*ones(1,2+N*4/8),4*ones(1,N*4/8)]';
36 % Guardamos la solucion para t=0 en la matriz solucion
37 sol(1,:)=u0'; 
38 U=u0; % Inicializacion
39 
40 % Calculamos U_n+1 a partir de U_n con un bucle for
41 for j=1:round(((length(t)-1)/T))
42 U=(eye(N+1)+dt*K)\(U+dt*F); % Euler Implicito
43 sol(j+1,:)= U'; % Guardamos la solucion
44 end
45 U(1)=[];
46 for j=round(((length(t)-1)/T)):length(t)-1 
47 U=(eye(N)+dt*K2)\(U+dt*F2); % Euler Implicito
48 sol(j+1,:)= [0,U']; % Guardamos la solucion
49 end
50 % Dibujamos la solucion
51 [xx,tt]=meshgrid(x,t);
52 surf(xx,tt,sol);
53  
54  % Calculamos cual es el error del 5%
55 limite=0.05*(23.550/8);
56 % Buscamos para que tiempo se alcanza ese valor
57 for i=1:length(t)
58 dist=(sqrt(sum(sol(i,:).^2)/(N+1)));
59 if dist<limite
60 tiempo=(i-1)*dt
61 break
62 end
63 end