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

Last change on this file since 992 was 948, checked in by Argimiro Resende Secchi, 10 years ago

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