function tm7
%TM7 resuelve el problema 7 de transporte de materia.
%   En un reactor Berty, que puede suponerse de mezcla completa, de 300 cm3
%   de volumen, se transforma anisol (metil-fenil-éter) en p-cresol a 723K 
%   y 101330 N/m2 de presión sobre un catalizador de fosfato de aluminio.
%   Suponiendo que la transformación se desarrolla según una reacción de 
%   primer orden, calcular el perfil de concentraciones de anisol dentro 
%   del catalizador, la velocidad de desaparición del anisol y el factor de
%   eficacia del catalizador.
%   
%   Datos:
%
%      - Catalizador: partículas cilíndricas de 1.6 mm de radio y 3 mm de
%        altura, con una porosidad de 0.2
%      - Porosidad del lecho: 0.6
%      - Difusividad efectiva del anisol en las partículas:DA=3.5e(-9) m2/s
%      - Constante de la velocidad de reacción: k1 = 0.003 (1/s)
%      - Concentración del anisol en la superficie externa del catalizador:
%        CA1 = 7.2e(-2) (kg/m3)
%
%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=concentración (en base másica) y x=radio.

VR = 3e-4; % volumen del reactor
T = 723;
P = 101330;
R1 = 1.6e-3;
R0 = 1e-12;
% no es posible usar R0=0 porque da error de valor singular (es 
% consecuencia de que en la ecuación [79] aparece 2/r, por lo que si usamos
% 0 dará indeterminación de infinito); por tanto para evitar este problema
% utilizamos un número suficientemente pequeño y próximo a cero.
L = 3e-3;
alfa_part = 0.2;
alfa_lecho = 0.6;
DA = 3.5e-9;
k1 = 3e-3;
CA1 = 7.2e-2;
% Definimos los datos del problema

global re;
re = 3*pi*R1^2*L/(2*pi*R1*L+2*pi*R1^2);
% Calculamos el radio equivalente de las partículas de catalizador
% cilíndricas, re(m).

function dydx = odefun_tm7(x,y)
%ODEFUN_TM7  Define el sistema de ODEs del problema TM7. 
%   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 TM7 quedaría definido con la siguiente 
%   expresión:
dydx =  [ y(2)
          k1/DA*y(1)-2/x*y(2) ];
end
     
 function res = bcfun_tm7(ya,yb)
%BCFUN_TM7  Define las condiciones límite del problema TM7.
%   RES = BCFUN_TM7 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=R0, b=re). 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). 
res = [ ya(2)
        yb(1)-CA1];
end
    
function yinit = yinit_tm7(x)
%YINIT_TM7 Solución aproximada (guess) para el problema TM7.
%   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=R0), o bien la 
%   expresión de una función que aproxime el comportamiento de la solución 
%   buscada; en el problema TM7 hemos optado por esta última opción 
%   sabiendo que el perfil de concentraciones del anisol es parabólico y 
%   presenta su mínimo en el centro de la partícula, por lo que una 
%   ecuación polinómica de segundo orden con su mínimo en x=R0=0, es decir 
%   "x^2", 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^2
          2*x   ];
end

solinit_tm7 = bvpinit(linspace(R0,re,10),@yinit_tm7);
%BVPINIT Es una función auxiliar para crear la estructura inicial de la
%solución aproximada.
%   Una vez tenemos YINIT_TM7 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_TM7. 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 R0 y re.
%           -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_TM7. 

sol = bvp4c(@odefun_tm7,@bcfun_tm7,solinit_tm7);
%BVP4C es la herramienta de resolucion de Matlab 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_TM7.
%           - bcfun: evaluación de los residuos en las condiciones límite;
%                en nuestro caso BCFUN_TM7.
%           - solinit: estructa de la solución inicial aproximada; en
%                nuestro caso SOLINIT_TM7.
%
%   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(R0,re,101);
[Sxint, SPxint] = deval(sol,xint);
%PERFIL CONCENTRACIONES: Para construir el perfil de concentraciones 
%   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 R0 y re en subintervalos de tamaño
%   0.000001565 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).

NAr_part = -DA*SPxint(1,:);
%PERFIL DENSIDAD DE FLUJO MASICO POR TRANSPORTE TOTAL PARA UNA PARTÍCULA,
%   NAr_part: 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 NAr_part para nuestro caso particular.

global FAr_part_re n_part FAr_lecho_re_gh eficacia;
FAr_part = 4*pi*re^2*alfa_part*NAr_part;
%FAr_part: velocidad de desaparición del anisol para una partícula
FAr_part_re = FAr_part(1,101);
%FAr_part_re: FAr_part en la superficie externa de la partícula (r=re)
n_part = (1-alfa_lecho)*VR/(pi*R1^2*L);
%n_part: número total de partículas en el lecho
FAr_lecho = FAr_part*n_part;
%FAr_lecho: velocidad de desaparición del anisol en todo el lecho
FAr_lecho_re_gh = FAr_lecho(1,101)*3600*1000;
%FAr_lecho_re: FAr_lecho en la superficie externa de todo el lecho en g/h
FAT = -4/3*pi*re^3*alfa_part*k1*CA1;
%FAT: velocidad de desaparción de anisol si la resistencia de la difusión es
%despreciable.
eficacia = FAr_part_re/FAT;

plot(xint,Sxint(1,:));
grid
title('Perfil Concentraciones')
ylabel('CA (kg/m3)')
xlabel('r (m)')
figure(2)
plot(xint,FAr_lecho)
grid
title('Velocidad de desaparición del anisol')
ylabel('FAr-lecho (kg/s)')
xlabel('r (m)')
global RESULTADOS7;
RESULTADOS7=[ xint' Sxint(1,:)' FAr_lecho'];
RESULTADOS7_VAR=[ '      r(m)      ' 'CA(kg/m3) ' 'FAr_lecho(kg/s)  '];
format short e
disp(RESULTADOS7_VAR)
disp(RESULTADOS7)
re
FAr_part_re 
n_part
FAr_lecho_re_gh
eficacia
%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 RESULTADOS7. Además finalmente mostramos por pantalla
%los valores de re(m), FAr_partR1(mol/s*part), n_part, FAr_lechoR1(mol/s),
%eficacia.
end
%
%EJECUCIÓN DE TM7: Ejecutar en la ventana de comandos:
%
%   global RESULTADOS7 re FAr_part_re n_part FAr_lecho_re_gh eficacia; tm7
%
% 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 RESULTADOS7. así como los valores de otros cálculos también 
% almacenados como variables globales, y que ahora tenemos disponibles 
% 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
