feat: Use stability test in flash (#3206)

src/coreComponents/constitutive/CMakeLists.txt

@@ -74,6 +74,7 @@ set( constitutive_headers
      fluid/multifluid/compositional/models/ComponentProperties.hpp
      fluid/multifluid/compositional/models/CompositionalDensity.hpp
      fluid/multifluid/compositional/models/ConstantViscosity.hpp
+     fluid/multifluid/compositional/models/CriticalVolume.hpp
      fluid/multifluid/compositional/models/EquationOfState.hpp
      fluid/multifluid/compositional/models/FunctionBase.hpp
      fluid/multifluid/compositional/models/LohrenzBrayClarkViscosity.hpp
@@ -218,6 +219,7 @@ set( constitutive_sources
      fluid/multifluid/CO2Brine/functions/WaterDensity.cpp
      fluid/multifluid/compositional/models/CompositionalDensity.cpp
      fluid/multifluid/compositional/models/ConstantViscosity.cpp
+     fluid/multifluid/compositional/models/CriticalVolume.cpp
      fluid/multifluid/compositional/models/LohrenzBrayClarkViscosity.cpp
      fluid/multifluid/compositional/models/NegativeTwoPhaseFlashModel.cpp
      fluid/multifluid/compositional/CompositionalMultiphaseFluid.cpp

src/coreComponents/constitutive/docs/CompositionalMultiphaseFluid.rst

@@ -119,6 +119,77 @@ where :math:`V_{ci}` and :math:`T_{ci}` are respectively the critical volume and
 :math:`i`. This is compared to the current temperature :math:`T` such that if :math:`T_{cp}<T` then the mixture
 is labeled as vapor and as liquid otherwise.
 
+Negative two-phase flash
+------------------------
+When a cell is identified as having an unstable mixture, it is necessary to determine the amounts in the liquid
+and vapor phases through phase splitting. This phase split is calculated by ensuring that the two phases are
+in thermodynamic equilibrium. For a system to be in thermodynamic equilibrium, the fugacities of each component
+in both the liquid and vapor phases must be equal:
+
+.. math::
+    \phi_{iL} = \phi_{iV} \hspace{1cm} i=1,2,3,\ldots,N_c
+
+where :math:`\phi_{iL}` is the fugacity of component :math:`i` in the liquid phase and :math:`\phi_{iL}` is the
+fugacity of component :math:`i` in the vapor phase.
+
+Fugacities are functions of temperature, pressure, and composition:
+
+.. math::
+    \phi_{iL} = \phi_{iL}(T, p, x_i) \hspace{1cm} i=1,2,3,\ldots,N_c
+
+and 
+
+.. math::
+    \phi_{iV} = \phi_{iV}(T, p, y_i) \hspace{1cm} i=1,2,3,\ldots,N_c
+
+and are calculated directly from an equation of state.
+
+Equilibrium constants, also known as K-values, are defined for each component as:
+
+.. math::
+    K_i = \frac{y_i}{x_i}
+
+where :math:`x_i` is the mole fraction of component :math:`i` in the liquid phase and :math:`y_i` is the
+mole fraction of component :math:`i` in the vapor phase. If we denote :math:`V` as the mole fraction
+of the vapor phase, the material balance indicates that the mole fractions of each component in the liquid
+and vapor phases are given by:
+
+.. math::
+    x_i = \frac{z_i}{1 + (K_i - 1)V}
+
+and
+
+.. math::
+    y_i = \frac{K_i z_i}{1 + (K_i - 1)V}
+
+The value of :math:`V` corresponding to a given set of K-values is determined by solving the
+so called Rachford and-Rice equation:
+
+.. math::
+    F(V) = \sum_{i=1}^{N_c} \left(x_i - y_i\right) = \sum_{i=1}^{N_c} \frac{z_i(1 - K_i)}{1 + (K_i - 1)V} = 0
+
+The flash calculation process is as follows:
+
+#. Once the mixture is confirmed to be stable, an initial set of K-values is chosen, typically using Wilson's formula.
+
+#. Given :math:`z_i` and :math:`K_i`, the Rachford-Rice equation is solved to determine the molar fraction of vapor,  :math:`V`. This is initially solved using successive substitution, followed by Newton iterations once the residual is sufficiently reduced.
+
+#. After  :math:`V` is calculated, the corresponding liquid and vapor mole fractions, :math:`x_i` and :math:`y_i`, are computed.
+
+#. These phase compositions are then used to calculate the component fugacities :math:`\phi_{iL}` and :math:`\phi_{iV}` in the liquid and vapor phases using the equation of state.
+
+#. Convergence is reached when the fugacities are equal for all components. The convergence criterion is defined as:
+
+   .. math::
+       \sum_{i=1}^{N_c} \left( \phi_{iL} - \phi_{iV} \right)^2 < \varepsilon
+   
+   where :math:`\varepsilon` is the convergence tolerance.
+
+#. If convergence is not achieved, successive substitution is used to update the set of K-values for the next iteration. The new K-values at iteration  :math:`t+1` are given by:
+
+   .. math::
+       K_i^{(t+1)} = K_i^{(t)} \frac{\phi_{iL}}{\phi_{iV}}
+
 Parameters
 =========================
 

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

@@ -196,6 +196,8 @@ CompositionalMultiphaseFluid< FLASH, PHASE1, PHASE2, PHASE3 >::deliverClone( str
                                                                              Group * const parent ) const
 {
   std::unique_ptr< ConstitutiveBase > clone = MultiFluidBase::deliverClone( name, parent );
+  CompositionalMultiphaseFluid & newFluid = dynamicCast< CompositionalMultiphaseFluid & >( *clone );
+  newFluid.createModels();
   return clone;
 }
 

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

@@ -64,39 +64,6 @@ struct KValueInitialization
     }
   }
 
-  /**
-   * @brief Calculate gas-liquid k-values near the convergence pressure
-   * @param[in] numComps number of components
-   * @param[in] pressure pressure
-   * @param[in] temperature temperature
-   * @param[in] componentProperties The compositional component properties
-   * @param[out] kValues the calculated k-values
-   **/
-  template< integer USD >
-  GEOS_HOST_DEVICE
-  GEOS_FORCE_INLINE
-  static void
-  computeConstantLiquidKvalue( integer const numComps,
-                               real64 const pressure,
-                               real64 const temperature,
-                               ComponentProperties::KernelWrapper const & componentProperties,
-                               arraySlice1d< real64, USD > const & kValues )
-  {
-    GEOS_UNUSED_VAR( pressure, temperature );
-    arrayView1d< real64 const > const & criticalPressure = componentProperties.m_componentCriticalPressure;
-    real64 averagePressure = 0.0; // Average pressure
-    for( integer ic = 0; ic < numComps; ++ic )
-    {
-      averagePressure += criticalPressure[ic];
-    }
-    averagePressure /= numComps;
-    constexpr real64 kValueGap = 0.01;
-    for( integer ic = 0; ic < numComps; ++ic )
-    {
-      kValues[ic] = criticalPressure[ic] < averagePressure ? 1.0/(1.0 + kValueGap) : 1.0/(1.0 - kValueGap);
-    }
-  }
-
 /**
  * @brief Calculate water-gas k-value
  * @param[in] pressure pressure

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

@@ -257,9 +257,6 @@ bool NegativeTwoPhaseFlash::compute( integer const numComps,
     }
   }
 
-  bool kValueReset = true;
-  constexpr real64 boundsTolerance = 1.0e-3;
-
   if( needInitialisation )
   {
     KValueInitialization::computeWilsonGasLiquidKvalue( numComps,
@@ -271,8 +268,6 @@ bool NegativeTwoPhaseFlash::compute( integer const numComps,
 
   auto const presentComponents = componentIndices.toSliceConst();
 
-  real64 const initialVapourFraction = RachfordRice::solve( kVapourLiquid.toSliceConst(), composition, presentComponents );
-
   bool converged = false;
   for( localIndex iterationCount = 0; iterationCount < MultiFluidConstants::maxSSIIterations; ++iterationCount )
   {
@@ -301,23 +296,9 @@ bool NegativeTwoPhaseFlash::compute( integer const numComps,
     }
 
     // Update K-values
-    if( (vapourPhaseMoleFraction < -boundsTolerance || 1.0-vapourPhaseMoleFraction < -boundsTolerance)
-        && 0.2 < LvArray::math::abs( vapourPhaseMoleFraction-initialVapourFraction )
-        && !kValueReset )
-    {
-      KValueInitialization::computeConstantLiquidKvalue( numComps,
-                                                         pressure,
-                                                         temperature,
-                                                         componentProperties,
-                                                         kVapourLiquid );
-      kValueReset =  true;
-    }
-    else
+    for( integer ic = 0; ic < numComps; ++ic )
     {
-      for( integer ic = 0; ic < numComps; ++ic )
-      {
-        kVapourLiquid[ic] *= exp( fugacityRatios[ic] );
-      }
+      kVapourLiquid[ic] *= exp( fugacityRatios[ic] );
     }
   }
 

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

@@ -52,7 +52,7 @@ struct StabilityTest
    * @param[out] kValues the k-values estimated from the stationary points
    * @return a flag indicating that 2 stationary points have been found
    */
-  template< integer USD1 >
+  template< integer USD1, integer USD2 >
   GEOS_HOST_DEVICE
   static bool compute( integer const numComps,
                        real64 const pressure,
@@ -61,25 +61,23 @@ struct StabilityTest
                        ComponentProperties::KernelWrapper const & componentProperties,
                        EquationOfStateType const & equationOfState,
                        real64 & tangentPlaneDistance,
-                       arraySlice1d< real64 > const & kValues )
+                       arraySlice1d< real64, USD2 > const & kValues )
   {
     constexpr integer numTrials = 2;    // Trial compositions
-    stackArray1d< real64, maxNumComps > logFugacity( numComps );
-    stackArray1d< real64, maxNumComps > normalizedComposition( numComps );
-    stackArray2d< real64, numTrials *maxNumComps > trialComposition( numTrials, numComps );
-    stackArray1d< real64, maxNumComps > logTrialComposition( numComps );
-    stackArray1d< real64, maxNumComps > hyperplane( numComps );     // h-parameter
+    stackArray2d< real64, 4*maxNumComps > workSpace( 4, numComps );
+    arraySlice1d< real64 > logFugacity = workSpace[0];
+    arraySlice1d< real64 > normalizedComposition = workSpace[1];
+    arraySlice1d< real64 > logTrialComposition = workSpace[2];
+    arraySlice1d< real64 > hyperplane = workSpace[3]; // h-parameter
     stackArray1d< integer, maxNumComps > availableComponents( numComps );
 
     calculatePresentComponents( numComps, composition, availableComponents );
     auto const presentComponents = availableComponents.toSliceConst();
 
+    LvArray::forValuesInSlice( workSpace.toSlice(), []( real64 & a ){ a = 0.0; } );
+
     // Calculate the hyperplane parameter
     // h_i = log( z_i ) + log( phi_i )
-    for( integer ic = 0; ic < numComps; ++ic )
-    {
-      hyperplane[ic] = 0.0;
-    }
     FugacityCalculator::computeLogFugacity( numComps,
                                             pressure,
                                             temperature,
@@ -99,43 +97,21 @@ struct StabilityTest
                                                         componentProperties,
                                                         kValues );
 
-    for( integer ic = 0; ic < numComps; ++ic )
-    {
-      trialComposition( 0, ic ) = composition[ic] / kValues[ic];
-      trialComposition( 1, ic ) = composition[ic] * kValues[ic];
-    }
-
-    integer numberOfStationaryPoints = 0;
     tangentPlaneDistance = LvArray::NumericLimits< real64 >::max;
     for( integer trialIndex = 0; trialIndex < numTrials; ++trialIndex )
     {
-      for( integer ic = 0; ic < numComps; ++ic )
-      {
-        normalizedComposition[ic] = trialComposition( trialIndex, ic );
-      }
-      normalizeComposition( numComps, normalizedComposition.toSlice() );
-
-      FugacityCalculator::computeLogFugacity( numComps,
-                                              pressure,
-                                              temperature,
-                                              normalizedComposition.toSliceConst(),
-                                              componentProperties,
-                                              equationOfState,
-                                              logFugacity );
+      // Initialise next sample
+      real64 const alpha = static_cast< real64 >(trialIndex)/(numTrials-1);
       for( integer const ic : presentComponents )
       {
-        logTrialComposition[ic] = LvArray::math::log( trialComposition( trialIndex, ic ) );
+        normalizedComposition[ic] = composition[ic]*(alpha * kValues[ic] + (1.0 - alpha) / kValues[ic]);
+        logTrialComposition[ic] = LvArray::math::log( normalizedComposition[ic] );
       }
+
       for( localIndex iterationCount = 0; iterationCount < MultiFluidConstants::maxSSIIterations; ++iterationCount )
       {
-        for( integer const ic : presentComponents )
-        {
-          logTrialComposition[ic] = hyperplane[ic] - logFugacity[ic];
-          trialComposition( trialIndex, ic ) = LvArray::math::exp( logTrialComposition[ic] );
-          normalizedComposition[ic] = trialComposition( trialIndex, ic );
-        }
-        normalizeComposition( numComps, normalizedComposition.toSlice() );
-
+        // Normalise the composition and calculate the fugacity
+        real64 const totalMoles = normalizeComposition( numComps, normalizedComposition );
         FugacityCalculator::computeLogFugacity( numComps,
                                                 pressure,
                                                 temperature,
@@ -144,39 +120,47 @@ struct StabilityTest
                                                 equationOfState,
                                                 logFugacity );
 
+        // Calculate the TPD
+        real64 tpd = 0.0;
+        for( integer const ic : presentComponents )
+        {
+          tpd += composition[ic] + totalMoles * normalizedComposition[ic] * (logTrialComposition[ic] + logFugacity[ic] - hyperplane[ic] - 1.0);
+        }
+        if( tpd < tangentPlaneDistance )
+        {
+          tangentPlaneDistance = tpd;
+        }
+        if( tangentPlaneDistance < -MultiFluidConstants::fugacityTolerance )
+        {
+          break;
+        }
+
+        // Check stationarity
         real64 error = 0.0;
         for( integer const ic : presentComponents )
         {
           real64 const dG =  logTrialComposition[ic] + logFugacity[ic] - hyperplane[ic];
           error += (dG*dG);
         }
         error = LvArray::math::sqrt( error );
-
         if( error < MultiFluidConstants::fugacityTolerance )
         {
-          // Calculate modified tangent plane distance (Michelsen, 1982b) of trial composition relative to input composition
-          real64 tpd = 1.0;
-          for( integer const ic : presentComponents )
-          {
-            tpd += trialComposition( trialIndex, ic ) * (logTrialComposition[ic] + logFugacity[ic] - hyperplane[ic] - 1.0);
-          }
-          if( tpd < tangentPlaneDistance )
-          {
-            tangentPlaneDistance = tpd;
-          }
-          numberOfStationaryPoints++;
           break;
         }
+
+        // Update to next step
+        for( integer const ic : presentComponents )
+        {
+          logTrialComposition[ic] = hyperplane[ic] - logFugacity[ic];
+          normalizedComposition[ic] = LvArray::math::exp( logTrialComposition[ic] );
+        }
       }
-    }
-    if( numberOfStationaryPoints == numTrials )
-    {
-      for( integer const ic : presentComponents )
+      if( tangentPlaneDistance < -MultiFluidConstants::fugacityTolerance )
       {
-        kValues[ic] = trialComposition( 1, ic ) / trialComposition( 0, ic );
+        break;
       }
     }
-    return numberOfStationaryPoints == numTrials;
+    return true;
   }
 
 private: