Branch data     Line data    Source code

       1                 :            : // *****************************************************************************
       2                 :            : /*!
       3                 :            :   \file      src/PDE/MultiSpecies/Mixture.cpp
       4                 :            :   \copyright 2012-2015 J. Bakosi,
       5                 :            :              2016-2018 Los Alamos National Security, LLC.,
       6                 :            :              2019-2021 Triad National Security, LLC.
       7                 :            :              All rights reserved. See the LICENSE file for details.
       8                 :            :   \brief     Multispecies mixture function
       9                 :            :   \details   This file declares functions for computing mixture flow quantities
      10                 :            : */
      11                 :            : // *****************************************************************************
      12                 :            : 
      13                 :            : #include "MultiSpecies/Mixture/Mixture.hpp"
      14                 :            : #include "MultiSpecies/MultiSpeciesIndexing.hpp"
      15                 :            : 
      16                 :            : using inciter::Mixture;
      17                 :            : 
      18                 :    7606896 : Mixture::Mixture(
      19                 :            :   std::size_t nspec,
      20                 :            :   const std::vector< tk::real >& ugp,
      21                 :    7606896 :   const std::vector< EOS >& mat_blk) :
      22                 :            :   m_nspec(nspec),
      23                 :            :   m_mix_density(0),
      24                 :            :   m_mix_R(0),
      25         [ +  - ]:    7606896 :   m_Ys(nspec, 0)
      26                 :            : // *************************************************************************
      27                 :            : //  Constructor (use during timestepping)
      28                 :            : //! \brief Initialize a mixture class using the state vector
      29                 :            : //! \param[in] nspec Number of species in mixture
      30                 :            : //! \param[in] ugp State vector
      31                 :            : //! \param[in] mat_blk EOS material block
      32                 :            : // *************************************************************************
      33                 :            : {
      34                 :            :   // Compute total density
      35                 :    7606896 :   m_mix_density = 0.;
      36         [ +  + ]:   19017240 :   for (std::size_t k=0; k<m_nspec; ++k)
      37                 :   11410344 :     m_mix_density += ugp[multispecies::densityIdx(m_nspec, k)];
      38                 :            : 
      39                 :            :   // Compute mass fractions
      40         [ +  + ]:   19017240 :   for (std::size_t k=0; k<m_nspec; ++k)
      41                 :   11410344 :     m_Ys[k] = ugp[multispecies::densityIdx(m_nspec, k)] / m_mix_density;
      42                 :            : 
      43                 :            :   // Compute mixture gas constant
      44                 :    7606896 :   m_mix_R = 0.;
      45         [ +  + ]:   19017240 :   for (std::size_t k = 0; k < m_nspec; k++)
      46         [ +  - ]:   11410344 :     m_mix_R += m_Ys[k] * mat_blk[k].compute< EOS::gas_constant >();
      47                 :    7606896 : }
      48                 :            : 
      49                 :     243294 : Mixture::Mixture(
      50                 :            :   std::size_t nspec,
      51                 :            :   const std::vector< tk::real >& Ys,
      52                 :            :   tk::real mix_pressure,
      53                 :            :   tk::real temperature,
      54                 :     243294 :   const std::vector< EOS >& mat_blk) :
      55                 :            :   m_nspec(nspec),
      56                 :            :   m_mix_density(0),
      57                 :            :   m_mix_R(0),
      58                 :     243294 :   m_Ys(Ys)
      59                 :            : // *************************************************************************
      60                 :            : //  Constructor (use during initialization)
      61                 :            : //! \brief Initialize a mixture class using the mixture thermodynamics and
      62                 :            : //!   known mass fractions.
      63                 :            : //! \param[in] nspec Number of species in mixture
      64                 :            : //! \param[in] Ys Mass fractions
      65                 :            : //! \param[in] mix_pressure Mixture pressure
      66                 :            : //! \param[in] temperature Temperature
      67                 :            : //! \param[in] mat_blk EOS material block
      68                 :            : // *************************************************************************
      69                 :            : {
      70                 :            :   // Compute mixture gas constant
      71                 :     243294 :   m_mix_R = 0.;
      72         [ +  + ]:     607110 :   for (std::size_t k = 0; k < m_nspec; k++)
      73         [ +  - ]:     363816 :     m_mix_R += m_Ys[k] * mat_blk[k].compute< EOS::gas_constant >();
      74                 :            : 
      75                 :            :   // Compute total density (via ideal gas EOS)
      76                 :     243294 :   m_mix_density = mix_pressure / (m_mix_R * temperature);
      77                 :     243294 : }
      78                 :            : 
      79                 :            : tk::real
      80                 :    3504600 : Mixture::frozen_soundspeed(
      81                 :            :   tk::real mix_density,
      82                 :            :   tk::real mix_temp,
      83                 :            :   const std::vector< EOS >& mat_blk) const
      84                 :            : // *************************************************************************
      85                 :            : //! \brief Calculate frozen speed of sound based on the mixture composition
      86                 :            : //!   and species parameters.
      87                 :            : //! \param[in] mix_density Mixture density (sum of species density)
      88                 :            : //! \param[in] mix_temp Mixture temperature (provided at call-site, since
      89                 :            : //!   it is reconstructed separately)
      90                 :            : //! \param[in] mat_blk EOS material block
      91                 :            : //! \return Mixture speed of sound using the ideal gas EoS
      92                 :            : // *************************************************************************
      93                 :            : {
      94                 :            :   // Clip pressure
      95                 :    3504600 :   auto mix_peff = std::max( 1.0e-15, pressure(mix_density, mix_temp) );
      96                 :            : 
      97                 :            :   // Compute beta, mixture parameters for sound speed calc.
      98                 :    3504600 :   tk::real mix_Cv = 0.;
      99         [ +  + ]:    8761500 :   for (std::size_t k = 0; k < m_nspec; k++) {
     100                 :    5256900 :     mix_Cv += mat_blk[k].compute< EOS::cv >(mix_temp) * m_Ys[k];
     101                 :            :   }
     102                 :    3504600 :   tk::real beta = m_mix_R / mix_Cv;
     103                 :            : 
     104                 :            :   // Compute speed of sound
     105                 :    3504600 :   tk::real a_sq = (1. + beta) * mix_peff / mix_density;
     106                 :    3504600 :   return std::sqrt(a_sq);
     107                 :            : }
     108                 :            : 
     109                 :            : tk::real
     110                 :    3654294 : Mixture::totalenergy(
     111                 :            :   tk::real mix_density,
     112                 :            :   tk::real u,
     113                 :            :   tk::real v,
     114                 :            :   tk::real w,
     115                 :            :   tk::real mix_temp,
     116                 :            :   const std::vector< EOS >& mat_blk) const
     117                 :            : // *************************************************************************
     118                 :            : //! \brief Calculate total energy based on the mixture composition
     119                 :            : //!   and species parameters.
     120                 :            : //! \param[in] mix_density Mixture density (sum of species density)
     121                 :            : //! \param[in] u Velocity component
     122                 :            : //! \param[in] v Velocity component
     123                 :            : //! \param[in] w Velocity component
     124                 :            : //! \param[in] mix_temp Mixture temperature (provided at call-site, since
     125                 :            : //!   it is reconstructed separately
     126                 :            : //! \param[in] mat_blk EOS material block
     127                 :            : //! \return Total energy
     128                 :            : // *************************************************************************
     129                 :            : {
     130                 :            :   // Compute mixture internal energy
     131                 :    3654294 :   tk::real mix_e = 0.;
     132         [ +  + ]:    9134610 :   for (std::size_t k = 0; k < m_nspec; k++) {
     133                 :    5480316 :     mix_e += m_Ys[k] * mat_blk[k].compute< EOS::internalenergy >(mix_temp);
     134                 :            :   }
     135                 :            : 
     136                 :            :   // Compute total energy
     137                 :    3654294 :   return mix_density * (mix_e + 0.5 * (u*u + v*v + w*w));
     138                 :            : }
     139                 :            : 
     140                 :            : tk::real
     141                 :    7009200 : Mixture::pressure(
     142                 :            :   tk::real mix_density,
     143                 :            :   tk::real mix_temp ) const
     144                 :            : // *************************************************************************
     145                 :            : //! \brief Calculate mixture pressure based on the mixture composition
     146                 :            : //!   and species parameters.
     147                 :            : //! \param[in] mix_density Mixture density (sum of species density)
     148                 :            : //! \param[in] mix_temp Mixture temperature
     149                 :            : //! \return Mixture pressure
     150                 :            : // *************************************************************************
     151                 :            : {
     152                 :            :   // Compute pressure based on the ideal gas EOS
     153                 :    7009200 :   return mix_density * m_mix_R * mix_temp;
     154                 :            : }
     155                 :            : 
     156                 :            : tk::real
     157                 :     691296 : Mixture::temperature(
     158                 :            :   tk::real mix_density,
     159                 :            :   tk::real u,
     160                 :            :   tk::real v,
     161                 :            :   tk::real w,
     162                 :            :   tk::real rhoE,
     163                 :            :   const std::vector< EOS >& mat_blk) const
     164                 :            : // *************************************************************************
     165                 :            : //! \brief Calculate temperature based on the mixture composition
     166                 :            : //!   and species parameters.
     167                 :            : //! \param[in] mix_density Mixture density (sum of species density)
     168                 :            : //! \param[in] u Velocity component
     169                 :            : //! \param[in] v Velocity component
     170                 :            : //! \param[in] w Velocity component
     171                 :            : //! \param[in] rhoE Total energy of the mixture
     172                 :            : //! \param[in] mat_blk EOS material block
     173                 :            : //! \return Mixture pressure
     174                 :            : // *************************************************************************
     175                 :            : {
     176                 :            :   // Compute internal energy
     177                 :     691296 :   tk::real e = rhoE / mix_density - 0.5 * (u*u + v*v + w*w);
     178                 :            : 
     179                 :            :   // Solve for temperature -- Newton's method
     180                 :     691296 :   tk::real temp = 1500; // Starting guess
     181                 :     691296 :   tk::real tol = std::max(1e-8, 1e-8 * e); // Stopping condition
     182                 :            :   tk::real err;
     183                 :     691296 :   std::size_t maxiter = 10;
     184                 :     691296 :   std::size_t i(0);
     185         [ +  - ]:    2243910 :   while (i < maxiter) {
     186                 :            :     // Construct f(T) = e(temp) - e
     187                 :    2243910 :     tk::real f_T = 0.;
     188         [ +  + ]:    5179116 :     for (std::size_t k = 0; k < m_nspec; k++) {
     189         [ +  - ]:    2935206 :       f_T += m_Ys[k] * mat_blk[k].compute< EOS::internalenergy >(temp);
     190                 :            :     }
     191                 :    2243910 :     f_T -= e;
     192                 :            : 
     193                 :            :     // Construct f'(T) = cv(temp)
     194                 :    2243910 :     tk::real fp_T = 0.;
     195         [ +  + ]:    5179116 :     for (std::size_t k = 0; k < m_nspec; k++) {
     196         [ +  - ]:    2935206 :       fp_T += m_Ys[k] * mat_blk[k].compute< EOS::cv >(temp);
     197                 :            :     }
     198                 :            : 
     199                 :            :     // Calculate next guess
     200                 :    2243910 :     temp = temp - f_T / fp_T;
     201                 :            : 
     202                 :            :     // Check stopping conditions
     203                 :    2243910 :     err = abs(f_T);
     204         [ +  + ]:    2243910 :     if (err <= tol) break;
     205                 :    1552614 :     i++;
     206         [ -  + ]:    1552614 :     if ( i == maxiter ) {
     207 [ -  - ][ -  - ]:          0 :       Throw("Mixture Newton's Method for temperature failed to converge after iterations "
         [ -  - ][ -  - ]
     208                 :            :       + std::to_string(i));
     209                 :            :     }
     210                 :            :   }
     211                 :            : 
     212                 :     691296 :   return temp;
     213                 :            : }
     214                 :            :