GasTransport.h

Go to the documentation of this file.

/**
 * @file GasTransport.h
 */
 
// This file is part of Cantera. See License.txt in the top-level directory or
// at https://cantera.org/license.txt for license and copyright information.
 
#ifndef CT_GAS_TRANSPORT_H
#define CT_GAS_TRANSPORT_H
 
#include "Transport.h"
#include "cantera/numerics/DenseMatrix.h"
 
namespace Cantera
{
 
class MMCollisionInt;
 
//! Class GasTransport implements some functions and properties that are
//! shared by the MixTransport and MultiTransport classes.
//!
//! For details, see Kee, et al. @cite kee2003 and @cite kee2017.
//!
//! @ingroup tranprops
class GasTransport : public Transport
{
public:
    //! Get the viscosity [Pa·s] of the mixture
    /*!
     * The viscosity is computed using the Wilke mixture rule, as given in Poling et al.
     * @cite poling2001, Eqs. (9-5.13 and 9-5.14):
     * @f[
     *     \eta = \sum_k \frac{\eta_k X_k}{\sum_j \Phi_{k,j} X_j}.
     * @f]
     *
     * Here, @f$ \eta_k @f$ is the viscosity of pure species *k* and the weighting
     * function @f$ \Phi_{k,j} @f$ is:
     * @f[
     *     \Phi_{k,j} = \frac{ \left[ 1 + \left( \eta_k / \eta_j \right)^{1/2}
     *                               \left( M_j / M_k \right)^{1/4} \right]^2 }
     *                      {\left[ 8 \left( 1 + M_k / M_j \right) \right]^{1/2}}
     * @f]
     * @see updateViscosity_T()
     */
    double viscosity() override;
 
    //! Get the pure-species viscosities [Pa·s]
    void getSpeciesViscosities(double* const visc) override {
        update_T();
        updateViscosity_T();
        std::copy(m_visc.begin(), m_visc.end(), visc);
    }
 
    void getBinaryDiffCoeffs(const size_t ld, double* const d) override;
 
    //! Returns the Mixture-averaged diffusion coefficients [m²/s].
    /*!
     * Returns the mixture averaged diffusion coefficients for a gas,
     * appropriate for calculating the mass averaged diffusive flux with respect
     * to the mass averaged velocity using gradients of the mole fraction.
     * Note, for the single species case or the pure fluid case the routine
     * returns the self-diffusion coefficient. This is needed to avoid a Nan
     * result in the formula below.
     *
     * This is Eqn. 12.180 from "Chemically Reacting Flow"
     *
     * @f[
     *     D_{km}' = \frac{\left( \bar{M} - X_k M_k \right)}{ \bar{\qquad M \qquad } }
     *               {\left( \sum_{j \ne k} \frac{X_j}{\mathcal{D}_{kj}} \right) }^{-1}
     * @f]
     *
     * @param[out] d  Vector of mixture diffusion coefficients, @f$ D_{km}' @f$ ,
     *     for each species; length is the number of species.
     */
    void getMixDiffCoeffs(double* const d) override;
 
    //! Returns the mixture-averaged diffusion coefficients [m²/s].
    //! These are the coefficients for calculating the molar diffusive fluxes
    //! from the species mole fraction gradients, computed according to
    //! Eq. 12.176 in "Chemically Reacting Flow":
    //!
    //! @f[  D_{km}^* = \frac{1-X_k}{\sum_{j \ne k}^K X_j/\mathcal{D}_{kj}} @f]
    //!
    //! @param[out] d vector of mixture-averaged diffusion coefficients for
    //!     each species, length #m_nsp.
    void getMixDiffCoeffsMole(double* const d) override;
 
    //! Returns the mixture-averaged diffusion coefficients [m²/s].
    /*!
     * These are the coefficients for calculating the diffusive mass fluxes
     * from the species mass fraction gradients, computed according to
     * Eq. 12.178 in "Chemically Reacting Flow":
     *
     * @f[
     *     \frac{1}{D_{km}} = \sum_{j \ne k}^K \frac{X_j}{\mathcal{D}_{kj}} +
     *     \frac{X_k}{1-Y_k} \sum_{j \ne k}^K \frac{Y_j}{\mathcal{D}_{kj}}
     * @f]
     *
     * @param[out] d vector of mixture-averaged diffusion coefficients for
     *     each species, length #m_nsp.
     */
    void getMixDiffCoeffsMass(double* const d) override;
 
    //! Return the polynomial fits to the viscosity of species `i`.
    //! @see fitProperties()
    void getViscosityPolynomial(size_t i, double* coeffs) const override;
 
    //! Return the temperature fits of the heat conductivity of species `i`.
    //! @see fitProperties()
    void getConductivityPolynomial(size_t i, double* coeffs) const override;
 
    //! Return the polynomial fits to the binary diffusivity of species pair (i, j)
    //! @see fitDiffCoeffs()
    void getBinDiffusivityPolynomial(size_t i, size_t j, double* coeffs) const override;
 
    //! Return the polynomial fits to the collision integral of species pair (i, j)
    //! @see fitCollisionIntegrals()
    void getCollisionIntegralPolynomial(size_t i, size_t j,
                                        double* astar_coeffs,
                                        double* bstar_coeffs,
                                        double* cstar_coeffs) const override;
 
    //! Modify the polynomial fits to the viscosity of species `i`
    //! @see fitProperties()
    void setViscosityPolynomial(size_t i, double* coeffs) override;
 
    //! Modify the temperature fits of the heat conductivity of species `i`
    //! @see fitProperties()
    void setConductivityPolynomial(size_t i, double* coeffs) override;
 
    //! Modify the polynomial fits to the binary diffusivity of species pair (i, j)
    //! @see fitDiffCoeffs()
    void setBinDiffusivityPolynomial(size_t i, size_t j, double* coeffs) override;
 
    //! Modify the polynomial fits to the collision integral of species pair (i, j)
    //! @see fitCollisionIntegrals()
    void setCollisionIntegralPolynomial(size_t i, size_t j,
                                        double* astar_coeffs,
                                        double* bstar_coeffs,
                                        double* cstar_coeffs, bool actualT) override;
 
    void init(ThermoPhase* thermo, int mode=0) override;
 
    bool CKMode() const override {
        return m_mode == CK_Mode;
    }
 
    void invalidateCache() override;
 
protected:
    GasTransport();
 
    virtual void update_T();
    virtual void update_C() = 0;
 
    //! Update the temperature-dependent viscosity terms.
    /**
     * Updates the array of pure species viscosities, and the weighting
     * functions in the viscosity mixture rule. The flag #m_visc_ok is set to true.
     * @see viscosity()
     */
    virtual void updateViscosity_T();
 
    //! Update the pure-species viscosities. These are evaluated from the
    //! polynomial fits of the temperature and are assumed to be independent
    //! of pressure.
    virtual void updateSpeciesViscosities();
 
    //! Update the binary diffusion coefficients
    /*!
     * These are evaluated from the polynomial fits of the temperature at the
     * unit pressure of 1 Pa.
     */
    virtual void updateDiff_T();
 
    //! @name Initialization
    //! @{
 
    //! Setup parameters for a new kinetic-theory-based transport manager for
    //! low-density gases
    virtual void setupCollisionParameters();
 
    //! Setup range for polynomial fits to collision integrals of
    //! Monchick & Mason @cite monchick1961
    void setupCollisionIntegral();
 
    //! Read the transport database
    /*!
     * Read transport property data from a file for a list of species. Given the
     * name of a file containing transport property parameters and a list of
     * species names.
     */
    virtual void getTransportData();
 
    //! Corrections for polar-nonpolar binary diffusion coefficients
    /*!
     * Calculate corrections to the well depth parameter and the diameter for use in
     * computing the binary diffusion coefficient of polar-nonpolar pairs. For more
     * information about this correction, see Dixon-Lewis @cite dixon-lewis1968.
     *
     * @param i        Species one - this is a bimolecular correction routine
     * @param j        species two - this is a bimolecular correction routine
     * @param f_eps    Multiplicative correction factor to be applied to epsilon(i,j)
     * @param f_sigma  Multiplicative correction factor to be applied to diam(i,j)
     */
    void makePolarCorrections(size_t i, size_t j, double& f_eps,
                              double& f_sigma);
 
    //! Generate polynomial fits to collision integrals
    /*!
     * @param integrals interpolator for the collision integrals
     */
    void fitCollisionIntegrals(MMCollisionInt& integrals);
 
    //! Generate polynomial fits to the viscosity @f$ \eta @f$ and conductivity
    //! @f$ \lambda @f$.
    /*!
     * If CKMode(), then the fits are of the form
     * @f[
     *      \ln \eta(i) = \sum_{n=0}^3 a_n(i) \, (\ln T)^n
     * @f]
     * and
     * @f[
     *      \ln \lambda(i) = \sum_{n=0}^3 b_n(i) \, (\ln T)^n
     * @f]
     * Otherwise the fits are of the form
     * @f[
     *      \left(\eta(i)\right)^{1/2} = T^{1/4} \sum_{n=0}^4 a_n(i) \, (\ln T)^n
     * @f]
     * and
     * @f[
     *      \lambda(i) = T^{1/2} \sum_{n=0}^4 b_n(i) \, (\ln T)^n
     * @f]
     *
     * @param integrals interpolator for the collision integrals
     */
    virtual void fitProperties(MMCollisionInt& integrals);
 
    //! Generate polynomial fits to the binary diffusion coefficients
    /*!
     * If CKMode(), then the fits are of the form
     * @f[
     *      \ln \mathcal{D}(i,j) = \frac{1}{p} \sum_{n=0}^3 c_n(i,j) \, (\ln T)^n
     * @f]
     * Otherwise the fits are of the form
     * @f[
     *      \mathcal{D}(i,j) = \frac{T^{3/2}}{p} \sum_{n=0}^4 c_n(i,j) \, (\ln T)^n
     * @f]
     *
     * @param integrals interpolator for the collision integrals
     */
    virtual void fitDiffCoeffs(MMCollisionInt& integrals);
 
    //! Second-order correction to the binary diffusion coefficients
    /*!
     * Calculate second-order corrections to binary diffusion coefficient pair
     * (dkj, djk). At first order, the binary diffusion coefficients are
     * independent of composition, and d(k,j) = d(j,k). But at second order,
     * there is a weak dependence on composition, with the result that d(k,j) !=
     * d(j,k). This method computes the multiplier by which the first-order
     * binary diffusion coefficient should be multiplied to produce the value
     * correct to second order. The expressions here are taken from Marerro and
     * Mason @cite marrero1972.
     *
     * @param t   Temperature (K)
     * @param integrals interpolator for the collision integrals
     * @param k   index of first species
     * @param j   index of second species
     * @param xk  Mole fraction of species k
     * @param xj  Mole fraction of species j
     * @param fkj multiplier for d(k,j)
     * @param fjk multiplier for d(j,k)
     *
     * @note This method is not used currently.
     */
    void getBinDiffCorrection(double t, MMCollisionInt& integrals, size_t k,
                              size_t j, double xk, double xj,
                              double& fkj, double& fjk);
 
    //! @}
 
    //! Vector of species mole fractions. These are processed so that all mole
    //! fractions are >= #Tiny. Length = #m_nsp.
    vector<double> m_molefracs;
 
    //! Internal storage for the viscosity of the mixture [Pa·s]
    double m_viscmix = 0.0;
 
    //! Update boolean for mixture rule for the mixture viscosity
    bool m_visc_ok = false;
 
    //! Update boolean for the weighting factors for the mixture viscosity
    bool m_viscwt_ok = false;
 
    //! Update boolean for the species viscosities
    bool m_spvisc_ok = false;
 
    //! Update boolean for the binary diffusivities at unit pressure
    bool m_bindiff_ok = false;
 
    //! Type of the polynomial fits to temperature. `CK_Mode` means Chemkin mode.
    //! Any other value means to use %Cantera's preferred fitting functions.
    int m_mode = 0;
 
    //! Viscosity weighting function. size = #m_nsp * #m_nsp
    DenseMatrix m_phi;
 
    //! work space length = #m_nsp
    vector<double> m_spwork;
 
    //! vector of species viscosities [Pa·s]. These are used in Wilke's
    //! rule to calculate the viscosity of the solution. length = #m_nsp.
    vector<double> m_visc;
 
    //! Polynomial fits to the viscosity of each species. `m_visccoeffs[k]` is
    //! the vector of polynomial coefficients for species `k` that fits the
    //! viscosity as a function of temperature.
    vector<vector<double>> m_visccoeffs;
 
    //! Local copy of the species molecular weights.
    vector<double> m_mw;
 
    //! Holds square roots of molecular weight ratios
    /*!
     *  @code
     *  m_wratjk(j,k)  = sqrt(mw[j]/mw[k])        j < k
     *  m_wratjk(k,j)  = sqrt(sqrt(mw[j]/mw[k]))  j < k
     *  @endcode
     */
    DenseMatrix m_wratjk;
 
    //! Holds square roots of molecular weight ratios
    /*!
     *  `m_wratjk1(j,k)  = sqrt(1.0 + mw[k]/mw[j])        j < k`
     */
    DenseMatrix m_wratkj1;
 
    //! vector of square root of species viscosities. These are used in Wilke's rule to
    //! calculate the viscosity of the solution. length = #m_nsp.
    vector<double> m_sqvisc;
 
    //! Powers of the ln temperature, up to fourth order
    vector<double> m_polytempvec;
 
    //! Current value of the temperature [K] at which the properties in this object
    //! are calculated.
    double m_temp = -1.0;
 
    //! Current value of Boltzmann constant times the temperature [J]
    double m_kbt = 0.0;
 
    //! current value of temperature to 1/2 power
    double m_sqrt_t = 0.0;
 
    //! Current value of the log of the temperature
    double m_logt = 0.0;
 
    //! Current value of temperature to 1/4 power
    double m_t14 = 0.0;
 
    //! Polynomial fits to the binary diffusivity of each species
    /*!
     * `m_diffcoeff[ic]` is the vector of polynomial coefficients that fits the binary
     * diffusion coefficient for species `i` and `j`. The relationship between `i`, `j`,
     * and `ic` is determined from the following algorithm:
     * @code
     * int ic = 0;
     * for (i = 0; i < m_nsp; i++) {
     *    for (j = i; j < m_nsp; j++) {
     *      ic++;
     *    }
     * }
     * @endcode
     */
    vector<vector<double>> m_diffcoeffs;
 
    //! Matrix of binary diffusion coefficients at the reference pressure and
    //! the current temperature Size is #m_nsp x #m_nsp.
    DenseMatrix m_bdiff;
 
    //! temperature fits of the heat conduction
    /*!
     *  Dimensions are number of species (#m_nsp) and polynomial order of the collision
     *  integral fit (degree+1).
     */
    vector<vector<double>> m_condcoeffs;
 
    //! Indices for the (i,j) interaction in collision integral fits
    /*!
     *  `m_poly[i][j]` contains the index for (i,j) interactions in
     *  #m_omega22_poly, #m_astar_poly, #m_bstar_poly, and #m_cstar_poly.
     */
    vector<vector<int>> m_poly;
 
    //! Fit for omega22 collision integral
    /*!
     * `m_omega22_poly[m_poly[i][j]]` is the vector of polynomial coefficients
     * (length degree+1) for the collision integral fit for the species pair (i,j).
     */
    vector<vector<double>> m_omega22_poly;
 
    //! Flag to indicate for which (i,j) interaction pairs the
    //! actual temperature is used instead of the reduced temperature
    vector<vector<int>> m_star_poly_uses_actualT;
 
    //! Fit for astar collision integral
    /*!
     * `m_astar_poly[m_poly[i][j]]` is the vector of polynomial coefficients
     * (length degree+1) for the collision integral fit for the species pair (i,j).
     */
    vector<vector<double>> m_astar_poly;
 
    //! Fit for bstar collision integral
    /*!
     * `m_bstar_poly[m_poly[i][j]]` is the vector of polynomial coefficients
     * (length degree+1) for the collision integral fit for the species pair (i,j).
     */
    vector<vector<double>> m_bstar_poly;
 
    //! Fit for cstar collision integral
    /*!
     * `m_bstar_poly[m_poly[i][j]]` is the vector of polynomial coefficients
     * (length degree+1) for the collision integral fit for the species pair (i,j).
     */
    vector<vector<double>> m_cstar_poly;
 
    //! Rotational relaxation number for each species
    /*!
     * length is the number of species in the phase. units are dimensionless
     */
    vector<double> m_zrot;
 
    //! Dimensionless rotational heat capacity of each species
    /*!
     * These values are 0, 1 and 1.5 for single-molecule, linear, and nonlinear
     * species respectively. length is the number of species in the phase.
     * Dimensionless  (Cr / R)
     */
    vector<double> m_crot;
 
    //! Vector of booleans indicating whether a species is a polar molecule
    /*!
     * Length is #m_nsp.
     */
    vector<bool> m_polar;
 
    //! Polarizability [m³] of each species in the phase
    vector<double> m_alpha;
 
    //! Lennard-Jones well-depth [J] of the species in the current phase
    /*!
     * length is the number of species in the phase. Note this is not J/kmol; this is a
     * per molecule amount.
     */
    vector<double> m_eps;
 
    //! Lennard-Jones diameter [m] of the species in the current phase
    vector<double> m_sigma;
 
    //! This is the reduced mass [kg] of the interaction between species i and j
    /*!
     * @f[
     *     m_{ij} = \frac{M_i M_j}{N_A (M_i + M_j)}
     * @f]
     *  Length #m_nsp * #m_nsp. This is a symmetric matrix
     */
    DenseMatrix m_reducedMass;
 
    //! hard-sphere diameter [m] for (i,j) collision
    /*!
     * @f[
     *     \sigma_{ij} = \frac{\sigma_i + \sigma_j}{2}
     * @f]
     *  Length #m_nsp * #m_nsp. This is a symmetric matrix.
     */
    DenseMatrix m_diam;
 
    //! The effective well depth [J] for (i,j) collisions
    /*!
     * @f[
     *     \epsilon_{ij} = \sqrt{\epsilon_i \epsilon_j}
     * @f]
     * Length #m_nsp * #m_nsp. This is a symmetric matrix.
     */
    DenseMatrix m_epsilon;
 
    //! The effective dipole moment [Coulomb·m] for (i,j) collisions
    /*!
     * @f[
     *    \mu_{ij} = \sqrt{\mu_i \mu_j}
     * @f]
     *
     * Dipole moments are conventionally given in Debye. The conversion factor to
     * Coulomb·m is @f$ 10^{-21} / c \approx 3.335 \times 10^{-30} @f$.
     *
     * Length #m_nsp * #m_nsp. This is a symmetric matrix.
     */
    DenseMatrix m_dipole;
 
    //! Reduced dipole moment of the interaction between two species
    /*!
     * @f[
     *     \tilde{\delta}^*_{ij} =
     *         \frac{ \mu_{ij}^2 }{ 8 \pi \varepsilon_0 \epsilon_{ij} \sigma_{ij}^3 }
     * @f]
     *  Length #m_nsp * #m_nsp. This is a symmetric matrix
     */
    DenseMatrix m_delta;
 
    //! Pitzer acentric factor [dimensionless]
    /*!
     * Length is the number of species in the phase.
     */
    vector<double> m_w_ac;
 
    //! Dispersion coefficient normalized by the square of the elementary charge [m⁵].
    vector<double> m_disp;
 
    //! Quadrupole polarizability
    vector<double> m_quad_polar;
};
 
} // namespace Cantera
 
#endif