feat: Add compositional enthalpy model

src/coreComponents/constitutive/CMakeLists.txt

@@ -86,6 +86,7 @@ set( constitutive_headers
      fluid/multifluid/CO2Brine/functions/PVTFunctionHelpers.hpp
      fluid/multifluid/CO2Brine/functions/SpanWagnerCO2Density.hpp
      fluid/multifluid/CO2Brine/functions/WaterDensity.hpp
+     fluid/multifluid/compositional/functions/CompositionalMultiphaseFluidUtilities.hpp
      fluid/multifluid/compositional/functions/CompositionalProperties_impl.hpp
      fluid/multifluid/compositional/functions/CompositionalProperties.hpp
      fluid/multifluid/compositional/functions/CubicEOSPhaseModel_impl.hpp
@@ -102,6 +103,7 @@ set( constitutive_headers
      fluid/multifluid/compositional/functions/StabilityTest.hpp
      fluid/multifluid/compositional/functions/Utilities.hpp
      fluid/multifluid/compositional/models/CompositionalDensity.hpp
+     fluid/multifluid/compositional/models/CompositionalEnthalpy.hpp
      fluid/multifluid/compositional/models/ConstantViscosity.hpp
      fluid/multifluid/compositional/models/FunctionBase.hpp
      fluid/multifluid/compositional/models/ImmiscibleWaterDensity.hpp
@@ -267,6 +269,7 @@ set( constitutive_sources
      fluid/multifluid/CO2Brine/functions/PureWaterProperties.cpp
      fluid/multifluid/CO2Brine/functions/WaterDensity.cpp
      fluid/multifluid/compositional/models/CompositionalDensity.cpp
+     fluid/multifluid/compositional/models/CompositionalEnthalpy.cpp
      fluid/multifluid/compositional/models/ConstantViscosity.cpp
      fluid/multifluid/compositional/models/ImmiscibleWaterDensity.cpp
      fluid/multifluid/compositional/models/ImmiscibleWaterFlashModel.cpp

src/coreComponents/constitutive/fluid/multifluid/compositional/CompositionalMultiphaseFluid.hpp

@@ -29,7 +29,6 @@
 #include "constitutive/fluid/multifluid/compositional/models/KValueFlashModel.hpp"
 #include "constitutive/fluid/multifluid/compositional/models/LohrenzBrayClarkViscosity.hpp"
 #include "constitutive/fluid/multifluid/compositional/models/NegativeTwoPhaseFlashModel.hpp"
-#include "constitutive/fluid/multifluid/compositional/models/NullModel.hpp"
 #include "constitutive/fluid/multifluid/compositional/models/PhaseModel.hpp"
 #include "constitutive/fluid/multifluid/compositional/models/PhillipsBrineDensity.hpp"
 #include "constitutive/fluid/multifluid/compositional/models/PhillipsBrineViscosity.hpp"
@@ -50,14 +49,15 @@ namespace constitutive
 template< typename FLASH, typename PHASE1, typename PHASE2, typename PHASE3 = compositional::NullPhaseModel >
 class CompositionalMultiphaseFluid : public MultiFluidBase
 {
+  using Traits = compositional::CompositionalMultiphaseFluidTraits< FLASH, PHASE1, PHASE2, PHASE3 >;
 public:
   using FlashModel = FLASH;
   using Phase1Model = PHASE1;
   using Phase2Model = PHASE2;
   using Phase3Model = PHASE3;
 
   // Get the number of phases
-  static constexpr integer NUM_PHASES = FlashModel::KernelWrapper::getNumberOfPhases();
+  static constexpr integer NUM_PHASES = Traits::getNumberOfPhases();
   // Currently restrict to 2 or 3 phases
   static_assert( NUM_PHASES == 2 || NUM_PHASES == 3 );
 

src/coreComponents/constitutive/fluid/multifluid/compositional/CompositionalMultiphaseFluidUpdates.hpp

@@ -21,6 +21,7 @@
 #define GEOS_CONSTITUTIVE_FLUID_MULTIFLUID_COMPOSITIONAL_COMPOSITIONALMULTIPHASEFLUIDUPDATES_HPP_
 
 #include "constitutive/fluid/multifluid/compositional/parameters/ComponentProperties.hpp"
+#include "constitutive/fluid/multifluid/compositional/functions/CompositionalMultiphaseFluidUtilities.hpp"
 
 #include "constitutive/fluid/multifluid/MultiFluidBase.hpp"
 #include "constitutive/fluid/multifluid/MultiFluidUtils.hpp"
@@ -40,6 +41,8 @@ namespace constitutive
 template< typename FLASH, typename PHASE1, typename PHASE2, typename PHASE3 >
 class CompositionalMultiphaseFluidUpdates final : public MultiFluidBase::KernelWrapper
 {
+  using Traits = compositional::CompositionalMultiphaseFluidTraits< FLASH, PHASE1, PHASE2, PHASE3 >;
+  static constexpr integer NUM_PHASES = Traits::getNumberOfPhases();
 public:
   CompositionalMultiphaseFluidUpdates( compositional::ComponentProperties const & componentProperties,
                                        FLASH const & flash,
@@ -191,8 +194,7 @@ CompositionalMultiphaseFluidUpdates< FLASH, PHASE1, PHASE2, PHASE3 >::compute(
   MultiFluidBase::FluidProp::SliceType const totalDensity ) const
 {
   integer constexpr maxNumComp = MultiFluidBase::MAX_NUM_COMPONENTS;
-  integer constexpr maxNumPhase = MultiFluidBase::MAX_NUM_PHASES - 1;
-  MultiFluidBase::PhaseComp::StackValueType< maxNumPhase *maxNumComp > kValues( 1, 1, numPhases() - 1, numComponents() );
+  MultiFluidBase::PhaseComp::StackValueType< NUM_PHASES *maxNumComp > kValues( 1, 1, numPhases() - 1, numComponents() );
 
   LvArray::forValuesInSlice( kValues[0][0], setZero );   // Force initialisation of k-Values
 
@@ -230,7 +232,6 @@ CompositionalMultiphaseFluidUpdates< FLASH, PHASE1, PHASE2, PHASE3 >::compute(
 {
   integer constexpr maxNumComp = MultiFluidBase::MAX_NUM_COMPONENTS;
   integer constexpr maxNumDof = MultiFluidBase::MAX_NUM_COMPONENTS + 2;
-  integer constexpr maxNumPhase = MultiFluidBase::MAX_NUM_PHASES;
   integer const numComp = numComponents();
   integer const numPhase = numPhases();
   integer const numDof = numComp + 2;
@@ -283,7 +284,7 @@ CompositionalMultiphaseFluidUpdates< FLASH, PHASE1, PHASE2, PHASE3 >::compute(
                             phaseMassDensity.value[m_phaseOrder[1]],
                             phaseMassDensity.derivs[m_phaseOrder[1]],
                             m_useMass );
-  if constexpr (2 < FLASH::KernelWrapper::getNumberOfPhases())
+  if constexpr (2 < NUM_PHASES)
   {
     m_phase3.density.compute( m_componentProperties,
                               pressure,
@@ -315,7 +316,7 @@ CompositionalMultiphaseFluidUpdates< FLASH, PHASE1, PHASE2, PHASE3 >::compute(
                               phaseViscosity.value[m_phaseOrder[1]],
                               phaseViscosity.derivs[m_phaseOrder[1]],
                               m_useMass );
-  if constexpr (2 < FLASH::KernelWrapper::getNumberOfPhases())
+  if constexpr (2 < NUM_PHASES)
   {
     m_phase3.viscosity.compute( m_componentProperties,
                                 pressure,
@@ -328,9 +329,38 @@ CompositionalMultiphaseFluidUpdates< FLASH, PHASE1, PHASE2, PHASE3 >::compute(
                                 m_useMass );
   }
 
-  // 5. Convert derivatives from phase composition to total composition
+  // 5. Calculate the phase enthapies
+  if constexpr (Traits::isThermalType())
+  {
+    m_phase1.enthalpy.compute( m_componentProperties,
+                               pressure,
+                               temperature,
+                               phaseCompFrac.value[m_phaseOrder[0]].toSliceConst(),
+                               phaseEnthalpy.value[m_phaseOrder[0]],
+                               phaseEnthalpy.derivs[m_phaseOrder[0]],
+                               m_useMass );
+    m_phase2.enthalpy.compute( m_componentProperties,
+                               pressure,
+                               temperature,
+                               phaseCompFrac.value[m_phaseOrder[1]].toSliceConst(),
+                               phaseEnthalpy.value[m_phaseOrder[1]],
+                               phaseEnthalpy.derivs[m_phaseOrder[1]],
+                               m_useMass );
+    if constexpr (2 < NUM_PHASES)
+    {
+      m_phase3.enthalpy.compute( m_componentProperties,
+                                 pressure,
+                                 temperature,
+                                 phaseCompFrac.value[m_phaseOrder[2]].toSliceConst(),
+                                 phaseEnthalpy.value[m_phaseOrder[2]],
+                                 phaseEnthalpy.derivs[m_phaseOrder[2]],
+                                 m_useMass );
+    }
+  }
+
+  // 6. Convert derivatives from phase composition to total composition
   stackArray1d< real64, maxNumDof > workSpace( numDof );
-  for( integer ip = 0; ip < FLASH::KernelWrapper::getNumberOfPhases(); ++ip )
+  for( integer ip = 0; ip < NUM_PHASES; ++ip )
   {
     convertDerivativesToTotalMoleFraction( numComp,
                                            phaseComponentFraction.derivs[ip].toSliceConst(),
@@ -344,9 +374,16 @@ CompositionalMultiphaseFluidUpdates< FLASH, PHASE1, PHASE2, PHASE3 >::compute(
                                            phaseComponentFraction.derivs[ip].toSliceConst(),
                                            phaseViscosity.derivs[ip],
                                            workSpace );
+    if constexpr (Traits::isThermalType())
+    {
+      convertDerivativesToTotalMoleFraction( numComp,
+                                             phaseCompFrac.derivs[ip].toSliceConst(),
+                                             phaseEnthalpy.derivs[ip],
+                                             workSpace );
+    }
   }
 
-  // 6. if mass variables used instead of molar, perform the conversion
+  // 7. if mass variables used instead of molar, perform the conversion
   if( m_useMass )
   {
     real64 phaseMolecularWeight[maxNumPhase]{};
@@ -390,7 +427,7 @@ CompositionalMultiphaseFluidUpdates< FLASH, PHASE1, PHASE2, PHASE3 >::compute(
     }
   }
 
-  // 7. Compute total fluid mass/molar density and derivatives
+  // 8. Compute total fluid mass/molar density and derivatives
 
   computeTotalDensity( phaseFraction,
                        phaseDensity,

src/coreComponents/constitutive/fluid/multifluid/compositional/functions/CompositionalMultiphaseFluidUtilities.hpp

@@ -0,0 +1,60 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2024 TotalEnergies
+ * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2023-2024 Chevron
+ * Copyright (c) 2019-     GEOS/GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file CompositionalMultiphaseFluidUtilities.hpp
+ */
+
+#ifndef GEOS_CONSTITUTIVE_FLUID_MULTIFLUID_COMPOSITIONAL_FUNCTIONS_COMPOSITIONALMULTIPHASEFLUIDUTILITIES_HPP_
+#define GEOS_CONSTITUTIVE_FLUID_MULTIFLUID_COMPOSITIONAL_FUNCTIONS_COMPOSITIONALMULTIPHASEFLUIDUTILITIES_HPP_
+
+#include "constitutive/fluid/multifluid/compositional/models/NullModel.hpp"
+
+namespace geos
+{
+namespace constitutive
+{
+namespace compositional
+{
+
+template< typename FLASH, typename PHASE1, typename PHASE2, typename PHASE3 >
+struct CompositionalMultiphaseFluidTraits
+{
+  static constexpr integer getNumberOfPhases()
+  {
+    return FLASH::KernelWrapper::getNumberOfPhases();
+  }
+
+  static constexpr bool isThermalType()
+  {
+    if constexpr (getNumberOfPhases() == 3)
+    {
+      return !( std::is_same_v< typename PHASE1::Enthalpy, compositional::NullModel > ||
+                std::is_same_v< typename PHASE2::Enthalpy, compositional::NullModel > ||
+                std::is_same_v< typename PHASE3::Enthalpy, compositional::NullModel > );
+    }
+    else
+    {
+      return !( std::is_same_v< typename PHASE1::Enthalpy, compositional::NullModel > ||
+                std::is_same_v< typename PHASE2::Enthalpy, compositional::NullModel > );
+    }
+  }
+};
+
+} //namespace compositional
+} // namespace constitutive
+} // namespace geos
+
+#endif //GEOS_CONSTITUTIVE_FLUID_MULTIFLUID_COMPOSITIONAL_FUNCTIONS_COMPOSITIONALMULTIPHASEFLUIDUTILITIES_HPP_

src/coreComponents/constitutive/fluid/multifluid/compositional/models/CompositionalEnthalpy.cpp

@@ -0,0 +1,91 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2024 TotalEnergies
+ * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2023-2024 Chevron
+ * Copyright (c) 2019-     GEOS/GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file CompositionalEnthalpy.cpp
+ */
+
+#include "CompositionalEnthalpy.hpp"
+#include "constitutive/fluid/multifluid/compositional/parameters/HeatCapacityCoefficients.hpp"
+#include "common/PhysicsConstants.hpp"
+
+namespace geos
+{
+namespace constitutive
+{
+namespace compositional
+{
+
+CompositionalEnthalpyUpdate::CompositionalEnthalpyUpdate( EquationOfStateType const equationOfState,
+                                                          arrayView1d< real64 const > const & referenceEnthalpy,
+                                                          arrayView2d< real64 const > const & coefficients )
+  : m_equationOfState( equationOfState ),
+  m_referenceEnthalpy( referenceEnthalpy ),
+  m_coefficients( coefficients )
+{}
+
+CompositionalEnthalpy::CompositionalEnthalpy( string const & name,
+                                              ComponentProperties const & componentProperties,
+                                              integer const phaseIndex,
+                                              ModelParameters const & modelParameters )
+  : FunctionBase( name, componentProperties ),
+  m_heatCapacityCoefficients ( modelParameters.get< HeatCapacityCoefficients >() )
+{
+  EquationOfState const * equationOfState = modelParameters.get< EquationOfState >();
+  string const eosName = equationOfState->m_equationsOfStateNames[phaseIndex];
+  m_equationOfState = EnumStrings< EquationOfStateType >::fromString( eosName );
+
+  const auto & heatCapacityCoefficients = m_heatCapacityCoefficients->m_coefficients;
+  integer const numComps = heatCapacityCoefficients.size( 1 );
+  integer const numCoeffs = heatCapacityCoefficients.size( 2 );
+  m_referenceEnthalpy.resize( numComps );
+  m_coefficients.resize( numComps, numCoeffs );
+  constexpr real64 R = constants::gasConstant;
+  for( integer ic = 0; ic < numComps; ic++ )
+  {
+    for( integer j = 0; j < numCoeffs; j++ )
+    {
+      m_coefficients( ic, j ) = R * heatCapacityCoefficients( phaseIndex, ic, j );
+    }
+    // Calculate the enthalpy at the reference temperature
+    real64 refEnthalpy = 0.0;
+    real64 refHeatCapacity = 0.0;
+    KernelWrapper::evaluatePolynomial( m_heatCapacityCoefficients->m_referenceTemperature[ic],
+                                       m_coefficients[ic],
+                                       refEnthalpy,
+                                       refHeatCapacity );
+    m_referenceEnthalpy[ic] = m_heatCapacityCoefficients->m_referenceEnthalpy( phaseIndex, ic ) - refEnthalpy;
+  }
+}
+
+CompositionalEnthalpy::KernelWrapper
+CompositionalEnthalpy::createKernelWrapper() const
+{
+  return KernelWrapper( m_equationOfState,
+                        m_referenceEnthalpy.toViewConst(),
+                        m_coefficients.toViewConst());
+}
+
+std::unique_ptr< ModelParameters >
+CompositionalEnthalpy::createParameters( std::unique_ptr< ModelParameters > parameters )
+{
+  auto params = EquationOfState::create( std::move( parameters ) );
+  params = HeatCapacityCoefficients::create( std::move( params ) );
+  return params;
+}
+
+} // namespace compositional
+} // namespace constitutive
+} // namespace geos