Tratamiento de residuos mediante análisis compartimental (Grupo 15B)

De MateWiki
Revisión del 12:25 4 mar 2014 de TrabajoEDO15B (Discusión | contribuciones) (Resolución numérica)

Saltar a: navegación, buscar
Warning.png Este artículo está en versión beta. El autor de este artículo no lo ha terminado todavía, por favor no lo edites hasta que elimine este mensaje.

1 Introducción

Tenemos tres pantanos A, B y C con 300 Hm3 de agua cada uno y con capacidad máxima de 500Hm3. Los pantanos A y C pueden recibir agua procedente de un rio y también pueden desalojarla por canales artificiales, pero tanto la entrada como salida de agua está regulada. Al pantano B puede fluir agua del rio, pero este pantano no tiene canal de desalojo de agua. Se produce un vertido tóxico al rio y al activar la entrada y salida a los pantanos del agua procedente del río, al pantano B empieza a entrar agua contaminada a razón de 2Hm3/día con 10kg/Hm de contaminante, mientras que:

hacia el pantano A entra agua contaminada a razón de 2Hm3/día con 5kg/Hm3 de contaminante y el agua contaminada sale por el canal a razón de 3Hm3/día; hacia el pantano C entra agua contaminada a razón de 2Hm3/día con 0,05kg/Hm3 de contaminante y el agua contaminada sale por el canal a razón de 1Hm3=día;

Cuando el pantano B se ha llenado, las autoridades cortan el acceso de agua al pantano B y filtran la entrada de agua a los pantanos A y C, de tal manera que deja de entrar agua contaminada.

En este momento, si designamos por xA(0), xB(0), xC(0) y VA, VB y VC la cantidad de contaminante y el volumen de agua respectivamente que hay en los pantanos A, B y C, determine estas cantidades.

2 Condiciones iniciales

Buscamos definir el aumento de la masa del contaminante respecto al tiempo en cada uno de los pantanos, incomunicados entre si. Para ello, calcularemos el desarrollo de la concentración, definida como la masa del soluto entre el volumen de la disolución.

Situación inicial para el cálculo analítico

La cual aplicada al desarrollo del volumen respecto al tiempo nos permite obtener las ecuaciones diferenciales correspondientes.

2.1 Pantano A:

Dados los datos iniciales del pantano:

[math]X_{entrada}= 5 Kg/Hm^3[/math]
[math]V_{entrada}= 2 Hm^3/día[/math]
[math]V_{salida}= 3 Hm^3/día[/math]
[math]V_{entrada soluto} (t)= 2 Hm^3/día \times 5 Kg/Hm^3=10 Kg/día[/math]
[math]V_{inicial}=300 Hm^3[/math]

Podemos desarrollar las ecuaciones que definen las variables respecto al tiempo:


[math]V(t)=V(0)+(V_{entrada}-V_{salida})t=300+(2-3)t[/math]
[math]C(t)=X(t)/(300+(2-3)t)[/math]
[math]V_{salida soluto} (t)=3 Hm^3/día \times Q(t)/(300-t)[/math]

Y así, el incremento diferencial de masa respecto al tiempo será:

[math]X_A’=10-3X_A(t)/(300-t)[/math]

Esta ecuación diferencial lineal se resuelve analíticamente, obteniendo una solución tal que, con la condición inicial de [math]Q(0)=0[/math] podemos reducir el haz a una solución final:

[math]X_A(t)= 5(300-t)-5,55 \times 10^{-5} \times(300-t)[/math]

2.2 Pantano B

En este caso el problema que se plantea es más sencillo, ya que este pantano no tiene salida antes de alcanzar las condiciones iniciales. Así pues:

[math]X_{entrada} = 2 Hm^3/día[/math]
[math]V_{max}= 500Hm^3[/math]
[math]V_{inicial}= 300 Hm^3[/math]
[math]V_{salida}= 0 Hm^3[/math]
[math]X_{entrada}=10 Kg/Hm^3[/math]

Al no vaciarse, tendrá un incremento directo de volumen:

[math] \Delta V= 500-300= 200 Hm^3[/math]

Con lo que calculamos el tiempo que tarda en llenarse:

[math]t= 200 Hm^3/ 2 Hm^3/día = 100 días[/math]

2.3 Pantano C

Este caso es idéntico al Pantano A:

[math]V_{entrada}= 2 Hm^3/día[/math]
[math]V_{max}= 500 Hm^3[/math]
[math]V_{inicial}= 300 Hm^3[/math]
[math]X_{entrada}= 0,05 Kg/Hm^3[/math]
[math]V_{salida}=1 Hm^3/día[/math]
[math]V_{entrada soluto}= 0,05 Kg/Hm^3 \times 2 Hm^3/día= 0,1 Kg/día[/math]
[math]V(t)=V(0)+(V_{entrada}-V_{salida})t=300+(2-3)t[/math]
[math]C(t)=X(t)/(300+(2-1)t)[/math]
[math]V_{salida}= 300 Hm^3/díax X_C/(300+t)[/math]

La ecuación será:

[math]X_C'= 0,1-3X_C/(300+t)[/math]

Siendo su solución en [math]Q(0)=0[/math] la siguiente:

[math]X_C(t)= (300+t)/40-300^4/40(300+t)^3[/math]

2.4 Obtención de los Valores Iniciales

Obtenidas las ecuaciones que definen los 3 pantanos, basta con aplicar la condición de t=100 días para obtener las condiciones iniciales de volumen de disolución y masa del soluto en cada uno de ellos:

Pantano A:

[math]X_A=556 Kg[/math]
[math]V_A=200 Hm^3[/math]

Pantano B:

[math]X_B=2000 Kg[/math]
[math]V_B=500 Hm^3[/math]

Pantano C:

[math]X_C=6, 8359 Kg[/math]
[math]V_C=400 Hm^3[/math]

3 Resolución numérica

3.1 Planteamiento

Para intentar eliminar el contaminante de los pantanos, se bombea agua entre estos, quedando la situación como sigue: al pantano A sigue entrando agua desde el rio a razón de 2Hm3/día, pero sin contaminar, y sale por el canal artificial 3Hm3/día de agua contaminada. Del pantano B al A entra agua a razón de 3Hm3/día y de A hacia B a razón de 2Hm3/día; al pantano C sigue entrando agua desde el rio a razón de2Hm3/día, pero sin contaminar y sale por el canal artificial 1Hm3/día de agua contaminada. Del pantano B recibe 1Hm3/día y sale hacia B 2Hm3/día;

Resolvemos numéricamente este problema de valor inicial con N: 50, 100, 300 (medimos en días) utilizando el método de Euler, método trapezoidal y método de Runge-Kuta, con una longitud de paso h = 1 y h = 0:1.

Con la nueva dinámica que producen los trasbases, pasamos a trabajar con los siguientes datos, además de los datos iniciales anteriormente obtenidos:

Situación inicial para el cálculo numérico

A partir de ahora los pantanos no trabajan de forma independiente, sino que el desarrollo de cada uno es dependiente de los otros dos. Así pues, trabajaremos con el siguiente sistema de ecuaciones diferenciales:

Sistema a resolver por los distintos métodos numéricos

3.2 Método de Euler

%Ejercicio de mezclas por el método de Euler
t0=0;
tf=50;
y0=[556 2000 6.8359]';
A=[-1/40 3/500 0; 1/100 -1/125 1/200; 0 1/500 -3/400]; 
%matriz de coeficientes
%Datos de la discretizacion
h=0.1;
N=(tf-t0)/h; 
%Construir de los vectores de tiempos y solucion aproximada
t=t0:h:tf;
y=zeros(3,N+1); 
%Hacemos una matriz con todas las soluciones conjuntas
%Inicialización
y(:,1)=y0; 
%la primera columna de la matriz de soluciones vale el valor de y0
yy=y0;
%Iteraciones y(tn)>>y(tn+1)
for n=1:N
  yy=yy+h*(A*yy);  
  y(:,n+1)=yy;
 end
clf
plot(t,y)
xlabel('Tiempo(días)')
ylabel('Cantidad soluto (kg)')
title('Metodo Euler h=0,1 N=50')
legend('Pantano A','Pantano B', 'Pantano C')


3.3 Método Trapezoidal

%Ejercicio de mezclas por el método del Trapecio
%Datos del problema
t0=0;
tf=50;
y0=[556 2000 6.8359]';
A=[-1/40 3/500 0; 1/100 -1/125 1/200; 0 1/500 -3/400];
%Matriz de coeficientes
%Datos de la discretizacion
h=0.1;
N=(tf-t0)/h; 
%Construir de los vectores de tiempos y solucion aproximada
t=t0:h:tf;
y=zeros(3,N+1); 
%Hacemos una matriz con todas las soluciones conjuntas
%Inicializacion
y(:,1)=y0; 
%La primera columna de la matriz de soluciones vale el valor de y0
yy=y0;
%Iteraciones y(tn)>>y(tn+1)
for n=1:N
  yy=(eye(3)-h/2*A)\((eye(3)+h/2*A)*yy);
  y(:,n+1)=yy;
 end
clf
plot(t,y)
xlabel('Tiempo(días)')
ylabel('Cantidad soluto (kg)')
title('Metodo Trapecio con h=0.1 y N=50')
legend('Pantano A','Pantano B','Pantano C')
legend('Pantano A','Pantano B','Pantano C')

3.4 Método de Runge-Kuta

%Ejercicio de mezclas por el método Runge Kutta 4
%Datos del problema
t0=0;
tf=50;
y0=[556 2000 6.8359]';
A=[-1/40 3/500 0; 1/100 -1/125 1/200; 0 1/500 -3/400];
%Datos de la discretizacion
h=0.1;
N=tf;
%Vectores de tiempos y aproximaciones
t=t0:h:tf;
y=zeros(3,N+1);
%metemos una matriz de ceros que iremos rellenando con el bucle 
%Inicializacion
y(:,1)=y0; 
%Para guardar la aproxmacion en cada tiempo vamos a usar yy es un valor numerico q va cambiando
yy=y0; 
%Ahora ya hacemos el bucle que vaya calculando el valor siguiente en cada instante las aproximaciones nos permite pasar del y(tn) al y(tn+1).
%Tendremos que hacer tantas iteraciones como intervalos, para eso tenemos N
for n=1:N
  %ahora depende del metodo numerico, antes de hacer la media ponderada tenemos q hacer la constante
  %Calculo k1=f(tn,yn)
    k1=[-1/40*yy(1)+3/500*yy(2),0.01*yy(1)-1/125*yy(2)+1/200*yy(3),1/500*yy(2)-3/400*yy(3)]';
  %Hemos metido el f(tn,yn) como una matriz en la q cada columna es una %de las dos ecuciacones del sistema y1 e y2 para k2 calculamos el argumento en el 
  %tiempo y en el valor de y
  yp=yy+0.5*h*k1;
  k2=[-1/40*yp(1)+3/500*yp(2),0.01*yp(1)-1/125*yp(2)+1/200*yp(3),1/500*yp(2)-3/400*yp(3)]';  
  %para k3 el argumento en el tiempo es el mismo q en k2 pero cambia el argumento en y
  yp=yy+1/2*h*k2;
  k3=[-1/40*yp(1)+3/500*yp(2),0.01*yp(1)-1/125*yp(2)+1/200*yp(3),1/500*yp(2)-3/400*yp(3)]';  %calculo de k3
   %para k4 cambiamos ambos argumentos
  yp=yy+h*k3;
  k4=[-1/40*yp(1)+3/500*yp(2),0.01*yp(1)-1/125*yp(2)+1/200*yp(3),1/500*yp(2)-3/400*yp(3)]';
    yy=yy+h/6*(k1+2*k2+2*k3+k4);
  y(:,n+1)=yy;
end

 
%Dibujo la solucion aproximada
figure(1)
clf
plot(t,y)
xlabel('Tiempo(días)')
ylabel('Cantidad soluto (kg)')
title('Metodo Runge Kutta 4 con h=0.1 y N=50')
legend('Pantano A','Pantano B','Pantano C')