Changeset 414
 Timestamp:
 Nov 25, 2007, 2:15:43 PM (15 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/eml/reactors/gibbs.mso
r403 r414 17 17 * 18 18 * 19 * 19 * Description: 20 21 * Thermodynamic equilibrium modeling of a reactor using Gibbs 22 * free energy minimization approach. 23 24 * 25 26 * Assumptions: 27 28 * * thermodynamic equilibrium 29 30 * * steadystate 31 32 * 33 34 * Specify: 35 36 * * inlet stream 37 * * number of elements related to components 38 39 * * matrix of elements by components 40 41 * * equilibrium temperature 20 42 * 21 43 * … … 30 52 * only vapour phase 31 53 **# 32 33 54 Model gibbs_vap as tank_vap 34 55 ATTRIBUTES 35 Pallete 36 Icon 37 Brief 38 Info 56 Pallete = true; 57 Icon = "icon/cstr"; 58 Brief = "Model of a generic vapourphase Gibbs CSTR"; 59 Info = " 39 60 Requires the information of: 40 61 * number of elements 41 * matrix of elements (elements , compounds)62 * matrix of elements (elements by compoments) 42 63 "; 43 64 … … 47 68 na(NElem,NComp) as Real (Brief="Number of elements per component"); 48 69 fs(NComp) as pressure (Brief="Fugacity in standard state", Default=1, DisplayUnit='atm'); 49 T 0as temperature (Brief="Reference temperature", Default=298.15);70 To as temperature (Brief="Reference temperature", Default=298.15); 50 71 51 72 VARIABLES 52 73 out Outlet as vapour_stream; # Outlet stream 53 74 54 G s(NComp) as energy_mol (Brief="Gibbs energy in standard state");75 G(NComp) as energy_mol (Brief="Gibbs freeenergy change of formation"); 55 76 lambda(NElem) as energy_mol (Brief="Lagrangian multiplier"); 77 activ(NComp) as Real (Brief="Activity", Lower=1e20); 56 78 phi(NComp) as fugacity (Brief="Fugacity coefficient", Default=1); 57 activ(NComp) as Real (Brief="Activity", Lower=1e20); 58 79 59 80 rate(NComp) as reaction_mol (Brief="Overall component rate of reaction"); 60 81 conv(NComp) as Real (Brief="Fractional conversion of component", Default=0); … … 80 101 sumt(Fi*na) = sumt(Outletm.F*Outletm.z*na); 81 102 82 "Gibbs Energy of Formation"83 G s= PP.IdealGasGibbsOfFormation(Outlet.T);84 85 # "Gibbs Energy of Formation without Cp correction"86 # G s = PP.IdealGasGibbsOfFormationAt25C()*Outlet.T/T0+PP.IdealGasEnthalpyOfFormationAt25C()*(1Outlet.T/T0);103 "Gibbs freeenergy of formation" 104 G = PP.IdealGasGibbsOfFormation(Outlet.T); 105 106 # "Gibbs freeenergy of formation without Cp correction" 107 # G = PP.IdealGasGibbsOfFormationAt25C()*Outlet.T/To+PP.IdealGasEnthalpyOfFormationAt25C()*(1Outlet.T/To); 87 108 88 109 for i in [1:NComp] 89 110 "Lagrangian multiplier" 90 G s(i) + sumt(lambda*na(:,i)) = Rg*Outlet.T*ln(activ(i));111 G(i) + sumt(lambda*na(:,i)) = Rg*Outlet.T*ln(activ(i)); 91 112 92 113 if (Outletm.z(i) > 0) then … … 110 131 activ = phi*Outlet.P*Outlet.z/fs; 111 132 end 133 134 135 #* 136 * only liquid phase 137 **# 138 Model gibbs_liq as tank_liq 139 ATTRIBUTES 140 Pallete = true; 141 Icon = "icon/cstr"; 142 Brief = "Model of a generic liquidphase Gibbs CSTR"; 143 Info = " 144 Requires the information of: 145 * number of elements 146 * matrix of elements (elements by compoments) 147 "; 148 149 PARAMETERS 150 outer 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 157 out Outlet as liquid_stream; # Outlet stream 158 159 G(NComp) as energy_mol (Brief="Gibbs freeenergy change of formation"); 160 lambda(NElem) as energy_mol (Brief="Lagrangian multiplier"); 161 activ(NComp) as Real (Brief="Activity", Lower=0); 162 gamma(NComp) as fugacity (Brief="Activity coefficient", Default=1); 163 164 rate(NComp) as reaction_mol (Brief="Overall component rate of reaction"); 165 conv(NComp) as Real (Brief="Fractional conversion of component", Default=0); 166 Fi(NComp) as flow_mol (Brief="Component molar flow rate"); 167 168 EQUATIONS 169 "Outlet stream" 170 Outlet.F*Outlet.z = Outletm.F*Outletm.z + rate*V; 171 172 "Mechanical equilibrium" 173 Outlet.P = Outletm.P; 174 175 "Steadystate" 176 Outlet.F = sum(Fi); 177 178 "Component molar flow rate" 179 Fi = Outlet.F*Outlet.z; 180 181 "Energy balance" 182 Outlet.F*Outlet.h = Outletm.F*Outletm.h; 183 184 "Element balance" 185 sumt(Fi*na) = sumt(Outletm.F*Outletm.z*na); 186 187 "Gibbs freeenergy of formation" 188 G = PP.IdealGasGibbsOfFormation(Outlet.T); 189 190 # "Gibbs freeenergy of formation without Cp correction" 191 # G = PP.IdealGasGibbsOfFormationAt25C()*Outlet.T/To+PP.IdealGasEnthalpyOfFormationAt25C()*(1Outlet.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) > 0) 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 209 end 210 211 "Activity coefficient" 212 gamma = PP.LiquidFugacityCoefficient(Outlet.T,Outlet.P,Outlet.z); 213 214 "Activity" 215 activ = gamma*Outlet.z*exp(PP.LiquidVolume(Outlet.T,Outlet.P,Outlet.z)*(Outlet.P  Ps)/Rg/Outlet.T); 216 end
Note: See TracChangeset
for help on using the changeset viewer.