source: trunk/eml/reactors/equil.mso @ 471

Last change on this file since 471 was 471, checked in by Rodolfo Rodrigues, 15 years ago

Updated equil and gibbs models

File size: 6.7 KB
Line 
1#*---------------------------------------------------------------------
2* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
3*
4* This LIBRARY is free software; you can distribute it and/or modify
5* it under the therms of the ALSOC FREE LICENSE as available at
6* http://www.enq.ufrgs.br/alsoc.
7*
8* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
9* from http://www.rps.eng.br Copyright (C) 2002-2004.
10* All rights reserved.
11*
12* EMSO is distributed under the therms of the ALSOC LICENSE as
13* available at http://www.enq.ufrgs.br/alsoc.
14*
15*----------------------------------------------------------------------
16* Model of an equilibrium reactor
17*----------------------------------------------------------------------
18*
19*   Description:
20*       Thermodynamic equilibrium modeling of a reactor based on
21*       equilibrium constants approach.
22*
23*   Assumptions:
24*               * single-phases involved
25*       * thermodynamic equilibrium
26*               * steady-state
27*
28*       Specify:
29*               * inlet stream
30*               * stoichiometric matrix
31*               * equilibrium temperature
32*
33*----------------------------------------------------------------------
34* Author: Rodolfo Rodrigues
35* $Id$
36*--------------------------------------------------------------------*#
37
38using "tank_basic";
39
40
41#*---------------------------------------------------------------------
42*       only vapour phase
43*--------------------------------------------------------------------*#
44Model equil_vap as tank_vap
45        ATTRIBUTES
46        Pallete         = true;
47        Icon            = "icon/cstr";
48        Brief           = "Model of a generic vapour-phase equilibrium CSTR";
49        Info            = "
50== Assumptions ==
51* only vapour-phase
52* thermodynamic equilibrium
53* steady-state
54
55== Specify ==
56* inlet stream
57* stoichiometric matrix
58* equilibrium temperature
59";
60
61        PARAMETERS
62        NReac           as Integer              (Brief="Number of reactions", Default=1);
63    stoic(NComp,NReac) as Real  (Brief="Stoichiometric matrix", Symbol="\nu");
64        Rg                      as Real                 (Brief="Universal gas constant", Unit='J/mol/K', Default=8.314);
65        fs(NComp)       as pressure     (Brief="Fugacity in standard state", Default=1, DisplayUnit='atm');
66        To                      as temperature  (Brief="Reference temperature", Default=298.15);
67
68        VARIABLES
69out Outlet              as vapour_stream(Brief="Outlet stream", PosX=1, PosY=1, Symbol="_{out}");
70
71        G(NComp)        as energy_mol   (Brief="Gibbs free-energy of formation");
72        K(NReac)        as Real                 (Brief="Equillibrium constant", Lower=0, Default=1.5);
73        activ(NComp)as Real             (Brief="Activity", Symbol="\hat{a}", Lower=0, Default=0.2);
74       
75        rate(NComp) as reaction_mol (Brief="Overall component rate of reaction");
76        extent(NReac) as flow_mol       (Brief="Extent of reaction", Symbol="\xi");
77        conv(NComp) as Real             (Brief="Fractional conversion of component", Symbol="X", Default=0); # Lower=-1e3, Upper=1e3);
78       
79        EQUATIONS
80        "Outlet stream"
81        Outlet.F*Outlet.z = Outletm.F*Outletm.z + rate*Tank.V;
82       
83        "Mechanical equilibrium"
84        Outlet.P = Outletm.P;
85       
86        "Energy balance"
87        Outlet.F*Outlet.h = Outletm.F*Outletm.h;
88       
89        "Steady-state"
90        Outlet.F = Inlet.F + sum(sumt(stoic*extent));
91
92        "Gibbs free-energy of formation"
93        G = PP.IdealGasGibbsOfFormation(Outlet.T);
94
95#       "Gibbs free-energy of formation without Cp correction"
96#       G = PP.IdealGasGibbsOfFormationAt25C()*Outlet.T/To
97#               + PP.IdealGasEnthalpyOfFormationAt25C()*(1 - Outlet.T/To);
98
99        for j in [1:NReac]
100        "Gibbs free energy of reaction"
101                sumt(G*stoic(:,j)) = -Rg*Outlet.T*ln(K(j));
102#               K(j) = exp(-sumt(G*stoic(:,j))/(Rg*Outlet.T));
103
104        "Equilibrium constant"
105                K(j) = prod(activ^stoic(:,j));
106        end
107
108        for i in [1:NComp]
109        "Outlet molar fraction"
110                Outlet.F*Outlet.z(i) = (Inlet.F*Inlet.z(i) + sumt(stoic(i,:)*extent));
111        end     
112
113        for i in [1:NComp]
114          if (Outletm.z(i) > 0) then
115            "Molar conversion"
116            Outlet.F*Outlet.z(i) = Outletm.F*Outletm.z(i)*(1 - conv(i));
117          else if (Outlet.z(i) > 0) then
118                        "Molar conversion"
119                                conv(i) = 1;    # ?
120                        else
121                        "Molar conversion"
122                                conv(i) = 0;    # ?
123                        end
124          end
125        end
126
127        "Activity"
128        activ = PP.VapourFugacityCoefficient(Outlet.T,Outlet.P,Outlet.z)*Outlet.P*Outlet.z/fs;
129end
130
131
132#*---------------------------------------------------------------------
133*       only liquid-phase
134*--------------------------------------------------------------------*#
135Model equil_liq as tank_liq
136        ATTRIBUTES
137        Pallete         = true;
138        Icon            = "icon/cstr";
139        Brief           = "Model of a generic liquid-phase equilibrium CSTR";
140        Info            = "
141== Assumptions ==
142* only liquid-phase
143* thermodynamic equilibrium
144* steady-state
145
146== Specify ==
147* inlet stream
148* stoichiometric matrix
149* equilibrium temperature
150";
151
152        PARAMETERS
153        NReac           as Integer              (Brief="Number of reactions", Default=1);
154    stoic(NComp,NReac) as Real  (Brief="Stoichiometric matrix", Symbol="\nu");
155        Rg                      as Real                 (Brief="Universal gas constant", Unit='J/mol/K', Default=8.314);
156        Ps                      as pressure     (Brief="Standard pressure", Default=1, DisplayUnit='bar');
157        To                      as temperature  (Brief="Reference temperature", Default=298.15);
158
159        VARIABLES
160out Outlet              as liquid_stream(Brief="Outlet stream", PosX=1, PosY=1, Symbol="_{out}");
161
162        G(NReac)        as enth_mol     (Brief="Gibbs free-energy of formation");
163        K(NReac)        as fraction             (Brief="Equillibrium constant");
164        activ(NComp)as Real             (Brief="Activity", Symbol="\hat{a}");
165       
166        rate(NComp) as reaction_mol (Brief="Overall component rate of reaction");
167        extent(NReac)as flow_mol        (Brief="Extent of reaction", Symbol="\xi");
168        conv(NComp) as Real             (Brief="Fractional conversion of component", Symbol="X", Default=0);
169       
170        EQUATIONS
171        "Outlet stream"
172        Outlet.F*Outlet.z = Outletm.F*Outletm.z + rate*Tank.V;
173       
174        "Mechanical equilibrium"
175        Outlet.P = Outletm.P;
176       
177        "Energy balance"
178        Outlet.F*Outlet.h = Outletm.F*Outletm.h;
179       
180        "Steady-state"
181        Outlet.F = Inlet.F + sum(sumt(stoic*extent));
182       
183        "Gibbs free-energy of formation"
184        G = PP.IdealGasGibbsOfFormation(Outlet.T);
185
186#       "Gibbs free-energy of formation without Cp correction"
187#       G = PP.IdealGasGibbsOfFormationAt25C()*Outlet.T/To
188#               + PP.IdealGasEnthalpyOfFormationAt25C()*(1 - Outlet.T/To);
189
190        "Gibbs free energy of reaction"
191        sumt(G*stoic) = -Rg*Outlet.T*ln(K);
192#       K = exp(-sumt(G*stoic)/(Rg*Outlet.T));
193       
194        for j in [1:NReac]
195        "Equilibrium constant"
196                K(j) = prod(activ^stoic(:,j));
197        end
198       
199        for i in [1:NComp]
200        "Outlet molar fraction"
201                Outlet.F*Outlet.z(i) = (Inlet.F*Inlet.z(i) + sumt(stoic(i,:)*extent));
202        end     
203
204        for i in [1:NComp]
205          if (Outletm.z(i) > 0) then
206            "Molar conversion"
207            Outlet.F*Outlet.z(i) = Outletm.F*Outletm.z(i)*(1 - conv(i));
208          else if (Outlet.z(i) > 0) then
209                        "Molar conversion"
210                                conv(i) = 1;    # ?
211                        else
212                        "Molar conversion"
213                                conv(i) = 0;    # ?
214                        end
215          end
216        end
217       
218        "Activity"
219        activ = PP.LiquidFugacityCoefficient(Outlet.T,Outlet.P,Outlet.z)*Outlet.z
220                *exp(PP.LiquidVolume(Outlet.T,Outlet.P,Outlet.z)*(Outlet.P - Ps)/Rg/Outlet.T);
221end
Note: See TracBrowser for help on using the repository browser.