[sheaker]

I am just an amateur openFOAM user but I always use Cantera for chemical calculations. So I recommend You to calculate Your own coefficients for Your case. Firstly - You will get some experience for further research, secondly - You will be able to compare Your output to two previous.

In case You wont calculate new coefficients I think that Your CP values are more realistic in lower range of temperature but immortality values seems more realistic in higher range of temperature. (Do not mix them!)

Have a nice day!
Sheaker

[er99]

Originally Posted by sheaker

I am just an amateur openFOAM user but I always use Cantera for chemical calculations. So I recommend You to calculate Your own coefficients for Your case. Firstly - You will get some experience for further research, secondly - You will be able to compare Your output to two previous.

In case You wont calculate new coefficients I think that Your CP values are more realistic in lower range of temperature but immortality values seems more realistic in higher range of temperature. (Do not mix them!)

Have a nice day!
Sheaker

Hi,

Thanks for your reply. I did calculate my own values, but I cited values from paper (see my previous post) and computed the high and low values from the air composition percentages. I did not use Cantera.

My simulations won't go above 350 K - so you said at lower temperatures mine look more realistic. I was wondering what made you make this assertion? Why do you think mine is more realistic?

Thanks

[sheaker]

Hello.
At first I though Your CP(T) was more realistic because I thought that there is nothing special in air properties in range T=223.16K and T=323.16K so heat capacity should be monotonic as Yours. Also I remember some (not for air) diagrams of CP(T) from old "Citrination" site where CP -> 0 for temperature -> 0; that's why I had typed Your CP as more realistic.
But according to Engineering Toolbox I was wrong.
Have a nice day.
Oskar

[er99]

Originally Posted by sheaker

Hello.
At first I though Your CP(T) was more realistic because I thought that there is nothing special in air properties in range T=223.16K and T=323.16K so heat capacity should be monotonic as Yours. Also I remember some (not for air) diagrams of CP(T) from old "Citrination" site where CP -> 0 for temperature -> 0; that's why I had typed Your CP as more realistic.
But according to Engineering Toolbox I was wrong.
Have a nice day.
Oskar

Could you finish your thought? What does that mean? Mine is wrong?

[sheaker]

I don't know. It is hard to say. "Immortality" output is similar to the one from Engineering Toolbox and Your output is similar to results from Cantera. (see attachment)

It is up to You to decide which one is proper.

I wish You best!
Oskar

[er99]

Originally Posted by sheaker

I don't know. It is hard to say. "Immortality" output is similar to the one from Engineering Toolbox and Your output is similar to results from Cantera. (see attachment)

It is up to You to decide which one is proper.

I wish You best!
Oskar

So from my references the O2 and N2 values are valid from 200 to 6000K.

References:
https://ntrs.nasa.gov/archive/nasa/c...9940013151.pdf
http://garfield.chem.elte.hu/Burcat/THERM.DAT

[er99]

Originally Posted by sheaker

I don't know. It is hard to say. "Immortality" output is similar to the one from Engineering Toolbox and Your output is similar to results from Cantera. (see attachment)

It is up to You to decide which one is proper.

I wish You best!
Oskar

Hi,

Have you plotted Cp(T) for my JANAF values? I'm getting something odd (see attached). For room temperature ~290K I seem to get a CP of 985 J/K and it should be 1005. Did you get that too using my values?

Thanks

[sheaker]

Dear er10

My plot from post #45 comes from Cantera 2.0.

I am sorry but I haven't checked Your values because they are divided by gas constant (and I was to lazy to check them :/ ). But I have recommended to You to do Your own calculations with Cantera to compare them with previous plots!

[er99]

Originally Posted by sheaker

Dear er10

My plot from post #45 comes from Cantera 2.0.

I am sorry but I haven't checked Your values because they are divided by gas constant (and I was to lazy to check them :/ ). But I have recommended to You to do Your own calculations with Cantera to compare them with previous plots!

Hi Sheaker!

I didn't do the correct calculations! My values actually work See attached! The "literature" markers are from the "Fundamentals of Heat and Mass Transfer" book, specifically Appendix A.4. I'm pretty sure my values are correct for the full range and their tables don't account for high temperature effects and that's where the difference comes from above a certain T. I'm only interested in 290 < T < 350 so I'm happy with my values.

Thanks for your kind and prompt replies!

[sheaker]

Now it looks good!

I wish You quick progress in Your simulation!

Have a nice day.
Sheaker

[arvindpj]

Originally Posted by sheaker

Hello. I think I have some experience with OLD! JANAF and openFOAM.

To calculate all a0 - a6 and b0 - b6 coefficient we need to use some additional tools like Cantera.

With Cantera installed we can use cpPolynomials.py:
http://chomikuj.pl/Sheaker/Studia/Me...,5632310992.py (I hope You can download it for free without registration) to calculate coefficient of Cp(T) of defined gas mixture for specific temperature ranges. Remember that Cantera outputs Cp(T) = [J/molK] while openFoam needs Cp(T)/R = [-].

To calculate b5 and b6 coefficient make sure Your H(T=273,15K) is equal to Standard Enthalpy of Formation and Your S(T=273,15K) is equal to Standard Entropy. Values of Standard Enthalpy of Formation and Standard Entropy comes from Cantera basic output.

To calculate a5 and a6 just make Your functions tangent near T_common.

I hope it put some light on JANAF coefficient.

Could you explain a little bit more? I am trying to simulate the buoyant flow of flue gas mimicked as a mixture [ CO2 - 13%; H20 -11%; N2 -76%].

I was using these fixed values.

 

t	ρ	cp	μ *106	ν *106
[O C]	[kg/m3]	[kJ/kgK] [Pas]	[m2/s]
0	1.295	1.042	15.8	12.2
100	0.95	1.068	20.4	21.54
200	0.748	1.097	24.5	32.8
300	0.617	1.122	28.2	45.81
400	0.525	1.151	31.7	60.38
500	0.457	1.185	34.8	76.3
600	0.405	1.214	37.9	93.61
700	0.363	1.239	40.7	112.1
800	0.33	1.264	43.4	131.8
900	0.301	1.29	45.9	152.5
1000	0.275	1.306	48.4	174.3
1100	0.257	1.323	50.7	197.1
1200	0.24	1.34	53	221

Now, I want to simulate using temperature varying properties. As a first step, I just migrated your code for cantera 2.2.1.

from cantera import *
import sys
import numpy as np
from numpy import *
import math
import matplotlib.pyplot as plt


def cpPolynomials(P1,q,mech,tempRange1,tempRange2,step):

   # gas = importPhase(mech);
    gas = Solution(mech);
    # specific heat of fresh mixture
    
    a = tempRange1[0]
    b = tempRange1[1]

    Tv1 = np.arange(a,b,step);
    cp1 = [];
        
    for i in range(a,b,step):
        #gas.set(T = i, P = P1, X = q);
        gas.TPX = i,P1,q
        #cp1.append(gas.cp_mass());
        #print gas.cp_mole
        cp1.append(gas.cp_mole);

    CpPoly1 = np.polyfit(Tv1, cp1, 4)

    #print 'Polynomials temperature range:'+str(a)+'-'+str(b)+''
    #print CpPoly1

    
    c = tempRange2[0]
    d = tempRange2[1]

    Tv2 = np.arange(c,d,step);
    cp2 = [];
        
    for i in range(c,d,step):
        #gas.set(T = i, P = P1, X = q);
        gas.TPX = i,P1,q
        #cp2.append(gas.cp_mass());
        cp2.append(gas.cp_mole);

    CpPoly2 = np.polyfit(Tv2, cp2, 4)
    
    #print 'Polynomials temperature range:'+str(c)+'-'+str(d)+''
    #print CpPoly2

    print 'Function cpPolynomials teturn list Tv1, cp1, CpPoly1, Tv2, cp2, CpPoly2'
    return [Tv1, cp1, CpPoly1, Tv2, cp2, CpPoly2]

######################################################################

P1 = 100000;
T1 = 298;
q = 'CH4:1'
# xCh4 + (1-x)H2 + (1,5x+0,5)O2+3,76
mech = 'gri30.cti'
tempRange1 = [290, 1000]
tempRange2 = [1000, 5000]
step = 100
out = cpPolynomials(P1,q,mech,tempRange1,tempRange2, step)
print out[2]
print out[5]

[arvindpj]

Originally Posted by sheaker

.....

I found my solution very similar to openFoam predefine coefficients for air as oxidant, octane/propane as fuel and stoichiometric burn products of those fuels.

Have a nice day.
Oskar

Are you computing the JanaF values for burn products using the cantera python code you had mentioned in the earlier post?

Could you explain a bit more?
Thanks.

[sheaker]

Hello.
I'm going to show You how You could deal with it.

1. Calculate Equilibrium state of Your fuel+oxidizer mixture as in example below.

For constant pressure:

import Cantera as ct
import numpy as np
import sys

gas = ct.importPhase('gri30.cti')
gas.set(T=689,P = 1838000,X='C3H8:1 O2:5 N2:18.8')
gas.equilibrate('HP')
print gas.temperature

or for constant volume:

import Cantera as ct
import numpy as np
import sys

gas = ct.importPhase('gri30.cti')
gas.set(T=689,P = 1838000,X='C3H8:1 O2:5 N2:18.8')
gas.equilibrate('UV')
print gas.temperature

2. Check for the important burn products You want to include in Your burn products in openFoam case. Keep in mind that Your chemical mechanism works only for specific species.

3. Prepare Your cpPolynomials.py as in example below.

from Cantera import *
import sys
import numpy as np
from numpy import *
import math
import matplotlib.pyplot as plt


def cpPolynomials(P1,q,mech,tempRange1,tempRange2,step):

    gas = importPhase(mech);
    # specific heat of fresh mixture
    
    a = tempRange1[0]
    b = tempRange1[1]

    Tv1 = np.arange(a,b,step);
    cp1 = [];
        
    for i in range(a,b,step):
        gas.set(T = i, P = P1, X = q);
        cp1.append(gas.cp_mass());

    CpPoly1 = np.polyfit(Tv1, cp1, 4)

    #print 'Polynomials temperature range:'+str(a)+'-'+str(b)+''
    #print CpPoly1

    
    c = tempRange2[0]
    d = tempRange2[1]

    Tv2 = np.arange(c,d,step);
    cp2 = [];
        
    for i in range(c,d,step):
        gas.set(T = i, P = P1, X = q);
        cp2.append(gas.cp_mass());

    CpPoly2 = np.polyfit(Tv2, cp2, 4)
    
    #print 'Polynomials temperature range:'+str(c)+'-'+str(d)+''
    #print CpPoly2

    print 'Function cpPolynomials teturn list Tv1, cp1, CpPoly1, Tv2, cp2, CpPoly2'
    return [Tv1, cp1, CpPoly1, Tv2, cp2, CpPoly2]

######################################################################

P1 = 100000;
T1 = 298;
q = 'O2:7, H2O:95, N2:719, CO:14, CO2:158, NO:4'
# xCh4 + (1-x)H2 + (1,5x+0,5)O2+3,76
mech = 'gri30.cti'
tempRange1 = [250, 1000]
tempRange2 = [1000, 5000]
step = 100
out = cpPolynomials(P1,q,mech,tempRange1,tempRange2, step)
print out[2]
print out[5]

Output:

[ -4.47755028e-10   1.07819449e-06  -9.08567169e-04   6.16503092e-01  9.10859105e+02]
[ -3.14356305e-12   4.25127042e-08  -2.23818856e-04   5.67985575e-01 8.65698780e+02]

For further step lets shorten it to:

[a   b   c   d   e]
[f   g   h   i   j]

4. You need to use tool like excel to make this part easier.
4a. Draw Your polynomials to check if everything is good (if both polynomials looks physical and match each other in T_Common).

For low temperature range

CP= (((a*T+b)*T+c)*T+d)*T+e

For high temperature range

CP= (((f*T+g)*T+h)*T+i)*T+j

4b.Remember that openFoam calculate CP in a little different way. Firstly You need to reverse coefficients, secondly You need to divide solution by R. Check:
https://cfd.direct/openfoam/user-guide/thermophysical/ part (7.4)

4c. I think it could be done in better way but I have not enough time to check it so I'm going to show how I have dealed with it.

4d. Create list of CP(T). Divide values by R. Plot it in X-Y Cartesian coordinate system.

4e. Create a tendline with polynomial of 4th order separately for low temperature range and for high temperature range. Those are values of a0 - a4 and b0 - b4 for openFOAM. (reverse them for openFOAM)

5. You need to check gasProperties to calculate enthalpy and entropy offsets. I'm sharing with You entire code but You only need those values of enthalpy and entropy. They will be at the beginning of the output.

######################
#Basic gas properties#
#Oskar Zuranski, 2017#
#Warsaw, Poland      #
######################
from Cantera import *
import numpy as np
import sys

gas = importPhase('gri30.cti')
gas.set(T=273.16,P=101325,X='O2:7, H2O:95, N2:719, CO:14, CO2:158, NO:4')

gasnames = gas.speciesNames()
gasmass = gas.massFractions()
gasmole = gas.moleFractions()
n=0
print ""
print "BASIC GAS PROPERTIES"
print ""
print ""
print gas
print ""
print ""
print "Temperature:\t", gas.temperature(), "K"
print "Pressure:\t", gas.pressure(), "Pa"

print ""
print "Molar fractions:"
for i in range(gas.nSpecies()):
    if gasmole[i]!=0:
        n+=1    
        if gasmole[i]>=0.1: print n, gasnames[i], "\t{0:.10f}".format(100.0*gasmole[i]), "%"
        else:
            print n, gasnames[i], "\t {0:.10f}".format(100.0*gasmole[i]), "%"

n=0
print ""
print "Mass fractions:"
for i in range(gas.nSpecies()):
    if gasmass[i]!=0:
        n+=1    
        if gasmass[i]>=0.1: print n, gasnames[i], "\t{0:.10f}".format(100.0*gasmass[i]), "%"
        else:
            print n, gasnames[i], "\t {0:.10f}".format(100.0*gasmass[i]), "%"    

#kg/m^3\t{0:.15f}".format(gasmole[i]), "mol/m^3"
print ""
print "Mean molar mass:"
print "M  =", gas.meanMolarMass(), "kg/mol"

print ""
print "Density:"
print "rho =", gas.density(), "kg/m^3"
print "rho =", gas.molarDensity()*1000.0, "mol/m^3"

print ""
print "Specific gas constant:"
print "R  =",GasConstant/gas.meanMolarMass(), "J/kg/K"

print ""
print "Molar heat capacities:"
print "Cp =", gas.cp_mole()/1000.0, "J/mol/K"
print "Cv =", gas.cv_mole()/1000.0, "J/mol/K"


print ""
print "Mass heat capacities:"
print "Cp =", gas.cp_mass(), "J/kg/K"
print "Cv =", gas.cv_mass(), "J/kg/K"

print ""
print "Heat capacity ratio"
print "kappa =", gas.cp_mole()/gas.cv_mole()

In this example it is:

enthalpy    -2.94802e+006      -8.736e+007     J
entropy         6735.43       1.996e+005     J/K

6. Calculate Your CP at temperature of 298.15 K.

7. Calculate offsets of enthalpy and entropy for low temperature:
Enthalpy (from cantera) = CP(T=298.15K) + a5
Entropy (from cantera) = CP(T=298.15K) + a6

7. Calculate offsets of enthalpy and entropy for high temperature:
To match polynomials values in T_Common.


I'm not sure if excel part (4.) is clear enough. If You have any further question do not hesitate to ask.

[arvindpj]

Thanks a ton for your quick detailed reply.
Cheers :-)
Jay

Sent from my SAMSUNG-SM-G920A using CFD Online Forum mobile app

[arvindpj]

Hi,

Here is the complete python code to generate JANAF coeffs for mixtures in OpenFOAM format. I had to make some changes while computing the enthalpy and entropy offsets.
Thanks for your hints and code.

Cheers,
Jay

 
from cantera import *  #Version 2.3.0 (stable)
import sys
import numpy as np
import matplotlib.pyplot as plt
import numpy.polynomial.polynomial as poly


def cpPolynomials(P1,q,mech,tempRange1,tempRange2,step):

    gas = Solution(mech);
    # specific heat of fresh mixture
    
    a = tempRange1[0]
    b = tempRange1[1]

    Tv1 = np.arange(a,b,step);
    cp1 = [];
        
    for i in range(a,b,step):
        gas.TPX = i,P1,q
        cp1.append(gas.cp_mass);  #[J/kg/K]

    CpPoly1 = poly.polyfit(Tv1, cp1, 4)
    #print CpPoly1

    
    c = tempRange2[0]
    d = tempRange2[1]

    Tv2 = np.arange(c,d,step);
    cp2 = [];
        
    for i in range(c,d,step):
        gas.TPX = i,P1,q
        cp2.append(gas.cp_mass);  #[J/kg/K]        

    CpPoly2 = poly.polyfit(Tv2, cp2, 4)
    #print CpPoly2

    #print 'Function cpPolynomials teturn list Tv1, cp1, CpPoly1, Tv2, cp2, CpPoly2'
    return [Tv1, cp1, CpPoly1, Tv2, cp2, CpPoly2]

######################################################################

P1 = 101325;
T1 = 298.15;


q = 'CH4:1'
#q = 'CO2:1'
#q = 'O2:7, H2O:95, N2:719, CO:14, CO2:158, NO:4'
# xCh4 + (1-x)H2 + (1,5x+0,5)O2+3,76


######################################################################


mech = 'gri30.cti'

Tref = 298.15 #K
Tlow = 200
Tcommon = 1000
Thigh = 3500

tempRange1 = [Tlow, Tcommon]
tempRange2 = [Tcommon, Thigh]
step = 100
out = cpPolynomials(P1,q,mech,tempRange1,tempRange2, step)

gas = Solution(mech);
P1 = 101325;
T1 = 298.15;
gas.TPX = T1,P1,q
print gas()


print ""
print "Specific gas constant:"
print "R  =",gas_constant/gas.mean_molecular_weight, "J/kg/K"  
R = gas_constant/gas.mean_molecular_weight
print ""

npoints = 50
temRL = np.linspace(200, 1200, npoints)
temRH = np.linspace(800, 5200, npoints)


cpL = np.zeros(npoints)
cpH = np.zeros(npoints)

cpL = poly.polyval(temRL, out[2])
cpH = poly.polyval(temRH, out[5])

#For OpenFoam
cpLbyR = np.zeros(npoints)
cpHbyR = np.zeros(npoints)

cpLbyR = poly.polyval(temRL, out[2]/R)
cpHbyR = poly.polyval(temRH, out[5]/R)

# From C02 only (for Verification)
#oLbyR = np.zeros(npoints)
#oHbyR = np.zeros(npoints)
#oH = [ 3.85746, 0.00441437, -2.21481e-06, 5.2349e-10, -4.72084e-14 ]
#oL = [ 2.35677, 0.0089846, -7.12356e-06, 2.45919e-09, -1.437e-13 ]
#oLbyR = poly.polyval(temRL, oL)
#oHbyR = poly.polyval(temRH, oH)
#plt.plot(temRL,oLbyR, lw=2)
#plt.plot(temRH,oHbyR, lw=2)


plt.plot(temRL,cpLbyR, lw=2)
plt.plot(temRH,cpHbyR, lw=2)


plt.xlabel('Temperature (K)')
plt.ylabel('$C_p / R$ ')
plt.grid(color='b', alpha=0.5, linewidth=0.5)

#plt.savefig("cp_vs_T.png")
plt.show()


Lcof = out[2]/R
Hcof = out[5]/R

hLoff = Lcof[0] + Lcof[1]*((Tref**1)/2) + Lcof[2]*((Tref**2)/3) + Lcof[3]*((Tref**3)/4) + Lcof[4]*((Tref**4)/5)
sLoff = Lcof[0]*np.log(Tref) + Lcof[1]*((Tref**1)/1) + Lcof[2]*((Tref**2)/2) + Lcof[3]*((Tref**3)/3) + Lcof[4]*((Tref**4)/4)

low_enthalpy_offset = gas.enthalpy_mass/R - hLoff*Tref
low_entropy_offset = gas.entropy_mass/R - sLoff


hHoff = Hcof[0] + Hcof[1]*((Tref**1)/2) + Hcof[2]*((Tref**2)/3) + Hcof[3]*((Tref**3)/4) + Hcof[4]*((Tref**4)/5)
sHoff = Hcof[0]*np.log(Tref)  + Hcof[1]*((Tref**1)/1) + Hcof[2]*((Tref**2)/2) + Hcof[3]*((Tref**3)/3) + Hcof[4]*((Tref**4)/4)

high_enthalpy_offset = gas.enthalpy_mass/R - hHoff*Tref
high_entropy_offset = gas.entropy_mass/R - sHoff

name = 'mixture'

print ""
print "%s" % name
print "{"
print "    specie"
print "    {"
print "        nMoles          1;"
print "        molWeight       %s;" % gas.mean_molecular_weight
print "    }"
print "    thermodynamics"
print "    {"
print "        Tlow            %f;" % Tlow
print "        Thigh           %f;" % Thigh
print "        Tcommon         %f;" % Tcommon
print "        highCpCoeffs    ( %g %g %g %g %g %g %g );" % (Hcof[0], Hcof[1], Hcof[2], Hcof[3], Hcof[4], high_enthalpy_offset, high_entropy_offset)
print "        lowCpCoeffs     ( %g %g %g %g %g %g %g );" % (Lcof[0], Lcof[1], Lcof[2], Lcof[3], Lcof[4], low_enthalpy_offset, low_entropy_offset)
print "    }"
print "    transport"
print "    {"
print "        As              1.67212e-06;"
print "        Ts              170.672;"
print "    }"
print "}"

[sheaker]

Good to see You working with JANAF.
Unfortunately Your script isn't working for me now (different versions of Cantera and Python). I will try to fix it for my usage and post it here As Soon As Possible. For now there is one thing You should check:

out = cpPolynomials(P1,q,mech,tempRange1,tempRange2, step) 

gas = Solution(mech);

It looks like You run cpPolynomials function before redefine gas mechanism. For me this led to usage of default composition with only one mole of hydrogen. I wonder if You had cheked Your output of air composition with the one from openFoam?

Have a nice day.
sheaker

[arvindpj]

Originally Posted by sheaker

Good to see You working with JANAF.
Unfortunately Your script isn't working for me now (different versions of Cantera and Python). I will try to fix it for my usage and post it here As Soon As Possible. For now there is one thing You should check:

out = cpPolynomials(P1,q,mech,tempRange1,tempRange2, step) 

gas = Solution(mech);

It looks like You run cpPolynomials function before redefine gas mechanism. For me this led to usage of default composition with only one mole of hydrogen. I wonder if You had cheked Your output of air composition with the one from openFoam?

Have a nice day.
sheaker

Hi Sheaker.

My code does work and I did verify it with the gri-mech data in the reactionFoam tutorials.

I had used this code in a ipython notebook and cantera 2.3
Please do let me know, if you have any difficulties.

Cheers,
Jay

 

from cantera import *
import sys
import numpy as np
import matplotlib.pyplot as plt
import numpy.polynomial.polynomial as poly


def cpPolynomials(P1,q,mech,tempRange1,tempRange2,step):

    gas = Solution(mech);
    # specific heat of fresh mixture
    
    a = tempRange1[0]
    b = tempRange1[1]

    Tv1 = np.arange(a,b,step);
    cp1 = [];
        
    for i in range(a,b,step):
        gas.TPX = i,P1,q
        cp1.append(gas.cp_mass);  #[J/kg/K]

    CpPoly1 = poly.polyfit(Tv1, cp1, 4)
    #print CpPoly1
    
    c = tempRange2[0]
    d = tempRange2[1]

    Tv2 = np.arange(c,d,step);
    cp2 = [];
        
    for i in range(c,d,step):
        gas.TPX = i,P1,q
        cp2.append(gas.cp_mass);  #[J/kg/K]        

    CpPoly2 = poly.polyfit(Tv2, cp2, 4)
    #print CpPoly2

    #print 'Function cpPolynomials teturn list Tv1, cp1, CpPoly1, Tv2, cp2, CpPoly2'
    return [Tv1, cp1, CpPoly1, Tv2, cp2, CpPoly2]

######################################################################

P1 = 101325;
T1 = 298.15;

#q = 'CO2:1'
# xCh4 + (1-x)H2 + (1,5x+0,5)O2+3,76

#q = 'CH4:1'
#q = 'O2:7, H2O:95, N2:719, CO:14, CO2:158, NO:4'  #Mole fraction
#gas.TPX = 273.16,101325,q

qmass = 'N2:.76, CO2:.13, H2O:.11'  #Flue gas
gas.TPY = 273.16,101325,qmass

######################################################################

mech = 'gri30.cti'

Tref = 298.15 #K
Tlow = 200
Tcommon = 1000
Thigh = 3500

tempRange1 = [Tlow, Tcommon]
tempRange2 = [Tcommon, Thigh]
step = 100
out = cpPolynomials(P1,q,mech,tempRange1,tempRange2, step)

######################################################################

print gas()

gasmass = gas.Y
gasmole = gas.X
gasnames = gas.species_names

print ""
print "Temperature:\t", gas.T, "K"
print "Pressure:\t", gas.P, "Pa"

n=0
print ""
print "Molar fractions:"
for i in range(gas.n_species):
    if gasmole[i]!=0:
        n+=1    
        if gasmole[i]>=0.1: print n, gasnames[i], "\t{0:.10f}".format(100.0*gasmole[i]), "%"
        else:
            print n, gasnames[i], "\t {0:.10f}".format(100.0*gasmole[i]), "%"

n=0
print ""
print "Mass fractions:"
for i in range(gas.n_species):
    if gasmass[i]!=0:
        n+=1    
        if gasmass[i]>=0.1: print n, gasnames[i], "\t{0:.10f}".format(100.0*gasmass[i]), "%"
        else:
            print n, gasnames[i], "\t {0:.10f}".format(100.0*gasmass[i]), "%"

print ""
print "Specific gas constant:"
print "R  =",gas_constant/gas.mean_molecular_weight, "J/kg/K"  #CHheck units ??????????
R = gas_constant/gas.mean_molecular_weight
print ""

npoints = 50
temRL = np.linspace(200, 1200, npoints)
temRH = np.linspace(800, 5200, npoints)

cpL = np.zeros(npoints)
cpH = np.zeros(npoints)

cpL = poly.polyval(temRL, out[2])
cpH = poly.polyval(temRH, out[5])

#For OpenFoam
cpLbyR = np.zeros(npoints)
cpHbyR = np.zeros(npoints)

cpLbyR = poly.polyval(temRL, out[2]/R)
cpHbyR = poly.polyval(temRH, out[5]/R)

# From C02 only (for Verification)
#oLbyR = np.zeros(npoints)
#oHbyR = np.zeros(npoints)
#oH = [ 3.85746, 0.00441437, -2.21481e-06, 5.2349e-10, -4.72084e-14 ]
#oL = [ 2.35677, 0.0089846, -7.12356e-06, 2.45919e-09, -1.437e-13 ]
#oLbyR = poly.polyval(temRL, oL)
#oHbyR = poly.polyval(temRH, oH)
#plt.plot(temRL,oLbyR, lw=2)
#plt.plot(temRH,oHbyR, lw=2)

plt.plot(temRL,cpLbyR, lw=2)
plt.plot(temRH,cpHbyR, lw=2)

plt.xlabel('Temperature (K)')
plt.ylabel('$C_p / R$ ')
plt.grid(color='b', alpha=0.5, linewidth=0.5)

#plt.savefig("cp_vs_T.png")
plt.show()

Lcof = out[2]/R
Hcof = out[5]/R

hLoff = Lcof[0] + Lcof[1]*((Tref**1)/2) + Lcof[2]*((Tref**2)/3) + Lcof[3]*((Tref**3)/4) + Lcof[4]*((Tref**4)/5)
sLoff = Lcof[0]*np.log(Tref) + Lcof[1]*((Tref**1)/1) + Lcof[2]*((Tref**2)/2) + Lcof[3]*((Tref**3)/3) + Lcof[4]*((Tref**4)/4)

low_enthalpy_offset = gas.enthalpy_mass/R - hLoff*Tref
low_entropy_offset = gas.entropy_mass/R - sLoff


hHoff = Hcof[0] + Hcof[1]*((Tref**1)/2) + Hcof[2]*((Tref**2)/3) + Hcof[3]*((Tref**3)/4) + Hcof[4]*((Tref**4)/5)
sHoff = Hcof[0]*np.log(Tref)  + Hcof[1]*((Tref**1)/1) + Hcof[2]*((Tref**2)/2) + Hcof[3]*((Tref**3)/3) + Hcof[4]*((Tref**4)/4)

high_enthalpy_offset = gas.enthalpy_mass/R - hHoff*Tref
high_entropy_offset = gas.entropy_mass/R - sHoff

name = 'mixture'

print '/*--------------------------------*- C++ -*----------------------------------*'+ '/'
print "|=========                 |                                                 |"
print "| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |"
print "|  \\    /   O peration     | Version:  dev                                   |"
print "|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |"
print "|    \\/     M anipulation  |                                                 |"
print "\*---------------------------------------------------------------------------*/"
print "FoamFile"
print "{"
print "    version     2.0;"
print "    format      ascii;"
print "    class       dictionary;"
print '    location    "constant";'
print "    object      thermophysicalProperties;"
print "}"
print "// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //"
print "thermoType"
print "{"
print "    type            heRhoThermo;"
print "    mixture         pureMixture;"
print "    transport       sutherland;"
print "    thermo          janaf;"
print "    equationOfState perfectGas;"
print "    specie          specie;"
print "    energy          sensibleEnthalpy;"
print "}"
print ""
print "%s" % name
print "{"
print "    specie"
print "    {"
print "        nMoles          1;"
print "        molWeight       %s;" % gas.mean_molecular_weight
print "    }"
print "    thermodynamics"
print "    {"
print "        Tlow            %f;" % Tlow
print "        Thigh           %f;" % Thigh
print "        Tcommon         %f;" % Tcommon
print "        highCpCoeffs    ( %g %g %g %g %g %g %g );" % (Hcof[0], Hcof[1], Hcof[2], Hcof[3], Hcof[4], high_enthalpy_offset, high_entropy_offset)
print "        lowCpCoeffs     ( %g %g %g %g %g %g %g );" % (Lcof[0], Lcof[1], Lcof[2], Lcof[3], Lcof[4], low_enthalpy_offset, low_entropy_offset)
print "    }"
print "    transport"
print "    {"
print "        As              1.67212e-06;"
print "        Ts              170.672;"
print "    }"
print "}"

Here is the output:

 

  gri30:

       temperature          273.16  K
          pressure          101325  Pa
           density         1.23277  kg/m^3
  mean mol. weight         27.6322  amu

                          1 kg            1 kmol
                       -----------      ------------
          enthalpy     -2.6664e+06       -7.368e+07     J
   internal energy     -2.7486e+06       -7.595e+07     J
           entropy          7100.7        1.962e+05     J/K
    Gibbs function      -4.606e+06       -1.273e+08     J
 heat capacity c_p          1097.8        3.033e+04     J/K
 heat capacity c_v          796.87        2.202e+04     J/K

                           X                 Y          Chem. Pot. / RT
                     -------------     ------------     ------------
               H2O        0.16872             0.11         -130.982
               CO2      0.0816225             0.13         -201.497
                N2       0.749657             0.76         -23.3349
     [  +50 minor]              0                0

None

Temperature:	273.16 K
Pressure:	101325.0 Pa

Molar fractions:
1 H2O 	16.8720458895 %
2 CO2 	 8.1622527076 %
3 N2 	74.9657014029 %

Mass fractions:
1 H2O 	11.0000000000 %
2 CO2 	13.0000000000 %
3 N2 	76.0000000000 %

Specific gas constant:
R  = 300.897153097 J/kg/K


/*--------------------------------*- C++ -*----------------------------------*/
|=========                 |                                                 |
| \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \    /   O peration     | Version:  dev                                   |
|   \  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      thermophysicalProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
thermoType
{
    type            heRhoThermo;
    mixture         pureMixture;
    transport       sutherland;
    thermo          janaf;
    equationOfState perfectGas;
    specie          specie;
    energy          sensibleEnthalpy;
}

mixture
{
    specie
    {
        nMoles          1;
        molWeight       27.6322391702;
    }
    thermodynamics
    {
        Tlow            200.000000;
        Thigh           3500.000000;
        Tcommon         1000.000000;
        highCpCoeffs    ( 0.128925 0.0230647 -9.87434e-06 2.10638e-09 -1.75368e-13 -9841.96 16.4077 );
        lowCpCoeffs     ( 8.87021 -0.0235471 8.47083e-05 -8.34927e-08 2.87116e-11 -11056.5 -23.004 );
    }
    transport
    {
        As              1.67212e-06;
        Ts              170.672;
    }
}

Here is similar output for CO2

gri30:

       temperature          273.16  K
          pressure          101325  Pa
           density         1.96343  kg/m^3
  mean mol. weight         44.0098  amu

                          1 kg            1 kmol
                       -----------      ------------
          enthalpy     -8.9621e+06       -3.944e+08     J
   internal energy     -9.0137e+06       -3.967e+08     J
           entropy            4785        2.106e+05     J/K
    Gibbs function     -1.0269e+07       -4.519e+08     J
 heat capacity c_p          817.81        3.599e+04     J/K
 heat capacity c_v          628.89        2.768e+04     J/K

                           X                 Y          Chem. Pot. / RT
                     -------------     ------------     ------------
               CO2              1                1         -198.991
     [  +52 minor]              0                0

None

Temperature:	273.16 K
Pressure:	101325.0 Pa

Molar fractions:
1 CO2 	100.0000000000 %

Mass fractions:
1 CO2 	100.0000000000 %

Specific gas constant:
R  = 188.92296943 J/kg/K

/*--------------------------------*- C++ -*----------------------------------*/
|=========                 |                                                 |
| \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \    /   O peration     | Version:  dev                                   |
|   \  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      thermophysicalProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
thermoType
{
    type            heRhoThermo;
    mixture         pureMixture;
    transport       sutherland;
    thermo          janaf;
    equationOfState perfectGas;
    specie          specie;
    energy          sensibleEnthalpy;
}

mixture
{
    specie
    {
        nMoles          1;
        molWeight       44.0098;
    }
    thermodynamics
    {
        Tlow            200.000000;
        Thigh           3500.000000;
        Tcommon         1000.000000;
        highCpCoeffs    ( 3.85746 0.00441437 -2.21481e-06 5.2349e-10 -4.72084e-14 -48765.8 2.12717 );
        lowCpCoeffs     ( 2.35677 0.0089846 -7.12356e-06 2.45919e-09 -1.437e-13 -48481.9 9.51614 );
    }
    transport
    {
        As              1.67212e-06;
        Ts              170.672;
    }
}

Thermo data from GRI-Mech for CO2 (for verification):

CO2
{
    specie
    {
        molWeight       44.01;
    }
    thermodynamics
    {
        Tlow            200;
        Thigh           3500;
        Tcommon         1000;
        highCpCoeffs    ( 3.85746 0.00441437 -2.21481e-06 5.2349e-10 -4.72084e-14 -48759.2 2.27164 );
        lowCpCoeffs     ( 2.35677 0.0089846 -7.12356e-06 2.45919e-09 -1.437e-13 -48372 9.90105 );
    }
    transport
    {
        As              1.67212e-06;
        Ts              170.672;
    }
    elements
    {
        C               1;
        O               2;
    }
}

[sheaker]

I have managed to change Your code for Cantera 2.0.

Here is the script:

from Cantera import *
import sys
import numpy as np


def cpPolynomials(P1,q,mech,tempRange1,tempRange2,step):

    gas = importPhase(mech);
    
    a = tempRange1[0]
    b = tempRange1[1]

    Tv1 = np.arange(a,b,step);
    cp1 = [];
        
    for i in range(a,b,step):
        gas.set(T = i, P = P1, X = q);
        cp1.append(gas.cp_mass());

    CpPoly1 = np.polyfit(Tv1, cp1, 4)
    
    c = tempRange2[0]
    d = tempRange2[1]

    Tv2 = np.arange(c,d,step);
    cp2 = [];
        
    for i in range(c,d,step):
        gas.set(T = i, P = P1, X = q);
        cp2.append(gas.cp_mass());

    CpPoly2 = np.polyfit(Tv2, cp2, 4)

    return [Tv1, cp1, CpPoly1, Tv2, cp2, CpPoly2]

#######################################################################
# USER EDITABLE PART
#######################################################################
mech = 'gri30_highT.cti'

P1 = 100000;
T1 = 298.15;
q = 'CH4:1'

Tlow = 200
Thigh = 6000
Tcommon = 1000
Tref = 298.15

#######################################################################
# END OF USER EDITABLE PART
#######################################################################

gas = importPhase(mech)
gas.set(T=T1,P=P1,X=q)

tempRange1 = [Tlow, Tcommon]
tempRange2 = [Tcommon, Thigh]
step = 100
out = cpPolynomials(P1,q,mech,tempRange1,tempRange2, step)

R = GasConstant/gas.meanMolarMass()

Lcof_rev = out[2]/R
Hcof_rev = out[5]/R

Lcof = list(reversed(Lcof_rev))
Hcof = list(reversed(Hcof_rev))


hLoff = Lcof[0] + Lcof[1]*((Tref**1)/2) + Lcof[2]*((Tref**2)/3) + Lcof[3]*((Tref**3)/4) + Lcof[4]*((Tref**4)/5)
sLoff = Lcof[0]*np.log(Tref) + Lcof[1]*((Tref**1)/1) + Lcof[2]*((Tref**2)/2) + Lcof[3]*((Tref**3)/3) + Lcof[4]*((Tref**4)/4)

low_enthalpy_offset = gas.enthalpy_mass()/R - hLoff*Tref
low_entropy_offset = gas.entropy_mass()/R - sLoff


hHoff = Hcof[0] + Hcof[1]*((Tref**1)/2) + Hcof[2]*((Tref**2)/3) + Hcof[3]*((Tref**3)/4) + Hcof[4]*((Tref**4)/5)
sHoff = Hcof[0]*np.log(Tref)  + Hcof[1]*((Tref**1)/1) + Hcof[2]*((Tref**2)/2) + Hcof[3]*((Tref**3)/3) + Hcof[4]*((Tref**4)/4)

high_enthalpy_offset = gas.enthalpy_mass()/R - hHoff*Tref
high_entropy_offset = gas.entropy_mass()/R - sHoff

print "    specie"
print "    {"
print "        nMoles          1;"
print "        molWeight       %s;" % gas.meanMolarMass()
print "    }"
print "    thermodynamics"
print "    {"
print "        Tlow            %f;" % Tlow
print "        Thigh           %f;" % Thigh
print "        Tcommon         %f;" % Tcommon
print  "        highCpCoeffs    ( %g %g %g %g %g %g %g );" % (Hcof[0],  Hcof[1], Hcof[2], Hcof[3], Hcof[4], high_enthalpy_offset,  high_entropy_offset)
print "        lowCpCoeffs     ( %g %g %g %g %g  %g %g );" % (Lcof[0], Lcof[1], Lcof[2], Lcof[3], Lcof[4],  low_enthalpy_offset, low_entropy_offset)
print "    }"
print "    transport"
print "    {"
print "        As              1.67212e-06;"
print "        Ts              170.672;"
print "    }"
print "}"

I have verified results with CH4 thermophysicalProperties from
OpenFOAM-1.6-ext\tutorials\combustion\fireFoam\les\smallPoolFir e2D\constant\thermophysicalProperties and they match perfectly. (see attachment)

I wonder if You found any more universal reaction mechanism than gri30?
There is nasa.cti thermo data file available but it doesn't contain reaction kinetics.


Have a nice day!
Sheaker

[sheaker]

Curran detailed mechanism includes no less than 6864 reactions and 874 species...
http://www.cerfacs.fr/cantera/docs/m...ran/curran.cti

[foadsf]

There are already some examples:

• Here you may find the thermophysicalProperties of compressible/rhoCentralFoam/biconic25-55Run35/ turorial
• and here another example for multiphase/reactingTwoPhaseEulerFoam/laminar/bubbleColumnEvaporating