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

Contents of /trunk/models/thermodynamics.a4l

Parent Directory Parent Directory | Revision Log Revision Log


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

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