Modelos de mezclas (grupo 5)

1 Introducción

Tenemos dos pantanos A y B, con 100 [math]Hm^3[/math] de agua cada uno, que están unidos por una presa que deja pasar agua de A a B. El pantano A recibe 3 [math]Hm^3/día[/math] de agua limpia proveniente de ríos y el B 1.5 [math]Hm^3/día[/math]. Para mantener el nivel de los pantanos estable la presa de A a B deja pasar una media de 3 [math]Hm^3/día[/math] mientras que la presa al final de B desaloja 4.5 [math]Hm^3/día[/math]. Se produce un vertido tóxico en el pantano A que deja 20 toneladas de un cierto contaminante. Supongamos que se dan las siguientes hipótesis::

• El contaminante está disuelto de forma homogénea en el agua de los pantanos.
• Al entrar o salir agua en un pantano, esta se mezcla con el agua del pantano de forma inmediata creando una mezcla homogénea.
• La variación de contaminante en un lago es la diferencia entre el contaminante que entra y sale en el lago, es decir si denotamos por [math] x_A (t) [/math] la cantidad de contaminante en el lago A, se tiene:: [math] \frac{dx_A}{dt} = velocidad de entrada – velocidad de salida[/math]

2 Planteamiento

2.1 Sin plan de limpieza

Sean [math] x_A [/math] y [math] x_B [/math] las cantidades de contaminante en cada uno de los lagos, la cantidad de contaminante en un instante de tiempo determinado se obtiene como diferencia entre la velocidad de entrada y salida de agua contaminada. En el lago A sólo entra agua limpia mientras que la velocidad de salida se obtiene como el producto de la concentración de contaminante en el lago A [math] \frac{X_A}{100} [/math] por el caudal de agua que sale 3 [math]^3/día [/math]. Así mismo, en el lago B la velocidad de entrada de agua contaminada es la misma que la velocidad de salida de agua contaminada del lago A, ya que el lago A vierte el agua contaminada en B, y la velocidad de salida del agua de B es el producto de su concentración [math] \frac{X_B}{100} [/math] por el caudal de agua que sale 4.5 [math] Hm^3/día [/math].

Para modelizar el problema, se plantea el siguiente sistema dos ecuaciones diferenciales de primer orden:

[math] \left\{\begin{matrix}\frac{\mathrm{d}X_A}{\mathrm{d} t}=-\frac{3X_A}{100}\\\frac{\mathrm{d} X_B}{\mathrm{d} t}=\frac{3X_A}{100}-\frac{4.5X_B}{100}\\X_A(t_{0})=20, \quad X_B(t_{0})=0\end{matrix}\right. [/math]

2.2 Con plan de limpieza

A continuación se propone activar un plan de limpieza basado en bombear 1[math] Hm^3/día [/math] de agua del pantano B al A, por lo que para mantener el volumen de agua entre los lagos constante se tiene que modificar el flujo entre ambos.

De este modo, si al lago A le llegan 3[math]Hm^3/día[/math] provenientes de los ríos y 1[math]Hm^3/día[/math] proveniente de B, el lago A deberá verter a B 4[math]Hm^3/día [/math] en lugar de los '3 [math] Hm^3/día [/math] que estaba vertiendo antes de activar el plan de limpieza. El lago B recibe 1.5 [math] Hm^3/día [/math] provenientes de los ríos y, como se ha visto, 4 [math] Hm^3/día [/math] provenientes del lago A. Tras activar el plan de limpieza bombea 1 [math] Hm^3/día [/math] a A, por lo que el vertido que desaloja B aguas abajo se mantendrá constante igual a 4.5 [math] Hm^3/día [/math] para mantener el volumen de agua. Así se plantea un sistema de ecuaciones diferenciales similar al anterior:: [math] \left\{\begin{matrix}\frac{\mathrm{d}X_A}{\mathrm{d}t}=-\frac{4X_A}{100}+\frac{1X_B}{100}\\\frac{\mathrm{d}X_B}{\mathrm{d} t}=\frac{4X_A}{100}-\frac{5.5X_B}{100}\\X_A(t_{0})=20, \quad X_B(t_{0})=0\end{matrix}\right. [/math]

2.3 Unión de un tercer pantano

Si hubiese un tercer pantano, al que se le llamará C, unido al anterior que desalojase 6 [math]Hm^3/día [/math] y recibiese 1.5 [math]Hm^3/día [/math], la velocidad de entrada de agua sería igual a la de salida de B, es decir, 4.5 [math] Hm^3/día [/math] ,mientras que la de salida serían estos 4.5 [math]Hm^3/día [/math] más 1.5 [math] Hm^3/día [/math] provenientes de los ríos, dando lugar a los 6 [math]Hm^3/día [/math] para que se mantenga constante el volumen. Así, el sistema de ecuaciones diferenciales quedará de la siguiente forma::

[math] \left\{\begin{matrix}\frac{\mathrm{d}X_A}{\mathrm{d}t}=\frac{-3X_A}{100}\\\frac{\mathrm{d}X_B}{\mathrm{d} t}=\frac{3X_A}{100}-\frac{4.5X_B}{100}\\\frac{\mathrm{d}X_C}{\mathrm{d} t}=\frac{4.5X_B}{100}-\frac{6X_C}{100}\\X_A(t_{0})=20, \quad X_B(t_{0})=0\end{matrix}\right. [/math]

3 Resolución numérica

A partir de este momento nos olvidamos del pantano C. Para resolver estos sistemas se utilizará un programade análisis numérico (en este caso Matlab®).

3.1 Método de Euler

Se utilizan los siguientes códigos::

Sin plan de limpieza:

clear all

N=100000; %número de iteraciones
t0=0;tN=325; %tiempo inicial y final
h=(tN-t0)/N; %paso: espacio entre dos iteraciones
x0=[20 0]'; %condiciones iniciales
A=[-0.03 0;0.03 -0.045];
x=x0;
xA(1)=x(1);
xB(1)=x(2);

for n=1:N    
  x=x+h*A*x;    
  xA(n+1)=x(1);    
  xB(n+1)=x(2);
end
y=t0:h:tN;

%dibujamos la gráfica
hold on
plot(y,xA,'b') %contaminante en pantano A: azul
plot(y,xB,'r') %contaminante en pantano B: rojo
title('Sin plan de limpieza')
xlabel('Tiempo (días)')
ylabel('Contaminante (T)')
hold off

Con plan de limpieza:

clear all

N=100000; %número de iteraciones
t0=0;tN=325; %tiempo inicial y final
h=(tN-t0)/N; %paso: espacio entre dos iteraciones
x0=[20 0]'; %condiciones iniciales
B=[-0.04 0.01;0.04 -0.055];
X=x0;
Xa(1)=X(1);
Xb(1)=X(2);

for n=1:N    
  X=X+h*B*X;    
  Xa(n+1)=X(1);    
  Xb(n+1)=X(2);
end
y=t0:h:tN;

%dibujamos la gráfica
hold on
plot(y,Xa,'g') %contaminante en pantano A: verde
plot(y,Xb,'m') %contaminante en pantano B: magenta
title('Con plan de limpieza')
xlabel('Tiempo (días)')
ylabel('Contaminante (T)')
hold off

Si queremos ver la diferencia de tiempo que tarda en desaparecer la mitad de contaminante en el pantano A activando el plan de limpieza o sin activarlo no tenemos más que observar las gráficas.

....

3.2 Método de Runge-Kutta

Aplicando el método deRunge-Kutta de cuarto orden se utilizan el siguiente código:

clear all

N=10000; %número de interacciones
t0=0; tN=325;  %tiempo inicial y final
x0=[20 0]'; %condiciones iniciales
h=(tN-t0)/N; %paso: espacio entre dos iteraciones


%sin plan de limpieza
A=[-0.03 0;0.03 -0.045];
x=x0;
xA(1)=x(1);
xB(1)=x(2);
for n=1:N     
   k1=A*x;    
   k2=A*(x+(1/2*k1*h));    
   k3=A*(x+(1/2*k2*h));    
   k4=A*(x+(k3*h));    
   x=x+h/6*(k1+2*k2+2*k3+k4);    
   xA(n+1)=x(1);    
   xB(n+1)=x(2);
end


%Con plan de limpieza
B=[-0.04 0.01;0.04 -0.055];
X=x0;
Xa(1)=X(1);
Xb(1)=X(2);
for n=1:N     
   k1=B*X;    
   k2=B*(X+(1/2*k1*h));    
   k3=B*(X+(1/2*k2*h));    
   k4=B*(X+(k3*h));    
   X=X+h/6*(k1+2*k2+2*k3+k4);    
   Xa(n+1)=X(1);    
   Xb(n+1)=X(2);
end


y=t0:h:tN;


%dibujamos las dos gráficas
hold on
figure(1) %sin plan de limpieza
plot(y,xA,'g'); %contaminante pantano A: verde
plot(y,xB,'m'); %contaminante pantano B: magenta
title('Sin plan de limpieza')
xlabel('Tiempo (dias)')
ylabel('Contaminante (T)')
figure(2) %con plan de limpieza
plot(y,Xa,'b') %contaminante pantano A: azul
plot(y,Xb,'r') %contaminante pantano B: rojo
title('Con plan de limpieza')
xlabel('Tiempo (dias)')
ylabel('Contaminante (T)')
hold off

Ambas soluciones quedan reflejadas en las siguientes gráficas:

...

3.3 Comparación método Euler con método Runge-Kutta

Para comparar ambos métodos basta con comparar sus gráficas. Realizamos el procedimiento con ambos sistemas. La gráfica del primer sistema es la siguiente: ...

Como puede verse la diferencia entre ambos métodos en tan pequeña que las gráficas aparecen superpuestas. Para poder apreciar la diferencia entre ellas ampliamos la grafica cogiendo dos intervalos pequeños detiempo. Runge-Kutta aparece representado en azul y Euler en verde. ...

Realizamos el mismo proceso con el segundo sistema y como es de esperar sucede lo mismo: ...

Una vez comparados ambos métodos se puede afirmar que la diferencia entre ellos es casi inapreciable.

4 Contaminante inicial desconocido

...