source: trunk/sample/miscellaneous/sample_gibbs_reactor_simple.mso @ 623

#*---------------------------------------------------------------------
* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
*
* This LIBRARY is free software; you can distribute it and/or modify
* it under the therms of the ALSOC FREE LICENSE as available at
* http://www.enq.ufrgs.br/alsoc.
*
* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
* from http://www.rps.eng.br Copyright (C) 2002-2004.
* All rights reserved.
*
* EMSO is distributed under the therms of the ALSOC LICENSE as
* available at http://www.enq.ufrgs.br/alsoc.
*
*----------------------------------------------------------------------
* Model of a simplified Gibbs reactor.
* This model requires VRTherm (www.vrtech.com.br) to run.
*----------------------------------------------------------------------
*
*
*
*----------------------------------------------------------------------
* Author: Rafael de Pelegrini Soares
* $Id$
*--------------------------------------------------------------------*#

using "types";

Model gibbs_reactor_simple
        ATTRIBUTES
        Info =
"Model for reactor in thermodynamic equilibrium based on Gibbs free
energy.

The user should specify T, P and ni.

There is a correction for G0 as a function of the temperature.
This correction considers that H0 does not depend on the temperature.
This is a good approximation for most cases (less than 2% of error)";
        
        PARAMETERS
outer PP                as Plugin (Brief="External physical properties", Type="PP");
outer NComp     as Integer (Brief="Number of components", Default=1);

        nu(NComp) as Real(Symbol="\nu_i");

        R  as Real(Brief="Universal gas constant", Unit='J/mol/K', Default=8.314);
        T0 as temperature(Default = 298.15, Symbol="T_{298}");
        P0 as Real(Default=1, Unit='bar', Symbol="P^0");
        g0(NComp) as energy_mol         (Brief="Gibbs energy in standard state", Symbol="g^0_{298,i}");
        h0(NComp) as energy_mol         (Brief="Enthalpy in standard state", Symbol="h^0_{298,i}");
        
        SET
        g0 = PP.IdealGasGibbsOfFormationAt25C();
        h0 = PP.IdealGasEnthalpyOfFormationAt25C();

        VARIABLES
        T as temperature;
        P as pressure;
        n0(NComp) as positive(Brief="Initial number of mols", Unit='mol', Symbol="n_{i,0}");
        n(NComp)  as positive(Brief="Number of mols at equilibrium", Unit='mol', Symbol="n_i");
        advance   as Real(Unit='mol', Symbol="\epsilon");

        K  as positive(Brief="Reaction equilibrium constant at T");
        K0 as positive(Brief="Reaction equilibrium constant at T0", Symbol="K_{298}");
        
        phi(NComp)      as fugacity(Brief="Fugacity coefficient", Default=1, Symbol="\phi_i");
        
        EQUATIONS
        "Equilibrium constant at 298 K"
        K0 = exp(-sum(nu*g0)/(R*T0));
        "Equilibrium constant at T"
        K = K0 * exp(-sum(nu*h0)/R*(1/T - 1/T0));

        "Equilibrium rule"
        K = prod( (n/sum(n)*phi*P/P0) ^ nu);

        "Reaction advance"
        n = n0 + nu * advance;
        
        "Fugacity coefficient"
        phi = PP.VapourFugacityCoefficient(T, P, n/sum(n));
end

# Ethane decomposition to produce ethylene and hydrogen at 1000 degC
FlowSheet gibbs_reactor_simple_sample as gibbs_reactor_simple
        PARAMETERS
        PP as Plugin(Brief="Physical Properties",
                Type="PP",
                Components = ["ethane", "ethylene", "hydrogen"],
                LiquidModel = "PR",
                VapourModel = "PR"
        );
        NComp as Integer;

        SET
        NComp = PP.NumberOfComponents;
        nu = [-1, 1, 1];

        SPECIFY
        T = (1000 + 273.15) * 'K';
        P = 1 * 'atm';
        n0 = [1, 0, 0] * 'mol';
        
        # Expected results:
        # advance = 0.9;
        # n = [0.1, 0.9, 0.9]
        
        OPTIONS
        Dynamic = false;
end


# Ethanol production by ethylene hydratation at 250 degC and 35 bar
FlowSheet gibbs_reactor_simple_sample2 as gibbs_reactor_simple
        PARAMETERS
        PP as Plugin(Brief="Physical Properties",
                Type="PP",
                Components = ["ethylene", "water", "ethanol"],
                LiquidModel = "IdealLiquid",
                VapourModel = "PR"
        );
        NComp as Integer;

        SET
        NComp = PP.NumberOfComponents;
        nu = [-1, -1, 1];

        SPECIFY
        T = (250 + 273.15) * 'K';
        P = 35 * 'atm';
        n0 = [1, 5, 0] * 'mol';
        
        # Expected results:
        # advance = 0.21;
        
        OPTIONS
        Dynamic = false;
end

# Water-gas-shift T = 1100 K, P = 1 bar
FlowSheet gibbs_reactor_simple_sample3 as gibbs_reactor_simple
        PARAMETERS
        PP as Plugin(Brief="Physical Properties",
                Type="PP",
                Components = ["carbon monoxide", "water", "carbon dioxide", "hydrogen"],
                LiquidModel = "IdealLiquid",
                VapourModel = "PR"
        );
        NComp as Integer;

        SET
        NComp = PP.NumberOfComponents;
        nu = [-1, -1, 1, 1];

        SPECIFY
        T = 1100 * 'K';
        P = 1 * 'bar';
        n0 = [1, 1, 0, 0] * 'mol';
        
        # Expected results:
        # advance = 0.41;
        
        OPTIONS
        Dynamic = false;
end


# Water-gas-shift T = 1100 K, P = 10 bar
FlowSheet gibbs_reactor_simple_sample4 as gibbs_reactor_simple
        PARAMETERS
        PP as Plugin(Brief="Physical Properties",
                Type="PP",
                Components = ["carbon monoxide", "water", "carbon dioxide", "hydrogen"],
                LiquidModel = "IdealLiquid",
                VapourModel = "PR"
        );
        NComp as Integer;

        SET
        NComp = PP.NumberOfComponents;
        nu = [-1, -1, 1, 1];

        SPECIFY
        T = 1100 * 'K';
        P = 1 * 'bar';
        n0 = [1, 1, 0, 0] * 'mol';
        
        # Expected results:
        # advance = 0.41;
        
        OPTIONS
        Dynamic = false;
end


# Water-gas-shift T = 1100 K, P = 1 bar excess of water
FlowSheet gibbs_reactor_simple_sample5 as gibbs_reactor_simple
        PARAMETERS
        PP as Plugin(Brief="Physical Properties",
                Type="PP",
                Components = ["carbon monoxide", "water", "carbon dioxide", "hydrogen"],
                LiquidModel = "IdealLiquid",
                VapourModel = "PR"
        );
        NComp as Integer;

        SET
        NComp = PP.NumberOfComponents;
        nu = [-1, -1, 1, 1];

        SPECIFY
        T = 1100 * 'K';
        P = 1 * 'bar';
        n0 = [1, 2, 0, 0] * 'mol';
        
        # Expected results:
        # advance = 0.56;
        
        OPTIONS
        Dynamic = false;
end

# Ammonia synthesis from nitrogen and hydrogen T = 500 degC, P = 1 bar
FlowSheet gibbs_reactor_simple_sample6 as gibbs_reactor_simple
        PARAMETERS
        PP as Plugin(Brief="Physical Properties",
                Type="PP",
                Components = ["nitrogen", "hydrogen", "ammonia"],
                LiquidModel = "IdealLiquid",
                VapourModel = "PR"
        );
        NComp as Integer;

        SET
        NComp = PP.NumberOfComponents;
        nu = [-1, -1, 2];

        SPECIFY
        T = (500 + 273.15) * 'K';
        P = 1 * 'bar';
        n0 = [1, 3, 0] * 'mol';
        
        # Expected results:
        # K = 6e-5;
        
        OPTIONS
        Dynamic = false;
end

# Ammonia synthesis from nitrogen and hydrogen T = 500 degC, P = 300 bar
FlowSheet gibbs_reactor_simple_sample7 as gibbs_reactor_simple
        PARAMETERS
        PP as Plugin(Brief="Physical Properties",
                Type="PP",
                Components = ["nitrogen", "hydrogen", "ammonia"],
                LiquidModel = "IdealLiquid",
                VapourModel = "PR"
        );
        NComp as Integer;

        SET
        NComp = PP.NumberOfComponents;
        nu = [-1, -3, 2];

        SPECIFY
        T = (500 + 273.15) * 'K';
        P = 300 * 'bar';
        n0 = [1, 3, 0] * 'mol';
        
        # Expected results:
        # advance = 0.54;
        
        OPTIONS
        Dynamic = false;
end