Pressure used in equilibrium constant

[JB Scoggins]

Hello,

I have a question concerning how the equilibrium constants for each reaction are computed.  It is my understanding that the equation for Kc should be 

K_c = \frac{P_atm}{Ru T}^(\Delta \nu) exp(- \frac{\Delta G}{Ru T})

where the P_atm term represents the standard atmosphere (101325 Pa).  However, in looking at the actual source code of Cantera, it appears that the total pressure of the mixture is used in place of P_atm.  Is this correct?  To be honest, I have seen it both ways in various chemical kinetics texts but I thought that P_atm is the correct pressure because the term \frac{P_atm}{Ru T}^(\Delta \nu) is in place to convert the units of K_p to concentration based units for K_c.  It is also my understanding that K_c is only a function of temperature and thus the mixture pressure should not show up in the expression.  I would really appreciate it if someone could clarify this issue for me.

Warm Regards,

JB Scoggins

[Ray Speth]

Hi JB,

You're certainly right that K_c should be a function of temperature only. At first glance, it does look like Cantera is using the mixture pressure when calculating K_c, but it turns out not affect the resulting equilibrium constants. For example, I tried this in Python (the reaction index is for H + H + H2 <-> H2 + H2, i.e. one where the number of moles changes):

In [1]: import Cantera as ct

In [2]: gas = ct.GRI30()

In [3]: gas.set(T=1700, P=101325, X='H2:1.0, N2:3.0')

In [4]: gas.equilibriumConstants()[39]

Out[4]: 6418881501.4493847

In [5]: gas.set(T=1700, P=101325*5, X='H2:1.0, N2:3.0')

In [6]: gas.equilibriumConstants()[39]

Out[6]: 6418881501.4493618

I'm not entirely sure, but I think Cantera does have a problem here, even though the correct equilibrium constants are being computed.

The function GasKinetics::getEquilibriumConstants uses Gibbs free energies computed using the function IdealGasPhase::getStandardChemPotentials in a formula for Kc that involves the mixture pressure (but shouldn't). However, the getStandardChePotentials function actually calculates the chemical potential at the current pressure. Using this value exactly cancels out the effect of using the mixture pressure in the formula for Kc.

I believe that IdealGasPhase::getStandardChemPotentials and IdealGasPhase::getChemPotentials are both wrong (with the latter effectively adding the RT*log(P) term twice), but I'd appreciate it if someone else could confirm this interpretation before I make any changes to the code.

Ray

[Moffat, Harry K]

Hi,

 

  The function getEquilibriumConstants() from InterfaceKinetics returns an equilibrium constant that is pressure dependent in the general case, but is pressure independent for reactions involving ideal gases. This is consistent with what’s done with solution thermodynamics (see Denbigh chaps 9 and 10). getEquilibriumConstants() doesn’t return K_p nor K_c per se, because these really aren’t defined well or aren’t meaningful for general non-ideal phases, which may be gases, solutions, solids, or a combination of those.

 

Cantera’s implementation of equilibrium is consistent with the activity coefficient model for presentation of chemical potentials with additions to also make it consistent with the kinetics operator. Let me  explain. I’ll leave out complications due to charge transfer reactions. The chemical potential of a species can be separated into the standard state chemical potential and the activity.

 

   chemPot_i(T, P, X_j) =  chemPot_i^ss(T,P)  + RT ln (a_i)

 

 

chemPot_i^ss(T,P)  is the standard state chemical potential of species i at the solution T and P. 

a_i is the activity of species i.

 

For a single reaction the sum of the product of the chemical potentials multiplied by the stoichiometric coefficient must be zero. The equilibrium constant returned by getEquilibriumConstants() is based on the sum of the product of [chemPot_i^ss(T,P) - RT ln(co_i)]    and the stoichiometric coefficients.

 

For ideal gases,

        chemPot_i^ss(T,P)  = chemPot_i^0(T,Po)  + RT ln (P / Po)

 

where Po is the reference pressure. Usually 1 bar is used in the JANAF tables. In other words, the standard state for gases is defined as a pure ideal gas species at the temperature and pressure of the solution. The reference state for gases is defined as a pure ideal gas species at the reference pressure of 1 bar and the temperature of the solution. (it’s the reference state thermo that can be looked up in the JANAF tables, btw).

 

For mass action kinetics, we multiply the activity by the standard concentration for species i to obtain an activity concentration, which is then used in kinetics expressions involving the phase.

 

    ac_i = (co_i) (a_i)

 

For gases, the standard concentration is defined as P/RT   , where P is the system pressure and T is the system temperature.  Note, for solids and liquids and surface ligands, usually the standard concentration is defined as unity, however. Then the activity concentration for species i, ac_i, is defined for gases as

 

   ac_i =    P X_i gamma_i   / RT

 

which is basically the concentration of gases, which is what you want for gas-phase mass action kinetics. (gamma_i is the activity coefficient of species i, which for ideal gases is equal to one)

 

Therefore, we can manipulate terms to get:

 

   chemPot_i(T, P, X_j) =  chemPot_i^ss(T,P)  + RT ln (ac_i) – RT ln(co_i)

 

 

For ideal gases, the terms left over after separating out the concentration looks a lot like the input to the concentration equilibrium coefficient.

 

         chemPot_i^ss(T,P)  – RT ln(co_i) =  chemPot_i^o(T,Po)  + RT ln (P / Po) - RT ln (RT / P)

                                                                           =  chemPot_i^o(T,Po)  + RT ln (RT / Po)

 

This is what’s used in the equilibrium constant. Therefore, as Ray’s example shows the equilibrium constant returned by Cantera is independent of P for gas phase reaction that are not equimolar, but only because there is a fortuitous cancelling of terms for the ideal gas case.

 

For non-ideal solutions (e.g., brines at elevated temperatures and pressures) the situation gets complicated, and there is no other general starting point possible other than the activity coefficient representation of the chemical potential. In particular reactions which lower the molar volumes are constituents are favored when the pressure goes up. This effect is captured by the pressure dependence of the standard state gibbs free energies of the species and is reflected in a change in the equilibrium constant as a function of pressure. Cantera’s implementation has been checked against multiple problems and is quite solid.

 

In particular IdealGasPhase::getStandardChemPotentials() adds the + RT ln (P / Po) term as appropriate to change the chemical potential from the reference pressure to the current pressure. IdealGasPhase::getChemPotentials() just adds the RTln(X_i) term to the resulting standard state chemical potential.

 

 

Best wishes,
Harry

--
You received this message because you are subscribed to the Google Groups "Cantera User's Group" group.
To view this discussion on the web visit https://groups.google.com/d/msg/cantera-users/-/SoEsUaxx-30J.
To post to this group, send email to canter...@googlegroups.com.
To unsubscribe from this group, send email to cantera-user...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/cantera-users?hl=en.

[JB Scoggins]

Hi Harry,

Thank you very much for your detailed reply.  I assumed Cantera is well vetted in this case but my own understanding was lacking.  My background is in Aerospace Engineering so obviously I only had training in gas phase kinetics.  Your discussion has cleared up my confusion though and now the Cantera source code is much more clear.

Thank you,

JB

[Manqi Zhu]

Hi Harry,

Thank you for your explanation about how Cantera cacultates the equilibrium constant in a reversible reaction but I am still a little confused, could you please help me to better understand it?

In fact, I am working on the reversible reaction 2C2H5 <=> C4H10 in a chemical schema containing 149 reactions.

I have run Cantera for 1 iteration and output the variable like dH_0, dS_0, dG_0, Equilibrium Constants, etc...

I am able to output the dG_0 by writing delta_G0() in Cantera, which given -1.6984417E+8 J/kmol, correspondent to what I get by myself (calculate manually dG_0=dH_0-TdS_0 using the thermal data of C2H5 and C4H10, and T is the initial temperature that I have setted as 1100K).

Then, I have first used this equation for calculating manually the equilibrium constant Kp:

Kp=exp(-dG_0/RT) 

where R=8314 J/kmol is the gas constant, T is my initial temperature setted as 1100K, I get manually Keq=1.16E+8.

Cantera gives me 4.858E+9 by writing equilibriumConstants(), which is said to be the "Equilibrium constants in concentration units for all reactions." on the website http://cantera.github.io/docs/sphinx/html/python/kinetics.html

So it should be Kc=k_fwd/k_rev to me, or something not far away...

So, according to Kuo (Principle of Combustions):

Kc=Kp*(P_0/RT)^(mu_{C4H10}-mu_{C2H5}) 

where P_0 is the reference pressure usually equal to 100 000Pa, mu_{C4H10}=2 and mu_{C2H5}=1 are the stoechiometric coefficient of C4H10 and C2H5 respectively. 

Finally I get manually Kc=1.06E+7… Comparing to 4.858E+9 in Cantera I am really confused…

From what you have explained, I know that Cantera does not calculate Kp or Kc at all… but I couldn't understand how it calculates exactly by the equation you have written:

chemPot_i^o(T,Po) + RT ln (RT / Po)

The chemPot_i^o(T,Po) is equal to the Gibbs free energy dG_0?

Could you explain more to me about that? And if you could provide me an equation like Kc=exp(-dG_0/RT)*(P_0/RT)^(mu_{C4H10}-mu_{C2H5}) allowing me to recalculated manually the equilibrium constant as in Cantera, that will help me a lot to understand.

 

Thanks in advance!

Manqi

[Ray Speth]

Hi Manqi,

I'm not sure what thermodynamic data you're using, but it doesn't give results that agree with what I get using data from one of the LLNL mechanisms. The only catch in these calculations, as far as Cantera goes, is the meaning of 'standard', which in some literature implies a value at the reference pressure, but in Cantera means a value at the current pressure that only excludes the specific composition effects. I have attached a sample input file which goes with the following script (written for the 'new' Python module) that should show you that all the methods for calculating the equilibrium constant are equivalent and give consistent results using Cantera:

import cantera as ct

import numpy as np

gas = ct.Solution('sample.cti')

gas.TPX = 1100, 5e5, 'c2h5:0.9, c4h10:0.1'

gas.equilibrate('TP')

p_c4h10, p_c2h5 = gas.P * gas['c4h10', 'c2h5'].X

RT = ct.gas_constant * gas.T

p0 = 1e5

nu_net = 1

# Kc

Kc_1 = gas.equilibrium_constants[0]

Kc_2 = ((p_c4h10/RT)**1) * ((p_c2h5/RT)**-2)

print 'Kc: {0:e}, {1:e}'.format(Kc_1, Kc_2)
# = 2.038503e+09

# Kp

Kp_1 = ((p_c4h10/p0)**1) * ((p_c2h5/p0)**-2)

# need delta G at p0

gas.TP = None, p0

Kp_2 = np.exp(-gas.delta_standard_gibbs[0] / RT)

Kp_3 = Kc_1 * (p0/RT)**nu_net

print 'Kp: {0:e}, {1:e}, {2:e}'.format(Kp_1, Kp_2, Kp_3)
# 2.228868e+07

Please let me know if you have any further questions.

Regards,

Ray

[Ray Speth]

Sorry, I forgot the attachment.

Ray

[Manqi Zhu]

Thank you so much Ray!!

With your explication and your script, now I much better understand how Cantera calculate the equilibrium constant.

On the other hand, one of my colleagues also showed me the approach in Cantera, I'd like to share it here:

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

In Cantera, in the formula Keq=exp(-dG_0/RT), with dG_0=dH_0-TdS_0, but the dS_0 is defined as the dS at pressure P instead of P_0 (we use usually the reference pressure here), where P is the real pressure of the system, as Ray said.

dH only depends on temperature but not on pressure so there will be no difference.

dS depends on temperature and pressure, so dS_0(P)=dS_0(P_0)-R*ln(P_0/P)^(nu_net)

First, Kp=exp(-dG_0/RT)=exp(-(dH_0-TdS_0(P))/RT) = exp(-(dH_0-TdS_0(P_0)+T*R*ln(P_0/P)^(nu_net))/RT) = exp(-(dH_0-TdS_0(P_0))/RT) * ((P_0/P)^(nu_net))

Then Keq=Kc=Kp*(P/RT)^(nu_net) 

= exp(-(dH_0-TdS_0(P_0))/RT) * ((P_0/P)^(nu_net)) * (P/RT)^(nu_net) 

= exp(-(dH_0-TdS_0(P_0))/RT) * (P_0/RT)^(nu_net)

here we refind the usual definition of the Kc with the reference pressure.

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

PS to Ray, the thermo data I used was from UGent, especially for a butane cracking reduced chemistry schema, maybe there are some difference from the LLNL data you used.

Thanks a lot!

--
You received this message because you are subscribed to a topic in the Google Groups "Cantera Users' Group" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/cantera-users/65bMdBD0dAI/unsubscribe.
To unsubscribe from this group and all its topics, send an email to cantera-user...@googlegroups.com.

To post to this group, send email to canter...@googlegroups.com.

Visit this group at http://groups.google.com/group/cantera-users.
For more options, visit https://groups.google.com/d/optout.