equilibrium between gas and solids using gas.equilibrate method

[Dong Wang]

Hi everyone, 

I'm using cantera to compute the equilibrium abundances of a mixture with both gas and condensed species. I found an obviously wrong result I cannot understand. So here is my simple test case. I consider water ice, liquid water and water vapor, at 1000 K, water vapor should dominate at the equilibrium state, however, my calculation shows water ice dominate. My .cti file looks like this: 

units(length = "cm", time = "s", quantity = "molec", act_energy = "K")

ideal_gas(name = "test_case", 

      elements =  "H O Ar", 

      species = """ Ar H2O H2O(s) H2O(L) """,

      initial_state = state(temperature = 300.0,

                        pressure = OneAtm)  

   )

species(name = "H2O",

    atoms = " H:2  O:1 ",

    thermo = (

       NASA( [  200.00,  1000.00], [  4.198640560E+00,  -2.036434100E-03, 

                6.520402110E-06,  -5.487970620E-09,   1.771978170E-12,

               -3.029372670E+04,  -8.490322080E-01] ),

       NASA( [ 1000.00,  6000.00], [  2.677037870E+00,   2.973183290E-03, 

               -7.737696900E-07,   9.443366890E-11,  -4.269009590E-15,

               -2.988589380E+04,   6.882555710E+00] )

             )

    # note = "L 8/89"

       )

species(name = "H2O(L)",

    atoms = " H:2  O:1 ",

    thermo = (

       NASA( [  273.15,  600.00], [  7.255750050E+01,  -6.624454020E-01, 

                2.561987460E-03,  -4.365919230E-06,   2.781789810E-09,

               -4.188654990E+04,  -2.882801370E+02] )

             )

    # note = "L 8/89"

       )

species(name = "H2O(s)",

    atoms = " H:2  O:1 ",

    thermo = (

       NASA( [  200.00,  273.15], [  5.296779700E+00,  -6.757492470E-02, 

                5.169421090E-04,  -1.438533600E-06,   1.525647940E-09,

               -3.622665570E+04,  -1.792204280E+01] )

             )

    # note = "L 8/89"

       )

species(name = "Ar",

    atoms = " Ar:1 ",

    thermo = (

       NASA( [  200.00,  1000.00], [  2.500000000E+00,   0.000000000E+00, 

                0.000000000E+00,   0.000000000E+00,   0.000000000E+00,

               -7.453750000E+02,   4.379674910E+00] ),

       NASA( [ 1000.00,  6000.00], [  2.500000000E+00,   0.000000000E+00, 

                0.000000000E+00,   0.000000000E+00,   0.000000000E+00,

               -7.453750000E+02,   4.379674910E+00] )

             )

    # note = "L 6/88"

       )

The python script looks like this: 

import numpy as np

import cantera as ct

ct.add_directory('/Users/wangdong/Documents/2015_summer/research/equilibrium_solver/debug')

gas1 = ct.Solution('test.cti')

composition = 'Ar:1.0,H2O:0.01'

gas1.X = composition

gas1.TP = 1000., 1.e5

# compute equilibrium abundances

gas1.equilibrate('TP')

gas1()

The result is 

       temperature            1000  K

          pressure          100000  Pa

           density        0.477852  kg/m^3

  mean mol. weight         39.7308  amu

 

                          1 kg            1 kmol

                       -----------      ------------

          enthalpy      4.7119e+05        1.872e+07     J

   internal energy      2.6192e+05        1.041e+07     J

           entropy          4728.6        1.879e+05     J/K

    Gibbs function     -4.2574e+06       -1.692e+08     J

 heat capacity c_p          1625.5        6.458e+04     J/K

 heat capacity c_v          1416.2        5.627e+04     J/K

 

                           X                 Y          Chem. Pot. / RT    

                     -------------     ------------     ------------

                Ar       0.990099         0.995511         -19.9176

               H2O      0.0001168       5.2961e-05         -63.0172

            H2O(s)     0.00972337        0.0044089         -63.0172

            H2O(L)    6.08195e-05      2.75776e-05         -63.0172

This is very strange because H2O(s) shouldn't be the dominate product! Please let me know where I get wrong. I appreciate it!! 

best

-Dong

[Dong Wang]

Just a follow up, I did some test, it seems that cantera automatically does extrapolation when the temperature you set is beyond the range of the temperature defined in the NASA polynomial. In this case, the extrapolation leads to wrong results. 

-Dong

[Ray Speth]

Dong,

While the extrapolation can be a problem in some cases, there's a more fundamental problem here: you are defining an ideal gas to contain species that are not ideal gases, which is nonsensical. You need to define the solid, liquid, and vapor as separate phases, of the types stoichiometric_solid, stoichiometric_liquid, and ideal_gas, respectively, and then use a Mixture object to combine them and find the equilibrium state. For example, consider the following example using the phase definitions from the included water.cti file:

import cantera as ct

liquid = ct.Solution('water.cti', 'liquid_water')
solid = ct.Solution('water.cti', 'ice')
gas = ct.Solution('h2o2.cti')

for T in [280, 1000]:
    gas.TPX = 300, 101325, 'AR:0.98, H2O:0.02'
    mix = ct.Mixture([(gas, 1.0), (liquid, 0.0), (solid, 0.0)])
    mix.T = T
    mix.P = 101325
    mix.equilibrate('TP', solver='gibbs')

    print('T = {}'.format(T))
    mix()

Which will show an equilibrium between liquid and vapor at 280 K (with the partial pressure of water in the gas phase being the vapor pressure at that temperature) and a pure gas-phase mixture at 1000 K. The other multiphase equilibrium solver (which can be specified by solver='vcs') exhibits incorrect behavior due to extrapolating the thermodynamic properties for ice.

Regards,

Ray

[Dong Wang]

Hi Ray, 

Thank you very much for your help! 

-Dong