[HenrikJohansson]

Hi,

I am working for quite some time to simulate an axial turbine with SES36 as real gas.
SES36 is alot like HFC-365mfc.
I think there is something wrong with my Nasa Polynomials because the rothalpy is not correct. Im running foam-extend 4.0.
Below are is my thermiphysicalProperties.

Is this the correct way to calculate the nasa polynomials for a fluid?
I have entropy, enthalpy and constant heat for a range of temperatures.
I then use least square to fit these to the nasa polynomials.

thermoType      realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<redlichKwong>>>>>;
mixture         SES36 1 184.53 28.5e5 450.7 2111670.0 -9762.24 -83.3023 0.845058 -0.00236983 3.14876e-06 -1.62616e-09 1.42087e-06 0.983296;

Description of coefficients

// ************************************************************************************************************* //
//   Coefficient:	[J/kmol*K]
//   SES36                		--> Name
//   1							--> number of Moles
//   184.53               		--> Molar mass, (Volume, old) [g/mol]
//   28.5e5	          			--> critical pressure	[Pa]
//   450.7               		--> critical temperatur	[K]
//   99263.04465                --> critical density (only for aungierRedlichKwong)	[Molar Volume] = [Molar Mass] / [Critical Density]
//   0.351             			--> acentric factor  (not needed for redlichKwong but for pengRobinson, aungierRedlichKwong and soaveRedlichKwong )
//   1329.5               		--> constant cp of perfect gas (used for constantHeatCapacity) [J/K]
//   49436.5054 --> 2.849677801e-13     --> 7 heat capacity polynomial coefficents --> NASA Heat Capacity Polynomial
//   1.42087e-6 0.983296                   --> two coefficents for sutherlandRealGasTransport, https://www.cfd-online.com/Forums/openfoam-pre-processing/62178-sutherlandtransport-setting-air.html

I used the following code to get my Nasa Polynomials for different pressures.

import numpy as np
from scipy.optimize import leastsq
import matplotlib.pyplot as plt
from tabulate import tabulate
np.set_printoptions(threshold=np.nan)
from SES36_Data import R, M, Import_Data, Get_Data_Names
'''
    ***   INPUT  ***
'''
print("Avalibale SES36 data is:")
names = Get_Data_Names()
print(names)
idx = int(input("Select data by number order, etc. 1,2,3...\r\nSelect: "))-1
print(names[idx] + " selected") 
T, Cp, s, h, Delta_Hf_298_15, h_298_15 = Import_Data(names[idx])

''' Sutherland Transport, Calculate the two, (As, Ts), coefficients '''
mu1 = 15.8257e-6  # Pa*s or kg/m*s
T1 = 126          # K @ 11.121bar
mu2 = 14.2570e-6  # Pa*s or kg/m*s
T2 = 102.62       # K @ 3,74bar

'''
    ***   Calculation  ***
'''

''' Sutherland Transport, Calculate the two coefficients '''
def Sutherland_Transport(mu1, T1, mu2, T2):
    ''' Sutherland Transport equation '''
    ''' mu = As * T^(3/2) / (T + Ts) '''
    rootT1 = np.sqrt(T1)
    mu1rootT2 = mu1*np.sqrt(T2)
    mu2rootT1 = mu2*rootT1
    # Calculate TS:
    Ts = (mu2rootT1 - mu1rootT2)/(mu1rootT2/T1 - mu2rootT1/T2)
    # Calculate As:
    As = mu1*(1.0 + Ts/T1)/rootT1
    return Ts, As  

def NASA_POLYNOM_PLOT_9(pars, T):
    a1, a2, a3, a4, a5, a6, a7, b1, b2 = pars
    y1 = a1/T**2 + a2/T + a3 + a4*T + a5*T**2 + a6*T**3 + a7*T**4
    y2 = -a1/T + a2*np.log(T) + a3*T + a4*T**2/2 + a5*T**3/3 + a6*T**4/4 + a7*T**5/5 + b1
    y3 = -a1/(T**2*2) - a2/T + a3*np.log(T) + a4*T + a5*T**2/2 + a6*T**3/3 + a7*T**4/4 + b2
    return np.array([y1, y2, y3])

def NASA_POLYNOM_9(pars, T, Cp, h, s):
    a1, a2, a3, a4, a5, a6, a7, b1, b2 = pars
    F1 = a1/T**2 + a2/T + a3 + a4*T + a5*T**2 + a6*T**3 + a7*T**4 - Cp
    F2 = -a1/T + a2*np.log(T) + a3*T + a4*T**2/2 + a5*T**3/3 + a6*T**4/4 + a7*T**5/5 + b1 - h
    F3 = -a1/(T**2*2) - a2/T + a3*np.log(T) + a4*T + a5*T**2/2 + a6*T**3/3 + a7*T**4/4 + b2 - s
    return np.concatenate((F1, F2, F3))

def Round(var):
    form = '5g'
    if var.size > 1:
        for i in range(0,var.size):
            var[i] = format(var[i], form)
    else:
        var = format(var, form)
    return var

''' Solve the Sutherland Transport equation '''

Ts, As = Sutherland_Transport(mu1, T1, mu2, T2)
print('\r\nSutherland Transport constants:\r\n\r\nAs = {0}\r\nTs = {1}\r\n'.format(Round(As),Round(Ts)))
#print(As * T1**(3/2) / (T1 + Ts))
#print(As * T2**(3/2) / (T2 + Ts))

''' Solve NASA Polynomials '''

# initial values 
par_init = np.zeros(9)
h = Delta_Hf_298_15 + h - h_298_15    # H(T) = Delta Hf(298) + [ H(T) - H(298) ]
best, cov, info, message, ier = leastsq(NASA_POLYNOM_9,
                                        par_init, args=(T, Cp/R, h/R, s/R),
                                        full_output=True)

print('\n\r{0}\r\n'.format(message))

'''
    ***   Print Result  ***
'''

y_calc1, y_calc2, y_calc3 = NASA_POLYNOM_PLOT_9(best, T) * R
best = Round(best)
print(tabulate( [best], headers=['a1', 'a2', 'a3', 'a4', 'a5', 'a6', 'a7', 'b1', 'b2']))
print('\r\n{0} {1} {2} {3} {4} {5} {6}\r\n'.format(best[0], best[1], best[2], best[3], best[4], best[5], best[6]))
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(12, 6), dpi=80, facecolor='w', edgecolor='k')
fig.subplots_adjust(left=0.05, bottom=None, right=0.98, top=None, wspace=None, hspace=None)
ax1.plot(T, Cp,'or', label = 'Input') 
ax1.plot(T, y_calc1, label = 'Leastsq')
ax1.set_title('Cp')
ax1.legend(loc=4)
ax2.plot(T, h, 'or', label = 'Input')
ax2.plot(T, y_calc2, label = 'Leastsq')
ax2.set_title('h')
ax2.legend(loc=4)
ax3.plot(T, s, 'or', label = 'Input')
ax3.plot(T, y_calc3, label = 'Leastsq')
ax3.set_title('s')
ax3.legend(loc=4)
plt.show()

SES36_Data.py file

import numpy as np
# -----------------------------------------------------------------------------------
'''
    ***   11.12 bar Isobar  ***
'''
# -----------------------------------------------------------------------------------
def Isobar11_12(select):
	if select == "data":
		T = np.arange(125, 230, 5) + 273.15 # K

		Cp = np.array([1.3293,1.3107,1.2943,1.283,1.2753,1.2702,1.267,1.2652,
                     1.2647,1.265,1.266,1.2676,1.2697,1.2723,1.2751,1.2782,
                     1.2815,1.285,1.2887,1.2925,1.2964
                     ]) * 1e-3 * M  # J/kmol*K

		s = np.array([1.6568,1.67,1.6861,1.7018,1.7172,1.7323,1.7472,1.7619,
                    1.7764,1.7908,1.805,1.819,1.8329,1.8467,1.8604,1.874,
                    1.8874,1.9008,1.914,1.9272,1.9402
                    ]) * 1e-3 * M # J/kmol*K

		h = np.array([439.02,444.32,450.83,457.27,463.67,470.03,476.37,482.7,
                    489.03,495.35,501.68,508.01,514.35,520.71,527.08,533.46,
                    539.86,546.28,552.71,559.16,565.64
                    ]) *1e-3 * M # J/kmol

		Delta_Hf_298_15 = 96.48 *1e-3 * M     # J/kmol  <-- Heat of formation
		h_298_15 = 435.12 *1e-3 * M            # J/kmol
		return T, Cp, s, h, Delta_Hf_298_15, h_298_15
	elif select == "name":
		return "Isobar 11.12"
	else:
		print("Somthing went wrong with getting data/ name from" + Isobar11_12("name"))
# -----------------------------------------------------------------------------------
'''
    ***   10 bar Isobar  ***
'''
# -----------------------------------------------------------------------------------
def Isobar10(select):
	if select == "data":
		T = np.arange(120, 230, 5) + 273.15 # K

		Cp = np.array([1.286,1.2726,1.2615,1.254,1.2492,1.2464,1.2451,1.2449,
                     1.2456,1.2471,1.2491,1.2516,1.2544,1.2576,1.2611,1.2647,
                     1.2686,1.2726,1.2767,1.2809,1.2852,1.2896
                     ]) * 1e-3 * M  # J/kmol*K

		s = np.array([1.6505,1.664,1.6798,1.6953,1.7105,1.7256,1.7404,1.755,
                    1.7694,1.7837,1.7979,1.8119,1.8258,1.8396,1.8533,1.8669,
                    1.8803,1.8937,1.9069,1.9201,1.9332,1.9462
                    ]) * 1e-3 * M # J/kmol*K

		h = np.array([435.11,440.47,446.8,453.09,459.35,465.59,471.81,478.04,
                    484.26,490.5,496.74,502.99,509.25,515.53,521.83,528.14,
                    534.48,540.83,547.2,553.6,560.01,566.45
                    ]) *1e-3 * M # J/kmol

		Delta_Hf_298_15 = 96.48 *1e-3 * M     # J/kmol  <-- Heat of formation
		h_298_15 = 435.12 *1e-3 * M            # J/kmol
		return T, Cp, s, h, Delta_Hf_298_15, h_298_15
	elif select == "name":
		return "Isobar10"
	else:
		print("Somthing went wrong with getting data/ name from" + Isobar10("name"))

# -----------------------------------------------------------------------------------
'''
    ***   5 bar Isobar  ***
'''
# -----------------------------------------------------------------------------------
def Isobar5(select):
	if select == "data":
		T = np.arange(90, 230, 5) + 273.15 # K


		Cp = np.array([1.1051, 1.1081, 1.1121, 1.1167, 1.1217, 1.1271, 1.1328, 1.1386, 1.1446,
					   1.1507, 1.1569, 1.1632, 1.1695, 1.1758, 1.1822, 1.1885, 1.1949, 1.2012,
					   1.2075, 1.2137, 1.2199, 1.2261, 1.2323, 1.2384, 1.2444, 1.2504, 1.2564,
					   1.2623]) * 1e-3 * M  # J/kmol*K

		s = np.array([1.6115, 1.6257, 1.6407, 1.6555, 1.6702, 1.6848, 1.6993, 1.7136, 1.7279,
					  1.7420, 1.7561, 1.7700, 1.7839, 1.7977, 1.8114, 1.8250, 1.8385, 1.8519,
					  1.8653, 1.8786, 1.8918, 1.9049, 1.9180, 1.9310, 1.9439, 1.9567, 1.9695,
					  1.9822]) * 1e-3 * M # J/kmol*K

		h = np.array([410.95, 416.16, 421.71, 427.28, 432.88, 438.50, 444.15, 449.83, 455.54,
					  461.28, 467.05, 472.85, 478.68, 484.54, 490.44, 496.36, 502.32, 508.31,
					  514.33, 520.39, 526.47, 532.58, 538.73, 544.91, 551.11, 557.35, 563.62,
					  569.92]) *1e-3 * M # J/kmol

		Delta_Hf_298_15 = 111.57 *1e-3 * M     # J/kmol  <-- Heat of formation
		h_298_15 = 410.95 *1e-3 * M            # J/kmol
		return T, Cp, s, h, Delta_Hf_298_15, h_298_15
	elif select == "name":
		return "Isobar 5"
	else:
		print("Somthing went wrong with getting data/ name from" + Isobar5("name"))

# -----------------------------------------------------------------------------------
'''
    ***   1 bar Isobar  ***
'''
# -----------------------------------------------------------------------------------
def Isobar1(select):
	if select == "data":
		T = np.arange(40, 230, 5) + 273.15 # K

		Cp = np.array([0.9165, 0.9278, 0.9390, 0.9500, 0.9608, 0.9715, 0.9821,
				       0.9925, 1.0027, 1.0128, 1.0228, 1.0326, 1.0422, 1.0517,
				       1.0611, 1.0704, 1.0795, 1.0884, 1.0973, 1.1060, 1.1146,
				       1.1231, 1.1315, 1.1397, 1.1479, 1.1559, 1.1638, 1.1716,
				       1.1793, 1.1869, 1.1943, 1.2017, 1.2089, 1.2161, 1.2231,
				       1.2301, 1.2369, 1.2437]) * 1e-3 * M  # J/kmol*K

		s = np.array([1.5559, 1.5706, 1.5851, 1.5996, 1.6141, 1.6285, 1.6428,
				      1.6571, 1.6713, 1.6855, 1.6996, 1.7136, 1.7276, 1.7416,
				      1.7554, 1.7693, 1.7830, 1.7967, 1.8103, 1.8239, 1.8374,
				      1.8509, 1.8643, 1.8776, 1.8909, 1.9041, 1.9173, 1.9304,
				      1.9435, 1.9564, 1.9694, 1.9822, 1.9950, 2.0078, 2.0205,
				      2.0331, 2.0457, 2.0582]) * 1e-3 * M # J/kmol*K

		h = np.array([370.03, 374.65, 379.31, 384.04, 388.81, 393.64, 398.53,
				      403.46, 408.45, 413.49, 418.58, 423.72, 428.90, 434.14,
				      439.42, 444.75, 450.13, 455.55, 461.01, 466.52, 472.07,
				      477.66, 483.30, 488.98, 494.70, 500.46, 506.26, 512.10,
				      517.97, 523.89, 529.84, 535.83, 541.86, 547.92, 554.02,
				      560.15, 566.32, 572.52]) *1e-3 * M # J/kmol

		Delta_Hf_298_15 = 129.41 *1e-3 * M     # J/kmol  <-- Heat of formation
		h_298_15 = 365.71 *1e-3 * M            # J/kmol
		return T, Cp, s, h, Delta_Hf_298_15, h_298_15
	elif select == "name":
		return "Isobar 1"
	else:
		print("Somthing went wrong with getting data/ name from" + Isobar1("name"))

# -----------------------------------------------------------------------------------
'''
    ***   0.1 bar Isobar  ***
'''
# -----------------------------------------------------------------------------------
def Isobar0_1(select):
	if select == "data":
		T = np.arange(-20, 230, 5) + 273.15 # K

		Cp = np.array([0.7498, 0.7635, 0.777, 0.7904, 0.8035, 0.8165, 0.8294, 0.842, 0.8545,
				       0.8668, 0.879, 0.891, 0.9028, 0.9144, 0.9259, 0.9372, 0.9484, 0.9594,
				       0.9703, 0.981, 0.9915, 1.0019, 1.0122, 1.0223, 1.0323, 1.0421, 1.0518,
				       1.0613, 1.0708, 1.08, 1.0892, 1.0982, 1.1071, 1.1159, 1.1245, 1.133,
				       1.1414, 1.1497, 1.1578, 1.1658, 1.1737, 1.1815, 1.1892, 1.1968, 1.2042,
				       1.2116, 1.2188, 1.2259, 1.2329, 1.2398]) * 1e-3 * M  # J/kmol*K

		s = np.array([1.4887, 1.5035, 1.5183, 1.533, 1.5477, 1.5624, 1.5771, 1.5917, 1.6063,
				      1.6209, 1.6354, 1.6499, 1.6643, 1.6787, 1.693, 1.7074, 1.7216, 1.7358,
				      1.75, 1.7641, 1.7782, 1.7922, 1.8061, 1.82, 1.8339, 1.8477, 1.8615,
				      1.8752, 1.8888, 1.9024, 1.9159, 1.9294, 1.9428, 1.9562, 1.9695, 1.9828,
				      1.996, 2.0091, 2.0222, 2.0353, 2.0482, 2.0612, 2.074, 2.0868, 2.0996,
				      2.1123, 2.1249, 2.1375, 2.15, 2.1625]) * 1e-3 * M # J/kmol*K

		h = np.array([322.43, 326.21, 330.06, 333.98, 337.97, 342.02, 346.13, 350.31, 354.55,
				      358.86, 363.22, 367.65, 372.13, 376.67, 381.27, 385.93, 390.65, 395.42,
				      400.24, 405.12, 410.05, 415.03, 420.07, 425.15, 430.29, 435.48, 440.71,
				      445.99, 451.32, 456.7, 462.13, 467.59, 473.11, 478.66, 484.27, 489.91,
				      495.6, 501.32, 507.09, 512.9, 518.75, 524.64, 530.57, 536.53, 542.53,
				      548.57, 554.65, 560.76, 566.91, 573.09]) *1e-3 * M # J/kmol

		Delta_Hf_298_15 = 142.49 *1e-3 * M     # J/kmol  <-- Heat of formation
		h_298_15 = 321.16 *1e-3 * M            # J/kmol
		return T, Cp, s, h, Delta_Hf_298_15, h_298_15
	elif select == "name":
		return "Isobar 0.1"
	else:
		print("Somthing went wrong with getting data/ name from" + Isobar0_1("name"))
# -----------------------------------------------------------------------------------
'''
    ***   Import Data  ***
'''
# -----------------------------------------------------------------------------------

def Import_Data(name):
   if name == Isobar11_12("name"):
       return Isobar11_12("data")
   elif name == Isobar10("name"):
       return Isobar10("data")
   elif name == Isobar5("name"):
       return Isobar5("data")
   elif name == Isobar1("name"):
       return Isobar1("data")
   elif name == Isobar0_1("name"):
       return Isobar0_1("data")
   else:
       print("\n\r Error!   Name does not exist! Could not import data.\n\r")

def Get_Data_Names():
	return Isobar11_12("name"), Isobar10("name"), Isobar5("name"), Isobar1("name"), Isobar0_1("name")

'''
    ***   Gas Properties  ***
'''
R = 8.314472e-3   # J/kmol*K
M = 184.53     # kg/kmol

[.bastian]

Originally Posted by HenrikJohansson

[...]
I have entropy, enthalpy and constant heat for a range of temperatures.
[...]

Where did you take the data from?
Because CoolProp gives different values.
Using the example of the heat capacity for the 1 bar isobar:

Your values:

Cp = array([0.16941502, 0.17150383, 0.17357415, 0.1756075 , 0.17760388,
            0.17958177, 0.18154119, 0.18346363, 0.18534909, 0.18721608,
            0.18906458, 0.19087611, 0.19265067, 0.19440674, 0.19614433,
            0.19786344, 0.19954557, 0.20119074, 0.20283591, 0.2044441 ,
            0.20603381, 0.20760503, 0.20915777, 0.21067354, 0.21218931,
            0.21366811, 0.21512843, 0.21657026, 0.21799361, 0.21939847,
            0.22076635, 0.22213424, 0.22346517, 0.22479608, 0.22609004,
              0.22738399, 0.22864097, 0.22989795]) # J/kmol*K

CoolProp (CP.PropsSI("CPMOLAR", "T", T, "P", 1e5, "SES36")*1e-3):

Cp = array([0.14436495, 0.14673565, 0.14921751, 0.15180839, 0.15450591,
            0.15730744, 0.1602101 , 0.16321083, 0.16630641, 0.16949345,
            0.17276842, 0.17612771, 0.17956763, 0.1830844 , 0.18667425,
            0.19033333, 0.19405783, 0.19784394, 0.20168785, 0.20558583,
            0.20953415, 0.21352918, 0.21756733, 0.22164509, 0.22575902,
            0.22990578, 0.23408211, 0.23828485, 0.24251091, 0.24675733,
            0.25102123, 0.25529983, 0.25959046, 0.26389054, 0.26819758,
             0.27250922, 0.27682317, 0.28113724]) # J/kmol*K

Accordingly, this could already be the source of the error.

[HenrikJohansson]

Hi .bastian,

Thank you for the input.

The data is from Solkane 9. Software supplied by the manufacturerer Solkatherm of SES36.
The rothalpy problem as I now understand it is a problem from the Mixing Plane.

I have started this old project again and will update here when I have found a solution or any other problem regarding the NASA Polynomials