Skip to content

Commit

Permalink
feat: Add immiscible water flash model (#3237)
Browse files Browse the repository at this point in the history
* Add immiscible water flash model
  • Loading branch information
dkachuma authored and rrsettgast committed Sep 17, 2024
1 parent 917c9e2 commit d04c70b
Show file tree
Hide file tree
Showing 11 changed files with 862 additions and 10 deletions.
4 changes: 4 additions & 0 deletions src/coreComponents/constitutive/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -77,6 +77,8 @@ set( constitutive_headers
fluid/multifluid/compositional/models/CriticalVolume.hpp
fluid/multifluid/compositional/models/EquationOfState.hpp
fluid/multifluid/compositional/models/FunctionBase.hpp
fluid/multifluid/compositional/models/ImmiscibleWaterFlashModel.hpp
fluid/multifluid/compositional/models/ImmiscibleWaterParameters.hpp
fluid/multifluid/compositional/models/LohrenzBrayClarkViscosity.hpp
fluid/multifluid/compositional/models/LohrenzBrayClarkViscosityImpl.hpp
fluid/multifluid/compositional/models/NegativeTwoPhaseFlashModel.hpp
Expand Down Expand Up @@ -222,6 +224,8 @@ set( constitutive_sources
fluid/multifluid/compositional/models/CompositionalDensity.cpp
fluid/multifluid/compositional/models/ConstantViscosity.cpp
fluid/multifluid/compositional/models/CriticalVolume.cpp
fluid/multifluid/compositional/models/ImmiscibleWaterFlashModel.cpp
fluid/multifluid/compositional/models/ImmiscibleWaterParameters.cpp
fluid/multifluid/compositional/models/LohrenzBrayClarkViscosity.cpp
fluid/multifluid/compositional/models/NegativeTwoPhaseFlashModel.cpp
fluid/multifluid/compositional/CompositionalMultiphaseFluid.cpp
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -174,7 +174,7 @@ struct NegativeTwoPhaseFlash
* @param[out] fugacityRatios the fugacity rations
* @return The error
*/
template< integer USD >
template< integer USD1, integer USD2 >
GEOS_HOST_DEVICE
static real64 computeFugacityRatio(
integer const numComps,
Expand All @@ -184,11 +184,11 @@ struct NegativeTwoPhaseFlash
ComponentProperties::KernelWrapper const & componentProperties,
EquationOfStateType const liquidEos,
EquationOfStateType const vapourEos,
arraySlice1d< real64 const, USD > const & kValues,
arraySlice1d< real64 const, USD1 > const & kValues,
arraySlice1d< integer const > const & presentComponents,
real64 & vapourPhaseMoleFraction,
arraySlice1d< real64, USD > const & liquidComposition,
arraySlice1d< real64, USD > const & vapourComposition,
arraySlice1d< real64, USD2 > const & liquidComposition,
arraySlice1d< real64, USD2 > const & vapourComposition,
arraySlice1d< real64 > const & logLiquidFugacity,
arraySlice1d< real64 > const & logVapourFugacity,
arraySlice1d< real64 > const & fugacityRatios );
Expand Down Expand Up @@ -487,7 +487,7 @@ void NegativeTwoPhaseFlash::computeDerivatives(
}
}

template< integer USD >
template< integer USD1, integer USD2 >
GEOS_HOST_DEVICE
real64 NegativeTwoPhaseFlash::computeFugacityRatio(
integer const numComps,
Expand All @@ -497,11 +497,11 @@ real64 NegativeTwoPhaseFlash::computeFugacityRatio(
ComponentProperties::KernelWrapper const & componentProperties,
EquationOfStateType const liquidEos,
EquationOfStateType const vapourEos,
arraySlice1d< real64 const, USD > const & kValues,
arraySlice1d< real64 const, USD1 > const & kValues,
arraySlice1d< integer const > const & presentComponents,
real64 & vapourPhaseMoleFraction,
arraySlice1d< real64, USD > const & liquidComposition,
arraySlice1d< real64, USD > const & vapourComposition,
arraySlice1d< real64, USD2 > const & liquidComposition,
arraySlice1d< real64, USD2 > const & vapourComposition,
arraySlice1d< real64 > const & logLiquidFugacity,
arraySlice1d< real64 > const & logVapourFugacity,
arraySlice1d< real64 > const & fugacityRatios )
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,7 @@ class ComponentProperties final
/**
* Data accessors
*/
arrayView1d< string > const & getComponentName() const { return m_componentNames; }
arrayView1d< real64 > const & getComponentMolarWeight() const { return m_componentMolarWeight; }
arrayView1d< real64 > const & getComponentCriticalPressure() const { return m_componentCriticalPressure; }
arrayView1d< real64 > const & getComponentCriticalTemperature() const { return m_componentCriticalTemperature; }
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,104 @@
/*
* ------------------------------------------------------------------------------------------------------------
* SPDX-License-Identifier: LGPL-2.1-only
*
* Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
* Copyright (c) 2018-2024 Total, S.A
* Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
* Copyright (c) 2018-2024 Chevron
* Copyright (c) 2019- GEOS/GEOSX Contributors
* All rights reserved
*
* See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
* ------------------------------------------------------------------------------------------------------------
*/

/**
* @file ImmiscibleWaterFlashModel.cpp
*/

#include "ImmiscibleWaterFlashModel.hpp"
#include "ImmiscibleWaterParameters.hpp"
#include "EquationOfState.hpp"

namespace geos
{

namespace constitutive
{

namespace compositional
{

// Naming conventions
string ImmiscibleWaterFlashModel::catalogName()
{
return "ThreePhase";
}

ImmiscibleWaterFlashModel::ImmiscibleWaterFlashModel( string const & name,
ComponentProperties const & componentProperties,
ModelParameters const & modelParameters ):
FunctionBase( name, componentProperties ),
m_parameters( modelParameters )
{
m_waterComponentIndex = ImmiscibleWaterParameters::getWaterComponentIndex( componentProperties );
}

ImmiscibleWaterFlashModel::KernelWrapper
ImmiscibleWaterFlashModel::createKernelWrapper() const
{
constexpr integer liquidIndex = 0;
constexpr integer vapourIndex = 1;
constexpr integer aqueousIndex = 2;
EquationOfState const * equationOfState = m_parameters.get< EquationOfState >();
EquationOfStateType const liquidEos = EnumStrings< EquationOfStateType >::fromString( equationOfState->m_equationsOfStateNames[liquidIndex] );
EquationOfStateType const vapourEos = EnumStrings< EquationOfStateType >::fromString( equationOfState->m_equationsOfStateNames[vapourIndex] );

array1d< real64 > componentCriticalVolume( m_componentProperties.getNumberOfComponents());

return KernelWrapper( m_componentProperties.getNumberOfComponents(),
liquidIndex,
vapourIndex,
aqueousIndex,
m_waterComponentIndex,
liquidEos,
vapourEos,
componentCriticalVolume );
}

ImmiscibleWaterFlashModelUpdate::ImmiscibleWaterFlashModelUpdate(
integer const numComponents,
integer const liquidIndex,
integer const vapourIndex,
integer const aqueousIndex,
integer const waterComponentIndex,
EquationOfStateType const liquidEos,
EquationOfStateType const vapourEos,
arrayView1d< real64 const > const componentCriticalVolume ):
m_twoPhaseModel( numComponents,
liquidIndex,
vapourIndex,
liquidEos,
vapourEos,
componentCriticalVolume ),
m_numComponents( numComponents ),
m_liquidIndex( liquidIndex ),
m_vapourIndex( vapourIndex ),
m_aquoesIndex( aqueousIndex ),
m_waterComponentIndex( waterComponentIndex )
{}

std::unique_ptr< ModelParameters >
ImmiscibleWaterFlashModel::createParameters( std::unique_ptr< ModelParameters > parameters )
{
auto params = NegativeTwoPhaseFlashModel::createParameters( std::move( parameters ) );
params = ImmiscibleWaterParameters::create( std::move( params ) );
return params;
}

} // end namespace compositional

} // namespace constitutive

} // end namespace geos
Original file line number Diff line number Diff line change
@@ -0,0 +1,212 @@
/*
* ------------------------------------------------------------------------------------------------------------
* SPDX-License-Identifier: LGPL-2.1-only
*
* Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
* Copyright (c) 2018-2024 Total, S.A
* Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
* Copyright (c) 2018-2024 Chevron
* Copyright (c) 2019- GEOS/GEOSX Contributors
* All rights reserved
*
* See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
* ------------------------------------------------------------------------------------------------------------
*/

/**
* @file ImmiscibleWaterFlashModel.hpp
*/

#ifndef GEOS_CONSTITUTIVE_FLUID_MULTIFLUID_COMPOSITIONAL_MODELS_IMMISCIBLEWATERFLASHMODEL_HPP_
#define GEOS_CONSTITUTIVE_FLUID_MULTIFLUID_COMPOSITIONAL_MODELS_IMMISCIBLEWATERFLASHMODEL_HPP_

#include "FunctionBase.hpp"
#include "EquationOfState.hpp"

#include "constitutive/fluid/multifluid/Layouts.hpp"
#include "constitutive/fluid/multifluid/MultiFluidUtils.hpp"
#include "NegativeTwoPhaseFlashModel.hpp"

namespace geos
{

namespace constitutive
{

namespace compositional
{

class ModelParameters;

class ImmiscibleWaterFlashModelUpdate final : public FunctionBaseUpdate
{
private:
static constexpr integer maxNumComps = MultiFluidConstants::MAX_NUM_COMPONENTS;
public:

using PhaseProp = NegativeTwoPhaseFlashModelUpdate::PhaseProp;
using PhaseComp = NegativeTwoPhaseFlashModelUpdate::PhaseComp;
using Deriv = multifluid::DerivativeOffset;

ImmiscibleWaterFlashModelUpdate( integer const numComponents,
integer const liquidIndex,
integer const vapourIndex,
integer const aqueousIndex,
integer const waterComponentIndex,
EquationOfStateType const liquidEos,
EquationOfStateType const vapourEos,
arrayView1d< real64 const > const componentCriticalVolume );

// Mark as a 3-phase flash
GEOS_HOST_DEVICE
static constexpr integer getNumberOfPhases() { return 3; }

template< int USD1, int USD2 >
GEOS_HOST_DEVICE
void compute( ComponentProperties::KernelWrapper const & componentProperties,
real64 const & pressure,
real64 const & temperature,
arraySlice1d< real64 const, USD1 > const & compFraction,
arraySlice2d< real64, USD2 > const & kValues,
PhaseProp::SliceType const phaseFraction,
PhaseComp::SliceType const phaseCompFraction ) const;

private:
template< int USD >
GEOS_FORCE_INLINE
GEOS_HOST_DEVICE
void convertCompositionDerivatives( real64 const hcMoleFraction,
arraySlice1d< real64 const > const & composition,
arraySlice1d< real64, USD > const & derivatives ) const
{
real64 dvdzi = 0.0;
for( integer ic = 0; ic < m_numComponents; ++ic )
{
dvdzi += derivatives[Deriv::dC+ic] * composition[ic];
}
for( integer ic = 0; ic < m_numComponents; ++ic )
{
derivatives[Deriv::dC+ic] /= hcMoleFraction;
}
derivatives[Deriv::dC+m_waterComponentIndex] = dvdzi / hcMoleFraction;
}

private:
NegativeTwoPhaseFlashModel::KernelWrapper const m_twoPhaseModel;
integer const m_numComponents;
integer const m_liquidIndex;
integer const m_vapourIndex;
integer const m_aquoesIndex;
integer const m_waterComponentIndex;
};

class ImmiscibleWaterFlashModel : public FunctionBase
{
public:
ImmiscibleWaterFlashModel( string const & name,
ComponentProperties const & componentProperties,
ModelParameters const & modelParameters );

static string catalogName();

FunctionType functionType() const override
{
return FunctionType::FLASH;
}

/// Type of kernel wrapper for in-kernel update
using KernelWrapper = ImmiscibleWaterFlashModelUpdate;

/**
* @brief Create an update kernel wrapper.
* @return the wrapper
*/
KernelWrapper createKernelWrapper() const;

// Create parameters unique to this model
static std::unique_ptr< ModelParameters > createParameters( std::unique_ptr< ModelParameters > parameters );

private:
ModelParameters const & m_parameters;
integer m_waterComponentIndex{-1};
};

template< int USD1, int USD2 >
GEOS_HOST_DEVICE
void ImmiscibleWaterFlashModelUpdate::compute( ComponentProperties::KernelWrapper const & componentProperties,
real64 const & pressure,
real64 const & temperature,
arraySlice1d< real64 const, USD1 > const & compFraction,
arraySlice2d< real64, USD2 > const & kValues,
PhaseProp::SliceType const phaseFraction,
PhaseComp::SliceType const phaseCompFraction ) const
{
LvArray::forValuesInSlice( phaseFraction.value, setZero );
LvArray::forValuesInSlice( phaseFraction.derivs, setZero );
LvArray::forValuesInSlice( phaseCompFraction.value, setZero );
LvArray::forValuesInSlice( phaseCompFraction.derivs, setZero );

// Water phase
phaseFraction.value[m_aquoesIndex] = compFraction[m_waterComponentIndex];
phaseFraction.derivs( m_aquoesIndex, Deriv::dC + m_waterComponentIndex ) = 1.0;
phaseCompFraction.value( m_aquoesIndex, m_waterComponentIndex ) = 1.0;

// Total hydrocarbon mole fraction
real64 const z_hc = 1.0 - compFraction[m_waterComponentIndex];

if( z_hc < MultiFluidConstants::minForSpeciesPresence )
{
// Single phase water
real64 const constantComposition = 1.0 / (m_numComponents - 1);
for( integer ic = 0; ic < m_numComponents; ++ic )
{
phaseCompFraction.value( m_liquidIndex, ic ) = constantComposition;
phaseCompFraction.value( m_vapourIndex, ic ) = constantComposition;
}
phaseCompFraction.value( m_liquidIndex, m_waterComponentIndex ) = 0.0;
phaseCompFraction.value( m_vapourIndex, m_waterComponentIndex ) = 0.0;
}
else
{
// Hydrocarbon phases

// Calculate normalised hyrdocarbon composition
stackArray1d< real64, maxNumComps > composition( m_numComponents );
for( integer ic = 0; ic < m_numComponents; ++ic )
{
composition[ic] = compFraction[ic] / z_hc;
}
composition[m_waterComponentIndex] = 0.0;

// Perform negative two-phase flash
m_twoPhaseModel.compute( componentProperties,
pressure,
temperature,
composition.toSliceConst(),
kValues,
phaseFraction,
phaseCompFraction );

for( integer const phaseIndex : {m_liquidIndex, m_vapourIndex} )
{
real64 const v = phaseFraction.value[phaseIndex];
phaseFraction.value[phaseIndex] *= z_hc;
LvArray::forValuesInSlice( phaseFraction.derivs[phaseIndex], [&]( real64 & a ){ a *= z_hc; } );
convertCompositionDerivatives( z_hc, composition.toSliceConst(), phaseFraction.derivs[phaseIndex] );
phaseFraction.derivs( phaseIndex, Deriv::dC+m_waterComponentIndex ) = -v;

for( integer ic = 0; ic < m_numComponents; ++ic )
{
convertCompositionDerivatives( z_hc, composition.toSliceConst(), phaseCompFraction.derivs[phaseIndex][ic] );
}
}
}
}

} // end namespace compositional

} // end namespace constitutive

} // end namespace geos

#endif //GEOS_CONSTITUTIVE_FLUID_MULTIFLUID_COMPOSITIONAL_MODELS_IMMISCIBLEWATERFLASHMODEL_HPP_
Loading

0 comments on commit d04c70b

Please sign in to comment.