 Timestamp:
 Mar 20, 2009, 4:59:07 PM (14 years ago)
 Location:
 branches/gui
 Files:

 1 added
 2 edited
Legend:
 Unmodified
 Added
 Removed

branches/gui/eml/reactors/pfr.mso
r574 r745 35 35 PARAMETERS 36 36 37 outer PP as Plugin (Brief = "External Physical Properties", Type="PP");37 outer PP as Plugin (Brief = "External Physical Properties", Type="PP"); 38 38 outer NComp as Integer (Brief="Number of components"); 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"); 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.572E5, Symbol = "\varepsilon"); 45 51 46 52 SET 47 53 48 Mw = PP.MolecularWeight(); 54 Mw = PP.MolecularWeight(); 55 Across = 0.25*pi*Dpipe^2; 56 dx = L/NDisc; 57 dv = Across*dx; 49 58 50 59 VARIABLES 51 60 52 in Inlet as stream (Brief = "Inlet Stream", PosX=0, PosY=0.5076, Symbol="_{in}");61 in Inlet as stream (Brief = "Inlet Stream", PosX=0, PosY=0.5076, Symbol="_{in}"); 53 62 out Outlet as stream (Brief = "Outlet Stream", PosX=1, PosY=0.5236, Symbol="_{out}"); 54 63 55 str(NDisc+1) as streamPH; 56 vel(NDisc+1) as velocity; 57 vol(NDisc+1) as vol_mol; 58 rho(NDisc+1) as dens_mass; 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}"); 59 71 60 q(NDisc) as heat_rate; 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=1e6); 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; 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 = 1E6); 78 mu(NDisc+1) as viscosity (Brief = "Viscosity Profile" , Symbol = "\mu"); 79 fns(NDisc+1) as fricfactor (Brief = "No Slip Friction Factor"); 67 80 68 81 EQUATIONS … … 81 94 Outlet.h = str(NDisc+1).h; 82 95 Outlet.v = str(NDisc+1).v; 96 97 "Reactor Initial Length" 98 Lincr(1) = 0*'m'; 99 100 C(:,1)*vol(1) = Inlet.z; 83 101 84 for z in [1:NDisc] do 85 for c in [1:NComp1] do 86 "Component Molar Balance" 87 diff(M(c,z)) = (str(z).F*str(z).z(c)  str(z+1).F*str(z+1).z(c)) 88 + sum(stoic(c,:)*r(:, z)) * Across*L/NDisc; 89 end 102 for i in [1:NDisc] do 90 103 91 "Energy Balance" 92 diff(E(z)) = str(z).F*str(z).h  str(z+1).F*str(z+1).h + 93 sum(Hr(:,z)*r(:,z)) * Across*L/NDisc  q(z); 104 for c in [1:NComp] do 94 105 95 "Energy Holdup" 96 E(z) = Mt(z)*str(z+1).h  str(z+1).P*Across*L/NDisc; 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; 97 109 98 "mass flow is considered constant" 99 str(z+1).F*vol(z+1) = str(z).F*vol(z); # FIXME: is this correct? No (constant velocity: only for equimolar) 100 #rho(z+1)*vel(z+1) = rho(z)*vel(z); # FIXME: this is correct! But does not converge. 110 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 end 140 141 for i in [1:NDisc+1] do 142 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; 101 160 102 "Molar concentration" 103 C(:,z) * Across*L/NDisc = M(:,z); 104 105 "Sum of M" 106 Mt(z) = sum(M(:,z)); 107 108 "Geometrical constraint" 109 Across*L/NDisc = Mt(z) * vol(z); 161 "Velocity" 162 str(i).F = vel(i)*Across*rhom(i); 110 163 111 "Molar fraction" 112 str(z+1).z * Mt(z) = M(:,z); 113 end 164 if Re(i) > 2300 165 166 then 167 "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 else 172 "Friction Factor for Pressure Drop  laminar Flow" 173 fns(i)*Re(i) = 16; 114 174 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); 175 end 121 176 122 "Velocity"123 vel(z)*Across = str(z).F*vol(z);124 end125 177 end 178 179 end 
branches/gui/sample/reactors/sample_pfr.mso
r585 r745 49 49 DEVICES 50 50 51 s1 as s ource;51 s1 as simple_source; 52 52 Reac as pfr; 53 53 … … 56 56 57 57 s1.ValidPhases = "VapourLiquid"; 58 s1.CompositionBasis = "Molar";58 #s1.CompositionBasis = "Molar"; 59 59 60 60 Reac.NDisc = 10; 61 Reac. Across = 0.7 * 'in^2';62 Reac.L = 2.28* 'm';61 Reac.Dpipe = 3 * 'in'; 62 Reac.L = 15 * 'm'; 63 63 Reac.NReac = 1; 64 Reac.Roughness = 4.572E5*'m'; 64 65 65 66 # Reaction 1: A > B + C … … 76 77 Reac.Hr(1,z) = 80.77 * 'kJ/mol'; 77 78 78 "Pressure Drop (no pressure drop)"79 Reac.str(z+1).P = Reac.str(z).P;80 81 79 end 82 80 … … 85 83 s1.T = 1035 * 'K'; 86 84 s1.P = 1.6 * 'atm'; 87 s1. Composition = [1, 0, 0];85 s1.MolarComposition = [1, 0, 0]; 88 86 89 87 Reac.q = 0 * 'J/s'; … … 98 96 99 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 100 Reac.str(z).z(1:NComp1) = Reac.Inlet.z(1:NComp1); 101 101 102 102 end 103 103 104 VARIABLES 105 106 SOMA as Real; 107 108 EQUATIONS 109 110 SOMA = sum(Reac.Outlet.z); 111 104 112 OPTIONS 105 113 # This model can be solved in both steady or dynamic … … 107 115 TimeStep = 0.05; 108 116 TimeEnd = 3; 117 DAESolver(File="sundials", RelativeAccuracy=1e2); 118 NLASolver(RelativeAccuracy=1e4); 119 #InitialFile = "PFR_AceticAnhydride"; 109 120 #Dynamic = false; 110 121 #DAESolver = "dassl";
Note: See TracChangeset
for help on using the changeset viewer.