Real Gas Equation of state implementation into thermodynamic properties in 1D flame simulations

1,705 views
Skip to first unread message

Refik Alper Tuncer

unread,
Oct 8, 2017, 8:57:08 AM10/8/17
to Cantera Users' Group
Hi All,

I have several questions regarding the real gas implementations into Cantera.

I need to implement a real gas equation of state instead of using ideal gas. ( Only need to change the thermodynamic properties such as enthalpy, internal energy vs, transport properties are not important at first stage )

Is it possible to do this in C++ version of Cantera ? Or MATLAB ? 

From what i remember, 1D flames only works with ideal gas. Is it still that way ?

I would like to start coding at once if real gas implementation is possible through C++ or MATLAB.

Best Regards.

Alper Tuncer

Ray Speth

unread,
Oct 9, 2017, 6:02:14 PM10/9/17
to Cantera Users' Group
Alper,

Cantera does have some real gas capabilities, but this is an area where there is still some work to be done. You can see an example of a reactor simulation where real gas effects are accounted for at https://github.com/Cantera/cantera-jupyter/blob/master/reactors/NonIdealShockTube.ipynb.

In the 1D flame model in Cantera, however, the energy equation is written in terms of the temperature assuming that ideal gas relationships between the enthalpy and temperature hold. You would need to rederive and reimplement the energy equation in a more general form in order to use the solver with a real gas model (and I think this would be a useful addition to Cantera if someone wanted to undertake this implementation).

Regards,
Ray

Chris Neal

unread,
Mar 20, 2018, 5:58:47 PM3/20/18
to Cantera Users' Group
Hey Ray,

I've been studying the relevant portions of the 1D solver with respect to the energy equation. With Steven D. working on the Peng-Robinson EOS, I think it's time to try to generalize that energy equation. I've done some working with trying to start with Kee's generalized equations and get down to the temperature form. I also plan to compare the final form to what Poinsot & Veynante derive for their general temperature form of an energy equation. I'll try to get some Latex formulas posted in this thread tomorrow to get your thoughts on the form the equation & how to best go about adjusting the current implementation of the energy equation.

-Chris N.

Chris Neal

unread,
Mar 21, 2018, 6:00:54 PM3/21/18
to Cantera Users' Group
So Kee starts with the following equation(3.203) in the 1st edition of Chemically Reacting Flow.

Equation:  \rho \frac{Dh}{Dt} = \frac{Dp}{Dt} + \nabla \cdot (\lambda \nabla T) - \sum_{k=1}^{K}\nabla \cdot h_k \vec{j}_k + \Phi

Would you agree that this is an appropriate starting point? 

I understand this enthalpy to be the sum of the sensible and chemical energy based off (eq 1.60) in Poinsot's book(below). Where Kee has neglected volumetric heating, body forces, and as expanded the heat flux vector into two components of a Fourier's law heat flux and a term for the diffusion of species with different enthalpies.

Equation:  \rho \frac{Dh}{Dt} = \frac{Dp}{Dt}  -\nabla \cdot \vec{q} + Q + \Phi  + \rho\sum_{k=1}^{K}Y_kf_{k,i}V_{k,i}


Poinsot's heat flux vector looks like the following.

Equation:  q_i = -\lambda \nabla T +  \rho\sum_{k=1}^{K} h_k Y_k V_{k,i}


It looks like Kee chose to keep the mass flux generalized to the variable j, whereas Poinsot has made the following substitution: j_i = rho * Y_k * V_k,i . 

Equation:  Kee \rightarrow  \sum_{k=1}^{K} h_k \vec{j}_i   .......  Poinsot \rightarrow   \sum_{k=1}^{K} h_k \rho Y_k \vec{V_{k}} 



I don't think this is something that we should do, and Kee probably has the more general form of the energy equation in that respect.

Steven DeCaluwe

unread,
Mar 21, 2018, 6:07:14 PM3/21/18
to canter...@googlegroups.com
Hi Chris,

That is correct. Keeping it in terms of a mass flux times the per-unit-mass enthalpy is advantageous, in that mass is a conserved quantity, and so this lets one keep the conservation equations more generally applicable. Decomposing it into a function of velocity assumes that you have some method of solving for the velocity field (or are at least interested in doing so) which may not always be the case. 

Best,
Steven

Sent from my iPhone
--
You received this message because you are subscribed to the Google Groups "Cantera Users' Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cantera-user...@googlegroups.com.
To post to this group, send email to canter...@googlegroups.com.
Visit this group at https://groups.google.com/group/cantera-users.
For more options, visit https://groups.google.com/d/optout.

Ray Speth

unread,
Mar 21, 2018, 8:12:21 PM3/21/18
to Cantera Users' Group

Chris,

One question that arises for me, even from the first equation, is what is the definition of h_k?. In a non-ideal mixture, it is not necessarily true that h = \sum h_k Y_k with h_k being the pure species enthalpies. However, it looks at first glance that this idea is essentially built into this equation where the h_k j_k term first appears (Eq. 3.159 of Kee). Perhaps there is a definition for h_k for which this holds, but I’m not sure what it is.

Regards,
Ray

To unsubscribe from this group and stop receiving emails from it, send an email to cantera-users+unsubscribe@googlegroups.com.

Steven DeCaluwe

unread,
Mar 22, 2018, 12:51:09 PM3/22/18
to <cantera-users@googlegroups.com>
Hi Ray,

I do not necessarily see a problem with h = sum (h_k*Y_k).  By my read, h_k is not the pure species enthalpy, but rather includes interaction potential terms, no?  This would be equivalent to the partial molar enthalpies, divided by the species molecular weight.  For an ideal gas they are obviously the same, but for RedlichKwongMFTP (lines 292-335) for example, they include interaction terms from the EoS.

Haven’t chased down the flame code to see what is actually called, but another option is to call the mixture-averaged enthalpy (again, in RedlichKwongMFTP it is enthalpy_mole(), which one would divide by the mmw), which also contains the interaction potentials.

Thanks,
Steven




To unsubscribe from this group and stop receiving emails from it, send an email to cantera-user...@googlegroups.com.

Chris Neal

unread,
Mar 22, 2018, 1:22:38 PM3/22/18
to Cantera Users' Group
That is a great point Ray. I believe the assumption of an ideal mixture(not the same as ideal gas assumption) is common. It's a neglection of the heat of mixing. So in this case we would be deriving an energy equation for real gases that are assumed to behave like ideal solutions. I think we would be essentially stuck if we went down the path of non-ideal solutions because the enthalpy of mixing is a very complex parameter to estimate and is often empirically tabulated for a finite set of species.


To unsubscribe from this group and stop receiving emails from it, send an email to cantera-user...@googlegroups.com.

Steven DeCaluwe

unread,
Mar 22, 2018, 2:53:20 PM3/22/18
to <cantera-users@googlegroups.com>
Again, I would disagree.  Unless the text says something explicit to define h_k to the contrary, h_k is just the enthalpy per unit mass of that species *at the current state*.  It is derived based on the equation of state, is a function of T, P, and Y_k (or other state variables), but is not the same as the pure species enthalpy.  

Chris Neal

unread,
Mar 22, 2018, 3:30:07 PM3/22/18
to Cantera Users' Group
Kee's book says the following before making the mixture enthalpy definition:

"The thermodynamic properties of the mixture that appear in the energy equation are given as mass-weighted averages of the individual species properties. Specifically of interest here is the enthalpy" 


It doesn't give a reference, but if you want to look it up, it's on page 114 of the 1st edition.

Refik Alper Tuncer

unread,
Mar 22, 2018, 3:37:30 PM3/22/18
to canter...@googlegroups.com
Dear Chris,

1. Take a look at the paper in the attachment. It provides the derivation of partial-mass species enthalpy values as the summation of ideal values and the departure functions.
2. Using enthalpy in energy equation instead of temperature is very wise because it will be less expensive computationally.

Best Regards.

Alper

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/96myuNgXN-o/unsubscribe.
To unsubscribe from this group and all its topics, send an email to cantera-users+unsubscribe@googlegroups.com.
A UNIFIED TREARMENT OF GENERAL FLUID THERMODYNAMICS.pdf

Ray Speth

unread,
Mar 22, 2018, 4:25:42 PM3/22/18
to Cantera Users' Group

Steven,

Perhaps it’s better to think of this in the reverse direction. That version of the energy equation requires a definition of the individual species enthalpies for which h=\sum X_k \hat{h}_k / \bar{W}. However, if you only have an equation of state which defines h=\sum X_k \hat{h}^0_k / \bar{W} + h_{mix}, then you’re in a bit of trouble. This is in fact how the Redlich-Kwong method in Cantera appears to be implemented (i.e. see the definitions of enthalpy_mole and hresid). However, based on a quick test it seems that the definition of partialMolarEnthalpies does apportion the enthalpy of mixing to the individual species in the desired way, even though this is not obvious from looking at the code. So as long as you’re working an equation of state where this is implemented this way, it should be fine, but I do think it is something that needs to be kept in mind.

Regards,
Ray

To unsubscribe from this group and stop receiving emails from it, send an email to cantera-users+unsubscribe@googlegroups.com.

Steven DeCaluwe

unread,
Mar 22, 2018, 4:43:01 PM3/22/18
to <cantera-users@googlegroups.com>
Hi Ray,

Yes, the second definition would get you into trouble, although I would contend that an EoS which does not give you access to the partial molar or mass properties for individual species is going to cause problems in a range of places, in terms of capabilities that a user would want and in terms of solving the energy equation in a reacting flow.  

In particular (assuming I’m thinking about this correctly): if you have diffusive fluxes, you need the individual h_k values.  A mixture averaged enthalpy will not be sufficient.

While I agree that it bears keeping in mind, I would contend that it sets a standard for what is considered a fully-implemented EoS, rather than indicating some problem with the governing equations of the 1D solver.

Best,
Steven

Chris Neal

unread,
Mar 22, 2018, 5:21:13 PM3/22/18
to Cantera Users' Group

Assuming I can set aside the issue of the definition of h_k. Would it be impractical/impossible to use an enthalpy-form of an energy equation? I imagine that this would complicate the boundary condition specifications. It seems like both Poinsot and Kee resort to making a perfect gas assumption when writing the total derivative of enthalpy in terms of temperature.

Equation:  \rho \frac{Dh_s}{Dt} = \sum_{k=1}^{K}h_{s,k} \rho \frac{DY_k}{Dt} + \rho C_p \frac{DT}{Dt}


Where the following definition is used to expand the total derivative.

Equation:  h_s =  \sum_{k=1}^{K}h_{s,k} Y_k 



The derivation chain starts on equation 3.206 in Kee's book. for reference. All of this is done to get to a temperature form of the sensible energy equation. In section 3.10.5 Kee derives an energy equation for a general EOS( a bulk modulus, beta, term appears), but this is only for a single species. And I'm not sure how such a derivation would be made for multi-species.

Steven DeCaluwe

unread,
Mar 22, 2018, 5:36:20 PM3/22/18
to <cantera-users@googlegroups.com>

Hi Chris,

I don’t have the text (an oversight I need to correct), but: are you saying that the equation 


makes a perfect gas assumption?  The only problem I can see is in the first term on the RHS:  D(h_{s,k}Y_k)/Dt is not necessarily equal to  h_{s,k} D(Y_k)/Dt.  Because h_k is a function of Y_k, a finite D(Y_k)/Dt also implies a finite D(h_{s,k})/Dt. Is this what you’re referring to?

Steven



On Mar 22, 2018, at 3:21 PM, Chris Neal <chri...@snumerics.com> wrote:

Assuming I can set aside the issue of the definition of h_k. Would it be impractical/impossible to use an enthalpy-form of an energy equation? I imagine that this would complicate the boundary condition specifications. It seems like both Poinsot and Kee resort to making a perfect gas assumption when writing the total derivative of enthalpy in terms of temperature.


Chris Neal

unread,
Mar 22, 2018, 5:45:23 PM3/22/18
to Cantera Users' Group

I'm just trying to wrap my head around this enthalpy issue. For a real fluid, the thermodynamic property of enthalpy can be split into an ideal component and a departure from that ideal value, which together give the correct real-fluid enthalpy.

Equation:  h_k = h^{ig}_{k} + h_{k, departure}


For a mixture of different species the mixture enthalpy can be expressed as:

Equation:  h_{mixture} =   \sum_{k = 1}^{K} x_k h_k  + \Delta h_{mixing} 



Where the enthalpy of mixing accounts for non-ideal behaviors due to things like species attractions, etc. An EOS object like a cubic equation of state would be able to give you the first estimate for the enthalpy of a species, or it's excess enthalpy, and maybe any other thermodynamic properties that can be expressed like that(internal energy, etc).  And an EOS object might or might not return the mixing enthalpy for the mixture enthalpy.  It would be safe to query the EOS for h_k, but there would possibly be uncertainty in querying h_mixture.  Is any of what I'm saying relevant to what is being discussed? lol

Chris Neal

unread,
Mar 22, 2018, 5:54:08 PM3/22/18
to Cantera Users' Group
Hey Steven,

The perfect gas assumption appears to sneak in with expressing the second term using dh=cpdT. The expansion is just a product rule

Equation:   \rho \frac{Dh_s}{Dt} = \rho \frac{D}{Dt}\left (\sum_{k=1}^{K}h_{s,k}Y_k  \right )   =   \rho \sum_{k=1}^{K}\left (h_{s,k}  \frac{DY_k}{Dt} + Y_k  \frac{Dh_{s,k}}{Dt}  \right ) 

So the book makes the assumption that dh =cpdT and thus the second term becomes a temperature derivative. At least that is my current understanding of the derivation.

Ray Speth

unread,
Mar 22, 2018, 5:56:04 PM3/22/18
to Cantera Users' Group

Steven,

I think that equation does make a perfect gas assumption, but the problem is in the second term, not the first. If we have a definition for h=\sum h_k Y_k, which I think we’ve settled is a requirement for the equation of state, then using the product rule to write

\frac{Dh}{Dt} = \frac{D}{Dt}\left( \sum h_k Y_k \right) = \sum \left( h_k\frac{DY_k}{Dt} + Y_k\frac{Dh_k}{dt} \right)

is perfectly fine, no matter what h_k depends on. However, in order to make the temperature appear, the route that is taken is use the chain rule to write

\frac{Dh_k}{Dt} = \frac{\partial h_k}{\partial T} \frac{DT}{Dt} = c_{p,k} \frac{DT}{Dt}

But this assumes that h_k is a function of temperature only. If h_k is also a function of Y_i, then there should be a bunch of additional terms involving those partial derivatives of the enthalpy as well. 

Regards,
Ray

To unsubscribe from this group and stop receiving emails from it, send an email to cantera-users+unsubscribe@googlegroups.com.

Steven DeCaluwe

unread,
Mar 22, 2018, 6:22:58 PM3/22/18
to canter...@googlegroups.com
Thanks, Ray - I wasn’t thinking that second term all the way through, in terms of where that came from. 

Side note: how are we defining “perfect gas?”  My textbook def says  “ideal gas” + “constant specific heat.” 

Of course, if we’re trying to be super general, this representation completely falls apart when phase change is involved. 

Steven

——————————————————————————————
Steven DeCaluwe, Ph.D.
Assistant Professor of Mechanical Engineering
Colorado School of Mines
G. Brown Hall, W410 B
Golden, CO, 80401

email: deca...@mines.edu

Ray Speth

unread,
Mar 22, 2018, 6:57:05 PM3/22/18
to Cantera Users' Group
Steven,

Sorry, you're right, perfect gas is not what I meant, as letting cp be a function of T is of course fine here. In this instance, the issue is the existence of any dependence of h_k on Y_i (for any i and k). So I suppose you could potentially have some form of non-ideality that doesn't involve such terms.

And Chris, I think this issue is an argument for keeping enthalpy as the variable rather than converting to temperature. The boundary conditions shouldn't be too much of a problem, at least in most cases. If you have boundary conditions for temperature and mass fractions, you should be able to combine those to get the boundary condition for enthalpy. The only case where I think it turns into a mess is if you have surface reactions, where you no longer get to have zero-gradient conditions for the species. I do worry a bit about getting a good Jacobian matrix in this case, in terms of calculating dh/dYi and how you decide to handle the constraint on the sum of mass fractions.

Regards,
Ray

Nick Curtis

unread,
Mar 23, 2018, 9:44:08 AM3/23/18
to Cantera Users' Group
Ray,

Re:

how you decide to handle the constraint on the sum of mass fractions.

What I've done in the past is to drop the mass fraction of the last species from the set of equations, then you have:

\Y_{ns} = 1 - \sum_{i=1}^{ns - 1} Y_{i}

This can complicate derivatives however, as now:

\frac{\partial \Y_{ns}}{\partial Y_{j}} \ne 0 \;, \forall j = 1 \dots ns - 1

Nick Curtis

unread,
Mar 23, 2018, 9:48:45 AM3/23/18
to Cantera Users' Group
Errr, those last two equations were:



and 


Steven DeCaluwe

unread,
Mar 23, 2018, 10:57:00 AM3/23/18
to <cantera-users@googlegroups.com>
I agree, Nick - this seems problematic to me, particularly for any problem involving diffusion.  Minimizing the resulting error requires some assumptions about the mechanism file (iirc, namely that the species with the highest concentration is listed as the final species.  That way the relative error introduced is small).

My thinking lately is that the ‘no norm’ option should be used on the mass fractions, but the implementation of no norm currently has issues that I’ve been meaning to correct.

In the long run, and for my own in-house codes, I prefer to store rho*Y_k for all species as my dependent variables (you don’t need to store rho, and so the # of variables is the same as if you store rho (or P) and Y_k for k=1 to (n-1)).  But this represents a significant re-working for any pre-existing codes.

Steven



——————————————————————————————————
Steven DeCaluwe, PhD
Assistant Professor of Mechanical Engineering
Colorado School of Mines
Brown Building W410B
Golden, CO 80401

Twitter: @DeCaluweGroup

Nick Curtis

unread,
Mar 23, 2018, 2:21:19 PM3/23/18
to Cantera Users' Group
Steven,

One thing you can do (and relatively simply with the python interface) is to simply re-order the mechanism such that the bath gas is placed at the end of the mechanism, e.g.:

# find the last species
n2_ind
= gas.species_index('N2')
# update the gas
specs
= gas.species()[:]
gas
= ct.Solution(thermo='IdealGas', kinetics='GasKinetics',
                  species
=[specs[x] for x in range(gas.n_species) if x != n2_ind] + [specs[n2_ind]],
                  reactions
=gas.reactions())

It can absolutely introduce a lot of additional error if you don't fiddle with the species ordering though (in pyJac we explicitly let the user specify which species to place at the end, defaulting to N2).
It does have the nice property of explicitly conserving mass / satisfying the mass fraction constraint however, but as you note requires some significant reworking of existing code

Nick

Gandhali Kogekar

unread,
Mar 23, 2018, 5:15:52 PM3/23/18
to Cantera Users' Group
Hi Chris,

I ran few quick tests to check whether h_{mix} = \sum {h_k x_k}, and indeed these two terms are different in case of non ideal EoS. The residual enthalpy calculation contains non linear terms and can not be directly summed up to obtain mixture enthalpy.

Steven DeCaluwe

unread,
Mar 23, 2018, 5:50:59 PM3/23/18
to <cantera-users@googlegroups.com>
It looks like the h_k calculations are the partial molar enthalpies, which is of course different from the molar-averaged species enthalpy.  I’ll need to think about how to calculate the h_k that we *do* want…


——————————————————————————————————
Steven DeCaluwe, PhD
Assistant Professor of Mechanical Engineering
Colorado School of Mines
Brown Building W410B
Golden, CO 80401

Twitter: @DeCaluweGroup

Chris Neal

unread,
Mar 26, 2018, 12:18:36 PM3/26/18
to Cantera Users' Group

Thanks Gandhali. I'm trying to tie all of these conversations back down to the energy equation. What your results suggest is that we should avoid performing an analytical substitution in the energy equation of the form:

Equation:  h_{mixture} =   \sum_{k = 1}^{K} x_k h_k 


Because doing so would neglect the behavior of the EOS object(which includes the mixing enthalpy). When you say 'residual' enthalpy, is that the mixing enthalpy? Am I interpreting your results correctly?

(To everyone) I know that all models are incorrect, so for this derivation I suppose the goal is to reduce the ideal gas assumption error as much as possible in the energy equation. Do we have any consensus about how to analytically proceed from the following equation:

Equation:  \rho \frac{Dh}{Dt} = \frac{Dp}{Dt} + \nabla \cdot (\lambda \nabla T) - \sum_{k=1}^{K}\nabla \cdot h_k \vec{j}_k + \Phi



Ray, the form of the enthalpy diffusion is presented in equation 3.155 in Kee's book. There isn't a discussion about where this form originated from, or any references for how it was determined. I believe this is what started us down the rabbit-hole of whether the Kee equation already had an ideal gas assumption buried further up the derivation chain.  I'm trying to figure out what would be the course of action here. 

Equation:  \left ( \frac{dQ}{dt} \right )_{species} = - \sum_{k=1}^{K} \int_{cs} h_k \vec{j}_k \cdot \vec{n} dA 







Are we trying to ensure that the derivation's handling of the enthalpy matches what the Cantera EOS objects provide?

Steven DeCaluwe

unread,
Mar 26, 2018, 7:32:05 PM3/26/18
to canter...@googlegroups.com
Hi Chris,

I would disagree with your first conclusion.  The trick is just getting the *correct* h_k from Cantera.

1.  The equation from Kee’s book is a rather uncontroversial finite-volume approach to the energy equation.  The enthalpy brought into the volume by the mass flux of a species k (j_k) is equal to j_k*h_k, where h_k is the average enthalpy per unit mass of that species.

a. For any situation where you have convection only, you can replace sum(j_k*h_k) with V*rho*h_mix (where rho is the phase’s mass density, kg/m3).

b. But if you have diffusive fluxes of any species, then the assumption upon which this solution relies (j_k = rho*V*Y_k) does not hold.  

2.  So as long as we stop before we convert enthalpy to a function of temperature, the approach in Kee’s book is generally applicable and free from assumptions, and the equation:
h_{mixture} =   \sum_{k = 1}^{K} y_k h_k 
is also free from assumption. It just assumes a specific definition of h_k (see item 4 for what this h_k is; note that I have changed to a per-unit-mass basis, relative to your equation).  The trick is getting Cantera to provide us with the required h_k.  It currently does not, but that is not a knock on the equation, it just means that the current cantera functionality is not that which we need, for this particular application.

3. What Cantera currently calculates is the partial molar enthalpy.  This is essentially the slope of total mixture enthalpy vs. moles of k:

h_k_partial_molar = dH/dn_k)const T,P.

It varies with the compositions, and says “For the current composition, what is the change in enthalpy if I add *one more mole* of species k?”  In other words, the partial molar enthalpy of the first mole of k added (when X_k = 0) will be different from that which makes X_k = 0.5, etc.

4. What we want is the average enthalpy per unit mass of species k.  I.e. assuming we could count the total enthalpy stored in species k in our finite volume, and divide it by the total mass of k in the volume.  This most certainly accounts for any enthalpy of mixing.

The partial molar enthalpy is, in theory, related to the average species enthalpy:

h_k = integral(h_k_partial_molar dn_k) / n_k_tot

The thing I don’t know is whether this relationship is at all helpful to us, or in general whether the h_k that we want is easily derivable from our EoS.  From Gandhali’s previous email we can at least assume that it will not be simple, due to non-linear terms in the EoS.

5.  But regardless, as I stated before, the mixture enthalpy will not be sufficient for any application where there are diffusive fluxes. 

6.  The next step is to look at the EoS and determine if we can calculate the required h_k.


Thanks,
Steven



——————————————————————————————
Steven DeCaluwe, Ph.D.
Assistant Professor of Mechanical Engineering
Colorado School of Mines
G. Brown Hall, W410 B
Golden, CO, 80401

email: deca...@mines.edu

Ray Speth

unread,
Mar 26, 2018, 11:49:23 PM3/26/18
to Cantera Users' Group

Hi all,

At least in the case of the Redlich-Kwong equation of state, I think Cantera does already calculate the needed values for \hat{h}_k such that \sum \hat{h}_k X_k = \hat{h} = h * \bar{W}, and these is what is returned as the “partial molar enthalpies”. For example, using the test input file for the R-K EOS, we can calculate (a) the weighted sum of the partial molar enthalpies, (b) the molar enthalpy, and (c), the weighted sum of the pure species enthalpies. At low pressures, where the non-ideal effects are small, these all give essentially the same values:

def h_pure(g):
    h = 0.0
    save = g.state
    for species,Xk in g.mole_fraction_dict().items():
        g.TPX = g.T, g.P, {species: 1.0}
        h += g.enthalpy_mole * Xk
    g.state = save
    return h

g.TPX = 500, 1e5, 'CO2:0.3, H2O:0.6, H2:0.1'
print('(a1)', sum(g.partial_molar_enthalpies * g.X))
print('(b1)', g.enthalpy_mole)
print('(c1)', h_pure(g))
(a1) -255933851.418
(b1) -255933851.41823298
(c1) -255938514.87131906

Closer to the critical point, the first two calculations agree, while the third diverges, as expected:

g.TPX = 480, 1e7, 'CO2:0.3, H2O:0.6, H2:0.1'
print('(a2)', sum(g.partial_molar_enthalpies * g.X))
print('(b2)', g.enthalpy_mole)
print('(c2)', h_pure(g))
(a2) -259392417.488
(b2) -259392417.48808572
(c2) -273891076.90921915

Gandhali, I’d be interested in seeing the tests that you did where you found that this summation does not work.

Steven, do you think the above is consistent with the definition of h_k_partial_molar that you mentioned above? How does the above definition of h_k avoid path dependencies?

Also, Steven, I thought j_k were usually defined to be only the diffusive mass fluxes, such that sum(j_k) = 0, which would fit with the net convective flux of enthalpy being captured in the total derivative on the LHS of the equation. Whether it makes sense to necessarily define a diffusion velocity or just directly relate the diffusive mass flux to the relevant gradients is I think independent of the questions about the enthalpy.

Regards,
Ray

On Monday, March 26, 2018 at 7:32:05 PM UTC-4, S. DeCaluwe wrote:

Hi Chris,

Steven DeCaluwe

unread,
Mar 27, 2018, 1:13:08 PM3/27/18
to <cantera-users@googlegroups.com>
Hi Ray,

1.  Interesting.  I am guessing Gandhali’s tests were in Matlab, so there may be something there to follow up on. 

2. Yeah, looking more carefully, the definition of partial molar enthalpy *is* what we want:
I think I was just over-reaching/over-thinking to try and come up with an explanation as to why Gandhali’s test failed.  

3.  As for j_k being the diffusive fluxes, I am sure you’re right.  That email yesterday was on my phone from an airport, so I really don’t have all the equations at my disposal, and even keeping track all the various equations in this thread is beyond me, at this point.  But in Chris’s first email on this, I think, there was some comparison between the sum(h_k*j_k) in Kee’s equation and something along the lines of (V_k*Y_k*rho) in another presentation of the governing equations.

I just wanted to make it clear that the reason the latter is insufficient is because of diffusive fluxes.  I agree that the representation of the diffusive fluxes is independent of the question of the enthalpies.


4.  Regarding diffusive fluxes vs. the relative gradients, this would be the next can of worms.  
    a.  For generalizing beyond ideal gas, the fluxes would be as a function of the chemical potential gradient.
    b.  It also assumes that you can lookup diffusion coefficients which are generally applicable.

The first item of course affects the code structure.  Writing the code in a generalized manner shouldn’t be too tough, but I do wonder how much/whether this slows down simulations.

The second item is not really relevant to the code structuring: it is not the code’s concern whether somebody’s diffusion coefficients are correct, and the software’s job is not to prevent people from making any and all mistakes, but rather “given the inputs you supplied, here are the results.”  
The only reason I bring it up is that a few years ago I put in some capabilities for high-pressure gas diffusion coefficients, using some curve-fitting and the theory of corresponding states.  But this method, if I remember correctly, is range limited - the correlations are only valid over a certain range of conditions, and I think the model returns an error outside this range.  I’ll have to check to see if I’m remembering correctly.  
Anyway, just mentioning it because it might require some thought as to how the simulation would handle a case where that exception is thrown.  Or modifications to the ‘highPressureGas’ transport model.

Thanks,
Steven

——————————————————————————————————
Steven DeCaluwe, PhD
Assistant Professor of Mechanical Engineering
Colorado School of Mines
Brown Building W410B
Golden, CO 80401

Twitter: @DeCaluweGroup

Gandhali Kogekar

unread,
Mar 27, 2018, 1:13:34 PM3/27/18
to Cantera Users' Group
Hi Ray,

I performed exactly same calculations, as you did. Basically, I compared mean enthalpy of the mixture with \sum h_k x_k. However, I was using Matlab instead of Python. And now it seems that Matlab subroutine for partial molar enthalpies always returns ideal gas values and is not implemented for non-ideal EoS. This explains the inconsistency between our results.

Steven DeCaluwe

unread,
Mar 27, 2018, 1:16:54 PM3/27/18
to <cantera-users@googlegroups.com>
Ha - I think our emails crossed paths in the ethernet cables.  

Interesting/disconcerting.  This is a good catch, and a significant error.  We’ll get this fixed, in short order.

The Matlab toolbox: the gift that keeps on giving! :)



——————————————————————————————————
Steven DeCaluwe, PhD
Assistant Professor of Mechanical Engineering
Colorado School of Mines
Brown Building W410B
Golden, CO 80401

Twitter: @DeCaluweGroup

Chris Neal

unread,
Mar 29, 2018, 5:54:12 PM3/29/18
to Cantera Users' Group

Going to take a shot at deriving the energy equation in temperature form without assuming ideal gas relationships for enthalpy.

Starting with Kee's general energy equation.

Equation:  \rho \frac{Dh}{Dt} = \frac{Dp}{Dt} + \nabla \cdot (\lambda \nabla T) - \sum_{k=1}^{K}\nabla \cdot h_k \vec{j}_k + \Phi


Substituting the following definition into the total derivative for enthalpy, 

Equation:  h =  \sum_{k=1}^{K} Y_k h_k 


using the product rule to expand the enthalpy total derivative gives,

Equation:  \rho \frac{Dh}{Dt} = \rho \frac{D}{Dt}\left (\sum_{k=1}^{K}h_{k}Y_k  \right )   =   \rho \sum_{k=1}^{K}\left (h_{k}  \frac{DY_k}{Dt} + Y_k  \frac{Dh_{k}}{Dt}  \right ) 


The total derivative of the enthalpy of species k can be expanded using the chain rule under the assumption that for an individual species, h_k(T,P).

Equation:  \frac{Dh_k}{Dt} = \frac{Dh_k}{DT}  \frac{DT}{Dt}  + \frac{Dh_k}{DP}  \frac{DP}{Dt} 


The above is essentially eq.3.217, and using equation 3.221, we can write the above relation as follows.

Equation:  \frac{Dh_k}{Dt} =  \frac{D}{Dt} \left [ c_{p,k} DT + \frac{1}{\rho_k}(1-T\beta_k)DP  \right ]   =  c_{p,k} \frac{DT}{Dt} + \frac{1}{\rho_k}(1-T\beta_k) \frac{DP}{Dt}  


If I take the above relation and substitute it into that product rule expansion for the total enthalpy derivative, the expression becomes,

Equation:    \rho \sum_{k=1}^{K}\left (h_{k}  \frac{DY_k}{Dt} + Y_k  c_{p,k} \frac{DT}{Dt} + \frac{Y_k}{\rho_k}(1-T\beta_k) \frac{DP}{Dt}   \right ) 


If I pull out total derivatives of P and T from the summation and distribute the summation across the terms, the following expression is obtained.

Equation:   \rho \sum_{k=1}^{K} h_{k}  \frac{DY_k}{Dt} + \rho \frac{DT}{Dt} \sum_{k=1}^{K} Y_k  c_{p,k}  + \rho \frac{DP}{Dt}  \sum_{k=1}^{K} \frac{Y_k}{\rho_k}    -  \rho \frac{DP}{Dt} T \sum_{k=1}^{K} Y_k \frac{\beta_k}{\rho_k} 


Substituting 3.124(species continuity) in for the total derivative of the species mass fraction, distributing the summation over the resulting expression, and re-introducing the right-hand-side of the equation gives,

Equation:   \sum_{k=1}^{K} h_{k}  \dot{\omega}W_k -  \sum_{k=1}^{K} h_{k} \nabla \cdot \vec{j}_k  + \rho \frac{DT}{Dt} \sum_{k=1}^{K} Y_k  c_{p,k}  + \rho \frac{DP}{Dt}  \sum_{k=1}^{K} \frac{Y_k}{\rho_k}    -  \rho \frac{DP}{Dt} T \sum_{k=1}^{K} Y_k \frac{\beta_k}{\rho_k}  =  \frac{Dp}{Dt} + \nabla \cdot (\lambda \nabla T) - \sum_{k=1}^{K}\nabla \cdot h_k \vec{j}_k + \Phi


Use the following identity to expand the enthalpy flux on the right-hand-side of the above equation,

Equation:   \nabla \cdot h_k \vec{j}_k = h_{k} \nabla \cdot \vec{j}_k  +  \vec{j}_k \cdot \nabla h_{k} 


Performing the substitution and cancelling the common terms on the right and left sides gives,


Equation:   \left (  \rho \sum_{k=1}^{K} Y_k  c_{p,k} \right ) \frac{DT}{Dt} + \left (\rho \sum_{k=1}^{K} \frac{Y_k}{\rho_k}  \right )  \frac{DP}{Dt}   -  \left ( \rho T \sum_{k=1}^{K} Y_k \frac{\beta_k}{\rho_k}  \right ) \frac{DP}{Dt} =  \\ \frac{Dp}{Dt} + \nabla \cdot (\lambda \nabla T)  - \sum_{k=1}^{K} \vec{j}_k \cdot \nabla h_{k}   + \sum_{k=1}^{K} h_{k} \dot{\omega}W_k + \Phi


At this point we have several summation terms that are individual species properties weighted by the species mass fractions. Whether to leave these as summations or to make the jump back to the mixture property is where I'm currently stuck. I'm leaning towards it being an ok thing to do, because we started the derivation by using the same assumption. Then there is also the gradient of the individual species enthalpies, when I'm unsure about what can be done with the term. It resembles the original temperature energy equation 3.213.


On Tuesday, March 27, 2018 at 1:16:54 PM UTC-4, S. DeCaluwe wrote:
Ha - I think our emails crossed paths in the ethernet cables.  

Interesting/disconcerting.  This is a good catch, and a significant error.  We’ll get this fixed, in short order.

The Matlab toolbox: the gift that keeps on giving! :)


——————————————————————————————————
Steven DeCaluwe, PhD
Assistant Professor of Mechanical Engineering
Colorado School of Mines
Brown Building W410B
Golden, CO 80401

Steven DeCaluwe

unread,
Mar 29, 2018, 6:40:13 PM3/29/18
to <cantera-users@googlegroups.com>
Hi Chris,

Thanks - I think this comprehensive approach is a much more useful way to go about this conversation - have all equations in one monster thread so we can keep track of it all!

So my current thoughts:
1. It would be helpful if we could also have our own numbering for these equations.  I’ve added some, below, in red.

2. I currently only get to equations 4 and 5 before I hit issues (partially issues with the equations, partially limitations in my own knowledge/understanding):
a. For equation 4, h_k is also a function of the composition.  If I understand correctly, that means there should be a third term on this sum, correct?  

b. Can you either include eq. 3.221 or explain the variable beta_k for me?  I don’t know what it represents.

c. I feel like the substitution Dh_k/DT = Cp_k makes an ideal gas assumption:
Strictly speaking, c_p = \frac{\partial{h}}{\partial{T}}\right)_{const. P}.  For an ideal gas, h is a function of T only, the ‘constant P’ condition becomes
irrelevant, and  c_p = Dh/DT.  For our purposes, h_k is going to be a function of P, T, and composition, so we cannot make this substitution.

d. Moreover, using c_p is wholly unsuitable if we will ever hit saturation / phase change.  Since h changes discontinuously during phase change, c_p is undefined, here.  That said, I’m not sure if this is something we want/need to account for, or whether we just stipulate that the model assumes a single phase throughout 
the simulation, and slap a big “CAVEAT EMPTOR” on the whole thing...

e. It seems like there is something not equivalent in going from eq. 4 to eq. 5: should the preceding term in the middle equation be 1/Dt, rather than D/Dt?

3. General thought for the overall derivation: many of the real-gas equations of state will be writing in terms of molar volume v and mole fractions X_k, rather than P and Y_k.  I’m not sure if this will complicate calculation of the necessary terms in the equations below (which are in terms of P and Y_k), but it is something on my mind.

Thanks for all your work, here - glad to be moving this derivation forward.

Steven

——————————————————————————————————
Steven DeCaluwe, PhD
Assistant Professor of Mechanical Engineering
Colorado School of Mines
Brown Building W410B
Golden, CO 80401

Twitter: @DeCaluweGroup
On Mar 29, 2018, at 3:54 PM, Chris Neal <chri...@snumerics.com> wrote:

Going to take a shot at deriving the energy equation in temperature form without assuming ideal gas relationships for enthalpy.

Starting with Kee's general energy equation.

Eq. 1

Equation:  \rho \frac{Dh}{Dt} = \frac{Dp}{Dt} + \nabla \cdot (\lambda \nabla T) - \sum_{k=1}^{K}\nabla \cdot h_k \vec{j}_k + \Phi


Substituting the following definition into the total derivative for enthalpy, 


Eq. 2

Equation:  h =  \sum_{k=1}^{K} Y_k h_k 


using the product rule to expand the enthalpy total derivative gives,

Eq. 3

Equation:  \rho \frac{Dh}{Dt} = \rho \frac{D}{Dt}\left (\sum_{k=1}^{K}h_{k}Y_k  \right )   =   \rho \sum_{k=1}^{K}\left (h_{k}  \frac{DY_k}{Dt} + Y_k  \frac{Dh_{k}}{Dt}  \right ) 


The total derivative of the enthalpy of species k can be expanded using the chain rule under the assumption that for an individual species, h_k(T,P).

Eq. 4 

Equation:  \frac{Dh_k}{Dt} = \frac{Dh_k}{DT}  \frac{DT}{Dt}  + \frac{Dh_k}{DP}  \frac{DP}{Dt} 


The above is essentially eq.3.217, and using equation 3.221, we can write the above relation as follows.

Eq. 5

Equation:  \frac{Dh_k}{Dt} =  \frac{D}{Dt} \left [ c_{p,k} DT + \frac{1}{\rho_k}(1-T\beta_k)DP  \right ]   =  c_{p,k} \frac{DT}{Dt} + \frac{1}{\rho_k}(1-T\beta_k) \frac{DP}{Dt}  


If I take the above relation and substitute it into that product rule expansion for the total enthalpy derivative, the expression becomes,

Eq. 6 

Equation:    \rho \sum_{k=1}^{K}\left (h_{k}  \frac{DY_k}{Dt} + Y_k  c_{p,k} \frac{DT}{Dt} + \frac{Y_k}{\rho_k}(1-T\beta_k) \frac{DP}{Dt}   \right ) 


If I pull out total derivatives of P and T from the summation and distribute the summation across the terms, the following expression is obtained.

Eq. 7 

Equation:   \rho \sum_{k=1}^{K} h_{k}  \frac{DY_k}{Dt} + \rho \frac{DT}{Dt} \sum_{k=1}^{K} Y_k  c_{p,k}  + \rho \frac{DP}{Dt}  \sum_{k=1}^{K} \frac{Y_k}{\rho_k}    -  \rho \frac{DP}{Dt} T \sum_{k=1}^{K} Y_k \frac{\beta_k}{\rho_k} 


Substituting 3.124(species continuity) in for the total derivative of the species mass fraction, distributing the summation over the resulting expression, and re-introducing the right-hand-side of the equation gives,

Eq. 8 

Equation:   \sum_{k=1}^{K} h_{k}  \dot{\omega}W_k -  \sum_{k=1}^{K} h_{k} \nabla \cdot \vec{j}_k  + \rho \frac{DT}{Dt} \sum_{k=1}^{K} Y_k  c_{p,k}  + \rho \frac{DP}{Dt}  \sum_{k=1}^{K} \frac{Y_k}{\rho_k}    -  \rho \frac{DP}{Dt} T \sum_{k=1}^{K} Y_k \frac{\beta_k}{\rho_k}  =  \frac{Dp}{Dt} + \nabla \cdot (\lambda \nabla T) - \sum_{k=1}^{K}\nabla \cdot h_k \vec{j}_k + \Phi


Use the following identity to expand the enthalpy flux on the right-hand-side of the above equation,

Eq. 9 

Equation:   \nabla \cdot h_k \vec{j}_k = h_{k} \nabla \cdot \vec{j}_k  +  \vec{j}_k \cdot \nabla h_{k} 


Performing the substitution and cancelling the common terms on the right and left sides gives,

Eq. 10 

Ray Speth

unread,
Mar 29, 2018, 11:37:31 PM3/29/18
to Cantera Users' Group

Steven,

I agree that the first issue is with the expansion of \frac{Dh_k}{Dt} (Eq. 4), and it’s the same problem that I previously mentioned with the application of the chain rule to the total enthalpy. I believe that the correct expansion, accounting for the fact that h_k = h_k(T, P, Y_1, \ldots, Y_K) (with K species) would be

\frac{Dh_k}{Dt} = \left,\frac{\partial h_k}{\partial T}\right|_{P,Y} \frac{DT}{Dt} + \left,\frac{\partial h_k}{\partial P}\right|_{T,Y} \frac{DP}{Dt} +\sum_{i=1}^K  \left,\frac{\partial h_k}{\partial Y_i}\right|_{T,P,Y}\frac{DY_k}{Dt}

(noting that the total derivative only makes sense for the D/Dt terms — the rest are just regular partial derivatives). I think the term multiplying DT/Dt is in fact c_{p,k}. The second term can be neglected if we make the usual zero Mach number assumption. The final term is the one that’s a mess — I’m pretty sure that Cantera doesn’t currently calculate these partial derivatives. And in fact, you get the fun case of having those partial derivatives depend on your choice of how the \sum Y_k = 1 condition is enforced.

In the case of a phase change, I think you’re going to run into problems either way. If you try to write the equations in terms of temperature, you have the problem of infinite cp, and if you try to keep using enthalpy, you have to allow for it to be discontinuous at the phase boundary. I would imagine that you would have to do explicitly model the phase boundary, in the absence of some really clever tricks.

I’m pretty convinced that keeping enthalpy as the state variable is the way to go.

Regards,
Ray

On Thursday, March 29, 2018 at 6:40:13 PM UTC-4, S. DeCaluwe wrote:

Hi Chris,

Chris Neal

unread,
Apr 2, 2018, 1:53:37 PM4/2/18
to Cantera Users' Group
Thanks for the feedback Steven.

1.) I agree that having a coherent derivation that isn't piecemeal split across many different threads/responses. Along those lines I've moved the derivation onto the ShareLatex online Latex document editor. Anyone with the Link can edit the document.

2.) 
   a.) This is where my understanding is shaky. The consideration of the enthalpy of the mixture as a weighted sum of the constituent enthalpies of the species in the mixture (equation 2 in the ShareLatex document) is why I assumed that when dealing with the h_k quantity there was no consideration of the composition i.e. h_k(T,P), whereas h(T,P,Yk)->Where Yk is introduced as the weighting factor.

  b.) I updated the document for explicitly show the equations 3.217 and 3.221 as equations (5) and (6) in the document.

  c.) Your observation gives evidence to Ray's assertion that sticking with enthalpy might be the way to do in order to avoid all of these unpleasant assumptions in getting to a workable temperature form of the energy equation.

  d.) I agree that this would be a situation that would have to be either handled or have the limitation documented.

  e.) I cleaned up that equation. If you divide both sides of (6) by Dt and compare to (4), then that is how (7) is determined.

3.) I had not realized that. And I'm a little unsure about how to go about "re-phrasing" the equations in terms of these alternative quantities.


On Thursday, March 29, 2018 at 6:40:13 PM UTC-4, S. DeCaluwe wrote:
Hi Chris,

Thanks - I think this comprehensive approach is a much more useful way to go about this conversation - have all equations in one monster thread so we can keep track of it all!

So my current thoughts:
1. It would be helpful if we could also have our own numbering for these equations.  I’ve added some, below, in red.

2. I currently only get to equations 4 and 5 before I hit issues (partially issues with the equations, partially limitations in my own knowledge/understanding):
a. For equation 4, h_k is also a function of the composition.  If I understand correctly, that means there should be a third term on this sum, correct?  

b. Can you either include eq. 3.221 or explain the variable beta_k for me?  I don’t know what it represents.

c. I feel like the substitution Dh_k/DT = Cp_k makes an ideal gas assumption:
Strictly speaking, c_p = \frac{\partial{h}}{\partial{T}}\right)_{const. P}.  For an ideal gas, h is a function of T only, the ‘constant P’ condition becomes
irrelevant, and  c_p = Dh/DT.  For our purposes, h_k is going to be a function of P, T, and composition, so we cannot make this substitution.

d. Moreover, using c_p is wholly unsuitable if we will ever hit saturation / phase change.  Since h changes discontinuously during phase change, c_p is undefined, here.  That said, I’m not sure if this is something we want/need to account for, or whether we just stipulate that the model assumes a single phase throughout 
the simulation, and slap a big “CAVEAT EMPTOR” on the whole thing...

e. It seems like there is something not equivalent in going from eq. 4 to eq. 5: should the preceding term in the middle equation be 1/Dt, rather than D/Dt?

3. General thought for the overall derivation: many of the real-gas equations of state will be writing in terms of molar volume v and mole fractions X_k, rather than P and Y_k.  I’m not sure if this will complicate calculation of the necessary terms in the equations below (which are in terms of P and Y_k), but it is something on my mind.

Thanks for all your work, here - glad to be moving this derivation forward.

Steven
——————————————————————————————————
Steven DeCaluwe, PhD
Assistant Professor of Mechanical Engineering
Colorado School of Mines
Brown Building W410B
Golden, CO 80401


<div style="letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;wo

Chris Neal

unread,
Apr 2, 2018, 3:06:54 PM4/2/18
to Cantera Users' Group
Ray,  Just to reiterate(for myself).  h_k = f(T, P, Yk) and not f(T,P) is because h_k is not the enthalpy of species k, but rather the partial enthalpy(derivative of enthalpy w.r.t. each species) which has an explicit dependence of the mass fraction. Something like this.

Ray Speth

unread,
Apr 2, 2018, 9:52:00 PM4/2/18
to Cantera Users' Group
Chris,

Equation 3.217 from Kee et al. is, as the section title states, only true for a single-component fluid -- the immediately preceding statement is the assumption that h=h(p, T).

Equation 4 is not correct because h_k, which is the enthalpy per unit mass of species k, depends on (potentially) all of the mass fractions in addition to the temperature and pressure. The h_k appearing in Equation 2 are not the pure species properties.

Regards,
Ray
<span style="

Chris Neal

unread,
Apr 3, 2018, 2:08:03 PM4/3/18
to Cantera Users' Group
Thanks Ray. So it looks to me like equation (1) in the Latex writeup is just about as far as we should go with respect to having a governing equation that has no ideal gas assumptions & does not run into infinite cp issues for phase changes? There is the choice of using the enthalpy(sensible+chemical) equation or the sensible enthalpy equation(where the chemical source terms are explicitly included & the enthalpy diffusion heat flux term is also present). I'm adding some equations to the writeup to illustrate the different forms.
<span style="text-align

Ray Speth

unread,
Apr 4, 2018, 11:38:35 AM4/4/18
to Cantera Users' Group
Chris,

Cantera works almost exclusively with enthalpies including both the sensible and chemical terms, so I would expect keeping that convention to make the most sense from a consistency standpoint. This is the approach taken in the ConstPressureReactor class, for example.

Regards,
Ray
<span styl

Chris Neal

unread,
Apr 17, 2018, 5:14:53 PM4/17/18
to Cantera Users' Group
Thanks Ray,

I took a look at the ConstPressureReactor documentation(https://www.cantera.org/docs/sphinx/html/reactors.html#governing-equations-for-single-reactors), and I see the enthalpy form of the energy equation that you mentioned. That enthalpy is referred to as the 'total enthalpy', which I understand to mean the sensible + chemical + kinetic. I believe the energy equation in Kee's book(3.203) that we started with is the sensible + kinetic i.e. not the total enthalpy. Unless there is just a syntactical difference between what is meant in the documentation. If we wanted to use the total enthalpy form, we'd have to get from 3.187 to the equivalent total enthalpy form of that equation.  So for consistency between the 0D and 1D reactors, should the total enthalpy form of the energy equation be used?

Thanks,
Chris
<span style="white-sp

Ray Speth

unread,
Apr 17, 2018, 7:22:54 PM4/17/18
to canter...@googlegroups.com
Chris,

I'm not sure why you think Eq. 3.203 is not an equation for the total enthalpy. The following paragraphs explain that this is why there is no term which explicitly contains the reaction rates.

Regards,
Ray

--

Chris N

unread,
Apr 18, 2018, 12:33:36 PM4/18/18
to Cantera Users' Group
I'm probably confused between my different references. It's not the source terms though that are causing me an issue with understanding. In Poinsot's book they derive a total enthalpy equation:


Equation:  \rho \frac{Dh_t}{Dt} = \frac{\partial p}{\partial t} - \frac{\partial q_i}{\partial x_i} +  \frac{\partial}{\partial x_j} (\tau_{ij} u_i)

Once the kinetic energy is removed from the total enthalpy equation the following form is given by Poinsot.


Equation:  \rho \frac{Dh}{Dt} = \frac{Dp}{Dt} - \frac{\partial q_i}{\partial x_i} + \tau_{ij} \frac{\partial u_i}{\partial x_j} 


The second form looks like 3.203 to me, so I assumed that this wasn't total enthalpy in the sense that total enthalpy = sensible + chemical + kinetic. Am I glossing over some sort of "poh-tay-toe", "poh-tah-toe" situation here where they are equivalent?

Ray Speth

unread,
Apr 24, 2018, 9:51:59 PM4/24/18
to Cantera Users' Group
Chris,

I think the confusion here is over the introduction of the kinetic energy, which we were not previously discussing. I am using the term 'total enthalpy' to refer to chemical + sensible enthalpy, although I do recognize that this term is sometimes used to refer to sensible + kinetic enthalpy. The formulation for the 1D flame used in Cantera makes use of the low Mach number approximation, whereby the thermodynamic pressure can be assumed to be constant through the domain, and I think it would make sense to continue making that approximation while generalizing to a real gas equation of state. Similarly, I would almost immediately discard the term which represents the viscous heat generation, which is obviously going to be much smaller than the chemical heat release rate.

Regards,
Ray
<blockquote class="gmail_quote" style="margin:0;margin-left:0.8ex;bor

Chris N

unread,
Apr 25, 2018, 4:06:49 PM4/25/18
to Cantera Users' Group
That clears things up for me Ray.  So equation 3.203 it is then.  I suppose now we can move on to discussion/implementation of the energy equation into the 1D solver. Ground zero for the energy equation change looks to be the StFlow class. I haven't yet found any other classes that would be affected by the change yet.
<blockquote class="gmail_quote" style

CatDog

unread,
May 10, 2018, 11:54:04 PM5/10/18
to Cantera Users' Group
Hi everyone,

I would suggest you implement a real gas compressible flow solver first if you really need to consider real gas effects in flow seriously. Because for low mach number flow, the absolute pressure is almost the same and we can neglect its effects in EoS (in most cases, in certain condition small variation of pressure is critical), so we got a "ideal gas"-like behavior: X=X(T,N_i;P=P0) = sum(N_i*X_i(T; P=P0)) + X_real. We can The non-ideal gas effects is in the high order term. However, UV explosion is an exception due to large absolute pressure variation.

Another problem is the fugacity, the pressure of real gas does not equal fugacity. When we need to calculate chemical equilibirum, fugacity is needed to replace pressure. But fugacity is the partial molar pressure. Partial molar quantity is difficult to calculate.

About thermodynamic derivatives, Bridgman's table would be very helpful. According to the table, not all thermodynamic derivatives are independent.

About the partial molar quantities, they are very useful in formula derivation, but they are very difficult to understand and calculate. According to Euler's integral , energy and entropy is extensive quantity, so they can be represented as sum of the corresponding partial molar quantities times mass or mole.

About the energy equation, it is more convenient to use total enthalpy (= sensible + chemical + mixing) and mass fraction in compressible CFD solver, which is usually density based. But we can still use Temperature as solution variable in FVM solver. (ref: Fluent Theory Guide)

I also have some suggestion about the implementation of non-ideal gas EoS in Cantera. According to some literature I read about Refprop of NIST and other real gas, they think the future of EoS implementation is Helmholtz form A=A(T,v,X_i). Because of thermodynamic completeness, if we have a explicit expression of Helmholtz energy, we can get Gibbs, Internal Energy, Enthalpy, Entropy, Chemical potential, etc. in explicit form through derivatives of Helmholtz energy. And Helmholtz energy is also very useful in phase equilibrium calculation (aka. flash, dew, bubbling point). Helmholtz energy is also very useful in CFD because the specific volume v in its formula is related to inverse of density, which is usually used as solution variable in compressible solver. It is more convenient than pressure in compressible solver. Further, all cubic form EoS (Peng-Robinson, VdW, Redlich-Kwong, etc. ) can be converted to Helmholtz form directly. So we can use auto-differentiation to implement all other thermodynamic quantities with accuracy and easiness.

Chris Neal

unread,
Jun 4, 2018, 1:53:25 PM6/4/18
to Cantera Users' Group
Thank you for your suggestions CatDog. Some of the suggestions go over my head(like where the fugacity issue would arise). The Bridgman table is very helpful.  The conclusion of the discussion on the thread is in agreement with your suggestion to use the total enthalpy equation for the form of the energy equation.

I'm ignorant about the Helmholtz EoS topic. I know that there is some work that could be done to have Cantera interface to CoolProp(open source REFPROP), and I believe that CoolProp probably does things close to how REFPROP does. So maybe that issue could be handled through the use of the CoolProp interface(I don't know though, just typing my thoughts).

Christopher Neal

unread,
Jun 6, 2018, 6:19:06 PM6/6/18
to Cantera Users' Group
Ray,

I've updated that PDF with the notes on the energy equation. I'm currently working to document the discretization for the enthalpy form of the energy equation. One issue that I have run into is how to handle the temperature when it shows up in the enthalpy energy equation in the form of the diffusion term. One thought that I had would be to re-define the "T" method that used to return the solution vector with one that queries the solution vector for the enthalpy at that point, and then uses an EOS call to determine the temperature at that point. I've attached a pdf of my notes, and anyone can edit the file here if they want to: https://www.sharelatex.com/3151828574zdbhmtwjbrxz

Thanks,
Chris
<blockquote class="gmail_quote" style
Cantera_Energy_Equation.pdf

raghavki...@gmail.com

unread,
Oct 20, 2020, 3:20:49 PM10/20/20
to Cantera Users' Group
Apologize for reviving this thread, but does Cantera have the capability of real gas 1D flame models? If not, what is a good alternative to do such a simulation? Chemkin?

Thanks
Raghav Sharma
Reply all
Reply to author
Forward
0 new messages