source: trunk/eml/stage_separators/tray.mso @ 706

Last change on this file since 706 was 706, checked in by Argimiro Resende Secchi, 14 years ago

Removing temporary variable from tray model.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 24.5 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: Paula B. Staudt
17* $Id: tray.mso 706 2009-01-26 22:08:28Z arge $
18*--------------------------------------------------------------------*#
19
20using "streams";
21
22Model trayBasic
23        ATTRIBUTES
24        Pallete         = false;
25        Icon            = "icon/Tray";
26        Brief           = "Basic equations of a tray column model.";
27        Info            =
28"This model contains only the main equations of a column tray equilibrium model without
29the hidraulic equations.
30       
31== Assumptions ==
32* both phases (liquid and vapour) exists all the time;
33* thermodymanic equilibrium with Murphree plate efficiency;
34* no entrainment of liquid or vapour phase;
35* no weeping;
36* the dymanics in the downcomer are neglected.
37";
38       
39        PARAMETERS
40outer PP as Plugin(Brief = "External Physical Properties", Type="PP");
41outer NComp as Integer;
42        V as volume(Brief="Total Volume of the tray");
43        Q as heat_rate (Brief="Rate of heat supply");
44        Ap as area (Brief="Plate area = Atray - Adowncomer");
45       
46        VARIABLES
47in      Inlet as stream (Brief="Feed stream", PosX=0, PosY=0.4932, Symbol="_{in}");
48in      InletL as stream (Brief="Inlet liquid stream", PosX=0.5195, PosY=0, Symbol="_{inL}");
49in      InletV as stream (Brief="Inlet vapour stream", PosX=0.4994, PosY=1, Symbol="_{inV}");
50out     OutletL as liquid_stream (Brief="Outlet liquid stream", PosX=0.8277, PosY=1, Symbol="_{outL}");
51out     OutletV as vapour_stream (Brief="Outlet vapour stream", PosX=0.8043, PosY=0, Symbol="_{outV}");
52
53        M(NComp) as mol (Brief="Molar Holdup in the tray");
54        ML as mol (Brief="Molar liquid holdup");
55        MV as mol (Brief="Molar vapour holdup");
56        E as energy (Brief="Total Energy Holdup on tray");
57        vL as volume_mol (Brief="Liquid Molar Volume");
58        vV as volume_mol (Brief="Vapour Molar volume");
59        Level as length (Brief="Height of clear liquid on plate");
60        yideal(NComp) as fraction;
61        Emv as Real (Brief = "Murphree efficiency");
62       
63        EQUATIONS
64        "Component Molar Balance"
65        diff(M)=Inlet.F*Inlet.z + InletL.F*InletL.z + InletV.F*InletV.z
66                - OutletL.F*OutletL.z - OutletV.F*OutletV.z;
67       
68        "Energy Balance"
69        diff(E) = ( Inlet.F*Inlet.h + InletL.F*InletL.h + InletV.F*InletV.h
70                - OutletL.F*OutletL.h - OutletV.F*OutletV.h + Q );
71       
72        "Molar Holdup"
73        M = ML*OutletL.z + MV*OutletV.z;
74       
75        "Energy Holdup"
76        E = ML*OutletL.h + MV*OutletV.h - OutletL.P*V;
77       
78        "Mol fraction normalisation"
79        sum(OutletL.z)= 1.0;
80        sum(OutletL.z)= sum(OutletV.z);
81       
82        "Liquid Volume"
83        vL = PP.LiquidVolume(OutletL.T, OutletL.P, OutletL.z);
84        "Vapour Volume"
85        vV = PP.VapourVolume(OutletV.T, OutletV.P, OutletV.z);
86       
87        "Chemical Equilibrium"
88        PP.LiquidFugacityCoefficient(OutletL.T, OutletL.P, OutletL.z)*OutletL.z =
89                PP.VapourFugacityCoefficient(OutletV.T, OutletV.P, yideal)*yideal;
90
91        "Murphree Efficiency"
92        OutletV.z = Emv * (yideal - InletV.z) + InletV.z;
93       
94        "Thermal Equilibrium"
95        OutletV.T = OutletL.T;
96       
97        "Mechanical Equilibrium"
98        OutletV.P = OutletL.P;
99       
100        "Geometry Constraint"
101        V = ML* vL + MV*vV;
102       
103        "Level of clear liquid over the weir"
104        Level = ML*vL/Ap;
105end
106
107Model tray as trayBasic
108        ATTRIBUTES
109        Pallete         = false;
110        Icon            = "icon/Tray";
111        Brief           = "Complete model of a column tray.";
112        Info            =
113"== Specify ==
114* the Feed stream
115* the Liquid inlet stream
116* the Vapour inlet stream
117* the Vapour outlet flow (OutletV.F)
118       
119== Initial ==
120* the plate temperature (OutletL.T)
121* the liquid height (Level) OR the liquid flow OutletL.F
122* (NoComps - 1) OutletL compositions
123
124== Options ==
125You can choose the equation for the liquid outlet flow and the vapour
126inlet flow calculation through the VapourFlowModel and LiquidFlowModel
127switchers.
128
129== References ==
130* ELGUE, S.; PRAT, L.; CABASSUD, M.; LANN, J. L.; CéZERAC, J. Dynamic models for start-up operations of batch distillation columns with experimental validation. Computers and Chemical Engineering, v. 28, p. 2735-2747, 2004.
131* FEEHERY, W. F. Dynamic Optimization with Path Constraints. Tese (Doutorado) - Massachusetts Institute of Technology, June 1998.
132* KLINGBERG, A. Modeling and Optimization of Batch Distillation. Dissertação (Mestrado) - Department of Automatic Control, Lund Institute of Technology, Lund, Sweden, fev. 2000.
133* OLSEN, I.; ENDRESTOL, G. O.; SIRA, T. A rigorous and efficient distillation column model for engineering and training simulators. Computers and Chemical Engineering,v. 21, n. Suppl, p. S193-S198, 1997.
134* REEPMEYER, F.; REPKE, J.-U.; WOZNY, G. Analysis of the start-up process for reactive distillation. Chemical Engineering Technology, v. 26, p. 81-86, 2003.
135* ROFFEL, B.; BETLEM, B.; RUIJTER, J. de. First principles dynamic modeling and multivariable control of a cryogenic distillation column process. Computers and Chemical Engineering, v. 24, p. 111-123, 2000.
136* WANG, L.; LI, P.; WOZNY, G.; WANG, S. A start-up model for simulation of batch distillation starting from a cold state. Computers and Chemical Engineering, v. 27, p.1485-1497, 2003.
137";     
138
139        PARAMETERS
140        Ah as area (Brief="Total holes area");
141        lw as length (Brief="Weir length");
142        g as acceleration (Default=9.81);
143        hw as length (Brief="Weir height");
144        beta as fraction (Brief="Aeration fraction");
145        alfa as fraction (Brief="Dry pressure drop coefficient");
146        w as Real (Brief="Feehery correlation coefficient", Unit='1/m^4', Default=1);
147        btray as Real (Brief="Elgue correlation coefficient", Unit='kg/m/mol^2', Default=1);
148        fw as Real      (Brief = "Olsen correlation coefficient", Default=1);
149        Np as Real      (Brief = "Number of liquid passes in the tray", Default=1);
150        Mw(NComp)       as molweight    (Brief = "Component Mol Weight");
151       
152        VapourFlow as Switcher(Valid = ["on", "off"], Default = "on");
153        LiquidFlow as Switcher(Valid = ["on", "off"], Default = "on");
154        VapourFlowModel as Switcher(Valid = ["Reepmeyer", "Feehery_Fv", "Roffel_Fv", "Klingberg", "Wang_Fv", "Elgue"], Default = "Reepmeyer");
155        LiquidFlowModel as Switcher(Valid = ["default", "Wang_Fl", "Olsen", "Feehery_Fl", "Roffel_Fl"], Default = "default");
156
157        SET
158        Mw = PP.MolecularWeight();
159       
160        VARIABLES
161        rhoL as dens_mass;
162        rhoV as dens_mass;
163       
164        EQUATIONS
165        "Liquid Density"
166        rhoL = PP.LiquidDensity(OutletL.T, OutletL.P, OutletL.z);
167        "Vapour Density"
168        rhoV = PP.VapourDensity(InletV.T, InletV.P, InletV.z);
169
170        switch LiquidFlow
171                case "on":
172                        switch LiquidFlowModel
173                                case "default":
174                                "Francis Equation"
175                                OutletL.F*vL = 1.84*'1/s'*lw*((Level-(beta*hw))/(beta))^2;
176                       
177                                case "Wang_Fl":
178                                OutletL.F*vL = 1.84*'m^0.5/s'*lw*((Level-(beta*hw))/(beta))^1.5;
179                       
180                                case "Olsen":
181                                OutletL.F / 'mol/s'= lw*Np*rhoL/sum(Mw*OutletV.z)/(0.665*fw)^1.5 * ((ML*sum(Mw*OutletL.z)/rhoL/Ap)-hw)^1.5 * 'm^0.5/mol';
182                       
183                                case "Feehery_Fl":
184                                OutletL.F = lw*rhoL/sum(Mw*OutletL.z) * ((Level-hw)/750/'mm')^1.5 * 'm^2/s';
185                       
186                                case "Roffel_Fl":
187                                OutletL.F = 2/3*rhoL/sum(Mw*OutletL.z)*lw*(ML*sum(Mw*OutletL.z)/(Ap*1.3)/rhoL)^1.5*sqrt(2*g/
188                                                        (2*(1 - 0.3593/'Pa^0.0888545'*(OutletV.F*sum(Mw*OutletV.z)/(Ap*1.3)/sqrt(rhoV))^0.177709)-1)); 
189                        end
190                when Level < (beta * hw) switchto "off";
191               
192                case "off":
193                "Low level"
194                OutletL.F = 0 * 'mol/h';
195                when Level > (beta * hw) + 1e-6*'m' switchto "on";
196        end
197       
198        switch VapourFlow
199                case "on":
200                        switch VapourFlowModel
201                                case "Reepmeyer":
202                                InletV.F*vV = sqrt((InletV.P - OutletV.P)/(rhoV*alfa))*Ah;
203                       
204                                case "Feehery_Fv":
205                                InletV.F = rhoV/Ap/w/sum(Mw*OutletV.z) * sqrt(((InletV.P - OutletV.P)-(rhoV*g*ML*vL/Ap))/rhoV);
206                       
207                                case "Roffel_Fv":
208                                InletV.F^1.08 * 0.0013 * 'kg/m/mol^1.08/s^0.92*1e5' = (InletV.P - OutletV.P)*1e5 - (beta*sum(M*Mw)/(Ap*1.3)*g*1e5) * (rhoV*Ah/sum(Mw*OutletV.z))^1.08 * 'm^1.08/mol^1.08';
209                       
210                                case "Klingberg":
211                                InletV.F * vV = Ap * sqrt(((InletV.P - OutletV.P)-rhoL*g*Level)/rhoV);
212                       
213                                case "Wang_Fv":
214                                InletV.F * vV = Ap * sqrt(((InletV.P - OutletV.P)-rhoL*g*Level)/rhoV*alfa);
215                               
216                                case "Elgue":
217                                InletV.F  = sqrt((InletV.P - OutletV.P)/btray);
218                        end
219                when InletV.F < 1e-6 * 'kmol/h' switchto "off";
220               
221                case "off":
222                InletV.F = 0 * 'mol/s';
223                when InletV.P > OutletV.P + Level*g*rhoL + 1e-1 * 'atm' switchto "on";
224        end
225
226end
227
228#*-------------------------------------------------------------------
229* Model of a tray with reaction
230*-------------------------------------------------------------------*#
231Model trayReact
232        ATTRIBUTES
233        Pallete         = false;
234        Icon            = "icon/Tray";
235        Brief           = "Model of a tray with reaction.";
236        Info            =
237"== Assumptions ==
238* both phases (liquid and vapour) exists all the time;
239* thermodymanic equilibrium with Murphree plate efficiency;
240* no entrainment of liquid or vapour phase;
241* no weeping;
242* the dymanics in the downcomer are neglected.
243       
244== Specify ==
245* the Feed stream;
246* the Liquid inlet stream;
247* the Vapour inlet stream;
248* the Vapour outlet flow (OutletV.F);
249* the reaction related variables.
250       
251== Initial ==
252* the plate temperature (OutletL.T)
253* the liquid height (Level) OR the liquid flow OutletL.F
254* (NoComps - 1) OutletL compositions
255";
256
257        PARAMETERS
258        outer PP as Plugin(Type="PP");
259        outer NComp as Integer;
260        V as volume(Brief="Total Volume of the tray");
261        Q as power (Brief="Rate of heat supply");
262        Ap as area (Brief="Plate area = Atray - Adowncomer");
263       
264        Ah as area (Brief="Total holes area");
265        lw as length (Brief="Weir length");
266        g as acceleration (Default=9.81);
267        hw as length (Brief="Weir height");
268        beta as fraction (Brief="Aeration fraction");
269        alfa as fraction (Brief="Dry pressure drop coefficient");
270
271        stoic(NComp) as Real(Brief="Stoichiometric matrix");
272        Hr as energy_mol;
273        Pstartup as pressure;
274       
275        VapourFlow as Switcher(Valid = ["on", "off"], Default = "off");
276        LiquidFlow as Switcher(Valid = ["on", "off"], Default = "off");
277       
278        VARIABLES
279in      Inlet as stream (Brief="Feed stream", PosX=0, PosY=0.4932, Symbol="_{in}");
280in      InletL as stream (Brief="Inlet liquid stream", PosX=0.5195, PosY=0, Symbol="_{inL}");
281in      InletV as stream (Brief="Inlet vapour stream", PosX=0.4994, PosY=1, Symbol="_{inV}");
282out     OutletL as liquid_stream (Brief="Outlet liquid stream", PosX=0.8277, PosY=1, Symbol="_{outL}");
283out     OutletV as vapour_stream (Brief="Outlet vapour stream", PosX=0.8043, PosY=0, Symbol="_{outV}");
284
285        yideal(NComp) as fraction;
286        Emv as Real (Brief = "Murphree efficiency");
287
288        M(NComp) as mol (Brief="Molar Holdup in the tray");
289        ML as mol (Brief="Molar liquid holdup");
290        MV as mol (Brief="Molar vapour holdup");
291        E as energy (Brief="Total Energy Holdup on tray");
292        vL as volume_mol (Brief="Liquid Molar Volume");
293        vV as volume_mol (Brief="Vapour Molar volume");
294        Level as length (Brief="Height of clear liquid on plate");
295        Vol as volume;
296       
297        rhoL as dens_mass;
298        rhoV as dens_mass;
299        r3 as reaction_mol (Brief = "Reaction resulting ethyl acetate", DisplayUnit = 'mol/l/s');
300        C(NComp) as conc_mol (Brief = "Molar concentration", Lower = -1); #, Unit = "mol/l");
301       
302        EQUATIONS
303        "Molar Concentration"
304        OutletL.z = vL * C;
305       
306        "Reaction"
307        r3 = exp(-7150*'K'/OutletL.T)*(4.85e4*C(1)*C(2) - 1.23e4*C(3)*C(4))*'l/mol/s';
308       
309        "Component Molar Balance"
310        diff(M)=Inlet.F*Inlet.z + InletL.F*InletL.z + InletV.F*InletV.z
311                - OutletL.F*OutletL.z - OutletV.F*OutletV.z + stoic*r3*ML*vL;
312       
313        "Energy Balance"
314        diff(E) = ( Inlet.F*Inlet.h + InletL.F*InletL.h + InletV.F*InletV.h
315                - OutletL.F*OutletL.h - OutletV.F*OutletV.h + Q ) + Hr * r3 * vL*ML;
316       
317        "Molar Holdup"
318        M = ML*OutletL.z + MV*OutletV.z;
319       
320        "Energy Holdup"
321        E = ML*OutletL.h + MV*OutletV.h - OutletL.P*V;
322       
323        "Mol fraction normalisation"
324        sum(OutletL.z)= 1.0;
325       
326        "Liquid Volume"
327        vL = PP.LiquidVolume(OutletL.T, OutletL.P, OutletL.z);
328        "Vapour Volume"
329        vV = PP.VapourVolume(OutletV.T, OutletV.P, OutletV.z);
330
331        "Thermal Equilibrium"
332        OutletV.T = OutletL.T;
333       
334        "Mechanical Equilibrium"
335        OutletV.P = OutletL.P;
336       
337        "Level of clear liquid over the weir"
338        Level = ML*vL/Ap;
339
340        Vol = ML*vL;
341       
342        "Liquid Density"
343        rhoL = PP.LiquidDensity(OutletL.T, OutletL.P, OutletL.z);
344        "Vapour Density"
345        rhoV = PP.VapourDensity(InletV.T, InletV.P, InletV.z);
346
347        switch LiquidFlow
348                case "on":
349                "Francis Equation"
350                OutletL.F*vL = 1.84*'1/s'*lw*((Level-(beta*hw)+1e-6*'m')/(beta))^2;
351                when Level < (beta * hw) switchto "off";
352               
353                case "off":
354                "Low level"
355                OutletL.F = 0 * 'mol/h';
356                when Level > (beta * hw) + 1e-6*'m' switchto "on";
357        end
358
359        switch VapourFlow
360                case "on":
361                #InletV.P = OutletV.P + Level*g*rhoL + rhoV*alfa*(InletV.F*vV/Ah)^2;
362                InletV.F*vV = sqrt((InletV.P - OutletV.P - Level*g*rhoL + 1e-8 * 'atm')/(rhoV*alfa))*Ah;
363                when InletV.P < OutletV.P + Level*g*rhoL switchto "off";
364               
365                case "off":
366                InletV.F = 0 * 'mol/s';
367                when InletV.P > OutletV.P + Level*g*rhoL + 3e-2 * 'atm' switchto "on";
368                #when InletV.P > OutletV.P + Level*beta*g*rhoL + 1e-2 * 'atm' switchto "on";
369        end
370
371        "Chemical Equilibrium"
372        PP.LiquidFugacityCoefficient(OutletL.T, OutletL.P, OutletL.z)*OutletL.z =
373                PP.VapourFugacityCoefficient(OutletV.T, OutletV.P, yideal)*yideal;
374       
375        OutletV.z = Emv * (yideal - InletV.z) + InletV.z;
376       
377        sum(OutletL.z)= sum(OutletV.z);
378       
379        "Geometry Constraint"
380        V = ML* vL + MV*vV;
381end
382
383#*-------------------------------------
384* Model of a packed column stage
385-------------------------------------*#
386Model packedStage
387        ATTRIBUTES
388        Pallete         = false;
389        Icon            = "icon/PackedStage";
390        Brief           = "Complete model of a packed column stage.";
391        Info            =
392"== Specify ==
393* the Feed stream
394* the Liquid inlet stream
395* the Vapour inlet stream
396* the stage pressure drop (deltaP)
397       
398== Initial ==
399* the plate temperature (OutletL.T)
400* the liquid molar holdup ML
401* (NoComps - 1) OutletL compositions
402";     
403       
404        PARAMETERS
405outer PP as Plugin(Brief = "External Physical Properties", Type="PP");
406outer NComp as Integer;
407        PPwater as Plugin(Brief="Physical Properties",
408                Type="PP",
409                Components = [ "water" ],
410                LiquidModel = "PR",
411                VapourModel = "PR"
412        );
413
414        V as volume(Brief="Total Volume of the tray");
415        Q as heat_rate (Brief="Rate of heat supply");
416        d as length (Brief="Column diameter"); 
417
418        a as Real (Brief="surface area per packing volume", Unit='m^2/m^3');
419        g as acceleration;
420        e as Real (Brief="Void fraction of packing, m^3/m^3");
421        Cpo as Real (Brief="Constant for resitance equation"); # Billet and Schultes, 1999.
422        Mw(NComp)       as molweight    (Brief = "Component Mol Weight");
423        hs as length (Brief="Height of the packing stage");
424        Qsil as positive (Brief="Resistance coefficient on the liquid load", Default=1);
425
426        VARIABLES
427in      Inlet as stream (Brief="Feed stream", PosX=0, PosY=0.4932, Symbol="_{in}");
428in      InletL as stream (Brief="Inlet liquid stream", PosX=0.5195, PosY=0, Symbol="_{inL}");
429in      InletV as stream (Brief="Inlet vapour stream", PosX=0.4994, PosY=1, Symbol="_{inV}");
430out     OutletL as liquid_stream (Brief="Outlet liquid stream", PosX=0.8277, PosY=1, Symbol="_{outL}");
431out     OutletV as vapour_stream (Brief="Outlet vapour stream", PosX=0.8043, PosY=0, Symbol="_{outV}");
432
433        M(NComp) as mol (Brief="Molar Holdup in the tray", Default=0.01, Lower=0, Upper=100);
434        ML as mol (Brief="Molar liquid holdup", Default=0.01, Lower=0, Upper=100);
435        MV as mol (Brief="Molar vapour holdup", Default=0.01, Lower=0, Upper=100);
436        E as energy (Brief="Total Energy Holdup on tray", Default=-500);
437        vL as volume_mol (Brief="Liquid Molar Volume");
438        vV as volume_mol (Brief="Vapour Molar volume");
439       
440        miL as viscosity (Brief="Liquid dynamic viscosity", DisplayUnit='kg/m/s');
441        miV as viscosity (Brief="Vapor dynamic viscosity", DisplayUnit='kg/m/s');
442        rhoL as dens_mass;
443        rhoV as dens_mass;
444       
445        deltaP as pressure;
446       
447        uL as velocity (Brief="volume flow rate of liquid, m^3/m^2/s", Lower=-10, Upper=100);
448        uV as velocity (Brief="volume flow rate of vapor, m^3/m^2/s", Lower=-10, Upper=100);
449        dp as length (Brief="Particle diameter", Default=1e-3, Lower=0, Upper=10);
450        invK as positive (Brief="Wall factor", Default=1, Upper=10);
451        Rev as Real (Brief="Reynolds number of the vapor stream", Default=4000);
452        Al as area (Brief="Area occupied by the liquid", Default=0.001, Upper=1);
453        hl as positive (Brief="Column holdup", Unit='m^3/m^3', Default=0.01,Upper=10);
454
455        SET
456        Mw = PP.MolecularWeight();
457
458        EQUATIONS
459        "Component Molar Balance"
460        diff(M)=Inlet.F*Inlet.z + InletL.F*InletL.z + InletV.F*InletV.z
461                - OutletL.F*OutletL.z - OutletV.F*OutletV.z;
462
463        "Energy Balance"
464        diff(E) = ( Inlet.F*Inlet.h + InletL.F*InletL.h + InletV.F*InletV.h
465                - OutletL.F*OutletL.h - OutletV.F*OutletV.h + Q );
466       
467        "Molar Holdup"
468        M = ML*OutletL.z + MV*OutletV.z;
469       
470        "Energy Holdup"
471        E = ML*OutletL.h + MV*OutletV.h - OutletL.P*V;
472       
473        "Mol fraction normalisation"
474        sum(OutletL.z)= 1.0;
475       
476        "Liquid Volume"
477        vL = PP.LiquidVolume(OutletL.T, OutletL.P, OutletL.z);
478        "Vapour Volume"
479        vV = PP.VapourVolume(OutletV.T, OutletV.P, OutletV.z);
480       
481        "Chemical Equilibrium"
482        PP.LiquidFugacityCoefficient(OutletL.T, OutletL.P, OutletL.z)*OutletL.z =
483                PP.VapourFugacityCoefficient(OutletV.T, OutletV.P, OutletV.z)*OutletV.z;
484       
485        "Thermal Equilibrium"
486        OutletV.T = OutletL.T;
487       
488        "Mechanical Equilibrium"
489        OutletL.P = OutletV.P;
490       
491        "Geometry Constraint"
492        V*e = ML*vL + MV*vV;
493       
494        "Liquid Density"
495        rhoL = PP.LiquidDensity(OutletL.T, OutletL.P, OutletL.z);
496        "Vapour Density"
497        rhoV = PP.VapourDensity(InletV.T, InletV.P, InletV.z);
498        "Liquid viscosity"
499        miL = PP.LiquidViscosity(OutletL.T, OutletL.P, OutletL.z);
500        "Vapour viscosity"
501        miV = PP.VapourViscosity(InletV.T, InletV.P, InletV.z);
502
503        "Area occupied by the liquid"
504        Al = ML*vL/hs;
505
506        "Volume flow rate of liquid, m^3/m^2/s"
507        uL * Al = OutletL.F * vL;
508        "Volume flow rate of vapor, m^3/m^2/s"
509        uV * ((d^2*3.14159/4)*e - Al) = OutletV.F * vV;
510       
511        "Liquid holdup"
512        hl = ML*vL/V/e;
513       
514        "Particle diameter"
515        dp = 6 * (1-e)/a;
516       
517        "Wall Factor"
518        invK = (1 + (2*dp/(3*d*(1-e))));
519       
520        "Reynolds number of the vapor stream"
521        Rev*invK = dp*uV*rhoV / (miV*(1-e));
522       
523        deltaP = InletV.P - OutletV.P;
524       
525        "Pressure drop and Vapor flow"
526        deltaP/hs  = Qsil*a*uV^2*rhoV*invK / (2*(e-hl)^3);
527
528        "Liquid holdup"
529        hl = (12*miL*a^2*uL/rhoL/g)^1/3;
530end
531
532#*-------------------------------------
533* Nonequilibrium Model
534-------------------------------------*#
535Model interface
536       
537        ATTRIBUTES
538        Pallete         = false;
539        Icon            = "icon/Tray";
540        Brief           = "Descrition of variables of the equilibrium interface.";
541        Info            =
542"This model contains only the variables of the equilibrium interface.";
543
544        PARAMETERS
545outer PP as Plugin(Brief = "External Physical Properties", Type="PP");
546outer NComp as Integer;
547outer NC1 as Integer;
548       
549        VARIABLES
550        NL(NComp) as flow_mol_delta     (Brief = "Stream Molar Rate on Liquid Phase");
551        NV(NComp) as flow_mol_delta     (Brief = "Stream Molar Rate on Vapour Phase");
552        T as temperature                (Brief = "Stream Temperature");
553        P as pressure                   (Brief = "Stream Pressure");
554        x(NComp) as fraction    (Brief = "Stream Molar Fraction on Liquid Phase");
555        y(NComp) as fraction    (Brief = "Stream Molar Fraction on Vapour Phase");
556        a as area                           (Brief = "Interface Area");
557        htL as heat_trans_coeff (Brief = "Heat Transference Coefficient on Liquid Phase");
558        htV as heat_trans_coeff (Brief = "Heat Transference Coefficient on Vapour Phase");     
559        E_liq as heat_rate      (Brief = "Liquid Energy Rate at interface");
560    E_vap as heat_rate      (Brief = "Vapour Energy Rate at interface");       
561        hL as enth_mol          (Brief = "Liquid Molar Enthalpy");
562        hV as enth_mol          (Brief = "Vapour Molar Enthalpy");
563        kL(NC1,NC1) as velocity (Brief = "Mass Transfer Coefficients");
564        kV(NC1,NC1) as velocity (Brief = "Mass Transfer Coefficients");
565       
566        EQUATIONS
567        "Liquid Enthalpy"
568        hL = PP.LiquidEnthalpy(T, P, x);
569       
570        "Vapour Enthalpy"
571        hV = PP.VapourEnthalpy(T, P, y);
572
573end
574
575Model trayRateBasic
576        ATTRIBUTES
577        Pallete         = false;
578        Icon            = "icon/Tray";
579        Brief           = "Basic equations of a tray rate column model.";
580        Info            =
581"This model contains only the main equations of a column tray nonequilibrium model without
582the hidraulic equations.
583       
584== Assumptions ==
585* both phases (liquid and vapour) exists all the time;
586* no entrainment of liquid or vapour phase;
587* no weeping;
588* the dymanics in the downcomer are neglected.
589";
590       
591        PARAMETERS
592outer PP as Plugin(Brief = "External Physical Properties", Type="PP");
593outer NComp as Integer;
594    NC1 as Integer;
595        V as volume(Brief="Total Volume of the tray");
596        Q as heat_rate (Brief="Rate of heat supply");
597        Ap as area (Brief="Plate area = Atray - Adowncomer");
598       
599        VARIABLES
600in      Inlet as stream (Brief="Feed stream", PosX=0, PosY=0.4932, Symbol="_{in}");
601in      InletFV as stream (Brief="Feed stream", PosX=0, PosY=0.4932, Symbol="_{in}");
602in      InletL as stream (Brief="Inlet liquid stream", PosX=0.5195, PosY=0, Symbol="_{inL}");
603in      InletV as stream (Brief="Inlet vapour stream", PosX=0.4994, PosY=1, Symbol="_{inV}");
604out     OutletL as liquid_stream (Brief="Outlet liquid stream", PosX=0.8277, PosY=1, Symbol="_{outL}");
605out     OutletV as vapour_stream (Brief="Outlet vapour stream", PosX=0.8043, PosY=0, Symbol="_{outV}");
606
607        M_liq(NComp) as mol (Brief="Liquid Molar Holdup in the tray");
608        M_vap(NComp) as mol (Brief="Vapour Molar Holdup in the tray");
609        ML as mol (Brief="Molar liquid holdup");
610        MV as mol (Brief="Molar vapour holdup");
611        E_liq as energy (Brief="Total Liquid Energy Holdup on tray");
612        E_vap as energy (Brief="Total Vapour Energy Holdup on tray");
613        vL as volume_mol (Brief="Liquid Molar Volume");
614        vV as volume_mol (Brief="Vapour Molar volume");
615        Level as length (Brief="Height of clear liquid on plate");
616        interf as interface;   
617
618        SET   
619        NC1=NComp-1;
620
621        EQUATIONS
622        "Component Molar Balance"
623        diff(M_liq)=Inlet.F*Inlet.z + InletL.F*InletL.z
624        - OutletL.F*OutletL.z + interf.NL;
625       
626        diff(M_vap)=InletFV.F*InletFV.z + InletV.F*InletV.z
627        - OutletV.F*OutletV.z - interf.NV;
628       
629        "Energy Balance"
630        diff(E_liq) = Inlet.F*Inlet.h + InletL.F*InletL.h
631                - OutletL.F*OutletL.h  + Q + interf.E_liq;
632       
633        diff(E_vap) = InletFV.F*InletFV.h + InletV.F*InletV.h
634                - OutletV.F*OutletV.h  - interf.E_vap;
635       
636        "Molar Holdup"
637        M_liq = ML*OutletL.z;
638       
639        M_vap = MV*OutletV.z;
640       
641        "Energy Holdup"
642        E_liq = ML*(OutletL.h - OutletL.P*vL);
643       
644        E_vap = MV*(OutletV.h - OutletV.P*vV);
645       
646        "Energy Rate through the interface"
647        interf.E_liq = interf.htL*interf.a*(interf.T-OutletL.T)+sum(interf.NL)*interf.hL;       
648       
649        interf.E_vap = interf.htV*interf.a*(OutletV.T-interf.T)+sum(interf.NV)*interf.hV;
650       
651        "Mass Conservation"
652        interf.NL = interf.NV;
653       
654        "Energy Conservation"
655        interf.E_liq = interf.E_vap;
656       
657        "Mol fraction normalisation"
658        sum(OutletL.z)= 1.0;
659        sum(OutletL.z)= sum(OutletV.z);
660        sum(interf.x)=1.0;
661        sum(interf.x)=sum(interf.y);
662       
663        "Liquid Volume"
664        vL = PP.LiquidVolume(OutletL.T, OutletL.P, OutletL.z);
665        "Vapour Volume"
666        vV = PP.VapourVolume(OutletV.T, OutletV.P, OutletV.z);
667       
668        "Chemical Equilibrium"
669        PP.LiquidFugacityCoefficient(interf.T, interf.P, interf.x)*interf.x =
670                PP.VapourFugacityCoefficient(interf.T, interf.P, interf.y)*interf.y;
671
672        "Geometry Constraint"
673        V = ML*vL + MV*vV;
674       
675        "Level of clear liquid over the weir"
676        Level = ML*vL/Ap;
677
678        "Total Mass Transfer Rates"
679        interf.NL(1:NC1)=interf.a*sumt(interf.kL*(interf.x(1:NC1)-OutletL.z(1:NC1)))/vL+
680                OutletL.z(1:NC1)*sum(interf.NL);
681
682#       interf.NL(1:NC1)=0.01*'kmol/s';
683       
684        interf.NV(1:NC1)=interf.a*sumt(interf.kV*(OutletV.z(1:NC1)-interf.y(1:NC1)))/vV+
685                OutletV.z(1:NC1)*sum(interf.NV);
686
687        "Mechanical Equilibrium"
688        OutletV.P = OutletL.P;
689        interf.P=OutletL.P;
690end
691
692Model trayRate as trayRateBasic
693        ATTRIBUTES
694        Pallete         = false;
695        Icon            = "icon/Tray";
696        Brief           = "Complete rate model of a column tray.";
697        Info            =
698"== Specify ==
699* the Feed stream
700* the Liquid inlet stream
701* the Vapour inlet stream
702* the Vapour outlet flow (OutletV.F)
703       
704== Initial ==
705* the plate temperature of both phases (OutletL.T and OutletV.T)
706* the liquid height (Level) OR the liquid flow holdup (ML)
707* the vapor holdup (MV)
708* (NoComps - 1) OutletL compositions
709";
710
711        PARAMETERS
712        Ah as area (Brief="Total holes area");
713        lw as length (Brief="Weir length");
714        g as acceleration (Default=9.81);
715        hw as length (Brief="Weir height");
716        beta as fraction (Brief="Aeration fraction");
717        alfa as fraction (Brief="Dry pressure drop coefficient");
718       
719        VapourFlow as Switcher(Valid = ["on", "off"], Default = "on");
720        LiquidFlow as Switcher(Valid = ["on", "off"], Default = "on");
721       
722        VARIABLES
723        rhoL as dens_mass;
724        rhoV as dens_mass;
725
726        EQUATIONS
727        "Liquid Density"
728        rhoL = PP.LiquidDensity(OutletL.T, OutletL.P, OutletL.z);
729        "Vapour Density"
730        rhoV = PP.VapourDensity(InletV.T, InletV.P, InletV.z);
731
732        switch LiquidFlow
733                case "on":
734                "Francis Equation"
735#               OutletL.F*vL = 1.84*'m^0.5/s'*lw*((Level-(beta*hw))/(beta))^1.5;
736                OutletL.F*vL = 1.84*'1/s'*lw*((Level-(beta*hw))/(beta))^2;
737                when Level < (beta * hw) switchto "off";
738               
739                case "off":
740                "Low level"
741                OutletL.F = 0 * 'mol/h';
742                when Level > (beta * hw) + 1e-6*'m' switchto "on";
743        end
744
745        switch VapourFlow
746                case "on":
747                InletV.F*vV = sqrt((InletV.P - OutletV.P)/(rhoV*alfa))*Ah;
748                when InletV.F < 1e-6 * 'kmol/h' switchto "off";
749               
750                case "off":
751                InletV.F = 0 * 'mol/s';
752                when InletV.P > OutletV.P + Level*g*rhoL + 1e-1 * 'atm' switchto "on";
753        end     
754end
Note: See TracBrowser for help on using the repository browser.