Skip to content

Commit

Permalink
Catching up to Cantera with viscosity polyl fitting scaling and disco…
Browse files Browse the repository at this point in the history
…vering causes of discrepencies in thermo properties between ct and pg.
  • Loading branch information
kaschau committed Jun 22, 2024
1 parent 4d5310c commit 4f71dbf
Show file tree
Hide file tree
Showing 3 changed files with 48 additions and 14 deletions.
8 changes: 5 additions & 3 deletions src/compute/transport/kineticTheory.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -93,7 +93,8 @@ void kineticTheory(block_ &b, const thtrdat_ &th, const int &nface,

// Evaluate all property polynomials
const double logT = log(T);
const double sqrt_T = exp(0.5 * logT);
const double sqrt_T = sqrt(T);
const double sqrtsqrt_T = sqrt(sqrt_T);
double logT_n[deg + 1];
logT_n[0] = 1.0;
for (int ply = 1; ply <= deg; ply++) {
Expand All @@ -112,9 +113,10 @@ void kineticTheory(block_ &b, const thtrdat_ &th, const int &nface,
}

// Set to the correct dimensions
mu_sp(n) *= sqrt_T;
mu_sp(n) *= sqrtsqrt_T;
mu_sp(n) *= mu_sp(n);
kappa_sp(n) *= sqrt_T;
const double T_3o2 = pow(T, 1.5);
const double T_3o2 = T * sqrt_T;
for (int n2 = n; n2 <= ns - 1; n2++) {
Dij(n, n2) *= T_3o2;
Dij(n2, n) = Dij(n, n2);
Expand Down
36 changes: 28 additions & 8 deletions src/peregrinepy/thermoTransport/kineticTheoryPoly.py
Original file line number Diff line number Diff line change
Expand Up @@ -30,15 +30,25 @@ def cp_R(cp0, poly, T, MW):

deg = 4
# Maximum and minumum temperatures to generate poly'l
# NOTE: These ranges vary by input file in Cantera. It seems to set the
# minTemp and maxTemp based on the min/max ranges of the NASA7 poly'l
# data. In testing, this seems to explain the errors we sometimes get
# in thermodynamic testing against Cantera.
# (i.e. takes error from 1% to 0.001%). For now we just use sensible values
# here.
Tmin = 200
Tmax = 5000
# Generate range of temperatures
npts = 50
Ts = np.linspace(Tmin, Tmax, npts)

# Collision integral interpolations
intrp_o22 = intrp.RectBivariateSpline(tstar22, delta, omega22_table, kx=5, ky=5)
intrp_Astar = intrp.RectBivariateSpline(tstar, delta, astar_table, kx=5, ky=5)
intrp_o22 = intrp.RectBivariateSpline(
tstar22, delta, omega22_table, kx=5, ky=5
)
intrp_Astar = intrp.RectBivariateSpline(
tstar, delta, astar_table, kx=5, ky=5
)

# Get molecular mass
MW = completeSpecies("MW", usersp, refsp)
Expand Down Expand Up @@ -135,7 +145,9 @@ def cp_R(cp0, poly, T, MW):
Tstar = T * kb / well
omga22 = intrp_o22(Tstar, r_deltastar.diagonal(), grid=False)
visc[i, :] = (
(5.0 / 16.0) * np.sqrt(np.pi * mass * kb * T) / (np.pi * diam**2 * omga22)
(5.0 / 16.0)
* np.sqrt(np.pi * mass * kb * T)
/ (np.pi * diam**2 * omga22)
)

##########################################
Expand Down Expand Up @@ -221,13 +233,16 @@ def cp_R(cp0, poly, T, MW):
logTs = np.log(Ts)
sqrtTs = np.sqrt(Ts)

# We fit the visc pol'y to the sqrtT as visc is proportional to sqrtT
# We fit the sqrt(sqrt(T)/visc) pol'y to the log(T) as visc is proportional to sqrtT
# we also reverse the numpy poly'l so lowest order is first
visc = visc / sqrtTs[:, None]
visc = np.sqrt(visc / sqrtTs[:, None])
w = 1.0 / (visc**2)
muPoly = np.flip(
np.array(
[list(np.polyfit(logTs, visc[:, k], deg=deg, w=w[:, k])) for k in range(ns)]
[
list(np.polyfit(logTs, visc[:, k], deg=deg, w=w[:, k]))
for k in range(ns)
]
),
-1,
)
Expand All @@ -238,7 +253,10 @@ def cp_R(cp0, poly, T, MW):
w = 1.0 / (cond**2)
kappaPoly = np.flip(
np.array(
[list(np.polyfit(logTs, cond[:, k], deg=deg, w=w[:, k])) for k in range(ns)]
[
list(np.polyfit(logTs, cond[:, k], deg=deg, w=w[:, k]))
for k in range(ns)
]
),
-1,
)
Expand All @@ -248,7 +266,9 @@ def cp_R(cp0, poly, T, MW):
w = 1.0 / (diff**2)
for k in range(ns):
for j in range(k, ns):
Dij.append(list(np.polyfit(logTs, diff[:, k, j], deg=deg, w=w[:, k, j])))
Dij.append(
list(np.polyfit(logTs, diff[:, k, j], deg=deg, w=w[:, k, j]))
)

DijPoly = np.flip(np.array(Dij), -1)

Expand Down
18 changes: 15 additions & 3 deletions tests/transport/test_kineticTheory.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,12 @@
##############################################
# Test kinetic theory
##############################################
# NOTE: These ranges vary by input file in Cantera. It seems to set the
# minTemp and maxTemp based on the min/max ranges of the NASA7 poly'l
# data. In testing, this seems to explain the errors we sometimes get
# in thermodynamic testing against Cantera.
# (i.e. takes error from 1% to 0.001%). For now we just use sensible values
# here.

pytestmark = pytest.mark.parametrize(
"ctfile,thfile",
Expand All @@ -30,7 +36,9 @@

def test_kineticTheory(my_setup, ctfile, thfile):
relpath = str(Path(__file__).parent)
ct.add_directory(relpath + "/../../src/peregrinepy/thermoTransport/database/source")
ct.add_directory(
relpath + "/../../src/peregrinepy/thermoTransport/database/source"
)

gas = ct.Solution(ctfile)
p = np.random.uniform(low=10000, high=1000000)
Expand Down Expand Up @@ -85,12 +93,16 @@ def print_diff(name, c, p):
pd.append(print_diff("T", gas.T, pgprim[4]))
for i, n in enumerate(gas.species_names[0:-1]):
pd.append(print_diff(n, gas.Y[i], pgprim[5 + i]))
pd.append(print_diff(gas.species_names[-1], gas.Y[-1], 1.0 - np.sum(pgprim[5::])))
pd.append(
print_diff(gas.species_names[-1], gas.Y[-1], 1.0 - np.sum(pgprim[5::]))
)
print("Mixture Properties")
pd.append(print_diff("mu", gas.viscosity, pgtrns[0]))
pd.append(print_diff("kappa", gas.thermal_conductivity, pgtrns[1]))
for i, n in enumerate(gas.species_names):
pd.append(print_diff(f"D_{n}", gas.mix_diff_coeffs_mass[i], pgtrns[2 + i]))
pd.append(
print_diff(f"D_{n}", gas.mix_diff_coeffs_mass[i], pgtrns[2 + i])
)

passfail = np.all(np.array(pd) < 1.0)
assert passfail

0 comments on commit 4f71dbf

Please sign in to comment.