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             : #include "NumericalFluxGasMixHLLC.h"
      11             : #include "VaporMixtureFluidProperties.h"
      12             : #include "THMIndicesGasMix.h"
      13             : #include "Numerics.h"
      14             : 
      15             : registerMooseObject("ThermalHydraulicsApp", NumericalFluxGasMixHLLC);
      16             : 
      17             : InputParameters
      18          45 : NumericalFluxGasMixHLLC::validParams()
      19             : {
      20          45 :   InputParameters params = NumericalFluxGasMixBase::validParams();
      21          45 :   params += NaNInterface::validParams();
      22          90 :   params.addRequiredParam<UserObjectName>("fluid_properties",
      23             :                                           "Name for fluid properties user object");
      24          45 :   params.addClassDescription("Computes internal side flux for the 1-D, 1-phase, variable-area "
      25             :                              "Euler equations using the HLLC approximate Riemann solver.");
      26          45 :   return params;
      27           0 : }
      28             : 
      29          24 : NumericalFluxGasMixHLLC::NumericalFluxGasMixHLLC(const InputParameters & parameters)
      30             :   : NumericalFluxGasMixBase(parameters),
      31             :     NaNInterface(this),
      32          24 :     _fp(getUserObject<VaporMixtureFluidProperties>("fluid_properties"))
      33             : {
      34          24 : }
      35             : 
      36             : void
      37       11736 : NumericalFluxGasMixHLLC::calcFlux(const std::vector<ADReal> & UL,
      38             :                                   const std::vector<ADReal> & UR,
      39             :                                   const RealVectorValue & nLR,
      40             :                                   const RealVectorValue & t1,
      41             :                                   const RealVectorValue & t2,
      42             :                                   std::vector<ADReal> & FL,
      43             :                                   std::vector<ADReal> & FR) const
      44             : {
      45             :   // extract the conserved variables and area
      46             : 
      47       11736 :   const ADReal xirhoAL = UL[THMGasMix3D::XIRHOA];
      48       11736 :   const ADReal rhoAL = UL[THMGasMix3D::RHOA];
      49       11736 :   const ADReal rhouAL = UL[THMGasMix3D::RHOUA];
      50       11736 :   const ADReal rhovAL = UL[THMGasMix3D::RHOVA];
      51       11736 :   const ADReal rhowAL = UL[THMGasMix3D::RHOWA];
      52       11736 :   const ADReal rhoEAL = UL[THMGasMix3D::RHOEA];
      53       11736 :   const ADReal AL = UL[THMGasMix3D::AREA];
      54             : 
      55       11736 :   const ADReal xirhoAR = UR[THMGasMix3D::XIRHOA];
      56       11736 :   const ADReal rhoAR = UR[THMGasMix3D::RHOA];
      57       11736 :   const ADReal rhouAR = UR[THMGasMix3D::RHOUA];
      58       11736 :   const ADReal rhovAR = UR[THMGasMix3D::RHOVA];
      59       11736 :   const ADReal rhowAR = UR[THMGasMix3D::RHOWA];
      60       11736 :   const ADReal rhoEAR = UR[THMGasMix3D::RHOEA];
      61       11736 :   const ADReal AR = UR[THMGasMix3D::AREA];
      62             : 
      63             :   // compute the primitive variables
      64             : 
      65             :   const ADReal rhoL = rhoAL / AL;
      66       11736 :   const ADRealVectorValue uvecL(rhouAL / rhoAL, rhovAL / rhoAL, rhowAL / rhoAL);
      67       11736 :   const ADReal unL = uvecL * nLR;
      68       11736 :   const ADReal ut1L = uvecL * t1;
      69       11736 :   const ADReal ut2L = uvecL * t2;
      70             :   const ADReal rhoEL = rhoEAL / AL;
      71       23472 :   const ADReal vL = 1.0 / rhoL;
      72             :   const ADReal EL = rhoEAL / rhoAL;
      73       23472 :   const ADReal eL = EL - 0.5 * uvecL * uvecL;
      74             :   const ADReal xiL = xirhoAL / rhoAL;
      75       11736 :   const ADReal pL = _fp.p_from_v_e(vL, eL, {xiL});
      76       11736 :   const ADReal TL = _fp.T_from_v_e(vL, eL, {xiL});
      77       23472 :   const ADReal cL = _fp.c_from_p_T(pL, TL, {xiL});
      78             : 
      79             :   const ADReal rhoR = rhoAR / AR;
      80       11736 :   const ADRealVectorValue uvecR(rhouAR / rhoAR, rhovAR / rhoAR, rhowAR / rhoAR);
      81       11736 :   const ADReal unR = uvecR * nLR;
      82       11736 :   const ADReal ut1R = uvecR * t1;
      83       11736 :   const ADReal ut2R = uvecR * t2;
      84             :   const ADReal rhoER = rhoEAR / AR;
      85       23472 :   const ADReal vR = 1.0 / rhoR;
      86             :   const ADReal ER = rhoEAR / rhoAR;
      87       23472 :   const ADReal eR = ER - 0.5 * uvecR * uvecR;
      88             :   const ADReal xiR = xirhoAR / rhoAR;
      89       11736 :   const ADReal pR = _fp.p_from_v_e(vR, eR, {xiR});
      90       11736 :   const ADReal TR = _fp.T_from_v_e(vR, eR, {xiR});
      91       23472 :   const ADReal cR = _fp.c_from_p_T(pR, TR, {xiR});
      92             : 
      93             :   // compute wave speeds
      94       11736 :   const ADReal sL = std::min(unL - cL, unR - cR);
      95       11736 :   const ADReal sR = std::max(unL + cL, unR + cR);
      96       11736 :   const ADReal sm = (rhoR * unR * (sR - unR) - rhoL * unL * (sL - unL) + pL - pR) /
      97       11736 :                     (rhoR * (sR - unR) - rhoL * (sL - unL));
      98             : 
      99             :   // compute Omega_L, Omega_R
     100       23472 :   const ADReal omegL = 1.0 / (sL - sm);
     101       23472 :   const ADReal omegR = 1.0 / (sR - sm);
     102             : 
     103             :   // compute p^*
     104       11736 :   const ADReal ps = rhoL * (sL - unL) * (sm - unL) + pL;
     105             : 
     106             :   // compute U_L^*, U_R^*
     107             : 
     108       11736 :   const ADReal rhoLs = omegL * (sL - unL) * rhoL;
     109       11736 :   const ADReal rhounLs = omegL * ((sL - unL) * rhoL * unL + ps - pL);
     110       11736 :   const ADReal rhoELs = omegL * ((sL - unL) * rhoEL - pL * unL + ps * sm);
     111             : 
     112       11736 :   const ADReal rhoRs = omegR * (sR - unR) * rhoR;
     113       11736 :   const ADReal rhounRs = omegR * ((sR - unR) * rhoR * unR + ps - pR);
     114       11736 :   const ADReal rhoERs = omegR * ((sR - unR) * rhoER - pR * unR + ps * sm);
     115             : 
     116       11736 :   std::vector<ADReal> UL_1d(THMGasMix1D::N_FLUX_INPUTS);
     117       11736 :   UL_1d[THMGasMix1D::XIRHOA] = UL[THMGasMix3D::XIRHOA];
     118       11736 :   UL_1d[THMGasMix1D::RHOA] = UL[THMGasMix3D::RHOA];
     119       11736 :   UL_1d[THMGasMix1D::RHOUA] = UL[THMGasMix3D::RHOUA];
     120       11736 :   UL_1d[THMGasMix1D::RHOEA] = UL[THMGasMix3D::RHOEA];
     121       11736 :   UL_1d[THMGasMix1D::AREA] = UL[THMGasMix3D::AREA];
     122             : 
     123       11736 :   std::vector<ADReal> UR_1d(THMGasMix1D::N_FLUX_INPUTS);
     124       11736 :   UR_1d[THMGasMix1D::XIRHOA] = UR[THMGasMix3D::XIRHOA];
     125       11736 :   UR_1d[THMGasMix1D::RHOA] = UR[THMGasMix3D::RHOA];
     126       11736 :   UR_1d[THMGasMix1D::RHOUA] = UR[THMGasMix3D::RHOUA];
     127       11736 :   UR_1d[THMGasMix1D::RHOEA] = UR[THMGasMix3D::RHOEA];
     128       11736 :   UR_1d[THMGasMix1D::AREA] = UR[THMGasMix3D::AREA];
     129             : 
     130       11736 :   const ADReal A_flow = computeFlowArea(UL_1d, UR_1d);
     131             : 
     132             :   // compute the fluxes
     133       11736 :   FL.resize(THMGasMix3D::N_FLUX_OUTPUTS);
     134       11736 :   if (sL > 0.0)
     135             :   {
     136           1 :     FL[THMGasMix3D::SPECIES] = unL * xiL * rhoL * A_flow;
     137           1 :     FL[THMGasMix3D::MASS] = unL * rhoL * A_flow;
     138           1 :     FL[THMGasMix3D::MOM_NORM] = (unL * rhoL * unL + pL) * A_flow;
     139           1 :     FL[THMGasMix3D::MOM_TAN1] = rhoL * unL * ut1L * A_flow;
     140           1 :     FL[THMGasMix3D::MOM_TAN2] = rhoL * unL * ut2L * A_flow;
     141           1 :     FL[THMGasMix3D::ENERGY] = unL * (rhoEL + pL) * A_flow;
     142             : 
     143           1 :     _last_region_index = 0;
     144             :   }
     145       11735 :   else if (sL <= 0.0 && sm > 0.0)
     146             :   {
     147        1113 :     FL[THMGasMix3D::SPECIES] = sm * xiL * rhoLs * A_flow;
     148        1113 :     FL[THMGasMix3D::MASS] = sm * rhoLs * A_flow;
     149        1113 :     FL[THMGasMix3D::MOM_NORM] = (sm * rhounLs + ps) * A_flow;
     150        1113 :     FL[THMGasMix3D::MOM_TAN1] = rhounLs * ut1L * A_flow;
     151        1113 :     FL[THMGasMix3D::MOM_TAN2] = rhounLs * ut2L * A_flow;
     152        1113 :     FL[THMGasMix3D::ENERGY] = sm * (rhoELs + ps) * A_flow;
     153             : 
     154        1113 :     _last_region_index = 1;
     155             :   }
     156       10622 :   else if (sm <= 0.0 && sR >= 0.0)
     157             :   {
     158       10621 :     FL[THMGasMix3D::SPECIES] = sm * xiR * rhoRs * A_flow;
     159       10621 :     FL[THMGasMix3D::MASS] = sm * rhoRs * A_flow;
     160       10621 :     FL[THMGasMix3D::MOM_NORM] = (sm * rhounRs + ps) * A_flow;
     161       10621 :     FL[THMGasMix3D::MOM_TAN1] = rhounRs * ut1R * A_flow;
     162       10621 :     FL[THMGasMix3D::MOM_TAN2] = rhounRs * ut2R * A_flow;
     163       10621 :     FL[THMGasMix3D::ENERGY] = sm * (rhoERs + ps) * A_flow;
     164             : 
     165       10621 :     _last_region_index = 2;
     166             :   }
     167           1 :   else if (sR < 0.0)
     168             :   {
     169           1 :     FL[THMGasMix3D::SPECIES] = unR * xiR * rhoR * A_flow;
     170           1 :     FL[THMGasMix3D::MASS] = unR * rhoR * A_flow;
     171           1 :     FL[THMGasMix3D::MOM_NORM] = (unR * rhoR * unR + pR) * A_flow;
     172           1 :     FL[THMGasMix3D::MOM_TAN1] = rhoR * unR * ut1R * A_flow;
     173           1 :     FL[THMGasMix3D::MOM_TAN2] = rhoR * unR * ut2R * A_flow;
     174           1 :     FL[THMGasMix3D::ENERGY] = unR * (rhoER + pR) * A_flow;
     175             : 
     176           1 :     _last_region_index = 3;
     177             :   }
     178             :   else
     179           0 :     std::fill(FL.begin(), FL.end(), getNaN());
     180             : 
     181       11736 :   FR = FL;
     182             : 
     183             :   const ADReal A_wall_L = AL - A_flow;
     184       11736 :   FL[THMGasMix3D::MOM_NORM] += pL * A_wall_L;
     185             : 
     186             :   const ADReal A_wall_R = AR - A_flow;
     187       11736 :   FR[THMGasMix3D::MOM_NORM] += pR * A_wall_R;
     188       82152 : }
     189             : 
     190             : ADReal
     191       11736 : NumericalFluxGasMixHLLC::computeFlowArea(const std::vector<ADReal> & UL,
     192             :                                          const std::vector<ADReal> & UR) const
     193             : {
     194       11736 :   return std::min(UL[THMGasMix1D::AREA], UR[THMGasMix1D::AREA]);
     195             : }