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

Last change on this file since 69 was 69, checked in by Rafael de Pelegrini Soares, 16 years ago

Updated PFR model

File size: 3.0 KB
Line 
1#*---------------------------------------------------------------------
2* This file is property of the author and cannot be used, copyed
3* or modified without permission.
4*
5* Copyright (C) 2004  the author
6*----------------------------------------------------------------------
7* Author: Rafael de P. Soares
8*         Paula B. Staudt
9* $Id: cstr.mso 1 2006-06-20 17:33:53Z rafael $
10*--------------------------------------------------------------------*#
11
12using "streams";
13
14#*
15* Generic PFR model with constant mass holdup
16*
17* Requires the information of:
18* - Reaction values
19* - Heat of reaction
20* - Pressure profile
21*#
22Model pfr
23        PARAMETERS
24ext PP   as CalcObject (Brief = "External Physical Properties");
25ext     NComp as Integer(Brief="Number of components");
26        NReac as Integer(Brief="Number of reactions");
27        stoic(NComp, NReac) as Real (Brief = "Stoichiometric Matrix");
28        NDisc as Integer(Brief="Number of points of discretization", Default=10);
29        Mw(NComp)  as molweight (Brief="Component Mol Weight");
30       
31        L as length(Brief="Reactor Length");
32        Across as area(Brief="Cross section area");
33
34        SET
35        Mw   = PP.MolecularWeight();
36
37        VARIABLES
38in      Inlet  as stream;
39out     Outlet as stream;
40        str(NDisc+1) as stream_therm;
41        vel(NDisc+1) as velocity;
42        vol(NDisc+1) as vol_mol;
43        rho(NDisc+1) as dens_mass;
44       
45        q(NDisc)    as heat_rate;
46        M(NComp, NDisc)    as mol (Brief = "Molar holdup");
47        Mt(NDisc)    as mol (Brief = "Molar holdup");
48        C(NComp, NDisc)  as conc_mol(Brief="Components concentration", Lower=-1e-6);
49        E(NDisc) as energy (Brief="Total Energy Holdup on element");
50       
51        r(NReac, NDisc) as reaction_mol;
52        Hr(NReac, NDisc) as heat_reaction;
53
54        EQUATIONS
55        "Vapourisation Fraction"
56        str.v = Inlet.v;
57
58        "Inlet boundary"
59        str(1).F = Inlet.F;
60        str(1).T = Inlet.T;
61        str(1).P = Inlet.P;
62        str(1).z = Inlet.z;
63
64        "Outlet boundary"
65        Outlet.F = str(NDisc+1).F;
66        Outlet.T = str(NDisc+1).T;
67        Outlet.P = str(NDisc+1).P;
68        Outlet.z = str(NDisc+1).z;
69        Outlet.h = str(NDisc+1).h;
70        Outlet.v = str(NDisc+1).v;
71       
72        for z in [1:NDisc]
73                for c in [1:NComp-1]
74                        "Component Molar Balance"
75                        diff(M(c,z)) = (str(z).F*str(z).z(c) - str(z+1).F*str(z+1).z(c))
76                                + sum(stoic(c,:)*r(:, z)) * Across*L/NDisc;
77                end
78
79                "Energy Balance"
80                diff(E(z)) = str(z).F*str(z).h - str(z+1).F*str(z+1).h +
81                        sum(Hr(:,z)*r(:,z)) * Across*L/NDisc - q(z);
82
83                "Energy Holdup"
84                E(z) = Mt(z)*str(z+1).h - str(z+1).P*Across*L/NDisc;
85
86                "mass flow is considered constant"
87                str(z+1).F*vol(z+1) = str(z).F*vol(z); # FIXME: is this correct? No (constant velocity: only for equimolar)
88                #rho(z+1)*vel(z+1) = rho(z)*vel(z);  # FIXME: this is correct! But does not converge.
89       
90                "Molar concentration"
91                C(:,z) * Across*L/NDisc = M(:,z);
92               
93                "Sum of M"
94                Mt(z) = sum(M(:,z));
95               
96                "Geometrical constraint"
97                Across*L/NDisc = Mt(z) * vol(z);
98
99                "Molar fraction"
100                str(z+1).z * Mt(z) = M(:,z);
101        end
102
103        for z in [1:NDisc+1]
104                "Specific Volume"
105                vol(z) = PP.VapourVolume(str(z).T, str(z).P, str(z).z);
106               
107                "Specific Mass"
108                rho(z) = PP.VapourDensity(str(z).T, str(z).P, str(z).z);
109
110                "Velocity"
111                vel(z)*Across = str(z).F*vol(z);
112        end
113end
Note: See TracBrowser for help on using the repository browser.