source: branches/gui/eml/reactors/pfr.mso @ 745

Last change on this file since 745 was 745, checked in by gerson bicca, 14 years ago

PFR reactor updates

  • Property svn:keywords set to Id
File size: 4.9 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* Author: Rafael de P. Soares and Paula B. Staudt
17* $Id: pfr.mso 745 2009-03-20 19:59:07Z bicca $
18*--------------------------------------------------------------------*#
19
20using "streams";
21
22Model pfr
23
24ATTRIBUTES
25        Pallete         = true;
26        Brief           = "Model of a Generic PFR with constant mass holdup";
27        Icon            = "icon/pfr";
28        Info            =
29"== Requires the information of ==
30* Reaction values
31* Heat of reaction
32* Pressure profile
33";
34
35PARAMETERS
36
37outer   PP                                                      as Plugin               (Brief = "External Physical Properties", Type="PP");
38outer   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");
51
52SET
53       
54        Mw              = PP.MolecularWeight();
55        Across  = 0.25*pi*Dpipe^2;
56        dx                      = L/NDisc;
57        dv                      = Across*dx;
58
59VARIABLES
60
61in              Inlet   as stream       (Brief = "Inlet Stream", PosX=0, PosY=0.5076, Symbol="_{in}");
62out     Outlet  as stream       (Brief = "Outlet Stream", PosX=1, PosY=0.5236, Symbol="_{out}");
63
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}");
71       
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");
80
81EQUATIONS
82
83"Inlet boundary"
84        str(1).F = Inlet.F;
85        str(1).T = Inlet.T;
86        str(1).P = Inlet.P;
87        str(1).z = Inlet.z;
88
89"Outlet boundary"
90        Outlet.F = str(NDisc+1).F;
91        Outlet.T = str(NDisc+1).T;
92        Outlet.P = str(NDisc+1).P;
93        Outlet.z = str(NDisc+1).z;
94        Outlet.h = str(NDisc+1).h;
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;
101       
102for i in [1:NDisc] do
103
104for c in [1:NComp] do
105
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
110end
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
139end
140
141for 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;
160       
161"Velocity"
162        str(i).F = vel(i)*Across*rhom(i);
163
164if 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;
174
175end
176
177end
178
179end
Note: See TracBrowser for help on using the repository browser.