source: branches/gui/eml/stage_separators/flash.mso @ 823

Last change on this file since 823 was 823, checked in by mamuller, 13 years ago

Added the static head contribution to the OutletLiquid? Pressure for the Dynamic Flash Vessel model.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 13.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* Author: Paula B. Staudt
16* $Id: flash.mso 823 2009-08-11 14:12:13Z mamuller $
17*--------------------------------------------------------------------*#
18
19using "streams";
20
21Model flash
22
23ATTRIBUTES
24        Pallete         = true;
25        Icon            = "icon/Flash";
26        Brief           = "Model of a Dynamic Flash Vessel.";
27        Info            =
28"== ASSUMPTIONS ==
29* perfect mixing of both phases;
30* thermodynamics equilibrium.
31
32== SET ==
33*Orientation: vessel position - vertical or horizontal;
34*Heads (bottom and top heads are identical)
35**elliptical: 2:1 elliptical heads (25% of vessel diameter);
36**hemispherical: hemispherical heads (50% of vessel diameter);
37**flat: flat heads (0% of vessel diameter);
38*Diameter: Vessel diameter;
39*Lenght: Side length of the cylinder shell;
40       
41== SPECIFY ==
42* the Inlet stream;
43* the outlet flows: OutletVapour.F and OutletLiquid.F;
44* the InletQ (the model requires an energy stream, also you can use a controller for setting the heat duty using the heat_flow model).
45
46== OPTIONAL ==
47* the Flash model has three control ports
48** TI OutletLiquid Temperature Indicator;
49** PI OutletLiquid Pressure Indicator;
50** LI Level Indicator;
51
52== INITIAL CONDITIONS ==
53* Initial_Temperature :  the Flash temperature (OutletLiquid.T);
54* Initial_Level : the Flash liquid level (Level);
55* Initial_Composition : (NoComps) OutletLiquid compositions.
56";
57       
58PARAMETERS
59outer PP                as Plugin       (Brief = "External Physical Properties", Type="PP");
60outer NComp     as Integer      (Brief = "Number of components", Lower = 1);
61
62        Mw(NComp)               as molweight    (Brief="Mol Weight", Hidden=true);
63        pi                      as positive             (Brief="Pi value", Default=3.141593,Hidden=true, Symbol="\pi");
64        g                               as acceleration         (Brief="Gravity Acceleration",Default=9.81,Hidden=true);
65       
66        Orientation     as Switcher     (Valid=["vertical","horizontal"],Default="vertical");
67        Heads                   as Switcher     (Valid=["elliptical","hemispherical","flat"],Default="flat");
68        Diameter                as length               (Brief="Vessel diameter", Symbol="D_{i}");
69        Lenght                  as length               (Brief="Side length of the cylinder shell", Symbol="L_{vessel}");
70       
71        Vhead_elliptical                as volume               (Brief="Elliptical Head Total Volume",Hidden=true, Symbol="V_{head}^{elliptical}");
72        Vhead_hemispherical     as volume               (Brief="Hemispherical Head Total Volume",Hidden=true, Symbol="V_{head}^{hemispherical}");
73        Vcylinder                               as volume               (Brief="Cylinder Total Volume",Hidden=true, Symbol="V_{cylinder}");
74        radius                                  as length               (Brief="Vessel radius",Hidden=true, Symbol="R_{cylinder}");
75       
76        Levelpercent_Initial                    as positive     (Brief="Initial liquid height in Percent", Default = 0.70);
77        Temperature_Initial                             as temperature  (Brief="Initial Liquid Temperature", Default = 330);
78        Composition_Initial(NComp)              as fraction             (Brief="Initial Composition", Default = 0.10);
79
80SET
81
82        Mw=PP.MolecularWeight();
83
84        Vhead_elliptical        = (pi*Diameter^3)/12;
85        Vhead_hemispherical = (pi*Diameter^3)/6;
86        Vcylinder = 0.25*(pi*Diameter^2)*Lenght;
87        radius = 0.5*Diameter;
88        g = 9.81 * 'm/(s^2)';
89
90VARIABLES
91
92in      Inlet                   as stream                       (Brief="Feed Stream", PosX=0, PosY=0.48, Symbol="_{in}");
93out     OutletLiquid    as liquid_stream        (Brief="Liquid outlet stream", PosX=0.43, PosY=1, Symbol="_{out}^{Liquid}");
94out     OutletVapour    as vapour_stream        (Brief="Vapour outlet stream", PosX=0.43, PosY=0, Symbol="_{out}^{Vapour}");
95in      InletQ                  as power                        (Brief="Heat Duty", PosX=1, PosY=0.81, Protected =true,Symbol="Q_{in}");
96
97        Vtotal                  as volume                       (Brief="Vessel total volume",Protected=true, Symbol="V_{total}");
98        Vfilled                 as volume                       (Brief="Vessel volume content",Protected=true, Symbol="V_{filled}");
99
100        TotalHoldup(NComp)              as mol  (Brief="Molar Holdup in the Vessel", Protected=true);
101        LiquidHoldup                    as mol  (Brief="Molar liquid holdup", Protected=true);
102        VapourHoldup                    as mol  (Brief="Molar vapour holdup", Protected=true);
103       
104        E                       as energy               (Brief="Total Energy Holdup in the Vessel", Protected=true);
105        vL                      as volume_mol   (Brief="Liquid Molar Volume", Protected=true);
106        vV                      as volume_mol   (Brief="Vapour Molar volume", Protected=true);
107        Level           as length               (Brief="liquid height", Protected=true);
108        Across          as area                 (Brief="Vessel cylinder shell Cross section area", Hidden=true, Symbol="A_{cross}");
109        vfrac           as positive     (Brief="Vapourization fraction", Symbol="\phi", Protected=true);
110        Pratio          as positive             (Brief = "Pressure Ratio", Symbol ="P_{ratio}", Protected=true);       
111        Pdrop           as press_delta  (Brief = "Pressure Drop", DisplayUnit = 'kPa', Symbol ="\Delta P", Protected=true);
112
113        Peq                     as pressure             (Brief="Equilibrium pressure on the liquid surface", Protected=true, DisplayUnit = 'kPa', Symbol="\Delta P_ª{eq}");
114        Pstatic         as pressure             (Brief="Static head at the bottom of the tank", Protected = true, DisplayUnit = 'kPa', Symbol="P_{static}^{Liquid}");
115
116out     TI as control_signal    (Brief="Temperature Indicator", PosX=1, PosY=0.39, Protected=true);
117out     PI as control_signal    (Brief="Pressure Indicator", PosX=1, PosY=0.21, Protected=true);
118out     LI as control_signal    (Brief="Level Indicator", PosX=1, PosY=0.59, Protected=true);
119
120INITIAL
121
122"Initial level Percent"
123        LI = Levelpercent_Initial;
124       
125"Initial Outlet Liquid Temperature"
126        OutletLiquid.T = Temperature_Initial;
127       
128"Initial Outlet Liquid Composition Normalized"
129        OutletLiquid.z(1:NComp - 1) = Composition_Initial(1:NComp - 1)/sum(Composition_Initial);
130
131EQUATIONS
132
133switch Orientation
134
135case "vertical":
136
137"Vessel Cross Section Area"
138        Across = 0.25*(pi*Diameter^2);
139
140switch Heads
141
142case "elliptical":
143
144"Vessel Total Volume"
145        Vtotal = Vhead_elliptical +     Vcylinder;
146
147if Level < 0.25*Diameter then
148
149"Vessel Filled Volume"
150        Vfilled = 0.25*pi*(((Diameter*Level)/(0.25*Diameter))^2)*(0.25*Diameter-Level/3);
151
152else
153
154"Vessel Filled Volume"
155        Vfilled = 0.25*pi*(Diameter^2)*(Level - 0.25*Diameter/3);
156
157end
158
159case "hemispherical":
160
161"Vessel Total Volume"
162        Vtotal = Vhead_hemispherical + Vcylinder;
163
164if Level < 0.5*Diameter then
165
166"Vessel Filled Volume"
167        Vfilled = 0.25*pi*(Level^2)*(2*Diameter-4*Level/3);
168
169else
170
171"Vessel Filled Volume"
172        Vfilled = 0.25*pi*((2/3)*((0.5*Diameter)^3) - (0.25*(Diameter)^3) + Level*Diameter^2);
173
174end
175
176case "flat":
177
178"Vessel Total Volume"
179        Vtotal = Vcylinder;
180
181"Vessel Filled Volume"
182        Vfilled = Across*Level;
183
184end
185
186case "horizontal":
187
188"Vessel Cross Section Area"
189        Across = (radius^2)*acos((radius-Level)/radius)-(radius-Level)*sqrt((2*radius*Level-Level^2));
190
191switch Heads
192
193case "elliptical":
194
195"Vessel Total Volume"
196        Vtotal = Vhead_elliptical +     Vcylinder;
197
198"Vessel Filled Volume"
199        Vfilled = 0.5236*Level^2*(1.5*Diameter-Level) + Across*Lenght;
200
201case "hemispherical":
202
203"Vessel Total Volume"
204        Vtotal = Vhead_hemispherical + Vcylinder;
205
206"Vessel Filled Volume"
207        Vfilled = 1.0472*Level^2*(1.5*Diameter-Level) + Across*Lenght;
208
209case "flat":
210
211"Vessel Total Volume"
212        Vtotal = Vcylinder;
213
214"Vessel Filled Volume"
215        Vfilled = Across*Lenght;
216
217end
218
219end
220
221"Component Molar Balance"
222        diff(TotalHoldup)=Inlet.F*Inlet.z - OutletLiquid.F*OutletLiquid.z - OutletVapour.F*OutletVapour.z;
223       
224"Energy Balance"
225        diff(E) = Inlet.F*Inlet.h - OutletLiquid.F*OutletLiquid.h - OutletVapour.F*OutletVapour.h + InletQ;
226       
227"Molar Holdup"
228        TotalHoldup = LiquidHoldup*OutletLiquid.z + VapourHoldup*OutletVapour.z;
229       
230"Energy Holdup"
231        E = LiquidHoldup*OutletLiquid.h + VapourHoldup*OutletVapour.h - OutletLiquid.P*Vtotal;
232       
233"Mol fraction normalisation"
234        sum(OutletLiquid.z)=1.0;
235
236"Mol fraction normalisation"
237        sum(OutletLiquid.z)=sum(OutletVapour.z);
238
239"Vaporization Fraction"
240        OutletVapour.F = Inlet.F * vfrac;
241
242"Liquid Volume"
243        vL = PP.LiquidVolume(OutletLiquid.T, Peq, OutletLiquid.z);
244
245"Vapour Volume"
246        vV = PP.VapourVolume(OutletVapour.T, Peq, OutletVapour.z);
247       
248"Chemical Equilibrium"
249        PP.LiquidFugacityCoefficient(OutletLiquid.T, Peq, OutletLiquid.z)*OutletLiquid.z =
250                PP.VapourFugacityCoefficient(OutletVapour.T, Peq, OutletVapour.z)*OutletVapour.z;
251       
252"Thermal Equilibrium"
253        OutletVapour.T = OutletLiquid.T;
254
255"Mechanical Equilibrium for the Vapour Phase"
256        OutletVapour.P = Peq;
257       
258"Static Head"   
259        Pstatic = PP.LiquidDensity(OutletLiquid.T, Peq, OutletLiquid.z) * g * Level;
260
261"Mechanical Equilibrium for the Liquid Phase"
262        OutletLiquid.P = Peq + Pstatic;
263
264"Pressure Drop"
265        OutletLiquid.P  = Inlet.P - Pdrop;
266
267"Pressure Ratio"
268        OutletLiquid.P = Inlet.P * Pratio;
269
270"Geometry Constraint"
271        Vtotal = LiquidHoldup * vL + VapourHoldup * vV;
272
273"Temperature indicator"
274        TI * 'K' = OutletLiquid.T;
275
276"Pressure indicator"
277        PI * 'atm' = OutletLiquid.P;
278
279"Level indicator"
280        LI*Vtotal= Vfilled;
281
282"Liquid Level"
283        LiquidHoldup * vL = Vfilled;
284
285end
286
287Model flash_steady
288
289ATTRIBUTES
290        Pallete         = true;
291        Icon            = "icon/Flash";
292        Brief           = "Model of a static PH flash.";
293        Info            =
294"This model is for using the flashPH routine available on VRTherm.
295
296== ASSUMPTIONS ==
297* perfect mixing of both phases;
298* thermodynamics equilibrium.
299* static model.
300
301== SPECIFY ==
302* The Inlet stream;
303* The heat duty;
304* The outlet pressure.
305";     
306
307PARAMETERS
308
309outer PP                as Plugin(Brief = "External Physical Properties", Type="PP");
310outer NComp     as Integer;
311
312VARIABLES
313
314in      Inlet                   as stream                       (Brief="Feed Stream", PosX=0, PosY=0.48, Symbol="_{in}");
315out     OutletLiquid    as liquid_stream        (Brief="Liquid outlet stream", PosX=0.43, PosY=1, Symbol="_{out}^{Liquid}");
316out     OutletVapour    as vapour_stream        (Brief="Vapour outlet stream", PosX=0.43, PosY=0, Symbol="_{out}^{Vapour}");
317in      InletQ                  as power                        (Brief="Heat Duty", PosX=1, PosY=0.81, Protected =true,Symbol="Q_{in}");
318
319        vfrac           as fraction             (Brief="Vaporization fraction", Symbol="\phi", Protected =true);
320        h                       as enth_mol             (Brief="Mixture enthalpy", Hidden =true);
321        Pratio          as positive     (Brief = "Pressure Ratio", Symbol ="P_{ratio}", Protected =true);       
322        Pdrop           as press_delta  (Brief = "Pressure Drop", DisplayUnit = 'kPa', Symbol ="\Delta P", Protected =true);
323
324EQUATIONS
325
326if vfrac > 0 and vfrac <1
327
328then
329"The flash calculation"
330        [vfrac, OutletLiquid.z, OutletVapour.z] = PP.Flash(OutletVapour.T, OutletVapour.P, Inlet.z);
331
332else
333"Chemical equilibrium"
334        [vfrac,OutletLiquid.z,OutletVapour.z]=PP.FlashPH(OutletLiquid.P,h,Inlet.z);
335
336end
337
338"Global Molar Balance"
339        Inlet.F = OutletVapour.F + OutletLiquid.F;
340
341"Vapour Fraction"
342        OutletVapour.F = Inlet.F * vfrac;
343
344"Energy Balance"
345        Inlet.F*(h - Inlet.h) = InletQ;
346        Inlet.F*h = Inlet.F*(1-vfrac)*OutletLiquid.h + Inlet.F*vfrac*OutletVapour.h;
347
348"Thermal Equilibrium"
349        OutletVapour.T = OutletLiquid.T;
350       
351"Mechanical Equilibrium"
352        OutletVapour.P = OutletLiquid.P;
353
354"Pressure Drop"
355        OutletLiquid.P  = Inlet.P - Pdrop;
356
357"Pressure Ratio"
358        OutletLiquid.P = Inlet.P * Pratio;
359
360end
361
362Model FlashPHSteady
363        ATTRIBUTES
364        Pallete         = false;
365        Icon            = "icon/Flash";
366        Brief           = "Another model of a static PH flash.";
367        Info            =
368"This model shows how to model a pressure enthalpy flash
369directly with the EMSO modeling language.
370
371This model is for demonstration purposes only, the flashPH
372routine available on VRTherm is much more robust.
373
374== Assumptions ==
375* perfect mixing of both phases;
376
377== Specify ==
378* the feed stream;
379* the heat duty;
380* the outlet pressure.
381";     
382       
383        PARAMETERS
384outer PP as Plugin(Brief = "External Physical Properties", Type="PP");
385outer NComp as Integer;
386        B as Real(Default=1000, Brief="Regularization Factor");
387
388        VARIABLES
389in      Inlet as stream(Brief="Feed Stream", PosX=0, PosY=0.5421, Symbol="_{in}");
390out     OutletL as liquid_stream(Brief="Liquid outlet stream", PosX=0.4790, PosY=1, Symbol="_{outL}");
391out     OutletV as vapour_stream(Brief="Vapour outlet stream", PosX=0.4877, PosY=0, Symbol="_{outV}");
392in      InletQ as power (Brief="Rate of heat supply", PosX=1, PosY=0.7559, Symbol="_{in}");
393        vfrac as fraction(Brief="Vaporization fraction", Symbol="\phi");
394        vsat as Real(Lower=-0.1, Upper=1.1, Brief="Vaporization fraction if saturated", Symbol="\phi_{sat}");
395        Tsat as temperature(Lower=173, Upper=1473, Brief="Temperature if saturated");
396        xsat(NComp) as Real(Lower=0, Upper=1, Brief="Liquid composition if saturated");
397        ysat(NComp) as Real(Lower=0, Upper=1, Brief="Vapour composition if saturated");
398        Pratio as positive (Brief = "Pressure Ratio", Symbol ="P_{ratio}");     
399        Pdrop as press_delta (Brief = "Pressure Drop", DisplayUnit = 'kPa', Symbol ="\Delta P");
400       
401        zero_one as fraction(Brief="Regularization Variable");
402        one_zero as fraction(Brief="Regularization Variable");
403
404        EQUATIONS
405        "Chemical equilibrium"
406        PP.LiquidFugacityCoefficient(Tsat, OutletL.P, xsat)*xsat =
407                PP.VapourFugacityCoefficient(Tsat, OutletV.P, ysat)*ysat;
408
409        "Global Molar Balance"
410        Inlet.F = OutletV.F + OutletL.F;
411        OutletV.F = Inlet.F * vfrac;
412
413        "Component Molar Balance"
414        Inlet.F*Inlet.z = OutletL.F*xsat + OutletV.F*ysat;
415        sum(xsat) = sum(ysat);
416
417        "Energy Balance if saturated"
418        Inlet.F*Inlet.h  + InletQ =
419                Inlet.F*(1-vsat)*PP.LiquidEnthalpy(Tsat, OutletL.P, xsat) +
420                Inlet.F*vsat*PP.VapourEnthalpy(Tsat, OutletV.P, ysat);
421
422        "Real Energy Balance"
423        Inlet.F*Inlet.h  + InletQ =
424                Inlet.F*(1-vfrac)*OutletL.h + Inlet.F*vfrac*OutletV.h;
425
426        "Thermal Equilibrium"
427        OutletV.T = OutletL.T;
428       
429        "Mechanical Equilibrium"
430        OutletV.P = OutletL.P;
431
432        "Pressure Drop"
433        OutletL.P  = Inlet.P - Pdrop;
434
435        "Pressure Ratio"
436        OutletL.P = Inlet.P * Pratio;
437
438        # regularization functions
439        zero_one = (1 + tanh(B * vsat))/2;
440        one_zero = (1 - tanh(B * (vsat - 1)))/2;
441       
442        vfrac = zero_one * one_zero * vsat + 1 - one_zero;
443        OutletL.z = zero_one*one_zero*xsat + (1-zero_one*one_zero)*Inlet.z;
444        OutletV.z = zero_one*one_zero*ysat + (1-zero_one*one_zero)*Inlet.z;
445end
Note: See TracBrowser for help on using the repository browser.