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

Diff of /trunk/models/thermodynamics.a4l

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2065 by johnpye, Tue May 9 03:42:08 2006 UTC revision 2066 by jpye, Tue Aug 11 03:50:20 2009 UTC
# Line 1  Line 1 
1  REQUIRE "components.a4l";  (*  ASCEND Modeling Library.
 (* => components.a4l, atoms.a4l, measures.a4l, system.a4l, basemodel.a4l *)  
 REQUIRE "phases.a4l";  
 (* => phases.a4l, system.a4l, basemodel.a4l *)  
 PROVIDE "thermodynamics.a4l";  
   
 (*  
  *  thermodynamics.a4l  
  *  by Arthur W. Westerberg, Ben Allan, and Vicente Rico-Ramirez  
  *  Part of the ASCEND Library  
  *  $Date: 1998/06/17 19:33:42 $  
  *  $Revision: 1.6 $  
  *  $Author: mthomas $  
  *  $Source: /afs/cs.cmu.edu/project/ascend/Repository/models/thermodynamics.a4l,v $  
  *  
  *  This file is part of the ASCEND Modeling Library.  
  *  
2   *  Copyright (C) 1998 Carnegie Mellon University   *  Copyright (C) 1998 Carnegie Mellon University
3   *   *
4   *  The ASCEND Modeling Library is free software; you can redistribute   *  The ASCEND Modeling Library is free software; you can redistribute
# Line 31  PROVIDE "thermodynamics.a4l"; Line 15  PROVIDE "thermodynamics.a4l";
15   *  along with the program; if not, write to the Free Software   *  along with the program; if not, write to the Free Software
16   *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139 USA.  Check   *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139 USA.  Check
17   *  the file named COPYING.   *  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;  MODEL td_model;
27  (* This MODEL is simply a root and documentation container for the  (* This MODEL is simply a root and documentation container for the
# Line 40  MODEL td_model; Line 30  MODEL td_model;
30   * the boundwidth which is standard for this library.   * the boundwidth which is standard for this library.
31   *)   *)
32  NOTES  NOTES
33  'licensing' SELF {This library is subject to the GNU Public License v.2.0}  'license' SELF {This library is subject to the GNU Public License v.2.0}
34    
35  'sacredness' SELF {  'sacredness' SELF {
36  There is nothing sacred or mysterious about either  There is nothing sacred or mysterious about either
# Line 63  single-component degenerate cases correc Line 53  single-component degenerate cases correc
53    
54  'author' SELF {Arthur W. Westerberg, Ben Allan, and Vicente Rico-Ramirez}  'author' SELF {Arthur W. Westerberg, Ben Allan, and Vicente Rico-Ramirez}
55  'created' SELF {November 9, 1997}  'created' SELF {November 9, 1997}
 'last revised' SELF {$Revision: 1.6 $}  
56    
57  'new users' SELF {  'new users' SELF {
58  The enduser MODEL in this library, if there is one, is  The end-user MODEL in this library, if there is one, is
59  "thermodynamics" which has a relatively simple interface.  "thermodynamics" which has a relatively simple interface.
60  Adding new correlations for properties is relatively easy. To  Adding new correlations for properties is relatively easy. To
61  understand how, examine by way of example the two models  understand how, examine by way of example the two models
62  "ideal_vapor_mixture" and "Pitzer_vapor_mixture".  "ideal_vapor_mixture" and "Pitzer_vapor_mixture".
63  }  }
64    
65  'telephone support' SELF {  'support' SELF {See our website at http://ascend.cheme.cmu.edu/}
 You must be kidding. CAPD consortium sponsors can contact the authors  
 for further consultancy. Everyone else can submit questions, suggestions,  
 and bug reports via the World Wide Web.  
 }  
   
 'URL' SELF {http://www.cs.cmu.edu/~ascend}  
66    
67  'design goal' SELF {  'design goal' SELF {
68  The approach used by this form for the physical property package is  The approach used by this form for the physical property package is
# Line 150  MODEL thermodynamic_properties( Line 133  MODEL thermodynamic_properties(
133      h "species molar enthalpy"  IS_A molar_energy;      h "species molar enthalpy"  IS_A molar_energy;
134      g "species molar Gibbs free energy" IS_A molar_energy;      g "species molar Gibbs free energy" IS_A molar_energy;
135    
   
136  NOTES  NOTES
137  'purpose' SELF {  'purpose' SELF {
138  The base thermodynamic property model.  Two other models in this  The base thermodynamic property model.  Two other models in this
# Line 161  enthalpy and molar Gibbs free energy of Line 143  enthalpy and molar Gibbs free energy of
143  END NOTES;  END NOTES;
144    
145  METHODS  METHODS
   
146  (* inherits ClearAll, reset *)  (* inherits ClearAll, reset *)
   
147  METHOD default_all;  METHOD default_all;
148      (* all the ATOM default values are fine so far as we know *)      (* all the ATOM default values are fine so far as we know *)
149      RUN default_self;      RUN default_self;
# Line 207  METHOD scale_self; Line 187  METHOD scale_self;
187  END scale_self;  END scale_self;
188    
189  METHOD bound_all;  METHOD bound_all;
190          (* resets most bounds to the interval centered at the current point.      (* resets most bounds to the interval centered at the current point.
191       * (variable.value) +/- boundwidth * variable.nominal       * (variable.value) +/- boundwidth * variable.nominal
192           *)      *)
193      T.lower_bound := 1.0e-8{K};      T.lower_bound := 1.0e-8{K};
194      P.lower_bound := 1.0e-8{atm};      P.lower_bound := 1.0e-8{atm};
195      T.upper_bound := T + boundwidth*T.nominal;      T.upper_bound := T + boundwidth*T.nominal;
# Line 320  of pressure and temperature for a single Line 300  of pressure and temperature for a single
300  the Pitzer vapor model.  the Pitzer vapor model.
301  }  }
302  END NOTES;  END NOTES;
303  Pitzer_component_v:      Pitzer_component_v:
304          P*v/(1{GAS_C})/data.Tc = T/data.Tc + (P/data.Pc)*              P*v/(1{GAS_C})/data.Tc = T/data.Tc + (P/data.Pc)*
     (0.083 - 0.422*(data.Tc/T)^1.6 + data.omega*  
         (0.139 - 0.172*(data.Tc/T)^4.2));  
   
 Pitzer_component_h:  
     h/(1{GAS_C})/data.Tc = data.H0/(1{GAS_C})/data.Tc +  
     data.cpvapa*(T - data.T0)/(1{GAS_C})/data.Tc +  
     data.cpvapb*(T^2 - data.T0^2)/2/(1{GAS_C})/data.Tc +  
     data.cpvapc*(T^3 - data.T0^3)/3/(1{GAS_C})/data.Tc +  
     data.cpvapd*(T^4 - data.T0^4)/4/(1{GAS_C})/data.Tc +  
     (P/data.Pc)*(0.083 - 1.097*(data.Tc/T)^1.6 + data.omega*  
         (0.139 - 0.894*(data.Tc/T)^4.2));  
   
 Pitzer_component_g:  
         g/(1{GAS_C})/data.Tc = data.G0/(1{GAS_C})/data.Tc -  
         (data.H0 - data.G0)*(T/data.T0 - 1)/(1{GAS_C})/data.Tc -  
         data.cpvapa*(T*lnm(T/data.T0) - T + data.T0)/1{GAS_C}/data.Tc -  
         data.cpvapb*(T^2 - 2*T*data.T0 + data.T0^2)/  
         (2*1{GAS_C}*data.Tc) -  
         data.cpvapc*(T^3/2 - 3*T*data.T0^2/2 + data.T0^3)/  
         (3*1{GAS_C}*data.Tc) -  
         data.cpvapd*(T^4/3 - 4*T*data.T0^3/3 + data.T0^4)/  
         (4*1{GAS_C}*data.Tc) +  
         T*lnm(P/data.P0)/data.Tc +  
         (P/data.Pc)*  
305          (0.083 - 0.422*(data.Tc/T)^1.6 + data.omega*          (0.083 - 0.422*(data.Tc/T)^1.6 + data.omega*
306          (0.139 - 0.172*(data.Tc/T)^4.2));              (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  METHODS
333  (* Inherited methods: ClearAll, check_all, check_self,  (* Inherited methods: ClearAll, check_all, check_self,
# Line 380  END Pitzer_vapor_component; Line 360  END Pitzer_vapor_component;
360    
361  (* ****************************************************************** *)  (* ****************************************************************** *)
362  MODEL ideal_vapor_component(  MODEL ideal_vapor_component(
363      P WILL_BE pressure;      P WILL_BE pressure;
364      T WILL_BE temperature;      T WILL_BE temperature;
365      data WILL_BE td_component_constants;      data WILL_BE td_component_constants;
366  ) REFINES pure_component;  ) REFINES pure_component;
367    
368  NOTES  NOTES
# Line 395  the ideal gas (vapor) model. Line 375  the ideal gas (vapor) model.
375  'FIXME' SELF {Needs to be independently checked. Did I get h and g correct?}  'FIXME' SELF {Needs to be independently checked. Did I get h and g correct?}
376  END NOTES;  END NOTES;
377    
378  ideal_component_v: P*v/(1{GAS_C})/data.Tc = T/data.Tc;      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_h:      ideal_component_g:
389      h/(1{GAS_C})/data.Tc =          g/(1{GAS_C})/data.Tc =
390      data.H0/(1{GAS_C})/data.Tc +          data.G0/(1{GAS_C})/data.Tc -
391      data.cpvapa*(T - data.T0)/(1{GAS_C})/data.Tc +          (data.H0 - data.G0)*(T/data.T0 - 1)/(1{GAS_C})/data.Tc -
392      data.cpvapb*(T^2 - data.T0^2)/2/(1{GAS_C})/data.Tc +          data.cpvapa*(T*lnm(T/data.T0) - T + data.T0)/1{GAS_C}/data.Tc -
393      data.cpvapc*(T^3 - data.T0^3)/3/(1{GAS_C})/data.Tc +          data.cpvapb*(T^2 - 2*T*data.T0 + data.T0^2)/(2*1{GAS_C}*data.Tc) -
394      data.cpvapd*(T^4 - data.T0^4)/4/(1{GAS_C})/data.Tc;          data.cpvapc*(T^3/2 - 3*T*data.T0^2/2 + data.T0^3)/
395              (3*1{GAS_C}*data.Tc) -
396  ideal_component_g:          data.cpvapd*(T^4/3 - 4*T*data.T0^3/3 + data.T0^4)/
397      g/(1{GAS_C})/data.Tc =            (4*1{GAS_C}*data.Tc) +
398      data.G0/(1{GAS_C})/data.Tc -          T*lnm(P/data.P0)/data.Tc;
     (data.H0 - data.G0)*(T/data.T0 - 1)/(1{GAS_C})/data.Tc -  
     data.cpvapa*(T*lnm(T/data.T0) - T + data.T0)/1{GAS_C}/data.Tc -  
     data.cpvapb*(T^2 - 2*T*data.T0 + data.T0^2)/(2*1{GAS_C}*data.Tc) -  
     data.cpvapc*(T^3/2 - 3*T*data.T0^2/2 + data.T0^3)/  
       (3*1{GAS_C}*data.Tc) -  
     data.cpvapd*(T^4/3 - 4*T*data.T0^3/3 + data.T0^4)/  
       (4*1{GAS_C}*data.Tc) +  
     T*lnm(P/data.P0)/data.Tc;  
399    
400  METHODS  METHODS
401  (* Inherited methods: ClearAll, check_all, check_self,  (* Inherited methods: ClearAll, check_all, check_self,
# Line 449  END ideal_vapor_component; Line 429  END ideal_vapor_component;
429  MODEL Rackett_liquid_component(  MODEL Rackett_liquid_component(
430      P WILL_BE pressure;      P WILL_BE pressure;
431      T WILL_BE temperature;      T WILL_BE temperature;
432          data WILL_BE td_component_constants;      data WILL_BE td_component_constants;
433  ) REFINES pure_component;  ) REFINES pure_component;
434    
435  NOTES  NOTES
# Line 589  the phase -- which can be used later to Line 569  the phase -- which can be used later to
569  for multi-phase mixtures.  for multi-phase mixtures.
570  }  }
571  END NOTES;  END NOTES;
572      P IS_A pressure;      P IS_A pressure;
573      T IS_A temperature;      T IS_A temperature;
574    
575      v_y "mixture molar volume" IS_A molar_volume;      v_y "mixture molar volume" IS_A molar_volume;
576      h_y "mixture molar enthalpy" IS_A molar_energy;      h_y "mixture molar enthalpy" IS_A molar_energy;
577      g_y "mixture molar Gibbs free energy" IS_A molar_energy;      g_y "mixture molar Gibbs free energy" IS_A molar_energy;
578    
579      components ALIASES cd.components;      components ALIASES cd.components;
580      partial[components] IS_A partial_component(P, T);      partial[components] IS_A partial_component(P, T);
581      alpha[components] IS_A relative_volatility;      alpha[components] IS_A relative_volatility;
582    
583  NOTES 'scaling' v_mix,h_mix,g_mix {  NOTES 'scaling' v_mix,h_mix,g_mix {
584  DIMENSIONLESS well-scaled equations are easier to converge from anywhere,  DIMENSIONLESS well-scaled equations are easier to converge from anywhere,
# Line 607  v*Pcref/(R*Tcref) = Sum(yi*vi |i) * Pcre Line 587  v*Pcref/(R*Tcref) = Sum(yi*vi |i) * Pcre
587  comes from v = Sum(yi*vi);  comes from v = Sum(yi*vi);
588  } END NOTES;  } END NOTES;
589    
590      v_mix: v_y * (cd.data[cd.reference].Pc /      v_mix: v_y * (cd.data[cd.reference].Pc /
591                 (cd.data[cd.reference].Tc * 1{GAS_C}))                 (cd.data[cd.reference].Tc * 1{GAS_C}))
592               = SUM[y[i]*partial[i].v | i IN cd.components] *               = SUM[y[i]*partial[i].v | i IN cd.components] *
593                 (cd.data[cd.reference].Pc /                 (cd.data[cd.reference].Pc /
594                 (cd.data[cd.reference].Tc * 1{GAS_C}));                 (cd.data[cd.reference].Tc * 1{GAS_C}));
595    
596      (* likewise for h and g *)      (* likewise for h and g *)
597      h_mix: h_y / (cd.data[cd.reference].Tc * 1{GAS_C})      h_mix: h_y / (cd.data[cd.reference].Tc * 1{GAS_C})
598               = SUM[y[i]*partial[i].h | i IN cd.components] /               = SUM[y[i]*partial[i].h | i IN cd.components] /
599                 (cd.data[cd.reference].Tc * 1{GAS_C});                 (cd.data[cd.reference].Tc * 1{GAS_C});
600    
601      g_mix: g_y / (cd.data[cd.reference].Tc * 1{GAS_C})      g_mix: g_y / (cd.data[cd.reference].Tc * 1{GAS_C})
602               = SUM[y[i]*partial[i].g | i IN cd.components] /               = SUM[y[i]*partial[i].g | i IN cd.components] /
603                 (cd.data[cd.reference].Tc * 1{GAS_C});                 (cd.data[cd.reference].Tc * 1{GAS_C});
604    
605      y[components] "mixture species mole fractions" IS_A mole_fraction;      y[components] "mixture species mole fractions" IS_A mole_fraction;
606    
607  NOTES 'slack explanation' slack_PhaseDisappearance {  NOTES 'slack explanation' slack_PhaseDisappearance {
608  When slack_PhaseDisappearance != 0.0, then the phase  When slack_PhaseDisappearance != 0.0, then the phase
609  does not exist in this formulation of multi-phase  does not exist in this formulation of multi-phase
610  equilibrium thermodynamics.  equilibrium thermodynamics.
611  } END NOTES;  } END NOTES;
612      slack_PhaseDisappearance "complementarity slack variable" IS_A fraction;      slack_PhaseDisappearance "complementarity slack variable" IS_A fraction;
613    
614      sum_y: SUM[y[components]] + slack_PhaseDisappearance  = 1.0;      sum_y: SUM[y[components]] + slack_PhaseDisappearance  = 1.0;
615    
616  METHODS  METHODS
617    
618  (* Inherited methods: ClearAll, reset.  *)  (* Inherited methods: ClearAll, reset.  *)
619  METHOD default_self;  METHOD default_self;
620          T.lower_bound := 1.0{K};      T.lower_bound := 1.0{K};
621          P.lower_bound := 1.0{Pa};      P.lower_bound := 1.0{Pa};
622          v_y := 24{liter/mole}; (* presupposes a vapor *)      v_y := 24{liter/mole}; (* presupposes a vapor *)
623          RUN partial[components].default_self;      RUN partial[components].default_self;
624  END default_self;  END default_self;
625    
626  METHOD default_all;  METHOD default_all;
# Line 664  METHOD check_self; Line 644  METHOD check_self;
644          STOP {Phase has negative molar volume?!};          STOP {Phase has negative molar volume?!};
645      END IF;      END IF;
646      (* Don't know how to check g_y,h_y *)      (* Don't know how to check g_y,h_y *)
647          RUN partial[components].check_self;      RUN partial[components].check_self;
648  END check_self;  END check_self;
649    
650  METHOD scale_all;  METHOD scale_all;
# Line 679  METHOD scale_self; Line 659  METHOD scale_self;
659      v_y.nominal := v_y;      v_y.nominal := v_y;
660      h_y.nominal := abs(h_y);      h_y.nominal := abs(h_y);
661      g_y.nominal := abs(g_y);      g_y.nominal := abs(g_y);
662          RUN partial[components].scale_self;      RUN partial[components].scale_self;
663  END scale_self;  END scale_self;
664    
665  METHOD bound_all;  METHOD bound_all;
# Line 726  METHOD specify_partials_seqmod; Line 706  METHOD specify_partials_seqmod;
706  END specify_partials_seqmod;  END specify_partials_seqmod;
707    
708  METHOD specify_partials_computed;  METHOD specify_partials_computed;
709         FIX P;      FIX P;
710         FIX T;      FIX T;
711         FIX alpha[components];      FIX alpha[components];
712         RUN partial[components].specify_partials_computed;      RUN partial[components].specify_partials_computed;
713         FIX y[cd.other_components];      FIX y[cd.other_components];
714         FIX slack_PhaseDisappearance;      FIX slack_PhaseDisappearance;
715         slack_PhaseDisappearance := 0.0;      slack_PhaseDisappearance := 0.0;
716  END specify_partials_computed;  END specify_partials_computed;
717    
718  METHOD specify;  METHOD specify;
719         FIX P;      FIX P;
720         FIX T;      FIX T;
721         FIX alpha[components];      FIX alpha[components];
722         RUN partial[components].specify;      RUN partial[components].specify;
723         FIX y[cd.other_components];      FIX y[cd.other_components];
724         FIX slack_PhaseDisappearance;      FIX slack_PhaseDisappearance;
725         slack_PhaseDisappearance := 0.0;      slack_PhaseDisappearance := 0.0;
726  END specify;  END specify;
727    
728  END phase_partials;  END phase_partials;
# Line 750  END phase_partials; Line 730  END phase_partials;
730  (* ****************************************************************** *)  (* ****************************************************************** *)
731    
732  MODEL ideal_vapor_mixture(  MODEL ideal_vapor_mixture(
733      cd WILL_BE components_data;      cd WILL_BE components_data;
734  ) REFINES phase_partials();  ) REFINES phase_partials();
735    
736  NOTES 'purpose' SELF {  NOTES 'purpose' SELF {
# Line 760  phase for molar volume, molar enthalpy a Line 740  phase for molar volume, molar enthalpy a
740  based on the ideal gas model for a mixture.  based on the ideal gas model for a mixture.
741  } END NOTES;  } END NOTES;
742    
743      data ALIASES cd.data;      data ALIASES cd.data;
744    
745      FOR i IN components CREATE      FOR i IN components CREATE
746          pure[i] "pure species ideal gas properties for the i-th species"          pure[i] "pure species ideal gas properties for the i-th species"
747          IS_A ideal_vapor_component(P, T, data[i]);              IS_A ideal_vapor_component(P, T, data[i]);
748      END FOR;      END FOR;
749    
750      FOR i IN components CREATE      FOR i IN components CREATE
   
751      ideal_mixture_v[i]:      ideal_mixture_v[i]:
752          (partial[i].v - pure[i].v) *          (partial[i].v - pure[i].v) *
753          (data[i].Pc/(1{GAS_C}*data[i].Tc)) = 0;          (data[i].Pc/(1{GAS_C}*data[i].Tc)) = 0;
# Line 776  based on the ideal gas model for a mixtu Line 755  based on the ideal gas model for a mixtu
755      ideal_mixture_h[i]:      ideal_mixture_h[i]:
756          (partial[i].h - pure[i].h)/((1{GAS_C})*data[i].Tc) = 0;          (partial[i].h - pure[i].h)/((1{GAS_C})*data[i].Tc) = 0;
757    
758          ideal_mixture_g[i]:      ideal_mixture_g[i]:
759          (partial[i].g - pure[i].g - (1{GAS_C})*T*lnm(y[i])) /          (partial[i].g - pure[i].g - (1{GAS_C})*T*lnm(y[i])) /
760          (1{GAS_C}*data[i].Tc) = 0;          (1{GAS_C}*data[i].Tc) = 0;
761      END FOR;      END FOR;
# Line 815  END ideal_vapor_mixture; Line 794  END ideal_vapor_mixture;
794  (* ****************************************************************** *)  (* ****************************************************************** *)
795    
796  MODEL Pitzer_vapor_mixture(  MODEL Pitzer_vapor_mixture(
797      cd WILL_BE components_data;      cd WILL_BE components_data;
798  ) REFINES phase_partials();  ) REFINES phase_partials();
799    
800  NOTES 'purpose' SELF {  NOTES 'purpose' SELF {
# Line 826  based on the Pitzer vapor model for a mi Line 805  based on the Pitzer vapor model for a mi
805  }  }
806  END NOTES;  END NOTES;
807    
808      data ALIASES cd.data;      data ALIASES cd.data;
809    
810      FOR i IN components CREATE      FOR i IN components CREATE
811          pure[i]          pure[i] "pure species properties for the i-th species using acentric factors"
812       "pure species properties for the i-th species using acentric factors"              IS_A Pitzer_vapor_component(P, T,  data[i]);
813          IS_A Pitzer_vapor_component(P, T,  data[i]);      END FOR;
     END FOR;  
814    
815      FOR i IN components CREATE      FOR i IN components CREATE
816    
817      Pitzer_mixture_v[i]:      Pitzer_mixture_v[i]:
818           (partial[i].v - pure[i].v)*(data[i].Pc/           (partial[i].v - pure[i].v)*(data[i].Pc/
# Line 862  END NOTES; Line 840  END NOTES;
840          (0.139 - 0.894*(sqrt(data[i].Tc*data[j].Tc)/T)^4.2))          (0.139 - 0.894*(sqrt(data[i].Tc*data[j].Tc)/T)^4.2))
841            | j IN components - [i]];            | j IN components - [i]];
842    
843          Pitzer_mixture_g[i]:      Pitzer_mixture_g[i]:
844          (partial[i].g - pure[i].g - (1{GAS_C})*T*lnm(y[i])) /          (partial[i].g - pure[i].g - (1{GAS_C})*T*lnm(y[i])) /
845          (1{GAS_C}*data[i].Tc) =          (1{GAS_C}*data[i].Tc) =
846          -(P/data[i].Pc)*          -(P/data[i].Pc)*
# Line 912  END Pitzer_vapor_mixture; Line 890  END Pitzer_vapor_mixture;
890    
891    
892  MODEL UNIFAC_liquid_mixture(  MODEL UNIFAC_liquid_mixture(
893      cd WILL_BE components_data;      cd WILL_BE components_data;
894  ) REFINES phase_partials();  ) REFINES phase_partials();
895    
896  NOTES 'purpose' SELF {  NOTES 'purpose' SELF {
# Line 922  phase for molar volume, molar enthalpy a Line 900  phase for molar volume, molar enthalpy a
900  based on the UNIFAC model for a liquid mixture.  based on the UNIFAC model for a liquid mixture.
901  } END NOTES;  } END NOTES;
902    
903      data ALIASES cd.data;      data ALIASES cd.data;
904    
905      FOR i IN components CREATE      FOR i IN components CREATE
906          pure[i]          pure[i] "pure species liquid properties for the i-th species by Rackett"
907        "pure species liquid properties for the i-th species by Rackett"              IS_A Rackett_liquid_component(P, T, data[i]);
         IS_A Rackett_liquid_component(P, T, data[i]);  
908      END FOR;      END FOR;
909    
910      subgroups       IS_A set OF symbol_constant;      subgroups       IS_A set OF symbol_constant;
# Line 967  based on the UNIFAC model for a liquid m Line 944  based on the UNIFAC model for a liquid m
944                  | k IN data[i].subgroups];                  | k IN data[i].subgroups];
945      END FOR;      END FOR;
946      FOR i IN components CREATE      FOR i IN components CREATE
947          UNIFAC_mixture_rv[i]:      UNIFAC_mixture_rv[i]:
948          rv[i] = Jv[i]*SUM[rv[j]*y[j] | j IN components];          rv[i] = Jv[i]*SUM[rv[j]*y[j] | j IN components];
949    
950          UNIFAC_mixture_qs[i]:      UNIFAC_mixture_qs[i]:
951          qs[i] = Ls[i]*SUM[qs[j]*y[j] | j IN components];          qs[i] = Ls[i]*SUM[qs[j]*y[j] | j IN components];
952          UNIFAC_partial_v[i]:      UNIFAC_partial_v[i]:
953          partial[i].v = pure[i].v;          partial[i].v = pure[i].v;
954    
955          UNIFAC_mixture_h_excess[i]:      UNIFAC_mixture_h_excess[i]:
956          (partial[i].h - pure[i].h)/(1{GAS_C})/data[i].Tc =          (partial[i].h - pure[i].h)/(1{GAS_C})/data[i].Tc =
957          SUM[0,theta[k]*          SUM[0,theta[k]*
958          SUM[0,SUM[0,theta[n]*              SUM[0,SUM[0,theta[n]*
959          ((uc.a[g][uc.group[k]] -                      ((uc.a[g][uc.group[k]] -
960          uc.a[uc.group[n]][uc.group[k]])/data[i].Tc)*                      uc.a[uc.group[n]][uc.group[k]])/data[i].Tc)*
961          exp(-(uc.a[g][uc.group[k]] +                      exp(-(uc.a[g][uc.group[k]] +
962          uc.a[uc.group[n]][uc.group[k]])/T)*                      uc.a[uc.group[n]][uc.group[k]])/T)*
963          SUM[0,data[i].nu[m]*uc.Q[m]                      SUM[0,data[i].nu[m]*uc.Q[m]
964          | m IN data[i].subgroups*uc.sub[g]]                          | m IN data[i].subgroups*uc.sub[g]
965          | g IN data[i].groups - [uc.group[n]]]                      ]
966          | n IN subgroups]/eta[k]/eta[k]                  | g IN data[i].groups - [uc.group[n]]]
967                | n IN subgroups]/eta[k]/eta[k]
968          | k IN subgroups] -          | k IN subgroups] -
969          SUM[0,(data[i].nu[k]*uc.Q[k]/(          SUM[0,(data[i].nu[k]*uc.Q[k]/(
970          SUM[0,data[i].nu[m]*uc.Q[m]              SUM[0,data[i].nu[m]*uc.Q[m]
971          | m IN data[i].subgroups*uc.sub[uc.group[k]]] +              | m IN data[i].subgroups*uc.sub[uc.group[k]]] +
972          SUM[0,SUM[0,data[i].nu[m]*uc.Q[m] *              SUM[0,SUM[0,data[i].nu[m]*uc.Q[m] *
973                  exp(-uc.a[g][uc.group[k]]/T)                      exp(-uc.a[g][uc.group[k]]/T)
974          | m IN data[i].subgroups*uc.sub[g]]                  | m IN data[i].subgroups*uc.sub[g]]
975          | g IN data[i].groups - [uc.group[k]]]))*              | g IN data[i].groups - [uc.group[k]]]))*
976          SUM[0,SUM[0,theta[n]*              SUM[0,SUM[0,theta[n]*
977          ((uc.a[g][uc.group[k]] -                      ((uc.a[g][uc.group[k]] -
978          uc.a[uc.group[n]][uc.group[k]])/data[i].Tc)*                      uc.a[uc.group[n]][uc.group[k]])/data[i].Tc)*
979          exp(-(uc.a[g][uc.group[k]] +                      exp(-(uc.a[g][uc.group[k]] +
980          uc.a[uc.group[n]][uc.group[k]])/T)*                      uc.a[uc.group[n]][uc.group[k]])/T)*
981          SUM[0,data[i].nu[m]*uc.Q[m]                      SUM[0,data[i].nu[m]*uc.Q[m]
982          | m IN data[i].subgroups*uc.sub[g]]                      | m IN data[i].subgroups*uc.sub[g]]
983          | g IN data[i].groups - [uc.group[n]]]                  | g IN data[i].groups - [uc.group[n]]]
984          | n IN subgroups]/eta[k]              | n IN subgroups]/eta[k]
985          | k IN data[i].subgroups];          | k IN data[i].subgroups];
986    
987          UNIFAC_mixture_g_excess[i]:      UNIFAC_mixture_g_excess[i]:
988          (partial[i].g - pure[i].g - (1{GAS_C})*T*lnm(y[i])) /          (partial[i].g - pure[i].g - (1{GAS_C})*T*lnm(y[i])) /
989          (1{GAS_C}*data[i].Tc) =          (1{GAS_C}*data[i].Tc) =
990          (1.0 - Jv[i] + lnm(Jv[i]) -          (1.0 - Jv[i] + lnm(Jv[i]) -
991          5.0*qs[i]*(1.0 - Jv[i]/Ls[i] + lnm(Jv[i]/Ls[i])) +          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 -          qs[i]*(1 - lnm(Ls[i])))*T/data[i].Tc -
993          SUM[0,theta[k]*(          SUM[0,theta[k]*(
994          SUM[0,data[i].nu[m]*uc.Q[m]              SUM[0,data[i].nu[m]*uc.Q[m]
995          | m IN data[i].subgroups*uc.sub[uc.group[k]]] +              | m IN data[i].subgroups*uc.sub[uc.group[k]]] +
996          SUM[0,SUM[0,data[i].nu[m] * uc.Q[m] *              SUM[0,SUM[0,data[i].nu[m] * uc.Q[m] *
997                  exp(-uc.a[g][uc.group[k]]/T)                      exp(-uc.a[g][uc.group[k]]/T)
998          | m IN data[i].subgroups*uc.sub[g]]                  | m IN data[i].subgroups*uc.sub[g]]
999          | g IN data[i].groups - [uc.group[k]]])/eta[k]              | g IN data[i].groups - [uc.group[k]]])/eta[k]
1000          | k IN subgroups]*T/data[i].Tc +          | k IN subgroups]*T/data[i].Tc +
1001          SUM[0,data[i].nu[k]*uc.Q[k]*lnm((          SUM[0,data[i].nu[k]*uc.Q[k]*lnm((
1002          SUM[0,data[i].nu[m]*uc.Q[m]              SUM[0,data[i].nu[m]*uc.Q[m]
1003          | m IN data[i].subgroups*uc.sub[uc.group[k]]] +              | m IN data[i].subgroups*uc.sub[uc.group[k]]] +
1004          SUM[0,SUM[0, data[i].nu[m] * uc.Q[m] *              SUM[0,SUM[0, data[i].nu[m] * uc.Q[m] *
1005                  exp(-uc.a[g][uc.group[k]]/T)                      exp(-uc.a[g][uc.group[k]]/T)
1006          | m IN data[i].subgroups*uc.sub[g]]                  | m IN data[i].subgroups*uc.sub[g]]
1007          | g IN data[i].groups - [uc.group[k]]])/eta[k])              | g IN data[i].groups - [uc.group[k]]])/eta[k])
1008          | k IN data[i].subgroups]*T/data[i].Tc;          | k IN data[i].subgroups]*T/data[i].Tc;
1009      END FOR;      END FOR;
1010    
# Line 1124  phase for molar volume, molar enthalpy a Line 1102  phase for molar volume, molar enthalpy a
1102  based on the Wilson model for a liquid mixture.  based on the Wilson model for a liquid mixture.
1103  } END NOTES;  } END NOTES;
1104    
1105      data ALIASES cd.data;      data ALIASES cd.data;
1106    
1107      FOR i IN components CREATE      FOR i IN components CREATE
1108          pure[i] IS_A Rackett_liquid_component(P, T,  data[i]);          pure[i] IS_A Rackett_liquid_component(P, T,  data[i]);
1109      END FOR;      END FOR;
1110    
1111      lambda[components][components] "Wilson binary parameters"      lambda[components][components] "Wilson binary parameters"
1112          IS_A factor;          IS_A factor;
# Line 1160  based on the Wilson model for a liquid m Line 1138  based on the Wilson model for a liquid m
1138          | s IN components]) -          | s IN components]) -
1139          (SUM[y[k]*lambda[k][i] | k IN components]*          (SUM[y[k]*lambda[k][i] | k IN components]*
1140          SUM[SUM[y[l]*lambda[m][l]*          SUM[SUM[y[l]*lambda[m][l]*
1141           data[m].del_ip[data[l].formula]                  data[m].del_ip[data[l].formula]
1142          | l IN components]              | l IN components]
1143          | m IN components]) / (          | m IN components]) / (
1144          SUM[SUM[y[n]*lambda[o][n]          SUM[SUM[y[n]*lambda[o][n]
1145          | n IN components]              | n IN components]
1146          | o IN components])^2 +          | o IN components])^2 +
1147          (SUM[y[p]*lambda[p][i]*          (SUM[y[p]*lambda[p][i]*
1148           data[p].del_ip[data[i].formula]               data[p].del_ip[data[i].formula]
1149           | p IN components]) /           | p IN components]) /
1150           (SUM[SUM[y[q]*lambda[r][q]           (SUM[SUM[y[q]*lambda[r][q]
1151           | q IN components]              | q IN components]
1152           | r IN components]);           | r IN components]);
1153      END FOR;      END FOR;
1154    
# Line 1190  and makes the molar volumes initial valu Line 1168  and makes the molar volumes initial valu
1168  (* The following initial bounds on lambda will only be useful  (* The following initial bounds on lambda will only be useful
1169    if meaninful bounds on pure[i].v have been set. *)    if meaninful bounds on pure[i].v have been set. *)
1170      FOR i IN components DO      FOR i IN components DO
1171          FOR j IN components DO          FOR j IN components DO
1172          IF (pure[i].v.upper_bound > 0{liter/mole}) THEN              IF (pure[i].v.upper_bound > 0{liter/mole}) THEN
1173              lambda[i][j].lower_bound :=                  lambda[i][j].lower_bound :=
1174              (pure[j].v.lower_bound / pure[i].v.upper_bound) *                      (pure[j].v.lower_bound / pure[i].v.upper_bound) *
1175              exp(-data[i].del_ip[data[j].formula] /                      exp(-data[i].del_ip[data[j].formula] /
1176              (1{GAS_C} * T));                      (1{GAS_C} * T));
1177          END IF;              END IF;
1178          IF (pure[i].v.lower_bound > 0{liter/mole}) THEN              IF (pure[i].v.lower_bound > 0{liter/mole}) THEN
1179              lambda[i][j].upper_bound :=                  lambda[i][j].upper_bound :=
1180              (pure[j].v.upper_bound / pure[i].v.lower_bound) *                      (pure[j].v.upper_bound / pure[i].v.lower_bound) *
1181              exp(-data[i].del_ip[data[j].formula] /                      exp(-data[i].del_ip[data[j].formula] /
1182              (1{GAS_C} * T));                      (1{GAS_C} * T));
1183          END IF;              END IF;
1184          END FOR;          END FOR;
1185      END FOR;      END FOR;
1186    
1187      (* liquid molar volumes are much smaller than vapor volumes *)      (* liquid molar volumes are much smaller than vapor volumes *)
# Line 1242  METHOD bound_self; Line 1220  METHOD bound_self;
1220      RUN phase_partials::bound_self;      RUN phase_partials::bound_self;
1221      RUN pure[components].bound_self;      RUN pure[components].bound_self;
1222      FOR i IN components DO      FOR i IN components DO
1223          FOR j IN components DO          FOR j IN components DO
1224          IF (pure[i].v.upper_bound > 0{liter/mole}) THEN              IF (pure[i].v.upper_bound > 0{liter/mole}) THEN
1225              lambda[i][j].lower_bound :=                  lambda[i][j].lower_bound :=
1226              (pure[j].v.lower_bound / pure[i].v.upper_bound) *                      (pure[j].v.lower_bound / pure[i].v.upper_bound) *
1227              exp(-data[i].del_ip[data[j].formula] /                      exp(-data[i].del_ip[data[j].formula] /
1228              (1{GAS_C} * T));                      (1{GAS_C} * T));
1229          END IF;              END IF;
1230          IF (pure[i].v.lower_bound > 0{liter/mole}) THEN              IF (pure[i].v.lower_bound > 0{liter/mole}) THEN
1231              lambda[i][j].upper_bound :=                  lambda[i][j].upper_bound :=
1232              (pure[j].v.upper_bound / pure[i].v.lower_bound) *                      (pure[j].v.upper_bound / pure[i].v.lower_bound) *
1233              exp(-data[i].del_ip[data[j].formula] /                      exp(-data[i].del_ip[data[j].formula] /
1234              (1{GAS_C} * T));                      (1{GAS_C} * T));
1235          END IF;              END IF;
1236          END FOR;          END FOR;
1237      END FOR;      END FOR;
1238  END bound_self;  END bound_self;
1239    
# Line 1269  END Wilson_liquid_mixture; Line 1247  END Wilson_liquid_mixture;
1247    
1248    
1249  MODEL thermodynamics(  MODEL thermodynamics(
1250     cd WILL_BE components_data;      cd WILL_BE components_data;
1251     pd WILL_BE phases_data;      pd WILL_BE phases_data;
1252     phase[pd.phases] WILL_BE phase_partials;      phase[pd.phases] WILL_BE phase_partials;
1253     equilibrated WILL_BE boolean;      equilibrated WILL_BE boolean;
1254  ) REFINES td_model;  ) REFINES td_model;
1255    
1256  NOTES  NOTES
# Line 1298  the phases (equilibrated := TRUE). Line 1276  the phases (equilibrated := TRUE).
1276  }  }
1277  END NOTES;  END NOTES;
1278    
1279     P ALIASES phase[pd.reference_phase].P;      P ALIASES phase[pd.reference_phase].P;
1280     T ALIASES phase[pd.reference_phase].T;      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     V "molar average volume over all species i, phases j" IS_A molar_volume;      FOR i IN cd.components CREATE
1321     H "molar average enthalpy over all species i, phases j" IS_A molar_energy;      overall_y[i]: y[i] = SUM[phase_fraction[j]*phase[j].y[i] | j IN pd.phases];
1322     G "molar average Gibbs free energy over all species i, phases j"      END FOR;
     IS_A molar_energy;  
    y[cd.components] "average mole fractions over all phases j" IS_A fraction;  
    phase_fraction[pd.phases] "division of moles into phases " IS_A fraction;  
    alpha_bar[pd.other_phases] "mole fraction weighted sum of volatilities"  
      IS_A relative_volatility;  
    number_of_other_phases "Np - 1" IS_A integer_constant;  
   
    FOR j IN [pd.phases] CREATE  
      slack_PhaseDisappearance[j] ALIASES phase[j].slack_PhaseDisappearance;  
      (* If you set slack_PhaseDisappearance[j].fixed := TRUE,  
       * (useful if you want to _force_ a phase to exist)  
       * then you must also set complementarity[j].included := FALSE;  
       * OTHERWISE the Jacobian matrix of this MODEL is singular.  
       * The only sensible value of slack_PhaseDisappearance[i]  
       * when it is fixed is 0.0.  
       *)  
    END FOR;  
   
    sum_phase_fractions: SUM[phase_fraction[pd.phases]] = 1.0;  
    overall_v:   V * (cd.data[cd.reference].Pc /  
                      (1{GAS_C} * cd.data[cd.reference].Tc))  
                   = SUM[phase_fraction[j]*phase[j].v_y | j IN pd.phases] *  
                      (cd.data[cd.reference].Pc /  
                      (1{GAS_C} * cd.data[cd.reference].Tc));  
    overall_h:   H / (1{GAS_C} * cd.data[cd.reference].Tc)  
                   = SUM[phase_fraction[j]*phase[j].h_y | j IN pd.phases] /  
                     (1{GAS_C} * cd.data[cd.reference].Tc);  
    overall_g:   G / (1{GAS_C} * cd.data[cd.reference].Tc)  
                   = SUM[phase_fraction[j]*phase[j].g_y | j IN pd.phases]/  
                     (1{GAS_C} * cd.data[cd.reference].Tc);  
   
    FOR i IN cd.components CREATE  
    overall_y[i]: y[i] = SUM[phase_fraction[j]*phase[j].y[i] | j IN pd.phases];  
    END FOR;  
1323    
1324     number_of_other_phases :== CARD[pd.other_phases];      number_of_other_phases :== CARD[pd.other_phases];
1325    
1326     SELECT (number_of_other_phases)      SELECT (number_of_other_phases)
1327      CASE 0:      CASE 0:
1328         (* Create Nothing *)          (* Create Nothing *)
1329      OTHERWISE:      OTHERWISE:
1330    
1331      (* Create equilibrium relationships *)          (* Create equilibrium relationships *)
1332    
1333      (* equilibrium *)          (* equilibrium *)
   
     (* pressure and temperature equilibrium.  
          * this condition of equilibrium is not UNIVERSAL  
          * to all equilibrium thermodynamic systems.  
          * It certainly applies to non-electrolytic Vapor-Liquid  
          * systems, however.  
          *)  
     FOR j IN [pd.other_phases] CREATE  
             equateP[j]: phase[j].P = phase[pd.reference_phase].P;  
             equateT[j]: phase[j].T = phase[pd.reference_phase].T;  
     END FOR;  
1334    
1335  NOTES          (* pressure and temperature equilibrium.
1336  'usage' alpha_bar, alphaeqn, phase[pd.phases].alpha[components] {          * this condition of equilibrium is not UNIVERSAL
1337  If using a relative volatility version for equilibrium, fix          * to all equilibrium thermodynamic systems.
1338  the relative volatilities alpha for all components in all          * It certainly applies to non-electrolytic Vapor-Liquid
1339  phases (except perhaps the reference phase) and compute          * systems, however.
1340  alpha_bar.          *)
1341            FOR j IN [pd.other_phases] CREATE
1342  Else activate the equations that equate the chemical              equateP[j]: phase[j].P = phase[pd.reference_phase].P;
1343  potential for each component in each phase to its value in              equateT[j]: phase[j].T = phase[pd.reference_phase].T;
1344  the reference phase, fix alpha_bar and set its value to one          END FOR;
 and compute alphas for all the components.  
1345    
1346  Note, there is a vector of alphas in the reference phase          NOTES
1347  that is unused.          'usage' alpha_bar, alphaeqn, phase[pd.phases].alpha[components] {
1348  }          If using a relative volatility version for equilibrium, fix
1349  'interpretation' alpha {          the relative volatilities alpha for all components in all
1350  When alpha_bar is fixed at 1.0, then the alpha[i] are essentially          phases (except perhaps the reference phase) and compute
1351  K-values. When alpha bar is computed from fixed alpha[i], the          alpha_bar.
1352  alpha[i] are essentially relative volatilities. The oldest of the  
1353  authors believes this note to be entirely redundant, but the          Else activate the equations that equate the chemical
1354  youngest has found it helps some people to relate the MODEL to          potential for each component in each phase to its value in
1355  their old unit operations textbooks.          the reference phase, fix alpha_bar and set its value to one
1356  }          and compute alphas for all the components.
1357  END NOTES;  
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          FOR j IN [pd.other_phases] CREATE
1380             FOR i IN cd.components CREATE              FOR i IN cd.components CREATE
1381                alphaeqn[j][i]:              td_eql[j][i]:
1382                  alpha_bar[j] * phase[j].y[i]                  phase[j].partial[i].g = phase[pd.reference_phase].partial[i].g;
1383                    = phase[j].alpha[i] * phase[pd.reference_phase].y[i];              END FOR;
1384             END FOR;          END FOR;
        END FOR;  
   
        FOR j IN [pd.other_phases] CREATE  
            FOR i IN cd.components CREATE  
            td_eql[j][i]:  
                   phase[j].partial[i].g  
             = phase[pd.reference_phase].partial[i].g;  
            END FOR;  
        END FOR;  
1385    
1386         WHEN (equilibrated)          WHEN (equilibrated)
1387             CASE TRUE:          CASE TRUE:
1388              USE td_eql;              USE td_eql;
1389             OTHERWISE:          OTHERWISE:
1390          (* relax free energy, but still insist T,P same over all *)              (* relax free energy, but still insist T,P same over all *)
1391         END WHEN;          END WHEN;
1392    
1393             (* Complementarity formulation *)          (* Complementarity formulation *)
1394    
1395         FOR j IN [pd.phases] CREATE          FOR j IN [pd.phases] CREATE
1396             complementarity[j]:          complementarity[j]:
1397          slack_PhaseDisappearance[j] * phase_fraction[j] = 0.0;              slack_PhaseDisappearance[j] * phase_fraction[j] = 0.0;
1398             END FOR;          END FOR;
1399    
1400      END SELECT;      END SELECT;
1401    
# Line 1481  METHOD specify; Line 1462  METHOD specify;
1462          RUN phase[pd.phases].specify_partials_seqmod_massonly;          RUN phase[pd.phases].specify_partials_seqmod_massonly;
1463      END IF;      END IF;
1464      (* The phase specify methods fix      (* The phase specify methods fix
1465       * all but one y variable per phase, and each phase's P and T.       * all but one y variable per phase, and each phase's P and T.
1466       * They also specify that the phases must not disappear.       * They also specify that the phases must not disappear.
1467       * They leave the alpha[i] free.       * They leave the alpha[i] free.
1468       *)       *)
1469      FREE phase[pd.other_phases].P;      FREE phase[pd.other_phases].P;
1470      FREE phase[pd.other_phases].T;      FREE phase[pd.other_phases].T;
1471      FREE phase[pd.phases].y[cd.components];      FREE phase[pd.phases].y[cd.components];
1472      (* Now only reference phase P,T fixed, unless a mass only MODEL. *)      (* Now only reference phase P,T fixed, unless a mass only MODEL. *)
1473    
1474      IF (number_of_other_phases != 0) THEN      IF (number_of_other_phases != 0) THEN
1475          (* allow phases to disappear *)          (* allow phases to disappear *)
1476          FREE slack_PhaseDisappearance[pd.phases];          FREE slack_PhaseDisappearance[pd.phases];
1477          slack_PhaseDisappearance[pd.phases] := 0.0;          slack_PhaseDisappearance[pd.phases] := 0.0;
# Line 1506  METHOD specify; Line 1487  METHOD specify;
1487          alpha_bar[pd.other_phases] := 1.0;          alpha_bar[pd.other_phases] := 1.0;
1488          (* When alphas are free, they are essentially K-values *)          (* When alphas are free, they are essentially K-values *)
1489          (* Solve for fixed phase fraction for equilibrium          (* Solve for fixed phase fraction for equilibrium
1490           * flash -- setting T may push flash out of multiphase region.           * flash -- setting T may push flash out of multiphase region.
1491           * Do we need to swap T,ph_fr here? Think so.           * Do we need to swap T,ph_fr here? Think so.
1492           *)           *)
1493      ELSE      ELSE
1494          (* Else solve as relative volatility model *)          (* Else solve as relative volatility model *)
1495          FIX phase_fraction[pd.other_phases];          FIX phase_fraction[pd.other_phases];
1496          FIX phase[pd.other_phases].alpha[cd.components];          FIX phase[pd.other_phases].alpha[cd.components];
1497      END IF;      END IF;
1498  END specify;  END specify;

Legend:
Removed from v.2065  
changed lines
  Added in v.2066

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