diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index 90caa852c88..c8eb7f9cf11 100644 --- a/BASELINE_NOTES.md +++ b/BASELINE_NOTES.md @@ -5,28 +5,25 @@ This file is designed to track changes to the integrated test baselines. Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining. These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD). +PR #4080 (2026-11-06) +Add thermal compositional models. + PR #4139 (2026-10-31) -===================== Add the mixed mimetic discretization: new `NumericalMethods/MixedMimetic` group in every restart file and new mixed MFD tests. PR #4169 (2026-10-30) -===================== Fix surface generator splitting for contact mechanics case with wells PR #4165 (2026-10-29) -===================== Add zero salinity, non-default compressibility pvt driver test PR #4164 (2026-10-26) -===================== Add a `Group` statistics layer for fracture surfaces. PR #4005 (2026-10-23) -===================== Added a `Group` statistics layer for sets statistics PR #4155 (2026-10-22) -===================== Add the missing restart baselines for the injection, production, and combined thermal compositional-well tests on one and four MPI ranks. In CI run 17777, all six simulations completed successfully, but their baseline checks failed diff --git a/inputFiles/compositionalMultiphaseFlow/benchmarks/jouleThomsonCooling/jouleThomsonCooling.ats b/inputFiles/compositionalMultiphaseFlow/benchmarks/jouleThomsonCooling/jouleThomsonCooling.ats new file mode 100644 index 00000000000..124f0388ecb --- /dev/null +++ b/inputFiles/compositionalMultiphaseFlow/benchmarks/jouleThomsonCooling/jouleThomsonCooling.ats @@ -0,0 +1,17 @@ +from geos.ats.test_builder import TestDeck, RestartcheckParameters, generate_geos_tests + +restartcheck_params = {} +restartcheck_params["atol"] = 1.0E-6 +restartcheck_params["rtol"] = 1.0E-5 + +decks = [ + TestDeck( + name="jouleThomsonCooling", + description="Joule-Thomson Cooling (thermal, 2-phase, cooling, radial)", + partitions=((1, 1, 1),), + restart_step=243, + check_step=320, + restartcheck_params=RestartcheckParameters(**restartcheck_params)), +] + +generate_geos_tests(decks) diff --git a/inputFiles/compositionalMultiphaseFlow/benchmarks/jouleThomsonCooling/jouleThomsonCooling.xml b/inputFiles/compositionalMultiphaseFlow/benchmarks/jouleThomsonCooling/jouleThomsonCooling.xml new file mode 100644 index 00000000000..b03182bfe8d --- /dev/null +++ b/inputFiles/compositionalMultiphaseFlow/benchmarks/jouleThomsonCooling/jouleThomsonCooling.xml @@ -0,0 +1,413 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/constitutiveDriver/constitutiveDriver.ats b/inputFiles/constitutiveDriver/constitutiveDriver.ats index 7c86d2bb31b..bf07cf8e15b 100644 --- a/inputFiles/constitutiveDriver/constitutiveDriver.ats +++ b/inputFiles/constitutiveDriver/constitutiveDriver.ats @@ -38,6 +38,7 @@ decks = [ create_pvt_test("testPVT_CO2BrineTables", "PVT Driver test for CO2 fluid models with solubility tables"), create_pvt_test("testPVT_PhaseComposition", "PVT Driver test with composition output"), create_pvt_test("testPVT_ThreePhaseCompositional", "PVT Driver test for three phase compositional models"), + create_pvt_test("testPVT_ThermalCompositional", "PVT Driver test for thermal compositional models"), create_pvt_test("testPVT_docExample", "PVT Driver documentation example"), create_relperm_test("testRelperm_Table_2_phase", "Relperm driver for 2 phase tables"), create_relperm_test("testRelperm_Table_3_phase", "Relperm driver for 3 phase tables"), diff --git a/inputFiles/constitutiveDriver/testPVT_ThermalCompositional.xml b/inputFiles/constitutiveDriver/testPVT_ThermalCompositional.xml new file mode 100644 index 00000000000..61baf42a2a1 --- /dev/null +++ b/inputFiles/constitutiveDriver/testPVT_ThermalCompositional.xml @@ -0,0 +1,131 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/constitutiveDriver/testPVT_data/kv_CO2.txt b/inputFiles/constitutiveDriver/testPVT_data/kv_CO2.txt new file mode 100644 index 00000000000..18b38db1084 --- /dev/null +++ b/inputFiles/constitutiveDriver/testPVT_data/kv_CO2.txt @@ -0,0 +1,65 @@ +2673.58 +1569.54 +921.641 +541.598 +318.844 +188.636 +113.161 +71.3553 +68.6664 +67.55 +64.5221 +60.3181 +55.2298 +6307.33 +3705.66 +2178.5 +1282.45 +756.968 +449.092 +269.132 +165.082 +107.296 +80.6751 +70.9355 +63.0822 +55.6087 +6889.77 +4048.08 +2379.8 +1400.98 +826.731 +490.356 +293.557 +179.514 +115.48 +84.0046 +71.4461 +62.7103 +54.8532 +7167.72 +4214.56 +2477.85 +1458.53 +860.685 +510.375 +305.381 +186.466 +119.421 +85.6759 +71.6673 +62.4604 +54.4108 +7695.74 +4532.4 +2666.22 +1569.45 +925.991 +548.848 +328.07 +199.751 +126.891 +88.9229 +72.0417 +61.8357 +53.4057 diff --git a/inputFiles/constitutiveDriver/testPVT_data/kv_H2O.txt b/inputFiles/constitutiveDriver/testPVT_data/kv_H2O.txt new file mode 100644 index 00000000000..ce3c58302b3 --- /dev/null +++ b/inputFiles/constitutiveDriver/testPVT_data/kv_H2O.txt @@ -0,0 +1,65 @@ +0.0098474 +0.00582094 +0.00345813 +0.00207293 +0.00126211 +0.000790421 +0.000521592 +0.000384888 +0.000793974 +0.000966344 +0.00107505 +0.00115461 +0.0012034 +0.188589 +0.111214 +0.0658019 +0.0391599 +0.0235417 +0.0144074 +0.00909634 +0.00608502 +0.00457049 +0.00447163 +0.00551906 +0.00617058 +0.00644326 +0.271182 +0.160069 +0.0946648 +0.0562932 +0.0337951 +0.0206318 +0.012967 +0.00859759 +0.00632712 +0.00581969 +0.00685352 +0.00764637 +0.00796688 +0.322712 +0.190546 +0.112675 +0.0669761 +0.0401834 +0.0245026 +0.0153674 +0.010148 +0.00740101 +0.00663886 +0.00762397 +0.00848277 +0.00882578 +0.447498 +0.265667 +0.157324 +0.0934601 +0.0560057 +0.0340784 +0.0212896 +0.0139536 +0.0100132 +0.00861264 +0.00940793 +0.010377 +0.0107577 diff --git a/inputFiles/constitutiveDriver/testPVT_data/kv_pressure.txt b/inputFiles/constitutiveDriver/testPVT_data/kv_pressure.txt new file mode 100644 index 00000000000..9b2759b2651 --- /dev/null +++ b/inputFiles/constitutiveDriver/testPVT_data/kv_pressure.txt @@ -0,0 +1,13 @@ +100000 +170400 +290400 +494900 +843400 +1.437e+06 +2.449e+06 +4.174e+06 +7.114e+06 +1.212e+07 +2.066e+07 +3.521e+07 +6e+07 diff --git a/inputFiles/constitutiveDriver/testPVT_data/kv_temperature.txt b/inputFiles/constitutiveDriver/testPVT_data/kv_temperature.txt new file mode 100644 index 00000000000..3c7730afdf9 --- /dev/null +++ b/inputFiles/constitutiveDriver/testPVT_data/kv_temperature.txt @@ -0,0 +1,5 @@ +283.15 +334.15 +342.15 +346.15 +354.15 diff --git a/src/coreComponents/constitutive/docs/compositional/Enthalpy.rst b/src/coreComponents/constitutive/docs/compositional/Enthalpy.rst index a14f4f35a79..0a8910c41c0 100644 --- a/src/coreComponents/constitutive/docs/compositional/Enthalpy.rst +++ b/src/coreComponents/constitutive/docs/compositional/Enthalpy.rst @@ -21,6 +21,9 @@ The ideal enthalpy represents the enthalpy of the mixture assuming it behaves as .. math:: C_{p,i}(T) = a_{i,0} + a_{i,1}\Delta T + a_{i,2}\Delta T^2 + a_{i,3}\Delta T^3 + a_{i,4}\Delta T^4 +.. note:: + The component coefficients :math:`a_{i,k}` are validated to ensure that :math:`C_{p,i}(T)>0` for T in the range :math:`[283.15, 473.15]` K (:math:`[10, 200]` °C). + To find the ideal enthalpy for a component, this heat capacity is integrated with respect to temperature from the reference temperature to the current temperature, and then added to a baseline reference enthalpy :math:`H_{ref,i}`: .. math:: diff --git a/src/coreComponents/constitutive/docs/compositional/Models.rst b/src/coreComponents/constitutive/docs/compositional/Models.rst index a7e01ac1c29..f0021abd335 100644 --- a/src/coreComponents/constitutive/docs/compositional/Models.rst +++ b/src/coreComponents/constitutive/docs/compositional/Models.rst @@ -54,4 +54,27 @@ The following table outlines the available compositional fluid models (specified * Possible applications: Saline aquifers and CO2 sequestration with stable compositions where pre-tabulated K-values are sufficient. * Phase split method: K-Value Flash. * Phase density models: Phillips brine density model; Compositional density. - * Phase viscosity models: Phillips brine viscosity model; Lohrenz-Bray-Clark viscosity. \ No newline at end of file + * Phase viscosity models: Phillips brine viscosity model; Lohrenz-Bray-Clark viscosity. + + * - ``CompositionalTwoPhaseThermalFluidLohrenzBrayClark`` :ref:`XML_CompositionalTwoPhaseThermalFluidLohrenzBrayClark` + + * Number of phases: 2. + * Possible applications: Non-isothermal miscible gas injection and thermal recovery processes involving drastic composition changes. + * Phase split method: Negative Two-Phase Flash. + * Phase density models: Compositional density. + * Phase viscosity models: Lohrenz-Bray-Clark viscosity. + + * - ``CompositionalTwoPhaseThermalFluidPhillipsBrine`` :ref:`XML_CompositionalTwoPhaseThermalFluidPhillipsBrine` + + * Number of phases: 2. + * Possible applications: Non-isothermal carbon sequestration (CO2-brine systems) or geothermal processes. + * Phase split method: Negative Two-Phase Flash. + * Phase density models: Phillips brine density model for the aqueous phase; Compositional density for the gas phase. + * Phase viscosity models: Phillips brine viscosity model for the aqueous phase; Lohrenz-Bray-Clark viscosity for the gas phase. + + * - ``CompositionalTwoPhaseKValueThermalFluidPhillipsBrine`` :ref:`XML_CompositionalTwoPhaseKValueThermalFluidPhillipsBrine` + + * Number of phases: 2. + * Possible applications: Non-isothermal saline aquifers and CO2 sequestration with stable compositions where pre-tabulated K-values are sufficient. + * Phase split method: K-Value Flash. + * Phase density models: Phillips brine density model; Compositional density. \ No newline at end of file diff --git a/src/coreComponents/constitutive/fluid/multifluid/MultiFluidBase.hpp b/src/coreComponents/constitutive/fluid/multifluid/MultiFluidBase.hpp index f84e49d1116..892a91eda6d 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/MultiFluidBase.hpp +++ b/src/coreComponents/constitutive/fluid/multifluid/MultiFluidBase.hpp @@ -1047,19 +1047,25 @@ MultiFluidBase::KernelWrapper:: { bool const phaseExists = (phaseFrac.value[ip] > 0); - if( !phaseExists ) + if( phaseExists ) { - continue; + real64 const densInv = 1.0 / phaseMassDens.value[ip]; + real64 const densInvSquared = densInv * densInv; + phaseInternalEnergy.value[ip] = phaseEnthalpy.value[ip] - pressure * densInv; + for( integer idof = 0; idof < numDOF; ++idof ) + { + phaseInternalEnergy.derivs[ip][idof] = phaseEnthalpy.derivs[ip][idof] + pressure * phaseMassDens.derivs[ip][idof] * densInvSquared; + } + phaseInternalEnergy.derivs[ip][multifluid::DerivativeOffset::dP] -= densInv; } - - real64 const densInv = 1.0 / phaseMassDens.value[ip]; - real64 const densInvSquared = densInv * densInv; - phaseInternalEnergy.value[ip] = phaseEnthalpy.value[ip] - pressure * densInv; - for( integer idof = 0; idof < numDOF; ++idof ) + else { - phaseInternalEnergy.derivs[ip][idof] = phaseEnthalpy.derivs[ip][idof] + pressure * phaseMassDens.derivs[ip][idof] * densInvSquared; + phaseInternalEnergy.value[ip] = 0.0; + for( integer idof = 0; idof < numDOF; ++idof ) + { + phaseInternalEnergy.derivs[ip][idof] = 0.0; + } } - phaseInternalEnergy.derivs[ip][multifluid::DerivativeOffset::dP] -= densInv; } } diff --git a/src/coreComponents/constitutive/fluid/multifluid/MultiFluidSelector.hpp b/src/coreComponents/constitutive/fluid/multifluid/MultiFluidSelector.hpp index cf809962701..c44e79cc6c9 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/MultiFluidSelector.hpp +++ b/src/coreComponents/constitutive/fluid/multifluid/MultiFluidSelector.hpp @@ -48,7 +48,10 @@ void constitutiveUpdatePassThru( constitutive::MultiFluidBase const & fluid, CompositionalTwoPhasePhillipsBrine, CompositionalKValuePhillipsBrine, CompositionalKValueLohrenzBrayClarkViscosity, - CompositionalTwoPhaseConstantViscosity + CompositionalTwoPhaseConstantViscosity, + CompositionalThermalTwoPhaseLohrenzBrayClarkViscosity, + CompositionalThermalTwoPhasePhillipsBrine, + CompositionalThermalKValuePhillipsBrine >::execute( fluid, std::forward< LAMBDA >( lambda ) ); } @@ -68,7 +71,10 @@ void constitutiveUpdatePassThru( constitutive::MultiFluidBase & fluid, CompositionalTwoPhasePhillipsBrine, CompositionalKValuePhillipsBrine, CompositionalKValueLohrenzBrayClarkViscosity, - CompositionalTwoPhaseConstantViscosity + CompositionalTwoPhaseConstantViscosity, + CompositionalThermalTwoPhaseLohrenzBrayClarkViscosity, + CompositionalThermalTwoPhasePhillipsBrine, + CompositionalThermalKValuePhillipsBrine >::execute( fluid, std::forward< LAMBDA >( lambda ) ); } diff --git a/src/coreComponents/constitutive/fluid/multifluid/compositional/CompositionalMultiphaseFluid.cpp b/src/coreComponents/constitutive/fluid/multifluid/compositional/CompositionalMultiphaseFluid.cpp index d8ab70a4a8f..0b88eaefaf5 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/compositional/CompositionalMultiphaseFluid.cpp +++ b/src/coreComponents/constitutive/fluid/multifluid/compositional/CompositionalMultiphaseFluid.cpp @@ -40,6 +40,7 @@ CompositionalMultiphaseFluid( string const & name, Group * const parent ) { using InputFlags = dataRepository::InputFlags; using RestartFlags = dataRepository::RestartFlags; + using PlotLevel = dataRepository::PlotLevel; getWrapperBase( viewKeyStruct::componentNamesString() ).setInputFlag( InputFlags::REQUIRED ); getWrapperBase( viewKeyStruct::componentMolarWeightString() ).setInputFlag( InputFlags::REQUIRED ); @@ -67,6 +68,17 @@ CompositionalMultiphaseFluid( string const & name, Group * const parent ) registerField< fields::multifluid::kValues >( &m_kValues ); + if constexpr ( isThermalType() ) + { + getField< fields::multifluid::phaseEnthalpy >(). + setPlotLevel( PlotLevel::LEVEL_0 ). + setRestartFlags( RestartFlags::WRITE_AND_READ ); + + getField< fields::multifluid::phaseInternalEnergy >(). + setPlotLevel( PlotLevel::LEVEL_0 ). + setRestartFlags( RestartFlags::WRITE_AND_READ ); + } + // Link parameters specific to each model m_parameters->registerParameters( this ); @@ -97,7 +109,9 @@ string CompositionalMultiphaseFluid< FLASH, PHASES... >::catalogName() // Use the first phase viscosity using FirstPhase = typename camp::at< camp::list< PHASES... >, camp::num< 0 > >::type; using ViscosityType = typename FirstPhase::Viscosity; - return GEOS_FMT( "Compositional{}Fluid{}", FLASH::catalogName(), ViscosityType::catalogName() ); + return GEOS_FMT( "Compositional{}{}Fluid{}", FLASH::catalogName(), + isThermalType() ? "Thermal": "", + ViscosityType::catalogName() ); } template< typename FLASH, typename ... PHASES > @@ -328,6 +342,18 @@ template class CompositionalMultiphaseFluid< compositional::KValueFlashModel< 2 >, compositional::PhaseModel< compositional::PhillipsBrineDensity, compositional::PhillipsBrineViscosity >, compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity > >; +template class CompositionalMultiphaseFluid< + compositional::NegativeTwoPhaseFlashModel, + compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity, compositional::CompositionalEnthalpy >, + compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity, compositional::CompositionalEnthalpy > >; +template class CompositionalMultiphaseFluid< + compositional::NegativeTwoPhaseFlashModel, + compositional::PhaseModel< compositional::PhillipsBrineDensity, compositional::PhillipsBrineViscosity, compositional::CompositionalEnthalpy >, + compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity, compositional::CompositionalEnthalpy > >; +template class CompositionalMultiphaseFluid< + compositional::KValueFlashModel< 2 >, + compositional::PhaseModel< compositional::PhillipsBrineDensity, compositional::PhillipsBrineViscosity, compositional::CompositionalEnthalpy >, + compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity, compositional::CompositionalEnthalpy > >; REGISTER_CATALOG_ENTRY( ConstitutiveBase, CompositionalTwoPhaseConstantViscosity, @@ -359,6 +385,21 @@ REGISTER_CATALOG_ENTRY( ConstitutiveBase, string const &, dataRepository::Group * const ) +REGISTER_CATALOG_ENTRY( ConstitutiveBase, + CompositionalThermalTwoPhaseLohrenzBrayClarkViscosity, + string const &, + dataRepository::Group * const ) + +REGISTER_CATALOG_ENTRY( ConstitutiveBase, + CompositionalThermalTwoPhasePhillipsBrine, + string const &, + dataRepository::Group * const ) + +REGISTER_CATALOG_ENTRY( ConstitutiveBase, + CompositionalThermalKValuePhillipsBrine, + string const &, + dataRepository::Group * const ) + } // namespace constitutive } // namespace geos diff --git a/src/coreComponents/constitutive/fluid/multifluid/compositional/CompositionalMultiphaseFluid.hpp b/src/coreComponents/constitutive/fluid/multifluid/compositional/CompositionalMultiphaseFluid.hpp index d116ad4ed73..cd85e5135d1 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/compositional/CompositionalMultiphaseFluid.hpp +++ b/src/coreComponents/constitutive/fluid/multifluid/compositional/CompositionalMultiphaseFluid.hpp @@ -23,6 +23,7 @@ #include "constitutive/fluid/multifluid/compositional/CompositionalMultiphaseFluidUpdates.hpp" #include "constitutive/fluid/multifluid/compositional/models/ConstantViscosity.hpp" #include "constitutive/fluid/multifluid/compositional/models/CompositionalDensity.hpp" +#include "constitutive/fluid/multifluid/compositional/models/CompositionalEnthalpy.hpp" #include "constitutive/fluid/multifluid/compositional/models/ImmiscibleWaterDensity.hpp" #include "constitutive/fluid/multifluid/compositional/models/ImmiscibleWaterFlashModel.hpp" #include "constitutive/fluid/multifluid/compositional/models/ImmiscibleWaterViscosity.hpp" @@ -88,6 +89,8 @@ class CompositionalMultiphaseFluid : public MultiFluidBase virtual void allocateConstitutiveData( dataRepository::Group & parent, localIndex const numPts ) override; + bool isThermal() const override { return isThermalType(); } + // TODO: This method should be implemented if an incorrect extrapolation of the pressure and temperature is encountered in the kernel /** * @copydoc MultiFluidBase::checkTablesParameters( real64 pressure, real64 temperature ) @@ -186,6 +189,18 @@ using CompositionalKValuePhillipsBrine = CompositionalMultiphaseFluid< compositional::KValueFlashModel< 2 >, compositional::PhaseModel< compositional::PhillipsBrineDensity, compositional::PhillipsBrineViscosity >, compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity > >; +using CompositionalThermalTwoPhaseLohrenzBrayClarkViscosity = CompositionalMultiphaseFluid< + compositional::NegativeTwoPhaseFlashModel, + compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity, compositional::CompositionalEnthalpy >, + compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity, compositional::CompositionalEnthalpy > >; +using CompositionalThermalTwoPhasePhillipsBrine = CompositionalMultiphaseFluid< + compositional::NegativeTwoPhaseFlashModel, + compositional::PhaseModel< compositional::PhillipsBrineDensity, compositional::PhillipsBrineViscosity, compositional::CompositionalEnthalpy >, + compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity, compositional::CompositionalEnthalpy > >; +using CompositionalThermalKValuePhillipsBrine = CompositionalMultiphaseFluid< + compositional::KValueFlashModel< 2 >, + compositional::PhaseModel< compositional::PhillipsBrineDensity, compositional::PhillipsBrineViscosity, compositional::CompositionalEnthalpy >, + compositional::PhaseModel< compositional::CompositionalDensity, compositional::LohrenzBrayClarkViscosity, compositional::CompositionalEnthalpy > >; } /* namespace constitutive */ diff --git a/src/coreComponents/constitutive/fluid/multifluid/compositional/CompositionalMultiphaseFluidUpdates.hpp b/src/coreComponents/constitutive/fluid/multifluid/compositional/CompositionalMultiphaseFluidUpdates.hpp index fb15d1a7c1c..510c1210435 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/compositional/CompositionalMultiphaseFluidUpdates.hpp +++ b/src/coreComponents/constitutive/fluid/multifluid/compositional/CompositionalMultiphaseFluidUpdates.hpp @@ -366,11 +366,14 @@ CompositionalMultiphaseFluidUpdates< FLASH, PHASES... >::compute( } // 5. Calculate the internal energy + // The enthalpy is in J/kg if m_useMass and in J/mol otherwise. At this stage phaseDensity + // is still the molar density (it is overwritten by the mass density in step 6), so the + // density consistent with the enthalpy units must be selected explicitly. if constexpr (isThermalType()) { computeInternalEnergy( pressure, phaseFraction, - phaseMassDensity, + m_useMass ? phaseMassDensity : phaseDensity, phaseEnthalpy, phaseInternalEnergy ); } diff --git a/src/coreComponents/constitutive/fluid/multifluid/compositional/models/CompositionalEnthalpy.cpp b/src/coreComponents/constitutive/fluid/multifluid/compositional/models/CompositionalEnthalpy.cpp index 96e51de6ba1..85d481f1be7 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/compositional/models/CompositionalEnthalpy.cpp +++ b/src/coreComponents/constitutive/fluid/multifluid/compositional/models/CompositionalEnthalpy.cpp @@ -30,14 +30,12 @@ namespace compositional CompositionalEnthalpyUpdate::CompositionalEnthalpyUpdate( integer const phaseIndex, integer const vapourIndex, - PhaseType const phaseType, EquationOfStateType const equationOfState, real64 const refTemperature, arrayView2d< real64 const > const & referenceEnthalpy, arrayView2d< real64 const > const & coefficients ): m_phaseIndex( phaseIndex ), m_vapourIndex( vapourIndex ), - m_phaseType( phaseType ), m_equationOfState( equationOfState ), m_refTemperature( refTemperature ), m_referenceEnthalpy( referenceEnthalpy ), @@ -58,10 +56,9 @@ CompositionalEnthalpy::CompositionalEnthalpy( string const & name, integer const numPhases = m_heatCapacityCoefficients->m_phaseTypes.size(); - m_phaseType = m_heatCapacityCoefficients->m_phaseTypes[phaseIndex]; for( integer ip = 0; ip < numPhases; ip++ ) { - if( m_heatCapacityCoefficients->m_phaseTypes[ip] == PhaseType::VAPOUR ) + if( isPhaseType( m_heatCapacityCoefficients->m_phaseTypes[ip], PhaseType::VAPOUR ) ) { m_vapourIndex = ip; } @@ -73,7 +70,6 @@ CompositionalEnthalpy::createKernelWrapper() const { return KernelWrapper( m_phaseIndex, m_vapourIndex, - m_phaseType, m_equationOfState, m_heatCapacityCoefficients->m_referenceTemperature, m_heatCapacityCoefficients->m_referenceEnthalpy.toViewConst(), diff --git a/src/coreComponents/constitutive/fluid/multifluid/compositional/models/CompositionalEnthalpy.hpp b/src/coreComponents/constitutive/fluid/multifluid/compositional/models/CompositionalEnthalpy.hpp index 24ac6fe88c8..36088fc8b40 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/compositional/models/CompositionalEnthalpy.hpp +++ b/src/coreComponents/constitutive/fluid/multifluid/compositional/models/CompositionalEnthalpy.hpp @@ -47,7 +47,6 @@ class CompositionalEnthalpyUpdate final : public FunctionBaseUpdate public: CompositionalEnthalpyUpdate( integer const phaseIndex, integer const vapourIndex, - PhaseType const phaseType, EquationOfStateType const equationOfState, real64 const refTemperature, arrayView2d< real64 const > const & referenceEnthalpy, @@ -102,13 +101,11 @@ class CompositionalEnthalpyUpdate final : public FunctionBaseUpdate real64 const & temperature, arraySlice1d< real64 const, USD1 > const & phaseComposition, real64 & enthalpy, - arraySlice1d< real64, USD2 > const & dEnthalpy, - PhaseType const & phaseType ) const; + arraySlice1d< real64, USD2 > const & dEnthalpy ) const; private: integer const m_phaseIndex{-1}; integer const m_vapourIndex{-1}; - PhaseType const m_phaseType; EquationOfStateType const m_equationOfState; real64 const m_refTemperature; arrayView2d< real64 const > const m_referenceEnthalpy; @@ -145,7 +142,6 @@ class CompositionalEnthalpy : public FunctionBase private: integer m_phaseIndex{-1}; integer m_vapourIndex{-1}; - PhaseType m_phaseType; EquationOfStateType m_equationOfState; arrayView1d< real64 const > m_criticalTemperature; HeatCapacityCoefficients const * const m_heatCapacityCoefficients{nullptr}; @@ -208,8 +204,7 @@ void CompositionalEnthalpyUpdate::compute( temperature, phaseComposition, hEosEnthalpy, - dhEosEnthalpy.toSlice(), - m_phaseType ); + dhEosEnthalpy.toSlice() ); } ); // 3. Combine the 2 values for the actual fluid enthalpy @@ -271,8 +266,7 @@ void CompositionalEnthalpyUpdate::calculateEquationOfStateEnthalpy( real64 const & temperature, arraySlice1d< real64 const, USD1 > const & phaseComposition, real64 & enthalpy, - arraySlice1d< real64, USD2 > const & dEnthalpy, - PhaseType const & phaseType ) const + arraySlice1d< real64, USD2 > const & dEnthalpy ) const { integer sizes[2] = {0, 0}; arraySlice2d< real64 const > derivs( nullptr, sizes, sizes ); @@ -280,8 +274,7 @@ void CompositionalEnthalpyUpdate::calculateEquationOfStateEnthalpy( componentProperties.m_componentBinaryCoeff.toSliceConst(), derivs.toSliceConst() ); - typename EOS::SelectedRoot const root = phaseType == PhaseType::LIQUID ? - EOS::SelectedRoot::MINIMUM : EOS::SelectedRoot::MAXIMUM; + typename EOS::SelectedRoot constexpr root = EOS::SelectedRoot::AUTO; EOS::initialiseStack( numComps, pressure, temperature, componentProperties, stack ); EOS::computeMixtureCoefficients( numComps, phaseComposition, stack ); diff --git a/src/coreComponents/constitutive/fluid/multifluid/compositional/parameters/HeatCapacityCoefficients.cpp b/src/coreComponents/constitutive/fluid/multifluid/compositional/parameters/HeatCapacityCoefficients.cpp index 72ba65a001a..db550a024a1 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/compositional/parameters/HeatCapacityCoefficients.cpp +++ b/src/coreComponents/constitutive/fluid/multifluid/compositional/parameters/HeatCapacityCoefficients.cpp @@ -23,7 +23,7 @@ #include "constitutive/fluid/multifluid/MultiFluidBase.hpp" #include "constitutive/fluid/multifluid/MultiFluidConstants.hpp" #include "dataRepository/InputFlags.hpp" -#include "common/PhysicsConstants.hpp" +#include "common/Units.hpp" namespace geos { @@ -63,11 +63,19 @@ void HeatCapacityCoefficients::registerParametersImpl( MultiFluidBase * fluid ) fluid->registerWrapper( viewKeyStruct::componentHeatCapacityCoefficientsString(), &m_coefficients ). setInputFlag( dataRepository::InputFlags::REQUIRED ). setDescription( "The polynomial coefficients for the specific heat capacity of each component." ); + + // Register extra wrappers to enable auto-cloning + fluid->registerWrapper( "heatCapacityPhaseTypes", &m_phaseTypes ). + setInputFlag( dataRepository::InputFlags::FALSE ). + setSizedFromParent( 0 ). + setRestartFlags( dataRepository::RestartFlags::NO_WRITE ); } void HeatCapacityCoefficients::postInputInitializationImpl( MultiFluidBase const * fluid, ComponentProperties const & componentProperties ) { + GEOS_UNUSED_VAR( componentProperties ); + integer const numPhases = fluid->numFluidPhases(); integer const numComps = fluid->numFluidComponents(); @@ -105,7 +113,7 @@ void HeatCapacityCoefficients::postInputInitializationImpl( MultiFluidBase const for( integer ip = 0; ip < numPhases; ip++ ) { PhaseType const phaseType = getPhaseTypeFromName( phaseNames[ip] ); - m_phaseTypes[ip] = phaseType; + m_phaseTypes[ip] = static_cast< integer >(phaseType); switch( phaseType ) { case PhaseType::LIQUID: @@ -176,25 +184,9 @@ void HeatCapacityCoefficients::postInputInitializationImpl( MultiFluidBase const m_referenceTemperature ), InputError ); - // Determine temperature range - constexpr real64 zeroC = constants::zeroDegreesCelsiusInKelvin; - real64 minTemperature = LvArray::math::min( zeroC, m_referenceTemperature ); - real64 maxTemperature = LvArray::math::max( zeroC, m_referenceTemperature ); - auto const criticalTemperature = componentProperties.getComponentCriticalTemperature(); - for( integer ic = 0; ic < numComps; ++ic ) - { - minTemperature = LvArray::math::min( minTemperature, criticalTemperature[ic] ); - maxTemperature = LvArray::math::max( maxTemperature, criticalTemperature[ic] ); - } - - // Extend interval by 10% in each direction - real64 const dt = LvArray::math::max( maxTemperature - minTemperature, 100.0 ); - minTemperature -= 0.1 * dt; - maxTemperature += 0.1 * dt; - - // Transform to reference temperature space - minTemperature -= m_referenceTemperature; - maxTemperature -= m_referenceTemperature; + // Follow the CO2 brine example and validate in the range [10,200] deg C + real64 const minTemperature = units::convertCToK( 10.0 ) - m_referenceTemperature; + real64 const maxTemperature = units::convertCToK( 200.0 ) - m_referenceTemperature; real64 negT = 0.0; real64 negHT = 0.0; diff --git a/src/coreComponents/constitutive/fluid/multifluid/compositional/parameters/HeatCapacityCoefficients.hpp b/src/coreComponents/constitutive/fluid/multifluid/compositional/parameters/HeatCapacityCoefficients.hpp index 9d4375c700b..42f369c478a 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/compositional/parameters/HeatCapacityCoefficients.hpp +++ b/src/coreComponents/constitutive/fluid/multifluid/compositional/parameters/HeatCapacityCoefficients.hpp @@ -51,7 +51,7 @@ class HeatCapacityCoefficients : public ModelParameters real64 m_referenceTemperature{298.15}; array2d< real64 > m_referenceEnthalpy; array2d< real64 > m_coefficients; - array1d< PhaseType > m_phaseTypes; + array1d< integer > m_phaseTypes; protected: void registerParametersImpl( MultiFluidBase * fluid ) override; diff --git a/src/coreComponents/constitutive/unitTests/CMakeLists.txt b/src/coreComponents/constitutive/unitTests/CMakeLists.txt index 94276e254da..d7a35a0754d 100644 --- a/src/coreComponents/constitutive/unitTests/CMakeLists.txt +++ b/src/coreComponents/constitutive/unitTests/CMakeLists.txt @@ -24,6 +24,7 @@ set( gtest_geosx_tests testModifiedCamClay.cpp testMultiFluidBlackOil.cpp testMultiFluidCO2Brine.cpp + testMultiFluidCompositionalThermalMultiphase.cpp testMultiFluidThreePhaseCompositionalMultiphase.cpp testMultiFluidTwoPhaseCompositionalMultiphase.cpp testMultiFluidDeadOil.cpp diff --git a/src/coreComponents/constitutive/unitTests/testCompositionalEnthalpy.cpp b/src/coreComponents/constitutive/unitTests/testCompositionalEnthalpy.cpp index bf52664e2ac..ffe75dce801 100644 --- a/src/coreComponents/constitutive/unitTests/testCompositionalEnthalpy.cpp +++ b/src/coreComponents/constitutive/unitTests/testCompositionalEnthalpy.cpp @@ -74,8 +74,11 @@ struct FluidData< 4 > } } coefficients->m_phaseTypes.resize( numPhases ); - coefficients->m_phaseTypes[static_cast< integer >(PhaseType::LIQUID)] = PhaseType::LIQUID; - coefficients->m_phaseTypes[static_cast< integer >(PhaseType::VAPOUR)] = PhaseType::VAPOUR; + for( PhaseType const phaseType : {PhaseType::LIQUID, PhaseType::VAPOUR} ) + { + integer const phase = static_cast< integer >(phaseType); + coefficients->m_phaseTypes[phase] = phase; + } } }; @@ -113,8 +116,11 @@ struct FluidData< 6 > } } coefficients->m_phaseTypes.resize( numPhases ); - coefficients->m_phaseTypes[static_cast< integer >(PhaseType::LIQUID)] = PhaseType::LIQUID; - coefficients->m_phaseTypes[static_cast< integer >(PhaseType::VAPOUR)] = PhaseType::VAPOUR; + for( PhaseType const phaseType : {PhaseType::LIQUID, PhaseType::VAPOUR} ) + { + integer const phase = static_cast< integer >(phaseType); + coefficients->m_phaseTypes[phase] = phase; + } } }; diff --git a/src/coreComponents/constitutive/unitTests/testHeatCapacityCoefficients.cpp b/src/coreComponents/constitutive/unitTests/testHeatCapacityCoefficients.cpp index fb04fc5314d..e86f5810e0b 100644 --- a/src/coreComponents/constitutive/unitTests/testHeatCapacityCoefficients.cpp +++ b/src/coreComponents/constitutive/unitTests/testHeatCapacityCoefficients.cpp @@ -196,12 +196,12 @@ TEST_F( HeatCapacityCoefficientsTestFixture, testHeatCapacityCoefficients ) setValues( {1.0, 0.0, 0.0, 0.0, 0.0} ); EXPECT_NO_THROW( m_parameters->postInputInitialization( m_fluid.get(), componentProperties ) ); - setValues( {9.7529018563e+00, -1.7379753331e-01, 1.3828200497e-03, -4.4470070713e-06, 4.3687436501e-09} ); + setValues( {9.7529018563e+00, -1.7379753331e-01, 1.3828200497e-03, -4.4470070713e-06, -2.35e-08} ); EXPECT_THROW( m_parameters->postInputInitialization( m_fluid.get(), componentProperties ), InputError ); setValues( {1.0, 0.0, 0.0, 0.0, 0.0} ); EXPECT_NO_THROW( m_parameters->postInputInitialization( m_fluid.get(), componentProperties ) ); - setValues( {8.440483444351720e-02, -8.156515860472481e-04, 2.921178619177309e-06, -4.649781358999949e-09, 2.775470116151693e-12} ); + setValues( {8.440483444351720e-02, -8.156515860472481e-04, 2.921178619177309e-06, -4.649781358999949e-09, -3.597e-10} ); EXPECT_THROW( m_parameters->postInputInitialization( m_fluid.get(), componentProperties ), InputError ); setValues( {1.0, 0.0, 0.0, 0.0, 0.0} ); EXPECT_NO_THROW( m_parameters->postInputInitialization( m_fluid.get(), componentProperties ) ); diff --git a/src/coreComponents/constitutive/unitTests/testMultiFluidCompositionalThermalMultiphase.cpp b/src/coreComponents/constitutive/unitTests/testMultiFluidCompositionalThermalMultiphase.cpp new file mode 100644 index 00000000000..6e0351ede15 --- /dev/null +++ b/src/coreComponents/constitutive/unitTests/testMultiFluidCompositionalThermalMultiphase.cpp @@ -0,0 +1,582 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 testMultiFluidCompositionalThermalMultiphase.cpp + */ + +#include "FluidModelTest.hpp" +#include "constitutive/fluid/multifluid/compositional/CompositionalMultiphaseFluid.hpp" +#include "constitutive/fluid/multifluid/compositional/parameters/HeatCapacityCoefficients.hpp" +#include "constitutive/fluid/multifluid/compositional/parameters/KValueFlashParameters.hpp" +#include "common/initializeEnvironment.hpp" + +using namespace geos::constitutive; +using namespace geos::constitutive::compositional; + +namespace geos +{ +namespace testing +{ + +enum class FLASH_TYPE : int { NEGATIVE, KVALUE }; +ENUM_STRINGS( FLASH_TYPE, "Negative", "KValue" ); + +enum class VISCOSITY_TYPE : int { LBC, PHILLIPS }; +ENUM_STRINGS( VISCOSITY_TYPE, "LohrenzBrayClark", "Phillips" ); + +template< FLASH_TYPE FLASH, VISCOSITY_TYPE VISCOSITY > +struct FluidModel {}; + +template<> +struct FluidModel< FLASH_TYPE::NEGATIVE, VISCOSITY_TYPE::LBC > +{ + using type = CompositionalThermalTwoPhaseLohrenzBrayClarkViscosity; +}; +template<> +struct FluidModel< FLASH_TYPE::NEGATIVE, VISCOSITY_TYPE::PHILLIPS > +{ + using type = CompositionalThermalTwoPhasePhillipsBrine; +}; +template<> +struct FluidModel< FLASH_TYPE::KVALUE, VISCOSITY_TYPE::PHILLIPS > +{ + using type = CompositionalThermalKValuePhillipsBrine; +}; + +template< typename FluidModel, integer NUM_COMP > +struct FluidData +{}; + +template< typename TEST_TYPE > +class MultiThermalFluidCompositionalMultiphaseTestFixture : public FluidModelTest< + typename FluidModel< std::tuple_element_t< 1, TEST_TYPE >::value, + std::tuple_element_t< 2, TEST_TYPE >::value >::type, + std::tuple_element_t< 3, TEST_TYPE >::value > +{ +public: + static constexpr EquationOfStateType EquationOfState = std::tuple_element_t< 0, TEST_TYPE >::value; + static constexpr FLASH_TYPE FlashModel = std::tuple_element_t< 1, TEST_TYPE >::value; + static constexpr VISCOSITY_TYPE ViscosityModel = std::tuple_element_t< 2, TEST_TYPE >::value; + using CompositionalMultiphaseFluid = typename FluidModel< FlashModel, ViscosityModel >::type; + using Base = FluidModelTest< CompositionalMultiphaseFluid, std::tuple_element_t< 3, TEST_TYPE >::value >; + static constexpr real64 relTol = 1.0e-4; + static constexpr real64 absTol = 1.0e-3; + + static constexpr integer pointCount = 10; + + static constexpr real64 minPressure = 0.99; + static constexpr real64 maxPressure = 600.1; + + static constexpr real64 minTemperature = 9.9; + static constexpr real64 maxTemperature = 151.0; + +public: + MultiThermalFluidCompositionalMultiphaseTestFixture() + { + Base::createFluid( getFluidName(), [&]( CompositionalMultiphaseFluid & fluid ){ + fillPhysicalProperties( fluid ); + if constexpr (ViscosityModel == VISCOSITY_TYPE::PHILLIPS) + { + generateInterpolationPoints( fluid ); + } + if constexpr (FlashModel == FLASH_TYPE::KVALUE) + { + generateTables( fluid ); + } + } ); + } + + ~MultiThermalFluidCompositionalMultiphaseTestFixture() override + { + for( string const & fileName : m_fileNames ) + { + removeFile( fileName ); + } + } + + void testNumericalDerivatives( const bool useMass ) + { + CompositionalMultiphaseFluid * fluid = this->getFluid( this->getFluidName() ); + + fluid->setMassFlag( useMass ); + + array2d< real64 > samples; + FluidData< CompositionalMultiphaseFluid, Base::numComp >::getSamples( samples ); + integer const sampleCount = samples.size( 0 ); + Feed< Base::numComp > sample; + + real64 constexpr eps = 1.0e-6; + + constexpr real64 pressures[] = { 10.0e5, 50.0e5, 100.0e5, 600.0e5 }; + constexpr real64 temperatures[] = { 15.5, 24.0, 40.0, 80.0 }; + + for( integer sampleIndex = 0; sampleIndex < sampleCount; ++sampleIndex ) + { + for( integer ic = 0; ic < Base::numComp; ++ic ) + { + sample[ic] = samples( sampleIndex, ic ); + } + for( real64 const pressure : pressures ) + { + for( real64 const temperature : temperatures ) + { + typename Base::TestPoint const data ( pressure, units::convertCToK( temperature ), sample ); + Base::testNumericalDerivatives( fluid, data, eps, relTol, absTol ); + } + } + } + } + + // Checks that u = h - p/rho for every existing phase. After the update, the phase density is + // in units consistent with the phase enthalpy (mass or molar), whatever the formulation. + void testInternalEnergyConsistency( const bool useMass ) + { + CompositionalMultiphaseFluid * fluid = this->getFluid( this->getFluidName() ); + + fluid->setMassFlag( useMass ); + + this->m_parent.resize( 1 ); + fluid->allocateConstitutiveData( this->m_parent, 1 ); + + array2d< real64 > samples; + FluidData< CompositionalMultiphaseFluid, Base::numComp >::getSamples( samples ); + integer const sampleCount = samples.size( 0 ); + + array2d< real64, compflow::LAYOUT_COMP > compositionArray( 1, Base::numComp ); + + constexpr real64 pressures[] = { 10.0e5, 50.0e5, 100.0e5, 600.0e5 }; + constexpr real64 temperatures[] = { 15.5, 24.0, 40.0, 80.0 }; + + string_array const & phaseNames = fluid->phaseNames(); + auto const fluidWrapper = fluid->createKernelWrapper(); + + for( integer sampleIndex = 0; sampleIndex < sampleCount; ++sampleIndex ) + { + for( integer ic = 0; ic < Base::numComp; ++ic ) + { + compositionArray( 0, ic ) = samples( sampleIndex, ic ); + } + auto const compositionView = compositionArray.toViewConst(); + + for( real64 const pressure : pressures ) + { + for( real64 const temperatureInC : temperatures ) + { + real64 const temperature = units::convertCToK( temperatureInC ); + + forAll< serialPolicy >( 1, [fluidWrapper, pressure, temperature, compositionView] ( localIndex const k ) + { + fluidWrapper.update( k, 0, pressure, temperature, compositionView[k] ); + } ); + + for( integer ip = 0; ip < Base::numPhase; ++ip ) + { + if( fluid->phaseFraction()( 0, 0, ip ) <= 0.0 ) + { + continue; + } + real64 const phaseEnthalpy = fluid->phaseEnthalpy()( 0, 0, ip ); + real64 const phaseDensity = fluid->phaseDensity()( 0, 0, ip ); + real64 const phaseInternalEnergy = fluid->phaseInternalEnergy()( 0, 0, ip ); + real64 const expectedInternalEnergy = phaseEnthalpy - pressure / phaseDensity; + checkRelativeError( phaseInternalEnergy, expectedInternalEnergy, relTol, absTol, + GEOS_FMT( "Phase internal energy ({}) at P = {}, T = {}, z = {}", + phaseNames[ip], pressure, temperature, + Base::toString( compositionView[0].toSliceConst() ) ) ); + } + } + } + } + } + + static string getFluidName(); + +private: + static void fillPhysicalProperties( CompositionalMultiphaseFluid & fluid ); + void generateInterpolationPoints( CompositionalMultiphaseFluid & fluid ); + void generateTables( CompositionalMultiphaseFluid & fluid ); + + static void writeToFile( string const & fileName, string const & content ); + + static void removeFile( string const & fileName ); +private: + string_array m_fileNames; +}; + +template< typename TEST_TYPE > +string MultiThermalFluidCompositionalMultiphaseTestFixture< TEST_TYPE >::getFluidName() +{ + return GEOS_FMT( "fluid_{}_{}_{}_{}", + Base::numComp, + EnumStrings< EquationOfStateType >::toString( EquationOfState ), + EnumStrings< FLASH_TYPE >::toString( FlashModel ), + EnumStrings< VISCOSITY_TYPE >::toString( ViscosityModel ) ); +} + +template< integer NUM_COMP > +static void fillBinaryCoeffs( array2d< real64 > & binaryCoeff, std::array< real64 const, NUM_COMP *(NUM_COMP-1)/2 > const data ) +{ + auto bic = data.begin(); + binaryCoeff.resize( NUM_COMP, NUM_COMP ); + for( integer i = 0; i < NUM_COMP; ++i ) + { + binaryCoeff( i, i ) = 0.0; + for( integer j = i+1; j < NUM_COMP; ++j ) + { + binaryCoeff( i, j ) = *bic++; + binaryCoeff( j, i ) = binaryCoeff( i, j ); + } + } +} + +template< integer NUM_COMP > +static void populateArray( arraySlice1d< real64 > array, std::array< real64 const, NUM_COMP > const data ) +{ + for( integer i = 0; i < NUM_COMP; ++i ) + { + array[i] = data[i]; + } +} + +template< typename FluidModel > +struct FluidData< FluidModel, 4 > +{ + static constexpr integer numComps = 4; + static constexpr integer numCoeffs = 5; + + static void fillProperties( dataRepository::Group & fluid ) + { + using Keys = typename FluidModel::viewKeyStruct; + + string_array & componentNames = fluid.getReference< string_array >( Keys::componentNamesString() ); + componentNames = {"N2", "C5", "C10", "H2O"}; + + array1d< real64 > & molarWeight = fluid.getReference< array1d< real64 > >( Keys::componentMolarWeightString() ); + TestFluid< numComps >::createArray( molarWeight, Feed< numComps >{28e-3, 134e-3, 142e-3, 18e-3} ); + + array1d< real64 > & criticalPressure = fluid.getReference< array1d< real64 > >( Keys::componentCriticalPressureString() ); + TestFluid< numComps >::createArray( criticalPressure, Feed< numComps >{34e5, 33.68e5, 21.03e5, 220.5e5} ); + array1d< real64 > & criticalTemperature = fluid.getReference< array1d< real64 > >( Keys::componentCriticalTemperatureString() ); + TestFluid< numComps >::createArray( criticalTemperature, Feed< numComps >{126.2, 469.7, 617.7, 647.0} ); + array1d< real64 > & acentricFactor = fluid.getReference< array1d< real64 > >( Keys::componentAcentricFactorString() ); + TestFluid< numComps >::createArray( acentricFactor, Feed< numComps >{0.04, 0.2510, 0.4884, 0.344} ); + array2d< real64 > & binaryCoeff = fluid.getReference< array2d< real64 > >( Keys::componentBinaryCoeffString() ); + fillBinaryCoeffs< numComps >( binaryCoeff, {0.0, 0.1, 0.0, 0.0, 0.0, 0.0} ); + + array2d< real64 > & referenceEnthalpy = fluid.getReference< array2d< real64 > >( HeatCapacityCoefficients::viewKeyStruct::referenceEnthalpyString() ); + referenceEnthalpy.resize( 2, numComps ); + referenceEnthalpy.zero(); + + array2d< real64 > & coefficients = fluid.getReference< array2d< real64 > >( HeatCapacityCoefficients::viewKeyStruct::componentHeatCapacityCoefficientsString() ); + coefficients.resize( numComps, numCoeffs ); + std::array< real64, numComps * numCoeffs > coefficientsData { + 2.89614125e+01, 2.03166820e-03, 2.47038198e-05, -9.54674193e-08, 1.15955044e-10, /* N2 */ + 1.21004526e+02, 2.47982404e-01, 7.27115529e-04, -2.64999262e-06, 2.52476069e-09, /* C5 */ + 3.11952000e+02, 5.26890000e-01, 8.16470000e-04, -4.03480000e-06, 9.60000000e-09, /* C10 */ + 2.69032982e+01, 5.51993739e-02, 6.31405927e-04, -3.66166336e-06, 4.56291585e-09, /* H2O */ + }; + for( int ic = 0; ic < numComps; ++ic ) + { + for( int j = 0; j < numCoeffs; ++j ) + { + coefficients( ic, j ) = coefficientsData[ic*numCoeffs + j]; + } + } + } + + static void getSamples( array2d< real64 > & samples ) + { + samples.resize( 3, numComps ); + populateArray< numComps >( samples[0], {0.099, 0.300, 0.600, 0.001} ); + populateArray< numComps >( samples[1], {0.100, 0.250, 0.350, 0.300} ); + populateArray< numComps >( samples[2], {0.050, 0.050, 0.100, 0.800} ); + } +}; + +template< typename FluidModel > +struct FluidData< FluidModel, 5 > +{ + static constexpr integer numComps = 5; + static constexpr integer numCoeffs = 5; + + static void fillProperties( dataRepository::Group & fluid ) + { + using Keys = typename FluidModel::viewKeyStruct; + + string_array & componentNames = fluid.getReference< string_array >( Keys::componentNamesString() ); + componentNames = {"CO2", "N2", "C1", "H2O", "C4"}; + + array1d< real64 > & molarWeight = fluid.getReference< array1d< real64 > >( Keys::componentMolarWeightString() ); + TestFluid< 5 >::createArray( molarWeight, Feed< 5 >{44.0098e-3, 28.0135e-3, 16.0428e-3, 18e-3, 82.4191e-3} ); + + array1d< real64 > & criticalPressure = fluid.getReference< array1d< real64 > >( Keys::componentCriticalPressureString() ); + TestFluid< 5 >::createArray( criticalPressure, Feed< 5 >{73.77300e5, 33.95800e5, 45.99200e5, 48.71800e5, 33.20710e5} ); + array1d< real64 > & criticalTemperature = fluid.getReference< array1d< real64 > >( Keys::componentCriticalTemperatureString() ); + TestFluid< 5 >::createArray( criticalTemperature, Feed< 5 >{304.1280, 126.1920, 190.5640, 305.3300, 504.2160} ); + array1d< real64 > & acentricFactor = fluid.getReference< array1d< real64 > >( Keys::componentAcentricFactorString() ); + TestFluid< 5 >::createArray( acentricFactor, Feed< 5 >{0.223000, 0.037200, 0.010400, 0.099100, 0.250274} ); + array1d< real64 > & volumeShift = fluid.getReference< array1d< real64 > >( Keys::componentVolumeShiftString() ); + TestFluid< 5 >::createArray( volumeShift, Feed< 5 >{1.845465e-01, -1.283880e-01, 9.225800e-02, 6.458060e-02, 0.000000e+00} ); + array2d< real64 > & binaryCoeff = fluid.getReference< array2d< real64 > >( Keys::componentBinaryCoeffString() ); + fillBinaryCoeffs< 5 >( binaryCoeff, {0.0, 0.1, 0.03, 0.139, 0.032, 0.0, 0.12, 0.03, 0.0, 0.0} ); + + array2d< real64 > & referenceEnthalpy = fluid.getReference< array2d< real64 > >( HeatCapacityCoefficients::viewKeyStruct::referenceEnthalpyString() ); + referenceEnthalpy.resize( 2, numComps ); + referenceEnthalpy.zero(); + + array2d< real64 > & coefficients = fluid.getReference< array2d< real64 > >( HeatCapacityCoefficients::viewKeyStruct::componentHeatCapacityCoefficientsString() ); + coefficients.resize( numComps, numCoeffs ); + std::array< real64, numComps * numCoeffs > coefficientsData { + 3.40933377e+01, 7.52742132e-02, 2.00271751e-04, -1.61900936e-06, 2.12199565e-09, /* CO2 */ + 2.89614125e+01, 2.03166820e-03, 2.47038198e-05, -9.54674193e-08, 1.15955044e-10, /* N2 */ + 3.54137115e+01, 4.34392612e-02, 1.42852598e-04, -4.70393972e-07, 4.80614694e-10, /* CH4 */ + 2.69032982e+01, 5.51993739e-02, 6.31405927e-04, -3.66166336e-06, 4.56291585e-09, /* H2O */ + 5.28630051e+01, 1.16281994e-01, 1.40990713e-04, -6.38108433e-07, 7.14166667e-10, /* C4H10 */ + }; + for( int ic = 0; ic < numComps; ++ic ) + { + for( int j = 0; j < numCoeffs; ++j ) + { + coefficients( ic, j ) = coefficientsData[ic*numCoeffs + j]; + } + } + } + + static void getSamples( array2d< real64 > & samples ) + { + samples.resize( 1, 5 ); + populateArray< 5 >( samples[0], {0.050, 0.150, 0.550, 0.150, 0.100} ); + } +}; + +template< typename TEST_TYPE > +void MultiThermalFluidCompositionalMultiphaseTestFixture< TEST_TYPE >::fillPhysicalProperties( CompositionalMultiphaseFluid & fluid ) +{ + string_array & phaseNames = fluid.template getReference< string_array >( MultiFluidBase::viewKeyStruct::phaseNamesString() ); + phaseNames = {"oil", "gas"}; + + string const eosName = EnumStrings< EquationOfStateType >::toString( EquationOfState ); + string_array & equationOfState = fluid.template getReference< string_array >( EquationOfState::viewKeyStruct::equationsOfStateString() ); + equationOfState = {eosName, eosName}; + if constexpr (EquationOfState == EquationOfStateType::SoreideWhitson) + { + phaseNames[0] = "water"; + equationOfState[1] = EnumStrings< EquationOfStateType >::toString( EquationOfStateType::PengRobinson ); + } + + dataRepository::Group & group = fluid; + FluidData< CompositionalMultiphaseFluid, Base::numComp >::fillProperties( group ); +} + +template< typename TEST_TYPE > +void MultiThermalFluidCompositionalMultiphaseTestFixture< TEST_TYPE >::generateInterpolationPoints( CompositionalMultiphaseFluid & fluid ) +{ + integer constexpr N = pointCount; + + array1d< real64 > & pressureCoordinates = fluid.template getReference< array1d< real64 > >( PressureTemperatureCoordinates::viewKeyStruct::pressureCoordinatesString() ); + array1d< real64 > & temperatureCoordinates = fluid.template getReference< array1d< real64 > >( PressureTemperatureCoordinates::viewKeyStruct::temperatureCoordinatesString() ); + pressureCoordinates.resize( N+1 ); + temperatureCoordinates.resize( N+1 ); + + real64 const dp = pow( maxPressure/minPressure, 1.0/N ); + real64 constexpr BAR_2_PA = 1.0e5; + pressureCoordinates[0] = BAR_2_PA * minPressure; + for( integer i = 1; i <= N; ++i ) + { + pressureCoordinates[i] = pressureCoordinates[i-1] * dp; + } + + real64 const dT = (maxTemperature-minTemperature)/N; + for( integer i = 0; i <= N; ++i ) + { + temperatureCoordinates[i] = minTemperature + i*dT + 273.15; + } +} + +template< typename TEST_TYPE > +void MultiThermalFluidCompositionalMultiphaseTestFixture< TEST_TYPE >::generateTables( CompositionalMultiphaseFluid & fluid ) +{ + string_array & kValueTables = fluid.template getReference< string_array >( KValueFlashParameters< 2 >::viewKeyStruct::kValueTablesString() ); + + string_array const & compNames = fluid.componentNames(); + string const fluidName = fluid.getName(); + + // Crookston correlations pressure (bar), temperature (K) + mapBase< string, std::array< real64 const, 5 >, std::false_type > crookstonCoefficients = { + {"H2O", {3.0620e+00, 8.9414e+02, 1.1912e-02, 5.3659e+02, 1.1951e+02}}, + {"CO2", {-1.8141e+00, 6.2655e+02, 6.7489e-03, 5.0732e+00, 2.5249e+02}}, + {"N2", {-3.5742e-01, 4.8660e+02, 5.7887e-03, 9.1910e+01, 1.8027e+02}}, + {"C1", {3.6577e+00, 7.1889e+02, 1.4154e-02, 5.7323e+02, 1.1088e+02}}, + {"C2", {8.2845e+00, 1.0144e+03, 5.2968e-02, 9.5800e+02, 9.3548e+01}}, + {"C4", {2.5068e+01, 2.2625e+03, 5.3888e-01, 1.7981e+03, 7.1948e+01}}, + {"C5", {3.6243e+01, 4.1985e+03, 1.9055e+00, 2.3575e+03, 5.8025e+01}}, + {"C10", {-2.3346e-01, 7.9356e-01, 6.8406e-03, 7.1217e+02, 1.3827e+02}} + }; + std::ostringstream content; + + // Create arrays of pressure and temperature + integer constexpr N = pointCount; + array1d< real64 > temperature( N+1 ); + array1d< real64 > pressure( N+1 ); + + real64 const dp = pow( maxPressure/minPressure, 1.0/N ); + content.str( "" ); + real64 constexpr BAR_2_PA = 1.0e5; + pressure[0] = minPressure; + content << BAR_2_PA * pressure[0] << "\n"; + for( integer i = 1; i <= N; ++i ) + { + pressure[i] = pressure[i-1] * dp; + content << BAR_2_PA * pressure[i] << "\n"; + } + string const pressureFileName = GEOS_FMT( "{}_PRESSURE.txt", fluidName ); + writeToFile( pressureFileName, content.str()); + m_fileNames.emplace_back( pressureFileName ); + + real64 const dT = (maxTemperature-minTemperature)/N; + content.str( "" ); + for( integer i = 0; i <= N; ++i ) + { + temperature[i] = minTemperature + i*dT + 273.15; + content << temperature[i] << "\n"; + } + string const temperatureFileName = GEOS_FMT( "{}_TEMP.txt", fluidName ); + writeToFile( temperatureFileName, content.str()); + m_fileNames.emplace_back( temperatureFileName ); + + FunctionManager & functionManager = FunctionManager::getInstance(); + + kValueTables.resize( Base::numComp ); + for( integer ic = 0; ic < Base::numComp; ++ic ) + { + auto const it = crookstonCoefficients.find( compNames[ic] ); + ASSERT_FALSE( it == crookstonCoefficients.end() ); + auto const [A, B, C, D, E] = std::apply( []( auto &&... args ) { return std::make_tuple( args ... ); }, it->second ); + string const tableName = GEOS_FMT( "{}_KVALUE_{}", fluidName, ic+1 ); + content.str( "" ); + for( integer pi = 0; pi <= N; ++pi ) + { + for( integer ti = 0; ti <= N; ++ti ) + { + real64 const kValue = (A + B/pressure[pi] + C*pressure[pi])*LvArray::math::exp( -D/(temperature[ti]-E)); + content << LvArray::math::max( kValue, 1.0e-8 ) << "\n"; + } + } + string const kValueFileName = GEOS_FMT( "{}.txt", tableName ); + writeToFile( kValueFileName, content.str()); + m_fileNames.emplace_back( kValueFileName ); + + TableFunction * tableFunction = dynamicCast< TableFunction * >( functionManager.createChild( "TableFunction", tableName ) ); + + string_array & inputVarNames = tableFunction->getWrapper< string_array >( dataRepository::keys::inputVarNames ).reference(); + inputVarNames.emplace_back( "pressure" ); + inputVarNames.emplace_back( "temperature" ); + + path_array & coordinateFiles = tableFunction->getWrapper< path_array >( TableFunction::viewKeyStruct::coordinateFilesString()) + .reference(); + coordinateFiles.emplace_back( Path() ); + coordinateFiles.emplace_back( Path() ); + dynamicCast< string & >( coordinateFiles[0] ) = pressureFileName; + dynamicCast< string & >( coordinateFiles[1] ) = temperatureFileName; + Path & voxelFile = tableFunction->getWrapper< Path >( TableFunction::viewKeyStruct::voxelFileString()) + .reference(); + dynamicCast< string & >( voxelFile ) = kValueFileName; + tableFunction->setInterpolationMethod( TableFunction::InterpolationType::Linear ); + tableFunction->initializeFunction(); + kValueTables[ic] = tableName; + } +} + +template< typename TEST_TYPE > +void MultiThermalFluidCompositionalMultiphaseTestFixture< TEST_TYPE >::writeToFile( string const & fileName, string const & content ) +{ + std::ofstream os( fileName ); + ASSERT_TRUE( os.is_open() ); + os << content; + os.close(); +} + +template< typename TEST_TYPE > +void MultiThermalFluidCompositionalMultiphaseTestFixture< TEST_TYPE >::removeFile( string const & fileName ) +{ + int const ret = std::remove( fileName.c_str() ); + ASSERT_TRUE( ret == 0 ); +} + +template< EquationOfStateType EOS, FLASH_TYPE FLASH, VISCOSITY_TYPE VISCOSITY, integer NUM_COMP > +using TestType = std::tuple< + std::integral_constant< EquationOfStateType, EOS >, + std::integral_constant< FLASH_TYPE, FLASH >, + std::integral_constant< VISCOSITY_TYPE, VISCOSITY >, + std::integral_constant< integer, NUM_COMP > + >; + +using TestTypes = ::testing::Types< + TestType< EquationOfStateType::PengRobinson, FLASH_TYPE::NEGATIVE, VISCOSITY_TYPE::LBC, 4 >, + TestType< EquationOfStateType::SoreideWhitson, FLASH_TYPE::NEGATIVE, VISCOSITY_TYPE::PHILLIPS, 4 >, + TestType< EquationOfStateType::PengRobinson, FLASH_TYPE::KVALUE, VISCOSITY_TYPE::PHILLIPS, 4 >, + TestType< EquationOfStateType::SoaveRedlichKwong, FLASH_TYPE::NEGATIVE, VISCOSITY_TYPE::LBC, 5 >, + TestType< EquationOfStateType::PengRobinson, FLASH_TYPE::NEGATIVE, VISCOSITY_TYPE::PHILLIPS, 5 >, + TestType< EquationOfStateType::PengRobinson, FLASH_TYPE::KVALUE, VISCOSITY_TYPE::PHILLIPS, 5 > + >; + +class NameGenerator +{ +public: + template< typename T > + static std::string GetName( int index ) + { + return GEOS_FMT( "CompositionalFluid_{}_{}", MultiThermalFluidCompositionalMultiphaseTestFixture< T >::getFluidName(), index ); + } +}; + +TYPED_TEST_SUITE( MultiThermalFluidCompositionalMultiphaseTestFixture, TestTypes, NameGenerator ); + +TYPED_TEST( MultiThermalFluidCompositionalMultiphaseTestFixture, numericalDerivativesMolar ) +{ + this->testNumericalDerivatives( false ); +} + +TYPED_TEST( MultiThermalFluidCompositionalMultiphaseTestFixture, numericalDerivativesMass ) +{ + this->testNumericalDerivatives( true ); +} + +TYPED_TEST( MultiThermalFluidCompositionalMultiphaseTestFixture, internalEnergyConsistencyMolar ) +{ + this->testInternalEnergyConsistency( false ); +} + +TYPED_TEST( MultiThermalFluidCompositionalMultiphaseTestFixture, internalEnergyConsistencyMass ) +{ + this->testInternalEnergyConsistency( true ); +} + +} // testing +} // geos + +int main( int argc, char * * argv ) +{ + ::testing::InitGoogleTest( &argc, argv ); + + geos::setupEnvironment( argc, argv ); + + int const result = RUN_ALL_TESTS(); + + geos::cleanupEnvironment(); + + return result; +} diff --git a/src/coreComponents/constitutiveDrivers/kernelSpecs.json b/src/coreComponents/constitutiveDrivers/kernelSpecs.json index 78ad42ef6ef..ccb975cb86e 100644 --- a/src/coreComponents/constitutiveDrivers/kernelSpecs.json +++ b/src/coreComponents/constitutiveDrivers/kernelSpecs.json @@ -38,6 +38,9 @@ "CO2BrinePhillipsThermalFluid", "CompositionalKValuePhillipsBrine", "CompositionalKValueLohrenzBrayClarkViscosity", + "CompositionalThermalTwoPhaseLohrenzBrayClarkViscosity", + "CompositionalThermalTwoPhasePhillipsBrine", + "CompositionalThermalKValuePhillipsBrine", "CompositionalThreePhaseLohrenzBrayClarkViscosity", "CompositionalTwoPhaseConstantViscosity", "CompositionalTwoPhaseLohrenzBrayClarkViscosity", diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernelSpecs.json b/src/coreComponents/physicsSolvers/fluidFlow/kernelSpecs.json index b3d00eda090..51e8043bce7 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernelSpecs.json +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernelSpecs.json @@ -89,6 +89,9 @@ "constitutive::CO2BrineEzrokhiThermalFluid", "constitutive::CompositionalKValuePhillipsBrine", "constitutive::CompositionalKValueLohrenzBrayClarkViscosity", + "constitutive::CompositionalThermalTwoPhaseLohrenzBrayClarkViscosity", + "constitutive::CompositionalThermalTwoPhasePhillipsBrine", + "constitutive::CompositionalThermalKValuePhillipsBrine", "constitutive::CompositionalThreePhaseLohrenzBrayClarkViscosity", "constitutive::CompositionalTwoPhaseConstantViscosity", "constitutive::CompositionalTwoPhaseLohrenzBrayClarkViscosity", @@ -125,7 +128,10 @@ "NCOMP": [ 2 ], "FLUID_TYPE": [ "constitutive::CO2BrinePhillipsThermalFluid", - "constitutive::CO2BrineEzrokhiThermalFluid" + "constitutive::CO2BrineEzrokhiThermalFluid", + "constitutive::CompositionalThermalTwoPhaseLohrenzBrayClarkViscosity", + "constitutive::CompositionalThermalTwoPhasePhillipsBrine", + "constitutive::CompositionalThermalKValuePhillipsBrine" ] }, "NCOMP_BOUNDS": { @@ -202,6 +208,9 @@ "constitutive::CO2BrineEzrokhiThermalFluid", "constitutive::CompositionalKValuePhillipsBrine", "constitutive::CompositionalKValueLohrenzBrayClarkViscosity", + "constitutive::CompositionalThermalTwoPhaseLohrenzBrayClarkViscosity", + "constitutive::CompositionalThermalTwoPhasePhillipsBrine", + "constitutive::CompositionalThermalKValuePhillipsBrine", "constitutive::CompositionalThreePhaseLohrenzBrayClarkViscosity", "constitutive::CompositionalTwoPhaseConstantViscosity", "constitutive::CompositionalTwoPhaseLohrenzBrayClarkViscosity", @@ -228,6 +237,9 @@ "constitutive::CO2BrineEzrokhiThermalFluid", "constitutive::CompositionalKValuePhillipsBrine", "constitutive::CompositionalKValueLohrenzBrayClarkViscosity", + "constitutive::CompositionalThermalTwoPhaseLohrenzBrayClarkViscosity", + "constitutive::CompositionalThermalTwoPhasePhillipsBrine", + "constitutive::CompositionalThermalKValuePhillipsBrine", "constitutive::CompositionalThreePhaseLohrenzBrayClarkViscosity", "constitutive::CompositionalTwoPhaseConstantViscosity", "constitutive::CompositionalTwoPhaseLohrenzBrayClarkViscosity", diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd index d1dd10ec6b0..b4eba69e6f2 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -739,6 +739,18 @@ + + + + + + + + + + + + @@ -7116,6 +7128,9 @@ Information output from lower logLevels is added with the desired log level + + + @@ -7826,6 +7841,168 @@ if the table is requested to be output in the log, and it is too large, a CSV fi + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/coreComponents/schema/schema.xsd.other b/src/coreComponents/schema/schema.xsd.other index 841191862d7..45d143d21b6 100644 --- a/src/coreComponents/schema/schema.xsd.other +++ b/src/coreComponents/schema/schema.xsd.other @@ -1663,6 +1663,9 @@ A field can represent a physical variable. (pressure, temperature, global compos + + + @@ -2428,6 +2431,168 @@ A field can represent a physical variable. (pressure, temperature, global compos + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/docs/sphinx/advancedExamples/validationStudies/carbonStorage/Index.rst b/src/docs/sphinx/advancedExamples/validationStudies/carbonStorage/Index.rst index 4005b439479..56a4f3dd76c 100644 --- a/src/docs/sphinx/advancedExamples/validationStudies/carbonStorage/Index.rst +++ b/src/docs/sphinx/advancedExamples/validationStudies/carbonStorage/Index.rst @@ -15,5 +15,7 @@ Carbon Storage isothermalHystInjection/Example + jouleThomsonCooling/Example + spe11b/Example diff --git a/src/docs/sphinx/advancedExamples/validationStudies/carbonStorage/jouleThomsonCooling/Example.rst b/src/docs/sphinx/advancedExamples/validationStudies/carbonStorage/jouleThomsonCooling/Example.rst new file mode 100644 index 00000000000..4060fd6de8c --- /dev/null +++ b/src/docs/sphinx/advancedExamples/validationStudies/carbonStorage/jouleThomsonCooling/Example.rst @@ -0,0 +1,434 @@ +.. |oslash| unicode:: U+00F8 + +.. |degree| unicode:: U+00B0 + +.. _jouleThomsonCoolingExample: + +########################################################################## +Cooling during CO2 Injection into a Gas Reservoir +########################################################################## + + +**Context** + +When a real gas expands from a high to a low pressure at constant enthalpy, its temperature changes. This is the Joule-Thomson (JT) effect. For CO2 at reservoir conditions the gas cools when it expands. During carbon sequestration, CO2 is injected into a depleted natural gas reservoir and the pressure drops steeply between the wellbore and the reservoir, so JT cooling may occur close to the injection well. If it were large, it could lead to CO2 or CH4 hydrate formation, to freezing of the residual water (with a loss of injectivity), or to thermal stresses able to fracture the formation. + +`(Oldenburg, 2007) `__ used the TOUGH2/EOS7C simulator to quantify the magnitude of this cooling. The model was first verified against the classical JT expansion experiment, and was then applied to an idealized one-dimensional radial gas reservoir with constant-rate CO2 injection. The study showed that the cooling grows when the permeability decreases (a larger pressure drop is required to inject the same mass of CO2), and that it is partly compensated by the heat released by the rock matrix and the residual water. + +This example reproduces the lowest-permeability case of the paper, :math:`k = 5\times 10^{-15}` m\ :sup:`2` (5 mD), whose results are shown in the paper. It is the case with the strongest cooling among the constant-rate injections with an injection temperature equal to the reservoir temperature, and therefore the most demanding one for the thermal formulation. It can serve as a validation of the thermal compositional multiphase flow solver of GEOS, in particular of the coupling between + +- a real-gas equation of state, which provides the enthalpy of the CO2-CH4-H2O mixture and thus the JT coefficient, +- the energy balance, including the advection of enthalpy, the heat stored in the rock and in the immobile water, and + conduction, +- thermal wells with prescribed mass and volumetric rates. + +------------------------------------------------------------------------ +Brief case description +------------------------------------------------------------------------ + +**Reservoir model of the paper** + +The paper considers an idealized depleted gas reservoir, with the geometry and the properties given in the paper. It is represented by a one-dimensional radially symmetric system, 50 m thick and 1130 m in radius, as illustrated below. + +.. code-block:: console + + injection well (r = 0) ---> CO2-rich gas | CH4-rich gas ---> production (r ~ 1100 m) + +The reservoir is initially filled with methane at 50 bar and 45 \ |degree|\ C. The water is immobile (the water saturation is equal to the residual saturation :math:`S_{lr} = 0.2`), and the only mobile phase is the gas. Pure CO2 is injected at a constant rate of 3 kg/s and at the reservoir temperature. Methane is produced at the outer boundary, at 0.56 kg/s. The system is closed: no fluid or heat crosses the outer boundary. Because the injection temperature is equal to the initial reservoir temperature, every temperature change in the reservoir is caused by the JT effect (attenuated by the heat provided by the rock and the water), which makes this setup well suited to isolate and validate this effect. + +The parameters of the problem are summarized in the following table. + +.. list-table:: Parameters of the problem + :widths: 30 70 + :header-rows: 1 + :class: align-top + + * - Parameter + - Value + * - Geometry + - 1D radial, :math:`R = 1130` m, thickness 50 m + + GEOS is modeled as a 1\ |degree|\ sector, :math:`r \in [0.1, 1130]` m. + * - Permeability + - :math:`5\times 10^{-15}` m\ :sup:`2` (5 mD) + * - Porosity + - 0.30 + * - Rock volumetric heat capacity + - :math:`2.6\times 10^{6}` J/m\ :sup:`3`/K + * - Thermal conductivity + - 2.51 W/m/K + * - Residual saturations + - :math:`S_{lr} = 0.2` (water), :math:`S_{gr} = 0.01` (gas) + * - Relative permeability + - Corey curves + * - Initial pressure + - 50 bar (:math:`5\times 10^{6}` Pa) + * - Initial temperature + - 45 \ |degree|\ C (318.15 K) + * - Initial fluids + - CH4 gas, water at :math:`S_l = S_{lr} = 0.2` + + GEOS initializes this via mass fractions: :math:`z_{CO_2} = 0`, :math:`z_{CH_4} = 0.115401`, :math:`z_{H_2O} = 0.884599`. + * - CO2 injection + - 3 kg/s of CO2 at 45 \ |degree|\ C + * - CH4 production + - 0.56 kg/s at the outer boundary (1100 m) + * - Simulated time + - 15 years + * - Equation of state + - Peng-Robinson for the gas phase + + The paper uses EOS7C/GasEOS, while GEOS explicitly uses S\ |oslash|\ reide-Whitson for the aqueous phase. + +------------------------------------------------------------------------ +Input base +------------------------------------------------------------------------ + +This benchmark test is based on the XML file located below: + +.. code-block:: console + + inputFiles/compositionalMultiphaseFlow/benchmarks/jouleThomsonCooling/jouleThomsonCooling.xml + +The mesh is a thin wedge of cylinder (``InternalWellbore``), 1\ |degree|\ wide (``theta`` from -0.5\ |degree|\ to 0.5\ |degree|\ ). The 109 radial cells are defined explicitly through the ``radius`` and ``nr`` arrays: the cell size grows geometrically with the distance to the injection well (the ratio between successive radii decreases from 1.2 close to the well to about 1.08 at the outer boundary), which refines the mesh where the pressure and temperature gradients are the strongest. The 360\ |degree|\ of the paper are accounted for by dividing all the well rates by 360. + +.. literalinclude:: ../../../../../../../inputFiles/compositionalMultiphaseFlow/benchmarks/jouleThomsonCooling/jouleThomsonCooling.xml + :language: xml + :start-after: + :end-before: + +The mesh is only one cell thick in the angular and vertical directions, so the flow is purely radial. The two wells are vertical and are located in the mid-plane of the wedge (:math:`\theta = 0^{\circ}`): the injector is at :math:`r = 0.101` m, in the first cell of the mesh, and the producer is at :math:`r = 1129.5` m, in the outermost cell, which ends at the outer boundary at 1130 m. As in the paper, the production is therefore taken from the outermost grid block. + +------------------------------------------------------------------------ +Constitutive models +------------------------------------------------------------------------ + +**Fluid model** + +The fluid is a mixture of three components (CO2, CH4 and H2O) distributed in two phases (an aqueous phase named ``water`` and a ``gas`` phase) through a two-phase flash. It is defined by the ``CompositionalTwoPhaseThermalFluidPhillipsBrine`` model: + +.. literalinclude:: ../../../../../../../inputFiles/compositionalMultiphaseFlow/benchmarks/jouleThomsonCooling/jouleThomsonCooling.xml + :language: xml + :start-after: + :end-before: + :emphasize-lines: 4,10-19 + +The main ingredients of this model are the following: + +- ``equationsOfState`` lists one equation of state for each phase, in the order of ``phaseNames``. The gas phase uses the Peng-Robinson (PR) cubic equation of state, as in the paper. The aqueous phase uses the S\ |oslash|\ reide-Whitson variant of PR, which modifies the attractive term for water and uses salinity-dependent interaction parameters, in order to represent the mutual solubilities of the gases and of water. The aqueous phase stays at its residual saturation in this problem and acts mostly as a heat reservoir, so that the solubility effects are weak, but the model remains fully compositional. +- the aqueous phase density and viscosity follow the brine correlations of Phillips et al., with ``salinity`` equal to 0 (fresh water). +- ``pressureCoordinates`` and ``temperatureCoordinates`` define the pressure-temperature grid on which the fluid properties are evaluated. They must bracket all the states visited during the simulation: here from 1 bar (the minimum bottom-hole pressure of the producer) to 300 bar (the maximum bottom-hole pressure of the injector), and from -10 \ |degree|\ C to 110 \ |degree|\ C. The states actually visited are between about 5 and 7 MPa, and between 31 and 45 \ |degree|\ C. + +.. _jouleThomsonCoolingHeatCapacity: + +**Heat capacity coefficients** + +The enthalpy of a phase is the sum of an ideal-gas contribution and of a departure (residual) contribution given by the equation of state, which accounts for the non-ideal behavior of the real gas. The ideal-gas part is obtained by integrating the isobaric ideal-gas heat capacity of each component, which is provided by ``componentHeatCapacityCoefficients`` (one row per component, in the same order as ``componentNames``): + +.. literalinclude:: ../../../../../../../inputFiles/compositionalMultiphaseFlow/benchmarks/jouleThomsonCooling/jouleThomsonCooling.xml + :language: xml + :start-after: + :end-before: + :lines: 15-19 + :dedent: 6 + +Each row contains the five coefficients :math:`a_0, \ldots, a_4` of a polynomial in the temperature difference +with respect to :math:`T_{ref} = 298.15` K (25 \ |degree|\ C): + +.. math:: + + c_p^{0}(T) = \sum_{i=0}^{4} a_i \, (T - T_{ref})^i \qquad [\mathrm{J/mol/K}], \quad T \text{ in K}. + +The first coefficient :math:`a_0` is therefore the ideal-gas molar heat capacity at 25 \ |degree|\ C (37.14 J/mol/K for CO2, 35.71 for CH4 and 33.59 for H2O). The following table gives the resulting values at the reservoir temperature, which are in agreement with the reference ideal-gas heat capacities of the NIST equations of state (`Lemmon et al., 1998 `__) within 0.01 J/mol/K over the whole temperature range of the fluid tables (-10 to 110 \ |degree|\ C). + +.. list-table:: Ideal-gas heat capacity :math:`c_p^0` at 318.15 K + :widths: 20 20 20 20 + :header-rows: 1 + + * - Component + - CO2 + - CH4 + - H2O + * - :math:`c_p^0` (J/mol/K) + - 38.04 + - 36.51 + - 33.69 + +These coefficients are important for two reasons. First, the JT coefficient is inversely proportional to the heat capacity, and the heat capacity of the gas, which is the sum of the ideal-gas part and of the real-gas departure, controls how much the gas cools for a given enthalpy variation. Second, the heat capacity sets the thermal inertia of the fluids against that of the rock, which determines how fast the cooling is attenuated by the matrix. + +.. _jouleThomsonBinaryCoeff: + +.. **Binary interaction parameters** + +.. The Peng-Robinson equation of state uses a mixing rule for the attractive term of the mixture in which the cross term of each pair of components :math:`(i,j)` is scaled by :math:`(1 - k_{ij})`. The binary interaction parameters :math:`k_{ij}` are given by ``componentBinaryCoeff`` as a symmetric matrix with a zero diagonal, in the order CO2, CH4, H2O (the same order as ``componentNames``), so that :math:`k_{CO_2-CH_4} = 0.1`, :math:`k_{CO_2-H_2O} = 0.1896` and :math:`k_{CH_4-H_2O} = 0.4850`. The CO2-CH4 coefficient is the one that matters most for this problem, since the gas phase is a CO2-CH4 mixture in the region where the CO2 front progresses and displaces the methane. The coefficients involving water are larger, and they only affect the small amount of water vapor in the gas phase at 45 \ |degree|\ C and the solubility of the gases in the aqueous phase, which stays at its residual saturation. + +**Rock, relative permeability and thermal conductivity** + +The rock is modeled as an incompressible porous medium (``compressibility`` = 0 in ``rockPorosity``) with a constant permeability and a constant porosity of 0.3. Its heat capacity is given as +a volumetric value, :math:`\rho_R C_R = 2600 \times 1000 = 2.6 \times 10^{6}` J/m\ :sup:`3`/K, which is the product +of the grain density and the grain heat capacity prescribed by `(Oldenburg, 2007) `__. The matrix contributes :math:`(1-\phi)\rho_R C_R` per unit of bulk volume to the energy accumulation, which is the term that attenuates the JT cooling. + +.. literalinclude:: ../../../../../../../inputFiles/compositionalMultiphaseFlow/benchmarks/jouleThomsonCooling/jouleThomsonCooling.xml + :language: xml + :start-after: + :end-before: + +The relative permeabilities are Brooks-Corey power laws with the residual saturations given by `(Oldenburg, 2007) `__. With a water residual saturation of 0.2 +(``phaseMinVolumeFraction``) and an initial water saturation equal to 0.2, the water relative permeability is zero and the water does not move: + +.. literalinclude:: ../../../../../../../inputFiles/compositionalMultiphaseFlow/benchmarks/jouleThomsonCooling/jouleThomsonCooling.xml + :language: xml + :start-after: + :end-before: + +The thermal conductivity of the formation is constant and isotropic: + +.. literalinclude:: ../../../../../../../inputFiles/compositionalMultiphaseFlow/benchmarks/jouleThomsonCooling/jouleThomsonCooling.xml + :language: xml + :start-after: + :end-before: + +.. note:: + The *GEOS* model is not strictly identical to the TOUGH2/EOS7C model of `(Oldenburg, 2007) `__. The main differences are: + + - the Corey relative permeabilities of the paper are approximated by Brooks-Corey laws; + - the capillary pressure (van Genuchten) and the molecular diffusion of the paper (Table 2) are not included in the + *GEOS* input; + - the aqueous phase uses the S\ |oslash|\ reide-Whitson equation of state, whereas EOS7C uses PR for the whole system; + - the producer is a well with a prescribed surface gas rate, located at 1129.5 m in the outermost grid block, instead of a source term in that block. + +------------------------------------------------------------------------ +Flow solver +------------------------------------------------------------------------ + +The flow is solved by ``CompositionalMultiphaseFVM``, with a TPFA discretization, on the reservoir region only: + +.. literalinclude:: ../../../../../../../inputFiles/compositionalMultiphaseFlow/benchmarks/jouleThomsonCooling/jouleThomsonCooling.xml + :language: xml + :start-after: + :end-before: + +.. literalinclude:: ../../../../../../../inputFiles/compositionalMultiphaseFlow/benchmarks/jouleThomsonCooling/jouleThomsonCooling.xml + :language: xml + :start-after: + :end-before: + +The ``isThermal`` flag activates the energy balance equation, so that the temperature is an unknown of the problem. The +``useMass`` flag selects a mass-based formulation: the compositions and the well rates are expressed as mass +fractions and mass rates (kg/s). The time step is controlled by the maximum relative changes of pressure (20 %), temperature (2 %) +and phase volume fraction (20 %) per time step, and ``maxCompFractionChange`` limits the update of the compositions during +the Newton iterations. + +The reservoir and the wells are coupled in a fully implicit way by ``CompositionalMultiphaseReservoir``. The problem is small +(one-dimensional), so a direct linear solver is used: + +.. literalinclude:: ../../../../../../../inputFiles/compositionalMultiphaseFlow/benchmarks/jouleThomsonCooling/jouleThomsonCooling.xml + :language: xml + :start-after: + :end-before: + +------------------------------------------------------------------------ +Initial and boundary conditions +------------------------------------------------------------------------ + +.. _jouleThomsonInitialState: + +**Initial state** + +The initial state is uniform over the whole reservoir and is specified with ``FieldSpecification`` tags: + +.. literalinclude:: ../../../../../../../inputFiles/compositionalMultiphaseFlow/benchmarks/jouleThomsonCooling/jouleThomsonCooling.xml + :language: xml + :start-after: + :end-before: + +It is chosen to match the initial conditions of the paper, for the following reasons. + +- **Pressure, 50 bar.** This is the initial pressure of the depleted reservoir. It is low compared to the pressure of + the CO2 supply, which is what creates a large pressure drop around the injector. The values are uniform (no + ``HydrostaticEquilibrium`` is required) since the model is one cell thick and the problem is + controlled by the radial flow. +- **Temperature, 45 °C (318.15 K).** This is the reservoir temperature as prescribed, and it is the same as the injection + temperature. The initial thermal state is therefore in equilibrium with the injected fluid, and any temperature + deviation observed afterwards is the result of the JT expansion and of its attenuation by the heat of the rock, + not of a thermal contrast between the injected fluid and the reservoir. +- **Composition.** The reservoir contains methane and water, with no CO2. In *GEOS*, the saturations are not primary inputs: the initial phase distribution is computed by a flash from the global composition, pressure and temperature. The paper + prescribes an initial water saturation equal to the residual saturation, :math:`S_w = 0.2`, so that the water is immobile and the only mobile phase is gas. The global mass fractions are calibrated to obtain exactly this saturation. At 50 bar and 45 \ |degree|\ C, the density of CH4 gas is :math:`\rho_g = 32.36` kg/m\ :sup:`3` and the density of water is :math:`\rho_w = 992.3` kg/m\ :sup:`3`, so that + + .. math:: + + z_{CH_4} = \frac{\rho_g (1 - S_w)}{\rho_g (1 - S_w) + \rho_w S_w} = \frac{32.36 \times 0.8}{32.36 \times 0.8 + 992.3 \times 0.2} \approx 0.1154, + \qquad z_{H_2O} = 1 - z_{CH_4} \approx 0.8846, + + which are the values of ``initialCH4`` and ``initialH2O`` (``initialCO2`` is the CO2 mass fraction and is equal to zero). The mutual solubility between the phases is neglected in this estimate, and its effect on the saturation is very small. + +.. note:: + These mass fractions are specific to the initial pressure and temperature. If these values are modified (for example + to reproduce the cases of the paper with a reservoir at 75 \ |degree|\ C), the densities change and the compositions must be + computed again to keep :math:`S_w = 0.2`. + +**Boundary conditions** + +No boundary condition is imposed: all the outer boundaries of the domain are closed for the mass and for the heat +(no-flow and no heat flux). This corresponds to the closed reservoir of the paper, in which the pressure +increases with the injection and the heat is exchanged only between the fluids and the rock. The only sources and sinks +are the two wells. + +.. _jouleThomsonCoolingWellControls: + +**Well controls** + +The paper specifies constant-rate injection of CO2 and constant-rate production of CH4. These are the controls that govern the +pressure drop between the wells and the reservoir and therefore the JT cooling. The injector is defined as follows: + +.. literalinclude:: ../../../../../../../inputFiles/compositionalMultiphaseFlow/benchmarks/jouleThomsonCooling/jouleThomsonCooling.xml + :language: xml + :start-after: + :end-before: + +and the producer as: + +.. literalinclude:: ../../../../../../../inputFiles/compositionalMultiphaseFlow/benchmarks/jouleThomsonCooling/jouleThomsonCooling.xml + :language: xml + :start-after: + :end-before: + +The correspondence with the paper is the following: + +.. list-table:: Well controls and their relation to the paper + :widths: 20 35 45 + :header-rows: 1 + + * - Item + - *GEOS* input + - Relation to the paper + * - Injector control + - ``control="massRate"`` with ``InjectionMassRateConstraint``, ``massRate`` = 8.3333e-3 kg/s + - Constant-rate injection of 3 kg/s. The rate is divided by 360 because the model is a 1\ |degree|\ sector + (:math:`3 \times 1/360`). + * - Injected fluid + - ``injectionStream`` = {0.9999, 0.0001, 0}, ``injectionTemperature`` = 318.15 K + - Pure CO2 at :math:`T_{inj}` = 45 \ |degree|\ C, equal to the reservoir temperature. The stream is + 99.99 % CO2 (in mass), with a trace of CH4 and no H2O. + * - Producer control + - ``control="phaseVolRate"`` for the ``gas`` phase, :math:`2.28872\times 10^{-3}` m\ :sup:`3`/s at surface conditions + (101325 Pa, 288.15 K) + - Constant production of 0.56 kg/s of CH4. In the sector, this is :math:`0.56/360 = 1.556 \times 10^{-3}` kg/s, + which is converted to a volumetric rate using the density of methane at surface conditions, + :math:`0.68` kg/m\ :sup:`3`. + * - Well status + - ``statusTableName="wellStatus"`` + - The ``wellStatus`` table function keeps both wells closed for :math:`t < 0` and opens them at :math:`t = 0`, which is the start of the injection. + * - Crossflow + - ``enableCrossflow="0"`` for both wells + - The injector can only inject and the producer can only produce. + +The mass rate of the producer in the paper and the surface volumetric rate used here are equivalent as long as the +produced gas is methane. A simple volumetric balance indicates that, after 15 years, the CO2 front remains at about 500-600 m from the injector, so the producer, at 1129.5 m, only produces methane during the simulation, and the rate is +constant in mass. + +The ``writeCSV`` flag exports the time history of the rates, of the bottom-hole pressure and of the +temperature of each well, which is useful to monitor the injection overpressure. + +------------------------------------------------------------------------ +Events and output +------------------------------------------------------------------------ + +The simulation covers 15 years. The events are defined in the ``Events`` block, whose ``minTime`` and ``maxTime`` are -1 day and slightly more than 15 years, so that the last output event is reached. They are of five kinds: + +- **Solver applications.** The coupled reservoir and well solver is applied by two periodic events. During the first 10 s the maximum time step is 1 s, to handle the sudden start of the injection and the initial pressure transient near the well. After that, the maximum time step is 30 days (2.592e6 s). In both cases the actual time step is further controlled by the targets defined in the flow solver (maximum relative changes of pressure, temperature and phase volume fraction). +- **VTK outputs.** Eight solo events write the full fields at the times of the curves shown in Figure 7 of the paper: the initial state, then 1, 10, 20, 30 and 90 days, 1 year and 15 years. +- **Statistics.** The ``CompositionalMultiphaseStatistics`` task (pressure, temperature, phase and component masses, CFL numbers) is executed at the initial time and then once per year. +- **Time history.** Around each of the seven output times after the initial state, a periodic event triggers the ``PackCollection`` tasks that collect the cell centers, the pressure, the temperature and the phase compositions. The collection is hourly, over a 3-hour window starting one hour before the output time, and the collected data are written once by the ``timeHistory`` output at 15 years. These files provide the radial profiles used for the comparison with the paper. +- **Restart.** Restart files are written every 5 years. + +.. note:: + The time tags in *GEOS* are in seconds. The ``maxTime`` of the events is slightly larger than 15 years, so that the + last output event is reached. + +------------------------------------------------------------------------ +Post-processing and reporting +------------------------------------------------------------------------ + +The results of the simulation can be post-processed and visualized using scripts and data located next to this documentation page at + +.. code-block:: console + + src/docs/sphinx/advancedExamples/validationStudies/carbonStorage/jouleThomsonCooling + +The profiles of the CO2 mass fraction in the gas phase, of the pressure and of the temperature along the radial direction are post-processed with the Python script ``plot_joule_thomson_profile.py``. It reads the time history written by the ``timeHistory`` output (``jouleThomsonCooling.hdf5``), selects the output times of the simulation (1, 10, 20, 30 and 90 days, 1 year and 15 years), and plots the profiles of the first row of the wedge on a logarithmic radial axis from 0.1 to 1000 m. For the temperature, the script also overlays the points digitized from Figure 7 of `(Oldenburg, 2007) `__ (``oldenburg.json``), shown as dots with the same color as the corresponding *GEOS* curve. + +The main features to observe are: + +- a cold region around the injector whose temperature decreases and whose extent increases with time. After 15 years the cooling is + about 14 \ |degree|\ C with respect to the initial 45 \ |degree|\ C (minimum temperature around 31 \ |degree|\ C), and the cooled region extends over + more than 100 m; +- a pressure increase near the well, of about 22 bar after 15 years, which decreases quickly with the distance to the well; +- a cooled and pressurized region that remains inside the region where the gas is almost pure CO2; +- a temperature that remains at 45 \ |degree|\ C farther from the well. + +**CO2 front** + +The CO2 mass fraction in the gas phase is the quantity that locates the injected CO2: it is equal to 1 behind the front and to 0 ahead of it, where the gas is methane. The front is already at about 9 m from the well after 1 day. It reaches about 25 m after 10 days and 550-600 m after 15 years, which is about half of the domain. The producer, at 1129.5 m, therefore only produces methane during the whole simulation, as assumed in the definition of the well controls. + +.. note:: + The front is sharp, with a transition from 1 to 0 over a factor of about 1.5 in radius. In the absence of molecular diffusion and dispersion in the *GEOS* model, this width is controlled by the numerical dispersion of the scheme and by the mesh size, and the paper, which includes molecular diffusion, can be expected to show a somewhat wider transition. The CO2 profile is important for the interpretation of the two other profiles, since the viscosity and the density of the gas change across the front. + +.. plot:: docs/sphinx/advancedExamples/validationStudies/carbonStorage/jouleThomsonCooling/plot_joule_thomson_profile.py plot_co2_fraction + :align: center + :caption: Radial profiles of the CO2 mass fraction in the gas phase at the output times. The sharp front moves from about 9 m (1 day) to 550-600 m (15 years). + +**Pressure** + +The pressure build-up is largest at the injector and is concentrated in the first tens of meters at early times. At the first cell (:math:`r \approx 0.1` m), the pressure is about 64 bar after 1 day, which is an increase of about 14 bar with respect to the initial pressure of 50 bar, and it rises slowly reaching 72 bar after 15 years (an increase of about 22 bar). In the region filled by CO2, the pressure decreases almost linearly with :math:`\log r`, with a drop of about 4-5 bar per decade of radius at 15 years. + +The pressure perturbation travels much faster than the CO2 front, since it propagates by diffusion of the pressure (the compressibility of the gas is the only storage mechanism) while the front is transported with the fluid. After 1 day, the pressure is already perturbed up to about 30-50 m while the front is at 9 m, and after 1 year the whole domain is affected. + +At the outer boundary, the pressure first decreases slightly below the initial value (about 49.7 bar after 1 year) because of the production, before the pressure wave of the injector arrives. It then increases to about 53 bar after 15 years: the system is closed and the injected mass (3 kg/s) is larger than the produced mass (0.56 kg/s), so that the average pressure of the reservoir increases. The pressure at the well is far below the 300 bar limit of the injector, which confirms that the injection stays at a constant rate. + +.. plot:: docs/sphinx/advancedExamples/validationStudies/carbonStorage/jouleThomsonCooling/plot_joule_thomson_profile.py plot_pressure + :align: center + :caption: Radial profiles of the pressure at the output times. The pressure is highest at the injector (about 72 bar after 15 years) and tends to the initial value of 50 bar far from the well at early times. + +**Temperature and comparison with the paper** + +The temperature profiles show the three zones that are expected for this problem. Near the well ($r=0.1$ m), the injected gas has not yet expanded significantly, maintaining an elevated temperature of approximately 43 \ |degree|\ C. Moving outward, JT expansion cools the gas and creates a temperature minimum. This minimum deepens and shifts farther from the well over time: after 1 day, the temperature drops to about 38.5 \ |degree|\ C near 0.7 m. After 15 years, the minimum falls to roughly 31 \ |degree|\ C at 30-40 m. Farther into the reservoir, the temperature recovers to the undisturbed baseline of 45 \ |degree|\ C. + +The cold region is much smaller than the region invaded by the CO2: after 15 years the temperature is back to about 45 \ |degree|\ C beyond 100-150 m, while the CO2 front is at 550-600 m. The thermal front is retarded with respect to the gas front by the heat capacity of the rock and of the immobile water, which have to be cooled before the cooling can move forward. This is the transient effect mentioned above: the matrix attenuates the cooling, but the thermal front keeps propagating, and the minimum temperature keeps decreasing because the pressure drop that drives the JT cooling increases with time. The weak undershoot below 45 \ |degree|\ C (a few tenths of a degree) at 15 years between the thermal front and the CO2 front is the cooling due to the small pressure drop in this region. + +The magnitude of the minimum temperature is larger than a simple estimate based on the JT coefficient at the initial conditions. The pressure difference between the well and the location of the minimum at 15 years is about 10 bar (72 bar at the well, 61-62 bar at 30-40 m), so that a cooling of 14 \ |degree|\ C requires an apparent JT coefficient of at least 1.4 K/bar, larger than the 0.9 K/bar of CO2 at 45 \ |degree|\ C and 50 bar. This is consistent with the strong increase of :math:`\mu_{JT}` when CO2 approaches its critical point (31 \ |degree|\ C and 73.8 bar), which is the case in the cold region, and it means that the cooling is sensitive to the thermodynamic model. + +.. plot:: docs/sphinx/advancedExamples/validationStudies/carbonStorage/jouleThomsonCooling/plot_joule_thomson_profile.py plot_temperature + :align: center + :caption: Radial profiles of the temperature at the output times (lines: *GEOS*; dots: results of Figure 7 of Oldenburg (2007), digitized). There is no reference curve for 1 day. The minimum temperature after 15 years (about 31 \ |degree|\ C at 30-40 m) is reproduced, while *GEOS* is warmer than the paper close to the well (r < 10 m). + +The comparison with the points taken from `(Oldenburg, 2007) `__ (dots) shows a good overall agreement, with some differences that are worth noting: + +- **Good agreement on the cold front and on the minimum temperature.** The position and the steepness of the thermal front are close to the digitized points at all times (for instance the rise from 35 to 45 \ |degree|\ C between about 10 and 30 m after 1 year, and from 35 to 40 \ |degree|\ C between about 60 and 100 m after 15 years), and the minimum temperature after 15 years, about 31 \ |degree|\ C at 30-40 m, is the same as in the paper. This minimum is the quantity that matters for the risk of hydrate or ice formation, and it remains far above the temperature of hydrate stability. The thermal front and the CO2 front also have the same relative position as in the paper (the cooled zone stays inside the CO2-rich region). +- **Warmer temperatures in the near-well region.** For :math:`r` smaller than about 10 m, the digitized points are lower than the *GEOS* curves. At early times (10 to 90 days), the paper gives 34-35 \ |degree|\ C at 1.5-3 m, about 2 \ |degree|\ C below the *GEOS* curves, and after 15 years the paper gives 35.5 \ |degree|\ C at 3 m and 33.8 \ |degree|\ C at 7 m, while *GEOS* is about 3-4 \ |degree|\ C warmer. The paper also shows a minimum that is more flat and more extended (about 31-34 \ |degree|\ C between 7 and 40 m at 15 years), whereas the *GEOS* profile decreases almost linearly with :math:`\log r` up to about 10 m and then drops more rapidly toward the minimum. *GEOS* therefore predicts a weaker cooling than the paper in the first meters around the injector. + +.. note:: + The cooling can be compared with the upper bound :math:`\mu_{JT}\,\Delta P`. The relevant pressure difference is the one between the well and the point of interest, and not the total pressure increase at the well, since the gas is injected at the reservoir temperature and cools as it expands toward lower pressures. As discussed above, the cooling observed in the simulation corresponds to an apparent JT coefficient larger than the one at the initial conditions, because of the proximity of the critical point, even though the heat released by the matrix, which is important in the first days and progressively exhausted, attenuates the cooling. + +------------------------------------------------------------------ +To go further +------------------------------------------------------------------ + +The other cases of the paper can be obtained from this input by changing a few parameters: + +- the permeability (``rockPerm``): :math:`10^{-13}` and :math:`10^{-14}` m\ :sup:`2` for Figures 5 and 6, and the summary of Figure 8; + the injection rate is kept unchanged, so the pressure drop and the JT cooling vary with the permeability; +- the porosity (``rockPorosity``), to study the compensation between the permeability and the heat supplied by the matrix (Figure 9); +- the injection temperature and the initial reservoir temperature (``injectionTemperature``, ``initialTemperature``) for the + cases of Figures 10 and 11 (20 \ |degree|\ C injection in a 75 \ |degree|\ C reservoir). The initial compositions must then be computed again, as explained above. + +**Feedback on this example** + +For any feedback on this example, please submit a `GitHub issue on the project's GitHub page `_. + +------------------------------------------------------------------ +References +------------------------------------------------------------------ + +- Oldenburg, C. M. 2007. Joule-Thomson cooling due to CO2 injection into natural gas reservoirs. *Energy Conversion and Management*, 48(6), 1808-1815. https://doi.org/10.1016/j.enconman.2007.01.010 +- Lemmon, E. W., McLinden, M. O., and Friend, D. G. 1998. Thermophysical Properties of Fluid Systems. In *NIST Chemistry WebBook*, NIST Standard Reference Database, eds. Linstrom, P. J. and Mallard, W. G., Vol. 69. Gaithersburg, Maryland, USA: National Institute of Standards and Technology. https://webbook.nist.gov/chemistry/fluid/ diff --git a/src/docs/sphinx/advancedExamples/validationStudies/carbonStorage/jouleThomsonCooling/jouleThomsonCooling.hdf5 b/src/docs/sphinx/advancedExamples/validationStudies/carbonStorage/jouleThomsonCooling/jouleThomsonCooling.hdf5 new file mode 100644 index 00000000000..92428786cbb Binary files /dev/null and b/src/docs/sphinx/advancedExamples/validationStudies/carbonStorage/jouleThomsonCooling/jouleThomsonCooling.hdf5 differ diff --git a/src/docs/sphinx/advancedExamples/validationStudies/carbonStorage/jouleThomsonCooling/oldenburg.json b/src/docs/sphinx/advancedExamples/validationStudies/carbonStorage/jouleThomsonCooling/oldenburg.json new file mode 100644 index 00000000000..c893d7fb4d3 --- /dev/null +++ b/src/docs/sphinx/advancedExamples/validationStudies/carbonStorage/jouleThomsonCooling/oldenburg.json @@ -0,0 +1 @@ +{"473364000": {"p": [[0.0, 5.0, 10.0, 20.0, 40.0, 60.0, 80.0, 100.0], [69.14, 66.44000000000001, 64.19000000000001, 62.25000000000001, 61.050000000000004, 60.00000000000001, 59.25000000000001, 58.95]], "t": [[0.0, 3.0, 7.0, 16.5, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0], [40.0, 35.5, 33.8, 31.7, 31.0, 31.0, 32.0, 33.0, 34.2, 35.6, 37.6, 40.0]], "y": [[0.0, 60.0, 80.0, 100.0], [1.0, 1.0, 0.97, 0.9]]}, "31557600": {"p": [[0.0, 10.0, 20.0, 40.0, 60.0, 80.0, 100.0], [65.39, 61.650000000000006, 59.25000000000001, 56.7, 55.510000000000005, 54.760000000000005, 54.160000000000004]], "t": [[0.0, 3.0, 7.0, 12.0, 15.0, 20.0, 26.0, 33.0, 40.0, 50.0, 60.0, 100.0], [40.0, 35.0, 33.5, 33.5, 35.0, 37.5, 41.4, 43.3, 44.4, 44.9, 45.0, 45.0]], "y": [[0.0, 100.0], [1.0, 0.99]]}, "7776000": {"p": [[0.0, 10.0, 20.0, 40.0, 60.0, 80.0, 100.0], [64.64, 59.400000000000006, 57.900000000000006, 55.660000000000004, 54.31, 53.56, 52.96]], "t": [[0.0, 2.0, 5.0, 8.0, 11.0, 15.0, 20.0, 25.0, 100.0], [40.0, 34.5, 34.9, 38.0, 40.5, 43.8, 44.7, 45.0, 45.0]], "y": [[0.0, 30.0, 40.0, 50.0, 60.0, 80.0, 100.0], [1.0, 1.0, 0.99, 0.92, 0.86, 0.59, 0.21]]}, "2592000": {"p": [[0.0, 10.0, 20.0, 40.0, 60.0, 80.0, 100.0], [63.900000000000006, 58.35, 56.85, 54.910000000000004, 53.410000000000004, 52.660000000000004, 52.06]], "t": [[0.0, 1.5, 4.0, 7.0, 10.0, 15.0, 20.0, 100.0], [40.0, 34.0, 35.5, 40.5, 43.0, 44.5, 45.0, 45.0]], "y": [[0.0, 15.0, 25.0, 40.0, 55.0, 65.0, 75.0, 80.0, 100.0], [1.0, 1.0, 0.88, 0.66, 0.25, 0.07, 0.02, 0.0, 0.0]]}, "1728000": {"p": [[0.0, 10.0, 20.0, 40.0, 60.0, 80.0, 100.0], [63.45, 57.75000000000001, 56.1, 54.010000000000005, 52.81, 52.06, 51.61000000000001]], "t": [[0.0, 1.5, 4.0, 7.0, 10.0, 15.0, 20.0, 100.0], [40.0, 34.5, 36.5, 41.0, 43.5, 44.8, 45.0, 45.0]], "y": [[0.0, 15.0, 25.0, 40.0, 50.0, 60.0, 70.0, 100.0], [1.0, 1.0, 0.74, 0.4, 0.2, 0.04, 0.0, 0.0]]}, "864000": {"p": [[0.0, 10.0, 20.0, 40.0, 60.0, 80.0, 100.0], [63.150000000000006, 56.85, 55.21000000000001, 53.260000000000005, 52.06, 51.46, 50.86000000000001]], "t": [[0.0, 1.5, 4.0, 7.0, 12.0, 100.0], [40.5, 34.5, 38.0, 44.3, 45.0, 45.0]], "y": [[0.0, 12.0, 18.0, 27.0, 40.0, 50.0, 60.0, 100.0], [1.0, 1.0, 0.86, 0.45, 0.06, 0.02, 0.0, 0.0]]}} \ No newline at end of file diff --git a/src/docs/sphinx/advancedExamples/validationStudies/carbonStorage/jouleThomsonCooling/plot_joule_thomson_profile.py b/src/docs/sphinx/advancedExamples/validationStudies/carbonStorage/jouleThomsonCooling/plot_joule_thomson_profile.py new file mode 100644 index 00000000000..361ba0b29a2 --- /dev/null +++ b/src/docs/sphinx/advancedExamples/validationStudies/carbonStorage/jouleThomsonCooling/plot_joule_thomson_profile.py @@ -0,0 +1,215 @@ +import h5py +import numpy +import json +import os +import argparse +import matplotlib.pyplot as plt +from matplotlib.lines import Line2D +from matplotlib.patches import Patch +import matplotlib.cm as cm +import typing + +def get_selected_times() -> dict: + """ + Defines the specific simulation times for which data will be extracted and plotted. + + Returns: + dict: A dictionary mapping time in seconds (int) to a human-readable label (str). + """ + DAY_2_SEC = 86400 + YEAR_2_SEC = int(DAY_2_SEC * 365.25 + 0.5) + return { + DAY_2_SEC: '1 day', + 10*DAY_2_SEC: '10 days', + 20*DAY_2_SEC: '20 days', + 30*DAY_2_SEC: '30 days', + 90*DAY_2_SEC: '90 days', + YEAR_2_SEC: '1 year', + 15*YEAR_2_SEC: '15 years', + } + +def get_time_indices(time_data: numpy.ndarray, selected_times: dict) -> dict: + """ + Finds the array indices closest to the requested target times. + + Args: + time_data (numpy.ndarray): The 1D array of time steps from the simulation data. + selected_times (dict): The dictionary of target times (keys are in seconds). + + Returns: + dict: A dictionary mapping the target time (in seconds) to its corresponding index in the data arrays. + """ + indices = {} + for time in selected_times: + # Find the index of the value in time_data that is closest to 'time' + indices[time] = numpy.absolute(time_data - time).argmin() + return indices + +def read_data(file_name: str, selected_times: dict) -> dict: + """ + Reads simulation results from an HDF5 file for specified time steps. + + Args: + file_name (str): The path to the HDF5 file containing the simulation data. + selected_times (dict): The dictionary defining the target times to extract. + + Returns: + dict: A nested dictionary containing the extracted data. + The top-level keys are the target times in seconds. + Each value is a dictionary containing arrays for radius ('r'), + pressure ('p'), temperature ('t'), and CO2 mass fraction ('y'). + """ + data = {} + with h5py.File(file_name, 'r') as h5: + # Assume that all the time arrays are the same, use elementCenter Time as reference + element_center_time = h5['elementCenter Time'][()].flatten() + time_indices = get_time_indices(element_center_time, selected_times) + + # Read the full arrays for element centers and properties + element_center = h5['elementCenter'][()] + pressure = h5['pressure'][()] + temperature = h5['temperature'][()] + phase_comp_fraction = h5['fluid_phaseCompFraction'][()] + + # Extract data for the closest matching time step for each selected time + for time in selected_times: + index = time_indices[time] + # Assuming radial geometry, take the first component (x-coordinate) as radius + r = element_center[index, :, 0] + # Convert Temperature from Kelvin to Celsius + t = temperature[index, :] - 273.15 + # Convert Pressure from Pascals to Bar + p = pressure[index, :] * 1.0e-5 + # Extract CO2 Mass Fraction (assuming it's the 4th component, index 3) + y = phase_comp_fraction[index, :, 3] + data[time] = {"r": r, "p": p, "t": t, "y": y} + + return data + +def plot_profile_data(item: str, data: dict, selected_times: dict, reference: dict = None): + """ + Plots the spatial profile of a specified property at different time steps. + + Args: + item (str): The property to plot ('p' for pressure, 't' for temperature, 'y' for mass fraction). + data (dict): The dictionary containing the simulation data extracted by read_data(). + selected_times (dict): The dictionary mapping times to labels, used for the legend. + reference (dict, optional): Reference data to plot alongside simulation results. Defaults to None. + """ + # Set global font size for matplotlib plots + plt.rcParams.update({'font.size': 16}) + + # Setup Colormap to distinguish different time steps + cmap = plt.colormaps['jet'] + num_times = len(data.keys()) + colors = [cmap((i+0.5)/num_times) for i in range(num_times)] + + # Configuration for different plot types + figs = { + "p" : { "t": 'Pressure', "u": 'bar', "l": (45, 75) }, + "t" : { "t": 'Temperature', "u": '°C', "l": (30, 50) }, + "y" : { "t": 'CO2 Mass Fraction', "u": '-' } + } + + # Check if the requested item is valid + if item not in figs: + print(f"Warning: Property '{item}' is not configured for plotting. Skipping.") + return + + fig, ax = plt.subplots(figsize=(9, 5)) + sorted_times = [t for t in sorted(data.keys())] + + # Plot data for each time step + for ti, time in enumerate(sorted_times): + r = data[time]["r"] + # Plot simulation data as a solid line + ax.plot(r, data[time][item], color=colors[ti], linewidth=2) + + # Plot reference data as points if available + if reference and str(time) in reference: + r0, v0 = reference[str(time)][item] + ax.plot(r0, v0, color=colors[ti], linewidth=0, linestyle="None", marker="o") + + # Create custom legend handles using Patch to match colors + custom_handles = [] + for ti, time in enumerate(sorted_times): + time_label = selected_times[time] + custom_handles.append(Patch(facecolor=colors[ti], edgecolor='black', label=time_label)) + + # Configure axes labels, title, and scales + ax.set_xlabel('r (m)') + ax.set_ylabel("%s (%s)" % (figs[item]["t"], figs[item]["u"])) + ax.set_title("%s Profile" % (figs[item]["t"])) + ax.set_xscale('log') + ax.set_xlim(0.1, 1000) + + # Apply y-axis limits if specified in the configuration + if "l" in figs[item]: + ax.set_ylim(figs[item]["l"]) + + ax.grid(True) + ax.legend(handles=custom_handles, loc='best', fontsize=10) + fig.tight_layout() + + plt.show(block=True) + +def plot_property(property: str, reference: bool = False): + # The data files are located next to this script + # (the Sphinx plot directive runs the script from its own directory and does not define __file__) + data_directory = os.path.dirname(os.path.abspath(__file__)) if '__file__' in globals() else os.getcwd() + + # Get target times and extract simulation data + selected_times = get_selected_times() + file_name = os.path.join(data_directory, "jouleThomsonCooling.hdf5") + + try: + data = read_data(file_name, selected_times) + except FileNotFoundError: + print(f"Error: Could not find data file '{file_name}'. Please ensure the file exists.") + exit(1) + + # Load reference data if requested + reference_data = None + if reference: + ref_file_name = os.path.join(data_directory, "oldenburg.json") + try: + with open(ref_file_name, mode="r") as f: + reference_data = json.load(f) + except FileNotFoundError: + print(f"Warning: Reference file '{ref_file_name}' not found. Plotting without reference data.") + + # Iterate through requested properties and generate plots + props_to_plot = property.lower() + for prop in ['p', 't', 'y']: + if prop in props_to_plot: + plot_profile_data(prop, data, selected_times, reference=reference_data) + +def plot_pressure(): + plot_property("p") + +def plot_temperature(): + plot_property("t", True) + +def plot_co2_fraction(): + plot_property("y") + +if __name__ == "__main__": + # Setup argparse for command line arguments + parser = argparse.ArgumentParser(description="Plot Joule-Thomson cooling profiles from HDF5 data.") + + parser.add_argument( + "-p", "--property", + type=str, + default="pty", + help="Properties to plot. Include 'p' for pressure, 't' for temperature, 'y' for CO2 mass fraction. (Default: 'pty')" + ) + + parser.add_argument( + "-r", "--reference", + action="store_true", + help="Flag to load and show reference data from 'oldenburg.json'. (Default: False)" + ) + + args = parser.parse_args() + + #plot_property(args.property, args.reference)