Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Incorrect activity coefficients for molality-based phases #1311

Open
speth opened this issue Jun 4, 2022 · 0 comments
Open

Incorrect activity coefficients for molality-based phases #1311

speth opened this issue Jun 4, 2022 · 0 comments

Comments

@speth
Copy link
Member

speth commented Jun 4, 2022

Problem description

The molar activity coefficients (returned by getActivityCoefficients) for phases that use molalities internally (those derived from the MolalityVPSSTP class) are not correct.

Steps to reproduce

For the Debye-Huckel model:

>>> p = ct.Solution('test/data/debye-huckel-all.yaml', 'debye-huckel-dilute')
>>> p()

  debye-huckel-dilute:

       temperature             300  K
          pressure          101325  Pa
           density         55.9229  kg/m^3
  mean mol. weight         21.3322  amu
         potential               0  V
                pH           8.861

                          1 kg            1 kmol
                       -----------      ------------
          enthalpy    -1.22199e+07       -2.607e+08     J
   internal energy    -1.22217e+07       -2.607e+08     J
           entropy         4194.41        8.948e+04     J/K
    Gibbs function    -1.34782e+07       -2.875e+08     J
 heat capacity c_p          158171        3.374e+06     J/K
 heat capacity c_v    <not implemented>

                           X           Molalities         Chem.Pot.    ChemPotSS    ActCoeffMolal
                                                          (J/kmol)      (J/kmol)
                     -------------    ------------     ------------  ------------    ------------
            H2O(L)      0.738115          55.5093     -3.07518e+08   -3.06815e+08        1.02227
               Na+      0.124393           9.3549     -2.57445e+08   -2.57956e+08       0.131209
               Cl-      0.124393           9.3549     -1.83605e+08   -1.84117e+08       0.131209
                H+   1.39607e-10       1.0499e-08     -5.08919e+07              0       0.131209
               OH-   1.83035e-08       1.3765e-06     -2.65584e+08   -2.26855e+08       0.131209
          NaCl(aq)     0.0130966          0.98492     -4.32958e+08    -4.3292e+08              1
          NaOH(aq)   5.16408e-08       3.8836e-06     -5.14827e+08   -4.83751e+08              1
          SiO2(aq)   4.01294e-07       3.0179e-05     -9.25875e+08   -8.99913e+08              1
      NaH3SiO4(aq)   9.14817e-07       6.8798e-05     -1.74563e+09   -1.72172e+09              1
           H3SiO4-   1.36043e-08       1.0231e-06     -1.50517e+09    -1.4657e+09       0.131209

>>> p.activities
array([7.54555306e-01, 1.22745026e+00, 1.22744995e+00, 1.37756687e-09,
       1.80609658e-07, 9.84920000e-01, 3.88360000e-06, 3.01790000e-05,
       6.87980000e-05, 1.34240277e-07])

>>> 0.131209 * 9.3549  # molal activity coefficient times molality for Na+
1.2274470741  # matches reported activity

>>> p.activity_coefficients * p.X  # should equal the activities, but doesn't
array([7.54555306e-01, 2.21125165e-02, 2.21125109e-02, 2.48168672e-11,
       3.25368299e-09, 1.77433338e-02, 6.99630540e-08, 5.43674685e-07,
       1.23939597e-06, 2.41833859e-09])

Likewise for IdealMolalSoln:

>>> p = ct.Solution('test/data/thermo-models.yaml', 'ideal-molal-aqueous')
>>> p()

       ideal-molal-aqueous:
     
            temperature          298.15  K
               pressure          101325  Pa
                density         12.0577  kg/m^3
       mean mol. weight           18.07  amu
              potential               0  V
     
                               1 kg            1 kmol
                            -----------      ------------
               enthalpy     0.000735017          0.01328     J
        internal energy        -8403.38       -1.518e+05     J
                entropy         5.18177            93.63     J/K
         Gibbs function    -2.15747e+06       -3.899e+07     J
      heat capacity c_p          1595.8        2.884e+04     J/K
      heat capacity c_v    <not implemented>
     
                                X           Molalities         Chem.Pot.    ChemPotSS    ActCoeffMolal
                                                               (J/kmol)      (J/kmol)
                          -------------    ------------     ------------  ------------    ------------
                 H2O(l)      0.997484          55.5093     -3.89686e+07   -3.89623e+07       0.999997
                CO2(aq)    0.00179697              0.1     -4.46703e+07   -3.89623e+07              1
                H2S(aq)    0.00053909             0.03     -4.76549e+07   -3.89623e+07              1
                CH4(aq)   0.000179697             0.01     -5.03783e+07   -3.89623e+07              1

>>> p.activities
array([0.99748108, 0.1       , 0.03      , 0.01      ])  # equals molal activity coeff times molality, except for solvent

>>> p.activity_coefficients
array([0.99999682, 1.0025221 , 1.0025221 , 1.0025221 ])  # not just returning molality based coefficients

>>> p.activity_coefficients * p.X
array([9.97481078e-01, 1.80150000e-03, 5.40450000e-04, 1.80150000e-04])  # should equal activities, but don't

System information

  • Cantera version: 2.6.0 or main at 1ab81ac

Additional context

This was discovered as part of #1299, which implements Cantera/enhancements#114.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
No open projects
Development

No branches or pull requests

1 participant