source: mso/eml/reactors/pfr.mso @ 68

#*---------------------------------------------------------------------
* This file is property of the author and cannot be used, copyed
* or modified without permission.
*
* Copyright (C) 2004  the author
*----------------------------------------------------------------------
* Author: Rafael de P. Soares
*         Paula B. Staudt
* $Id: cstr.mso 1 2006-06-20 17:33:53Z rafael $
*--------------------------------------------------------------------*#

using "streams";

#*
* Generic PFR model with constant mass holdup
*
* Requires the information of:
* - Reaction values
* - Heat of reaction
* - Pressure profile
*#
Model pfr
        PARAMETERS
ext PP   as CalcObject (Brief = "External Physical Properties");
ext     NComp as Integer(Brief="Number of components");
        NReac as Integer(Brief="Number of reactions");
        stoic(NComp, NReac) as Real (Brief = "Stoichiometric Matrix");
        NDisc as Integer(Brief="Number of points of discretization", Default=10);
        Mw(NComp)  as molweight (Brief="Component Mol Weight");
        
        L as length(Brief="Reactor Length");
        Across as area(Brief="Cross section area");

        SET
        Mw   = PP.MolecularWeight();

        VARIABLES
in      Inlet  as stream;
out     Outlet as stream;
        str(NDisc+1) as stream_therm;
        vel(NDisc+1) as velocity;
        vol(NDisc+1) as vol_mol;
        rho(NDisc+1) as dens_mass;
        
        q(NDisc)    as heat_rate;
        M(NComp, NDisc)    as mol (Brief = "Molar holdup");
        C(NComp, NDisc)  as conc_mol(Brief="Components concentration");
        E(NDisc) as energy (Brief="Total Energy Holdup on element");
        
        r(NReac, NDisc) as reaction_mol;
        Hr(NReac, NDisc) as heat_reaction;

        EQUATIONS
        "Vapourisation Fraction"
        str.v = Inlet.v;

        "Inlet boundary"
        str(1).F = Inlet.F;
        str(1).T = Inlet.T;
        str(1).P = Inlet.P;
        str(1).z = Inlet.z;

        "Outlet boundary"
        Outlet.F = str(NDisc+1).F;
        Outlet.T = str(NDisc+1).T;
        Outlet.P = str(NDisc+1).P;
        Outlet.z = str(NDisc+1).z;
        Outlet.h = str(NDisc+1).h;
        Outlet.v = str(NDisc+1).v;
        
        for z in [1:NDisc]
                for c in [1:NComp]
                        "Component Molar Balance"
                        diff(M(c,z)) = (str(z).F*str(z).z(c) - str(z+1).F*str(z+1).z(c))
                                + sum(stoic(c,:)*r(:, z)) * Across*L/NDisc;
                end

                "Energy Balance"
                diff(E(z)) = str(z).F*str(z).h - str(z+1).F*str(z+1).h +
                        sum(Hr(:,z)*r(:,z)) * Across*L/NDisc - q(z);

                "Energy Holdup"
                E(z) = sum(M(:,z))*str(z+1).h - str(z+1).P*Across*L/NDisc;

                "mass flow is considered constant"
                str(z+1).F*vol(z+1) = str(z).F*vol(z); # FIXME: is this correct? No (constant velocity: only for equimolar)
#               rho(z+1)*vel(z+1) = rho(z)*vel(z);  # FIXME: this is correct! But does not converge.
        
                "Molar concentration"
                C(:,z) * Across*L/NDisc = M(:,z);

                "Molar fraction"
                str(z+1).z = C(:,z) * vol(z+1);
        end     

        for z in [1:NDisc+1]
                "Specific Volume"
                vol(z) = PP.VapourVolume(str(z).T, str(z).P, str(z).z);
                
                "Specific Mass"
                rho(z) = sum(str(z).z*Mw)/vol(z);

                "Velocity"
                vel(z)*Across = str(z).F*vol(z);
        end
end