OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Added dummy water phase to compositional simulator.

Browse Source 
Note: compositional is now always three phase!
This commit is contained in:
Svenn Tveit 2025-01-08 13:08:17 +01:00
parent 51edfd6f95
commit df6bdbe700
8 changed files with 92 additions and 32 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
flowexperimental/comp/flowexp_comp.hpp
View File

@ -20,7 +20,7 @@
#define FLOWEXP_COMP_HPP
#include <opm/material/constraintsolvers/PTFlash.hpp>
#include <opm/material/fluidsystems/GenericOilGasFluidSystem.hpp>
#include <opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp>
#include <opm/models/discretization/common/baseauxiliarymodule.hh>
#include <opm/models/ptflash/flashmodel.hh>
@ -169,7 +169,7 @@ private:
static constexpr int num_comp = getPropValue<TypeTag, Properties::NumComp>();
public:
using type = Opm::GenericOilGasFluidSystem<Scalar, num_comp>;
using type = Opm::GenericOilGasWaterFluidSystem<Scalar, num_comp>;
};
template<class TypeTag, int NumComp>
struct EnableMech<TypeTag, TTag::FlowExpCompProblem<NumComp>> {

14
opm/models/ptflash/flashindices.hh
View File

@ -51,13 +51,16 @@ class FlashIndices
using EnergyIndices = Opm::EnergyIndices<PVOffset + numComponents, enableEnergy>;
public:
static constexpr bool waterEnabled = false;
//! All phases active (note: immiscible/"dummy" water phase)
static constexpr bool waterEnabled = true;
static constexpr bool gasEnabled = true;
static constexpr bool oilEnabled = true;
static constexpr int waterPhaseIdx = -1;
static constexpr int numPhases = 2;
//! number of active phases
static constexpr int numPhases = 3;
//! number of equations/primary variables
static const int numEq = numComponents + EnergyIndices::numEq_;
static const int numEq = numComponents + EnergyIndices::numEq_ + 1;
// Primary variable indices
@ -67,6 +70,9 @@ public:
//! Index of the molefraction of the first component
static constexpr int z0Idx = pressure0Idx + 1;
//! Index of water saturation
static constexpr int water0Idx = z0Idx + numComponents - 1;
// equation indices
//! Index of the mass conservation equation for the first

21
opm/models/ptflash/flashintensivequantities.hh
View File

@ -76,6 +76,7 @@ class FlashIntensiveQuantities
enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
enum { dimWorld = GridView::dimensionworld };
enum { pressure0Idx = Indices::pressure0Idx };
enum { water0Idx = Indices::water0Idx};
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
@ -215,20 +216,22 @@ public:
(R * fluidState_.temperature(FluidSystem::gasPhaseIdx));
// Update saturation
// \Note: the current implementation assume oil-gas system.
// Update oil/gas saturation; water saturation is a primary variable
Evaluation Sw = priVars.makeEvaluation(water0Idx, timeIdx);
Evaluation L = fluidState_.L();
Evaluation So = Opm::max((L * Z_L / ( L * Z_L + (1 - L) * Z_V)), 0.0);
Evaluation Sg = Opm::max(1 - So, 0.0);
Scalar sumS = Opm::getValue(So) + Opm::getValue(Sg);
Evaluation So = Opm::max((1 - Sw) * (L * Z_L / ( L * Z_L + (1 - L) * Z_V)), 0.0);
Evaluation Sg = Opm::max(1 - So - Sw, 0.0);
Scalar sumS = Opm::getValue(So) + Opm::getValue(Sg) + Opm::getValue(Sw);
So /= sumS;
Sg /= sumS;
Sw /= sumS;
fluidState_.setSaturation(0, So);
fluidState_.setSaturation(1, Sg);
fluidState_.setSaturation(FluidSystem::oilPhaseIdx, So);
fluidState_.setSaturation(FluidSystem::gasPhaseIdx, Sg);
fluidState_.setSaturation(FluidSystem::waterPhaseIdx, Sw);
fluidState_.setCompressFactor(0, Z_L);
fluidState_.setCompressFactor(1, Z_V);
fluidState_.setCompressFactor(FluidSystem::oilPhaseIdx, Z_L);
fluidState_.setCompressFactor(FluidSystem::gasPhaseIdx, Z_V);
// Print saturation
if (flashVerbosity >= 5) {

28
opm/models/ptflash/flashlocalresidual.hh
View File

@ -49,12 +49,16 @@ class FlashLocalResidual: public GetPropType<TypeTag, Properties::DiscLocalResid
using Indices = GetPropType<TypeTag, Properties::Indices>;
using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
enum { water0Idx = Indices::water0Idx };
enum { conti0EqIdx = Indices::conti0EqIdx };
enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
using DiffusionModule = Opm::DiffusionModule<TypeTag, enableDiffusion>;
@ -77,7 +81,16 @@ public:
const IntensiveQuantities& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
const auto& fs = intQuants.fluidState();
// compute storage term of all components within all phases
// compute water storage term
if (phaseIdx == waterPhaseIdx) {
unsigned eqIdx = conti0EqIdx + numComponents;
storage[eqIdx] =
Toolbox::template decay<LhsEval>(fs.density(phaseIdx))
* Toolbox::template decay<LhsEval>(fs.saturation(phaseIdx))
* Toolbox::template decay<LhsEval>(intQuants.porosity());
}
else {
// compute storage term of all components within oil/gas phases
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
unsigned eqIdx = conti0EqIdx + compIdx;
storage[eqIdx] +=
@ -86,6 +99,7 @@ public:
* Toolbox::template decay<LhsEval>(fs.saturation(phaseIdx))
* Toolbox::template decay<LhsEval>(intQuants.porosity());
}
}
EnergyModule::addPhaseStorage(storage, elemCtx.intensiveQuantities(dofIdx, timeIdx), phaseIdx);
}
@ -146,22 +160,34 @@ public:
up.fluidState().density(phaseIdx)
* extQuants.volumeFlux(phaseIdx);
if (phaseIdx == waterPhaseIdx) {
unsigned eqIdx = conti0EqIdx + numComponents;
flux[eqIdx] = tmp;
}
else {
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
flux[conti0EqIdx + compIdx] +=
tmp*up.fluidState().massFraction(phaseIdx, compIdx);
}
}
}
else {
Evaluation tmp =
Toolbox::value(up.fluidState().density(phaseIdx))
* extQuants.volumeFlux(phaseIdx);
if (phaseIdx == waterPhaseIdx) {
unsigned eqIdx = conti0EqIdx + numComponents;
flux[eqIdx] = tmp;
}
else {
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
flux[conti0EqIdx + compIdx] +=
tmp*Toolbox::value(up.fluidState().massFraction(phaseIdx, compIdx));
}
}
}
}
EnergyModule::addAdvectiveFlux(flux, elemCtx, scvfIdx, timeIdx);
}

6
opm/models/ptflash/flashnewtonmethod.hh
View File

@ -125,6 +125,12 @@ protected:
for (unsigned compIdx = 0; compIdx < numComponents - 1; ++compIdx) {
clampValue_(nextValue[z0Idx + compIdx], tol, 1-tol);
}
// limit change in water saturation to 0.2
constexpr Scalar dSwMax = 0.2;
if (update[Indices::water0Idx] > dSwMax) {
nextValue[Indices::water0Idx] = currentValue[Indices::water0Idx] - dSwMax;
}
}
private:
void clampValue_(Scalar& val, Scalar minVal, Scalar maxVal) const

16
opm/models/ptflash/flashprimaryvariables.hh
View File

@ -61,6 +61,12 @@ class FlashPrimaryVariables : public FvBasePrimaryVariables<TypeTag>
// primary variable indices
enum { z0Idx = Indices::z0Idx };
enum { pressure0Idx = Indices::pressure0Idx };
enum { water0Idx = Indices::water0Idx };
// phase indices
enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
@ -108,10 +114,15 @@ public:
// the energy module
EnergyModule::setPriVarTemperatures(*this, fluidState);
// assign components total fraction
for (int i = 0; i < numComponents - 1; ++i)
(*this)[z0Idx + i] = getValue(fluidState.moleFraction(i));
(*this)[pressure0Idx] = getValue(fluidState.pressure(0));
// assign pressure
(*this)[pressure0Idx] = getValue(fluidState.pressure(oilPhaseIdx));
// assign water saturation
(*this)[water0Idx] = getValue(fluidState.saturation(waterPhaseIdx));
}
/*!
@ -121,12 +132,13 @@ public:
*/
void print(std::ostream& os = std::cout) const
{
os << "(p_" << FluidSystem::phaseName(0) << " = "
os << "(p_" << FluidSystem::phaseName(FluidSystem::oilPhaseIdx) << " = "
<< this->operator[](pressure0Idx);
for (unsigned compIdx = 0; compIdx < numComponents - 2; ++compIdx) {
os << ", z_" << FluidSystem::componentName(compIdx) << " = "
<< this->operator[](z0Idx + compIdx);
}
os << ", S_w = " << this->operator[](water0Idx);
os << ")" << std::flush;
}
};

7
opm/simulators/flow/FlowProblemComp.hpp
View File

@ -362,6 +362,9 @@ public:
Dune::FieldVector<Scalar, numComponents> z(0.0);
Scalar sumMoles = 0.0;
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
if (phaseIdx == waterPhaseIdx){
continue;
}
const auto saturation = fs.saturation(phaseIdx);
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
Scalar tmp = fs.molarity(phaseIdx, compIdx) * saturation;
@ -527,6 +530,10 @@ protected:
- waterSaturationData[dofIdx]
- gasSaturationData[dofIdx]);
}
if (water_active) {
dofFluidState.setSaturation(FluidSystem::waterPhaseIdx,
waterSaturationData[dofIdx]);
}
//////
// set phase pressures

4
opm/simulators/flow/GenericOutputBlackoilModule.cpp
View File

@ -30,7 +30,7 @@
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
#include <opm/material/fluidsystems/BlackOilDefaultIndexTraits.hpp>
#include <opm/material/fluidsystems/GenericOilGasFluidSystem.hpp>
#include <opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp>
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/EclipseState/Runspec.hpp>
@ -1677,7 +1677,7 @@ INSTANTIATE_TYPE(float)
#endif
#define INSTANTIATE_COMP(NUM) \
template<class T> using FS##NUM = GenericOilGasFluidSystem<T, NUM>; \
template<class T> using FS##NUM = GenericOilGasWaterFluidSystem<T, NUM>; \
template class GenericOutputBlackoilModule<FS##NUM<double>>;
INSTANTIATE_COMP(0) // \Note: to register the parameter ForceDisableFluidInPlaceOutput