/[ascend]/trunk/models/thermodynamics.a4l
ViewVC logotype

Contents of /trunk/models/thermodynamics.a4l

Parent Directory Parent Directory | Revision Log Revision Log


Revision 576 - (show annotations) (download) (as text)
Tue May 9 03:42:08 2006 UTC (14 years, 2 months ago) by johnpye
File MIME type: text/x-ascend
File size: 45841 byte(s)
Changed all cases of *.fixed := {TRUE,FALSE} to 'FIX' and 'FREE' statements.
1 REQUIRE "components.a4l";
2 (* => components.a4l, atoms.a4l, measures.a4l, system.a4l, basemodel.a4l *)
3 REQUIRE "phases.a4l";
4 (* => phases.a4l, system.a4l, basemodel.a4l *)
5 PROVIDE "thermodynamics.a4l";
6
7 (*
8 * thermodynamics.a4l
9 * by Arthur W. Westerberg, Ben Allan, and Vicente Rico-Ramirez
10 * Part of the ASCEND Library
11 * $Date: 1998/06/17 19:33:42 $
12 * $Revision: 1.6 $
13 * $Author: mthomas $
14 * $Source: /afs/cs.cmu.edu/project/ascend/Repository/models/thermodynamics.a4l,v $
15 *
16 * This file is part of the ASCEND Modeling Library.
17 *
18 * Copyright (C) 1998 Carnegie Mellon University
19 *
20 * The ASCEND Modeling Library is free software; you can redistribute
21 * it and/or modify it under the terms of the GNU General Public
22 * License as published by the Free Software Foundation; either
23 * version 2 of the License, or (at your option) any later version.
24 *
25 * The ASCEND Modeling Library is distributed in hope that it
26 * will be useful, but WITHOUT ANY WARRANTY; without even the implied
27 * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
28 * See the GNU General Public License for more details.
29 *
30 * You should have received a copy of the GNU General Public License
31 * along with the program; if not, write to the Free Software
32 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139 USA. Check
33 * the file named COPYING.
34 *)
35
36 MODEL td_model;
37 (* This MODEL is simply a root and documentation container for the
38 * thermodynamic modeling library.
39 * It should contain nothing except notes on the library and
40 * the boundwidth which is standard for this library.
41 *)
42 NOTES
43 'licensing' SELF {This library is subject to the GNU Public License v.2.0}
44
45 'sacredness' SELF {
46 There is nothing sacred or mysterious about either
47 thermodynamics in general or this implementation of
48 property calculations in particular. Please experiment --
49 you really cannot break anything in such a way that you
50 cannot also fix it.
51
52 It's all really just a matter of consistent and correct
53 modeling practice combined with good coefficient data
54 supplied by thermodynamicists.
55 }
56
57 'purpose' SELF {
58 Thermodynamic-based physical property package for
59 expressing the state information about multiphase,
60 multicomponent streams and holdups. Handles single-phase,
61 single-component degenerate cases correctly.
62 }
63
64 'author' SELF {Arthur W. Westerberg, Ben Allan, and Vicente Rico-Ramirez}
65 'created' SELF {November 9, 1997}
66 'last revised' SELF {$Revision: 1.6 $}
67
68 'new users' SELF {
69 The enduser MODEL in this library, if there is one, is
70 "thermodynamics" which has a relatively simple interface.
71 Adding new correlations for properties is relatively easy. To
72 understand how, examine by way of example the two models
73 "ideal_vapor_mixture" and "Pitzer_vapor_mixture".
74 }
75
76 'telephone support' SELF {
77 You must be kidding. CAPD consortium sponsors can contact the authors
78 for further consultancy. Everyone else can submit questions, suggestions,
79 and bug reports via the World Wide Web.
80 }
81
82 'URL' SELF {http://www.cs.cmu.edu/~ascend}
83
84 'design goal' SELF {
85 The approach used by this form for the physical property package is
86 to have the modeler using it define only the components (see model
87 "components_data") in mixture and the phases and the choice of
88 mixture models (see model "phases_data") to be used for the vapor
89 and liquid phases present. The models contained herein will set up
90 all the needed pure component data models, the physical property
91 models and the equilibrium relationships from just this
92 information.
93
94 New mixing rules should be easily added (and are) by defining new
95 refinements of the phase_partials MODEL here and adding them
96 to the alternatives allowed in phases.a4l.
97 }
98
99 'history' SELF {
100 Contains previously created models for a pure Pitzer vapor
101 component and a pure Rackett liquid component, and for a Pitzer
102 vapor mixture, a UNIFAC liquid mixture and a Wilson liquid mixture
103 created by Joe Zaher, Bob Huss and Boyd Safrit.
104 New models for ideal gases added by Ben Allan.
105 }
106
107 'technique' SELF {
108 The models in this file follow the pattern for easy reuse of:
109 Given (MODEL name and parameter list are the givens)
110 Assumptions (These are enforced as needed by the WHERE statements)
111 Wanted (NOTES 'purpose' SELF are explain what the MODEL is to do)
112 Equations and variables.
113 METHODS
114 Starting default values for a typical problem
115 Checking functions that verify the solutions calculated if requested.
116 Bounding functions that limit the variables for nonlinear search.
117 Scaling functions that provide scaling values for variables and some eqns.
118 Specify and reset functions for typical degree of freedom configurations.
119 }
120
121 'naming conventions' SELF {
122 Note the following variable naming convention applies here.
123 If a variable is the simple sum of simple products,
124 such as Vs = Sum Vi * yi, then it is notated v_y,
125 the inner product of vector V dot vector y.
126 The same convention applies in streams and holdup models.
127 This convention does not generalize prettily to sums of
128 vector products.
129 For example Sum(phi[j]*phase[j].v_y) would be rather
130 ambiguously notated as phi_v_y. So there it gets
131 simplified to a capital letter V. We need superscript
132 underbars or curly underbars or some such.
133 }
134 END NOTES;
135
136 boundwidth "width factor for bound calculations" IS_A bound_width;
137
138 END td_model;
139
140 (* murphree_equilibrium_mixture added by FIXME!! fixme! *)
141
142 (* ****************************************************************** *)
143
144 MODEL thermodynamic_properties(
145 P WILL_BE pressure;
146 T WILL_BE temperature;
147 ) REFINES td_model;
148
149 v "species molar volume" IS_A molar_volume;
150 h "species molar enthalpy" IS_A molar_energy;
151 g "species molar Gibbs free energy" IS_A molar_energy;
152
153
154 NOTES
155 'purpose' SELF {
156 The base thermodynamic property model. Two other models in this
157 library are refinements of this model. Its purpose is to
158 create the variables for pressure, temperature, molar volume, molar
159 enthalpy and molar Gibbs free energy of a pure species.
160 }
161 END NOTES;
162
163 METHODS
164
165 (* inherits ClearAll, reset *)
166
167 METHOD default_all;
168 (* all the ATOM default values are fine so far as we know *)
169 RUN default_self;
170 END default_all;
171
172 METHOD default_self;
173 (* All the ATOM default values are fine so far as we know.
174 * Might want to change the boundwidth value, though, later..
175 *)
176 END default_self;
177
178 METHOD check_all;
179 IF (T < 0.0{K}) THEN
180 STOP {Negative absolute temperature?!};
181 END IF;
182 IF (P < 0.0{Pa}) THEN
183 STOP {Negative absolute pressure?!};
184 END IF;
185 RUN check_self;
186 END check_all;
187
188 METHOD check_self;
189 IF (v <= 0.0{liter/mole}) THEN
190 STOP {Negative molar volume?};
191 END IF;
192 (* Don't know how to check g,h *)
193 END check_self;
194
195 METHOD scale_all;
196 (* assuming all are not near 0.0, which is not necessarily so. *)
197 T.nominal := T;
198 P.nominal := P;
199 RUN scale_self;
200 END scale_all;
201
202 METHOD scale_self;
203 (* assuming all are not near 0.0, which is not necessarily so. *)
204 v.nominal := v;
205 h.nominal := abs(h);
206 g.nominal := abs(g);
207 END scale_self;
208
209 METHOD bound_all;
210 (* resets most bounds to the interval centered at the current point.
211 * (variable.value) +/- boundwidth * variable.nominal
212 *)
213 T.lower_bound := 1.0e-8{K};
214 P.lower_bound := 1.0e-8{atm};
215 T.upper_bound := T + boundwidth*T.nominal;
216 P.upper_bound := P + boundwidth*P.nominal;
217 RUN bound_self;
218 END bound_all;
219
220 METHOD bound_self;
221 v.upper_bound := v + boundwidth*v.nominal;
222 v.lower_bound := 1.0e-8{liter/g_mole};
223 h.upper_bound := h + boundwidth*h.nominal;
224 h.lower_bound := h - boundwidth*h.nominal;
225 g.upper_bound := g + boundwidth*g.nominal;
226 g.lower_bound := g - boundwidth*g.nominal;
227 END bound_self;
228
229 METHOD specify;
230 FIX T;
231 FIX P;
232 FIX v;
233 FIX h;
234 FIX g;
235 END specify;
236
237 END thermodynamic_properties;
238
239 (* ****************************************************************** *)
240
241
242 MODEL pure_component(
243 P WILL_BE pressure;
244 T WILL_BE temperature;
245 data WILL_BE td_component_constants;
246 ) REFINES thermodynamic_properties;
247
248 NOTES
249 'purpose' SELF {
250 This model provides the pressure, temperature, molar volume, molar
251 enthalpy and molar Gibbs free energy for a pure component.
252 }
253 END NOTES;
254
255 METHODS
256 (* Inherits: ClearAll, default_self, bound_self,
257 * check_self, scale_self, specify, reset.
258 *)
259 METHOD check_all;
260 RUN data.check_all;
261 RUN check_self;
262 END check_all;
263 METHOD default_all;
264 RUN data.default_all;
265 RUN default_self;
266 END default_all;
267 METHOD scale_all;
268 RUN data.scale_all;
269 RUN scale_self;
270 END scale_all;
271 METHOD bound_all;
272 RUN data.bound_all;
273 RUN bound_self;
274 END bound_all;
275
276 END pure_component;
277
278 (* ****************************************************************** *)
279
280
281 MODEL partial_component(
282 P WILL_BE pressure;
283 T WILL_BE temperature;
284 ) REFINES thermodynamic_properties;
285
286 NOTES
287 'purpose' SELF {
288 This model provides the pressure, temperature, partial
289 molar volume, partial molar enthalpy and partial molar
290 Gibbs free energy for a single component.
291 }
292 END NOTES;
293
294 METHODS
295 (* Inherits: ClearAll, default_all, default_self, bound_all, bound_self,
296 * check_all, check_self, scale_all, scale_self, specify, reset.
297 *)
298 METHOD specify_partials_computed;
299 (* assumes equations for v,h,g have been added by refinement. *)
300 FIX T;
301 FIX P;
302 END specify_partials_computed;
303
304 END partial_component;
305
306 (* ****************************************************************** *)
307
308
309 MODEL Pitzer_vapor_component(
310 P WILL_BE pressure;
311 T WILL_BE temperature;
312 data WILL_BE td_component_constants;
313 ) REFINES pure_component;
314
315 NOTES
316 'purpose' SELF {
317 This model provides the equations that compute molar
318 volume, molar enthalpy and molar Gibbs free energy in terms
319 of pressure and temperature for a single component based on
320 the Pitzer vapor model.
321 }
322 END NOTES;
323 Pitzer_component_v:
324 P*v/(1{GAS_C})/data.Tc = T/data.Tc + (P/data.Pc)*
325 (0.083 - 0.422*(data.Tc/T)^1.6 + data.omega*
326 (0.139 - 0.172*(data.Tc/T)^4.2));
327
328 Pitzer_component_h:
329 h/(1{GAS_C})/data.Tc = data.H0/(1{GAS_C})/data.Tc +
330 data.cpvapa*(T - data.T0)/(1{GAS_C})/data.Tc +
331 data.cpvapb*(T^2 - data.T0^2)/2/(1{GAS_C})/data.Tc +
332 data.cpvapc*(T^3 - data.T0^3)/3/(1{GAS_C})/data.Tc +
333 data.cpvapd*(T^4 - data.T0^4)/4/(1{GAS_C})/data.Tc +
334 (P/data.Pc)*(0.083 - 1.097*(data.Tc/T)^1.6 + data.omega*
335 (0.139 - 0.894*(data.Tc/T)^4.2));
336
337 Pitzer_component_g:
338 g/(1{GAS_C})/data.Tc = data.G0/(1{GAS_C})/data.Tc -
339 (data.H0 - data.G0)*(T/data.T0 - 1)/(1{GAS_C})/data.Tc -
340 data.cpvapa*(T*lnm(T/data.T0) - T + data.T0)/1{GAS_C}/data.Tc -
341 data.cpvapb*(T^2 - 2*T*data.T0 + data.T0^2)/
342 (2*1{GAS_C}*data.Tc) -
343 data.cpvapc*(T^3/2 - 3*T*data.T0^2/2 + data.T0^3)/
344 (3*1{GAS_C}*data.Tc) -
345 data.cpvapd*(T^4/3 - 4*T*data.T0^3/3 + data.T0^4)/
346 (4*1{GAS_C}*data.Tc) +
347 T*lnm(P/data.P0)/data.Tc +
348 (P/data.Pc)*
349 (0.083 - 0.422*(data.Tc/T)^1.6 + data.omega*
350 (0.139 - 0.172*(data.Tc/T)^4.2));
351
352 METHODS
353 (* Inherited methods: ClearAll, check_all, check_self,
354 * bound_all, bound_self, scale_all, scale_self, reset.
355 *)
356 METHOD default_all;
357 (* Rather arbitrary but ok limits for defaults.
358 * the default bounds on T, P are rather too large.
359 *)
360 T.lower_bound := 1.0{K};
361 P.lower_bound := 0.1{atm};
362 T.upper_bound := 1000{K};
363 P.upper_bound := 20{atm};
364 RUN data.default_all;
365 RUN default_self;
366 END default_all;
367
368 METHOD default_self;
369 v := 24{liter/mole};
370 (* ATOM defaults ok for rest. *)
371 END default_self;
372
373 METHOD specify;
374 (* release v,h,g since equations calculate them. *)
375 FIX T;
376 FIX P;
377 END specify;
378
379 END Pitzer_vapor_component;
380
381 (* ****************************************************************** *)
382 MODEL ideal_vapor_component(
383 P WILL_BE pressure;
384 T WILL_BE temperature;
385 data WILL_BE td_component_constants;
386 ) REFINES pure_component;
387
388 NOTES
389 'purpose' SELF {
390 This model provides the equations that compute molar
391 volume, molar enthalpy and molar Gibbs free energy in terms
392 of pressure and temperature for a single component based on
393 the ideal gas (vapor) model.
394 }
395 'FIXME' SELF {Needs to be independently checked. Did I get h and g correct?}
396 END NOTES;
397
398 ideal_component_v: P*v/(1{GAS_C})/data.Tc = T/data.Tc;
399
400 ideal_component_h:
401 h/(1{GAS_C})/data.Tc =
402 data.H0/(1{GAS_C})/data.Tc +
403 data.cpvapa*(T - data.T0)/(1{GAS_C})/data.Tc +
404 data.cpvapb*(T^2 - data.T0^2)/2/(1{GAS_C})/data.Tc +
405 data.cpvapc*(T^3 - data.T0^3)/3/(1{GAS_C})/data.Tc +
406 data.cpvapd*(T^4 - data.T0^4)/4/(1{GAS_C})/data.Tc;
407
408 ideal_component_g:
409 g/(1{GAS_C})/data.Tc =
410 data.G0/(1{GAS_C})/data.Tc -
411 (data.H0 - data.G0)*(T/data.T0 - 1)/(1{GAS_C})/data.Tc -
412 data.cpvapa*(T*lnm(T/data.T0) - T + data.T0)/1{GAS_C}/data.Tc -
413 data.cpvapb*(T^2 - 2*T*data.T0 + data.T0^2)/(2*1{GAS_C}*data.Tc) -
414 data.cpvapc*(T^3/2 - 3*T*data.T0^2/2 + data.T0^3)/
415 (3*1{GAS_C}*data.Tc) -
416 data.cpvapd*(T^4/3 - 4*T*data.T0^3/3 + data.T0^4)/
417 (4*1{GAS_C}*data.Tc) +
418 T*lnm(P/data.P0)/data.Tc;
419
420 METHODS
421 (* Inherited methods: ClearAll, check_all, check_self,
422 * bound_all, bound_self, scale_all, scale_self, reset.
423 *)
424 METHOD default_all;
425 (* Rather arbitrary but ok limits for defaults.
426 * the default bounds on T, P are rather too large.
427 *)
428 T.lower_bound := 1.0{K};
429 P.lower_bound := 0.1{atm};
430 T.upper_bound := 1000{K};
431 P.upper_bound := 10{atm};
432 RUN data.default_all;
433 RUN default_self;
434 END default_all;
435 METHOD default_self;
436 v := 24{liter/mole};
437 (* ATOM defaults ok for rest. *)
438 END default_self;
439 METHOD specify;
440 (* release v,h,g since equations calculate them. *)
441 FIX T;
442 FIX P;
443 END specify;
444 END ideal_vapor_component;
445
446 (* ****************************************************************** *)
447
448
449 MODEL Rackett_liquid_component(
450 P WILL_BE pressure;
451 T WILL_BE temperature;
452 data WILL_BE td_component_constants;
453 ) REFINES pure_component;
454
455 NOTES
456 'purpose' SELF {
457 This model provides the equations that compute vapor
458 pressure, molar volume, molar enthalpy and molar Gibbs free
459 energy in terms of pressure and temperature for a pure
460 component based on the Rackett liquid model.
461 }
462 END NOTES;
463
464 VP "pure species vapor pressure" IS_A pressure;
465
466 (*Rackett_component_VP:*)
467 SELECT (data.vp_correlation)
468 CASE 1:
469 lnm(VP/data.Pc)*T/data.Tc*abs(data.Tc - T) =
470 (data.vpa * abs(1.0 - T/data.Tc) +
471 data.vpb * abs(1.0 - T/data.Tc)^1.5 +
472 data.vpc * abs(1.0 - T/data.Tc)^3.0 +
473 data.vpd * abs(1.0 - T/data.Tc)^6.0) * (data.Tc - T);
474 CASE 2:
475 lnm(VP/1{bar}) / lnm(data.Pc/1{bar}) =
476 (data.vpa - data.vpb/T + data.vpc * lnm(T/1{K}) +
477 data.vpd * VP/sqr(T)) / lnm(data.Pc/1{bar});
478 CASE 3:
479 lnm(VP/1{bar}) / lnm(data.Pc/1{bar}) =
480 (data.vpa - data.vpb/(data.vpc + T)) / lnm(data.Pc/1{bar});
481 END SELECT;
482
483 Rackett_component_v: v/data.Vliq =
484 data.Zc^(abs(1.0 - T/data.Tc)^(2/7))/
485 data.Zc^(abs(1.0 - data.Tliq/data.Tc)^(2/7));
486
487 Rackett_component_h:
488 h / (1{GAS_C} * data.Tc) = data.H0/(1{GAS_C})/data.Tc +
489 (T - data.T0) * (data.cpvapa/(1{GAS_C})/data.Tc) +
490 (T^2 - data.T0^2) * (data.cpvapb/2/(1{GAS_C})/data.Tc) +
491 (T^3 - data.T0^3) * (data.cpvapc/3/(1{GAS_C})/data.Tc) +
492 (T^4 - data.T0^4) * (data.cpvapd/4/(1{GAS_C})/data.Tc) +
493 (VP/data.Pc)*
494 (0.083 - 1.097*(data.Tc/T)^1.6 + data.omega*
495 (0.139 - 0.894*(data.Tc/T)^4.2)) -
496 (data.Hv/(1{GAS_C})/data.Tc)*
497 abs((data.Tc-T)/(data.Tc-data.Tb))^0.38 +
498 (P - VP)*(data.Vliq/(1{GAS_C})/data.Tc)*
499 (data.Zc^(abs(1.0 - T/data.Tc)^(2/7))/
500 data.Zc^(abs(1.0 - data.Tliq/data.Tc)^(2/7)))*(1.0 -
501 (-2/7)*(T/data.Tc)*(abs(1 - T/data.Tc)^(-5/7))*lnm(data.Zc));
502
503 Rackett_component_g: g/(1{GAS_C}*data.Tc) =
504 data.G0/(1{GAS_C})/data.Tc -
505 (T/data.T0 - 1) * ((data.H0 - data.G0)/(1{GAS_C})/data.Tc) -
506 (T*lnm(T/data.T0) - T + data.T0) *
507 (data.cpvapa/(1{GAS_C})/data.Tc) -
508 (T^2 - T*(2*data.T0) + data.T0^2) *
509 (data.cpvapb/2/(1{GAS_C})/data.Tc) -
510 (T^3/2 - T*(3/2*data.T0^2) + data.T0^3) *
511 (data.cpvapc/3/(1{GAS_C})/data.Tc) -
512 (T^4/3 - T*(4/3*data.T0^3) + data.T0^4) *
513 (data.cpvapd/4/(1{GAS_C})/data.Tc) +
514 T*lnm(VP/data.P0)/data.Tc +
515 (VP/data.Pc) *
516 (0.083 - 0.422*(data.Tc/T)^1.6 + data.omega *
517 (0.139 - 0.172*(data.Tc/T)^4.2)) +
518 (P - VP)*(data.Vliq/(1{GAS_C})/data.Tc) *
519 (data.Zc^(abs(1.0 - T/data.Tc)^(2/7))/
520 data.Zc^(abs(1.0 - data.Tliq/data.Tc)^(2/7)));
521
522 METHODS
523 (* Inherited methods: ClearAll, check_all, check_self, reset. *)
524 (* we need to see if there are any qualitative checks of VP
525 * or the VP slope
526 * to be added in check_self. At present, we just
527 * inherit checking of v. fixme, FIXME.
528 *)
529 METHOD default_all;
530 T.lower_bound := 1.0{K};
531 P.lower_bound := 0.1{atm};
532 T.upper_bound := 1000{K};
533 P.upper_bound := 20{atm};
534 RUN data.default_all;
535 RUN default_self;
536 END default_all;
537
538 METHOD default_self;
539 VP.lower_bound := 1.0e-12{Pa};
540 v := 0.1{liter/mole};
541 END default_self;
542
543 METHOD scale_all;
544 RUN data.scale_all;
545 RUN scale_self;
546 END scale_all;
547
548 METHOD scale_self;
549 VP.nominal := VP;
550 RUN pure_component::scale_self;
551 END scale_self;
552
553 METHOD bound_all;
554 RUN data.bound_all;
555 RUN bound_self;
556 END bound_all;
557
558 METHOD bound_self;
559 VP.lower_bound := VP - boundwidth*VP.nominal;
560 IF (VP.lower_bound < 1.0e-12{Pa}) THEN
561 VP.lower_bound := 1.0e-12{Pa};
562 END IF;
563 VP.upper_bound := VP + boundwidth*VP.nominal;
564 RUN pure_component::bound_self;
565 END bound_self;
566
567 METHOD specify;
568 FIX T;
569 FIX P;
570 END specify;
571
572 END Rackett_liquid_component;
573
574 (* ****************************************************************** *)
575
576 MODEL phase_partials(
577 cd WILL_BE components_data;
578 ) REFINES td_model;
579
580 NOTES
581 'purpose' SELF {
582 This model is the base model for computing the mixture
583 thermodynamic properties of molar volume, molar enthalpy and molar
584 Gibbs free energy in terms of pressure and temperature for a single
585 phase, based on summing the partial molar quantities times the
586 corresponding mole fractions for the components in the mixture.
587 This base model also provides space for relative volatilities for
588 the phase -- which can be used later to computer phase equilibrium
589 for multi-phase mixtures.
590 }
591 END NOTES;
592 P IS_A pressure;
593 T IS_A temperature;
594
595 v_y "mixture molar volume" IS_A molar_volume;
596 h_y "mixture molar enthalpy" IS_A molar_energy;
597 g_y "mixture molar Gibbs free energy" IS_A molar_energy;
598
599 components ALIASES cd.components;
600 partial[components] IS_A partial_component(P, T);
601 alpha[components] IS_A relative_volatility;
602
603 NOTES 'scaling' v_mix,h_mix,g_mix {
604 DIMENSIONLESS well-scaled equations are easier to converge from anywhere,
605 even if simply bilinear in dimensional form.
606 v*Pcref/(R*Tcref) = Sum(yi*vi |i) * Pcref/(R*Tcref);
607 comes from v = Sum(yi*vi);
608 } END NOTES;
609
610 v_mix: v_y * (cd.data[cd.reference].Pc /
611 (cd.data[cd.reference].Tc * 1{GAS_C}))
612 = SUM[y[i]*partial[i].v | i IN cd.components] *
613 (cd.data[cd.reference].Pc /
614 (cd.data[cd.reference].Tc * 1{GAS_C}));
615
616 (* likewise for h and g *)
617 h_mix: h_y / (cd.data[cd.reference].Tc * 1{GAS_C})
618 = SUM[y[i]*partial[i].h | i IN cd.components] /
619 (cd.data[cd.reference].Tc * 1{GAS_C});
620
621 g_mix: g_y / (cd.data[cd.reference].Tc * 1{GAS_C})
622 = SUM[y[i]*partial[i].g | i IN cd.components] /
623 (cd.data[cd.reference].Tc * 1{GAS_C});
624
625 y[components] "mixture species mole fractions" IS_A mole_fraction;
626
627 NOTES 'slack explanation' slack_PhaseDisappearance {
628 When slack_PhaseDisappearance != 0.0, then the phase
629 does not exist in this formulation of multi-phase
630 equilibrium thermodynamics.
631 } END NOTES;
632 slack_PhaseDisappearance "complementarity slack variable" IS_A fraction;
633
634 sum_y: SUM[y[components]] + slack_PhaseDisappearance = 1.0;
635
636 METHODS
637
638 (* Inherited methods: ClearAll, reset. *)
639 METHOD default_self;
640 T.lower_bound := 1.0{K};
641 P.lower_bound := 1.0{Pa};
642 v_y := 24{liter/mole}; (* presupposes a vapor *)
643 RUN partial[components].default_self;
644 END default_self;
645
646 METHOD default_all;
647 RUN cd.default_all;
648 RUN default_self;
649 END default_all;
650
651 METHOD check_all;
652 RUN cd.check_all;
653 IF (T < 0.0{K}) THEN
654 STOP {Phase has negative absolute temperature?!};
655 END IF;
656 IF (P < 0.0{Pa}) THEN
657 STOP {Phase has negative absolute pressure?!};
658 END IF;
659 RUN check_self;
660 END check_all;
661
662 METHOD check_self;
663 IF (v_y <= 0.0{liter/mole}) THEN
664 STOP {Phase has negative molar volume?!};
665 END IF;
666 (* Don't know how to check g_y,h_y *)
667 RUN partial[components].check_self;
668 END check_self;
669
670 METHOD scale_all;
671 RUN cd.scale_all;
672 RUN scale_self;
673 END scale_all;
674
675 METHOD scale_self;
676 (* assuming all are not near 0.0, which is not necessarily so. *)
677 T.nominal := T;
678 P.nominal := P;
679 v_y.nominal := v_y;
680 h_y.nominal := abs(h_y);
681 g_y.nominal := abs(g_y);
682 RUN partial[components].scale_self;
683 END scale_self;
684
685 METHOD bound_all;
686 RUN cd.bound_all;
687 RUN bound_self;
688 END bound_all;
689
690 METHOD bound_self;
691 (* Resets most bounds to the interval centered at the current point.
692 * (variable.value) +/- boundwidth * variable.nominal
693 *)
694 T.upper_bound := T + boundwidth*T.nominal;
695 T.lower_bound := 1.0e-8{K};
696 P.upper_bound := P + boundwidth*P.nominal;
697 P.lower_bound := 1.0e-8{atm};
698 v_y.upper_bound := v_y + boundwidth*v_y.nominal;
699 v_y.lower_bound := 1.0e-8{liter/g_mole};
700 h_y.upper_bound := h_y + boundwidth*h_y.nominal;
701 h_y.lower_bound := h_y - boundwidth*h_y.nominal;
702 g_y.upper_bound := g_y + boundwidth*g_y.nominal;
703 g_y.lower_bound := g_y - boundwidth*g_y.nominal;
704 END bound_self;
705
706 METHOD specify_partials_seqmod_massonly;
707 (* like specify_partials_computed, except leaves alpha[i] free *)
708 FIX P;
709 FIX T;
710 FREE alpha[components];
711 RUN partial[components].specify;
712 FIX y[cd.other_components];
713 FIX slack_PhaseDisappearance;
714 slack_PhaseDisappearance := 0.0;
715 END specify_partials_seqmod_massonly;
716
717 METHOD specify_partials_seqmod;
718 (* like specify_partials_computed, except leaves alpha[i] free *)
719 FIX P;
720 FIX T;
721 FREE alpha[components];
722 RUN partial[components].specify_partials_computed;
723 FIX y[cd.other_components];
724 FIX slack_PhaseDisappearance;
725 slack_PhaseDisappearance := 0.0;
726 END specify_partials_seqmod;
727
728 METHOD specify_partials_computed;
729 FIX P;
730 FIX T;
731 FIX alpha[components];
732 RUN partial[components].specify_partials_computed;
733 FIX y[cd.other_components];
734 FIX slack_PhaseDisappearance;
735 slack_PhaseDisappearance := 0.0;
736 END specify_partials_computed;
737
738 METHOD specify;
739 FIX P;
740 FIX T;
741 FIX alpha[components];
742 RUN partial[components].specify;
743 FIX y[cd.other_components];
744 FIX slack_PhaseDisappearance;
745 slack_PhaseDisappearance := 0.0;
746 END specify;
747
748 END phase_partials;
749
750 (* ****************************************************************** *)
751
752 MODEL ideal_vapor_mixture(
753 cd WILL_BE components_data;
754 ) REFINES phase_partials();
755
756 NOTES 'purpose' SELF {
757 A refinement of the phase_partials model given just above, this
758 model computes the partial molar quantities for a vapor mixture
759 phase for molar volume, molar enthalpy and molar Gibbs free energy
760 based on the ideal gas model for a mixture.
761 } END NOTES;
762
763 data ALIASES cd.data;
764
765 FOR i IN components CREATE
766 pure[i] "pure species ideal gas properties for the i-th species"
767 IS_A ideal_vapor_component(P, T, data[i]);
768 END FOR;
769
770 FOR i IN components CREATE
771
772 ideal_mixture_v[i]:
773 (partial[i].v - pure[i].v) *
774 (data[i].Pc/(1{GAS_C}*data[i].Tc)) = 0;
775
776 ideal_mixture_h[i]:
777 (partial[i].h - pure[i].h)/((1{GAS_C})*data[i].Tc) = 0;
778
779 ideal_mixture_g[i]:
780 (partial[i].g - pure[i].g - (1{GAS_C})*T*lnm(y[i])) /
781 (1{GAS_C}*data[i].Tc) = 0;
782 END FOR;
783
784 METHODS
785
786 (* Inherited methods: ClearAll, default_all, check_all, bound_all,
787 * scale_all, reset.
788 *)
789 METHOD default_self;
790 RUN phase_partials::default_self;
791 RUN pure[components].default_self;
792 END default_self;
793
794 METHOD check_self;
795 RUN phase_partials::check_self;
796 RUN pure[components].check_self;
797 END check_self;
798
799 METHOD bound_self;
800 RUN phase_partials::bound_self;
801 RUN pure[components].bound_self;
802 END bound_self;
803
804 METHOD scale_self;
805 RUN phase_partials::scale_self;
806 RUN pure[components].scale_self;
807 END scale_self;
808
809 METHOD specify;
810 RUN phase_partials::specify_partials_computed;
811 END specify;
812
813 END ideal_vapor_mixture;
814
815 (* ****************************************************************** *)
816
817 MODEL Pitzer_vapor_mixture(
818 cd WILL_BE components_data;
819 ) REFINES phase_partials();
820
821 NOTES 'purpose' SELF {
822 A refinement of the phase_partials model given just above, this
823 model computes the partial molar quantities for a vapor mixture
824 phase for molar volume, molar enthalpy and molar Gibbs free energy
825 based on the Pitzer vapor model for a mixture.
826 }
827 END NOTES;
828
829 data ALIASES cd.data;
830
831 FOR i IN components CREATE
832 pure[i]
833 "pure species properties for the i-th species using acentric factors"
834 IS_A Pitzer_vapor_component(P, T, data[i]);
835 END FOR;
836
837 FOR i IN components CREATE
838
839 Pitzer_mixture_v[i]:
840 (partial[i].v - pure[i].v)*(data[i].Pc/
841 (1{GAS_C}*data[i].Tc)) =
842 -1.0*(0.083 - 0.422*(data[i].Tc/T)^1.6 + data[i].omega*
843 (0.139 - 0.172*(data[i].Tc/T)^4.2))*(1.0 - y[i])^2 +
844 0.50*(1.0 - y[i])*
845 SUM[y[j]*((1.0 + (data[j].Vc/data[i].Vc)^(1/3))^3/
846 (1.0 + data[j].Zc/data[i].Zc))*
847 (0.083 - 0.422*(sqrt(data[i].Tc*data[j].Tc)/T)^1.6 +
848 0.5*(data[i].omega + data[j].omega)*
849 (0.139 - 0.172*(sqrt(data[i].Tc*data[j].Tc)/T)^4.2))
850 | j IN components - [i]];
851
852 Pitzer_mixture_h[i]:
853 (partial[i].h - pure[i].h)/((1{GAS_C})*data[i].Tc) =
854 -(P/data[i].Pc)*
855 (0.083 - 1.097*(data[i].Tc/T)^1.6 + data[i].omega*
856 (0.139 - 0.894*(data[i].Tc/T)^4.2))*(1.0 - y[i])^2 +
857 (1.0 - y[i])*(P/(data[i].Pc*2))*
858 SUM[y[j]*((1.0 + (data[j].Vc/data[i].Vc)^(1/3))^3/
859 (1.0 + data[j].Zc/data[i].Zc))*
860 (0.083 - 1.097*(sqrt(data[i].Tc*data[j].Tc)/T)^1.6 +
861 0.5*(data[i].omega + data[j].omega)*
862 (0.139 - 0.894*(sqrt(data[i].Tc*data[j].Tc)/T)^4.2))
863 | j IN components - [i]];
864
865 Pitzer_mixture_g[i]:
866 (partial[i].g - pure[i].g - (1{GAS_C})*T*lnm(y[i])) /
867 (1{GAS_C}*data[i].Tc) =
868 -(P/data[i].Pc)*
869 (0.083 - 0.422*(data[i].Tc/T)^1.6 + data[i].omega*
870 (0.139 - 0.172*(data[i].Tc/T)^4.2))*(1.0 - y[i])^2 +
871 (1.0 - y[i])*(P/(data[i].Pc*2))*
872 SUM[y[j]*((1.0 + (data[j].Vc/data[i].Vc)^(1/3))^3/
873 (1.0 + data[j].Zc/data[i].Zc))*
874 (0.083 - 0.422*(sqrt(data[i].Tc*data[j].Tc)/T)^1.6 +
875 0.5*(data[i].omega + data[j].omega)*
876 (0.139 - 0.172*(sqrt(data[i].Tc*data[j].Tc)/T)^4.2))
877 | j IN components - [i]];
878 END FOR;
879
880 METHODS
881
882 (* Inherited methods: ClearAll, default_all, check_all, bound_all,
883 * scale_all, reset.
884 *)
885 METHOD default_self;
886 RUN phase_partials::default_self;
887 RUN pure[components].default_self;
888 END default_self;
889
890 METHOD check_self;
891 RUN phase_partials::check_self;
892 RUN pure[components].check_self;
893 END check_self;
894
895 METHOD bound_self;
896 RUN phase_partials::bound_self;
897 RUN pure[components].bound_self;
898 END bound_self;
899
900 METHOD scale_self;
901 RUN phase_partials::scale_self;
902 RUN pure[components].scale_self;
903 END scale_self;
904
905 METHOD specify;
906 RUN phase_partials::specify_partials_computed;
907 END specify;
908
909 END Pitzer_vapor_mixture;
910
911 (* ****************************************************************** *)
912
913
914 MODEL UNIFAC_liquid_mixture(
915 cd WILL_BE components_data;
916 ) REFINES phase_partials();
917
918 NOTES 'purpose' SELF {
919 A refinement of the phase_partials model given just above, this
920 model computes the partial molar quantities for a liquid mixture
921 phase for molar volume, molar enthalpy and molar Gibbs free energy
922 based on the UNIFAC model for a liquid mixture.
923 } END NOTES;
924
925 data ALIASES cd.data;
926
927 FOR i IN components CREATE
928 pure[i]
929 "pure species liquid properties for the i-th species by Rackett"
930 IS_A Rackett_liquid_component(P, T, data[i]);
931 END FOR;
932
933 subgroups IS_A set OF symbol_constant;
934 groups IS_A set OF integer_constant;
935 comps[subgroups] IS_A set OF symbol_constant;
936 rv[components] "computed molecular volumes, ri constant"
937 IS_A UNIFAC_size;
938 qs[components] "computed molecular areas, qi constant"
939 IS_A UNIFAC_size;
940 Jv[components] IS_A factor;
941 Ls[components] IS_A factor;
942 theta[subgroups] IS_A factor;
943 eta[subgroups] IS_A factor;
944 uc IS_A UNIFAC_constants;
945
946 subgroups :== UNION[data[i].subgroups | i IN components];
947 groups :== UNION[data[i].groups | i IN components];
948 FOR k IN subgroups CREATE
949 comps[k] :== [i IN components | k IN data[i].subgroups];
950 END FOR;
951 FOR k IN subgroups CREATE
952 UNIFAC_mixture_theta[k]:
953 theta[k] = uc.Q[k]*SUM[data[i].nu[k]*y[i] | i IN comps[k]];
954
955 UNIFAC_mixture_eta[k]:
956 eta[k] =
957 SUM[theta[m] | m IN subgroups*uc.sub[uc.group[k]]] +
958 SUM[SUM[theta[m]*exp(-uc.a[g][uc.group[k]]/T)
959 | m IN subgroups*uc.sub[g]]
960 | g IN groups - [uc.group[k]]];
961 END FOR;
962
963 FOR i IN components CREATE
964 rv[i] :== SUM[0, data[i].nu[k]*uc.R[k]
965 | k IN data[i].subgroups];
966 qs[i] :== SUM[0, data[i].nu[k]*uc.Q[k]
967 | k IN data[i].subgroups];
968 END FOR;
969 FOR i IN components CREATE
970 UNIFAC_mixture_rv[i]:
971 rv[i] = Jv[i]*SUM[rv[j]*y[j] | j IN components];
972
973 UNIFAC_mixture_qs[i]:
974 qs[i] = Ls[i]*SUM[qs[j]*y[j] | j IN components];
975 UNIFAC_partial_v[i]:
976 partial[i].v = pure[i].v;
977
978 UNIFAC_mixture_h_excess[i]:
979 (partial[i].h - pure[i].h)/(1{GAS_C})/data[i].Tc =
980 SUM[0,theta[k]*
981 SUM[0,SUM[0,theta[n]*
982 ((uc.a[g][uc.group[k]] -
983 uc.a[uc.group[n]][uc.group[k]])/data[i].Tc)*
984 exp(-(uc.a[g][uc.group[k]] +
985 uc.a[uc.group[n]][uc.group[k]])/T)*
986 SUM[0,data[i].nu[m]*uc.Q[m]
987 | m IN data[i].subgroups*uc.sub[g]]
988 | g IN data[i].groups - [uc.group[n]]]
989 | n IN subgroups]/eta[k]/eta[k]
990 | k IN subgroups] -
991 SUM[0,(data[i].nu[k]*uc.Q[k]/(
992 SUM[0,data[i].nu[m]*uc.Q[m]
993 | m IN data[i].subgroups*uc.sub[uc.group[k]]] +
994 SUM[0,SUM[0,data[i].nu[m]*uc.Q[m] *
995 exp(-uc.a[g][uc.group[k]]/T)
996 | m IN data[i].subgroups*uc.sub[g]]
997 | g IN data[i].groups - [uc.group[k]]]))*
998 SUM[0,SUM[0,theta[n]*
999 ((uc.a[g][uc.group[k]] -
1000 uc.a[uc.group[n]][uc.group[k]])/data[i].Tc)*
1001 exp(-(uc.a[g][uc.group[k]] +
1002 uc.a[uc.group[n]][uc.group[k]])/T)*
1003 SUM[0,data[i].nu[m]*uc.Q[m]
1004 | m IN data[i].subgroups*uc.sub[g]]
1005 | g IN data[i].groups - [uc.group[n]]]
1006 | n IN subgroups]/eta[k]
1007 | k IN data[i].subgroups];
1008
1009 UNIFAC_mixture_g_excess[i]:
1010 (partial[i].g - pure[i].g - (1{GAS_C})*T*lnm(y[i])) /
1011 (1{GAS_C}*data[i].Tc) =
1012 (1.0 - Jv[i] + lnm(Jv[i]) -
1013 5.0*qs[i]*(1.0 - Jv[i]/Ls[i] + lnm(Jv[i]/Ls[i])) +
1014 qs[i]*(1 - lnm(Ls[i])))*T/data[i].Tc -
1015 SUM[0,theta[k]*(
1016 SUM[0,data[i].nu[m]*uc.Q[m]
1017 | m IN data[i].subgroups*uc.sub[uc.group[k]]] +
1018 SUM[0,SUM[0,data[i].nu[m] * uc.Q[m] *
1019 exp(-uc.a[g][uc.group[k]]/T)
1020 | m IN data[i].subgroups*uc.sub[g]]
1021 | g IN data[i].groups - [uc.group[k]]])/eta[k]
1022 | k IN subgroups]*T/data[i].Tc +
1023 SUM[0,data[i].nu[k]*uc.Q[k]*lnm((
1024 SUM[0,data[i].nu[m]*uc.Q[m]
1025 | m IN data[i].subgroups*uc.sub[uc.group[k]]] +
1026 SUM[0,SUM[0, data[i].nu[m] * uc.Q[m] *
1027 exp(-uc.a[g][uc.group[k]]/T)
1028 | m IN data[i].subgroups*uc.sub[g]]
1029 | g IN data[i].groups - [uc.group[k]]])/eta[k])
1030 | k IN data[i].subgroups]*T/data[i].Tc;
1031 END FOR;
1032
1033 FOR i IN components CREATE
1034 rv[i] :== SUM[0,data[i].nu[k]*uc.R[k]
1035 | k IN data[i].subgroups];
1036 qs[i] :== SUM[0,data[i].nu[k]*uc.Q[k]
1037 | k IN data[i].subgroups];
1038 END FOR;
1039
1040 METHODS
1041
1042 (* inherited methods: ClearAll, default_all, check_all, bound_all,
1043 * scale_all, reset.
1044 *)
1045
1046 METHOD default_UNIFAC;
1047 NOTES 'purpose' SELF {
1048 This method should be run after partial and pure are defaulted.
1049 This takes care of the mess of local variables in UNIFAC
1050 and makes the molar volumes initial values to a liquid size.
1051 } END NOTES;
1052 Jv[components].lower_bound := 1.0e-8;
1053 Ls[components].lower_bound := 1.0e-8;
1054 theta[subgroups].lower_bound := 0.0;
1055 eta[subgroups].lower_bound := 0.0;
1056 (* liquid molar volumes are much smaller than vapor volumes *)
1057 v_y := 0.1{liter/mole};
1058 partial[components].v := 0.1{liter/mole};
1059 pure[components].v := 0.1{liter/mole};
1060 END default_UNIFAC;
1061
1062 METHOD default_self;
1063 RUN phase_partials::default_self;
1064 RUN pure[components].default_self;
1065 RUN default_UNIFAC;
1066 END default_self;
1067
1068 METHOD check_self;
1069 (* fixme FIXME. Are there limits, on eta,theta,Ls,Jv we need to check?*)
1070 RUN phase_partials::check_self;
1071 RUN pure[components].check_self;
1072 END check_self;
1073
1074 METHOD scale_self;
1075 RUN phase_partials::scale_self;
1076 RUN pure[components].scale_self;
1077 (* we may need to further investigate the scaling on these
1078 * when we approach small numbers for Ls, Jv.
1079 * fixme FIXME.
1080 *)
1081 FOR i IN components DO
1082 Ls[i].nominal := Ls[i];
1083 Jv[i].nominal := Jv[i];
1084 END FOR;
1085 FOR j IN subgroups DO
1086 theta[j].nominal := theta[j];
1087 eta[j].nominal := eta[j];
1088 END FOR;
1089 END scale_self;
1090
1091 METHOD bound_self;
1092 RUN phase_partials::bound_self;
1093 RUN pure[components].bound_self;
1094 (* we may need to further investigate the bounding on these
1095 * when we approach small numbers for Ls, Jv.
1096 * fixme FIXME.
1097 *)
1098 FOR i IN components DO
1099 Ls[i].upper_bound := Ls[i] + boundwidth*Ls[i].nominal;
1100 Jv[i].upper_bound := Jv[i] + boundwidth*Jv[i].nominal;
1101 END FOR;
1102 FOR j IN subgroups DO
1103 theta[j].upper_bound := theta[j] + boundwidth*theta[j].nominal;
1104 eta[j].upper_bound := eta[j] + boundwidth*eta[j].nominal;
1105 END FOR;
1106 END bound_self;
1107
1108 METHOD specify;
1109 RUN phase_partials::specify_partials_computed;
1110 END specify;
1111
1112 END UNIFAC_liquid_mixture;
1113
1114 (* ****************************************************************** *)
1115
1116 MODEL Wilson_liquid_mixture(
1117 cd WILL_BE components_data;
1118 ) REFINES phase_partials();
1119
1120 NOTES 'purpose' SELF {
1121 A refinement of the phase_partials model given just above, this
1122 model computes the partial molar quantities for a liquid mixture
1123 phase for molar volume, molar enthalpy and molar Gibbs free energy
1124 based on the Wilson model for a liquid mixture.
1125 } END NOTES;
1126
1127 data ALIASES cd.data;
1128
1129 FOR i IN components CREATE
1130 pure[i] IS_A Rackett_liquid_component(P, T, data[i]);
1131 END FOR;
1132
1133 lambda[components][components] "Wilson binary parameters"
1134 IS_A factor;
1135
1136 FOR i IN components CREATE
1137 FOR j IN components CREATE
1138 lambda[i][j] = (pure[j].v / pure[i].v) *
1139 exp(-data[i].del_ip[data[j].formula] /
1140 (1{GAS_C} * T));
1141 END FOR;
1142 END FOR;
1143
1144 FOR i IN components CREATE
1145 partial[i].v = pure[i].v; (* fix me. scaling. *)
1146
1147 Wilson_mixture_g_excess[i]:
1148 partial[i].g - pure[i].g - (1{GAS_C})*T*lnm(y[i]) =
1149 (1{GAS_C})*T*(-lnm( SUM[y[j]*lambda[i][j] | j IN components])
1150 + 1 -
1151 SUM[(y[k] * lambda[k][i]) / SUM[ y[j] * lambda[k][j]
1152 | j IN components] | k IN components]);
1153
1154 Wilson_mixture_h_excess[i]:
1155 partial[i].h - pure[i].h =
1156 (SUM[ y[j] * lambda[i][j] *
1157 data[i].del_ip[data[j].formula]
1158 | j IN components]) /
1159 (SUM[y[s]*lambda[i][s]
1160 | s IN components]) -
1161 (SUM[y[k]*lambda[k][i] | k IN components]*
1162 SUM[SUM[y[l]*lambda[m][l]*
1163 data[m].del_ip[data[l].formula]
1164 | l IN components]
1165 | m IN components]) / (
1166 SUM[SUM[y[n]*lambda[o][n]
1167 | n IN components]
1168 | o IN components])^2 +
1169 (SUM[y[p]*lambda[p][i]*
1170 data[p].del_ip[data[i].formula]
1171 | p IN components]) /
1172 (SUM[SUM[y[q]*lambda[r][q]
1173 | q IN components]
1174 | r IN components]);
1175 END FOR;
1176
1177 METHODS
1178
1179 (* inherited methods: ClearAll, default_all, check_all, bound_all,
1180 * scale_all, reset.
1181 *)
1182
1183 METHOD default_Wilson;
1184 NOTES 'purpose' SELF {
1185 This method should be run after partial and pure are defaulted.
1186 This takes care of the mess of local variables in Wilson
1187 and makes the molar volumes initial values to a liquid size.
1188 } END NOTES;
1189
1190 (* The following initial bounds on lambda will only be useful
1191 if meaninful bounds on pure[i].v have been set. *)
1192 FOR i IN components DO
1193 FOR j IN components DO
1194 IF (pure[i].v.upper_bound > 0{liter/mole}) THEN
1195 lambda[i][j].lower_bound :=
1196 (pure[j].v.lower_bound / pure[i].v.upper_bound) *
1197 exp(-data[i].del_ip[data[j].formula] /
1198 (1{GAS_C} * T));
1199 END IF;
1200 IF (pure[i].v.lower_bound > 0{liter/mole}) THEN
1201 lambda[i][j].upper_bound :=
1202 (pure[j].v.upper_bound / pure[i].v.lower_bound) *
1203 exp(-data[i].del_ip[data[j].formula] /
1204 (1{GAS_C} * T));
1205 END IF;
1206 END FOR;
1207 END FOR;
1208
1209 (* liquid molar volumes are much smaller than vapor volumes *)
1210 v_y := 0.1{liter/mole};
1211 partial[components].v := 0.1{liter/mole};
1212 pure[components].v := 0.1{liter/mole};
1213 END default_Wilson;
1214
1215 METHOD default_self;
1216 RUN phase_partials::default_self;
1217 RUN pure[components].default_self;
1218 RUN default_Wilson;
1219 END default_self;
1220
1221 METHOD check_self;
1222 (* fixme FIXME. Are there limits, on lambdaIJ we need to check?*)
1223 RUN phase_partials::check_self;
1224 RUN pure[components].check_self;
1225 END check_self;
1226
1227 METHOD scale_self;
1228 RUN phase_partials::scale_self;
1229 RUN pure[components].scale_self;
1230 (* we may need to further investigate the scaling on these
1231 * when we approach small numbers for lambdaIJ.
1232 * fixme FIXME.
1233 *)
1234 FOR i IN components DO
1235 FOR j IN components DO
1236 lambda[i][j].nominal := lambda[i][j];
1237 END FOR;
1238 END FOR;
1239 END scale_self;
1240
1241 METHOD bound_self;
1242 RUN phase_partials::bound_self;
1243 RUN pure[components].bound_self;
1244 FOR i IN components DO
1245 FOR j IN components DO
1246 IF (pure[i].v.upper_bound > 0{liter/mole}) THEN
1247 lambda[i][j].lower_bound :=
1248 (pure[j].v.lower_bound / pure[i].v.upper_bound) *
1249 exp(-data[i].del_ip[data[j].formula] /
1250 (1{GAS_C} * T));
1251 END IF;
1252 IF (pure[i].v.lower_bound > 0{liter/mole}) THEN
1253 lambda[i][j].upper_bound :=
1254 (pure[j].v.upper_bound / pure[i].v.lower_bound) *
1255 exp(-data[i].del_ip[data[j].formula] /
1256 (1{GAS_C} * T));
1257 END IF;
1258 END FOR;
1259 END FOR;
1260 END bound_self;
1261
1262 METHOD specify;
1263 RUN phase_partials::specify_partials_computed;
1264 END specify;
1265
1266 END Wilson_liquid_mixture;
1267
1268 (* ****************************************************************** *)
1269
1270
1271 MODEL thermodynamics(
1272 cd WILL_BE components_data;
1273 pd WILL_BE phases_data;
1274 phase[pd.phases] WILL_BE phase_partials;
1275 equilibrated WILL_BE boolean;
1276 ) REFINES td_model;
1277
1278 NOTES
1279 'purpose' SELF {
1280 This model is the target model for this library of models.
1281 It computes the state of a multiphase multicomponent
1282 mixture. The state of a multiphase multicomponent mixture
1283 is a sufficient set of intensive variables to completely
1284 characterize it when one assumes the phases are in
1285 thermodynamic equilibrium. Here a sufficient set of
1286 variables is: pressure, temperature, molar composition,
1287 molar volume, molar enthalpy and molar Gibbs free energy
1288 for each of the phases and overall for all the phases.
1289 }
1290
1291 'equilibrium' SELF {
1292 Equilibrium is expressed in this MODEL by equating
1293 temperature and pressure in all the phases and then by
1294 using either a relative volatility model (equilibrated :=
1295 FALSE) or by equating the partial molar Gibbs free energy
1296 (chemical potential) for each of the components in each of
1297 the phases (equilibrated := TRUE).
1298 }
1299 END NOTES;
1300
1301 P ALIASES phase[pd.reference_phase].P;
1302 T ALIASES phase[pd.reference_phase].T;
1303
1304 V "molar average volume over all species i, phases j" IS_A molar_volume;
1305 H "molar average enthalpy over all species i, phases j" IS_A molar_energy;
1306 G "molar average Gibbs free energy over all species i, phases j"
1307 IS_A molar_energy;
1308 y[cd.components] "average mole fractions over all phases j" IS_A fraction;
1309 phase_fraction[pd.phases] "division of moles into phases " IS_A fraction;
1310 alpha_bar[pd.other_phases] "mole fraction weighted sum of volatilities"
1311 IS_A relative_volatility;
1312 number_of_other_phases "Np - 1" IS_A integer_constant;
1313
1314 FOR j IN [pd.phases] CREATE
1315 slack_PhaseDisappearance[j] ALIASES phase[j].slack_PhaseDisappearance;
1316 (* If you set slack_PhaseDisappearance[j].fixed := TRUE,
1317 * (useful if you want to _force_ a phase to exist)
1318 * then you must also set complementarity[j].included := FALSE;
1319 * OTHERWISE the Jacobian matrix of this MODEL is singular.
1320 * The only sensible value of slack_PhaseDisappearance[i]
1321 * when it is fixed is 0.0.
1322 *)
1323 END FOR;
1324
1325 sum_phase_fractions: SUM[phase_fraction[pd.phases]] = 1.0;
1326 overall_v: V * (cd.data[cd.reference].Pc /
1327 (1{GAS_C} * cd.data[cd.reference].Tc))
1328 = SUM[phase_fraction[j]*phase[j].v_y | j IN pd.phases] *
1329 (cd.data[cd.reference].Pc /
1330 (1{GAS_C} * cd.data[cd.reference].Tc));
1331 overall_h: H / (1{GAS_C} * cd.data[cd.reference].Tc)
1332 = SUM[phase_fraction[j]*phase[j].h_y | j IN pd.phases] /
1333 (1{GAS_C} * cd.data[cd.reference].Tc);
1334 overall_g: G / (1{GAS_C} * cd.data[cd.reference].Tc)
1335 = SUM[phase_fraction[j]*phase[j].g_y | j IN pd.phases]/
1336 (1{GAS_C} * cd.data[cd.reference].Tc);
1337
1338 FOR i IN cd.components CREATE
1339 overall_y[i]: y[i] = SUM[phase_fraction[j]*phase[j].y[i] | j IN pd.phases];
1340 END FOR;
1341
1342 number_of_other_phases :== CARD[pd.other_phases];
1343
1344 SELECT (number_of_other_phases)
1345 CASE 0:
1346 (* Create Nothing *)
1347 OTHERWISE:
1348
1349 (* Create equilibrium relationships *)
1350
1351 (* equilibrium *)
1352
1353 (* pressure and temperature equilibrium.
1354 * this condition of equilibrium is not UNIVERSAL
1355 * to all equilibrium thermodynamic systems.
1356 * It certainly applies to non-electrolytic Vapor-Liquid
1357 * systems, however.
1358 *)
1359 FOR j IN [pd.other_phases] CREATE
1360 equateP[j]: phase[j].P = phase[pd.reference_phase].P;
1361 equateT[j]: phase[j].T = phase[pd.reference_phase].T;
1362 END FOR;
1363
1364 NOTES
1365 'usage' alpha_bar, alphaeqn, phase[pd.phases].alpha[components] {
1366 If using a relative volatility version for equilibrium, fix
1367 the relative volatilities alpha for all components in all
1368 phases (except perhaps the reference phase) and compute
1369 alpha_bar.
1370
1371 Else activate the equations that equate the chemical
1372 potential for each component in each phase to its value in
1373 the reference phase, fix alpha_bar and set its value to one
1374 and compute alphas for all the components.
1375
1376 Note, there is a vector of alphas in the reference phase
1377 that is unused.
1378 }
1379 'interpretation' alpha {
1380 When alpha_bar is fixed at 1.0, then the alpha[i] are essentially
1381 K-values. When alpha bar is computed from fixed alpha[i], the
1382 alpha[i] are essentially relative volatilities. The oldest of the
1383 authors believes this note to be entirely redundant, but the
1384 youngest has found it helps some people to relate the MODEL to
1385 their old unit operations textbooks.
1386 }
1387 END NOTES;
1388
1389 FOR j IN [pd.other_phases] CREATE
1390 FOR i IN cd.components CREATE
1391 alphaeqn[j][i]:
1392 alpha_bar[j] * phase[j].y[i]
1393 = phase[j].alpha[i] * phase[pd.reference_phase].y[i];
1394 END FOR;
1395 END FOR;
1396
1397 FOR j IN [pd.other_phases] CREATE
1398 FOR i IN cd.components CREATE
1399 td_eql[j][i]:
1400 phase[j].partial[i].g
1401 = phase[pd.reference_phase].partial[i].g;
1402 END FOR;
1403 END FOR;
1404
1405 WHEN (equilibrated)
1406 CASE TRUE:
1407 USE td_eql;
1408 OTHERWISE:
1409 (* relax free energy, but still insist T,P same over all *)
1410 END WHEN;
1411
1412 (* Complementarity formulation *)
1413
1414 FOR j IN [pd.phases] CREATE
1415 complementarity[j]:
1416 slack_PhaseDisappearance[j] * phase_fraction[j] = 0.0;
1417 END FOR;
1418
1419 END SELECT;
1420
1421 METHODS
1422
1423 (* inherited methods: ClearAll, reset. *)
1424 METHOD check_self;
1425 (* there really isn't much that can be checked here since
1426 * whether phases missing or not is an error we can't tell here.
1427 * STOP {MODEL "thermodynamics" needs a check method for complementarity.};
1428 *)
1429 END check_self;
1430 METHOD check_all;
1431 RUN cd.check_all;
1432 RUN pd.check_all;
1433 RUN phase[pd.phases].check_all;
1434 RUN check_self;
1435 END check_all;
1436 METHOD default_self;
1437 V := 30 {liter/mole};
1438 END default_self;
1439 METHOD default_all;
1440 RUN cd.default_all;
1441 RUN pd.default_all;
1442 RUN phase[pd.phases].default_all;
1443 equilibrated := FALSE;
1444 RUN default_self;
1445 END default_all;
1446 METHOD scale_self;
1447 (* fixme! assumed nonzero! always TRUE? *)
1448 V.nominal := V;
1449 G.nominal := abs(G);
1450 H.nominal := abs(H);
1451 (* should we rescale y, phase_fraction, alpha_bar here?? *)
1452 END scale_self;
1453 METHOD scale_all;
1454 RUN cd.scale_all;
1455 RUN pd.scale_all;
1456 RUN phase[pd.phases].scale_all;
1457 RUN scale_self;
1458 END scale_all;
1459
1460 METHOD bound_all;
1461 RUN cd.bound_all;
1462 RUN pd.bound_all;
1463 RUN phase[pd.phases].bound_all;
1464 RUN bound_self;
1465 END bound_all;
1466
1467 METHOD bound_self;
1468 (* bounds on alpha_bar? y, phase_fraction? *)
1469 V.upper_bound := V + boundwidth*V.nominal;
1470 V.lower_bound := 1.0e-8{liter/g_mole};
1471 H.upper_bound := H + boundwidth*H.nominal;
1472 H.lower_bound := H - boundwidth*H.nominal;
1473 G.upper_bound := G + boundwidth*G.nominal;
1474 G.lower_bound := G - boundwidth*G.nominal;
1475 END bound_self;
1476
1477 METHOD specify;
1478 IF (pd.phase_indicator != 'M') THEN
1479 RUN phase[pd.phases].specify_partials_seqmod;
1480 ELSE
1481 RUN phase[pd.phases].specify_partials_seqmod_massonly;
1482 END IF;
1483 (* The phase specify methods fix
1484 * all but one y variable per phase, and each phase's P and T.
1485 * They also specify that the phases must not disappear.
1486 * They leave the alpha[i] free.
1487 *)
1488 FREE phase[pd.other_phases].P;
1489 FREE phase[pd.other_phases].T;
1490 FREE phase[pd.phases].y[cd.components];
1491 (* Now only reference phase P,T fixed, unless a mass only MODEL. *)
1492
1493 IF (number_of_other_phases != 0) THEN
1494 (* allow phases to disappear *)
1495 FREE slack_PhaseDisappearance[pd.phases];
1496 slack_PhaseDisappearance[pd.phases] := 0.0;
1497 END IF;
1498
1499 (* Fix Np - 1 phase fractions and Nc - 1 overall mole fractions *)
1500 FIX y[cd.other_components];
1501
1502 IF (equilibrated) THEN
1503 (* use thermo instead of relative volatility. free alphas. *)
1504 FREE phase[pd.phases].alpha[cd.components];
1505 FIX alpha_bar[pd.other_phases];
1506 alpha_bar[pd.other_phases] := 1.0;
1507 (* When alphas are free, they are essentially K-values *)
1508 (* Solve for fixed phase fraction for equilibrium
1509 * flash -- setting T may push flash out of multiphase region.
1510 * Do we need to swap T,ph_fr here? Think so.
1511 *)
1512 ELSE
1513 (* Else solve as relative volatility model *)
1514 FIX phase_fraction[pd.other_phases];
1515 FIX phase[pd.other_phases].alpha[cd.components];
1516 END IF;
1517 END specify;
1518
1519 METHOD reset_to_massbalance;
1520 RUN ClearAll;
1521 equilibrated := FALSE;
1522 RUN specify;
1523 END reset_to_massbalance;
1524
1525 METHOD reset_to_fullthermo;
1526 RUN ClearAll;
1527 equilibrated := TRUE;
1528 RUN specify;
1529 END reset_to_fullthermo;
1530
1531 END thermodynamics;

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