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

Last change on this file since 530 was 530, checked in by Argimiro Resende Secchi, 14 years ago

Fix some examples.

File size: 6.5 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]
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]
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.