function tm5
%TM5 resuelve el problema 5 de transporte de materia.
%
%   Una columna de paredes mojadas térmicamente aislada se utiliza para la 
%   rectificación adiabática de una mezcla de benceno y tolueno cuyo 
%   comportamiento puede suponerse ideal y con calores latentes de 
%   vaporización de ambos componentes prácticamente iguales. Una caldera 
%   situada en su base le suministra una corriente vapor y un condensador 
%   situado en su parte superior devuelve a la columna como reflujo, en 
%   forma de capa líquida periférica, parte del vapor condensado. En una 
%   cierta sección transversal de la columna, la temperatura es de 363 K y 
%   las corrientes de vapor y líquido contienen 70 y 59 moles por 100 de 
%   benceno respectivamente. Suponiendo que la resistencia total a la 
%   transferencia de ambos componentes que ofrece el vapor equivale a la 
%   que ofrecería a su difusión una película estacionaria del mismo sobre 
%   la superficie líquida, de 0.2 mm de espesor, calcular los perfiles de 
%   fracciones molares en la película estacionaria y los caudales molares 
%   de intercambio de los dos componentes entre ambas fases por unidad de
%   superficie interfacial.
%   Datos:
%
%      - Presión en toda la columna: 101330 N/m2
%      - Presion de vapor del tolueno a 363 K: 53945 N/m2
%      - Difusividad mútua de ambos componentes: 5.1e(-6) m2/s
%
%Aplicaremos la herramienta de resolución (solver) BVP4C para resolver
%problemas de ecuaciones diferenciales ordinarias (Ordinary Differential
%Equations, ODEs) con condiciones de valores límite (Boundary Value 
%Problem, BVP). 

% NOTA: A nivel interno de cálculo, y dado que la estructura "solinit" que
% Matlab utiliza en BVP4C maneja como variables "x,y", tomaremos 
% y=fracción molar y x=z(distancia sobre el eje z).

yB1 = 0.7;
xB2 = 0.59;
DBT = 5.1e-6;
R = 8.314;
P = 101330;
PvT=53945;
T = 363;
z1= 0;
z2 = 0.0002;
% Definimos los datos del problema

yT1 = 1-yB1;
PT2 = PvT*(1-xB2);
yT2 = PT2/P;
% Calculamos las fracciones molares del tolueno en los límites de la
% película estacionaria

function dydx = odefun_tm5(x,y)
%ODEFUN_TM5  Define el sistema de ODEs del problema TM5. 
%   Matlab sólo trabaja con ecuaciones diferenciales de primer orden, por
%   lo que si aparecen ecuaciones diferenciales de orden superior es
%   preciso transformarlas en un sistema de ecuaciones difereciales de
%   primer orden. Para ello se realizará un cambio de variables llamando
%   y(1)=y, y(2)=y', y(3)=y'', y(4)=y''', y así sucesivamente; en nuestro 
%   caso como tenemos una ecuación diferencial de orden dos es preciso 
%   definir sólo y(1)=y, y(2)=y'. El sistema de ecuaciones diferenciales de
%   primer orden del problema TM5 quedaría definido con la siguiente 
%   expresión:
dydx =  [ y(2)
           0 ];
end
     
 function res = bcfun_tm5(ya,yb)
%BCFUN_TM5  Define las condiciones límite del problema TM5.
%   RES = BCFUN_TM5 devuelve el vector columna RES con los residuos
%   resultantes de las soluciones aproximadas YA e YB en los extremos del
%   intervalo definido por las condiciones límite [a b] (a=z1, b=z2). El
%   problema BVP está resuelto cuando RES=0 (es decir cuando se anula la
%   diferencia entre la solución aproximada y el valor de la solución en 
%   los puntos correspondientes a las condiciones límite). 
%   
%   Nota: 
%   z1=0 corresponde al comienzo de la película estacionaria y z2=0.0002 
%   corresponde al fin de dicha capa cuando se alcanza la superficie del 
%   líquido.

res = [ ya(1)-yT1
        yb(1)-yT2];
end
    
function yinit = yinit_tm5(x)
%YINIT_TM5 Solución aproximada (guess) para el problema TM5.
%   Para resolver problemas BVP es preciso introducir una aproximación a la
%   solución mediante un vector columna YINIT. La aproximación será bien
%   una aproximación al valor de la función y al valor de la primera 
%   derivada de la función en la condición límite xa (a=z1), o bien la 
%   expresión de una función que aproxime el comportamiento de la solución 
%   buscada; en el problema TM5 hemos optado por esta última opción 
%   sabiendo que el perfil de fracciones molares del tolueno es lineal 
%   descendente por lo que la ecuación polinómica de primer grado "-x" 
%   será adecuada. Tomando este criterio el primer elemento del vector 
%   columna YINIT será esta función y el segundo elemento su primera 
%   derivada:
yinit = [ -x
          -1   ];
end

solinit_tm5 = bvpinit(linspace(z1,z2,10),@yinit_tm5);
%BVPINIT Es una función auxiliar para crear la estructura inicial de la
%solución aproximada.
%   Una vez tenemos YINIT_TM5 Matlab precisa crear una matriz inicial de
%   valores aproximados de la solución para diferentes valores de x
%   dentro del intervalo definido por las condiciones límite. Esa matriz se
%   construye utilizando la función auxiliar BVPINIT, cuya solución
%   almacenaremos en SOLINIT_TM5. La sintaxis de BVPINIT es la siguiente:
%       
%       bvpinit(x,yinit)
%
%           - x: valores de x para la matriz inicial. Aquí es donde se
%                definen que valores de x comprendidos entre las 
%                condiciones límite [a b] ensayaremos para calcular la 
%                solución inicial aproximada. En nuestro caso nos hemos 
%                basado en el comando LINSPACE para crear 10 puntos 
%                linealmente espaciados entre z1 y z2.
%           -yinit: aproximación del valor de la solución en los valores de 
%                x de la matriz inicial. En nuestro caso la hemos definido 
%                en la función previa YINIT_TM5. 

sol = bvp4c(@odefun_tm5,@bcfun_tm5,solinit_tm5);
%BVP4C es la herramienta de resolucion de Matlaba para BVPs
%   La sintaxis de BVP4C es la siguiente:
%
%       bvp4c(odefun,bcfun,solinit)
%
%           - odefun: sistema de ecuaciones diferenciales de primer orden
%                a solucionar; en nuestro caso ODEFUN_TM5.
%           - bcfun: evaluación de los residuos en las condiciones límite;
%                en nuestro caso BCFUN_TM5.
%           - solinit: estructa de la solución inicial aproximada; en
%                nuestro caso SOLINIT_TM5.
%
%   La solución generada se almacena en la estructura SOL, donde:
%           - sol.x: son los valores de x seleccionados por BVP4C para
%                evaluar la solución.
%           - sol.y: la solución del problema en los puntos sol.x.
%           - sol.yp: evaluación de la primera derivada de la solución, y',
%                en los puntos sol.x.
%           - sol.solver: nombre de la herramienta utilizada para resolver
%                el problema, en nuestro caso 'bvp4c'.
%
xint = linspace(z1,z2,101);
[Sxint, SPxint] = deval(sol,xint);
NTz = -DBT*P/(R*T)*SPxint(1,:);
yB = 1-Sxint(1,:);
%PERFIL DE FRACCIONES MOLARES: Para construir el perfil de fracciones 
%   molares utilizaremos la función auxiliar DEVAL que nos devuelve el 
%   valor de la solución del problema para los valores de x que deseemos. 
%   En nuestro caso dividiremos el intervalo entre z1 y z2 en subintervalos
%   de tamaño 0.000002 m mediante el comando LINSPACE, guardando el valor 
%   de cada punto en XINT. Almacenaremos los valores de la solución en la 
%   variable SXINT (los valores se alamcenan en la primera fila de SXINT).
%PERFIL DENSIDAD DE FLUJO MOLAR POR TRANSPORTE TOTAL, NTz: La orden DEVAL
%   también nos permite obtener el valor de la primera derivada del
%   polinomio que interpolaría la solución entre los puntos de la
%   solución, resumiendo la primera derivada de la solución, la cual se
%   almacena en la variable SPXINT(los valores se alamcenan en la primera 
%   fila de SPXINT). A partir de SPXINT obtendremos el perfil 
%   de NTz para nuestro caso particular (es constante). 
plot(xint,Sxint(1,:),xint,yB,'r');
grid
title('Perfil Fracciones Molares')
ylabel('y (fracción molar)')
xlabel('z (m)')
legend('Tolueno','Benceno');
figure(2)
plot(xint,NTz)
grid
title('Perfil Densidad de Flujo Molar por Transporte Total')
ylabel('NTz=-NBz (mol/sm2)')
xlabel('z (m)')
global RESULTADOS5;
RESULTADOS5=[ xint' Sxint(1,:)' yB' NTz'];
RESULTADOS5_VAR=[ '      z(m)       ' 'yTolueno    ' 'yBenceno   ' 'NTz (mol/sm2)'];
format short e
disp(RESULTADOS5_VAR)
disp(RESULTADOS5)
%REPRESENTACIÓN Y OBTENCIÓN DE LOS RESULTADOS: La secuencia de comandos
%anterior permite la representación de gráficos de los perfiles, la
%obtención de los resultados en una tabla y el almacenamiento de estos en
%la variable global RESULTADOS5.
end
%
%EJECUCIÓN DE TM5: Ejecutar en la ventana de comandos:
%
%                global RESULTADOS5; tm5
%
% Obtenemos los resultados en una tabla en la ventana de comandos, los 
% gráficos de los perfiles y los resultados almacenados en la variable
% global RESULTADOS5 que ahora tenemos disponible en el espacio de trabajo 
% por si queremos hacer algún cálculo posterior con ellos.
%
% 000000000000000000000000000000000
% ---------------------------------
% Juan F. Martínez Gallegos
% Dpto. Ingeniería Química
% UGR
% ---------------------------------
% 000000000000000000000000000000000
