Thermodynamics with ASCEND

The ASCEND ModelLibrary includes a range of models for dealing with chemical and physical properties of materials, phase equilibrium, and reaction kinetics.

We have a brief 'HOW TO' guide for doing thermodynamics calculations with ASCEND, thanks for Pedro Vale Lima, here: howto_thermo.pdf. The following text is our partial (ongoing) transcription of that document into wikitext. The same ASCEND thermodynamics library is also covered in Chapter 11 of the ASCEND user's manual.

In addition to the methods described in this page, there are other ways of calculating thermodynamic properties using ASCEND. See FPROPS and freesteam for example.

You can access the models from the HOW TO in the ModelLibrary as models/thermodynamics_example.a4c.

Physical property data

Physical property data from the book Reid, Prausnitz & Polling, The Properties of Gases and Liquids, 4th Ed.,McGraw-Hill, 1987, are included in the ASCEND model file models/components.a4l.

To add support for a new fluid component, it is suggested that you create your own copy of that file, and add some new code. First add your component to the supported_components list:

supported_components :== [
   'hydrogen'
   ,'carbon_dioxide'
   ,(* add another one here *)
];

Then, add a CASE in the SWITCH statement to define the constants for your particular substance:

CASE 'my_new_substance':
   rpp_index :== 111;

   (* add your data here... *)

Here is the meaning of the constants that are used in the components models:

constant	meaning
rpp_index	Index in the Reid et al book
formula	Chemical formula (optional, informational only)
subgroups	Component subgroups for the UNIFAC model.
wilson_set	Wilson groups for the Wilson liquid model.
mw	Molecular weight
Tb	Boiling-point temperature (at what pressure??)
Tc	Temperature at critical point
Pc	Pressure at critical point
Vc	(Molar?) volume at critical point
Zc	Compressibility factor at critical point.
omega	Pitzer acentric factor
cpvapa..d	Constants for vapour-phase <math>c_p</math> calculation. <math>c_p = a + bT +c T^2 + d T^3</math>.
Hf	Enthalpy of formation
Gf	Gibbs free energy of formation
vpa..d	Constants for vapour pressure calculation (three equations can be used)
vp_correlation	Number that represents which correlation is used for calculation of vapour pressure
Hv	Enthalpy of vapourisation
lden	Liquid density
Tliq	Temperature at which the liquid density (lden) is specified.

Finally, add the name of your new component in the check list of components model:

FOR i IN components DO

   ASSERT (i IN ['hydrogen', 'carbon_dioxide' (* others here *)) == TRUE;

END FOR;

Calculating properties for a single component

The simples state to model is the one that represents a single component in a single phase. To elaborate this model, one needs an equation of state (EOS). As there is no universal equation of state, the thermodynamic libraries in ASCEND contain several models, and leave the decision about which one to use up to the user.

The ASCEND thermodynamics library has two models for vapour components (ideal_vapour_component and Pitzer_vapor_component) and one model for liquid components (Rackett_liquid_component). All these models are similar in concept; only the EOS changes. We'll use the ideal vapour model to illustrate the main concepts.

The model defines the equations needed to calculate the molar volume <math>{\bar v}</math>, enthalpy <math>{\bar h}</math> and Gibbs free energy <math>{\bar g}</math>, based on the temperature T, pressure p and component properties.

The first equation needed by the ideal gas model is the well known equation of state,

([[#equation_{{{3}}}|{{{3}}}]])

1

where <math>{\bar R}</math> is the ideal gas constant. To define <math>{\bar h}</math>, one can use the constant-pressure molar heat capacity, <math>{\bar c}_p</math>. The heat capacity is define implicitly as a polynomial function of temperature:

([[#equation_{{{3}}}|{{{3}}}]])

2

where a, b, c, and d are defined (as cpvapa..d) in the components data file. Remembering that, for an ideal gas, <math>{\bar h}</math> is defined as

([[#equation_{{{3}}}|{{{3}}}]])

3

results in the following explicit equation for <math>{\bar h}</math>:

([[#equation_{{{3}}}|{{{3}}}]])

4

Note that, if it were not an ideal gas, <math>{\bar h}</math> would also depend on p. Looking back again to your old thermodynamics textbook, it can then be proven that

([[#equation_{{{3}}}|{{{3}}}]])

{{{2}}}</math>. To make it easier to understand the goal of this representation let's see what these equations represent for a two phase system (V,L) of the two components (A,B).

([[#equation_{{{3}}}|{{{3}}}]])

24

([[#equation_{{{3}}}|{{{3}}}]])

25

in the relative volatility model you fix the <math>\alpha</math> values and compute the <math>\bar \alpha</math>. Note that <math>\alpha_A^V / {\bar \alpha}</math> will be a K-value. Before seeing the model in action let's do some accounting. To simplify, we suppose a two-phase system, since ASCEND is not yet able to compute liquid-liquid models.

FIXME: add the equation-counting table here.

Chemical potential equilibrium

In the full equilibrium model, the Gibbs free energy for each phase is equated. To choose full equilibrium, you should set the variable equilibrated to TRUE.

([[#equation_{{{3}}}|{{{3}}}]])

{{{2}}}</math>|26}}

for each partial component and each phase except the reference one.

When using this model <math>\bar \alpha</math> is fixed and equal to one and the <math>\alpha</math>'s are released and computed by equation 23. This equation becomes a computation of K-values. Now let's see the infamous accounting table for this model.

FIXME equation-counting for VL-equilibrium G-model.

Now let's look at some more examples.

Example 4

Calculate the change of enthalpy for a 50% water, 50% ethanol being heated from 300 K to 400 K at constant atmospheric pressure.

This is a simple problem but it's useful for introducing the correct way to use the thermodynamics model.

MODEL howto_thermo_ex8 REFINES cmumodel;

   cd IS_A components_data(['water','ethanol'], 'water');
   pd IS_A phases_data('VL', 'ideal_vapor_mixture', 'UNIFAC_liquid_mixture',

                       'none');
   equilibrated IS_A boolean;
   phases ALIASES pd.phases;

   FOR j IN phases CREATE
      smt[j] IS_A select_mixture_type(cd, pd.phase_type[j]);

   END FOR;
   FOR j IN phases CREATE
      phase[j] ALIASES smt[j].phase;

   END FOR;

   state IS_A thermodynamics(cd, pd,phase, equilibrated);

METHODS
(* ... *)
   METHOD values;
      state.P := 1 {atm};

      state.T := 300 {K};
      state.y['ethanol'] := 0.5;

      equilibrated := TRUE;
   END values;
(* ... *)
END howto_thermo_ex8;

Use of the thermodynamics model is slightly more complex. First, you should define the components data, then you choose the models for the phases with the phases_data model. This model's first parameter is a phase definition (it can be 'V', 'L' or 'VL'), then the name for the vapour model (or none if there's not vapour phase), and then the liquid model name. The last parameter is for a second liquid phase model when LL equilibrium is implemented; for now, you should just use none.

The select_mixture_type model builds the chosen models so that they can be sent to the main thermodynamic model as a parameter. This may seem complicated at first, but as you probable will never need to change that piece of code, you don't need to worry about it too much. Just peek at the library code when you feel brave enough.

The value definitions are very simple; just notice the newly introduced equilibrated variable.

When solving multi-phase models you should always check state.slack_PhaseDisappearance. Only if it's zero the phase exists. Sometimes you can get strange results if you're looking at properties of a non-existent phase.

If you run the model, you'll get only a liquid phase and <math>\bar h</math> = -282.351 kJ/mol. With <math>T</math> = 400 K, you'll get <math>\bar h</math> = -233.022 kJ/mol, and so <math>\Delta \bar h</math> = 49.329 kJ/mol.

Example 5

Build the VL equilibrium chart for water-ethanol at atmospheric pressure.

The model presented in the last example specifies <math>p</math> and <math>T</math>. To build an equilibrium chart, it's better to specify <math>p</math> and a mole fraction, and to free <math>T</math>. That way, we can be sure to have an equilibrium. The following method, added to the last example model, uses that trick. We'll also use the Pitzer vapour model, in the hope of improving accuracy.

MODEL howto_thermo_ex9 REFINES cmumodel;
   cd IS_A components_data(['water','ethanol'], 'water');

   pd IS_A phases_data('VL', 'Pitzer_vapor_mixture', 'UNIFAC_liquid_mixture',
                       'none');

   equilibrated IS_A boolean;
   phases ALIASES pd.phases;
(*...*)

   METHOD reset_Px;
      equilibrated := TRUE;
      RUN state.specify;

      state.T.fixed := FALSE;
      state.phase['liquid1'].y['water'].fixed := TRUE;

      state.P := 1 {atm};
      state.y['ethanol'] := 0.5;

      state.phase['liquid1'].y['water'] := 0.6;
   END reset_Px;

END howto_thermo_ex9;

Now, instead of calling reset and values just run reset_Px. We'll run the model several times and build a table like the following:

<math>y_w^L</math>	<math>y_w^V</math>	T / [K]
0.99	0.8786	369.764
0.95	0.652	362.471
0.9	0.5517	358.80
:	:	:

When building the table, you must take some care for a phase not to disappear. You can do that by changing the variable state.phase['liquid1'].y['water'].