%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Difusion en una Pelicula liquida descendente: Transferencia de 
%especies por conveccion forzada
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%Ultimo a las 03:09 am del jueves 22-06-2023
%
%1 - Analisis del problema:
%SUPUESTOS:
%       --Consideramos concretamente la absorción de un gas A por una 
%pelicula laminar descendente de liquido B. 
%       --La substancia A es solo ligeramente soluble en B, de forma que la 
%viscosidad del liquido no varia apreciablemente.
%       --Haremos la restriccion que la difusion en la pelicula liquida
%       tiene lugar tan lentamente que A "penetra" muy poco en B (la
%       distancia de penetracion es pequeña en comparacion con el espesor
%       de la pelicula).
%
%Consideracion especial:
%
%Para el calculo del perfil de concentracion de la especie A, se piensa que
%la substancia A "tiene la impresion" de que la pelicula se mueve toda ella
%con una velocidad igual a la maxima. Por otro lado, si la penetracion de A
%no es muy grande, no se estaria "enterando" de la existencia de la pared
%solida situada en x=delta. Por otro lado, si la pelicula tuviese un
%espesor infinito y se moviese con la velocidad maxima, el material que
%difunde no notaria la diferencia.
%
%2 - Definicion del problema:
%
%Datos:
%       --Velocidad maxima del perfil laminar en la interfase en m/s: 0.2
%       --Difusividad de la especie A en la mezcla binaria AB en m^2/s:
%       5.0*10^(-6)
%       --Constante c=D_{AB}/V_{max}
%       --Concentracion de la especie quimica A en mol/L o kmol/m^3: 10
%
%Salidas:
%       --Perfil de concentración de la especie A a traves del medio de B
%       --x y z, es decir, las distancias horizontales y verticales 
%       respectivamente conforme avanza la penetracion de especie gaseosa 
%       A en el medio liquido B.
%
%3 - Esquema del proceso (algoritmo y programa)
%
%El solver de PDE pdepe de MATLAB resuelve sistemas de varias PDE 
%parabolicas y elipticas en 1D.
%El coeficiente f(x,t,u,dudx) es un termino de flujo y s(x,t,u,dudx) es un 
%termino de origen.
%
function transfmatconvforzada
function [c,f,s]=pdeconvforz(x,z,u,dudx)
    %Variables de entrada:
    c=200000;
    f=dudx;
    s=0;
end

%Condición al inicio del dominio de z:
    function u0=icpde(x)
        u0=0;
    end
%Condiciones de Frontera en el dominio de x:
    function [pL,qL,pR,qR]=bcpde(xL,uL,xR,uR,z)
        pL=uL-10;% Se asume que la concentracion en x=0 es de 10 kmol/m^3
        qL=0;% Aqui "q a la izquierda" es cero por similitud con la 
             % ecuacion general para el calculo de condiciones de frontera.
        pR=uR;% Aqui la concentracion es 0 kmol/m^3 para cuando x es igual
              % al espesor de la pelicula liquida descendente "delta".
        qR=0;%Aqui "q a la izquierda" es cero por similitud con la ecuacion 
             %general para el calculo de condiciones de frontera.
    end
%Dominio de espacio en x: el espesor de la pelicula descendente es de
%aproximadamente 1 cm, sin embargo en las inmediaciones de la interfase
%gas-líquido es donde se genera la mayor diferencia de concentración o 
%fuerza impulsora, por lo que la concentración disminuirá notablemente
%para compensar esa diferencia de especies en ambas fases, sumado a la poca 
%solubilidad del gas A en el seno del liquido B.
x=linspace(0,0.05,50);
%Dominio de espacio en z: aquí la distancia es más bien arbitraria, dado
%que la distancia en dirección z es en donde los esfuerzas de gravedad 
%están influenciando. Por lo tanto para favorecer la película descendente
%intuitivamente la distancia debe ser mucho mayor comparada con la
%distancia en x.
z=linspace(0,5,50);
%Por ser geometría rectangular:
m=0;
%La solucion sera entonces:
sol=pdepe(m,@pdeconvforz,@icpde,@bcpde,x,z);%uso del comando "sol", donde 
                                            %pdepe devuelve la solución en 
                                            %una matriz 3-D sol.
u=sol(:,:,1);%Para este problema, u tiene solo un componente, por lo que 
             %sol es una matriz, pero en general puede extraer 
             %el k-ésimo componente de la solución con el comando 
             %u = sol(:,:,k). Entonces se extrae el primer componente de la
             %solucion "sol".
%
%Visualizacion de los resultados (evaluacion y aplicacion):
%
figure(1) 
surf(x,z,u)
axis([0 0.05 0 5 0 10]);
title ('Ecuacion de Difusion de especies');
xlabel('Distancia en x');
ylabel ('Distancia en z');
zlabel('Concentración de la especie A u(x,z)');
%
figure(2)
waterfall(x,z,u)
axis([0 0.05 0 5 0 10]);
title ('Ecuacion de Difusion de especies');
xlabel('Distancia en x');
ylabel ('Distancia en z');
zlabel('Concentración de la especie A u(x,z)');
%
figure(3)
mesh(x,z,u)
axis([0 0.05 0 5 0 10]);
title ('Ecuacion de Difusion de especies');
xlabel('Distancia en x');
ylabel ('Distancia en z');
zlabel('Concentración de la especie A u(x,z)');
%Curvas de Nivel (Az=0 Y El=90)
figure(4)
contour(x,z,u)
axis([0 0.05 0 5 0 10]);
title ('Ecuacion de Difusion de especies');
xlabel('Distancia en x');
ylabel ('Distancia en z');
zlabel('Concentración de la especie A u(x,z)');
end
