|
[Sponsors] |
how to get Sutherland and JANAF coefficients of air? |
![]() |
|
LinkBack | Thread Tools | Search this Thread | Display Modes |
![]() |
![]() |
#41 |
Senior Member
Oskar
Join Date: Nov 2015
Location: Poland
Posts: 184
Rep Power: 10 ![]() |
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 |
|
![]() |
![]() |
![]() |
![]() |
#42 | |
New Member
Join Date: Jan 2017
Posts: 24
Rep Power: 9 ![]() |
Quote:
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 |
||
![]() |
![]() |
![]() |
![]() |
#43 |
Senior Member
Oskar
Join Date: Nov 2015
Location: Poland
Posts: 184
Rep Power: 10 ![]() |
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 |
|
![]() |
![]() |
![]() |
![]() |
#44 | |
New Member
Join Date: Jan 2017
Posts: 24
Rep Power: 9 ![]() |
Quote:
|
||
![]() |
![]() |
![]() |
![]() |
#45 |
Senior Member
Oskar
Join Date: Nov 2015
Location: Poland
Posts: 184
Rep Power: 10 ![]() |
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 |
|
![]() |
![]() |
![]() |
![]() |
#46 | |
New Member
Join Date: Jan 2017
Posts: 24
Rep Power: 9 ![]() |
Quote:
References: https://ntrs.nasa.gov/archive/nasa/c...9940013151.pdf http://garfield.chem.elte.hu/Burcat/THERM.DAT |
||
![]() |
![]() |
![]() |
![]() |
#47 | |
New Member
Join Date: Jan 2017
Posts: 24
Rep Power: 9 ![]() |
Quote:
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 |
||
![]() |
![]() |
![]() |
![]() |
#48 |
Senior Member
Oskar
Join Date: Nov 2015
Location: Poland
Posts: 184
Rep Power: 10 ![]() |
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! ![]() |
|
![]() |
![]() |
![]() |
![]() |
#49 | |
New Member
Join Date: Jan 2017
Posts: 24
Rep Power: 9 ![]() |
Quote:
I didn't do the correct calculations! My values actually work ![]() ![]() Thanks for your kind and prompt replies! |
||
![]() |
![]() |
![]() |
![]() |
#50 |
Senior Member
Oskar
Join Date: Nov 2015
Location: Poland
Posts: 184
Rep Power: 10 ![]() |
Now it looks good!
I wish You quick progress in Your simulation! Have a nice day. Sheaker |
|
![]() |
![]() |
![]() |
![]() |
#51 | |
Member
Arvind Jay
Join Date: Sep 2012
Posts: 96
Rep Power: 15 ![]() |
Quote:
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 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] |
||
![]() |
![]() |
![]() |
![]() |
#52 | |
Member
Arvind Jay
Join Date: Sep 2012
Posts: 96
Rep Power: 15 ![]() |
Quote:
Could you explain a bit more? Thanks. |
||
![]() |
![]() |
![]() |
![]() |
#53 |
Senior Member
Oskar
Join Date: Nov 2015
Location: Poland
Posts: 184
Rep Power: 10 ![]() |
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 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 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] 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] Code:
[a b c d e] [f g h i j] 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 Code:
CP= (((f*T+g)*T+h)*T+i)*T+j 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() Code:
enthalpy -2.94802e+006 -8.736e+007 J entropy 6735.43 1.996e+005 J/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. |
|
![]() |
![]() |
![]() |
![]() |
#54 |
Member
Arvind Jay
Join Date: Sep 2012
Posts: 96
Rep Power: 15 ![]() |
Thanks a ton for your quick detailed reply.
Cheers :-) Jay Sent from my SAMSUNG-SM-G920A using CFD Online Forum mobile app |
|
![]() |
![]() |
![]() |
![]() |
#55 |
Member
Arvind Jay
Join Date: Sep 2012
Posts: 96
Rep Power: 15 ![]() |
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 "}" |
|
![]() |
![]() |
![]() |
![]() |
#56 |
Senior Member
Oskar
Join Date: Nov 2015
Location: Poland
Posts: 184
Rep Power: 10 ![]() |
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); Have a nice day. sheaker |
|
![]() |
![]() |
![]() |
![]() |
#57 | |
Member
Arvind Jay
Join Date: Sep 2012
Posts: 96
Rep Power: 15 ![]() |
Quote:
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; } } 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; } } |
||
![]() |
![]() |
![]() |
![]() |
#58 |
Senior Member
Oskar
Join Date: Nov 2015
Location: Poland
Posts: 184
Rep Power: 10 ![]() |
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 "}" 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 |
|
![]() |
![]() |
![]() |
![]() |
#59 |
Senior Member
Oskar
Join Date: Nov 2015
Location: Poland
Posts: 184
Rep Power: 10 ![]() |
Curran detailed mechanism includes no less than 6864 reactions and 874 species...
http://www.cerfacs.fr/cantera/docs/m...ran/curran.cti |
|
![]() |
![]() |
![]() |
![]() |
#60 |
Member
Foad
Join Date: Aug 2017
Posts: 58
Rep Power: 9 ![]() |
There are already some examples:
Last edited by foadsf; August 30, 2018 at 09:51. |
|
![]() |
![]() |
![]() |
Tags |
janaf, openfoam, sutherland |
|
|