/[ascend]/trunk/models/ben/benHGthermo.a4l
ViewVC logotype

Contents of /trunk/models/ben/benHGthermo.a4l

Parent Directory Parent Directory | Revision Log Revision Log


Revision 576 - (show annotations) (download) (as text)
Tue May 9 03:42:08 2006 UTC (18 years, 1 month ago) by johnpye
File MIME type: text/x-ascend
File size: 36883 byte(s)
Changed all cases of *.fixed := {TRUE,FALSE} to 'FIX' and 'FREE' statements.
1 REQUIRE "ben/benpropertyoptions.a4l";
2 (* --> measures,system,atoms,components *)
3 PROVIDE "benHGthermo.a4l";
4 (*********************************************************************\
5 H_G_thermodynamics.lib
6 by Joseph J. Zaher
7 Part of the Ascend Library
8
9 This file is part of the Ascend modeling library.
10
11 Copyright (C) 1994 -1997 Carnegie Mellon University
12
13 The Ascend modeling library is free software; you can redistribute
14 it and/or modify it under the terms of the GNU General Public License as
15 published by the Free Software Foundation; either version 2 of the
16 License, or (at your option) any later version.
17
18 The Ascend Language Interpreter is distributed in hope that it will be
19 useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
20 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
21 General Public License for more details.
22
23 You should have received a copy of the GNU General Public License along with
24 the program; if not, write to the Free Software Foundation, Inc., 675
25 Mass Ave, Cambridge, MA 02139 USA. Check the file named COPYING.
26
27 \*********************************************************************)
28
29 (**
30 **
31 ** T H E R M O D Y N A M I C S . L I B
32 ** ----------------------------------------------------
33 **
34 ** AUTHOR: Joseph J. Zaher
35 **
36 ** DATES: 07/91 - Original code.
37 ** 02/92 - Made compatible with new set version of ASCEND.
38 ** Scaled equations to enhance convergence, updated
39 ** vapor pressure correlation, added Pitzer extension
40 ** to vapor mixtures and UNIFAC extension to liquid
41 ** mixtures with help of Bob Huss.
42 ** 03/92 - Removed stream model. Library remains purely
43 ** intensive without any assumption to static or
44 ** dynamic modeling.
45 ** 07/92 - Structural changes to provide a common thermodynamic
46 ** properties root model as the library interface.
47 ** Modified the existing phase distribution model
48 ** to incorporate an intensive mass balance over the
49 ** phases. Residual quantities for pure vapor
50 ** components estimate corrections from ideal gas
51 ** behavior while residual quantities for pure liquid
52 ** components estimate corrections from incompressible
53 ** fluid behavior.
54 ** 08/92 - Allowed component names in mixtures to be freely
55 ** specified by user.
56 ** 03/94 - Made compatible with gnu-ascend.
57 ** 05/94 - Removed refinement link of models correction and
58 ** and partial_component which should not contain T,
59 ** P, and R anyway. The interface to the library
60 ** is now returned to model thermodynamic_properties
61 ** where refinement to pure_component,
62 ** homogeneous_mixture, or heterogeneous_mixture
63 ** is possible.
64 ** 06/94 - Changed canonical variables from V, H, and S to
65 ** V, H, and G. Also, liquid component model was
66 ** condensed, eliminating instance saturated.
67 **
68 ** 08/94 - Slight structural changes made by Bob Huss to
69 ** allow refinement of non-thermodynamic models,
70 ** and to include Wilson liquid mixture written
71 ** by Boyd Safrit.
72 **
73 ** 08/95 - Addition OF zeros by Jen Stokes to UNIFAC which
74 ** cleared up the wild dimensionality problem it
75 ** was having due to empty sets. However, a mixture
76 ** with ammonia, hydrogen, or carbon_dioxide alone
77 ** is still not feasible IN this version OF
78 ** UNIFAC because they DO not contain any groups or
79 ** subgroups and would CREATE completely NULL
80 ** systems. Also added names to the equations.
81 **
82 ** 03/97 - Addition of SELECT statements for correlations
83 ** by Jen (Stokes) Perry.
84 ** 03/97 - Euthanasia of mixture and alpha_mixture types
85 ** by Ben Allan.
86 **
87 **
88 **
89 ** CONTENTS: ASCEND structure for calculating the basic set of intensive
90 ** thermodynamic properties molar volume, enthalpy, and
91 ** entropy for single and multiple phase streams of pure and
92 ** mixed components. Specify METHODs are included which
93 ** have been designed to provide a means of calculating ideal
94 ** approximations when base models are used. For pure
95 ** component vapors, the ideal gas law can be obtained whereas
96 ** for pure component liquids, incompressibility can be
97 ** specified. Ideal vapor and liquid mixtures are maintained
98 ** by setting all partial molar excess properties to zero.
99 ** Distribution of components among multiple phases can be
100 ** ideally computed using constant relative volatilities.
101 **
102 ** For more rigorous non-ideal calculations, some generalized
103 ** refinements of the base models are provided. For pure
104 ** component vapors, a Pitzer correlation of the two term
105 ** virial equation allows a more accurate compressibility and
106 ** residual calculation. The widely used Rackett correlation
107 ** is accurate in estimating the effect of temperature on
108 ** liquid volumes. Non-ideal vapor mixtures are computed using
109 ** an extension of the Pitzer correlation where the exact
110 ** composition dependence of the second virial coefficient is
111 ** given by statistical mechanics. A reliable UNIFAC model
112 ** estimates non-ideal liquid mixtures. Phase equilibrium
113 ** can be enforced rigorously among multiple phases which
114 ** in turn will allow calculation of the true relative
115 ** volatilities.
116 **
117 **)
118
119 MODEL hgmodel() REFINES cmumodel();
120 END hgmodel;
121
122 MODEL thermodynamic_properties(
123 P WILL_BE pressure;
124 T WILL_BE temperature;
125 ) REFINES hgmodel();
126
127 V IS_A molar_volume;
128 H IS_A molar_energy;
129 G IS_A molar_energy;
130 boundwidth IS_A bound_width;
131
132
133 METHODS
134 METHOD clear;
135 FREE T;
136 FREE P;
137 FREE V;
138 FREE H;
139 FREE G;
140 END clear;
141 METHOD specify;
142 FIX T;
143 FIX P;
144 FIX V;
145 FIX H;
146 FIX G;
147 END specify;
148 METHOD reset;
149 RUN clear;
150 RUN specify;
151 END reset;
152 METHOD scale;
153 T.nominal := T;
154 P.nominal := P;
155 V.nominal := V;
156 H.nominal := sqrt(sqr(H));
157 G.nominal := sqrt(sqr(G));
158 T.lower_bound := 1.0e-8{K};
159 P.lower_bound := 1.0e-8{atm};
160 V.lower_bound := 1.0e-8{liter/g_mole};
161 H.lower_bound := H - boundwidth*H.nominal;
162 G.lower_bound := G - boundwidth*G.nominal;
163 T.upper_bound := T + boundwidth*T.nominal;
164 P.upper_bound := P + boundwidth*P.nominal;
165 V.upper_bound := V + boundwidth*V.nominal;
166 H.upper_bound := H + boundwidth*H.nominal;
167 G.upper_bound := G + boundwidth*G.nominal;
168 END scale;
169
170 END thermodynamic_properties;
171
172 MODEL pure_component(
173 P WILL_BE pressure;
174 T WILL_BE temperature;
175 data WILL_BE td_component_constants;
176 phase IS_A symbol_constant;
177 correlation WILL_BE symbol_constant;
178 ) REFINES thermodynamic_properties();
179 (* no body *)
180 END pure_component;
181
182
183 MODEL partial_component() REFINES hgmodel();
184
185 V IS_A molar_volume;
186 H IS_A molar_energy;
187 G IS_A molar_energy;
188 boundwidth IS_A bound_width;
189
190 METHODS
191 METHOD clear;
192 FREE V;
193 FREE H;
194 FREE G;
195 END clear;
196 METHOD specify;
197 FIX V;
198 FIX H;
199 FIX G;
200 END specify;
201 METHOD reset;
202 RUN clear;
203 RUN specify;
204 END reset;
205 METHOD scale;
206 V.nominal := V;
207 H.nominal := sqrt(sqr(H));
208 G.nominal := sqrt(sqr(G));
209 V.lower_bound := 1.0e-8{liter/g_mole};
210 H.lower_bound := H - boundwidth*H.nominal;
211 G.lower_bound := G - boundwidth*G.nominal;
212 V.upper_bound := V + boundwidth*V.nominal;
213 H.upper_bound := H + boundwidth*H.nominal;
214 G.upper_bound := G + boundwidth*G.nominal;
215 END scale;
216
217 END partial_component;
218
219 (* mixture gone away *)
220
221 MODEL homogeneous_mixture(
222 P WILL_BE pressure;
223 T WILL_BE temperature;
224 options WILL_BE single_phase_options;
225 ) REFINES thermodynamic_properties();
226
227 y[options.ds.components] IS_A mole_fraction;
228 mixture_y_def: SUM[y[i] | i IN options.ds.components] = 1.0;
229 phi[options.phase] IS_A constant;
230 phi[options.phase] :== 1;
231
232 METHODS
233 METHOD clear;
234 FREE y[options.ds.components];
235 END clear;
236 METHOD specify;
237 FIX y[options.ds.components-[options.ds.reference]];
238 END specify;
239 METHOD reset;
240 RUN clear;
241 RUN specify;
242 END reset;
243 METHOD scale;
244 (* P5
245 y[options.ds.components].nominal := 0.5;
246 *)
247 FOR i IN options.ds.components DO
248 y[i].nominal := y[i];
249 END FOR;
250 END scale;
251
252 END homogeneous_mixture;
253
254 MODEL two_phase_mixture(
255 P WILL_BE pressure;
256 T WILL_BE temperature;
257 light WILL_BE homogeneous_mixture;
258 heavy WILL_BE homogeneous_mixture;
259 ) WHERE (
260 P, heavy.P, light.P WILL_BE_THE_SAME;
261 T, heavy.T, light.T WILL_BE_THE_SAME;
262 light.options.phase != heavy.options.phase;
263 light.options.ds.reference == heavy.options.ds.reference;
264 (* constrain the reference component to be a mobile species.*)
265 INTERSECTION[ light.options.ds.components,
266 heavy.options.ds.components,
267 [light.options.ds.reference]
268 ] == [light.options.ds.reference];
269 ) REFINES thermodynamic_properties;
270 (* this MODEL could have species which are not able to
271 * exist in one of the phases.
272 *)
273 components "all the components found in at least one phase"
274 IS_A set OF symbol_constant;
275 components :== [light.options.ds.components,
276 heavy.options.ds.components];
277
278 y[components] "the mole fraction of each species in total mixture"
279 IS_A fraction;
280 mixture_y_def: SUM[y[i] | i IN components] = 1.0;
281
282 mobile "all the components moving between phases"
283 IS_A set OF symbol_constant;
284 mobile :==
285 INTERSECTION[light.options.ds.components,heavy.options.ds.components];
286
287 number_stuck "the number of nondistributing components"
288 IS_A integer_constant;
289 number_stuck :== CARD[components - mobile];
290
291 phases "the names of the phases" IS_A set OF symbol_constant;
292 phases :== [light.options.phase,heavy.options.phase];
293 phi[light.options.phase,heavy.options.phase]
294 "fraction of total moles in each phase" IS_A fraction;
295
296 reference IS_A symbol_constant;
297 reference :== heavy.options.phase;
298
299
300 SELECT (number_stuck)
301 CASE 0: (* this is the standard case *)
302 alpha[components] IS_A partition_coefficient;
303 ave_alpha IS_A partition_coefficient;
304 FOR i IN components CREATE
305 two_phase_mixture_alpha[i]:
306 ave_alpha * light.y[i] = alpha[i] * heavy.y[i];
307 two_phase_mixture_phi[i]:
308 y[i] = phi[light.options.phase] * light.y[i] +
309 phi[heavy.options.phase] * heavy.y[i];
310 END FOR;
311 OTHERWISE:
312 (* this cases needs to be rewritten to allow for species which
313 * are stuck in one phase or the other.
314 *)
315 END SELECT; (* select number_stuck *)
316
317 METHODS
318 METHOD clear;
319 RUN light.clear;
320 RUN heavy.clear;
321 FREE V;
322 FREE H;
323 FREE G;
324 FREE y[components];
325 FREE alpha[components];
326 FREE ave_alpha;
327 FREE phi[phases];
328 END clear;
329 METHOD specify;
330 (* this method buggy if number_stuck > 0 components *)
331 RUN light.specify;
332 FREE light.y[components];
333 RUN heavy.specify;
334 FIX V;
335 FIX H;
336 FIX G;
337 FIX alpha[components];
338 FIX y[light.options.ds.other];
339 FIX phi[light.options.phase];
340 END specify;
341 METHOD scale;
342 alpha[components].lower_bound := 1e-8;
343 (* P5
344 y[components].nominal := 0.5;
345 *)
346 FOR i IN components DO
347 alpha[i].nominal := alpha[i] +1;
348 alpha[i].upper_bound := alpha[i] + boundwidth * alpha[i].nominal;
349 y[i].nominal := y[i];
350 END FOR;
351 RUN light.scale;
352 RUN heavy.scale;
353 END scale;
354 END two_phase_mixture;
355
356 (* alpha_mixture gone away *)
357
358 MODEL td_homogeneous_mixture(
359 P WILL_BE pressure;
360 T WILL_BE temperature;
361 options WILL_BE single_phase_options;
362 ) REFINES homogeneous_mixture;
363
364 ds ALIASES options.ds;
365 FOR i IN ds.components CREATE
366 pure[i] IS_A
367 pure_component(P, T, ds.data[i], options.phase,
368 options.component_thermo_correlation);
369 END FOR;
370 partial[ds.components] IS_A partial_component;
371
372 td_homogeneous_mixture_V:
373 V * (ds.data[ds.reference].Pc / (1{GAS_C}*ds.data[ds.reference].Tc)) =
374 SUM[ y[i] * partial[i].V | i IN ds.components] *
375 (ds.data[ds.reference].Pc /
376 (1{GAS_C}*ds.data[ds.reference].Tc));
377 td_homogeneous_mixture_H:
378 H / (1{GAS_C} * ds.data[ds.reference].Tc) =
379 SUM[ y[i] * partial[i].H | i IN ds.components]/
380 (1{GAS_C} * ds.data[ds.reference].Tc);
381 td_homogeneous_mixture_G:
382 G / (1{GAS_C} * ds.data[ds.reference].Tc) =
383 SUM[ y[i] * partial[i].G | i IN ds.components] /
384 (1{GAS_C} * ds.data[ds.reference].Tc);
385
386
387 METHODS
388
389 METHOD defaults;
390 y[ds.components].lower_bound := 1.0e-12;
391 END defaults;
392 METHOD clear;
393 RUN pure[ds.components].clear;
394 RUN partial[ds.components].clear;
395 FREE y[ds.components];
396 FREE V;
397 FREE H;
398 FREE G;
399 END clear;
400 METHOD specify;
401 RUN pure[ds.components].specify;
402 RUN partial[ds.components].specify;
403 FIX y[ds.components];
404 FREE y[ds.reference];
405 END specify;
406 METHOD scale;
407 RUN pure[ds.components].scale;
408 RUN partial[ds.components].scale;
409 (* P5
410 y[ds.components].nominal := 0.5;
411 *)
412 FOR i IN ds.components DO
413 y[i].nominal := y[i];
414 END FOR;
415 T.nominal := T;
416 P.nominal := P;
417 V.nominal := V;
418 H.nominal := sqrt(sqr(H));
419 G.nominal := sqrt(sqr(G));
420 T.lower_bound := 1.0e-8{K};
421 P.lower_bound := 1.0e-8{atm};
422 V.lower_bound := 1.0e-8{liter/g_mole};
423 H.lower_bound := H - boundwidth*H.nominal;
424 G.lower_bound := G - boundwidth*G.nominal;
425 T.upper_bound := T + boundwidth*T.nominal;
426 P.upper_bound := P + boundwidth*P.nominal;
427 V.upper_bound := V + boundwidth*V.nominal;
428 H.upper_bound := H + boundwidth*H.nominal;
429 G.upper_bound := G + boundwidth*G.nominal;
430 END scale;
431
432 END td_homogeneous_mixture;
433
434 MODEL td_two_phase_mixture( (* was td_alpha_mixture *)
435 P WILL_BE pressure;
436 T WILL_BE temperature;
437 light WILL_BE td_homogeneous_mixture;
438 heavy WILL_BE td_homogeneous_mixture;
439 equilibrated WILL_BE boolean;
440 ) WHERE (
441 P, heavy.P, light.P WILL_BE_THE_SAME;
442 T, heavy.T, light.T WILL_BE_THE_SAME;
443 light.options.phase != heavy.options.phase;
444 light.options.ds.reference == heavy.options.ds.reference;
445 (* constrain the reference component to be a mobile species.*)
446 INTERSECTION[ light.options.ds.components,
447 heavy.options.ds.components,
448 [light.options.ds.reference]
449 ] == [light.options.ds.reference];
450 light.options.ds.components == heavy.options.ds.components;
451 ) REFINES two_phase_mixture;
452
453 data ALIASES heavy.options.ds.data;
454 ds ALIASES heavy.options.ds;
455
456 td_two_phase_mixture_V:
457 V * (data[ds.reference].Pc /
458 ((1{GAS_C})*data[ds.reference].Tc)) =
459 (phi[light.options.phase] * light.V +
460 phi[heavy.options.phase] * heavy.V) *
461 (data[ds.reference].Pc / ((1{GAS_C})*data[ds.reference].Tc));
462
463 td_two_phase_mixture_H:
464 H/((1{GAS_C})*data[ds.reference].Tc) =
465 (phi[light.options.phase] * light.H +
466 phi[heavy.options.phase] * heavy.H) /
467 ((1{GAS_C})*data[ds.reference].Tc);
468
469 td_two_phase_mixture_G:
470 G/((1{GAS_C})*data[ds.reference].Tc) =
471 (phi[light.options.phase] * light.G +
472 phi[heavy.options.phase] * heavy.G) /
473 ((1{GAS_C})*data[ds.reference].Tc);
474
475 FOR i IN ds.components CREATE
476 equil_condition[i]:
477 (light.partial[i].G - heavy.partial[i].G ) / (1E5{J/mole}) = 0;
478 (* oops, can't scale with G0 because hydrogen G0 = 0.0 *)
479 END FOR;
480
481 (*
482 NOTES 'default' FOR
483 equilibrated {
484 The equilibrium condition can be relaxed, effectively making this
485 model devolve to td_heterogeneous mixture. It is arguable whether
486 this type and its immediate ancestor should be merged.
487 };
488 END NOTES;
489 *)
490
491 WHEN (equilibrated)
492 CASE TRUE:
493 USE equil_condition;
494 END WHEN;
495
496
497
498 METHODS
499 METHOD defaults;
500 V := 30{liter/mole};
501 equilibrated := TRUE;
502 END defaults;
503 METHOD clear;
504 RUN two_phase_mixture::clear;
505 FREE V;
506 FREE H;
507 FREE G;
508 END clear;
509
510 METHOD specify;
511 FIX P;
512 FIX phi[light.options.phase];
513 FIX y[light.options.ds.other];
514 IF equilibrated THEN
515 ave_alpha := 1.0;
516 FIX ave_alpha;
517 ELSE
518 FIX T;
519 FIX alpha[ds.components];
520 END IF;
521 END specify;
522
523 METHOD reset_to_massbal;
524 RUN clear;
525 equilibrated := FALSE;
526 RUN specify;
527 END reset_to_massbal;
528
529 METHOD reset_to_fullthermo;
530 RUN clear;
531 equilibrated := TRUE;
532 RUN specify;
533 END reset_to_fullthermo;
534
535 METHOD reset;
536 (* reset for whatever value of equilibrated has been previously set *)
537 RUN clear;
538 RUN specify;
539 END reset;
540
541 METHOD scale;
542 RUN two_phase_mixture::scale;
543 T.nominal := T;
544 P.nominal := P;
545 V.nominal := V;
546 H.nominal := sqrt(sqr(H));
547 G.nominal := sqrt(sqr(G));
548 T.lower_bound := 1.0e-8{K};
549 P.lower_bound := 1.0e-8{atm};
550 V.lower_bound := 1.0e-8{liter/g_mole};
551 H.lower_bound := H - boundwidth*H.nominal;
552 G.lower_bound := G - boundwidth*G.nominal;
553 T.upper_bound := T + boundwidth*T.nominal;
554 P.upper_bound := P + boundwidth*P.nominal;
555 V.upper_bound := V + boundwidth*V.nominal;
556 H.upper_bound := H + boundwidth*H.nominal;
557 G.upper_bound := G + boundwidth*G.nominal;
558 END scale;
559
560 END td_two_phase_mixture;
561
562 (* td_equilibrium_mixture merged with td_two_phase *)
563
564 MODEL vapor_component(
565 P WILL_BE pressure;
566 T WILL_BE temperature;
567 data WILL_BE td_component_constants;
568 correlation WILL_BE symbol_constant;
569 ) REFINES pure_component(phase :== 'vapor';);
570 (* no body *)
571 END vapor_component;
572
573 MODEL liquid_component(
574 P WILL_BE pressure;
575 T WILL_BE temperature;
576 data WILL_BE td_component_constants;
577 correlation WILL_BE symbol_constant;
578 ) REFINES pure_component(phase :== 'liquid';);
579
580 VP IS_A pressure;
581 END liquid_component;
582
583 MODEL Pitzer_vapor_component(
584 P WILL_BE pressure;
585 T WILL_BE temperature;
586 data WILL_BE td_component_constants;
587 correlation WILL_BE symbol_constant WITH_VALUE 'Pitzer';
588 ) REFINES vapor_component;
589
590 Pitzer_component_V:
591 P*V/(1{GAS_C})/data.Tc = T/data.Tc + (P/data.Pc)*
592 (0.083 - 0.422*(data.Tc/T)^1.6 + data.omega*
593 (0.139 - 0.172*(data.Tc/T)^4.2));
594
595 Pitzer_component_H:
596 H/(1{GAS_C})/data.Tc = data.H0/(1{GAS_C})/data.Tc +
597 data.cpvapa*(T - data.T0)/(1{GAS_C})/data.Tc +
598 data.cpvapb*(T^2 - data.T0^2)/2/(1{GAS_C})/data.Tc +
599 data.cpvapc*(T^3 - data.T0^3)/3/(1{GAS_C})/data.Tc +
600 data.cpvapd*(T^4 - data.T0^4)/4/(1{GAS_C})/data.Tc +
601 (P/data.Pc)*(0.083 - 1.097*(data.Tc/T)^1.6 + data.omega*
602 (0.139 - 0.894*(data.Tc/T)^4.2));
603
604 Pitzer_component_G:
605 G/(1{GAS_C})/data.Tc = data.G0/(1{GAS_C})/data.Tc -
606 (data.H0 - data.G0)*(T/data.T0 - 1)/(1{GAS_C})/data.Tc -
607 data.cpvapa*(T*ln(T/data.T0) - T + data.T0)/1{GAS_C}/data.Tc -
608 data.cpvapb*(T^2 - 2*T*data.T0 + data.T0^2)/
609 (2*1{GAS_C}*data.Tc) -
610 data.cpvapc*(T^3/2 - 3*T*data.T0^2/2 + data.T0^3)/
611 (3*1{GAS_C}*data.Tc) -
612 data.cpvapd*(T^4/3 - 4*T*data.T0^3/3 + data.T0^4)/4/(1{GAS_C})/data.Tc +
613 T*ln(P/data.P0)/data.Tc +
614 (P/data.Pc)*
615 (0.083 - 0.422*(data.Tc/T)^1.6 + data.omega*
616 (0.139 - 0.172*(data.Tc/T)^4.2));
617
618 METHODS
619
620 METHOD defaults;
621 T.lower_bound := 1.0{K};
622 P.lower_bound := 0.1{atm};
623 T.upper_bound := 1000{K};
624 P.upper_bound := 20{atm};
625 V := 24{liter/mole};
626 END defaults;
627
628 METHOD specify;
629 FIX T;
630 FIX P;
631 END specify;
632 END Pitzer_vapor_component;
633
634 MODEL Rackett_liquid_component(
635 P WILL_BE pressure;
636 T WILL_BE temperature;
637 data WILL_BE td_component_constants;
638 correlation WILL_BE symbol_constant WITH_VALUE 'Rackett';
639 vp_correlation WILL_BE integer_constant;
640 ) REFINES liquid_component;
641 (*Rackett_component_VP:*)
642 SELECT (vp_correlation)
643 CASE 1:
644 ln(VP/data.Pc)*T/data.Tc*abs(data.Tc - T) =
645 (data.vpa * abs(1.0 - T/data.Tc) +
646 data.vpb * abs(1.0 - T/data.Tc)^1.5 +
647 data.vpc * abs(1.0 - T/data.Tc)^3.0 +
648 data.vpd * abs(1.0 - T/data.Tc)^6.0) * (data.Tc - T);
649 CASE 2:
650 ln(VP/1{bar}) / ln(data.Pc/1{bar}) =
651 (data.vpa - data.vpb/T + data.vpc * ln(T/1{K}) +
652 data.vpd * VP/sqr(T)) / ln(data.Pc/1{bar});
653 CASE 3:
654 ln(VP/1{bar}) / ln(data.Pc/1{bar}) =
655 (data.vpa - data.vpb/(data.vpc + T)) / ln(data.Pc/1{bar});
656 END SELECT;
657
658 Rackett_component_V: V/data.Vliq =
659 data.Zc^(abs(1.0 - T/data.Tc)^(2/7))/
660 data.Zc^(abs(1.0 - data.Tliq/data.Tc)^(2/7));
661
662 Rackett_component_H:
663 H / (1{GAS_C} * data.Tc) = data.H0/(1{GAS_C})/data.Tc +
664 (T - data.T0) * (data.cpvapa/(1{GAS_C})/data.Tc) +
665 (T^2 - data.T0^2) * (data.cpvapb/2/(1{GAS_C})/data.Tc) +
666 (T^3 - data.T0^3) * (data.cpvapc/3/(1{GAS_C})/data.Tc) +
667 (T^4 - data.T0^4) * (data.cpvapd/4/(1{GAS_C})/data.Tc) +
668 (VP/data.Pc)*
669 (0.083 - 1.097*(data.Tc/T)^1.6 + data.omega*
670 (0.139 - 0.894*(data.Tc/T)^4.2)) -
671 (data.Hv/(1{GAS_C})/data.Tc)*
672 abs((data.Tc-T)/(data.Tc-data.Tb))^0.38 +
673 (P - VP)*(data.Vliq/(1{GAS_C})/data.Tc)*
674 (data.Zc^(abs(1.0 - T/data.Tc)^(2/7))/
675 data.Zc^(abs(1.0 - data.Tliq/data.Tc)^(2/7)))*(1.0 -
676 (-2/7)*(T/data.Tc)*(abs(1 - T/data.Tc)^(-5/7))*ln(data.Zc));
677
678 Rackett_component_G: G/(1{GAS_C}*data.Tc) =
679 data.G0/(1{GAS_C})/data.Tc -
680 (T/data.T0 - 1) * ((data.H0 - data.G0)/(1{GAS_C})/data.Tc) -
681 (T*ln(T/data.T0) - T + data.T0) *
682 (data.cpvapa/(1{GAS_C})/data.Tc) -
683 (T^2 - T*(2*data.T0) + data.T0^2) *
684 (data.cpvapb/2/(1{GAS_C})/data.Tc) -
685 (T^3/2 - T*(3/2*data.T0^2) + data.T0^3) *
686 (data.cpvapc/3/(1{GAS_C})/data.Tc) -
687 (T^4/3 - T*(4/3*data.T0^3) + data.T0^4) *
688 (data.cpvapd/4/(1{GAS_C})/data.Tc) +
689 T*ln(VP/data.P0)/data.Tc +
690 (VP/data.Pc) *
691 (0.083 - 0.422*(data.Tc/T)^1.6 + data.omega *
692 (0.139 - 0.172*(data.Tc/T)^4.2)) +
693 (P - VP)*(data.Vliq/(1{GAS_C})/data.Tc) *
694 (data.Zc^(abs(1.0 - T/data.Tc)^(2/7))/
695 data.Zc^(abs(1.0 - data.Tliq/data.Tc)^(2/7)));
696
697
698 METHODS
699 METHOD defaults;
700 VP.lower_bound := 1.0e-12{Pa};
701 V := 0.1{liter/mole};
702 T.lower_bound := 1.0{K};
703 P.lower_bound := 0.1{atm};
704 T.upper_bound := 1000{K};
705 P.upper_bound := 20{atm};
706 END defaults;
707
708 METHOD clear;
709 RUN liquid_component::clear;
710 FREE VP;
711 END clear;
712 METHOD specify;
713 FIX T;
714 FIX P;
715 END specify;
716 METHOD scale;
717 RUN liquid_component::scale;
718 VP.nominal := VP;
719 VP.lower_bound := 1.0e-8{atm};
720 VP.upper_bound := VP + boundwidth*VP.nominal;
721 END scale;
722
723
724 END Rackett_liquid_component;
725
726 MODEL Pitzer_partials(
727 P WILL_BE pressure;
728 T WILL_BE temperature;
729 V WILL_BE molar_volume;
730 H WILL_BE molar_energy;
731 G WILL_BE molar_energy;
732 boundwidth WILL_BE bound_width;
733 options WILL_BE single_phase_options;
734 y[options.ds.components] WILL_BE mole_fraction;
735 pure[options.ds.components] WILL_BE Pitzer_vapor_component;
736 partial[options.ds.components] WILL_BE partial_component;
737 ) REFINES hgmodel();
738
739 ds ALIASES options.ds;
740
741 FOR i IN ds.components CREATE
742 Pitzer_mixture_V[i]:
743 (partial[i].V - pure[i].V)*(ds.data[i].Pc/
744 (1{GAS_C}*ds.data[i].Tc)) =
745 -1.0*(0.083 - 0.422*(ds.data[i].Tc/T)^1.6 + ds.data[i].omega*
746 (0.139 - 0.172*(ds.data[i].Tc/T)^4.2))*(1.0 - y[i])^2 +
747 0.50*(1.0 - y[i])*
748 SUM[y[j]*((1.0 + (ds.data[j].Vc/ds.data[i].Vc)^(1/3))^3/
749 (1.0 + ds.data[j].Zc/ds.data[i].Zc))*
750 (0.083 - 0.422*(sqrt(ds.data[i].Tc*ds.data[j].Tc)/T)^1.6 +
751 0.5*(ds.data[i].omega + ds.data[j].omega)*
752 (0.139 - 0.172*(sqrt(ds.data[i].Tc*ds.data[j].Tc)/T)^4.2))
753 | j IN ds.components - [i]];
754 Pitzer_mixture_H[i]:
755 (partial[i].H - pure[i].H)/((1{GAS_C})*ds.data[i].Tc) =
756 -(P/ds.data[i].Pc)*
757 (0.083 - 1.097*(ds.data[i].Tc/T)^1.6 + ds.data[i].omega*
758 (0.139 - 0.894*(ds.data[i].Tc/T)^4.2))*(1.0 - y[i])^2 +
759 (1.0 - y[i])*(P/(ds.data[i].Pc*2))*
760 SUM[y[j]*((1.0 + (ds.data[j].Vc/ds.data[i].Vc)^(1/3))^3/
761 (1.0 + ds.data[j].Zc/ds.data[i].Zc))*
762 (0.083 - 1.097*(sqrt(ds.data[i].Tc*ds.data[j].Tc)/T)^1.6 +
763 0.5*(ds.data[i].omega + ds.data[j].omega)*
764 (0.139 - 0.894*(sqrt(ds.data[i].Tc*ds.data[j].Tc)/T)^4.2))
765 | j IN ds.components - [i]];
766 Pitzer_mixture_G[i]:
767 (partial[i].G - pure[i].G - (1{GAS_C})*T*ln(y[i])) /
768 (1{GAS_C}*ds.data[i].Tc) =
769 -(P/ds.data[i].Pc)*
770 (0.083 - 0.422*(ds.data[i].Tc/T)^1.6 + ds.data[i].omega*
771 (0.139 - 0.172*(ds.data[i].Tc/T)^4.2))*(1.0 - y[i])^2 +
772 (1.0 - y[i])*(P/(ds.data[i].Pc*2))*
773 SUM[y[j]*((1.0 + (ds.data[j].Vc/ds.data[i].Vc)^(1/3))^3/
774 (1.0 + ds.data[j].Zc/ds.data[i].Zc))*
775 (0.083 - 0.422*(sqrt(ds.data[i].Tc*ds.data[j].Tc)/T)^1.6 +
776 0.5*(ds.data[i].omega + ds.data[j].omega)*
777 (0.139 - 0.172*(sqrt(ds.data[i].Tc*ds.data[j].Tc)/T)^4.2))
778 | j IN ds.components - [i]];
779 END FOR;
780
781 METHODS
782
783 METHOD defaults;
784 T.lower_bound := 1.0{K};
785 P.lower_bound := 1.0{Pa};
786 V := 24{liter/mole};
787 partial[ds.components].V := 24{liter/mole};
788 pure[ds.components].V := 24{liter/mole};
789 END defaults;
790
791 METHOD specify;
792 FREE partial[ds.components].V;
793 FREE partial[ds.components].H;
794 FREE partial[ds.components].G;
795 END specify;
796
797 END Pitzer_partials;
798
799 (*UNIFAC_partials is a model which helps to select mixing rules at the
800 * liquid mixture level, hiding unneccessary code from users.
801 *)
802 MODEL UNIFAC_partials(
803 P WILL_BE pressure;
804 T WILL_BE temperature;
805 V WILL_BE molar_volume;
806 H WILL_BE molar_energy;
807 G WILL_BE molar_energy;
808 boundwidth WILL_BE bound_width;
809 options WILL_BE single_phase_options;
810 y[options.ds.components] WILL_BE mole_fraction;
811 pure[options.ds.components] WILL_BE Rackett_liquid_component;
812 partial[options.ds.components] WILL_BE partial_component;
813 ) REFINES hgmodel();
814
815 ds ALIASES options.ds;
816 subgroups IS_A set OF symbol_constant;
817 groups IS_A set OF integer_constant;
818 comps[subgroups] IS_A set OF symbol_constant;
819 rv[ds.components] IS_A UNIFAC_size;
820 qs[ds.components] IS_A UNIFAC_size;
821 Jv[ds.components] IS_A factor;
822 Ls[ds.components] IS_A factor;
823 theta[subgroups] IS_A factor;
824 eta[subgroups] IS_A factor;
825 uc IS_A UNIFAC_constants;
826
827 subgroups :== UNION[ds.data[i].subgroups | i IN ds.components];
828 groups :== UNION[ds.data[i].groups | i IN ds.components];
829 FOR k IN subgroups CREATE
830 comps[k] :== [i IN ds.components | k IN ds.data[i].subgroups];
831 END FOR;
832 FOR k IN subgroups CREATE
833 UNIFAC_mixture_theta[k]:
834 theta[k] = uc.Q[k]*SUM[ds.data[i].nu[k]*y[i] | i IN comps[k]];
835
836 UNIFAC_mixture_eta[k]:
837 eta[k] =
838 SUM[theta[m] | m IN subgroups*uc.sub[uc.group[k]]] +
839 SUM[SUM[theta[m]*exp(-uc.a[g][uc.group[k]]/T)
840 | m IN subgroups*uc.sub[g]]
841 | g IN groups - [uc.group[k]]];
842 END FOR;
843
844 FOR i IN ds.components CREATE
845 rv[i] :== SUM[0, ds.data[i].nu[k]*uc.R[k]
846 | k IN ds.data[i].subgroups];
847 qs[i] :== SUM[0, ds.data[i].nu[k]*uc.Q[k]
848 | k IN ds.data[i].subgroups];
849 END FOR;
850 FOR i IN ds.components CREATE
851 UNIFAC_mixture_rv[i]:
852 rv[i] = Jv[i]*SUM[rv[j]*y[j] | j IN ds.components];
853
854 UNIFAC_mixture_qs[i]:
855 qs[i] = Ls[i]*SUM[qs[j]*y[j] | j IN ds.components];
856 partial[i].V = pure[i].V;
857
858 UNIFAC_mixture_H_excess[i]:
859 (partial[i].H - pure[i].H)/(1{GAS_C})/ds.data[i].Tc =
860 SUM[0,theta[k]*
861 SUM[0,SUM[0,theta[n]*
862 ((uc.a[g][uc.group[k]] -
863 uc.a[uc.group[n]][uc.group[k]])/ds.data[i].Tc)*
864 exp(-(uc.a[g][uc.group[k]] +
865 uc.a[uc.group[n]][uc.group[k]])/T)*
866 SUM[0,ds.data[i].nu[m]*uc.Q[m]
867 | m IN ds.data[i].subgroups*uc.sub[g]]
868 | g IN ds.data[i].groups - [uc.group[n]]]
869 | n IN subgroups]/eta[k]/eta[k]
870 | k IN subgroups] -
871 SUM[0,(ds.data[i].nu[k]*uc.Q[k]/(
872 SUM[0,ds.data[i].nu[m]*uc.Q[m]
873 | m IN ds.data[i].subgroups*uc.sub[uc.group[k]]] +
874 SUM[0,SUM[0,ds.data[i].nu[m]*uc.Q[m] *
875 exp(-uc.a[g][uc.group[k]]/T)
876 | m IN ds.data[i].subgroups*uc.sub[g]]
877 | g IN ds.data[i].groups - [uc.group[k]]]))*
878 SUM[0,SUM[0,theta[n]*
879 ((uc.a[g][uc.group[k]] -
880 uc.a[uc.group[n]][uc.group[k]])/ds.data[i].Tc)*
881 exp(-(uc.a[g][uc.group[k]] +
882 uc.a[uc.group[n]][uc.group[k]])/T)*
883 SUM[0,ds.data[i].nu[m]*uc.Q[m]
884 | m IN ds.data[i].subgroups*uc.sub[g]]
885 | g IN ds.data[i].groups - [uc.group[n]]]
886 | n IN subgroups]/eta[k]
887 | k IN ds.data[i].subgroups];
888
889 UNIFAC_mixture_G_excess[i]:
890 (partial[i].G - pure[i].G - (1{GAS_C})*T*ln(y[i])) /
891 (1{GAS_C}*ds.data[i].Tc) =
892 (1.0 - Jv[i] + ln(Jv[i]) -
893 5.0*qs[i]*(1.0 - Jv[i]/Ls[i] + ln(Jv[i]/Ls[i])) +
894 qs[i]*(1 - ln(Ls[i])))*T/ds.data[i].Tc -
895 SUM[0,theta[k]*(
896 SUM[0,ds.data[i].nu[m]*uc.Q[m]
897 | m IN ds.data[i].subgroups*uc.sub[uc.group[k]]] +
898 SUM[0,SUM[0,ds.data[i].nu[m] * uc.Q[m] *
899 exp(-uc.a[g][uc.group[k]]/T)
900 | m IN ds.data[i].subgroups*uc.sub[g]]
901 | g IN ds.data[i].groups - [uc.group[k]]])/eta[k]
902 | k IN subgroups]*T/ds.data[i].Tc +
903 SUM[0,ds.data[i].nu[k]*uc.Q[k]*ln((
904 SUM[0,ds.data[i].nu[m]*uc.Q[m]
905 | m IN ds.data[i].subgroups*uc.sub[uc.group[k]]] +
906 SUM[0,SUM[0, ds.data[i].nu[m] * uc.Q[m] *
907 exp(-uc.a[g][uc.group[k]]/T)
908 | m IN ds.data[i].subgroups*uc.sub[g]]
909 | g IN ds.data[i].groups - [uc.group[k]]])/eta[k])
910 | k IN ds.data[i].subgroups]*T/ds.data[i].Tc;
911 END FOR;
912
913
914 FOR i IN ds.components CREATE
915 rv[i] :== SUM[0,ds.data[i].nu[k]*uc.R[k]
916 | k IN ds.data[i].subgroups];
917 qs[i] :== SUM[0,ds.data[i].nu[k]*uc.Q[k]
918 | k IN ds.data[i].subgroups];
919 END FOR;
920
921 METHODS
922 METHOD defaults;
923 T.lower_bound := 1.0{K};
924 P.lower_bound := 1.0{Pa};
925 Jv[ds.components].lower_bound := 1.0e-8;
926 Ls[ds.components].lower_bound := 1.0e-8;
927 theta[subgroups].lower_bound := 0.0;
928 eta[subgroups].lower_bound := 0.0;
929 V := 0.1{liter/mole};
930 partial[ds.components].V := 0.1{liter/mole};
931 pure[ds.components].V := 0.1{liter/mole};
932 END defaults;
933 METHOD clear;
934 FREE Jv[ds.components];
935 FREE Ls[ds.components];
936 FREE theta[subgroups];
937 FREE eta[subgroups];
938 END clear;
939 METHOD specify;
940 FREE partial[ds.components].V;
941 FREE partial[ds.components].H;
942 FREE partial[ds.components].G;
943 END specify;
944
945 METHOD scale;
946 (* we may need to further investigate the scaling on these
947 * when we approach small numbers for Ls, Jv.
948 *)
949 FOR i IN ds.components DO
950 Ls[i].nominal := Ls[i];
951 Jv[i].nominal := Jv[i];
952 Ls[i].upper_bound := Ls[i]+ boundwidth*Ls[i].nominal;
953 Jv[i].upper_bound := Jv[i] + boundwidth*Jv[i].nominal;
954 END FOR;
955 FOR j IN subgroups DO
956 theta[j].nominal := theta[j];
957 eta[j].nominal := eta[j];
958 theta[j].upper_bound := theta[j] + boundwidth*theta[j].nominal;
959 eta[j].upper_bound := eta[j] + boundwidth*eta[j].nominal;
960 END FOR;
961 END scale;
962 END UNIFAC_partials;
963
964 MODEL Wilson_partials(
965 P WILL_BE pressure;
966 T WILL_BE temperature;
967 V WILL_BE molar_volume;
968 H WILL_BE molar_energy;
969 G WILL_BE molar_energy;
970 boundwidth WILL_BE bound_width;
971 options WILL_BE single_phase_options;
972 y[options.ds.components] WILL_BE mole_fraction;
973 pure[options.ds.components] WILL_BE Rackett_liquid_component;
974 partial[options.ds.components] WILL_BE partial_component;
975 ) REFINES hgmodel();
976
977 ds ALIASES options.ds;
978 lambda[ds.components][ds.components] IS_A factor;
979
980 FOR i IN ds.components CREATE
981 FOR j IN ds.components CREATE
982 lambda[i][j] = (pure[j].V / pure[i].V) *
983 exp(-pure[i].data.del_ip[pure[j].data.formula] /
984 (1{GAS_C} * T));
985 END FOR;
986 END FOR;
987
988 FOR i IN ds.components CREATE
989 partial[i].V = pure[i].V;
990
991 Wilson_mixture_G_excess[i]:
992 partial[i].G - pure[i].G - (1{GAS_C})*T*ln(y[i]) =
993 (1{GAS_C})*T*(-ln( SUM[y[j]*lambda[i][j] | j IN ds.components])
994 + 1 -
995 SUM[(y[k] * lambda[k][i]) / SUM[ y[j] * lambda[k][j]
996 | j IN ds.components] | k IN ds.components]);
997
998 Wilson_mixture_H_excess[i]:
999 partial[i].H - pure[i].H =
1000 (SUM[ y[j] * lambda[i][j] *
1001 pure[i].data.del_ip[pure[j].data.formula]
1002 | j IN ds.components]) /
1003 (SUM[y[s]*lambda[i][s]
1004 | s IN ds.components]) -
1005 (SUM[y[k]*lambda[k][i] | k IN ds.components]*
1006 SUM[SUM[y[l]*lambda[m][l]*
1007 pure[m].data.del_ip[pure[l].data.formula]
1008 | l IN ds.components]
1009 | m IN ds.components]) / (
1010 SUM[SUM[y[n]*lambda[o][n]
1011 | n IN ds.components]
1012 | o IN ds.components])^2 +
1013 (SUM[y[p]*lambda[p][i]*
1014 pure[p].data.del_ip[pure[i].data.formula]
1015 | p IN ds.components]) /
1016 (SUM[SUM[y[q]*lambda[r][q]
1017 | q IN ds.components]
1018 | r IN ds.components]);
1019 END FOR;
1020
1021 METHODS
1022 METHOD defaults;
1023 V := 0.1{liter/mole};
1024 partial[ds.components].V := 0.1{liter/mole};
1025 pure[ds.components].V := 0.1{liter/mole};
1026 END defaults;
1027 METHOD clear;
1028 FREE lambda[ds.components][ds.components];
1029 END clear;
1030 METHOD specify;
1031 FREE partial[ds.components].V;
1032 FREE partial[ds.components].G;
1033 FREE partial[ds.components].H;
1034 END specify;
1035 METHOD scale;
1036 FOR i IN ds.components DO
1037 FOR j IN ds.components DO
1038 lambda[i][j].nominal :=
1039 lambda[i][j];
1040 lambda[i][j].lower_bound := lambda[i][j] -
1041 boundwidth*lambda[i][j].nominal;
1042 lambda[i][j].upper_bound := lambda[i][j] +
1043 boundwidth*lambda[i][j].nominal;
1044 END FOR;
1045 END FOR;
1046 END scale;
1047
1048 END Wilson_partials;
1049
1050 MODEL vapor_mixture(
1051 P WILL_BE pressure;
1052 T WILL_BE temperature;
1053 options WILL_BE vapor_phase_options;
1054 ) REFINES td_homogeneous_mixture;
1055
1056 SELECT (options.mixture_thermo_correlation)
1057 CASE 'Pitzer':
1058 FOR i IN ds.components CREATE
1059 pure[i] IS_REFINED_TO
1060 Pitzer_vapor_component(P, T, ds.data[i],
1061 options.component_thermo_correlation);
1062 END FOR;
1063 Pitzer_mixing_rule IS_A
1064 Pitzer_partials(P, T, V, H, G, boundwidth, options,
1065 y, pure, partial);
1066 END SELECT;
1067
1068 METHODS
1069 METHOD specify;
1070 RUN td_homogeneous_mixture::specify;
1071 IF options.mixture_thermo_correlation = 'Pitzer' THEN
1072 RUN Pitzer_mixing_rule.specify;
1073 END IF;
1074 END specify;
1075
1076 END vapor_mixture;
1077
1078 MODEL liquid_mixture(
1079 P WILL_BE pressure;
1080 T WILL_BE temperature;
1081 options WILL_BE liquid_phase_options;
1082 ) REFINES td_homogeneous_mixture;
1083
1084 SELECT (options.mixture_thermo_correlation)
1085 CASE 'UNIFAC':
1086 FOR i IN ds.components CREATE
1087 pure[i] IS_REFINED_TO
1088 Rackett_liquid_component(P, T, ds.data[i],
1089 options.component_thermo_correlation,
1090 ds.data[i].vp_correlation);
1091 END FOR;
1092
1093 UNIFAC_mixing_rule IS_A
1094 UNIFAC_partials(P, T, V, H, G, boundwidth,
1095 options, y, pure, partial);
1096
1097 CASE 'Wilson':
1098 Wilson_mixing_rule IS_A
1099 Wilson_partials(P, T, V, H, G, boundwidth,
1100 options, y, pure, partial);
1101 END SELECT; (* select correlation *)
1102
1103 (* The Rackett_liquid_component may not
1104 be invariant across all future liquid_mixture options
1105 and may need to be moved inside the select IN the future.
1106 It is now outside to prevent warnings about undefined names
1107 due to refinements inside select statements not being
1108 recognized when a file is imported *)
1109 FOR i IN ds.components CREATE
1110 pure[i] IS_REFINED_TO
1111 Rackett_liquid_component(P, T, ds.data[i],
1112 options.component_thermo_correlation,
1113 ds.data[i].vp_correlation);
1114 END FOR;
1115 METHODS
1116 METHOD clear;
1117 RUN td_homogeneous_mixture::clear;
1118 IF options.mixture_thermo_correlation = 'UNIFAC' THEN
1119 RUN UNIFAC_mixing_rule.clear;
1120 END IF;
1121 IF options.mixture_thermo_correlation = 'Wilson' THEN
1122 RUN Wilson_mixing_rule.clear;
1123 END IF;
1124 END clear;
1125
1126 METHOD specify;
1127 RUN td_homogeneous_mixture::specify;
1128 SWITCH (options.mixture_thermo_correlation)
1129 CASE 'UNIFAC':
1130 FREE partial[ds.components].V;
1131 FREE partial[ds.components].H;
1132 FREE partial[ds.components].G;
1133 (* We ought to have break/fallthrough behavior here.
1134 * There's a sneaky way to get it without deep grammar changes.
1135 *)
1136 CASE 'Wilson':
1137 FREE partial[ds.components].V;
1138 FREE partial[ds.components].H;
1139 FREE partial[ds.components].G;
1140 END SWITCH;
1141 END specify;
1142 METHOD scale;
1143 RUN td_homogeneous_mixture::scale;
1144 IF options.mixture_thermo_correlation = 'UNIFAC' THEN
1145 RUN UNIFAC_mixing_rule.scale;
1146 END IF;
1147 IF options.mixture_thermo_correlation = 'Wilson' THEN
1148 RUN Wilson_mixing_rule.scale;
1149 END IF;
1150 END scale;
1151 END liquid_mixture;
1152
1153 MODEL td_VLE_mixture(
1154 P WILL_BE pressure;
1155 T WILL_BE temperature;
1156 light WILL_BE vapor_mixture;
1157 heavy WILL_BE liquid_mixture;
1158 equilibrated WILL_BE boolean;
1159 ) WHERE (
1160 P, heavy.P, light.P WILL_BE_THE_SAME;
1161 T, heavy.T, light.T WILL_BE_THE_SAME;
1162 light.options.phase != heavy.options.phase;
1163 light.options.ds.reference == heavy.options.ds.reference;
1164 (* constrain the reference component to be a mobile species.*)
1165 INTERSECTION[ light.options.ds.components,
1166 heavy.options.ds.components,
1167 [light.options.ds.reference]
1168 ] == [light.options.ds.reference];
1169 light.options.ds.components == heavy.options.ds.components;
1170 light.options.phase == 'vapor';
1171 heavy.options.phase == 'liquid';
1172 ) REFINES td_two_phase_mixture;
1173 (* no body *)
1174
1175 METHODS
1176 (* uses old ones in td_two_phase_mixture *)
1177 END td_VLE_mixture;
1178
1179 MODEL murphree_equilibrium_mixture(
1180 P WILL_BE pressure;
1181 T WILL_BE temperature;
1182 light WILL_BE vapor_mixture;
1183 heavy WILL_BE liquid_mixture;
1184 equilibrated WILL_BE boolean;
1185 ) WHERE (
1186 P, heavy.P, light.P WILL_BE_THE_SAME;
1187 T, heavy.T, light.T WILL_BE_THE_SAME;
1188 light.options.phase != heavy.options.phase;
1189 light.options.ds.reference == heavy.options.ds.reference;
1190 (* constrain the reference component to be a mobile species.*)
1191 INTERSECTION[ light.options.ds.components,
1192 heavy.options.ds.components,
1193 [light.options.ds.reference]
1194 ] == [light.options.ds.reference];
1195 light.options.ds.components == heavy.options.ds.components;
1196 light.options.phase == 'vapor';
1197 heavy.options.phase == 'liquid';
1198 ) REFINES td_VLE_mixture;
1199
1200 booleqnlabel: equilibrated == TRUE; (* logical constraint *)
1201
1202 alpha[ds.components],
1203 ave_alpha IS_REFINED_TO relative_volatility;
1204
1205 vap_eq IS_A vapor_mixture(light.P, light.T, light.options);
1206
1207 equil_alpha[ds.components] IS_A relative_volatility;
1208 ref_y[ds.components] IS_A fraction;
1209 murph_eff IS_A factor;
1210
1211 murphee_equil_ref_y_def:
1212 SUM[ref_y[ds.components]] = 1;
1213
1214 FOR i IN ds.components CREATE
1215 murphee_equil_vap_eq_y_def[i]:
1216 vap_eq.y[i] = equil_alpha[i] * heavy.y[i];
1217 (vap_eq.partial[i].G - heavy.partial[i].G)/1E+5 {J/mole} =0;
1218 END FOR;
1219
1220 FOR i IN ds.other CREATE
1221 murphee_equil_murph_eff_def[i]:
1222 murph_eff*(vap_eq.y[i] - ref_y[i]) = light.y[i] - ref_y[i];
1223 END FOR;
1224
1225 METHODS
1226 METHOD clear;
1227 RUN td_VLE_mixture::clear;
1228 RUN vap_eq.clear;
1229 FREE equil_alpha[ds.components];
1230 FREE ref_y[ds.components];
1231 FREE murph_eff;
1232 END clear;
1233 METHOD specify;
1234 RUN td_two_phase_mixture::specify;
1235 RUN vap_eq.specify;
1236 FREE vap_eq.y[ds.components];
1237 FREE equil_alpha[ds.components];
1238 FIX ref_y[ds.other];
1239 FIX murph_eff;
1240 FREE T;
1241 END specify;
1242
1243 END murphree_equilibrium_mixture;

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22