source: branches/gui/eml/reactors/equil.mso @ 574

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

Updated the models to work with some language constraints

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 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] do
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] do
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] do
114          if (Outletm.z(i) > 1e-16) 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] do
195        "Equilibrium constant"
196                K(j) = prod(activ^stoic(:,j));
197        end
198       
199        for i in [1:NComp] do
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] do
205          if (Outletm.z(i) > 1e-16) 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.