From 49eff2ee373b9e6c22f2c7c2248ce7b966f5dceb Mon Sep 17 00:00:00 2001 From: Kenneth Budzinski Date: Mon, 20 Aug 2018 10:04:25 -0400 Subject: [PATCH 1/4] Changing how we calculate the mixture averaged diffusion coefficients by using mole fractions instead of mass fractions, like cantera, which keeps the mixture molecular weight more bearable. As well as clipping mole fractions that were negative or over 1. --- .../mixture_averaged_transport_evaluator.h | 26 ++++++++++++++----- 1 file changed, 19 insertions(+), 7 deletions(-) diff --git a/src/transport/include/antioch/mixture_averaged_transport_evaluator.h b/src/transport/include/antioch/mixture_averaged_transport_evaluator.h index bd624d19..73b94670 100644 --- a/src/transport/include/antioch/mixture_averaged_transport_evaluator.h +++ b/src/transport/include/antioch/mixture_averaged_transport_evaluator.h @@ -579,12 +579,23 @@ namespace Antioch } case(MASS_FLUX_MASS_FRACTION): { - VectorStateType molar_fractions = zero_clone(mass_fractions); - - mixture.X(mixture.M(mass_fractions),mass_fractions,molar_fractions); - + +VectorStateType molar_fractions = zero_clone(mass_fractions); + mixture.X(mixture.M(mass_fractions),mass_fractions,molar_fractions); + + //Clipping unrealistic values, and calculating MW_Mixture through the new mole fractions + typename value_type::type MW_Mixture = zero_clone(mass_fractions[0]); + for(unsigned int s=0; s < D_vec.size(); s++) + { + if(molar_fractions[s] > 1) + molar_fractions[s] = 1; + if(molar_fractions[s] < 1e-16) + molar_fractions[s] = 1e-16; + MW_Mixture += molar_fractions[s]*mixture.M(s); + } + typename value_type::type one = constant_clone(mass_fractions[0],1); - + // term1 term2 // 1/D_s = (sum_{j\ne s} X_j/D_{s,j}) + X_s/(1-Y_s)\sum_{j\ne s} Y_j/D_{s,j} for(unsigned int s = 0; s < D_vec.size(); s++) @@ -599,15 +610,16 @@ namespace Antioch term1 += molar_fractions[j]/D_mat[s][j]; - term2 += mass_fractions[j]/D_mat[s][j]; + term2 += molar_fractions[j]*mixture.M(j)/D_mat[s][j]; } - term2 *= molar_fractions[s]/(one - mass_fractions[s]); + term2 *= molar_fractions[s]/(MW_Mixture - mixture.M(s)*molar_fractions[s]); D_vec[s] = one/(term1+term2); } break; } + default: { antioch_error_msg("ERROR: Invalid DiffusivityType in MixtureAveragedTransportEvaluator::diffusion_mixing_rule"); From 3ad5a4acb96f6c5fdb78ed0ef4011cadf77a6178 Mon Sep 17 00:00:00 2001 From: Kenneth Budzinski Date: Fri, 24 Aug 2018 14:57:42 -0400 Subject: [PATCH 2/4] Changing the way we calculate the mixture thermal conductivity from the species conductivities. The formula used here is taken from "Chemically Reacting Flow" by Kee, Coltrin, and Glarborg, page 519. --- .../include/antioch/mixture_averaged_transport_evaluator.h | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/transport/include/antioch/mixture_averaged_transport_evaluator.h b/src/transport/include/antioch/mixture_averaged_transport_evaluator.h index 73b94670..cc815ca3 100644 --- a/src/transport/include/antioch/mixture_averaged_transport_evaluator.h +++ b/src/transport/include/antioch/mixture_averaged_transport_evaluator.h @@ -344,6 +344,8 @@ namespace Antioch mu_mix = zero_clone(T); k_mix = zero_clone(T); + StateType k_sum1 = zero_clone(T); + StateType k_sum2 = zero_clone(T); D_vec = zero_clone(mass_fractions); VectorStateType mu = zero_clone(mass_fractions); @@ -398,11 +400,12 @@ namespace Antioch mu[s] ); } - k_mix += k[s]*chi[s]/phi_s; + k_sum1 += k[s]*chi[s]; + k_sum2 += chi[s]/k[s]; mu_mix += mu[s]*chi[s]/phi_s; } - + k_mix = .5*(k_sum1+1/k_sum2); if( DiffusionTraits::is_species_diffusion ) From 13f61791032bd769060434a28912b8b6ac90ffd5 Mon Sep 17 00:00:00 2001 From: Kenneth Budzinski Date: Fri, 24 Aug 2018 15:12:37 -0400 Subject: [PATCH 3/4] Fixing the unit test that was tripped when altering the calculation of the thermal conductivity --- test/wilke_transport_unit.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/wilke_transport_unit.C b/test/wilke_transport_unit.C index 2bfdfd87..578ac6d3 100644 --- a/test/wilke_transport_unit.C +++ b/test/wilke_transport_unit.C @@ -276,7 +276,7 @@ int tester() */ const Scalar mu_kt_long_double = 4.49877527305932602332e-05L; - const Scalar k_kt_long_double = 8.22050332419571328635e-02L; + const Scalar k_kt_long_double = 8.27884568936580333168e-02L; std::vector D_kt_long_double(5,0); D_kt_long_double[0] = 1.95418749838889089562e-04L; From aedef6f58e6efca8289c70c024134e3ff920d9fa Mon Sep 17 00:00:00 2001 From: Kenneth Budzinski Date: Tue, 28 Aug 2018 14:40:21 -0400 Subject: [PATCH 4/4] Instead of clipping the mole fractions, implemented a method that will calculate perturbed mass fractions that can be used to calculate the transport, this will make it easier to implement a switch for perturbing the mass fractions. Also added a little more documentation, but can still use more here. --- .../mixture_averaged_transport_evaluator.h | 92 ++++++++++++++----- 1 file changed, 71 insertions(+), 21 deletions(-) diff --git a/src/transport/include/antioch/mixture_averaged_transport_evaluator.h b/src/transport/include/antioch/mixture_averaged_transport_evaluator.h index cc815ca3..77ea9535 100644 --- a/src/transport/include/antioch/mixture_averaged_transport_evaluator.h +++ b/src/transport/include/antioch/mixture_averaged_transport_evaluator.h @@ -43,6 +43,9 @@ #include "antioch/diffusion_traits.h" #include "antioch/conductivity_traits.h" +// C++ +#include + namespace Antioch { //! Compute transport properties using ``mixture averaged" model @@ -87,6 +90,7 @@ namespace Antioch MASS_FLUX_MASS_FRACTION, MASS_FLUX_MOLE_FRACTION }; + //! Mixture diffusivities for each species, in [m^2/s] /*! Only valid for BinaryDiffusionBase species diffusion models. * Compile time error if otherwise. @@ -166,6 +170,17 @@ namespace Antioch DiffusivityType diff_type, VectorStateType & D_vec ) const; + //! Compute Perturbe mole fractions + /*! To avoid singularities from dividing by very small numbers in the vanishing mass fraction regions, we use the EGLIB approach + * (http://blanche.polytechnique.fr/www.eglib/manual.ps page 5) + * We Use a perturbed mole fraction that ensures that the mole fractions are between 0 and 1. + */ + template + void perturbed_molar_fractions( const ChemicalMixture & mixture, + const VectorStateType & mass_fractions, + VectorStateType & perturbed_mass_fractions) const; + + const MixtureAveragedTransportMixture& _mixture; const MixtureDiffusion& _diffusion; @@ -342,6 +357,8 @@ namespace Antioch DiffusionTraits::is_species_diffusion), "Incompatible thermal conductivity and diffusion models!" ); + + mu_mix = zero_clone(T); k_mix = zero_clone(T); StateType k_sum1 = zero_clone(T); @@ -351,21 +368,25 @@ namespace Antioch VectorStateType mu = zero_clone(mass_fractions); VectorStateType k = zero_clone(mass_fractions); VectorStateType chi = zero_clone(mass_fractions); + VectorStateType perturbed_mass_fractions = zero_clone(mass_fractions); + this->perturbed_molar_fractions( _mixture.chem_mixture(), + mass_fractions, + perturbed_mass_fractions); typedef typename Antioch::rebind::type MatrixStateType; MatrixStateType mu_mu_sqrt(mu.size()); Antioch::init_constant(mu_mu_sqrt,mu); - MatrixStateType D_mat(mass_fractions.size()); + MatrixStateType D_mat(perturbed_mass_fractions.size()); init_constant(D_mat,D_vec); - this->compute_mu_chi( T, mass_fractions, mu, chi ); + this->compute_mu_chi( T, perturbed_mass_fractions, mu, chi ); this->compute_mu_mu_sqrt( mu, mu_mu_sqrt); - const StateType molar_density = rho / _mixture.chem_mixture().M(mass_fractions); // total molar density + const StateType molar_density = rho / _mixture.chem_mixture().M(perturbed_mass_fractions); // total molar density // If we're using a binary diffusion model, compute D_mat, D_vec now if( DiffusionTraits::is_binary_diffusion ) @@ -373,7 +394,7 @@ namespace Antioch _diffusion.compute_binary_diffusion_matrix(T, molar_density, D_mat); this->diffusion_mixing_rule( _mixture.chem_mixture(), - mass_fractions, + perturbed_mass_fractions, D_mat, diff_type, D_vec ); @@ -382,7 +403,7 @@ namespace Antioch for( unsigned int s = 0; s < _mixture.transport_mixture().n_species(); s++ ) { StateType phi_s = this->compute_phi( mu_mu_sqrt, chi, s ); - + // To calculate the mixture_averaged_conductivity, we need to calculate the pure species conductivity. if( ConductivityTraits::requires_diffusion ) { k[s] = _conductivity.conductivity_with_diffusion( s, @@ -400,6 +421,9 @@ namespace Antioch mu[s] ); } + /// To calculate the mixture averaged thermal conductivity and viscosity, we use the model described in: + /// Kee, Coltrin, and Glarborg's "Chemically Reacting Flow:Theory and Practice", John Wiley & Sons, 2003. + k_sum1 += k[s]*chi[s]; k_sum2 += chi[s]/k[s]; mu_mix += mu[s]*chi[s]/phi_s; @@ -582,23 +606,14 @@ namespace Antioch } case(MASS_FLUX_MASS_FRACTION): { - -VectorStateType molar_fractions = zero_clone(mass_fractions); - mixture.X(mixture.M(mass_fractions),mass_fractions,molar_fractions); - - //Clipping unrealistic values, and calculating MW_Mixture through the new mole fractions - typename value_type::type MW_Mixture = zero_clone(mass_fractions[0]); - for(unsigned int s=0; s < D_vec.size(); s++) - { - if(molar_fractions[s] > 1) - molar_fractions[s] = 1; - if(molar_fractions[s] < 1e-16) - molar_fractions[s] = 1e-16; - MW_Mixture += molar_fractions[s]*mixture.M(s); - } - + + VectorStateType molar_fractions = zero_clone(mass_fractions); + typename value_type::type MW_Mixture = zero_clone(mass_fractions[0]); + MW_Mixture = mixture.M(mass_fractions); + mixture.X(MW_Mixture,mass_fractions,molar_fractions); + typename value_type::type one = constant_clone(mass_fractions[0],1); - + // term1 term2 // 1/D_s = (sum_{j\ne s} X_j/D_{s,j}) + X_s/(1-Y_s)\sum_{j\ne s} Y_j/D_{s,j} for(unsigned int s = 0; s < D_vec.size(); s++) @@ -630,6 +645,41 @@ VectorStateType molar_fractions = zero_clone(mass_fractions); } // switch(diff_type) } + template + template + void MixtureAveragedTransportEvaluator::perturbed_molar_fractions( const ChemicalMixture & mixture, + const VectorStateType & mass_fractions, + VectorStateType & perturbed_mass_fractions) const + { + // convenient + typedef typename value_type::type StateType; + VectorStateType perturbed_molar_fractions = zero_clone(mass_fractions); + + mixture.X(mixture.M(mass_fractions),mass_fractions,perturbed_molar_fractions); + + // EGlib uses eps = 1e-16 + typename raw_value_type::type eps(std::numeric_limits::epsilon() * 10); + StateType mol_frac_sum = zero_clone(mass_fractions[0]); + + // (i) evaluate the mixture molecular weight over the total number of species + for(unsigned int s = 0; s < perturbed_molar_fractions.size(); s++) + mol_frac_sum += perturbed_molar_fractions[s]; + mol_frac_sum /= mixture.n_species(); + + // (ii) evaluate the perturbed mole fractions + // (iii) evaluate the perturbed mean molecular weight + StateType M_perturbed = zero_clone(perturbed_mass_fractions[0]); + for(unsigned int s =0; s < perturbed_molar_fractions.size(); s++) + { + perturbed_molar_fractions[s] += eps * (mol_frac_sum - perturbed_molar_fractions[s]); + M_perturbed += perturbed_molar_fractions[s] * mixture.M(s); + } + + // (iv) Calculate the perturbed mass_fractions + for(unsigned int s=0; s < perturbed_molar_fractions.size(); s++) + perturbed_mass_fractions[s] = perturbed_molar_fractions[s]*mixture.M(s)/M_perturbed; + } + } // end namespace Antioch #endif // ANTIOCH_WILKE_TRANSPORT_EVALUATOR_H