Ticket #23: sample_equil.mso

File sample_equil.mso, 4.4 KB (added by Rodolfo Rodrigues, 16 years ago)
Line 
1using "streams";
2
3#*---------------------------------------------------------------------
4* Resolved example 15.5 - Smith et al. Introduction to Chemical
5* Engineering Thermodynamics. 5th ed., 1996 - p.527 (brazilian edition)
6*
7* System in gas fase contents CO, H2O, CO2, and H2.
8*
9*       Co + H2O -> CO2 + H2
10*
11* Initial conditions
12*   1bar, 1100K, 1mol of CO, and 1mol of H2O
13* Specify
14*       Outlet.T
15*       Outlet.P
16* Return
17*       Outlet.z
18----------------------------------------------------------------------*#
19
20
21FlowSheet sample_equil
22        PARAMETERS
23        PP                      as CalcObject   (Brief="External physical properties", File="vrpp");
24        NComp           as Integer;
25        NReac           as Integer;
26       
27        stoic(NComp,NReac) as Real      (Brief="Stoichiometric matrix");
28        Rg                      as Real                 (Brief="Universal gas constant", Unit="J/mol/K", Default=8.314);
29        fs                      as pressure     (Brief="Fugacity in standard state", Default=1, Unit="atm");
30       
31       
32        VARIABLES
33        Inlet           as streamTP; # Outlet stream
34        Outlet          as streamTP; # Outlet stream
35
36        Gs(NComp)       as energy_mol   (Brief="Gibbs energy in standard state");
37        K(NReac)        as Real                 (Brief="Equillibrium constant");
38        activ(NComp)as Real             (Brief="Activity");
39       
40        extent(NReac)as flow_mol        (Brief="Extent of reaction");
41       
42       
43        EQUATIONS
44        "Change time in T"
45        Outlet.T = time*"K/s"; 
46
47        "Mechanical equilibrium"
48        Outlet.P = Inlet.P;
49
50        "Steady-state"
51        Outlet.F = Inlet.F;
52       
53        "Gibbs energy in standard state"
54        sumt(Gs*stoic) = -Rg*Outlet.T*ln(K);
55       
56        for j in [1:NReac]
57        "Equilibrium constant"
58                K(j) = prod(activ^stoic(:,j));
59        end
60
61#       K = (activ(1)^stoic(1,1))*(activ(2)^stoic(2,1))*(activ(3)^stoic(3,1))*(activ(4)^stoic(4,1));
62       
63        "Activity"
64        activ = PP.VapourFugacityCoefficient(Outlet.T,Outlet.P,Outlet.z)*Outlet.P*Outlet.z/fs;
65
66        for i in [1:NComp]
67        "Outlet molar fraction"
68                (Inlet.F + sum(sumt(stoic*extent)))*Outlet.z(i) = (Inlet.F*Inlet.z(i) + sumt(stoic(i,:)*extent));
69        end
70
71       
72        SET
73        PP.Components = ["carbon monoxide","water","carbon dioxide","hydrogen"];
74        PP.LiquidModel = "IdealLiquid";
75        PP.VapourModel = "Ideal";
76        NComp = PP.NumberOfComponents;
77
78        NReac = 1;
79        stoic(:,1) = [-1.0, -1.0, 1.0, 1.0];
80       
81       
82        SPECIFY
83        Inlet.F = 2*"mol/s";
84        Inlet.P = 1*"bar";
85        Inlet.T = 1100*"K";
86        Inlet.z = [0.5, 0.5, 0.0, 0.0];
87       
88        Gs(1) = PP.VapourGibbsFreeEnergy(Outlet.T,Outlet.P,[1, 0, 0, 0]);
89        Gs(2) = PP.VapourGibbsFreeEnergy(Outlet.T,Outlet.P,[0, 1, 0, 0]);
90        Gs(3) = PP.VapourGibbsFreeEnergy(Outlet.T,Outlet.P,[0, 0, 1, 0]);
91        Gs(4) = PP.VapourGibbsFreeEnergy(Outlet.T,Outlet.P,[0, 0, 0, 1]);
92       
93       
94        OPTIONS
95        time = [500:10:2000];
96end
97
98
99
100FlowSheet sample_equil_SS
101        PARAMETERS
102        PP                      as CalcObject   (Brief="External physical properties", File="vrpp");
103        NComp           as Integer;
104        NReac           as Integer;
105       
106        stoic(NComp,NReac) as Real      (Brief="Stoichiometric matrix");
107        Rg                      as Real                 (Brief="Universal gas constant", Unit="J/mol/K", Default=8.314);
108        fs                      as pressure     (Brief="Fugacity in standard state", Default=1, Unit="atm");
109       
110       
111        VARIABLES
112        Inlet           as streamTP; # Outlet stream
113        Outlet          as streamTP; # Outlet stream
114
115        Gs(NComp)       as energy_mol   (Brief="Gibbs energy in standard state");
116        K(NReac)        as Real                 (Brief="Equillibrium constant");
117        activ(NComp)as Real             (Brief="Activity");
118       
119        extent(NReac)as flow_mol        (Brief="Extent of reaction");
120       
121       
122        EQUATIONS
123        "Mechanical equilibrium"
124        Outlet.P = Inlet.P;
125
126        "Steady-state"
127        Outlet.F = Inlet.F;
128       
129        "Gibbs energy in standard state"
130        sumt(Gs*stoic) = -Rg*Outlet.T*ln(K);
131       
132        for j in [1:NReac]
133        "Equilibrium constant"
134                K(j) = prod(activ^stoic(:,j));
135        end
136
137        "Activity"
138        activ = PP.VapourFugacityCoefficient(Outlet.T,Outlet.P,Outlet.z)*Outlet.P*Outlet.z/fs;
139
140        for i in [1:NComp]
141        "Outlet molar fraction"
142                (Inlet.F + sum(sumt(stoic*extent)))*Outlet.z(i) = (Inlet.F*Inlet.z(i) + sumt(stoic(i,:)*extent));
143        end
144
145       
146        SET
147        PP.Components = ["carbon monoxide","water","carbon dioxide","hydrogen"];
148        PP.LiquidModel = "IdealLiquid";
149        PP.VapourModel = "Ideal";
150        NComp = PP.NumberOfComponents;
151
152        NReac = 1;
153        stoic(:,1) = [-1.0, -1.0, 1.0, 1.0];
154       
155       
156        SPECIFY
157        Inlet.F = 2*"mol/s";
158        Inlet.P = 1*"bar";
159        Inlet.T = 1100*"K";
160        Inlet.z = [0.5, 0.5, 0.0, 0.0];
161
162        Outlet.T = 1100*"K";
163       
164        Gs(1) = PP.VapourGibbsFreeEnergy(Outlet.T,Outlet.P,[1, 0, 0, 0]);
165        Gs(2) = PP.VapourGibbsFreeEnergy(Outlet.T,Outlet.P,[0, 1, 0, 0]);
166        Gs(3) = PP.VapourGibbsFreeEnergy(Outlet.T,Outlet.P,[0, 0, 1, 0]);
167        Gs(4) = PP.VapourGibbsFreeEnergy(Outlet.T,Outlet.P,[0, 0, 0, 1]);
168       
169       
170        OPTIONS
171        mode = "steady";
172end
173