[sheaker]

Dear All.
I'm still working on engineFoam solver.
I would like to change fuel.
In cobustionProperites I have GuldersCoeffs for Methane, Propane and IsoOctane.
In thermophisicalProperites highCpCoeffs and lowCpCoeffs are calculated for IsoOctane.
I have function cpPolynomials.py for Cantera which returns n polynomial coefficients for isooctane in [J/kgK].

When I create CP(T=300) for IsoOctane from polynomials it is about 1670 J/kgK (and it is correct) but for lowCpCoeffs from thermodynamicalProperties is only about 4,3*R

I understand that for:
lowCpCoeffs [a0 a1 a2 a3 a4 a5 a6]

cp = R((((a4*T+a3)*T+a2)*T+a1)*T+a0)

Does R = 8.315? Then CP(T=300) = 35,7.

What do I wrong?


How to calculate a5 and a6 as Low temperature enthalpy offset and Low temperature entropy offset?

[KarenRei]

Haha, sounds like you're in the same boat as me:

http://www.cfd-online.com/Forums/ope...y-entropy.html

Well, at least we can be ignored together!

Now, that's concerning entropy and enthalpy - I'm still at a loss as to how they're using the offsets. But maybe I can help with the Cp, because I'm calculating reasonable Cp figures in the data I've tried so far. Could you post all of the raw data you're dealing with, such as where you're getting your data for what Cp of isoOctane should be at that temperature and what values you have for a0 through a4?

[sheaker]

Thank You for Your reply.

I'm using
https://github.com/datasqr/PythonSim...Polynomials.py

And I have two mechanism. Both gives same results:
isooctane.cti
octan_highT.cti

http://speedy.sh/cgA3T/mech.zip

I was trying to use nasa_gas.cti but it doesn't work for me.

Specific heat of isooctane:
http://webbook.nist.gov/cgi/cbook.cgi?ID=C111659&Mask=1
cp=188,7 [J/mol*K] = 1693,5 [J/kg*K]

There is mistake in my prev post. Specific heat calculated with cpPolynomials is about 1663 [J/kg*K].


low
-2,37927084e-27 1,09332465e-23 -3,22984875e-10 2,12832096e-06 -5,75993012e-03 8,11136326e+00 -3,06335488e+02

high
-2,14573662e-18 4,16500666e-14 -3,35702764e-10 1,47399791e-06 -3,86785086e-03 6,06659546e+00 4,80104221e+02


Probably Enthalpy and Entropy offsets must be calculated from thermodynamic equations.

Ok. Now I read your entire topic and know there is no answer:/

[KarenRei]

Where are you getting those a0 to a4 parameters from? They look wrong. a0 should be the largest and a4 should be the smallest. For example, for H2 the parameters (from a number of OpenFOAM example chemkin thermo.dat files), converted to OpenFOAM format:


highCpCoeffs ( 2.99142 0.000700064 -5.63383e-08 -9.23158e-12 1.58275e-15 -835.034 -1.35511 );
lowCpCoeffs ( 3.29812 0.000824944 -8.14302e-07 -9.47543e-11 4.13487e-13 -1012.52 -3.29409 );

Note that they start out large and shrink down - all of them are like that. But yours do just the opposite, they start out small and increase.

I'm still (fruitlessly) trying to figure out where on Earth they're getting the a5 and a6 (enthalpy and entropy) parameters. I have noticed a few things. One, while there's a wide range of values over which a5 can be, the high and low a5 coefficients are always similar to each other. Also, it seems that there may (or may not) be a positive correlation between a5 and molecular weight, and there may (or may not) also be a a negative correlation for a5 when the specie is in its reference state (monoatomic elements, diatomic gases).

Yeah, the documentation is terrible :Þ Not sure what to do about it. I'm trying to model aluminum combustion so I too have to add species....

[KarenRei]

I just set this up to try every combination I could think of:

zero = 0.0
low = 200.0
reference = 298.15
common = 1000.0
high = 5000.0

def calculate_enthalpy(T1, T2, a, substract_reference_enthalpy = False):
sub = 0
if subtract_reference_enthalpy:
sub = (a[4]*math.pow(reference, 5)/5 + a[3]*math.pow(reference, 4)/4 + a[2]*math.pow(reference, 3)/3 + a[1]*math.pow(reference, 2)/2 + a[0]*reference)
return R * (a[4]*math.pow(T1, 5)/5 + a[3]*math.pow(T1, 4)/4 + a[2]*math.pow(T1, 3)/3 + a[1]*math.pow(T1, 2)/2 + a[0]*T1 -
(a[4]*math.pow(T2, 5)/5 + a[3]*math.pow(T2, 4)/4 + a[2]*math.pow(T2, 3)/3 + a[1]*math.pow(T2, 2)/2 + a[0]*T2) -
sub
)

for b in [ True, False ]:
for i in [ (zero, zero_enthalpy), (low, low_enthalpy), (reference, reference_enthalpy), (common, common_enthalpy), (high, high_enthalpy) ]:
for j in [ (zero, zero_enthalpy), (low, low_enthalpy), (reference, reference_enthalpy), (common, common_enthalpy), (high, high_enthalpy) ]:
calculated_low_enthalpy = calculate_enthalpy(i[0], j[0], low_results[1], b)
calculated_high_enthalpy = calculate_enthalpy(i[0], j[0], high_results[1], b)

low_enthalpy_offset = ((i[1] - j[1]) - calculated_low_enthalpy) / R
high_enthalpy_offset = ((i[1] - j[1]) - calculated_high_enthalpy) / R

print "%d, %6.2f, %6.2f:\t%6.3f/%6.3f\t%6.3f/%6.3f" % (b, i[0], j[0], low_enthalpy_offset, -11918.1, high_enthalpy_offset, -12131.5)

And ran it on HCl, which has the thermodynamics parameters (from the examples):

Tlow 200;
Thigh 5000;
Tcommon 1000;
highCpCoeffs ( 2.75534 0.00147358 -4.97125e-07 8.10866e-11 -5.07206e-15 -11918.1 6.51512 );
lowCpCoeffs ( 3.33853 0.00126821 -3.66692e-06 4.70399e-09 -1.83601e-12 -12131.5 3.19355 );

So again, note that the values of high and low enthalpy match the print statement. So in theory at least one of those combinations should work, right? I mean, I'm trying to calculate the offset with every possible combination of differences in temperatures. So let's see what I get:

1, 0.00, 0.00: 1027.783/-11918.100 882.766/-12131.500
1, 0.00, 200.00: 2055.658/-11918.100 1804.832/-12131.500
1, 0.00, 298.15: 2054.526/-11918.100 1764.494/-12131.500
1, 0.00, 1000.00: 2054.614/-11918.100 1696.153/-12131.500
1, 0.00, 2000.00: 917.937/-11918.100 1696.457/-12131.500
1, 200.00, 0.00: -0.093/-11918.100 -39.299/-12131.500
1, 200.00, 200.00: 1027.783/-11918.100 882.766/-12131.500
1, 200.00, 298.15: 1026.651/-11918.100 842.428/-12131.500
1, 200.00, 1000.00: 1026.739/-11918.100 774.088/-12131.500
1, 200.00, 2000.00: -109.938/-11918.100 774.391/-12131.500
1, 298.15, 0.00: 1.039/-11918.100 1.039/-12131.500
1, 298.15, 200.00: 1028.915/-11918.100 923.105/-12131.500
1, 298.15, 298.15: 1027.783/-11918.100 882.766/-12131.500
1, 298.15, 1000.00: 1027.871/-11918.100 814.426/-12131.500
1, 298.15, 2000.00: -108.806/-11918.100 814.730/-12131.500
1, 1000.00, 0.00: 0.951/-11918.100 69.379/-12131.500
1, 1000.00, 200.00: 1028.826/-11918.100 991.445/-12131.500
1, 1000.00, 298.15: 1027.694/-11918.100 951.107/-12131.500
1, 1000.00, 1000.00: 1027.783/-11918.100 882.766/-12131.500
1, 1000.00, 2000.00: -108.895/-11918.100 883.070/-12131.500
1, 2000.00, 0.00: 1137.628/-11918.100 69.076/-12131.500
1, 2000.00, 200.00: 2165.504/-11918.100 991.142/-12131.500
1, 2000.00, 298.15: 2164.372/-11918.100 950.803/-12131.500
1, 2000.00, 1000.00: 2164.460/-11918.100 882.463/-12131.500
1, 2000.00, 2000.00: 1027.783/-11918.100 882.766/-12131.500
0, 0.00, 0.00: 0.000/-11918.100 0.000/-12131.500
0, 0.00, 200.00: 1027.875/-11918.100 922.066/-12131.500
0, 0.00, 298.15: 1026.743/-11918.100 881.727/-12131.500
0, 0.00, 1000.00: 1026.832/-11918.100 813.387/-12131.500
0, 0.00, 2000.00: -109.845/-11918.100 813.690/-12131.500
0, 200.00, 0.00: -1027.875/-11918.100 -922.066/-12131.500
0, 200.00, 200.00: 0.000/-11918.100 0.000/-12131.500
0, 200.00, 298.15: -1.132/-11918.100 -40.338/-12131.500
0, 200.00, 1000.00: -1.044/-11918.100 -108.679/-12131.500
0, 200.00, 2000.00: -1137.721/-11918.100 -108.375/-12131.500
0, 298.15, 0.00: -1026.743/-11918.100 -881.727/-12131.500
0, 298.15, 200.00: 1.132/-11918.100 40.338/-12131.500
0, 298.15, 298.15: 0.000/-11918.100 0.000/-12131.500
0, 298.15, 1000.00: 0.088/-11918.100 -68.340/-12131.500
0, 298.15, 2000.00: -1136.589/-11918.100 -68.037/-12131.500
0, 1000.00, 0.00: -1026.832/-11918.100 -813.387/-12131.500
0, 1000.00, 200.00: 1.044/-11918.100 108.679/-12131.500
0, 1000.00, 298.15: -0.088/-11918.100 68.340/-12131.500
0, 1000.00, 1000.00: 0.000/-11918.100 0.000/-12131.500
0, 1000.00, 2000.00: -1136.677/-11918.100 0.304/-12131.500
0, 2000.00, 0.00: 109.845/-11918.100 -813.690/-12131.500
0, 2000.00, 200.00: 1137.721/-11918.100 108.375/-12131.500
0, 2000.00, 298.15: 1136.589/-11918.100 68.037/-12131.500
0, 2000.00, 1000.00: 1136.677/-11918.100 -0.304/-12131.500
0, 2000.00, 2000.00: 0.000/-11918.100 0.000/-12131.500

All combinations were at least an order of magnitude below the value reported in the thermomdynamic data supplied with the OpenFOAM examples.

Sooo.... how exactly are we supposed to get these figures? :Þ I just don't get it at all.

Here's another thing: check out the thermodynamics data about H2 and HCl. First, Cp:

http://citrination.com/uploads/20/sa...city-data-type
http://citrination.com/uploads/20/sa...rimental-state

Now enthalpy:

http://citrination.com/uploads/20/sa...alpy-data-type
http://citrination.com/uploads/20/sa...alpy-data-type

Pretty similar, right? So why is the enthalpy offset a5 *12 to 15 times higher* in the example thermodynamics data for HCl than H2?

It's just maddening trying to figure out what they're actually doing here...

[sheaker]

Originally Posted by KarenRei

Where are you getting those a0 to a4 parameters from? They look wrong. a0 should be the largest and a4 should be the smallest. For example, for H2 the parameters (from a number of OpenFOAM example chemkin thermo.dat files), converted to OpenFOAM format:


highCpCoeffs ( 2.99142 0.000700064 -5.63383e-08 -9.23158e-12 1.58275e-15 -835.034 -1.35511 );
lowCpCoeffs ( 3.29812 0.000824944 -8.14302e-07 -9.47543e-11 4.13487e-13 -1012.52 -3.29409 );

Note that they start out large and shrink down - all of them are like that. But yours do just the opposite, they start out small and increase.

I'm still (fruitlessly) trying to figure out where on Earth they're getting the a5 and a6 (enthalpy and entropy) parameters. I have noticed a few things. One, while there's a wide range of values over which a5 can be, the high and low a5 coefficients are always similar to each other. Also, it seems that there may (or may not) be a positive correlation between a5 and molecular weight, and there may (or may not) also be a a negative correlation for a5 when the specie is in its reference state (monoatomic elements, diatomic gases).

Yeah, the documentation is terrible :Þ Not sure what to do about it. I'm trying to model aluminum combustion so I too have to add species....

It is because I show You results of np.polyfit and it returns polynomials coefficients it other directions than in lowCpCoeffs.


KarenRei, thank You for Your work on this offsets. It looks really strange. Why do they calculated it in so complicated way?

[KarenRei]

So, I dug into the source code. Here's the janaf equations:

template<class EquationOfState>
inline Foam::scalar Foam::janafThermo<EquationOfState>::ha(const scalar p, const scalar T) const
{
const coeffArray& a = coeffs(T);
return RR*(((((a[4]/5.0*T + a[3]/4.0)*T + a[2]/3.0)*T + a[1]/2.0)*T + a[0])*T + a[5]);
}

template<class EquationOfState>
inline Foam::scalar Foam::janafThermo<EquationOfState>::hs(const scalar p, const scalar T ) const
{
return ha(p, T) - hc();
}

template<class EquationOfState>
inline Foam::scalar Foam::janafThermo<EquationOfState>::hc() const
{
const coeffArray& a = lowCpCoeffs_;
return RR*(((((a[4]/5.0*Tstd + a[3]/4.0)*Tstd + a[2]/3.0)*Tstd + a[1]/2.0)*Tstd + a[0])*Tstd + a[5]);
}

template<class EquationOfState>
inline Foam::scalar Foam::janafThermo<EquationOfState>::s(const scalar p, const scalar T) const
{
const coeffArray& a = coeffs(T);
return RR*((((a[4]/4.0*T + a[3]/3.0)*T + a[2]/2.0)*T + a[1])*T + a[0]*log(T) + a[6]) + EquationOfState::s(p, T);
}


Great... except that it still doesn't explain what they're pulling these entropy / enthalpy values out of. It says how they're using them, but not how they're getting them, and thus still doesn't help me.

But you know what, I don't care. Because I can see how they're using them, I'm going to pick values that make sense in that context, regardless of how [i]they've[i] been picking numbers to use in there. a[5] (enthalpy) should be just the average difference between calculated enthalpy and real enthalpy in the range from low to common for the low params and from common to high for the high params, divided by R in each case. The exact same for enthalpy.

Ugh, I hope I'll choose the units the same as them....

[sheaker]

Originally Posted by KarenRei

So, I dug into the source code. Here's the janaf equations:

template<class EquationOfState>
inline Foam::scalar Foam::janafThermo<EquationOfState>::ha(const scalar p, const scalar T) const
{
const coeffArray& a = coeffs(T);
return RR*(((((a[4]/5.0*T + a[3]/4.0)*T + a[2]/3.0)*T + a[1]/2.0)*T + a[0])*T + a[5]);
}

template<class EquationOfState>
inline Foam::scalar Foam::janafThermo<EquationOfState>::hs(const scalar p, const scalar T ) const
{
return ha(p, T) - hc();
}

template<class EquationOfState>
inline Foam::scalar Foam::janafThermo<EquationOfState>::hc() const
{
const coeffArray& a = lowCpCoeffs_;
return RR*(((((a[4]/5.0*Tstd + a[3]/4.0)*Tstd + a[2]/3.0)*Tstd + a[1]/2.0)*Tstd + a[0])*Tstd + a[5]);
}

template<class EquationOfState>
inline Foam::scalar Foam::janafThermo<EquationOfState>::s(const scalar p, const scalar T) const
{
const coeffArray& a = coeffs(T);
return RR*((((a[4]/4.0*T + a[3]/3.0)*T + a[2]/2.0)*T + a[1])*T + a[0]*log(T) + a[6]) + EquationOfState::s(p, T);
}


Great... except that it still doesn't explain what they're pulling these entropy / enthalpy values out of. It says how they're using them, but not how they're getting, and it still doesn't help me figure out from what hat they're pulling these numbers.

But you know what, I don't care. Because I can see how they're using them, I'm going to pick values that make sense in that context, regardless of how [i]they've[i] been picking numbers to use in there. a[5] (enthalpy) should be just the average difference between calculated enthalpy and real enthalpy in the range from low to common for the low params and from common to high for the high, divided by R. The exact same for enthalpy.

Ugh, I hope I'll choose the units the same as them....

Thank You for this explanation. I'm going to study it and check if I can calculate those offsets for my own case.

Happy New Year!

[KarenRei]

You too!

[KarenRei]

So aaaaanyway.... I think I've written a python script to generate janaf parameters:

http://pastebin.com/85C17qTa

Check it out. Note the functions calculate_specific_heat, calculate_enthalpy and calculate_entropy. They're the exact same as is used in /src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H except for the a[5] and a[6] parameters, which it's calculating. Then note the verificaton code at the bottom adds the effect of the a[5] and a[6] parameters back in. So, for example, on a citrination datafile for H2 I get:

Errors: low=0.001820, high=0.002440
Specific heat checks, low: 27.52 27.45 / 30.20 30.20
Specific heat checks, high: 30.17 30.20 / 42.06 41.97
Enthalpy checks, low: -2772.68 -2774.00 / 20684.18 20680.00
Enthalpy checks, low: 20643.31 20680.00 / 208344.87 208341.00
Entropy checks, low: 119.42 119.41 / 166.22 166.22
Entropy checks, low: 166.18 166.22 / 230.32 230.32
H2
{
specie
{
nMoles 1;
molWeight 2.01588;
}
thermodynamics
{
Tlow 200.000000;
Thigh 6000.000000;
Tcommon 1000.000000;
highCpCoeffs ( 2.90499 0.000875347 -1.70915e-07 2.03727e-11 -1.03855e-15 -807.758 -0.876266 );
lowCpCoeffs ( 2.4187 0.00727126 -1.73593e-05 1.76502e-08 -6.34839e-12 -923.003 0.396006 );
}
}

So you can see that the figures it calculates come in very close to the actual, using the formulas OpenFOAM uses. And you can see, in the thermodynamics data that comes with OpenFOAM:

H2
{
specie
{
nMoles 1;
molWeight 2.01594;
}
thermodynamics
{
Tlow 200;
Thigh 6000;
Tcommon 1000;
highCpCoeffs ( 2.93287 0.000826608 -1.46402e-07 1.541e-11 -6.88805e-16 -813.066 -1.02433 );
lowCpCoeffs ( 2.34433 0.00798052 -1.94782e-05 2.01572e-08 -7.37612e-12 -917.935 0.68301 );
}
}


Very similar to what we generate! Now, with HCl we still get differences:

Errors: low=0.000015, high=0.000493
Specific heat checks, low: 29.13 29.12 / 31.63 31.63
Specific heat checks, high: 31.64 31.63 / 39.18 39.28
Enthalpy checks, low: -2859.19 -2859.00 / 21046.94 21046.00
Enthalpy checks, low: 21054.50 21046.00 / 206688.97 206693.00
Entropy checks, low: 175.27 175.27 / 222.90 222.90
Entropy checks, low: 222.91 222.90 / 287.57 287.57
HCl
{
specie
{
nMoles 1;
molWeight 36.46094;
}
thermodynamics
{
Tlow 200.000000;
Thigh 6000.000000;
Tcommon 1000.000000;
highCpCoeffs ( 2.78657 0.00140237 -4.48212e-07 6.81079e-11 -3.90719e-15 -822.32 6.36106 );
lowCpCoeffs ( 3.46745 0.000440939 -1.89332e-06 3.18209e-09 -1.39324e-12 -1042.33 2.65022 );
}
}

vs the supplied...

HCl
{
specie
{
nMoles 1;
molWeight 36.461;
}
thermodynamics
{
Tlow 200;
Thigh 6000;
Tcommon 1000;
highCpCoeffs ( 2.75758 0.00145387 -4.79647e-07 7.77909e-11 -4.79574e-15 -11913.8 6.52197 );
lowCpCoeffs ( 3.46376 0.000476484 -2.00301e-06 3.31714e-09 -1.44958e-12 -12144.4 2.66428 );
}
}

You'll see that we're very similar except that they for no obvious reason use huge enthalpy offsets here. From looking at the code, I think they're wrong. They're at the very least not consistent.

[KarenRei]

Here's an updated version - more elements supported, etc:

http://pastebin.com/1DW6nx2P

[sheaker]

Great job, KarenRei! Thank You again for Your investigation. Thanks to people like You openFoam become more popular!

I'm using openFoam 2.2.1 now but and I know that openfoam 2.1.1 is bugged (some thernodynamics equation are wrong). I think they could make other mistakes like with enthalpy offset for HCL.

I'm going to check Your code for others gases I found in tutorials like for CH4, C2H6,C3H8, O2+3,76N2.

[KarenRei]

I wouldn't be surprised if it's bugged. And note from their code, in calls of hs(), the a[5] parameter doesn't even matter - it gets cancelled out! It only ever affects ha(). But most people operate in "sensible enthalpy" mode (hs), that's what's recommended by all of the docs. And even in the case of ha being called, since it's the difference between enthalpies that matters, not absolute values, comparisons between ha() on the same material will also cancel out a[5]. So bad values of a[5] will often go unnoticed even on ha().

I'm just glad that I can finally move on to adding in my missing aluminum species and get on to debugging my model! Hope that the arrhenius parameters I pieced together are accurate... your case is probably easier than mine as long hydrocarbons by and large just break down to simpler hydrocarbons, and their reactions with air are well modeled; in my case no two papers can seem to agree on exactly what species of aluminum oxides/hydroxides occur when and how they react with each other, and so you can't just fill in the gaps in one set of parameters from one paper with data from another paper (it doesn't even burn completely in the real world as it entombs itself in oxide as it burns; a droplet's reaction stops as soon as it loses the ability to crack or otherwise seep out of its crust)

[sheaker]

I wish You luck and all You need!

[KarenRei]

You too!

[AgnesD]

Hi KarenRei,

Thank you for sharing your Python script. I am new to OpenFOAM and your script will be very useful as I am currently trying to simulate reactions involving Ar and some hydrocarbons. Your script seems to require an input file. I tried to work out how that's supposed to look like but it doesn't appear to be straightforward. Could you please share an example file? It would be extremely useful.

Many thanks in advance,

Agnes