function absorcion()

% Absorcion Isoterma. Fenomenos de Transporte 2004 - 2005
% Determinacion en regimen permanente de los flujos de difusion de
% sulfuro de hidrogeno, dioxido de carbono, metano y etano en una
% seccion de una columna de relleno.
% Por: Jose Luis de la Mata San Marcos 01305

% Sistema gaseoso:
%  - 1: Sulfuro de hidrogeno
%  - 2: Dioxido de carbono
%  - 3: Metano
%  - 4: Etano

% Sistema liquido:
%  - 1: Sulfuro de hidrogeno
%  - 2: Dioxido de carbono
%  - 3: Agua

% ********************************************************************

yI0 = [0.1 0.1 0.1];
yI = fsolve(@residuo,yI0,0);
yI
k_E = [60.92 185.6]; % Constante de equilibrio [sulfuro dioxido]
xI(1) = yI(1)/k_E(1);
xI(2) = yI(2)/k_E(2);
xI
FG = flujo_gas(yI)
FL = flujo_liquido(yI)

% ********************************************************************

function r = residuo(yI)

    FG = flujo_gas(yI);
    FL = flujo_liquido(yI);
    
    r(1) = FG(1)+FL(1);
    r(2) = FG(2)+FL(2);
    r(3) = FG(3);

% ********************************************************************

function FG = flujo_gas(yI)

% Parametros del sistema

P = 1.013e6; % Presion en Pa
T = 30 + 273; % Temperatura del sistema en K
y = [0.00674 0.16465 0.1067 0.72191]; % Composicion gases fraccion molar
g = 9.81; % Aceleracion de la gravedad en m/s^2
k_E = [60.92 185.6]; % Constante de equilibrio [sulfuro dioxido]
G = 125;    % Caudal del gas en kmol/h
D = 1.2; % Diametro de la columna en m
A = pi*D^2/4; % Seccion de la columna en m^2
visG = 0.0000085;    % Viscosidad en Pa·s


% Determinacion de las difusividades del sistema gaseoso

    PMG = [34 44 16 30]; % Pesos moleculares g/mol
    pmG = PMG*y'; % Peso molecular del gas
    VDG = [27.52 26.7 25.14 45.66]; % Volumenes de difusion

    % Difusividades de Maxwell-Stefan:
    for i = 1:4
        for j = 1:4
            % Peso molecular medio para la ecuacion de Fuller:
            PMM = 2/((1/PMG(i))+(1/PMG(j)));
            % Difusividades de Maxwell-Stefan:
            DMS(i,j) = 1.013e-2*T^1.75;
            DMS(i,j) = DMS(i,j)/(PMM^0.5*P*(VDG(i)^(1/3)+VDG(j)^(1/3))^2);
            % Valores de difusividades en m^2/s
        end
    end    
     

% Determinacion de los coeficientes de transferencia para el gas
% empleando la correlacion de Onda

    ap = 190; % Relacion superficie/volumen del relleno en m^2/m^3
    wG = G*pmG/(A*3600); % velocidad masica superficial en kg*m^-2*s^-1
    a = 0.76*ap; % En m^2/m^3
    dp = 25e-3; % Tamaño del relleno en m
    C = 5.23;   % Parametro de la correlacion dependiente de dp
    R = 8.314;  % Unidades SI
    denG = (pmG*P)/(R*T*1000);  % Densidad del gas en kg/m^3
    
    % Coeficientes de transferencia binarios
    for i = 1:4
        for j = 1:4
            kc(i,j) = DMS(i,j)*ap*C*(wG/(visG*ap))^0.7;
            kc(i,j) = kc(i,j)*(visG/(denG*DMS(i,j)))^(0.33);
            kc(i,j) = kc(i,j)*(dp*ap)^-2;
        end
    end
    
    % Matriz de coeficientes de transferencia
    
    for i = 1:3
        for j = 1:3
            if i == j
                R(i,j) = y(i)/kc(i,4);
                for k = 1:4
                    if i ~= k
                        R(i,j) = R(i,j) + y(k)/kc(i,k);
                    end
                end
            else
                R(i,j) = -y(i)*((1/kc(i,j))-(1/kc(i,4)));
            end
        end
    end
    K = R^-1;
    
    % Matriz beta
    
    for i = 1:3
        for j = 1:3
            if i == j
                beta(i,j) = 1 + y(i)/y(4);
            else
                beta(i,j) = y(i)/y(4);
            end
        end
    end
    
    c = denG/pmG;
    
    % Iteracion para el flujo
    cont = 0;
    dif = [1 1 1];
    tol = [1 1 1]*1e-9;
 
    for i = 1:3
         grad_y(i) = yI(i)-y(i);
    end
    grad_y = grad_y';
    
    N = c*beta*K*grad_y;
    N(4) = 0;
    Nt = N(1) + N(2) + N(3);
    
    while (dif > tol) & (cont < 1000)
    
        for i = 1:3
            for j = 1:3
                if i == j
                    Fi(i,j) = N(i)/(c*kc(i,4));
                    for k = 1:4
                        if i ~= k
                            Fi(i,j) = Fi(i,j) + N(k)/(c*kc(i,k));
                        end
                    end
                else
                    Fi(i,j) = -N(i)*((1/(c*kc(i,j))-(1/(c*kc(i,4)))));
                end
            end
        end
    
        I = [1 0 0;0 1 0;0 0 1];
        Xi = Fi*expm(Fi)*((expm(Fi)-I)^-1);
        J = c*K*Xi*grad_y;
        N1 = beta*J;
        N1(4) = 0;
        dif(1) = abs(N(1)-N1(1));
        dif(2) = abs(N(2)-N1(2));
        dif(3) = abs(N(3)-N1(3));
        N = N1;
        clear N1;
        cont = cont + 1;
    end
    
    clear grad_y;

    if cont == 100
        'No converge la iteracion'
        return
    end

FG = N;

% ********************************************************************

function FL = flujo_liquido(yI)

% Parametros del sistema

    T = 30 + 273; % Temperatura del sistema en K
    x = [0.00005454 0.0008495 0.9936965]; % Composicion liquido fraccion molar
    g = 9.81; % Aceleracion de la gravedad en m/s^2
    k_E = [60.92 185.6]; % Constante de equilibrio [sulfuro dioxido]
    L = 15500;  % Caudal del liquido en kmol/h
    D = 1.2;    % Diametro de la columna en m
    A = (pi*D^2)/4; % Seccion de la columna en m^2
    pmL = 18;   % Peso molecular del liquido en g/mol
    denL = 983.7905;  % Densidad del liquido en kg/m^3


% Determinacion de las difusividades del sistema liquido

    vis_L = 0.789e-3; % Viscosidad del agua en Pa·s a 30ºC
    vol = [32.9 34]; % Volumen molecular en el punto de ebullicion en cc/mol

    % Difusividad para el sulfuro de hidrogeno
    DFL(1) = 7.4e-8*T*((2.26*pmL)^0.5)/(vis_L*1000*vol(1)^0.6);
    DFL(1) = DFL(1)/10000; % Difusividad en m^2/s

    % Difusividad para el dioxido de carbono
    DFL(2) = 7.4e-8*T*((2.26*pmL)^0.5)/(vis_L*1000*vol(2)^0.6);
    DFL(2) = DFL(2)/10000; % Difusividad en m^2/s


% Determinacion de los coeficientes de transferencia para el liquido
% empleando la correlacion de Mohunta

    den_L = 983.7905; % Densidad del agua a 30ºC en kg/m^3
    ap = 190; % Relacion superficie/volumen del relleno en m^2/m^3
    wL = L*pmL/(A*3600); % velocidad masica superficial en kg*m^-2*s^-1
    a = 0.76*ap; % En m^2/m^3

    % Para el sulfuro de hidrogeno
    kL(1)=(2.5e-3)*(den_L*g/ap/vis_L)^(2/3)*(g^2*den_L/vis_L)^(1/9);
    kL(1)=kL(1)*(wL^3*vis_L*ap^3/g^2/den_L^4)^(1/4)*(vis_L/den_L/DFL(1))^-0.5;
    kL(1)=kL(1)/a; % En m/s
    
    % Para el dioxidod de carbono
    kL(2)=(2.5e-3)*(den_L*g/ap/vis_L)^(2/3)*(g^2*den_L/vis_L)^(1/9);
    kL(2)=kL(2)*(wL^3*vis_L*ap^3/g^2/den_L^4)^(1/4)*(vis_L/den_L/DFL(2))^-0.5;
    kL(2)=kL(2)/a; % En m/s


% Determinacion del flujo:

xI(1) = yI(1)/k_E(1);
xI(2) = yI(2)/k_E(2);

    % Matriz k
    k(1,1) = kL(1);
    k(2,2) = kL(2);
    k(1,2) = 0;
    k(2,1) = 0;
    
    % Matriz beta
    beta(1,1) = 1 + x(1)/x(3);
    beta(2,2) = 1 + x(2)/x(3);
    beta(1,2) = x(1)/x(3);
    beta(2,1) = x(2)/x(3);
    
    % Estimacion inicial del flujo
    grad_x(1) = xI(1)-x(1);
    grad_x(2) = xI(2)-x(2);
    grad_x = grad_x';
    c = denL/pmL;
    N = c*beta*k*grad_x;
    
    dif = [1;1];
    tol = [1;1]*1e-10;
    cont = 0;
    
    while (dif > tol) & (cont < 50)
        Fi(1) = (N(1)+N(2))/(c*kL(1));
        Fi(2) = (N(1)+N(2))/(c*kL(2));
    
        J(1) = c*kL(1)*(Fi(1)/(exp(Fi(1))-1))*grad_x(1);
        J(2) = c*kL(2)*(Fi(2)/(exp(Fi(2))-1))*grad_x(2);
        J = J';
    
        N1 = beta*J;
        clear J
    
        dif = abs(N-N1);
        N = N1;
        cont = cont + 1;        
    end
    
    clear grad_x;

    if cont == 50
        'No converge la iteracion'
        break
    end
    
FL = N;    