source: trunk/eml/reactors/gibbs.mso @ 972

Last change on this file since 972 was 574, checked in by Rafael de Pelegrini Soares, 14 years ago

Updated the models to work with some language constraints

File size: 6.3 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 a Gibbs reactor
17*----------------------------------------------------------------------
18*
19*   Description:
20*       Thermodynamic equilibrium modeling of a reactor using Gibbs
21*       free energy minimization approach.
22*
23*   Assumptions:
24*               * single-phases involved
25*       * thermodynamic equilibrium
26*               * steady-state
27*
28*       Specify:
29*               * inlet stream
30*               * number of elements related to components
31*               * matrix of elements by components
32*               * equilibrium temperature
33*
34*----------------------------------------------------------------------
35* Author: Rodolfo Rodrigues
36* $Id$
37*--------------------------------------------------------------------*#
38
39using "tank_basic";
40
41
42#*---------------------------------------------------------------------
43*       only vapour phase
44*--------------------------------------------------------------------*#
45Model gibbs_vap as tank_vap
46        ATTRIBUTES
47        Pallete = true;
48        Icon    = "icon/cstr";
49        Brief   = "Model of a generic vapour-phase Gibbs CSTR";
50        Info    = "
51== Assumptions ==
52* thermodynamic equilibrium
53* steady-state
54
55== Specify ==
56* inlet stream
57* number of elements related to components
58* matrix of elements by components
59* equilibrium temperature
60";
61
62        PARAMETERS
63outer NElem                     as Integer              (Brief="Number of elements", Default=1);
64        Rg                              as Real                 (Brief="Universal gas constant", Unit='J/mol/K', Default=8.314);
65        na(NElem,NComp) as Real                 (Brief="Number of elements per component");
66        fs(NComp)               as pressure     (Brief="Fugacity in standard state", Default=1, DisplayUnit='atm');
67        To                              as temperature  (Brief="Reference temperature", Default=298.15);
68       
69        VARIABLES
70out Outlet                      as vapour_stream(Brief="Outlet stream", PosX=1, PosY=1, Symbol="_{out}");
71
72        G(NComp)                as energy_mol   (Brief="Gibbs free-energy change of formation");
73        lambda(NElem)   as energy_mol   (Brief="Lagrangian multiplier", Symbol="\lambda");
74        activ(NComp)    as Real                 (Brief="Activity", Symbol="\hat{a}", Lower=0);
75       
76        rate(NComp)     as reaction_mol (Brief="Overall component rate of reaction");
77        conv(NComp)     as Real                 (Brief="Fractional conversion of component", Symbol="X", Default=0);
78        Fi(NComp)               as flow_mol             (Brief="Component molar flow rate");
79
80        EQUATIONS
81        "Outlet stream"
82        Outlet.F*Outlet.z = Outletm.F*Outletm.z + rate*Tank.V;
83       
84        "Mechanical equilibrium"
85        Outlet.P = Outletm.P;
86
87        "Steady-state"
88        Outlet.F = sum(Fi);
89
90        "Component molar flow rate"
91        Fi = Outlet.F*Outlet.z;
92       
93        "Energy balance"
94        Outlet.F*Outlet.h = Outletm.F*Outletm.h;
95
96        "Element balance"
97        sumt(Fi*na) = sumt(Outletm.F*Outletm.z*na);
98
99        "Gibbs free-energy of formation"
100        G = PP.IdealGasGibbsOfFormation(Outlet.T);
101
102#       "Gibbs free-energy of formation without Cp correction"
103#       G = PP.IdealGasGibbsOfFormationAt25C()*Outlet.T/To
104#               + PP.IdealGasEnthalpyOfFormationAt25C()*(1 - Outlet.T/To);
105
106        for i in [1:NComp] do
107        "Lagrangian multiplier"
108                G(i) + sumt(lambda*na(:,i)) = -Rg*Outlet.T*ln(activ(i));
109
110          if (Outletm.z(i) > 1e-16) then
111            "Molar conversion"
112            Fi(i) = Outletm.F*Outletm.z(i)*(1 - conv(i));
113          else if (Outlet.z(i) > 0) then
114                        "Molar conversion"
115                                conv(i) = 1;    # ?
116                        else
117                        "Molar conversion"
118                                conv(i) = 0;    # ?
119                        end
120          end
121        end
122       
123        "Activity"
124        activ = PP.VapourFugacityCoefficient(Outlet.T,Outlet.P,Outlet.z)
125                *Outlet.P*Outlet.z/fs;
126end
127
128
129#*---------------------------------------------------------------------
130*       only liquid phase
131*--------------------------------------------------------------------*#
132Model gibbs_liq as tank_liq
133        ATTRIBUTES
134        Pallete = true;
135        Icon    = "icon/cstr";
136        Brief   = "Model of a generic liquid-phase Gibbs CSTR";
137        Info    = "
138== Assumptions ==
139* thermodynamic equilibrium
140* steady-state
141
142== Specify ==
143* inlet stream
144* number of elements related to components
145* matrix of elements by components
146* equilibrium temperature
147";
148
149        PARAMETERS
150outer NElem                     as Integer              (Brief="Number of elements", Default=1);
151        Rg                              as Real                 (Brief="Universal gas constant", Unit='J/mol/K', Default=8.314);
152        na(NElem,NComp) as Real                 (Brief="Number of elements per component");
153        Ps                              as pressure             (Brief="Pressure of standard state", Default=1, DisplayUnit='atm');
154        To                              as temperature  (Brief="Reference temperature", Default=298.15);
155       
156        VARIABLES
157out Outlet                      as liquid_stream(Brief="Outlet stream", PosX=1, PosY=1, Symbol="_{out}");
158
159        G(NComp)                as energy_mol   (Brief="Gibbs free-energy change of formation");
160        lambda(NElem)   as energy_mol   (Brief="Lagrangian multiplier", Symbol="\lambda");
161        activ(NComp)    as Real                 (Brief="Activity", Symbol="\hat{a}", Lower=0);
162       
163        rate(NComp)     as reaction_mol (Brief="Overall component rate of reaction");
164        conv(NComp)     as Real                 (Brief="Fractional conversion of component", Symbol="X", Default=0);
165        Fi(NComp)               as flow_mol             (Brief="Component molar flow rate");
166
167        EQUATIONS
168        "Outlet stream"
169        Outlet.F*Outlet.z = Outletm.F*Outletm.z + rate*Tank.V;
170       
171        "Mechanical equilibrium"
172        Outlet.P = Outletm.P;
173
174        "Steady-state"
175        Outlet.F = sum(Fi);
176
177        "Component molar flow rate"
178        Fi = Outlet.F*Outlet.z;
179       
180        "Energy balance"
181        Outlet.F*Outlet.h = Outletm.F*Outletm.h;
182
183        "Element balance"
184        sumt(Fi*na) = sumt(Outletm.F*Outletm.z*na);
185       
186        "Gibbs free-energy of formation"
187        G = PP.IdealGasGibbsOfFormation(Outlet.T);
188
189#       "Gibbs free-energy of formation without Cp correction"
190#       G = PP.IdealGasGibbsOfFormationAt25C()*Outlet.T/To
191#               + PP.IdealGasEnthalpyOfFormationAt25C()*(1 - Outlet.T/To);
192
193        for i in [1:NComp] do
194        "Lagrangian multiplier"
195                G(i) + sumt(lambda*na(:,i)) = -Rg*Outlet.T*ln(activ(i));
196
197          if (Outletm.z(i) > 1e-16) then
198            "Molar conversion"
199            Fi(i) = Outletm.F*Outletm.z(i)*(1 - conv(i));
200          else if (Outlet.z(i) > 0) then
201                        "Molar conversion"
202                                conv(i) = 1;    # ?
203                        else
204                        "Molar conversion"
205                                conv(i) = 0;    # ?
206                        end
207          end
208        end
209       
210        "Activity"
211        activ = PP.LiquidFugacityCoefficient(Outlet.T,Outlet.P,Outlet.z)*Outlet.z
212                *exp(PP.LiquidVolume(Outlet.T,Outlet.P,Outlet.z)*(Outlet.P - Ps)/Rg/Outlet.T);
213end
Note: See TracBrowser for help on using the repository browser.