CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

how to get Sutherland and JANAF coefficients of air?

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree17Likes

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   May 10, 2017, 07:19
Default
  #41
Senior Member
 
sheaker's Avatar
 
Oskar
Join Date: Nov 2015
Location: Poland
Posts: 184
Rep Power: 10
sheaker is on a distinguished road
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
sheaker is offline   Reply With Quote

Old   May 22, 2017, 07:08
Default
  #42
New Member
 
Join Date: Jan 2017
Posts: 24
Rep Power: 9
er99 is on a distinguished road
Quote:
Originally Posted by sheaker View Post
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
er99 is offline   Reply With Quote

Old   May 22, 2017, 12:57
Default
  #43
Senior Member
 
sheaker's Avatar
 
Oskar
Join Date: Nov 2015
Location: Poland
Posts: 184
Rep Power: 10
sheaker is on a distinguished road
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
sheaker is offline   Reply With Quote

Old   May 22, 2017, 17:55
Default
  #44
New Member
 
Join Date: Jan 2017
Posts: 24
Rep Power: 9
er99 is on a distinguished road
Quote:
Originally Posted by sheaker View Post
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?
er99 is offline   Reply With Quote

Old   May 22, 2017, 18:26
Default
  #45
Senior Member
 
sheaker's Avatar
 
Oskar
Join Date: Nov 2015
Location: Poland
Posts: 184
Rep Power: 10
sheaker is on a distinguished road
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
Attached Images
File Type: png CPofAir.png (13.2 KB, 148 views)
sheaker is offline   Reply With Quote

Old   May 23, 2017, 08:53
Default
  #46
New Member
 
Join Date: Jan 2017
Posts: 24
Rep Power: 9
er99 is on a distinguished road
Quote:
Originally Posted by sheaker View Post
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 is offline   Reply With Quote

Old   May 31, 2017, 12:15
Default
  #47
New Member
 
Join Date: Jan 2017
Posts: 24
Rep Power: 9
er99 is on a distinguished road
Quote:
Originally Posted by sheaker View Post
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
Attached Images
File Type: jpg cp_air_zoomed.jpg (46.3 KB, 76 views)
er99 is offline   Reply With Quote

Old   May 31, 2017, 13:22
Default
  #48
Senior Member
 
sheaker's Avatar
 
Oskar
Join Date: Nov 2015
Location: Poland
Posts: 184
Rep Power: 10
sheaker is on a distinguished road
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!
sheaker is offline   Reply With Quote

Old   May 31, 2017, 13:26
Default
  #49
New Member
 
Join Date: Jan 2017
Posts: 24
Rep Power: 9
er99 is on a distinguished road
Quote:
Originally Posted by sheaker View Post
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!
Attached Images
File Type: jpg cp_air_full.jpg (54.5 KB, 93 views)
File Type: jpg cp_air_zoomed.jpg (49.5 KB, 64 views)
er99 is offline   Reply With Quote

Old   May 31, 2017, 14:29
Default
  #50
Senior Member
 
sheaker's Avatar
 
Oskar
Join Date: Nov 2015
Location: Poland
Posts: 184
Rep Power: 10
sheaker is on a distinguished road
Now it looks good!

I wish You quick progress in Your simulation!

Have a nice day.
Sheaker
sheaker is offline   Reply With Quote

Old   June 13, 2017, 11:19
Default
  #51
Member
 
Arvind Jay
Join Date: Sep 2012
Posts: 96
Rep Power: 14
arvindpj is on a distinguished road
Quote:
Originally Posted by sheaker View Post
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.

Code:
 

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.

Code:
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 is offline   Reply With Quote

Old   June 13, 2017, 11:23
Default
  #52
Member
 
Arvind Jay
Join Date: Sep 2012
Posts: 96
Rep Power: 14
arvindpj is on a distinguished road
Quote:
Originally Posted by sheaker View Post
.....

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.
arvindpj is offline   Reply With Quote

Old   June 13, 2017, 13:43
Default
  #53
Senior Member
 
sheaker's Avatar
 
Oskar
Join Date: Nov 2015
Location: Poland
Posts: 184
Rep Power: 10
sheaker is on a distinguished road
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:
Code:
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:
Code:
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.

Code:
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:
Code:
[ -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:
Code:
[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
Code:
CP= (((a*T+b)*T+c)*T+d)*T+e
For high temperature range
Code:
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.

Code:
######################
#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:
Code:
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 likes this.
sheaker is offline   Reply With Quote

Old   June 13, 2017, 13:51
Default
  #54
Member
 
Arvind Jay
Join Date: Sep 2012
Posts: 96
Rep Power: 14
arvindpj is on a distinguished road
Thanks a ton for your quick detailed reply.
Cheers :-)
Jay

Sent from my SAMSUNG-SM-G920A using CFD Online Forum mobile app
arvindpj is offline   Reply With Quote

Old   July 17, 2017, 15:34
Default JANAF Coeff for Mixtures
  #55
Member
 
Arvind Jay
Join Date: Sep 2012
Posts: 96
Rep Power: 14
arvindpj is on a distinguished road
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

Code:
 
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 "}"
kcjarvis56 and francescomarra like this.
arvindpj is offline   Reply With Quote

Old   July 17, 2017, 17:26
Default
  #56
Senior Member
 
sheaker's Avatar
 
Oskar
Join Date: Nov 2015
Location: Poland
Posts: 184
Rep Power: 10
sheaker is on a distinguished road
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:

Code:
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
sheaker is offline   Reply With Quote

Old   July 17, 2017, 18:15
Default JANAF Coeff for Flue gas
  #57
Member
 
Arvind Jay
Join Date: Sep 2012
Posts: 96
Rep Power: 14
arvindpj is on a distinguished road
Quote:
Originally Posted by sheaker View Post
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:

Code:
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


Code:
 

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:

Code:
 

  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

Code:
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):

Code:
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;
    }
}
Attached Images
File Type: jpg CO2_Compared_cp_vs_T.jpg (45.4 KB, 96 views)
arvindpj is offline   Reply With Quote

Old   July 18, 2017, 15:09
Default
  #58
Senior Member
 
sheaker's Avatar
 
Oskar
Join Date: Nov 2015
Location: Poland
Posts: 184
Rep Power: 10
sheaker is on a distinguished road
I have managed to change Your code for Cantera 2.0.

Here is the script:
Code:
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
Attached Images
File Type: png janafscript.png (21.8 KB, 79 views)
arvindpj likes this.
sheaker is offline   Reply With Quote

Old   July 27, 2017, 18:40
Default
  #59
Senior Member
 
sheaker's Avatar
 
Oskar
Join Date: Nov 2015
Location: Poland
Posts: 184
Rep Power: 10
sheaker is on a distinguished road
Curran detailed mechanism includes no less than 6864 reactions and 874 species...
http://www.cerfacs.fr/cantera/docs/m...ran/curran.cti
arvindpj and dzordz like this.
sheaker is offline   Reply With Quote

Old   August 28, 2018, 17:21
Default
  #60
Member
 
Foad
Join Date: Aug 2017
Posts: 58
Rep Power: 8
foadsf is on a distinguished road
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

Last edited by foadsf; August 30, 2018 at 09:51.
foadsf is offline   Reply With Quote

Reply

Tags
janaf, openfoam, sutherland

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On



All times are GMT -4. The time now is 08:12.