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

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.

3 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))

"La masa es constante en el tiempo"

"Variación de la masa en el punto medio"

6 Estudio de la concentración.

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

% 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 Uso de un limpiador de contaminante. 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 que este valor es 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
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 .

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