diff --git a/src/transport/include/antioch/mixture_averaged_transport_evaluator.h b/src/transport/include/antioch/mixture_averaged_transport_evaluator.h index bd624d19..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,28 +357,36 @@ 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); + StateType k_sum2 = zero_clone(T); D_vec = zero_clone(mass_fractions); 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 ) @@ -371,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 ); @@ -380,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, @@ -398,11 +421,15 @@ namespace Antioch mu[s] ); } - k_mix += k[s]*chi[s]/phi_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; } - + k_mix = .5*(k_sum1+1/k_sum2); if( DiffusionTraits::is_species_diffusion ) @@ -579,9 +606,11 @@ 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); + 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); @@ -599,15 +628,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"); @@ -615,6 +645,41 @@ namespace Antioch } // 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 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;