source: trunk/eml/stage_separators/flash.mso @ 848

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

Adding dew and bubble point models in flash.mso.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 19.0 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 844 2009-09-10 11:47:35Z arge $
17*--------------------------------------------------------------------*#
18
19using "streams";
20
21Model flash
22        ATTRIBUTES
23        Pallete         = true;
24        Icon            = "icon/Flash";
25        Brief           = "Model of a dynamic flash.";
26        Info            =
27"== Assumptions ==
28* both phases are perfectly mixed.
29       
30== Specify ==
31* the feed stream;
32* the outlet flows: OutletV.F and OutletL.F.
33
34== Initial Conditions ==
35* the flash initial temperature (OutletL.T);
36* the flash initial level (Level);
37* (NoComps - 1) OutletL (OR OutletV) compositions.
38";
39       
40        PARAMETERS
41outer PP as Plugin (Brief = "External Physical Properties", Type="PP");
42outer NComp as Integer (Brief = "Number of chemical components", Lower = 1);
43        V as volume (Brief="Total Volume of the flash");
44        Mw(NComp) as molweight;
45        orientation as Switcher (Valid=["vertical","horizontal"],Default="vertical");
46        diameter as length (Brief="Vessel diameter");
47
48        SET
49        Mw=PP.MolecularWeight();
50
51        VARIABLES
52in      Inlet as stream(Brief="Feed Stream", PosX=0, PosY=0.5421, Symbol="_{in}");
53out     OutletL as liquid_stream(Brief="Liquid outlet stream", PosX=0.4790, PosY=1, Symbol="_{outL}");
54out     OutletV as vapour_stream(Brief="Vapour outlet stream", PosX=0.4877, PosY=0, Symbol="_{outV}");
55in      InletQ as energy_stream (Brief="Rate of heat supply", PosX=1, PosY=0.7559, Symbol="_{in}");
56
57        M(NComp) as mol (Brief="Molar Holdup in the tray");
58        ML as mol (Brief="Molar liquid holdup");
59        MV as mol (Brief="Molar vapour holdup");
60        E as energy (Brief="Total Energy Holdup on tray");
61        vL as volume_mol (Brief="Liquid Molar Volume");
62        vV as volume_mol (Brief="Vapour Molar volume");
63        Level as length (Brief="liquid height");
64        Across as area (Brief="Flash Cross section area");
65        vMfrac as positive (Brief="Vapour Molar fraction", Symbol="\ksi");
66        vfrac as positive (Brief="Vapourization fraction", Symbol="\phi");
67        Pratio as positive      (Brief = "Pressure Ratio", Symbol ="P_{ratio}");       
68        Pdrop as press_delta (Brief = "Pressure Drop", DisplayUnit = 'kPa', Symbol ="\Delta P");
69
70        EQUATIONS
71        "Component Molar Balance"
72        diff(M)=Inlet.F*Inlet.z - OutletL.F*OutletL.z - OutletV.F*OutletV.z;
73       
74        "Energy Balance"
75        diff(E) = Inlet.F*Inlet.h - OutletL.F*OutletL.h - OutletV.F*OutletV.h + InletQ.Q;
76       
77        "Molar Holdup"
78        M = ML*OutletL.z + MV*OutletV.z;
79       
80        "Energy Holdup"
81        E = ML*OutletL.h + MV*OutletV.h - OutletL.P*V;
82       
83        "Mol fraction normalisation"
84        sum(OutletL.z)=1.0;
85
86        "Mol fraction normalisation"
87        sum(OutletL.z)=sum(OutletV.z);
88
89        if Inlet.F > 0 then
90                "Vaporization Ratio"
91                OutletV.F = Inlet.F * vfrac;
92        else
93                "Vaporization Ratio"
94                OutletV.F = (OutletV.F + OutletL.F) * vfrac;
95        end
96
97        "Vaporization Fraction"
98        MV = (ML + MV) * vMfrac;
99
100        "Liquid Volume"
101        vL = PP.LiquidVolume(OutletL.T, OutletL.P, OutletL.z);
102
103        "Vapour Volume"
104        vV = PP.VapourVolume(OutletV.T, OutletV.P, OutletV.z);
105       
106        "Chemical Equilibrium"
107        PP.LiquidFugacityCoefficient(OutletL.T, OutletL.P, OutletL.z)*OutletL.z =
108                PP.VapourFugacityCoefficient(OutletV.T, OutletV.P, OutletV.z)*OutletV.z;
109       
110        "Thermal Equilibrium"
111        OutletV.T = OutletL.T;
112       
113        "Mechanical Equilibrium"
114        OutletV.P = OutletL.P;
115
116        "Pressure Drop"
117        OutletL.P  = Inlet.P - Pdrop;
118
119        "Pressure Ratio"
120        OutletL.P = Inlet.P * Pratio;
121
122        "Geometry Constraint"
123        V = ML * vL + MV * vV;
124
125        switch orientation
126        case "vertical":
127        "Cross Section Area"
128                Across = 0.5 * asin(1) * diameter^2;
129       
130        "Liquid Level"
131                ML * vL = Across * Level;
132
133        case "horizontal":
134        "Cylindrical Side Area"
135                Across = 0.25*diameter^2 * (asin(1) - asin((diameter - 2*Level)/diameter)) +
136                                (Level - 0.5*diameter)*sqrt(Level*(diameter - Level));
137
138        "Liquid Level"
139                0.5 * asin(1) * diameter^2 * ML* vL = Across * V;
140        end
141end
142
143#*----------------------------------------------------------------------
144* Model of a  Steady State flash
145*---------------------------------------------------------------------*#
146Model flash_steady
147        ATTRIBUTES
148        Pallete         = true;
149        Icon            = "icon/Flash";
150        Brief           = "Model of a Steady State flash.";
151        Info            =
152"== Assumptions ==
153* both phases are perfectly mixed.
154       
155== Specify ==
156* the feed stream;
157* the outlet pressure (OutletV.P);
158* the outlet temperature OR the heat supplied.
159";
160       
161        PARAMETERS
162outer PP as Plugin(Brief = "External Physical Properties", Type="PP");
163       
164        VARIABLES
165in      Inlet as stream(Brief="Feed Stream", PosX=0, PosY=0.5421, Symbol="_{in}");
166out     OutletL as liquid_stream(Brief="Liquid outlet stream", PosX=0.4790, PosY=1, Symbol="_{outL}");
167out     OutletV as vapour_stream(Brief="Vapour outlet stream", PosX=0.4877, PosY=0, Symbol="_{outV}");
168in      InletQ as energy_stream (Brief="Rate of heat supply", PosX=1, PosY=0.7559, Symbol="_{in}");
169        vfrac as fraction (Brief="Vapourization fraction", Symbol="\phi");
170        Pratio as positive      (Brief = "Pressure Ratio", Symbol ="P_{ratio}");       
171        Pdrop as press_delta (Brief = "Pressure Drop", DisplayUnit = 'kPa', Symbol ="\Delta P");
172
173        EQUATIONS
174        "The flash calculation"
175        [vfrac, OutletL.z, OutletV.z] = PP.Flash(OutletV.T, OutletV.P, Inlet.z);
176       
177        "Global Molar Balance"
178        Inlet.F = OutletV.F + OutletL.F;
179       
180        "Vaporization Fraction"
181        OutletV.F = Inlet.F * vfrac;
182       
183        "Energy Balance"
184        Inlet.F*Inlet.h  + InletQ.Q = OutletL.F*OutletL.h + OutletV.F*OutletV.h;
185       
186        "Thermal Equilibrium"
187        OutletV.T = OutletL.T;
188       
189        "Mechanical Equilibrium"
190        OutletV.P = OutletL.P;
191
192        "Pressure Drop"
193        OutletL.P  = Inlet.P - Pdrop;
194
195        "Pressure Ratio"
196        OutletL.P = Inlet.P * Pratio;
197end
198
199#*----------------------------------------------------------------------
200* Model of a Steady State Bubble flash
201*---------------------------------------------------------------------*#
202Model bubble_steady
203        ATTRIBUTES
204        Pallete         = true;
205        Icon            = "icon/Flash";
206        Brief           = "Model of a Steady State flash.";
207        Info            =
208"== Assumptions ==
209* both phases are perfectly mixed.
210       
211== Specify ==
212* the feed stream;
213* the outlet pressure (OutletV.P);
214* the outlet temperature OR the heat supplied.
215";
216       
217        PARAMETERS
218outer PP as Plugin(Brief = "External Physical Properties", Type="PP");
219       
220        VARIABLES
221in      Inlet as stream(Brief="Feed Stream", PosX=0, PosY=0.5421, Symbol="_{in}");
222out     OutletL as liquid_stream(Brief="Liquid outlet stream", PosX=0.4790, PosY=1, Symbol="_{outL}");
223out     OutletV as vapour_stream(Brief="Vapour outlet stream", PosX=0.4877, PosY=0, Symbol="_{outV}");
224in      InletQ as energy_stream (Brief="Rate of heat supply", PosX=1, PosY=0.7559, Symbol="_{in}");
225        vfrac as fraction (Brief="Vapourization fraction", Symbol="\phi");
226        Pratio as positive      (Brief = "Pressure Ratio", Symbol ="P_{ratio}");       
227        Pdrop as press_delta (Brief = "Pressure Drop", DisplayUnit = 'kPa', Symbol ="\Delta P");
228
229        EQUATIONS
230        "The flash calculation"
231        PP.LiquidFugacityCoefficient(OutletL.T, OutletL.P, OutletL.z)*OutletL.z =
232                PP.VapourFugacityCoefficient(OutletV.T, OutletV.P, OutletV.z)*OutletV.z;
233
234        "Component Molar Balance"
235        Inlet.z = OutletL.z;
236        sum(OutletL.z) = sum(OutletV.z);
237
238        "Global Molar Balance"
239        Inlet.F = OutletV.F + OutletL.F;
240       
241        "Vaporization Fraction"
242        OutletV.F = Inlet.F * vfrac;
243
244        "Energy Balance"
245        Inlet.F*Inlet.h  + InletQ.Q = OutletL.F*OutletL.h + OutletV.F*OutletV.h;
246       
247        "Thermal Equilibrium"
248        OutletV.T = OutletL.T;
249       
250        "Mechanical Equilibrium"
251        OutletV.P = OutletL.P;
252
253        "Pressure Drop"
254        OutletL.P  = Inlet.P - Pdrop;
255
256        "Pressure Ratio"
257        OutletL.P = Inlet.P * Pratio;
258       
259        "Vapor fraction"
260        vfrac = 0;
261end
262
263#*----------------------------------------------------------------------
264* Model of a Steady State Dew flash
265*---------------------------------------------------------------------*#
266Model dew_steady
267        ATTRIBUTES
268        Pallete         = true;
269        Icon            = "icon/Flash";
270        Brief           = "Model of a Steady State flash.";
271        Info            =
272"== Assumptions ==
273* both phases are perfectly mixed.
274       
275== Specify ==
276* the feed stream;
277* the outlet pressure (OutletV.P);
278* the outlet temperature OR the heat supplied.
279";
280       
281        PARAMETERS
282outer PP as Plugin(Brief = "External Physical Properties", Type="PP");
283       
284        VARIABLES
285in      Inlet as stream(Brief="Feed Stream", PosX=0, PosY=0.5421, Symbol="_{in}");
286out     OutletL as liquid_stream(Brief="Liquid outlet stream", PosX=0.4790, PosY=1, Symbol="_{outL}");
287out     OutletV as vapour_stream(Brief="Vapour outlet stream", PosX=0.4877, PosY=0, Symbol="_{outV}");
288in      InletQ as energy_stream (Brief="Rate of heat supply", PosX=1, PosY=0.7559, Symbol="_{in}");
289        vfrac as fraction (Brief="Vapourization fraction", Symbol="\phi");
290        Pratio as positive      (Brief = "Pressure Ratio", Symbol ="P_{ratio}");       
291        Pdrop as press_delta (Brief = "Pressure Drop", DisplayUnit = 'kPa', Symbol ="\Delta P");
292
293        EQUATIONS
294        "The flash calculation"
295        PP.LiquidFugacityCoefficient(OutletL.T, OutletL.P, OutletL.z)*OutletL.z =
296                PP.VapourFugacityCoefficient(OutletV.T, OutletV.P, OutletV.z)*OutletV.z;
297
298        "Component Molar Balance"
299        Inlet.z = OutletV.z;
300        sum(OutletL.z) = sum(OutletV.z);
301
302        "Global Molar Balance"
303        Inlet.F = OutletV.F + OutletL.F;
304       
305        "Vaporization Fraction"
306        OutletV.F = Inlet.F * vfrac;
307
308        "Energy Balance"
309        Inlet.F*Inlet.h  + InletQ.Q = OutletL.F*OutletL.h + OutletV.F*OutletV.h;
310       
311        "Thermal Equilibrium"
312        OutletV.T = OutletL.T;
313       
314        "Mechanical Equilibrium"
315        OutletV.P = OutletL.P;
316
317        "Pressure Drop"
318        OutletL.P  = Inlet.P - Pdrop;
319
320        "Pressure Ratio"
321        OutletL.P = Inlet.P * Pratio;
322       
323        "Vapor fraction"
324        vfrac = 1;
325end
326
327#*----------------------------------------------------------------------
328* Model of a Steady State flash
329*---------------------------------------------------------------------*#
330Model flash_steady_full
331        ATTRIBUTES
332        Pallete         = true;
333        Icon            = "icon/Flash";
334        Brief           = "Model of a Steady State flash.";
335        Info            =
336"== Assumptions ==
337* both phases are perfectly mixed.
338       
339== Specify ==
340* the feed stream;
341* the outlet pressure (OutletV.P);
342* the outlet temperature OR the heat supplied.
343";
344       
345        PARAMETERS
346outer PP as Plugin(Brief = "External Physical Properties", Type="PP");
347       
348        VARIABLES
349in      Inlet as stream(Brief="Feed Stream", PosX=0, PosY=0.5421, Symbol="_{in}");
350out     OutletL as liquid_stream(Brief="Liquid outlet stream", PosX=0.4790, PosY=1, Symbol="_{outL}");
351out     OutletV as vapour_stream(Brief="Vapour outlet stream", PosX=0.4877, PosY=0, Symbol="_{outV}");
352in      InletQ as energy_stream (Brief="Rate of heat supply", PosX=1, PosY=0.7559, Symbol="_{in}");
353        vfrac as fraction (Brief="Vapourization fraction", Symbol="\phi");
354        Pratio as positive      (Brief = "Pressure Ratio", Symbol ="P_{ratio}");       
355        Pdrop as press_delta (Brief = "Pressure Drop", DisplayUnit = 'kPa', Symbol ="\Delta P");
356
357        EQUATIONS
358        "The flash calculation"
359        PP.LiquidFugacityCoefficient(OutletL.T, OutletL.P, OutletL.z)*OutletL.z =
360                PP.VapourFugacityCoefficient(OutletV.T, OutletV.P, OutletV.z)*OutletV.z;
361
362        "Component Molar Balance"
363        Inlet.F*Inlet.z = OutletL.F*OutletL.z + OutletV.F*OutletV.z;
364        sum(OutletL.z) = sum(OutletV.z);
365
366        "Global Molar Balance"
367        Inlet.F = OutletV.F + OutletL.F;
368       
369        "Vaporization Fraction"
370        OutletV.F = Inlet.F * vfrac;
371
372        "Energy Balance"
373        Inlet.F*Inlet.h  + InletQ.Q = OutletL.F*OutletL.h + OutletV.F*OutletV.h;
374       
375        "Thermal Equilibrium"
376        OutletV.T = OutletL.T;
377       
378        "Mechanical Equilibrium"
379        OutletV.P = OutletL.P;
380
381        "Pressure Drop"
382        OutletL.P  = Inlet.P - Pdrop;
383
384        "Pressure Ratio"
385        OutletL.P = Inlet.P * Pratio;
386end
387
388#*----------------------------------------------------------------------
389* Model of a Steady State flash
390*---------------------------------------------------------------------*#
391Model flash_steady_bd
392        ATTRIBUTES
393        Pallete         = true;
394        Icon            = "icon/Flash";
395        Brief           = "Model of a Steady State flash.";
396        Info            =
397"== Assumptions ==
398* both phases are perfectly mixed.
399       
400== Specify ==
401* the feed stream;
402* the outlet pressure (OutletV.P);
403* the outlet temperature OR the heat supplied.
404";
405       
406        PARAMETERS
407outer PP as Plugin(Brief = "External Physical Properties", Type="PP");
408outer NComp as Integer;
409
410        VARIABLES
411in      Inlet as stream(Brief="Feed Stream", PosX=0, PosY=0.5421, Symbol="_{in}");
412out     OutletL as liquid_stream(Brief="Liquid outlet stream", PosX=0.4790, PosY=1, Symbol="_{outL}");
413out     OutletV as vapour_stream(Brief="Vapour outlet stream", PosX=0.4877, PosY=0, Symbol="_{outV}");
414in      InletQ as energy_stream (Brief="Rate of heat supply", PosX=1, PosY=0.7559, Symbol="_{in}");
415        vfrac as fraction (Brief="Vapourization fraction", Symbol="\phi");
416        Pratio as positive      (Brief = "Pressure Ratio", Symbol ="P_{ratio}");       
417        Pdrop as press_delta (Brief = "Pressure Drop", DisplayUnit = 'kPa', Symbol ="\Delta P");
418        T_dew as temperature (Brief="Dew-point Temperature");
419        T_bubble as temperature (Brief="Bubble-point Temperature");
420        x_dew(NComp) as fraction (Brief="Dew-point liquid composition");
421        y_bubble(NComp) as fraction (Brief="Bubble-point Vapour composition");
422
423        EQUATIONS
424       
425        if OutletL.T > T_bubble and OutletL.T < T_dew then
426           "The flash calculation"
427                PP.LiquidFugacityCoefficient(OutletL.T, OutletL.P, OutletL.z)*OutletL.z =
428                PP.VapourFugacityCoefficient(OutletV.T, OutletV.P, OutletV.z)*OutletV.z;
429
430           "Composition constraint"
431                sum(OutletL.z) = sum(OutletV.z);
432
433           "Component Molar Balance"
434                Inlet.F*Inlet.z = OutletL.F*OutletL.z + OutletV.F*OutletV.z;
435
436        else if OutletL.T <= T_bubble then
437                        "Bubble-point result"
438                        OutletL.z = Inlet.z;
439                        OutletV.z = y_bubble;
440                        vfrac = 0;
441                 else
442                        "Dew-point result"
443                        OutletL.z = x_dew;
444                        OutletV.z = Inlet.z;
445                        vfrac = 1;
446                 end
447        end
448
449        "Dew-point equations"
450        PP.LiquidFugacityCoefficient(T_dew, OutletL.P, x_dew)*x_dew =
451                PP.VapourFugacityCoefficient(T_dew, OutletV.P, Inlet.z)*Inlet.z;
452        sum(x_dew) = 1;
453
454        "Bubble-point equations"
455        PP.LiquidFugacityCoefficient(T_bubble, OutletL.P, Inlet.z)*Inlet.z =
456                PP.VapourFugacityCoefficient(T_bubble, OutletV.P, y_bubble)*y_bubble;
457        sum(y_bubble) = 1;
458
459        "Global Molar Balance"
460        Inlet.F = OutletV.F + OutletL.F;
461       
462        "Vaporization Fraction"
463        OutletV.F = Inlet.F * vfrac;
464
465        "Energy Balance"
466        Inlet.F*Inlet.h  + InletQ.Q = OutletL.F*OutletL.h + OutletV.F*OutletV.h;
467       
468        "Thermal Equilibrium"
469        OutletV.T = OutletL.T;
470       
471        "Mechanical Equilibrium"
472        OutletV.P = OutletL.P;
473
474        "Pressure Drop"
475        OutletL.P  = Inlet.P - Pdrop;
476
477        "Pressure Ratio"
478        OutletL.P = Inlet.P * Pratio;
479end
480
481
482#*----------------------------------------------------------------------
483* Model of a steady-state PH flash.
484*---------------------------------------------------------------------*#
485Model FlashPHSteady
486        ATTRIBUTES
487        Pallete         = true;
488        Icon            = "icon/Flash";
489        Brief           = "Model of a static PH flash.";
490        Info            =
491"This model is for using the flashPH routine available on VRTherm.
492
493== Assumptions ==
494* perfect mixing of both phases;
495
496== Specify ==
497* the feed stream;
498* the heat duty;
499* the outlet pressure.
500";     
501
502        PARAMETERS
503outer PP as Plugin(Brief = "External Physical Properties", Type="PP");
504outer NComp as Integer;
505
506        VARIABLES
507in      Inlet as stream(Brief="Feed Stream", PosX=0, PosY=0.5421, Symbol="_{in}");
508out     OutletL as liquid_stream(Brief="Liquid outlet stream", PosX=0.4790, PosY=1, Symbol="_{outL}");
509out     OutletV as vapour_stream(Brief="Vapour outlet stream", PosX=0.4877, PosY=0, Symbol="_{outV}");
510in      InletQ as energy_stream (Brief="Rate of heat supply", PosX=1, PosY=0.7559, Symbol="_{in}");
511        vfrac as fraction(Brief="Vaporization fraction", Symbol="\phi");
512        h as enth_mol(Brief="Mixture enthalpy");
513        Pratio as positive (Brief = "Pressure Ratio", Symbol ="P_{ratio}");     
514        Pdrop as press_delta (Brief = "Pressure Drop", DisplayUnit = 'kPa', Symbol ="\Delta P");
515
516        EQUATIONS
517
518        "Chemical equilibrium"
519        [vfrac,OutletL.z,OutletV.z]=PP.FlashPH(OutletL.P,h,Inlet.z);
520
521        "Global Molar Balance"
522        Inlet.F = OutletV.F + OutletL.F;
523        OutletV.F = Inlet.F * vfrac;
524
525        "Energy Balance"
526        Inlet.F*(h - Inlet.h) = InletQ.Q;
527        Inlet.F*h = Inlet.F*(1-vfrac)*OutletL.h + Inlet.F*vfrac*OutletV.h;
528
529        "Thermal Equilibrium"
530        OutletV.T = OutletL.T;
531       
532        "Mechanical Equilibrium"
533        OutletV.P = OutletL.P;
534
535        "Pressure Drop"
536        OutletL.P  = Inlet.P - Pdrop;
537
538        "Pressure Ratio"
539        OutletL.P = Inlet.P * Pratio;
540end
541
542#*----------------------------------------------------------------------
543* Another model of a steady-state PH flash.
544* It is recommended to use [v,x,y]=PP.FlashPH(P,h,z) instead of.
545*---------------------------------------------------------------------*#
546Model FlashPHSteadyA
547        ATTRIBUTES
548        Pallete         = true;
549        Icon            = "icon/Flash";
550        Brief           = "Another model of a static PH flash.";
551        Info            =
552"This model shows how to model a pressure enthalpy flash
553directly with the EMSO modeling language.
554
555This model is for demonstration purposes only, the flashPH
556routine available on VRTherm is much more robust.
557
558== Assumptions ==
559* perfect mixing of both phases;
560
561== Specify ==
562* the feed stream;
563* the heat duty;
564* the outlet pressure.
565";     
566       
567        PARAMETERS
568outer PP as Plugin(Brief = "External Physical Properties", Type="PP");
569outer NComp as Integer;
570        B as Real(Default=1000, Brief="Regularization Factor");
571
572        VARIABLES
573in      Inlet as stream(Brief="Feed Stream", PosX=0, PosY=0.5421, Symbol="_{in}");
574out     OutletL as liquid_stream(Brief="Liquid outlet stream", PosX=0.4790, PosY=1, Symbol="_{outL}");
575out     OutletV as vapour_stream(Brief="Vapour outlet stream", PosX=0.4877, PosY=0, Symbol="_{outV}");
576in      InletQ as energy_stream (Brief="Rate of heat supply", PosX=1, PosY=0.7559, Symbol="_{in}");
577        vfrac as fraction(Brief="Vaporization fraction", Symbol="\phi");
578        vsat as Real(Lower=-0.1, Upper=1.1, Brief="Vaporization fraction if saturated", Symbol="\phi_{sat}");
579        Tsat as temperature(Lower=173, Upper=1473, Brief="Temperature if saturated");
580        xsat(NComp) as Real(Lower=0, Upper=1, Brief="Liquid composition if saturated");
581        ysat(NComp) as Real(Lower=0, Upper=1, Brief="Vapour composition if saturated");
582        Pratio as positive (Brief = "Pressure Ratio", Symbol ="P_{ratio}");     
583        Pdrop as press_delta (Brief = "Pressure Drop", DisplayUnit = 'kPa', Symbol ="\Delta P");
584       
585        zero_one as fraction(Brief="Regularization Variable");
586        one_zero as fraction(Brief="Regularization Variable");
587
588        EQUATIONS
589        "Chemical equilibrium"
590        PP.LiquidFugacityCoefficient(Tsat, OutletL.P, xsat)*xsat =
591                PP.VapourFugacityCoefficient(Tsat, OutletV.P, ysat)*ysat;
592
593        "Global Molar Balance"
594        Inlet.F = OutletV.F + OutletL.F;
595        OutletV.F = Inlet.F * vfrac;
596
597        "Component Molar Balance"
598        Inlet.F*Inlet.z = OutletL.F*xsat + OutletV.F*ysat;
599        sum(xsat) = sum(ysat);
600
601        "Energy Balance if saturated"
602        Inlet.F*Inlet.h  + InletQ.Q =
603                Inlet.F*(1-vsat)*PP.LiquidEnthalpy(Tsat, OutletL.P, xsat) +
604                Inlet.F*vsat*PP.VapourEnthalpy(Tsat, OutletV.P, ysat);
605
606        "Real Energy Balance"
607        Inlet.F*Inlet.h  + InletQ.Q =
608                Inlet.F*(1-vfrac)*OutletL.h + Inlet.F*vfrac*OutletV.h;
609
610        "Thermal Equilibrium"
611        OutletV.T = OutletL.T;
612       
613        "Mechanical Equilibrium"
614        OutletV.P = OutletL.P;
615
616        "Pressure Drop"
617        OutletL.P  = Inlet.P - Pdrop;
618
619        "Pressure Ratio"
620        OutletL.P = Inlet.P * Pratio;
621
622        # regularization functions
623        zero_one = (1 + tanh(B * vsat))/2;
624        one_zero = (1 - tanh(B * (vsat - 1)))/2;
625       
626        vfrac = zero_one * one_zero * vsat + 1 - one_zero;
627        OutletL.z = zero_one*one_zero*xsat + (1-zero_one*one_zero)*Inlet.z;
628        OutletV.z = zero_one*one_zero*ysat + (1-zero_one*one_zero)*Inlet.z;
629end
Note: See TracBrowser for help on using the repository browser.