Skip to content

Commit

Permalink
Compositional multiphase fluid properties (#2542)
Browse files Browse the repository at this point in the history
Creates a utility for the calculation of compositional multiphase properties
- Molar density
- Mass density
- Adds unit tests to compare against values generated by PVTPackage
  • Loading branch information
dkachuma authored Sep 20, 2023
1 parent 9001157 commit c4052de
Show file tree
Hide file tree
Showing 6 changed files with 826 additions and 4 deletions.
6 changes: 4 additions & 2 deletions src/coreComponents/constitutive/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -50,10 +50,11 @@ set( constitutive_headers
fluid/multifluid/CO2Brine/functions/CO2EOSSolver.hpp
fluid/multifluid/CO2Brine/functions/PureWaterProperties.hpp
fluid/multifluid/CO2Brine/functions/WaterDensity.hpp
fluid/multifluid/compositional/functions/CubicEOSPhaseModel.hpp
fluid/multifluid/compositional/functions/CompositionalProperties.hpp
fluid/multifluid/compositional/functions/CubicEOSPhaseModel.hpp
fluid/multifluid/compositional/functions/KValueInitialization.hpp
fluid/multifluid/compositional/functions/NegativeTwoPhaseFlash.hpp
fluid/multifluid/compositional/functions/RachfordRice.hpp
fluid/multifluid/compositional/functions/RachfordRice.hpp
fluid/multifluid/reactive/ReactiveBrineFluid.hpp
fluid/multifluid/reactive/ReactiveMultiFluid.hpp
fluid/multifluid/reactive/ReactiveMultiFluidFields.hpp
Expand Down Expand Up @@ -183,6 +184,7 @@ set( constitutive_sources
fluid/multifluid/CO2Brine/functions/CO2EOSSolver.cpp
fluid/multifluid/CO2Brine/functions/PureWaterProperties.cpp
fluid/multifluid/CO2Brine/functions/WaterDensity.cpp
fluid/multifluid/compositional/functions/CompositionalProperties.cpp
fluid/multifluid/reactive/ReactiveBrineFluid.cpp
fluid/multifluid/reactive/ReactiveMultiFluid.cpp
fluid/multifluid/reactive/ReactiveFluidDriver.cpp
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,146 @@
/*
* ------------------------------------------------------------------------------------------------------------
* SPDX-License-Identifier: LGPL-2.1-only
*
* Copyright (c) 2018-2020 Lawrence Livermore National Security LLC
* Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University
* Copyright (c) 2018-2020 TotalEnergies
* Copyright (c) 2019- GEOSX Contributors
* All rights reserved
*
* See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
* ------------------------------------------------------------------------------------------------------------
*/

/**
* @file CompositionalProperties.hpp
*/

#include "CompositionalProperties.hpp"

namespace geos
{

namespace constitutive
{

/*
* Calculate the molar volume and apply the Peneloux shift parameters. The parameters should be in
* dimensional form.
* Peneloux, A et al. 1982. Fluid phase equilibria, 8(1):7–23.
* https://doi.org/10.1016/0378-3812(82)80002-2
*/
GEOS_HOST_DEVICE
void CompositionalProperties::computeMolarDensity( integer const numComps,
real64 const pressure,
real64 const temperature,
arraySlice1d< real64 const > const & composition,
arraySlice1d< real64 const > const & volumeShift,
real64 const compressibilityFactor,
real64 & molarDensity )
{

real64 vEos = gasConstant * temperature * compressibilityFactor / pressure;
real64 vCorrected = vEos;

for( integer ic = 0; ic < numComps; ++ic )
{
vCorrected += composition[ic] * volumeShift[ic];
}

if( epsilon < vCorrected )
{
molarDensity = 1.0 / vCorrected;
}
else
{
molarDensity = 0.0;
}
}

GEOS_HOST_DEVICE
void CompositionalProperties::computeMolarDensity( integer const numComps,
real64 const pressure,
real64 const temperature,
arraySlice1d< real64 const > const & GEOS_UNUSED_PARAM ( composition ),
arraySlice1d< real64 const > const & volumeShift,
real64 const compressibilityFactor,
real64 const dCompressibilityFactor_dp,
real64 const dCompressibilityFactor_dt,
arraySlice1d< real64 const > const & dCompressibilityFactor_dz,
real64 const molarDensity,
real64 & dMolarDensity_dp,
real64 & dMolarDensity_dt,
arraySlice1d< real64 > const & dMolarDensity_dz )
{
if( molarDensity < epsilon )
{
dMolarDensity_dp = 0.0;
dMolarDensity_dt = 0.0;
for( integer ic = 0; ic < numComps; ++ic )
{
dMolarDensity_dz[ic] = 0.0;
}
return;
}

real64 dvCorrected_dx = 0.0;

// Pressure derivative
dvCorrected_dx = gasConstant * temperature * (dCompressibilityFactor_dp - compressibilityFactor / pressure) / pressure;
dMolarDensity_dp = -molarDensity * molarDensity * dvCorrected_dx;

// Temperature derivative
dvCorrected_dx = gasConstant * (temperature * dCompressibilityFactor_dt + compressibilityFactor) / pressure;
dMolarDensity_dt = -molarDensity * molarDensity * dvCorrected_dx;

// Composition derivative
for( integer ic = 0; ic < numComps; ++ic )
{
dvCorrected_dx = gasConstant * temperature * dCompressibilityFactor_dz[ic] / pressure + volumeShift[ic];
dMolarDensity_dz[ic] = -molarDensity * molarDensity * dvCorrected_dx;
}
}

GEOS_HOST_DEVICE
void CompositionalProperties::computeMassDensity( integer const numComps,
arraySlice1d< real64 const > const & composition,
arraySlice1d< real64 const > const & molecularWeight,
real64 const molarDensity,
real64 & massDensity )
{
massDensity = 0.0;
for( integer ic = 0; ic < numComps; ++ic )
{
massDensity += molecularWeight[ic] * composition[ic] * molarDensity;
}
}

GEOS_HOST_DEVICE
void CompositionalProperties::computeMassDensity( integer const numComps,
arraySlice1d< real64 const > const & molecularWeight,
real64 const molarDensity,
real64 const dMolarDensity_dp,
real64 const dMolarDensity_dt,
arraySlice1d< real64 const > const dMolarDensity_dz,
real64 const massDensity,
real64 & dMassDensity_dp,
real64 & dMassDensity_dt,
arraySlice1d< real64 > const & dMassDensity_dz )
{
// Pressure derivative
dMassDensity_dp = massDensity * dMolarDensity_dp / molarDensity;

// Temperature derivative
dMassDensity_dt = massDensity * dMolarDensity_dt / molarDensity;

// Composition derivative
for( integer ic = 0; ic < numComps; ++ic )
{
dMassDensity_dz[ic] = massDensity * dMolarDensity_dz[ic] / molarDensity + molecularWeight[ic] * molarDensity;
}
}

} // namespace constitutive

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

/**
* @file CompositionalProperties.hpp
*/

#ifndef GEOS_CONSTITUTIVE_FLUID_MULTIFLUID_COMPOSITIONAL_FUNCTIONS_COMPOSITIONALPROPERTIES_HPP_
#define GEOS_CONSTITUTIVE_FLUID_MULTIFLUID_COMPOSITIONAL_FUNCTIONS_COMPOSITIONALPROPERTIES_HPP_

#include "common/DataTypes.hpp"

namespace geos
{

namespace constitutive
{

struct CompositionalProperties
{
public:
/// Epsilon used in the calculations
static constexpr real64 epsilon = LvArray::NumericLimits< real64 >::epsilon;
/// Universal gas constant
static constexpr real64 gasConstant = 8.31446261815324;

/**
* @brief Compute the molar density of a mixture from the composition and the compressibility factor
* @param[in] numComps number of components
* @param[in] pressure pressure
* @param[in] temperature temperature
* @param[in] composition composition of the mixture
* @param[in] volumeShift dimensional volume shift parameters
* @param[in] compressibilityFactor compressibility factor (z-factor)
* @param[out] molarDensity the calculated molar density
* @note The volume shifts can result in a negative molar density which will be truncated to zero
*/
GEOS_HOST_DEVICE
static void computeMolarDensity( integer const numComps,
real64 const pressure,
real64 const temperature,
arraySlice1d< real64 const > const & composition,
arraySlice1d< real64 const > const & volumeShift,
real64 const compressibilityFactor,
real64 & molarDensity );

/**
* @brief Compute the molar density derivatives
* @param[in] numComps number of components
* @param[in] pressure pressure
* @param[in] temperature temperature
* @param[in] composition composition of the mixture
* @param[in] volumeShift dimensional volume shift parameters
* @param[in] compressibilityFactor compressibility factor (z-factor)
* @param[in] dCompressibilityFactor_dp derivative of the compressibility factor (z-factor) wrt pressure
* @param[in] dCompressibilityFactor_dp derivative of the compressibility factor (z-factor) wrt temperature
* @param[in] dCompressibilityFactor_dz derivative of the compressibility factor (z-factor) wrt composition
* @param[in] molarDensity the calculated molar density
* @param[out] dMolarDensity_dp derivative of the molar density wrt pressure
* @param[out] dMolarDensity_dt derivative of the molar density wrt temperature
* @param[out] dMolarDensity_dz derivative of the molar density wrt composition
*/
GEOS_HOST_DEVICE
static void computeMolarDensity( integer const numComps,
real64 const pressure,
real64 const temperature,
arraySlice1d< real64 const > const & composition,
arraySlice1d< real64 const > const & volumeShift,
real64 const compressibilityFactor,
real64 const dCompressibilityFactor_dp,
real64 const dCompressibilityFactor_dt,
arraySlice1d< real64 const > const & dCompressibilityFactor_dz,
real64 const molarDensity,
real64 & dMolarDensity_dp,
real64 & dMolarDensity_dt,
arraySlice1d< real64 > const & dMolarDensity_dz );

/**
* @brief Compute the mass density of a mixture from the composition and the molar density
* @param[in] numComps number of components
* @param[in] composition composition of the mixture
* @param[in] molecularWeight the component molecular weights
* @param[in] molarDensity the mixture molar density
* @param[out] massDensity the calculated mass density
*/
GEOS_HOST_DEVICE
static void computeMassDensity( integer const numComps,
arraySlice1d< real64 const > const & composition,
arraySlice1d< real64 const > const & molecularWeight,
real64 const molarDensity,
real64 & massDensity );

/**
* @brief Compute the mass density derivatives
* @param[in] numComps number of components
* @param[in] molecularWeight the component molecular weights
* @param[in] molarDensity the mixture molar density
* @param[in] dMolarDensity_dp derivative of the molar density wrt pressure
* @param[in] dMolarDensity_dt derivative of the molar density wrt temperature
* @param[in] dMolarDensity_dz derivative of the molar density wrt composition
* @param[in] massDensity mass density
* @param[out] dMassDensity_dp derivative of the mass density wrt pressure
* @param[out] dMassDensity_dt derivative of the mass density wrt temperature
* @param[out] dMassDensity_dz derivative of the mass density wrt composition
*/
GEOS_HOST_DEVICE
static void computeMassDensity( integer const numComps,
arraySlice1d< real64 const > const & molecularWeight,
real64 const molarDensity,
real64 const dMolarDensity_dp,
real64 const dMolarDensity_dt,
arraySlice1d< real64 const > const dMolarDensity_dz,
real64 const massDensity,
real64 & dMassDensity_dp,
real64 & dMassDensity_dt,
arraySlice1d< real64 > const & dMassDensity_dz );
};

} // namespace constitutive

} // namespace geos

#endif //GEOS_CONSTITUTIVE_FLUID_MULTIFLUID_COMPOSITIONAL_FUNCTIONS_COMPOSITIONALPROPERTIES_HPP_
1 change: 1 addition & 0 deletions src/coreComponents/constitutive/unitTests/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
#

set( gtest_geosx_tests
testCompositionalProperties.cpp
testDamageUtilities.cpp
testDruckerPrager.cpp
testElasticIsotropic.cpp
Expand Down
9 changes: 7 additions & 2 deletions src/coreComponents/constitutive/unitTests/TestFluid.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -47,8 +47,9 @@ struct Fluid
static constexpr integer Tc = 1;
static constexpr integer Ac = 2;
static constexpr integer Mw = 3;
static constexpr integer Vs = 4;

static std::array< real64, 44 > data;
static std::array< real64, 55 > data;
};

template< int NC >
Expand All @@ -67,6 +68,7 @@ class TestFluid
createArray( testFluid->criticalTemperature, components, Fluid::Tc, Fluid::data );
createArray( testFluid->acentricFactor, components, Fluid::Ac, Fluid::data );
createArray( testFluid->molecularWeight, components, Fluid::Mw, Fluid::data );
createArray( testFluid->volumeShift, components, Fluid::Vs, Fluid::data );
return testFluid;
}

Expand Down Expand Up @@ -107,7 +109,7 @@ class TestFluid
}
};

std::array< real64, 44 > Fluid::data = {
std::array< real64, 55 > Fluid::data = {
// -- Pc
2.2050e+07, 7.3750e+06, 3.4000e+06, 8.9630e+06, 1.2960e+06, 4.8721e+06,
4.2481e+06, 3.6400e+06, 4.5990e+06, 2.5300e+06, 1.4600e+06,
Expand All @@ -120,6 +122,9 @@ std::array< real64, 44 > Fluid::data = {
// -- Mw
1.8015e+01, 4.4010e+01, 2.8013e+01, 3.4100e+01, 1.6043e+01, 3.0070e+01,
4.4097e+01, 5.8124e+01, 7.2151e+01, 1.1423e+02, 1.4228e+02,
// -- Vs
0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
};

}// testing
Expand Down
Loading

0 comments on commit c4052de

Please sign in to comment.