Changeset 813
- Timestamp:
- Aug 5, 2009, 4:00:16 PM (14 years ago)
- Location:
- branches/gui
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/gui/eml/reactors/pfr.mso
r745 r813 35 35 PARAMETERS 36 36 37 outer 37 outer PP as Plugin (Brief = "External Physical Properties", Type="PP"); 38 38 outer NComp as Integer (Brief="Number of components"); 39 40 NReac as Integer (Brief="Number of reactions"); 41 stoic(NComp, NReac) as Real (Brief = "Stoichiometric Matrix"); 42 NDisc as Integer (Brief="Number of points of discretization", Default=10); 43 Mw(NComp) as molweight (Brief="Component Mol Weight"); 44 L as length (Brief="Reactor Length"); 45 Across as area (Brief="Cross section area"); 46 Dpipe as length (Brief="Reactor Inner Diameter"); 47 pi as Real (Brief="pi number",Default=3.141592, Symbol = "\pi"); 48 dx as length (Brief = "Incremental Length"); 49 dv as volume (Brief = "Incremental volume"); 50 Roughness as length (Brief="Reactor Tube Roughness", Default = 4.572E-5, Symbol = "\varepsilon"); 39 NReac as Integer (Brief="Number of reactions"); 40 stoic(NComp, NReac) as Real (Brief = "Stoichiometric Matrix"); 41 NDisc as Integer (Brief="Number of points of discretization", Default=10); 42 Mw(NComp) as molweight (Brief="Component Mol Weight"); 43 L as length (Brief="Reactor Length"); 44 Across as area (Brief="Cross section area"); 51 45 52 46 SET 53 47 54 Mw = PP.MolecularWeight(); 55 Across = 0.25*pi*Dpipe^2; 56 dx = L/NDisc; 57 dv = Across*dx; 48 Mw = PP.MolecularWeight(); 58 49 59 50 VARIABLES 60 51 61 in 52 in Inlet as stream (Brief = "Inlet Stream", PosX=0, PosY=0.5076, Symbol="_{in}"); 62 53 out Outlet as stream (Brief = "Outlet Stream", PosX=1, PosY=0.5236, Symbol="_{out}"); 63 54 64 str(NDisc+1) as streamPH; 65 vel(NDisc+1) as velocity; 66 vol(NDisc+1) as vol_mol; 67 rho(NDisc+1) as dens_mass; 68 rhom(NDisc+1) as dens_mol; 69 dP(NDisc+1) as press_delta (Brief = "Friction Pressure Drop"); 70 Lincr(NDisc+1) as length (Brief = "Length Points", Symbol = "L_{incr}"); 55 str(NDisc+1) as vapour_stream; 56 vol(NDisc+1) as vol_mol; 57 rho(NDisc+1) as dens_mass; 71 58 72 q(NDisc) as heat_rate; 73 C(NComp, NDisc+1) as conc_mol (Brief="Components concentration"); 74 E(NDisc) as energy (Brief="Total Energy Holdup on element"); 75 r(NReac, NDisc) as reaction_mol; 76 Hr(NReac, NDisc) as heat_reaction; 77 Re(NDisc+1) as Real (Brief = "Reynolds Number Profile",Lower = 1E-6); 78 mu(NDisc+1) as viscosity (Brief = "Viscosity Profile" , Symbol = "\mu"); 79 fns(NDisc+1) as fricfactor (Brief = "No Slip Friction Factor"); 59 q(NDisc) as heat_rate; 60 61 M(NComp, NDisc) as mol (Brief = "Molar holdup"); 62 Mt(NDisc) as mol (Brief = "Molar holdup"); 63 C(NComp, NDisc) as conc_mol (Brief="Components concentration", Lower=-1e-6); 64 E(NDisc) as energy (Brief="Total Energy Holdup on element"); 65 r(NReac, NDisc) as reaction_mol; 66 Hr(NReac, NDisc) as heat_reaction; 67 #Hf(NComp, NDisc) as heat_reaction; 80 68 81 69 EQUATIONS … … 94 82 Outlet.h = str(NDisc+1).h; 95 83 Outlet.v = str(NDisc+1).v; 84 85 for z in [1:NDisc] do 86 for c in [1:NComp] do 87 "Component Molar Balance" 88 diff(M(c,z)) = (str(z).F*str(z).z(c) - str(z+1).F*str(z+1).z(c)) 89 + sum(stoic(c,:)*r(:, z)) * Across*L/NDisc; 90 end 91 92 "Energy Balance" 93 diff(E(z)) = str(z).F*str(z).h - str(z+1).F*str(z+1).h + 94 sum(Hr(:,z)*r(:,z)) * Across*L/NDisc - q(z); 96 95 97 "Reactor Initial Length"98 Lincr(1) = 0*'m';96 "Energy Holdup" 97 E(z) = Mt(z)*str(z+1).h - str(z+1).P*Across*L/NDisc; 99 98 100 C(:,1)*vol(1) = Inlet.z; 99 "mass flow is considered constant" 100 str(z+1).F*vol(z+1)*rho(z+1) = str(z).F*vol(z)*rho(z); 101 101 102 for i in [1:NDisc] do 102 "Molar concentration" 103 C(:,z) * Across*L/NDisc = M(:,z); 104 105 "Geometrical constraint" 106 Across*L/NDisc = Mt(z) * vol(z); 103 107 104 for c in [1:NComp] do 108 "Molar fraction" 109 str(z+1).z * Mt(z) = M(:,z); 110 111 #Hf(:,z) = PP.IdealGasEnthalpyOfFormation(str(z+1).T); 112 #Hr(:,z) = -sum(stoic*Hf(:, z)); 113 end 105 114 106 "Component Molar Balance" 107 #diff(M(c,i)) = str(i).F*str(i).z(c) - str(i+1).F*str(i+1).z(c) + sum(stoic(c,:)*r(:, i)) * dv; 108 diff(C(c,i+1)*dv) = (vel(i)*C(c,i) - vel(i+1)*C(c,i+1) )*Across + sum(stoic(c,:)*r(:, i))*dv; 109 115 for z in [1:NDisc+1] do 116 "Specific Volume" 117 vol(z) = PP.VapourVolume(str(z).T, str(z).P, str(z).z); 118 119 "Specific Mass" 120 rho(z) = PP.VapourDensity(str(z).T, str(z).P, str(z).z); 121 end 110 122 end 111 112 "Constraint"113 sum(str(i+1).z ) = 1;114 # Mt(i) = sum(M(:,i));115 116 "Molar concentration"117 #C(:,i) = rhom(i)*str(i).z;118 str(i+1).z = C(:,i+1) * vol(i+1);119 120 #"Molar fraction"121 #str(i+1).z * Mt(i) = M(:,i);122 123 #"Geometrical constraint"124 #Mt(i) = dv * rhom(i);125 126 "Energy Balance"127 diff(E(i)) = str(i).F*str(i).h - str(i+1).F*str(i+1).h + sum(Hr(:,i)*r(:,i)) * dv - q(i);128 129 "Energy Holdup"130 #E(i) = Mt(i)*str(i+1).h - str(i+1).P*dv;131 E(i) = dv/vol(i+1) * str(i+1).h - str(i+1).P*dv;132 133 "Outlet Pressure"134 str(i+1).P =str(1).P - dP(i+1);135 136 "Incremental Length"137 Lincr(i+1) = i*dx;138 139 end140 141 for i in [1:NDisc+1] do142 143 "Incremental Pressure Drop"144 dP(i) = 0.5*fns(i)*Lincr(i)*rho(i)*vel(i)*vel(i)/Dpipe; # simple equation for pressure drop (Darcy Equation)145 146 "Specific Volume"147 vol(i) = PP.VapourVolume(str(i).T, str(i).P, str(i).z);148 149 "Specific Mass"150 rho(i) = PP.VapourDensity(str(i).T, str(i).P, str(i).z);151 152 "Molar Density"153 rhom(i)* vol(i) = 1;154 155 "Viscosity"156 mu(i) = 0.027*'cP';#PP.VapourViscosity(str(i).T,str(i).P,str(i).z); # VRTherm mu=16 cP !!!!!!157 158 "Reynolds Number"159 Re(i)*mu(i) = rho(i)*vel(i)*Dpipe;160 161 "Velocity"162 str(i).F = vel(i)*Across*rhom(i);163 164 if Re(i) > 2300165 166 then167 "Friction Factor for Pressure Drop - Turbulent Flow"168 #1/sqrt(fns(i))= -2*log(abs(Roughness/Dpipe/3.7+2.51/Re(i)/sqrt(fns(i))));169 1/sqrt(fns(i))= -2*log(Roughness/Dpipe/3.7+2.51/Re(i)/sqrt(fns(i)));170 171 else172 "Friction Factor for Pressure Drop - laminar Flow"173 fns(i)*Re(i) = 16;174 175 end176 177 end178 179 end -
branches/gui/sample/reactors/sample_pfr.mso
r745 r813 39 39 PARAMETERS 40 40 41 PP 41 PP as Plugin (Brief="Physical Properties", 42 42 Type = "PP", 43 Components = ["acetone", " acetic anhydride", "methane" ],43 Components = ["acetone", "ketene", "methane" ], 44 44 LiquidModel = "PR", 45 45 VapourModel = "PR" 46 46 ); 47 47 NComp as Integer; 48 SET 49 NComp = PP.NumberOfComponents; 48 50 49 51 DEVICES 52 s1 as source; 53 Reac as pfr; 50 54 51 s1 as simple_source; 52 Reac as pfr; 55 CONNECTIONS 56 s1.Outlet to Reac.Inlet; 57 58 SET 59 Reac.NDisc = 15; 60 Reac.L = 2.28 * 'm'; 61 # Area was chosen to match the original example data 62 Reac.Across = 1.27*'m^3'/Reac.L; 53 63 54 SET 55 NComp = PP.NumberOfComponents; 56 57 s1.ValidPhases = "Vapour-Liquid"; 58 #s1.CompositionBasis = "Molar"; 59 60 Reac.NDisc = 10; 61 Reac.Dpipe = 3 * 'in'; 62 Reac.L = 15 * 'm'; 64 # Only one reaction 63 65 Reac.NReac = 1; 64 Reac.Roughness = 4.572E-5*'m'; 66 # Reaction 1: acetone -> ketene + methane 67 Reac.stoic(:,1) = [-1, 1, 1]; 65 68 66 # Reaction 1: A -> B + C 67 Reac.stoic(:,1) = [-1, 1, 1]; 69 SPECIFY 70 # Feed specification 71 s1.Fw = 8000 * 'kg/h'; 72 s1.T = 1035 * 'K'; 73 s1.P = 1.6 * 'atm'; 74 s1.Composition = [0.98, 0.01, 0.01]; 75 76 # Adiabatic 77 Reac.q = 0 * 'J/s'; 68 78 69 79 EQUATIONS … … 77 87 Reac.Hr(1,z) = -80.77 * 'kJ/mol'; 78 88 89 "Pressure Drop (no pressure drop)" 90 Reac.str(z+1).P = Reac.str(z).P; 91 79 92 end 80 93 81 SPECIFY82 s1.F = 138/1000 * 'kmol/h';83 s1.T = 1035 * 'K';84 s1.P = 1.6 * 'atm';85 s1.MolarComposition = [1, 0, 0];86 87 Reac.q = 0 * 'J/s';88 89 CONNECTIONS90 91 s1.Outlet to Reac.Inlet;92 93 94 INITIAL 94 95 95 96 for z in [2:Reac.NDisc+1] do 97 Reac.str(z).T = Reac.Inlet.T; 98 Reac.str(z).z(1:NComp) = Reac.Inlet.z(1:NComp); 99 end 96 100 97 Reac.str(z).T = Reac.Inlet.T;98 Reac.vel(z) = Reac.vel(1);99 #Reac.str(z).F = Reac.str(1).F;100 Reac.str(z).z(1:NComp-1) = Reac.Inlet.z(1:NComp-1);101 102 end103 104 VARIABLES105 106 SOMA as Real;107 108 EQUATIONS109 110 SOMA = sum(Reac.Outlet.z);111 112 101 OPTIONS 113 # This model can be solved in both steady or dynamic114 102 Dynamic = true; 115 103 TimeStep = 0.05; 116 TimeEnd = 3; 117 DAESolver(File="sundials", RelativeAccuracy=1e-2); 118 NLASolver(RelativeAccuracy=1e-4); 104 TimeEnd = 10; 119 105 #InitialFile = "PFR_AceticAnhydride"; 120 #Dynamic = false; 121 #DAESolver = "dassl"; 106 #Dynamic = false; # Requires a GuessFile from a dynamic simulation 122 107 end
Note: See TracChangeset
for help on using the changeset viewer.