Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #pragma once
      11             : 
      12             : #include "VaporMixtureFluidProperties.h"
      13             : #include "SinglePhaseFluidProperties.h"
      14             : #include "THMIndicesGasMix.h"
      15             : 
      16             : class VaporMixtureFluidProperties;
      17             : 
      18             : namespace FlowModelGasMixUtils
      19             : {
      20             : 
      21             : /**
      22             :  * Computes the primitive solution vector from the conservative solution vector
      23             :  *
      24             :  * @param[in] U   Conservative solution vector
      25             :  * @param[in] fp  Fluid properties object
      26             :  */
      27             : template <bool is_ad>
      28             : std::vector<GenericReal<is_ad>>
      29        5609 : computePrimitiveSolution(const std::vector<GenericReal<is_ad>> & U,
      30             :                          const VaporMixtureFluidProperties & fp)
      31             : {
      32             :   const auto & xirhoA = U[THMGasMix1D::XIRHOA];
      33             :   const auto & rhoA = U[THMGasMix1D::RHOA];
      34             :   const auto & rhouA = U[THMGasMix1D::RHOUA];
      35             :   const auto & rhoEA = U[THMGasMix1D::RHOEA];
      36             :   const auto & A = U[THMGasMix1D::AREA];
      37             : 
      38        5608 :   const auto xi = xirhoA / rhoA;
      39        5608 :   const auto rho = rhoA / A;
      40        5608 :   const auto vel = rhouA / rhoA;
      41        5609 :   const auto v = 1.0 / rho;
      42        5609 :   const auto e = rhoEA / rhoA - 0.5 * vel * vel;
      43        5609 :   const auto p = fp.p_from_v_e(v, e, {xi});
      44        5609 :   const auto T = fp.T_from_v_e(v, e, {xi});
      45             : 
      46        5609 :   std::vector<GenericReal<is_ad>> W(THMGasMix1D::N_PRIM_VARS);
      47        5609 :   W[THMGasMix1D::MASS_FRACTION] = xi;
      48        5609 :   W[THMGasMix1D::PRESSURE] = p;
      49        5609 :   W[THMGasMix1D::VELOCITY] = vel;
      50        5609 :   W[THMGasMix1D::TEMPERATURE] = T;
      51             : 
      52        5609 :   return W;
      53           2 : }
      54             : 
      55             : /**
      56             :  * Computes the conservative solution vector from the primitive solution vector
      57             :  *
      58             :  * @param[in] W   Primitive solution vector
      59             :  * @param[in] A   Cross-sectional area
      60             :  * @param[in] fp  Fluid properties object
      61             :  */
      62             : template <bool is_ad>
      63             : std::vector<GenericReal<is_ad>>
      64           6 : computeConservativeSolution(const std::vector<GenericReal<is_ad>> & W,
      65             :                             const GenericReal<is_ad> & A,
      66             :                             const VaporMixtureFluidProperties & fp)
      67             : {
      68             :   const auto & xi = W[THMGasMix1D::MASS_FRACTION];
      69             :   const auto & p = W[THMGasMix1D::PRESSURE];
      70             :   const auto & T = W[THMGasMix1D::TEMPERATURE];
      71             :   const auto & vel = W[THMGasMix1D::VELOCITY];
      72             : 
      73          12 :   const auto rho = fp.rho_from_p_T(p, T, {xi});
      74             :   const auto rhoA = rho * A;
      75           6 :   const auto e = fp.e_from_p_rho(p, rho, {xi});
      76           6 :   const auto E = e + 0.5 * vel * vel;
      77             : 
      78           6 :   std::vector<GenericReal<is_ad>> U(THMGasMix1D::N_FLUX_INPUTS);
      79           6 :   U[THMGasMix1D::XIRHOA] = xi * rhoA;
      80           6 :   U[THMGasMix1D::RHOA] = rhoA;
      81           6 :   U[THMGasMix1D::RHOUA] = rhoA * vel;
      82           6 :   U[THMGasMix1D::RHOEA] = rhoA * E;
      83           6 :   U[THMGasMix1D::AREA] = A;
      84             : 
      85           6 :   return U;
      86          12 : }
      87             : 
      88             : /**
      89             :  * Computes the numerical flux vector from the primitive solution vector
      90             :  *
      91             :  * @param[in] W   Primitive solution vector
      92             :  * @param[in] A   Cross-sectional area
      93             :  * @param[in] fp  Fluid properties object
      94             :  */
      95             : template <bool is_ad>
      96             : std::vector<GenericReal<is_ad>>
      97           1 : computeFluxFromPrimitive(const std::vector<GenericReal<is_ad>> & W,
      98             :                          const GenericReal<is_ad> & A,
      99             :                          const VaporMixtureFluidProperties & fp)
     100             : {
     101             :   const auto & xi = W[THMGasMix1D::MASS_FRACTION];
     102             :   const auto & p = W[THMGasMix1D::PRESSURE];
     103             :   const auto & T = W[THMGasMix1D::TEMPERATURE];
     104             :   const auto & vel = W[THMGasMix1D::VELOCITY];
     105             : 
     106           1 :   const auto rho = fp.rho_from_p_T(p, T, {xi});
     107           1 :   const auto e = fp.e_from_p_rho(p, rho, {xi});
     108           1 :   const auto E = e + 0.5 * vel * vel;
     109             : 
     110           1 :   std::vector<ADReal> F(THMGasMix1D::N_FLUX_OUTPUTS, 0.0);
     111           1 :   F[THMGasMix1D::SPECIES] = xi * rho * vel * A;
     112           1 :   F[THMGasMix1D::MASS] = rho * vel * A;
     113           1 :   F[THMGasMix1D::MOMENTUM] = (rho * vel * vel + p) * A;
     114           1 :   F[THMGasMix1D::ENERGY] = vel * (rho * E + p) * A;
     115             : 
     116           1 :   return F;
     117           2 : }
     118             : 
     119             : /**
     120             :  * Gets the elemental conservative solution vector
     121             :  *
     122             :  * @param[in] elem         Element
     123             :  * @param[in] U_vars       Vector of conservative variable pointers
     124             :  * @param[in] is_implicit  Is implicit?
     125             :  */
     126             : template <bool is_ad>
     127             : std::vector<GenericReal<is_ad>>
     128           0 : getElementalSolutionVector(const Elem * elem,
     129             :                            const std::vector<MooseVariable *> & U_vars,
     130             :                            bool is_implicit)
     131             : {
     132             :   mooseAssert(elem, "The supplied element is a nullptr.");
     133             : 
     134           0 :   std::vector<GenericReal<is_ad>> U(THMGasMix1D::N_FLUX_INPUTS, 0.0);
     135             : 
     136           0 :   if (is_implicit)
     137             :   {
     138           0 :     for (const auto i : make_range(THMGasMix1D::N_FLUX_INPUTS))
     139             :     {
     140             :       mooseAssert(U_vars[i], "The supplied variable is a nullptr.");
     141           0 :       U[i] = U_vars[i]->getElementalValue(elem);
     142             :     }
     143             : 
     144             :     std::vector<dof_id_type> dof_indices;
     145             : 
     146           0 :     const std::vector<unsigned int> ind = {
     147             :         THMGasMix1D::XIRHOA, THMGasMix1D::RHOA, THMGasMix1D::RHOUA, THMGasMix1D::RHOEA};
     148           0 :     for (const auto j : index_range(ind))
     149             :     {
     150           0 :       const auto i = ind[j];
     151           0 :       U_vars[i]->dofMap().dof_indices(elem, dof_indices, U_vars[i]->number());
     152           0 :       Moose::derivInsert(U[i].derivatives(), dof_indices[0], 1.0);
     153             :     }
     154             :   }
     155             :   else
     156             :   {
     157           0 :     for (const auto i : make_range(THMGasMix1D::N_FLUX_INPUTS))
     158           0 :       U[i] = U_vars[i]->getElementalValueOld(elem);
     159             :   }
     160             : 
     161           0 :   return U;
     162             : }
     163             : 
     164             : template <bool is_ad>
     165             : GenericReal<is_ad>
     166       22540 : computeSecondaryMoleFraction(const GenericReal<is_ad> & xi_secondary,
     167             :                              const VaporMixtureFluidProperties & fp)
     168             : {
     169             :   mooseAssert(fp.getNumberOfSecondaryVapors() == 1,
     170             :               "This function assumes there is a single secondary fluid.");
     171       22540 :   const SinglePhaseFluidProperties & fp_primary = fp.getPrimaryFluidProperties();
     172       22540 :   const SinglePhaseFluidProperties & fp_secondary = fp.getSecondaryFluidProperties();
     173             : 
     174       22540 :   const GenericReal<is_ad> xi_primary = 1 - xi_secondary;
     175             : 
     176       22540 :   const GenericReal<is_ad> moles_primary = xi_primary / fp_primary.molarMass();
     177       45080 :   const GenericReal<is_ad> moles_secondary = xi_secondary / fp_secondary.molarMass();
     178             : 
     179       22540 :   return moles_secondary / (moles_primary + moles_secondary);
     180             : }
     181             : }