diff --git a/data/example_data b/data/example_data index 1a5d27e508..dbc42129c9 160000 --- a/data/example_data +++ b/data/example_data @@ -1 +1 @@ -Subproject commit 1a5d27e508a38b1791543e9fded80ffd5c5b8d75 +Subproject commit dbc42129c998c3a977dfa2a8b779851776fb2982 diff --git a/include/cantera/transport/HighPressureGasTransport.h b/include/cantera/transport/HighPressureGasTransport.h index 202a18adcd..e158b1480f 100644 --- a/include/cantera/transport/HighPressureGasTransport.h +++ b/include/cantera/transport/HighPressureGasTransport.h @@ -11,7 +11,6 @@ // Cantera includes #include "GasTransport.h" -#include "cantera/numerics/DenseMatrix.h" #include "cantera/transport/MixTransport.h" namespace Cantera @@ -45,7 +44,7 @@ class HighPressureGasTransport : public MixTransport public: string transportModel() const override { - return "HighPressureGas"; + return "high-pressure"; } double thermalConductivity() override; @@ -148,7 +147,7 @@ class ChungHighPressureGasTransport : public MixTransport public: string transportModel() const override { - return "ChungHighPressureGas"; + return "high-pressure-chung"; } double viscosity() override; diff --git a/src/transport/HighPressureGasTransport.cpp b/src/transport/HighPressureGasTransport.cpp index 8841e03a56..96992c31d4 100644 --- a/src/transport/HighPressureGasTransport.cpp +++ b/src/transport/HighPressureGasTransport.cpp @@ -7,7 +7,6 @@ // at https://cantera.org/license.txt for license and copyright information. #include "cantera/transport/HighPressureGasTransport.h" -#include "cantera/numerics/ctlapack.h" #include "cantera/base/utilities.h" #include "cantera/thermo/IdealGasPhase.h" #include "cantera/transport/TransportFactory.h" @@ -401,14 +400,6 @@ double HighPressureGasTransport::FQ_i(double Q, double Tr, double MW) *sign(Tr - 12.0)); } -// Sign function for use in the quantum correction term calculation. -template -int sign(T val) { - if (val > 0) return 1; - if (val < 0) return -1; - return 0; -} - // Calculates the pressure correction factor for the Lucas method that // is used to calculate high pressure pure fluid viscosity. This is equation // 9-4.18 in Poling et al. (2001): @@ -425,6 +416,12 @@ double HighPressureGasTransport::FP_i(double mu_r, double Tr, double Z_crit) double HighPressureGasTransport::compute_correction_factor(double Pr, double Tr) { + // In the low pressure limit, no correction is needed. Interpolate + // the value towards 1 as pressure drops below the 0.1 threshold. + if (Pr < 0.1) { + return 1.0; + } + // Data from Table 2 of Takahashi 1975 paper: const static double Pr_lookup[17] = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0, 4.0, 5.0}; @@ -442,25 +439,22 @@ double HighPressureGasTransport::compute_correction_factor(double Pr, double Tr) 14., 10.00900, 8.57519, 10.37483, 11.21674, 1., 6.19043, 1., 1.}; // Interpolate to obtain the value of the constants (DP)_R, A, B, C, E at - // the provided value of the reduced pressure (Pr): + // the provided value of the reduced pressure (Pr). int Pr_i = 0; double frac = 0.0; double A, B, C,E, DP_Rt; - if (Pr < 0.1) { // In the low pressure limit, no correction is needed: - frac = (Pr - Pr_lookup[0])/(Pr_lookup[1] - Pr_lookup[0]); - } else { - for (int j = 1; j < 17; j++) { - if (Pr_lookup[j] > Pr) { - frac = (Pr - Pr_lookup[j-1])/(Pr_lookup[j] - Pr_lookup[j-1]); - break; - } - Pr_i++; + for (int j = 1; j < 17; j++){ + if (Pr_lookup[j] > Pr) { + frac = (Pr - Pr_lookup[j-1])/(Pr_lookup[j] - Pr_lookup[j-1]); + break; } - // If this loop completes without finding a bounding value of Pr, use - // the final table value. - frac = 1.0; + Pr_i++; } + // If this loop completes without finding a bounding value of Pr, use + // the final table value. + frac = 1.0; + DP_Rt = DP_Rt_lookup[Pr_i]*(1.0 - frac) + DP_Rt_lookup[Pr_i+1]*frac; A = A_ij_lookup[Pr_i]*(1.0 - frac) + A_ij_lookup[Pr_i+1]*frac; @@ -478,16 +472,25 @@ double HighPressureGasTransport::compute_correction_factor(double Pr, double Tr) namespace Cantera { - +// This implements the high-pressure gas mixture viscosity model of Chung +// described in chapter 9-7 of Poling et al. (2001). double ChungHighPressureGasTransport::viscosity() { + ChungMixtureParameters params; + compute_mixture_parameters(params); + + // Compute T_star using equation 9-5.26, using the mixture parameters + double tKelvin = m_thermo->temperature(); + double T_star = tKelvin / params.epsilon_over_k_mix; + + double density = m_thermo->density(); // The density is required for high-pressure gases - // Vc needs to be in units of cm^3/mol - // Starting on page 9-7 in the "Method of Chung, et. al. (1984, 1988)" section - double epsilon_over_k = Tc / 1.2593; + // This result is in units of micropoise + double viscosity = high_pressure_viscosity(T_star, params.MW_mix, density, params.Vc_mix, params.Tc_mix, params.acentric_factor_mix, params.mu_r_mix, params.kappa_mix); - // Equation 9-4.1 - double T_star = T / epsilon_over_k; + double micropoise_to_pascals_second = 1e-7; + + return viscosity*micropoise_to_pascals_second; } @@ -577,23 +580,6 @@ void ChungHighPressureGasTransport::compute_mixture_parameters(ChungMixtureParam } -// This function is structured such that it can be used for pure species or mixtures, with the -// only difference being the values that are passed to the function (pure values versus mixture values). -double ChungHighPressureGasTransport::low_pressure_viscosity(double T, double T_star, double MW, double acentric_factor, double mu_r, double sigma, double kappa) -{ - // Equation 9-4.3 in Poling et al. (2001) - double omega = neufeld_collision_integral(T_star); - - // Molecular shapes and polarities factor, equation 9-4.11 - double Fc = 1 -0.2756*acentric_factor + 0.059035*pow(mu_r, 4.0) + kappa; - - // Equation 9-3.9 in Poling et al. (2001), multiplied by the Chung factor, Fc - // (another way of writing 9-4.10 that avoids explicit use of the critical volume in this method) - double viscosity = Fc* (26.69*sqrt(MW*T)/(sigma*sigma*omega)); - - return viscosity; -} - // Implementation of equation 9-4.3 in Poling et al. (2001). // Applicable over the range of 0.3 <= T_star <= 100. double neufeld_collision_integral(double T_star) @@ -610,6 +596,22 @@ double neufeld_collision_integral(double T_star) return omega; } +// This function is structured such that it can be used for pure species or mixtures, with the +// only difference being the values that are passed to the function (pure values versus mixture values). +double ChungHighPressureGasTransport::low_pressure_viscosity(double T, double T_star, double MW, double acentric_factor, double mu_r, double sigma, double kappa) +{ + // Equation 9-4.3 in Poling et al. (2001) + double omega = neufeld_collision_integral(T_star); + + // Molecular shapes and polarities factor, equation 9-4.11 + double Fc = 1 -0.2756*acentric_factor + 0.059035*pow(mu_r, 4.0) + kappa; + + // Equation 9-3.9 in Poling et al. (2001), multiplied by the Chung factor, Fc + // (another way of writing 9-4.10 that avoids explicit use of the critical volume in this method) + double viscosity = Fc* (26.69*sqrt(MW*T)/(sigma*sigma*omega)); + + return viscosity; +} // Computes the high-pressure viscosity using the Chung method (Equation 9-6.18). // This function is structured such that it can be used for pure species or mixtures, with the @@ -658,7 +660,7 @@ double ChungHighPressureGasTransport::high_pressure_viscosity(double T_star, dou double eta = eta_1 * 36.344 * sqrt(MW*Tc) / pow(Vc, 2.0/3.0); // Equation 9-6.18 - return eta + return eta; } diff --git a/src/transport/TransportFactory.cpp b/src/transport/TransportFactory.cpp index f721ca08e0..ceef2150dc 100644 --- a/src/transport/TransportFactory.cpp +++ b/src/transport/TransportFactory.cpp @@ -46,6 +46,7 @@ TransportFactory::TransportFactory() addDeprecatedAlias("water", "Water"); reg("high-pressure", []() { return new HighPressureGasTransport(); }); addDeprecatedAlias("high-pressure", "HighP"); + reg("high-pressure-chung", []() { return new ChungHighPressureGasTransport(); }); m_CK_mode["CK_Mix"] = m_CK_mode["mixture-averaged-CK"] = true; m_CK_mode["CK_Multi"] = m_CK_mode["multicomponent-CK"] = true; } diff --git a/test/python/test_transport.py b/test/python/test_transport.py index 513958c631..c1ea726b1d 100644 --- a/test/python/test_transport.py +++ b/test/python/test_transport.py @@ -147,6 +147,19 @@ def test_add_species_multi(self, cantera_data_path): assert gas1.thermal_conductivity == approx(gas2.thermal_conductivity) assert gas1.multi_diff_coeffs == approx(gas2.multi_diff_coeffs) + def test_high_pressure_chung_viscosity(self): + state = 520, 6e7, 'NH3:1.0' + + gas1 = ct.Solution('gri30.yaml') + gas1.transport_model = 'high-pressure-chung' + gas1.TPX = state + + # Values from Poling et al. (2001), example 9-12 for ammonia + experimental_viscosity = 466e-7 # 466 micropoise = 466e-7 Pa s + + error = abs(gas1.viscosity - experimental_viscosity) / experimental_viscosity + assert error <= 0.05 # Error should be less than 5% + def test_species_visosities(self, phase): for species_name in phase.species_names: # check that species viscosity matches overall for single-species