OilfieldToolsNet/opm-common
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

fluid systems: rename unexpressive phase and component indices

Browse Source 
These index names have been fully fluid system dependent for a while
and are supposed to be just used for convenience. This means that
phase names are now actual camelCase words.
This commit is contained in:
Andreas Lauser
2014-04-03 16:34:23 +02:00
parent 96b3d753e2
commit c0a7947357
27 changed files with 617 additions and 618 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
opm/material/binarycoefficients/Brine_CO2.hpp
View File

@@ -41,8 +41,8 @@ class Brine_CO2 {
typedef Opm::H2O<Scalar> H2O;
typedef Opm::CO2<Scalar, CO2Tables> CO2;
typedef Opm::IdealGas<Scalar> IdealGas;
static const int lPhaseIdx = 0; // index of the liquid phase
static const int gPhaseIdx = 1; // index of the gas phase
static const int liquidPhaseIdx = 0; // index of the liquid phase
static const int gasPhaseIdx = 1; // index of the gas phase
public:
/*!
@@ -119,7 +119,7 @@ public:
// if only liquid phase is present the mole fraction of CO2 in brine is given and
// and the virtual equilibrium mole fraction of water in the non-existing gas phase can be estimated
// with the mutual solubility function
if (knownPhaseIdx == lPhaseIdx) {
if (knownPhaseIdx == liquidPhaseIdx) {
ygH2O = A * (1 - xlCO2 - x_NaCl);
}
@@ -127,7 +127,7 @@ public:
// if only gas phase is present the mole fraction of water in the gas phase is given and
// and the virtual equilibrium mole fraction of CO2 in the non-existing liquid phase can be estimated
// with the mutual solubility function
if (knownPhaseIdx == gPhaseIdx) {
if (knownPhaseIdx == gasPhaseIdx) {
//y_H2o = fluidstate.
xlCO2 = 1 - x_NaCl - ygH2O / A;
}

42
opm/material/fluidmatrixinteractions/BrooksCorey.hpp
View File

@@ -92,8 +92,8 @@ public:
template <class Container, class FluidState>
static void capillaryPressures(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wPhaseIdx] = 0.0; // reference phase
values[Traits::nPhaseIdx] = pcnw(params, fs);
values[Traits::wettingPhaseIdx] = 0.0; // reference phase
values[Traits::nonWettingPhaseIdx] = pcnw(params, fs);
}
/*!
@@ -103,8 +103,8 @@ public:
template <class Container, class FluidState>
static void saturations(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wPhaseIdx] = Sw(params, fs);
values[Traits::nPhaseIdx] = 1 - values[Traits::wPhaseIdx];
values[Traits::wettingPhaseIdx] = Sw(params, fs);
values[Traits::nonWettingPhaseIdx] = 1 - values[Traits::wettingPhaseIdx];
}
/*!
@@ -120,8 +120,8 @@ public:
template <class Container, class FluidState>
static void relativePermeabilities(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wPhaseIdx] = krw(params, fs);
values[Traits::nPhaseIdx] = krn(params, fs);
values[Traits::wettingPhaseIdx] = krw(params, fs);
values[Traits::nonWettingPhaseIdx] = krn(params, fs);
}
/*!
@@ -134,10 +134,10 @@ public:
const FluidState &state,
int satPhaseIdx)
{
values[Traits::wPhaseIdx] = 0;
values[Traits::nPhaseIdx] = 0;
if (satPhaseIdx == Traits::wPhaseIdx)
values[Traits::nPhaseIdx] = dPcnw_dSw(params, state);
values[Traits::wettingPhaseIdx] = 0;
values[Traits::nonWettingPhaseIdx] = 0;
if (satPhaseIdx == Traits::wettingPhaseIdx)
values[Traits::nonWettingPhaseIdx] = dPcnw_dSw(params, state);
}
/*!
@@ -195,13 +195,13 @@ public:
const FluidState &state,
int satPhaseIdx)
{
if (satPhaseIdx == Traits::wPhaseIdx) {
values[Traits::wPhaseIdx] = twoPhaseSatDKrw_dSw(params, state.saturation(Traits::wPhaseIdx));
values[Traits::nPhaseIdx] = 0;
if (satPhaseIdx == Traits::wettingPhaseIdx) {
values[Traits::wettingPhaseIdx] = twoPhaseSatDKrw_dSw(params, state.saturation(Traits::wettingPhaseIdx));
values[Traits::nonWettingPhaseIdx] = 0;
}
else {
values[Traits::wPhaseIdx] = 0;
values[Traits::nPhaseIdx] = - twoPhaseSatDKrn_dSw(params, 1 - state.saturation(Traits::nPhaseIdx));
values[Traits::wettingPhaseIdx] = 0;
values[Traits::nonWettingPhaseIdx] = - twoPhaseSatDKrn_dSw(params, 1 - state.saturation(Traits::nonWettingPhaseIdx));
}
}
@@ -266,7 +266,7 @@ public:
template <class FluidState>
static Scalar pcnw(const Params &params, const FluidState &fs)
{
Scalar Sw = fs.saturation(Traits::wPhaseIdx);
Scalar Sw = fs.saturation(Traits::wettingPhaseIdx);
return twoPhaseSatPcnw(params, Sw);
}
@@ -292,7 +292,7 @@ public:
template <class FluidState>
static Scalar Sw(const Params &params, const FluidState &fs)
{
Scalar pc = fs.pressure(Traits::nPhaseIdx) - fs.pressure(Traits::wPhaseIdx);
Scalar pc = fs.pressure(Traits::nonWettingPhaseIdx) - fs.pressure(Traits::wettingPhaseIdx);
return twoPhaseSatSw(params, pc);
}
@@ -330,7 +330,7 @@ public:
template <class FluidState>
static Scalar dPcnw_dSw(const Params &params, const FluidState &fs)
{
Scalar Sw = fs.saturation(Traits::wPhaseIdx);
Scalar Sw = fs.saturation(Traits::wettingPhaseIdx);
return twoPhaseSatDPcnw_dSw(params, Sw);
}
@@ -351,7 +351,7 @@ public:
template <class FluidState>
static Scalar dSw_dpcnw(const Params &params, const FluidState &fs)
{
Scalar Sw = fs.saturation(Traits::wPhaseIdx);
Scalar Sw = fs.saturation(Traits::wettingPhaseIdx);
return twoPhaseSatDSw_dpcnw(params, Sw);
}
@@ -371,7 +371,7 @@ public:
*/
template <class FluidState>
static Scalar krw(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrw(params, fs.saturation(Traits::wPhaseIdx)); }
{ return twoPhaseSatKrw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatKrw(const Params &params, Scalar Sw)
{
@@ -406,7 +406,7 @@ public:
*/
template <class FluidState>
static Scalar krn(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrn(params, 1.0 - fs.saturation(Traits::nPhaseIdx)); }
{ return twoPhaseSatKrn(params, 1.0 - fs.saturation(Traits::nonWettingPhaseIdx)); }
static Scalar twoPhaseSatKrn(const Params &params, Scalar Sw)
{

35
opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp
View File

@@ -1,6 +1,5 @@
/*
Copyright (C) 2008-2013 by Andreas Lauser
Copyright (C) 2011 by Bernd Flemisch
This file is part of the Open Porous Media project (OPM).
@@ -83,10 +82,10 @@ public:
typedef typename Traits::Scalar Scalar;
static const int numPhases = 3;
static const int wPhaseIdx = Traits::wPhaseIdx;
static const int nPhaseIdx = Traits::nPhaseIdx;
static const int oPhaseIdx = Traits::nPhaseIdx;
static const int gPhaseIdx = Traits::gPhaseIdx;
static const int waterPhaseIdx = Traits::wettingPhaseIdx;
static const int nonWettingPhaseIdx = Traits::nonWettingPhaseIdx;
static const int oilPhaseIdx = Traits::nonWettingPhaseIdx;
static const int gasPhaseIdx = Traits::gasPhaseIdx;
//! Specify whether this material law implements the two-phase
//! convenience API
@@ -131,9 +130,9 @@ public:
const Params &params,
const FluidState &state)
{
values[gPhaseIdx] = pcgn(params, state);
values[oPhaseIdx] = 0;
values[wPhaseIdx] = - pcnw(params, state);
values[gasPhaseIdx] = pcgn(params, state);
values[oilPhaseIdx] = 0;
values[waterPhaseIdx] = - pcnw(params, state);
}
/*!
@@ -149,7 +148,7 @@ public:
static Scalar pcgn(const Params &params,
const FluidState &fs)
{
Scalar Sw = 1 - fs.saturation(gPhaseIdx);
Scalar Sw = 1 - fs.saturation(gasPhaseIdx);
return GasOilMaterial::twoPhaseSatPcnw(params.gasOilParams(), Sw);
}
@@ -166,7 +165,7 @@ public:
static Scalar pcnw(const Params &params,
const FluidState &fs)
{
Scalar Sw = fs.saturation(wPhaseIdx);
Scalar Sw = fs.saturation(waterPhaseIdx);
return OilWaterMaterial::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
}
@@ -231,9 +230,9 @@ public:
const Params &params,
const FluidState &fluidState)
{
values[wPhaseIdx] = krw(params, fluidState);
values[nPhaseIdx] = krn(params, fluidState);
values[gPhaseIdx] = krg(params, fluidState);
values[waterPhaseIdx] = krw(params, fluidState);
values[nonWettingPhaseIdx] = krn(params, fluidState);
values[gasPhaseIdx] = krg(params, fluidState);
}
/*!
@@ -243,7 +242,7 @@ public:
static Scalar krg(const Params &params,
const FluidState &fluidState)
{
Scalar Sw = 1 - fluidState.saturation(gPhaseIdx);
Scalar Sw = 1 - fluidState.saturation(gasPhaseIdx);
return GasOilMaterial::twoPhaseSatKrn(params.oilWaterParams(), Sw);
}
@@ -254,7 +253,7 @@ public:
static Scalar krw(const Params &params,
const FluidState &fluidState)
{
Scalar Sw = fluidState.saturation(wPhaseIdx);
Scalar Sw = fluidState.saturation(waterPhaseIdx);
return OilWaterMaterial::twoPhaseSatKrw(params.oilWaterParams(), Sw);
}
@@ -265,9 +264,9 @@ public:
static Scalar krn(const Params &params,
const FluidState &fluidState)
{
Scalar Sw = std::min(1.0, std::max(0.0, fluidState.saturation(wPhaseIdx)));
Scalar So = std::min(1.0, std::max(0.0, fluidState.saturation(oPhaseIdx)));
Scalar Sg = std::min(1.0, std::max(0.0, fluidState.saturation(gPhaseIdx)));
Scalar Sw = std::min(1.0, std::max(0.0, fluidState.saturation(waterPhaseIdx)));
Scalar So = std::min(1.0, std::max(0.0, fluidState.saturation(oilPhaseIdx)));
Scalar Sg = std::min(1.0, std::max(0.0, fluidState.saturation(gasPhaseIdx)));
// connate water. According to the Eclipse TD, this is
// probably only relevant if hysteresis is enabled...

16
opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp
View File

@@ -308,7 +308,7 @@ public:
static typename std::enable_if<implementsTwoPhaseSatApi, ScalarT>::type
twoPhaseSatPcnw(const Params &params, Scalar SwAbs)
{
Scalar SwEff = effectiveSaturation(params, SwAbs, Traits::wPhaseIdx);
Scalar SwEff = effectiveSaturation(params, SwAbs, Traits::wettingPhaseIdx);
return EffLaw::twoPhaseSatPcnw(params, SwEff);
}
@@ -331,12 +331,12 @@ public:
*/
template <class FluidState>
static Scalar Sw(const Params &params, const FluidState &fs)
{ return absoluteSaturation(params, EffLaw::Sw(params, fs), Traits::wPhaseIdx); }
{ return absoluteSaturation(params, EffLaw::Sw(params, fs), Traits::wettingPhaseIdx); }
template <class ScalarT = Scalar>
static typename std::enable_if<implementsTwoPhaseSatApi, ScalarT>::type
twoPhaseSatSw(const Params &params, Scalar Sw)
{ return absoluteSaturation(params, EffLaw::twoPhaseSatSw(params, Sw), Traits::wPhaseIdx); }
{ return absoluteSaturation(params, EffLaw::twoPhaseSatSw(params, Sw), Traits::wettingPhaseIdx); }
/*!
* \brief Calculate non-wetting liquid phase saturation given that
@@ -344,12 +344,12 @@ public:
*/
template <class FluidState>
static Scalar Sn(const Params &params, const FluidState &fs)
{ return absoluteSaturation(params, EffLaw::Sn(params, fs), Traits::nPhaseIdx); }
{ return absoluteSaturation(params, EffLaw::Sn(params, fs), Traits::nonWettingPhaseIdx); }
template <class ScalarT = Scalar>
static typename std::enable_if<implementsTwoPhaseSatApi, ScalarT>::type
twoPhaseSatSn(const Params &params, Scalar Sw)
{ return absoluteSaturation(params, EffLaw::twoPhaseSatSn(params, Sw), Traits::nPhaseIdx); }
{ return absoluteSaturation(params, EffLaw::twoPhaseSatSn(params, Sw), Traits::nonWettingPhaseIdx); }
/*!
* \brief Calculate gas phase saturation given that the rest of
@@ -360,7 +360,7 @@ public:
template <class FluidState, class ScalarT = Scalar>
static typename std::enable_if< (Traits::numPhases > 2), ScalarT>::type
Sg(const Params &params, const FluidState &fs)
{ return absoluteSaturation(params, EffLaw::Sg(params, fs), Traits::gPhaseIdx); }
{ return absoluteSaturation(params, EffLaw::Sg(params, fs), Traits::gasPhaseIdx); }
/*!
* \brief Returns the partial derivative of the capillary
@@ -436,7 +436,7 @@ public:
template <class ScalarT = Scalar>
static typename std::enable_if<implementsTwoPhaseSatApi, ScalarT>::type
twoPhaseSatKrw(const Params &params, Scalar Sw)
{ return EffLaw::twoPhaseSatKrw(params, effectiveSaturation(params, Sw, Traits::nPhaseIdx)); }
{ return EffLaw::twoPhaseSatKrw(params, effectiveSaturation(params, Sw, Traits::nonWettingPhaseIdx)); }
/*!
* \brief The relative permeability of the non-wetting phase.
@@ -464,7 +464,7 @@ public:
template <class ScalarT = Scalar>
static typename std::enable_if<implementsTwoPhaseSatApi, ScalarT>::type
twoPhaseSatKrn(const Params &params, Scalar Sw)
{ return EffLaw::twoPhaseSatKrn(params, effectiveSaturation(params, Sw, Traits::nPhaseIdx)); }
{ return EffLaw::twoPhaseSatKrn(params, effectiveSaturation(params, Sw, Traits::nonWettingPhaseIdx)); }
/*!
* \brief The relative permability of the gas phase

38
opm/material/fluidmatrixinteractions/LinearMaterial.hpp
View File

@@ -267,19 +267,19 @@ public:
template <class FluidState>
static Scalar pcnw(const Params &params, const FluidState &fs)
{
Scalar S = fs.saturation(Traits::wPhaseIdx);
Scalar S = fs.saturation(Traits::wettingPhaseIdx);
Valgrind::CheckDefined(S);
Scalar wPhasePressure =
S*params.pcMaxSat(Traits::wPhaseIdx) +
(1.0 - S)*params.pcMinSat(Traits::wPhaseIdx);
S*params.pcMaxSat(Traits::wettingPhaseIdx) +
(1.0 - S)*params.pcMinSat(Traits::wettingPhaseIdx);
S = fs.saturation(Traits::nPhaseIdx);
S = fs.saturation(Traits::nonWettingPhaseIdx);
Valgrind::CheckDefined(S);
Scalar nPhasePressure =
S*params.pcMaxSat(Traits::nPhaseIdx) +
(1.0 - S)*params.pcMinSat(Traits::nPhaseIdx);
S*params.pcMaxSat(Traits::nonWettingPhaseIdx) +
(1.0 - S)*params.pcMinSat(Traits::nonWettingPhaseIdx);
return nPhasePressure - wPhasePressure;
}
@@ -290,12 +290,12 @@ public:
twoPhaseSatPcnw(const Params &params, Scalar Sw)
{
Scalar wPhasePressure =
Sw*params.pcMaxSat(Traits::wPhaseIdx) +
(1.0 - Sw)*params.pcMinSat(Traits::wPhaseIdx);
Sw*params.pcMaxSat(Traits::wettingPhaseIdx) +
(1.0 - Sw)*params.pcMinSat(Traits::wettingPhaseIdx);
Scalar nPhasePressure =
(1.0 - Sw)*params.pcMaxSat(Traits::nPhaseIdx) +
Sw*params.pcMinSat(Traits::nPhaseIdx);
(1.0 - Sw)*params.pcMaxSat(Traits::nonWettingPhaseIdx) +
Sw*params.pcMinSat(Traits::nonWettingPhaseIdx);
return nPhasePressure - wPhasePressure;
}
@@ -342,7 +342,7 @@ public:
*/
template <class FluidState>
static Scalar krw(const Params &params, const FluidState &fs)
{ return std::max(0.0, std::min(1.0, fs.saturation(Traits::wPhaseIdx))); }
{ return std::max(0.0, std::min(1.0, fs.saturation(Traits::wettingPhaseIdx))); }
template <class ScalarT = Scalar>
static typename std::enable_if<Traits::numPhases == 2, ScalarT>::type
@@ -354,7 +354,7 @@ public:
*/
template <class FluidState>
static Scalar krn(const Params &params, const FluidState &fs)
{ return std::max(0.0, std::min(1.0, fs.saturation(Traits::nPhaseIdx))); }
{ return std::max(0.0, std::min(1.0, fs.saturation(Traits::nonWettingPhaseIdx))); }
template <class ScalarT = Scalar>
static typename std::enable_if<Traits::numPhases == 2, ScalarT>::type
@@ -369,7 +369,7 @@ public:
template <class FluidState, class ScalarT=Scalar>
static typename std::enable_if< (Traits::numPhases > 2), ScalarT>::type
krg(const Params &params, const FluidState &fs)
{ return std::max(0.0, std::min(1.0, fs.saturation(Traits::gPhaseIdx))); }
{ return std::max(0.0, std::min(1.0, fs.saturation(Traits::gasPhaseIdx))); }
/*!
* \brief The difference between the pressures of the gas and the non-wetting phase.
@@ -380,19 +380,19 @@ public:
static typename std::enable_if< (Traits::numPhases > 2), ScalarT>::type
pcgn(const Params &params, const FluidState &fs)
{
Scalar S = fs.saturation(Traits::nPhaseIdx);
Scalar S = fs.saturation(Traits::nonWettingPhaseIdx);
Valgrind::CheckDefined(S);
Scalar nPhasePressure =
S*params.pcMaxSat(Traits::nPhaseIdx) +
(1.0 - S)*params.pcMinSat(Traits::nPhaseIdx);
S*params.pcMaxSat(Traits::nonWettingPhaseIdx) +
(1.0 - S)*params.pcMinSat(Traits::nonWettingPhaseIdx);
S = fs.saturation(Traits::gPhaseIdx);
S = fs.saturation(Traits::gasPhaseIdx);
Valgrind::CheckDefined(S);
Scalar gPhasePressure =
S*params.pcMaxSat(Traits::gPhaseIdx) +
(1.0 - S)*params.pcMinSat(Traits::gPhaseIdx);
S*params.pcMaxSat(Traits::gasPhaseIdx) +
(1.0 - S)*params.pcMinSat(Traits::gasPhaseIdx);
return gPhasePressure - nPhasePressure;
}

50
opm/material/fluidmatrixinteractions/MaterialTraits.hpp
View File

@@ -51,7 +51,7 @@ public:
*
* \brief A generic traits class for two-phase material laws.
*/
template <class ScalarT, int wettingPhaseIdxV, int nonWettingPhaseIdxV>
template <class ScalarT, int wettinPhaseIdxV, int nonWettingasPhaseIdxV>
class TwoPhaseMaterialTraits
{
public:
@@ -62,19 +62,19 @@ public:
static const int numPhases = 2;
//! The index of the wetting phase
static const int wPhaseIdx = wettingPhaseIdxV;
static const int wettingPhaseIdx = wettinPhaseIdxV;
//! The index of the non-wetting phase
static const int nPhaseIdx = nonWettingPhaseIdxV;
static const int nonWettingPhaseIdx = nonWettingasPhaseIdxV;
// some safety checks...
static_assert(0 <= wPhaseIdx && wPhaseIdx < numPhases,
"wPhaseIdx is out of range");
static_assert(0 <= nPhaseIdx && nPhaseIdx < numPhases,
"nPhaseIdx is out of range");
static_assert(0 <= wettingPhaseIdx && wettingPhaseIdx < numPhases,
"wettingPhaseIdx is out of range");
static_assert(0 <= nonWettingPhaseIdx && nonWettingPhaseIdx < numPhases,
"nonWettingPhaseIdx is out of range");
static_assert(wPhaseIdx != nPhaseIdx,
"wPhaseIdx and nPhaseIdx must be different");
static_assert(wettingPhaseIdx != nonWettingPhaseIdx,
"wettingPhaseIdx and nonWettingPhaseIdx must be different");
};
/*!
@@ -82,7 +82,7 @@ public:
*
* \brief A generic traits class for three-phase material laws.
*/
template <class ScalarT, int wettingPhaseIdxV, int nonWettingPhaseIdxV, int gasPhaseIdxV>
template <class ScalarT, int wettingPhaseIdxV, int nonWettingasPhaseIdxV, int gasPhaseIdxV>
class ThreePhaseMaterialTraits
{
public:
@@ -93,28 +93,28 @@ public:
static const int numPhases = 3;
//! The index of the wetting liquid phase
static const int wPhaseIdx = wettingPhaseIdxV;
static const int wettingPhaseIdx = wettingPhaseIdxV;
//! The index of the non-wetting liquid phase
static const int nPhaseIdx = nonWettingPhaseIdxV;
static const int nonWettingPhaseIdx = nonWettingasPhaseIdxV;
//! The index of the gas phase (i.e., the least wetting phase)
static const int gPhaseIdx = gasPhaseIdxV;
static const int gasPhaseIdx = gasPhaseIdxV;
// some safety checks...
static_assert(0 <= wPhaseIdx && wPhaseIdx < numPhases,
"wPhaseIdx is out of range");
static_assert(0 <= nPhaseIdx && nPhaseIdx < numPhases,
"nPhaseIdx is out of range");
static_assert(0 <= gPhaseIdx && gPhaseIdx < numPhases,
"gPhaseIdx is out of range");
static_assert(0 <= wettingPhaseIdx && wettingPhaseIdx < numPhases,
"wettingPhaseIdx is out of range");
static_assert(0 <= nonWettingPhaseIdx && nonWettingPhaseIdx < numPhases,
"nonWettingPhaseIdx is out of range");
static_assert(0 <= gasPhaseIdx && gasPhaseIdx < numPhases,
"gasPhaseIdx is out of range");
static_assert(wPhaseIdx != nPhaseIdx,
"wPhaseIdx and nPhaseIdx must be different");
static_assert(wPhaseIdx != gPhaseIdx,
"wPhaseIdx and gPhaseIdx must be different");
static_assert(nPhaseIdx != gPhaseIdx,
"nPhaseIdx and gPhaseIdx must be different");
static_assert(wettingPhaseIdx != nonWettingPhaseIdx,
"wettingPhaseIdx and nonWettingPhaseIdx must be different");
static_assert(wettingPhaseIdx != gasPhaseIdx,
"wettingPhaseIdx and gasPhaseIdx must be different");
static_assert(nonWettingPhaseIdx != gasPhaseIdx,
"nonWettingPhaseIdx and gasPhaseIdx must be different");
};
} // namespace Opm

6
opm/material/fluidmatrixinteractions/NullMaterial.hpp
View File

@@ -296,7 +296,7 @@ public:
template <class FluidState, class ScalarT = Scalar>
static typename std::enable_if<(numPhases > 1), ScalarT>::type
krw(const Params &params, const FluidState &fs)
{ return std::max(0.0, std::min(1.0, fs.saturation(Traits::wPhaseIdx))); }
{ return std::max(0.0, std::min(1.0, fs.saturation(Traits::wettingPhaseIdx))); }
template <class ScalarT = Scalar>
static typename std::enable_if<numPhases == 2, ScalarT>::type
@@ -309,7 +309,7 @@ public:
template <class FluidState, class ScalarT=Scalar>
static typename std::enable_if<(numPhases > 1), ScalarT>::type
krn(const Params &params, const FluidState &fs)
{ return std::max(0.0, std::min(1.0, fs.saturation(Traits::nPhaseIdx))); }
{ return std::max(0.0, std::min(1.0, fs.saturation(Traits::nonWettingPhaseIdx))); }
template <class ScalarT = Scalar>
static typename std::enable_if<numPhases == 2, ScalarT>::type
@@ -324,7 +324,7 @@ public:
template <class FluidState, class ScalarT=Scalar>
static typename std::enable_if< (numPhases > 2), ScalarT>::type
krg(const Params &params, const FluidState &fs)
{ return std::max(0.0, std::min(1.0, fs.saturation(Traits::gPhaseIdx))); }
{ return std::max(0.0, std::min(1.0, fs.saturation(Traits::gasPhaseIdx))); }
/*!
* \brief The difference between the pressures of the gas and the non-wetting phase.

34
opm/material/fluidmatrixinteractions/ParkerLenhard.hpp
View File

@@ -299,7 +299,7 @@ public:
template <class FluidState>
static void update(Params &params, const FluidState &fs)
{
Scalar Sw = fs.saturation(Traits::wPhaseIdx);
Scalar Sw = fs.saturation(Traits::wettingPhaseIdx);
if (Sw > 1 - 1e-5) {
// if the absolute saturation is almost 1,
@@ -339,8 +339,8 @@ public:
template <class Container, class FluidState>
static void capillaryPressures(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wPhaseIdx] = 0.0; // reference phase
values[Traits::nPhaseIdx] = pcnw(params, fs);
values[Traits::wettingPhaseIdx] = 0.0; // reference phase
values[Traits::nonWettingPhaseIdx] = pcnw(params, fs);
}
/*!
@@ -358,8 +358,8 @@ public:
template <class Container, class FluidState>
static void relativePermeabilities(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wPhaseIdx] = krw(params, fs);
values[Traits::nPhaseIdx] = krn(params, fs);
values[Traits::wettingPhaseIdx] = krw(params, fs);
values[Traits::nonWettingPhaseIdx] = krn(params, fs);
}
/*!
@@ -372,10 +372,10 @@ public:
const FluidState &state,
int satPhaseIdx)
{
values[Traits::wPhaseIdx] = 0;
values[Traits::nPhaseIdx] = 0;
if (satPhaseIdx == Traits::wPhaseIdx)
values[Traits::nPhaseIdx] = twoPhaseSatDPcnw_dSw(params, state.saturation(Traits::wPhaseIdx));
values[Traits::wettingPhaseIdx] = 0;
values[Traits::nonWettingPhaseIdx] = 0;
if (satPhaseIdx == Traits::wettingPhaseIdx)
values[Traits::nonWettingPhaseIdx] = twoPhaseSatDPcnw_dSw(params, state.saturation(Traits::wettingPhaseIdx));
}
/*!
@@ -433,13 +433,13 @@ public:
const FluidState &state,
int satPhaseIdx)
{
if (satPhaseIdx == Traits::wPhaseIdx) {
values[Traits::wPhaseIdx] = twoPhaseSatDKrw_dSw(params, state.saturation(Traits::wPhaseIdx));
values[Traits::nPhaseIdx] = 0;
if (satPhaseIdx == Traits::wettingPhaseIdx) {
values[Traits::wettingPhaseIdx] = twoPhaseSatDKrw_dSw(params, state.saturation(Traits::wettingPhaseIdx));
values[Traits::nonWettingPhaseIdx] = 0;
}
else {
values[Traits::wPhaseIdx] = 0;
values[Traits::nPhaseIdx] = - twoPhaseSatDKrn_dSw(params, 1 - state.saturation(Traits::nPhaseIdx));
values[Traits::wettingPhaseIdx] = 0;
values[Traits::nonWettingPhaseIdx] = - twoPhaseSatDKrn_dSw(params, 1 - state.saturation(Traits::nonWettingPhaseIdx));
}
}
@@ -494,7 +494,7 @@ public:
*/
template <class FluidState>
static Scalar pcnw(const Params &params, const FluidState &fs)
{ return twoPhaseSatPcnw(params, fs.saturation(Traits::wPhaseIdx)); }
{ return twoPhaseSatPcnw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatPcnw(const Params &params, Scalar Sw)
{
@@ -560,7 +560,7 @@ public:
*/
template <class FluidState>
static Scalar krw(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrw(params, fs.saturation(Traits::wPhaseIdx)); }
{ return twoPhaseSatKrw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatKrw(const Params &params, Scalar Sw)
{
@@ -579,7 +579,7 @@ public:
*/
template <class FluidState>
static Scalar krn(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrn(params, fs.saturation(Traits::wPhaseIdx)); }
{ return twoPhaseSatKrn(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatKrn(const Params &params, Scalar Sw)
{

38
opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterial.hpp
View File

@@ -94,8 +94,8 @@ public:
template <class Container, class FluidState>
static void capillaryPressures(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wPhaseIdx] = 0.0; // reference phase
values[Traits::nPhaseIdx] = pcnw(params, fs);
values[Traits::wettingPhaseIdx] = 0.0; // reference phase
values[Traits::nonWettingPhaseIdx] = pcnw(params, fs);
}
/*!
@@ -112,8 +112,8 @@ public:
template <class Container, class FluidState>
static void relativePermeabilities(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wPhaseIdx] = krw(params, fs);
values[Traits::nPhaseIdx] = krn(params, fs);
values[Traits::wettingPhaseIdx] = krw(params, fs);
values[Traits::nonWettingPhaseIdx] = krn(params, fs);
}
@@ -127,10 +127,10 @@ public:
const FluidState &state,
int satPhaseIdx)
{
values[Traits::wPhaseIdx] = 0;
values[Traits::nPhaseIdx] = 0;
if (satPhaseIdx == Traits::wPhaseIdx) {
values[Traits::nPhaseIdx] = evalDeriv_(params.pcnwSamples(), state.saturation(Traits::wPhaseIdx));
values[Traits::wettingPhaseIdx] = 0;
values[Traits::nonWettingPhaseIdx] = 0;
if (satPhaseIdx == Traits::wettingPhaseIdx) {
values[Traits::nonWettingPhaseIdx] = evalDeriv_(params.pcnwSamples(), state.saturation(Traits::wettingPhaseIdx));
}
}
@@ -189,13 +189,13 @@ public:
const FluidState &state,
int satPhaseIdx)
{
if (satPhaseIdx == Traits::wPhaseIdx) {
values[Traits::wPhaseIdx] = evalDeriv_(params.krwSamples(), state.saturation(Traits::wPhaseIdx));
values[Traits::nPhaseIdx] = 0;
if (satPhaseIdx == Traits::wettingPhaseIdx) {
values[Traits::wettingPhaseIdx] = evalDeriv_(params.krwSamples(), state.saturation(Traits::wettingPhaseIdx));
values[Traits::nonWettingPhaseIdx] = 0;
}
else {
values[Traits::wPhaseIdx] = 0;
values[Traits::nPhaseIdx] = - evalDeriv_(params.krwSamples(), 1 - state.saturation(Traits::nPhaseIdx));
values[Traits::wettingPhaseIdx] = 0;
values[Traits::nonWettingPhaseIdx] = - evalDeriv_(params.krwSamples(), 1 - state.saturation(Traits::nonWettingPhaseIdx));
}
}
@@ -250,7 +250,7 @@ public:
template <class FluidState>
static Scalar pcnw(const Params &params, const FluidState &fs)
{
Scalar Sw = fs.saturation(Traits::wPhaseIdx);
Scalar Sw = fs.saturation(Traits::wettingPhaseIdx);
return twoPhaseSatPcnw(params, Sw);
}
@@ -288,7 +288,7 @@ public:
*/
template <class FluidState>
static Scalar dPcnw_dSw(const Params &params, const FluidState &fs)
{ return twoPhaseSatDPcnw_dSw(params, fs.saturation(Traits::wPhaseIdx)); }
{ return twoPhaseSatDPcnw_dSw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatDPcnw_dSw(const Params &params, Scalar Sw)
{
@@ -302,7 +302,7 @@ public:
*/
template <class FluidState>
static Scalar krw(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrw(params, fs.saturation(Traits::wPhaseIdx)); }
{ return twoPhaseSatKrw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatKrw(const Params &params, Scalar Sw)
{
@@ -321,7 +321,7 @@ public:
*/
template <class FluidState>
static Scalar dKrw_dSw(const Params &params, const FluidState &fs)
{ return twoPhaseSatDkrw_dSw(params, fs.saturation(Traits::wPhaseIdx)); }
{ return twoPhaseSatDkrw_dSw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatDKrw_dSw(const Params &params, Scalar Sw)
{
@@ -339,7 +339,7 @@ public:
*/
template <class FluidState>
static Scalar krn(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrn(params, 1.0 - fs.saturation(Traits::nPhaseIdx)); }
{ return twoPhaseSatKrn(params, 1.0 - fs.saturation(Traits::nonWettingPhaseIdx)); }
static Scalar twoPhaseSatKrn(const Params &params, Scalar Sw)
{
@@ -358,7 +358,7 @@ public:
*/
template <class FluidState>
static Scalar dKrn_dSw(const Params &params, const FluidState &fs)
{ return twoPhaseSatDkrn_dSw(params, fs.saturation(Traits::wPhaseIdx)); }
{ return twoPhaseSatDkrn_dSw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatDKrn_dSw(const Params &params, Scalar Sw)
{

38
opm/material/fluidmatrixinteractions/RegularizedBrooksCorey.hpp
View File

@@ -115,8 +115,8 @@ public:
template <class Container, class FluidState>
static void capillaryPressures(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wPhaseIdx] = 0.0; // reference phase
values[Traits::nPhaseIdx] = pcnw(params, fs);
values[Traits::wettingPhaseIdx] = 0.0; // reference phase
values[Traits::nonWettingPhaseIdx] = pcnw(params, fs);
}
/*!
@@ -126,8 +126,8 @@ public:
template <class Container, class FluidState>
static void saturations(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wPhaseIdx] = Sw(params, fs);
values[Traits::nPhaseIdx] = 1 - values[Traits::wPhaseIdx];
values[Traits::wettingPhaseIdx] = Sw(params, fs);
values[Traits::nonWettingPhaseIdx] = 1 - values[Traits::wettingPhaseIdx];
}
/*!
@@ -143,8 +143,8 @@ public:
template <class Container, class FluidState>
static void relativePermeabilities(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wPhaseIdx] = krw(params, fs);
values[Traits::nPhaseIdx] = krn(params, fs);
values[Traits::wettingPhaseIdx] = krw(params, fs);
values[Traits::nonWettingPhaseIdx] = krn(params, fs);
}
/*!
@@ -157,10 +157,10 @@ public:
const FluidState &state,
int satPhaseIdx)
{
values[Traits::wPhaseIdx] = 0;
values[Traits::nPhaseIdx] = 0;
if (satPhaseIdx == Traits::wPhaseIdx)
values[Traits::nPhaseIdx] = twoPhaseSatDPcnw_dSw(params, state.saturation(Traits::wPhaseIdx));
values[Traits::wettingPhaseIdx] = 0;
values[Traits::nonWettingPhaseIdx] = 0;
if (satPhaseIdx == Traits::wettingPhaseIdx)
values[Traits::nonWettingPhaseIdx] = twoPhaseSatDPcnw_dSw(params, state.saturation(Traits::wettingPhaseIdx));
}
/*!
@@ -218,13 +218,13 @@ public:
const FluidState &state,
int satPhaseIdx)
{
if (satPhaseIdx == Traits::wPhaseIdx) {
values[Traits::wPhaseIdx] = twoPhaseSatDKrw_dSw(params, state.saturation(Traits::wPhaseIdx));
values[Traits::nPhaseIdx] = 0;
if (satPhaseIdx == Traits::wettingPhaseIdx) {
values[Traits::wettingPhaseIdx] = twoPhaseSatDKrw_dSw(params, state.saturation(Traits::wettingPhaseIdx));
values[Traits::nonWettingPhaseIdx] = 0;
}
else {
values[Traits::wPhaseIdx] = 0;
values[Traits::nPhaseIdx] = - twoPhaseSatDKrn_dSw(params, 1 - state.saturation(Traits::nPhaseIdx));
values[Traits::wettingPhaseIdx] = 0;
values[Traits::nonWettingPhaseIdx] = - twoPhaseSatDKrn_dSw(params, 1 - state.saturation(Traits::nonWettingPhaseIdx));
}
}
@@ -299,7 +299,7 @@ public:
*/
template <class FluidState>
static Scalar pcnw(const Params &params, const FluidState &fs)
{ return twoPhaseSatPcnw(params, fs.saturation(Traits::wPhaseIdx)); }
{ return twoPhaseSatPcnw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatPcnw(const Params &params, Scalar Sw)
{
@@ -336,7 +336,7 @@ public:
template <class FluidState>
static Scalar Sw(const Params &params, const FluidState &fs)
{
Scalar pcnw = fs.pressure(Traits::nPhaseIdx) - fs.pressure(Traits::wPhaseIdx);
Scalar pcnw = fs.pressure(Traits::nonWettingPhaseIdx) - fs.pressure(Traits::wettingPhaseIdx);
return twoPhaseSatSw(params, pcnw);
}
@@ -455,7 +455,7 @@ public:
*/
template <class FluidState>
static Scalar krw(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrw(params, fs.saturation(Traits::wPhaseIdx)); }
{ return twoPhaseSatKrw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatKrw(const Params &params, Scalar Sw)
{
@@ -491,7 +491,7 @@ public:
*/
template <class FluidState>
static Scalar krn(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrn(params, 1.0 - fs.saturation(Traits::nPhaseIdx)); }
{ return twoPhaseSatKrn(params, 1.0 - fs.saturation(Traits::nonWettingPhaseIdx)); }
static Scalar twoPhaseSatKrn(const Params &params, Scalar Sw)
{

42
opm/material/fluidmatrixinteractions/RegularizedVanGenuchten.hpp
View File

@@ -113,8 +113,8 @@ public:
template <class Container, class FluidState>
static void capillaryPressures(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wPhaseIdx] = 0.0; // reference phase
values[Traits::nPhaseIdx] = pcnw(params, fs);
values[Traits::wettingPhaseIdx] = 0.0; // reference phase
values[Traits::nonWettingPhaseIdx] = pcnw(params, fs);
}
/*!
@@ -124,8 +124,8 @@ public:
template <class Container, class FluidState>
static void saturations(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wPhaseIdx] = Sw(params, fs);
values[Traits::nPhaseIdx] = 1 - values[Traits::wPhaseIdx];
values[Traits::wettingPhaseIdx] = Sw(params, fs);
values[Traits::nonWettingPhaseIdx] = 1 - values[Traits::wettingPhaseIdx];
}
/*!
@@ -135,8 +135,8 @@ public:
template <class Container, class FluidState>
static void relativePermeabilities(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wPhaseIdx] = krw(params, fs);
values[Traits::nPhaseIdx] = krn(params, fs);
values[Traits::wettingPhaseIdx] = krw(params, fs);
values[Traits::nonWettingPhaseIdx] = krn(params, fs);
}
@@ -150,10 +150,10 @@ public:
const FluidState &state,
int satPhaseIdx)
{
values[Traits::wPhaseIdx] = 0;
values[Traits::nPhaseIdx] = 0;
if (satPhaseIdx == Traits::wPhaseIdx)
values[Traits::nPhaseIdx] = twoPhaseSatDPcnw_dSw(params, state.saturation(Traits::wPhaseIdx));
values[Traits::wettingPhaseIdx] = 0;
values[Traits::nonWettingPhaseIdx] = 0;
if (satPhaseIdx == Traits::wettingPhaseIdx)
values[Traits::nonWettingPhaseIdx] = twoPhaseSatDPcnw_dSw(params, state.saturation(Traits::wettingPhaseIdx));
}
/*!
@@ -211,13 +211,13 @@ public:
const FluidState &state,
int satPhaseIdx)
{
if (satPhaseIdx == Traits::wPhaseIdx) {
values[Traits::wPhaseIdx] = twoPhaseSatDKrw_dSw(params, state.saturation(Traits::wPhaseIdx));
values[Traits::nPhaseIdx] = 0;
if (satPhaseIdx == Traits::wettingPhaseIdx) {
values[Traits::wettingPhaseIdx] = twoPhaseSatDKrw_dSw(params, state.saturation(Traits::wettingPhaseIdx));
values[Traits::nonWettingPhaseIdx] = 0;
}
else {
values[Traits::wPhaseIdx] = 0;
values[Traits::nPhaseIdx] = - twoPhaseSatDKrn_dSw(params, 1 - state.saturation(Traits::nPhaseIdx));
values[Traits::wettingPhaseIdx] = 0;
values[Traits::nonWettingPhaseIdx] = - twoPhaseSatDKrn_dSw(params, 1 - state.saturation(Traits::nonWettingPhaseIdx));
}
}
@@ -280,7 +280,7 @@ public:
*/
template <class FluidState>
static Scalar pcnw(const Params &params, const FluidState &fs)
{ return twoPhaseSatPcnw(params, fs.saturation(Traits::wPhaseIdx)); }
{ return twoPhaseSatPcnw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatPcnw(const Params &params, Scalar Sw)
{
@@ -336,7 +336,7 @@ public:
template <class FluidState>
static Scalar Sw(const Params &params, const FluidState &fs)
{
Scalar pC = fs.pressure(Traits::nPhaseIdx) - fs.pressure(Traits::wPhaseIdx);
Scalar pC = fs.pressure(Traits::nonWettingPhaseIdx) - fs.pressure(Traits::wettingPhaseIdx);
return twoPhaseSatSw(params, pC);
}
@@ -404,7 +404,7 @@ public:
*/
template <class FluidState>
static Scalar dPcnw_dSw(const Params &params, const FluidState &fs)
{ return twoPhaseSatDPcnw_dSw(params, fs.saturation(Traits::wPhaseIdx)); }
{ return twoPhaseSatDPcnw_dSw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatDPcnw_dSw(const Params &params, Scalar Sw)
{
@@ -440,7 +440,7 @@ public:
template <class FluidState>
static Scalar dSw_dpC(const Params &params, const FluidState &fs)
{
Scalar pC = fs.pressure(Traits::nPhaseIdx) - fs.pressure(Traits::wPhaseIdx);
Scalar pC = fs.pressure(Traits::nonWettingPhaseIdx) - fs.pressure(Traits::wettingPhaseIdx);
// calculate the saturation which corrosponds to the
// saturation in the non-regularized verision of van
@@ -483,7 +483,7 @@ public:
*/
template <class FluidState>
static Scalar krw(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrw(params, fs.saturation(Traits::wPhaseIdx)); }
{ return twoPhaseSatKrw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatKrw(const Params &params, Scalar Sw)
{
@@ -520,7 +520,7 @@ public:
*/
template <class FluidState>
static Scalar krn(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrn(params, 1.0 - fs.saturation(Traits::nPhaseIdx)); }
{ return twoPhaseSatKrn(params, 1.0 - fs.saturation(Traits::nonWettingPhaseIdx)); }
static Scalar twoPhaseSatKrn(const Params &params, Scalar Sw)
{

32
opm/material/fluidmatrixinteractions/ThreePhaseParkerVanGenuchten.hpp
View File

@@ -55,9 +55,9 @@ public:
typedef typename Traits::Scalar Scalar;
static const int numPhases = 3;
static const int wPhaseIdx = Traits::wPhaseIdx;
static const int nPhaseIdx = Traits::nPhaseIdx;
static const int gPhaseIdx = Traits::gPhaseIdx;
static const int wettingPhaseIdx = Traits::wettingPhaseIdx;
static const int nonWettingPhaseIdx = Traits::nonWettingPhaseIdx;
static const int gasPhaseIdx = Traits::gasPhaseIdx;
//! Specify whether this material law implements the two-phase
//! convenience API
@@ -99,9 +99,9 @@ public:
const Params &params,
const FluidState &fluidState)
{
values[gPhaseIdx] = pcgn(params, fluidState);
values[nPhaseIdx] = 0;
values[wPhaseIdx] = - pcnw(params, fluidState);
values[gasPhaseIdx] = pcgn(params, fluidState);
values[nonWettingPhaseIdx] = 0;
values[wettingPhaseIdx] = - pcnw(params, fluidState);
}
/*!
@@ -121,8 +121,8 @@ public:
// sum of liquid saturations
Scalar St =
fluidState.saturation(wPhaseIdx)
+ fluidState.saturation(nPhaseIdx);
fluidState.saturation(wettingPhaseIdx)
+ fluidState.saturation(nonWettingPhaseIdx);
Scalar Se = (St - params.Swrx())/(1. - params.Swrx());
@@ -171,7 +171,7 @@ public:
static Scalar pcnw(const Params &params,
const FluidState &fluidState)
{
Scalar Sw = fluidState.saturation(wPhaseIdx);
Scalar Sw = fluidState.saturation(wettingPhaseIdx);
Scalar Se = (Sw-params.Swr())/(1.-params.Snr());
Scalar PC_VG_REG = 0.01;
@@ -251,9 +251,9 @@ public:
const Params &params,
const FluidState &fluidState)
{
values[wPhaseIdx] = krw(params, fluidState);
values[nPhaseIdx] = krn(params, fluidState);
values[gPhaseIdx] = krg(params, fluidState);
values[wettingPhaseIdx] = krw(params, fluidState);
values[nonWettingPhaseIdx] = krn(params, fluidState);
values[gasPhaseIdx] = krg(params, fluidState);
}
/*!
@@ -270,7 +270,7 @@ public:
{
// transformation to effective saturation
Scalar Se = (fluidState.saturation(wPhaseIdx) - params.Swr()) / (1-params.Swr());
Scalar Se = (fluidState.saturation(wettingPhaseIdx) - params.Swr()) / (1-params.Swr());
// regularization
if(Se > 1.0) return 1.;
@@ -296,8 +296,8 @@ public:
static Scalar krn(const Params &params,
const FluidState &fluidState)
{
Scalar Sn = fluidState.saturation(nPhaseIdx);
Scalar Sw = fluidState.saturation(wPhaseIdx);
Scalar Sn = fluidState.saturation(nonWettingPhaseIdx);
Scalar Sw = fluidState.saturation(wettingPhaseIdx);
Scalar Swe = std::min((Sw - params.Swr()) / (1 - params.Swr()), 1.);
Scalar Ste = std::min((Sw + Sn - params.Swr()) / (1 - params.Swr()), 1.);
@@ -341,7 +341,7 @@ public:
static Scalar krg(const Params &params,
const FluidState &fluidState)
{
Scalar Sg = fluidState.saturation(gPhaseIdx);
Scalar Sg = fluidState.saturation(gasPhaseIdx);
Scalar Se = std::min(((1-Sg) - params.Sgr()) / (1 - params.Sgr()), 1.);
// regularization

44
opm/material/fluidmatrixinteractions/VanGenuchten.hpp
View File

@@ -108,8 +108,8 @@ public:
template <class Container, class FluidState>
static void capillaryPressures(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wPhaseIdx] = 0.0; // reference phase
values[Traits::nPhaseIdx] = pcnw(params, fs);
values[Traits::wettingPhaseIdx] = 0.0; // reference phase
values[Traits::nonWettingPhaseIdx] = pcnw(params, fs);
}
/*!
@@ -119,8 +119,8 @@ public:
template <class Container, class FluidState>
static void saturations(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wPhaseIdx] = Sw(params, fs);
values[Traits::nPhaseIdx] = 1 - values[Traits::wPhaseIdx];
values[Traits::wettingPhaseIdx] = Sw(params, fs);
values[Traits::nonWettingPhaseIdx] = 1 - values[Traits::wettingPhaseIdx];
}
/*!
@@ -136,8 +136,8 @@ public:
template <class Container, class FluidState>
static void relativePermeabilities(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wPhaseIdx] = krw(params, fs);
values[Traits::nPhaseIdx] = krn(params, fs);
values[Traits::wettingPhaseIdx] = krw(params, fs);
values[Traits::nonWettingPhaseIdx] = krn(params, fs);
}
@@ -151,10 +151,10 @@ public:
const FluidState &state,
int satPhaseIdx)
{
values[Traits::wPhaseIdx] = 0;
values[Traits::nPhaseIdx] = 0;
if (satPhaseIdx == Traits::wPhaseIdx)
values[Traits::nPhaseIdx] = dPcnw_dSw(params, state);
values[Traits::wettingPhaseIdx] = 0;
values[Traits::nonWettingPhaseIdx] = 0;
if (satPhaseIdx == Traits::wettingPhaseIdx)
values[Traits::nonWettingPhaseIdx] = dPcnw_dSw(params, state);
}
/*!
@@ -212,13 +212,13 @@ public:
const FluidState &state,
int satPhaseIdx)
{
if (satPhaseIdx == Traits::wPhaseIdx) {
values[Traits::wPhaseIdx] = twoPhaseSatDKrw_dSw(params, state.saturation(Traits::wPhaseIdx));
values[Traits::nPhaseIdx] = 0;
if (satPhaseIdx == Traits::wettingPhaseIdx) {
values[Traits::wettingPhaseIdx] = twoPhaseSatDKrw_dSw(params, state.saturation(Traits::wettingPhaseIdx));
values[Traits::nonWettingPhaseIdx] = 0;
}
else {
values[Traits::wPhaseIdx] = 0;
values[Traits::nPhaseIdx] = - twoPhaseSatDKrn_dSw(params, 1 - state.saturation(Traits::nPhaseIdx));
values[Traits::wettingPhaseIdx] = 0;
values[Traits::nonWettingPhaseIdx] = - twoPhaseSatDKrn_dSw(params, 1 - state.saturation(Traits::nonWettingPhaseIdx));
}
}
@@ -284,7 +284,7 @@ public:
template <class FluidState>
static Scalar pcnw(const Params &params, const FluidState &fs)
{
Scalar Sw = fs.saturation(Traits::wPhaseIdx);
Scalar Sw = fs.saturation(Traits::wettingPhaseIdx);
assert(0 <= Sw && Sw <= 1);
return twoPhaseSatPcnw(params, Sw);
@@ -322,7 +322,7 @@ public:
template <class FluidState>
static Scalar Sw(const Params &params, const FluidState &fs)
{
Scalar pC = fs.pressure(Traits::nPhaseIdx) - fs.pressure(Traits::wPhaseIdx);
Scalar pC = fs.pressure(Traits::nonWettingPhaseIdx) - fs.pressure(Traits::wettingPhaseIdx);
return twoPhaseSatSw(params, pC);
}
@@ -360,7 +360,7 @@ public:
*/
template <class FluidState>
static Scalar dPcnw_dSw(const Params &params, const FluidState &fs)
{ return twoPhaseSatDPcnw_dSw(params, fs.saturation(Traits::wPhaseIdx)); }
{ return twoPhaseSatDPcnw_dSw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatDPcnw_dSw(const Params &params, Scalar Sw)
{
@@ -384,7 +384,7 @@ public:
*/
template <class FluidState>
static Scalar krw(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrw(params, fs.saturation(Traits::wPhaseIdx)); }
{ return twoPhaseSatKrw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatKrw(const Params &params, Scalar Sw)
{
@@ -407,7 +407,7 @@ public:
*/
template <class FluidState>
static Scalar dKrw_dSw(const Params &params, const FluidState &fs)
{ return twoPhaseSatDkrw_dSw(params, fs.saturation(Traits::wPhaseIdx)); }
{ return twoPhaseSatDkrw_dSw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatDKrw_dSw(const Params &params, Scalar Sw)
{
@@ -430,7 +430,7 @@ public:
*/
template <class FluidState>
static Scalar krn(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrn(params, 1.0 - fs.saturation(Traits::nPhaseIdx)); }
{ return twoPhaseSatKrn(params, 1.0 - fs.saturation(Traits::nonWettingPhaseIdx)); }
static Scalar twoPhaseSatKrn(const Params &params, Scalar Sw)
{
@@ -454,7 +454,7 @@ public:
*/
template <class FluidState>
static Scalar dKrn_dSw(const Params &params, const FluidState &fs)
{ return twoPhaseSatDkrn_dSw(params, fs.saturation(Traits::wPhaseIdx)); }
{ return twoPhaseSatDkrn_dSw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatDKrn_dSw(const Params &params, Scalar Sw)
{

34
opm/material/fluidsystems/2pImmiscibleFluidSystem.hpp
View File

@@ -72,9 +72,9 @@ public:
static const int numPhases = 2;
//! Index of the wetting phase
static const int wPhaseIdx = 0;
static const int wettingPhaseIdx = 0;
//! Index of the non-wetting phase
static const int nPhaseIdx = 1;
static const int nonWettingPhaseIdx = 1;
//! \copydoc BaseFluidSystem::phaseName
static const char *phaseName(int phaseIdx)
@@ -93,7 +93,7 @@ public:
{
//assert(0 <= phaseIdx && phaseIdx < numPhases);
return
(phaseIdx == wPhaseIdx)
(phaseIdx == wettingPhaseIdx)
? WettingPhase::isLiquid()
: NonwettingPhase::isLiquid();
}
@@ -104,7 +104,7 @@ public:
//assert(0 <= phaseIdx && phaseIdx < numPhases);
return
(phaseIdx == wPhaseIdx)
(phaseIdx == wettingPhaseIdx)
? WettingPhase::isCompressible()
: NonwettingPhase::isCompressible();
}
@@ -116,7 +116,7 @@ public:
// let the fluids decide
return
(phaseIdx == wPhaseIdx)
(phaseIdx == wettingPhaseIdx)
? WettingPhase::isIdealGas()
: NonwettingPhase::isIdealGas();
}
@@ -138,16 +138,16 @@ public:
static const int numComponents = 2;
//! Index of the wetting phase's component
static const int wCompIdx = 0;
static const int wettingCompIdx = 0;
//! Index of the non-wetting phase's component
static const int nCompIdx = 1;
static const int nonWettingCompIdx = 1;
//! \copydoc BaseFluidSystem::componentName
static const char *componentName(int compIdx)
{
assert(0 <= compIdx && compIdx < numComponents);
if (compIdx == wCompIdx)
if (compIdx == wettingCompIdx)
return WettingPhase::name();
return NonwettingPhase::name();
}
@@ -159,7 +159,7 @@ public:
// let the fluids decide
return
(compIdx == wCompIdx)
(compIdx == wettingCompIdx)
? WettingPhase::molarMass()
: NonwettingPhase::molarMass();
}
@@ -172,7 +172,7 @@ public:
//assert(0 <= compIdx && compIdx < numComponents);
// let the fluids decide
return
(compIdx == wCompIdx)
(compIdx == wettingCompIdx)
? WettingPhase::criticalTemperature()
: NonwettingPhase::criticalTemperature();
}
@@ -185,7 +185,7 @@ public:
//assert(0 <= compIdx && compIdx < numComponents);
// let the fluids decide
return
(compIdx == wCompIdx)
(compIdx == wettingCompIdx)
? WettingPhase::criticalPressure()
: NonwettingPhase::criticalPressure();
}
@@ -198,7 +198,7 @@ public:
//assert(0 <= compIdx && compIdx < numComponents);
// let the fluids decide
return
(compIdx == wCompIdx)
(compIdx == wettingCompIdx)
? WettingPhase::acentricFactor()
: NonwettingPhase::acentricFactor();
}
@@ -225,7 +225,7 @@ public:
Scalar temperature = fluidState.temperature(phaseIdx);
Scalar pressure = fluidState.pressure(phaseIdx);
if (phaseIdx == wPhaseIdx)
if (phaseIdx == wettingPhaseIdx)
return WettingPhase::density(temperature, pressure);
return NonwettingPhase::density(temperature, pressure);
}
@@ -240,7 +240,7 @@ public:
Scalar temperature = fluidState.temperature(phaseIdx);
Scalar pressure = fluidState.pressure(phaseIdx);
if (phaseIdx == wPhaseIdx)
if (phaseIdx == wettingPhaseIdx)
return WettingPhase::viscosity(temperature, pressure);
return NonwettingPhase::viscosity(temperature, pressure);
}
@@ -274,7 +274,7 @@ public:
Scalar temperature = fluidState.temperature(phaseIdx);
Scalar pressure = fluidState.pressure(phaseIdx);
if (phaseIdx == wPhaseIdx)
if (phaseIdx == wettingPhaseIdx)
return WettingPhase::enthalpy(temperature, pressure);
return NonwettingPhase::enthalpy(temperature, pressure);
}
@@ -289,7 +289,7 @@ public:
Scalar temperature = fluidState.temperature(phaseIdx);
Scalar pressure = fluidState.pressure(phaseIdx);
if (phaseIdx == wPhaseIdx)
if (phaseIdx == wettingPhaseIdx)
return WettingPhase::thermalConductivity(temperature, pressure);
return NonwettingPhase::thermalConductivity(temperature, pressure);
}
@@ -304,7 +304,7 @@ public:
Scalar temperature = fluidState.temperature(phaseIdx);
Scalar pressure = fluidState.pressure(phaseIdx);
if (phaseIdx == wPhaseIdx)
if (phaseIdx == wettingPhaseIdx)
return WettingPhase::heatCapacity(temperature, pressure);
return NonwettingPhase::heatCapacity(temperature, pressure);
}

90
opm/material/fluidsystems/BlackOilFluidSystem.hpp
View File

@@ -63,11 +63,11 @@ public:
static const int numPhases = 3;
//! Index of the oil phase
static const int oPhaseIdx = 0;
static const int oilPhaseIdx = 0;
//! Index of the water phase
static const int wPhaseIdx = 1;
static const int waterPhaseIdx = 1;
//! Index of the gas phase
static const int gPhaseIdx = 2;
static const int gasPhaseIdx = 2;
/*!
* \copydoc BaseFluidSystem::init
@@ -113,9 +113,9 @@ public:
Scalar rhoWater,
Scalar rhoGas)
{
surfaceDensity_[oPhaseIdx] = rhoOil;
surfaceDensity_[wPhaseIdx] = rhoWater;
surfaceDensity_[gPhaseIdx] = rhoGas;
surfaceDensity_[oilPhaseIdx] = rhoOil;
surfaceDensity_[waterPhaseIdx] = rhoWater;
surfaceDensity_[gasPhaseIdx] = rhoGas;
}
/*!
@@ -223,18 +223,18 @@ public:
// calculate molar masses
// water is simple: 18 g/mol
molarMass_[wCompIdx] = 18e-3;
molarMass_[waterCompIdx] = 18e-3;
// for gas, we take the density at surface pressure and assume
// it to be ideal
Scalar p = 1.0135e5;
Scalar rho_g = surfaceDensity_[gPhaseIdx];
Scalar rho_g = surfaceDensity_[gasPhaseIdx];
Scalar T = 311;
molarMass_[gCompIdx] = Opm::Constants<Scalar>::R*T*rho_g / p;
molarMass_[gasCompIdx] = Opm::Constants<Scalar>::R*T*rho_g / p;
// finally, for oil phase, we take the molar mass from the
// spe9 paper
molarMass_[oCompIdx] = 175e-3; // kg/mol
molarMass_[oilCompIdx] = 175e-3; // kg/mol
// create the spline representing saturation pressure
// depending of the mass fraction in gas
@@ -290,7 +290,7 @@ public:
static bool isLiquid(const int phaseIdx)
{
assert(0 <= phaseIdx && phaseIdx < numPhases);
return phaseIdx != gPhaseIdx;
return phaseIdx != gasPhaseIdx;
}
/****************************************
@@ -301,11 +301,11 @@ public:
static const int numComponents = 3;
//! Index of the oil component
static const int oCompIdx = 0;
static const int oilCompIdx = 0;
//! Index of the water component
static const int wCompIdx = 1;
static const int waterCompIdx = 1;
//! Index of the gas component
static const int gCompIdx = 2;
static const int gasCompIdx = 2;
//! \copydoc BaseFluidSystem::componentName
static const char *componentName(const int compIdx)
@@ -350,10 +350,10 @@ public:
Scalar p = fluidState.pressure(phaseIdx);
switch (phaseIdx) {
case wPhaseIdx: return waterDensity_(p);
case gPhaseIdx: return gasDensity_(p);
case oPhaseIdx: {
Scalar pSat = oilSaturationPressure(fluidState.massFraction(oPhaseIdx, gCompIdx));
case waterPhaseIdx: return waterDensity_(p);
case gasPhaseIdx: return gasDensity_(p);
case oilPhaseIdx: {
Scalar pSat = oilSaturationPressure(fluidState.massFraction(oilPhaseIdx, gasCompIdx));
// retrieve the gas formation factor and the oil formation volume factor
Scalar Rs = gasDissolutionFactorSpline_.eval(pSat, /*extrapolate=*/true);
@@ -367,19 +367,19 @@ public:
// define the derivatives of oil regarding oil component
// mass fraction and pressure
Scalar drhoo_dXoO =
surfaceDensity_[oPhaseIdx]*
surfaceDensity_[oilPhaseIdx]*
(1 + oilCompressibility()*(p - 1.0135e5));
Scalar drhoo_dp =
oilCompressibility();
// Calculate the derivative of the density of saturated
// oil regarding pressure
Scalar drhoosat_dp = - surfaceDensity_[oPhaseIdx]*dBo_dp / (Bo * Bo);
Scalar drhoosat_dp = - surfaceDensity_[oilPhaseIdx]*dBo_dp / (Bo * Bo);
// calculate the derivative of the gas component mass
// fraction regarding pressure in saturated oil
Scalar dXoOsat_dp =
- surfaceDensity_[gPhaseIdx]/surfaceDensity_[oPhaseIdx]
- surfaceDensity_[gasPhaseIdx]/surfaceDensity_[oilPhaseIdx]
*(Bo * dRs_dp + Rs * dBo_dp);
// Using the previous derivatives, define a derivative
@@ -388,13 +388,13 @@ public:
drhoo_dXoO + (drhoo_dp - drhoosat_dp) / dXoOsat_dp;
// calculate the composition of saturated oil.
Scalar XoGsat = surfaceDensity_[gPhaseIdx]/surfaceDensity_[oPhaseIdx] * Rs * Bo;
Scalar XoGsat = surfaceDensity_[gasPhaseIdx]/surfaceDensity_[oilPhaseIdx] * Rs * Bo;
Scalar XoOsat = 1.0 - XoGsat;
Scalar rhoo =
surfaceDensity_[oPhaseIdx]/Bo*(1 + drhoo_dp*(p - pSat))
+ (XoOsat - fluidState.massFraction(oPhaseIdx, oCompIdx))*drhoo_dXoO
+ (XoGsat - fluidState.massFraction(oPhaseIdx, gCompIdx))*drhoo_dXoG;
surfaceDensity_[oilPhaseIdx]/Bo*(1 + drhoo_dp*(p - pSat))
+ (XoOsat - fluidState.massFraction(oilPhaseIdx, oilCompIdx))*drhoo_dXoO
+ (XoGsat - fluidState.massFraction(oilPhaseIdx, gasCompIdx))*drhoo_dXoG;
return rhoo;
}
@@ -415,9 +415,9 @@ public:
Scalar p = fluidState.pressure(phaseIdx);
switch (phaseIdx) {
case wPhaseIdx: return fugCoefficientInWater(compIdx, p);
case gPhaseIdx: return fugCoefficientInGas(compIdx, p);
case oPhaseIdx: return fugCoefficientInOil(compIdx, p);
case waterPhaseIdx: return fugCoefficientInWater(compIdx, p);
case gasPhaseIdx: return fugCoefficientInGas(compIdx, p);
case oilPhaseIdx: return fugCoefficientInOil(compIdx, p);
}
OPM_THROW(std::logic_error, "Unhandled phase or component index");
@@ -434,9 +434,9 @@ public:
Scalar p = fluidState.pressure(phaseIdx);
switch (phaseIdx) {
case oPhaseIdx: return oilViscositySpline_.eval(p, /*extrapolate=*/true);
case wPhaseIdx: return waterViscosity_(p);
case gPhaseIdx: return gasViscosity_(p);
case oilPhaseIdx: return oilViscositySpline_.eval(p, /*extrapolate=*/true);
case waterPhaseIdx: return waterViscosity_(p);
case gasPhaseIdx: return gasViscosity_(p);
}
OPM_THROW(std::logic_error, "Unhandled phase index " << phaseIdx);
@@ -501,9 +501,9 @@ public:
// pressure of water as a starting point. (we just set it to
// 30 kPa to ease interpreting the results.)
const Scalar pvWater = 30e3;
if (compIdx == oCompIdx)
if (compIdx == oilCompIdx)
return 1e3*pvWater / pressure;
else if (compIdx == gCompIdx)
else if (compIdx == gasCompIdx)
return 1e6*pvWater / pressure;
return pvWater / pressure;
@@ -534,9 +534,9 @@ public:
// pressure to ease physical interpretation of the results
Scalar phi_oO = 20e3/pressure;
if (compIdx == oCompIdx)
if (compIdx == oilCompIdx)
return phi_oO;
else if (compIdx == wCompIdx)
else if (compIdx == waterCompIdx)
// assume that the affinity of the water component to the
// oil phase is one million times smaller than that of the
// oil component
@@ -553,7 +553,7 @@ public:
// then, scale the gas component's gas phase fugacity
// coefficient, so that the oil phase ends up at the right
// composition if we were doing a flash experiment
Scalar phi_gG = fugCoefficientInGas(gCompIdx, pressure);
Scalar phi_gG = fugCoefficientInGas(gasCompIdx, pressure);
return phi_gG / x_oGf;
}
@@ -603,10 +603,10 @@ public:
{
// first, we calculate the total reservoir oil phase density
// [kg/m^3]
Scalar rho_gRef = surfaceDensity_[gPhaseIdx];
Scalar rho_gRef = surfaceDensity_[gasPhaseIdx];
Scalar B_o = oilFormationVolumeFactor(pressure);
Scalar rho_o = surfaceDensity(oPhaseIdx)/B_o;
Scalar rho_o = surfaceDensity(oilPhaseIdx)/B_o;
// then, we calculate the mass of the gas component [kg/m^3]
// in the oil phase. This is equivalent to the gas formation
@@ -629,8 +629,8 @@ public:
// which can be converted to mole fractions, given the
// components' molar masses
Scalar MG = molarMass(gCompIdx);
Scalar MO = molarMass(oCompIdx);
Scalar MG = molarMass(gasCompIdx);
Scalar MO = molarMass(oilCompIdx);
Scalar avgMolarMass = MO/(1 + XoG*(MO/MG - 1));
Scalar xoG = XoG*avgMolarMass/MG;
@@ -645,10 +645,10 @@ public:
Scalar Bo = oilFormationVolumeFactor(pressure);
// oil formation volume factor at standard pressure
Scalar BoRef = refFormationVolumeFactor_[oPhaseIdx];
Scalar BoRef = refFormationVolumeFactor_[oilPhaseIdx];
// surface density of oil
Scalar rhoRef = surfaceDensity_[oPhaseIdx];
Scalar rhoRef = surfaceDensity_[oilPhaseIdx];
// reservoir density is surface density scaled by the ratio of
// the volume formation factors
@@ -663,10 +663,10 @@ private:
// gas formation volume factor at standard pressure
Scalar BgRef = refFormationVolumeFactor_[gPhaseIdx];
Scalar BgRef = refFormationVolumeFactor_[gasPhaseIdx];
// surface density of gas
Scalar rhoRef = surfaceDensity_[gPhaseIdx];
Scalar rhoRef = surfaceDensity_[gasPhaseIdx];
// reservoir density is surface density scaled by the ratio of
// the volume formation factors
@@ -676,7 +676,7 @@ private:
static Scalar waterDensity_(Scalar pressure)
{
// compressibility of water times standard density
Scalar rhoRef = surfaceDensity_[wPhaseIdx];
Scalar rhoRef = surfaceDensity_[waterPhaseIdx];
return rhoRef *
(1 +
waterCompressibilityScalar_

40
opm/material/fluidsystems/BrineCO2FluidSystem.hpp
View File

@@ -81,9 +81,9 @@ public:
static const int numPhases = 2;
//! The index of the liquid phase
static const int lPhaseIdx = 0;
static const int liquidPhaseIdx = 0;
//! The index of the gas phase
static const int gPhaseIdx = 1;
static const int gasPhaseIdx = 1;
/*!
* \copydoc BaseFluidSystem::phaseName
@@ -106,7 +106,7 @@ public:
{
assert(0 <= phaseIdx && phaseIdx < numPhases);
return phaseIdx != gPhaseIdx;
return phaseIdx != gasPhaseIdx;
}
/*!
@@ -116,7 +116,7 @@ public:
{
assert(0 <= phaseIdx && phaseIdx < numPhases);
if (phaseIdx == gPhaseIdx)
if (phaseIdx == gasPhaseIdx)
return CO2::gasIsIdeal();
return false;
}
@@ -236,12 +236,12 @@ public:
Scalar temperature = fluidState.temperature(phaseIdx);
Scalar pressure = fluidState.pressure(phaseIdx);
if (phaseIdx == lPhaseIdx) {
if (phaseIdx == liquidPhaseIdx) {
// use normalized composition for to calculate the density
// (the relations don't seem to take non-normalized
// compositions too well...)
Scalar xlBrine = std::min(1.0, std::max(0.0, fluidState.moleFraction(lPhaseIdx, BrineIdx)));
Scalar xlCO2 = std::min(1.0, std::max(0.0, fluidState.moleFraction(lPhaseIdx, CO2Idx)));
Scalar xlBrine = std::min(1.0, std::max(0.0, fluidState.moleFraction(liquidPhaseIdx, BrineIdx)));
Scalar xlCO2 = std::min(1.0, std::max(0.0, fluidState.moleFraction(liquidPhaseIdx, CO2Idx)));
Scalar sumx = xlBrine + xlCO2;
xlBrine /= sumx;
xlCO2 /= sumx;
@@ -255,13 +255,13 @@ public:
return result;
}
assert(phaseIdx == gPhaseIdx);
assert(phaseIdx == gasPhaseIdx);
// use normalized composition for to calculate the density
// (the relations don't seem to take non-normalized
// compositions too well...)
Scalar xgBrine = std::min(1.0, std::max(0.0, fluidState.moleFraction(gPhaseIdx, BrineIdx)));
Scalar xgCO2 = std::min(1.0, std::max(0.0, fluidState.moleFraction(gPhaseIdx, CO2Idx)));
Scalar xgBrine = std::min(1.0, std::max(0.0, fluidState.moleFraction(gasPhaseIdx, BrineIdx)));
Scalar xgCO2 = std::min(1.0, std::max(0.0, fluidState.moleFraction(gasPhaseIdx, CO2Idx)));
Scalar sumx = xgBrine + xgCO2;
xgBrine /= sumx;
xgCO2 /= sumx;
@@ -287,7 +287,7 @@ public:
Scalar temperature = fluidState.temperature(phaseIdx);
Scalar pressure = fluidState.pressure(phaseIdx);
if (phaseIdx == lPhaseIdx) {
if (phaseIdx == liquidPhaseIdx) {
// assume pure brine for the liquid phase. TODO: viscosity
// of mixture
Scalar result = Brine::liquidViscosity(temperature, pressure);
@@ -295,7 +295,7 @@ public:
return result;
}
assert(phaseIdx == gPhaseIdx);
assert(phaseIdx == gasPhaseIdx);
Scalar result = CO2::gasViscosity(temperature, pressure);
Valgrind::CheckDefined(result);
return result;
@@ -313,7 +313,7 @@ public:
assert(0 <= phaseIdx && phaseIdx < numPhases);
assert(0 <= compIdx && compIdx < numComponents);
if (phaseIdx == gPhaseIdx)
if (phaseIdx == gasPhaseIdx)
// use the fugacity coefficients of an ideal gas. the
// actual value of the fugacity is not relevant, as long
// as the relative fluid compositions are observed,
@@ -366,10 +366,10 @@ public:
{
Scalar temperature = fluidState.temperature(phaseIdx);
Scalar pressure = fluidState.pressure(phaseIdx);
if (phaseIdx == lPhaseIdx)
if (phaseIdx == liquidPhaseIdx)
return BinaryCoeffBrineCO2::liquidDiffCoeff(temperature, pressure);
assert(phaseIdx == gPhaseIdx);
assert(phaseIdx == gasPhaseIdx);
return BinaryCoeffBrineCO2::gasDiffCoeff(temperature, pressure);
}
@@ -386,7 +386,7 @@ public:
Scalar temperature = fluidState.temperature(phaseIdx);
Scalar pressure = fluidState.pressure(phaseIdx);
if (phaseIdx == lPhaseIdx) {
if (phaseIdx == liquidPhaseIdx) {
Scalar XlCO2 = fluidState.massFraction(phaseIdx, CO2Idx);
Scalar result = liquidEnthalpyBrineCO2_(temperature,
pressure,
@@ -396,8 +396,8 @@ public:
return result;
}
else {
Scalar XCO2 = fluidState.massFraction(gPhaseIdx, CO2Idx);
Scalar XBrine = fluidState.massFraction(gPhaseIdx, BrineIdx);
Scalar XCO2 = fluidState.massFraction(gasPhaseIdx, CO2Idx);
Scalar XBrine = fluidState.massFraction(gasPhaseIdx, BrineIdx);
Scalar result = 0;
result += XBrine * Brine::gasEnthalpy(temperature, pressure);
@@ -416,7 +416,7 @@ public:
int phaseIdx)
{
// TODO way too simple!
if (phaseIdx == lPhaseIdx)
if (phaseIdx == liquidPhaseIdx)
return 0.6; // conductivity of water[W / (m K ) ]
// gas phase
@@ -440,7 +440,7 @@ public:
const ParameterCache &paramCache,
int phaseIdx)
{
if(phaseIdx == lPhaseIdx)
if(phaseIdx == liquidPhaseIdx)
return H2O::liquidHeatCapacity(fluidState.temperature(phaseIdx),
fluidState.pressure(phaseIdx));
else

52
opm/material/fluidsystems/H2OAirFluidSystem.hpp
View File

@@ -75,16 +75,16 @@ public:
static const int numPhases = 2;
//! The index of the liquid phase
static const int lPhaseIdx = 0;
static const int liquidPhaseIdx = 0;
//! The index of the gas phase
static const int gPhaseIdx = 1;
static const int gasPhaseIdx = 1;
//! \copydoc BaseFluidSystem::phaseName
static const char *phaseName(int phaseIdx)
{
switch (phaseIdx) {
case lPhaseIdx: return "liquid";
case gPhaseIdx: return "gas";
case liquidPhaseIdx: return "liquid";
case gasPhaseIdx: return "gas";
};
OPM_THROW(std::logic_error, "Invalid phase index " << phaseIdx);
}
@@ -93,14 +93,14 @@ public:
static bool isLiquid(int phaseIdx)
{
//assert(0 <= phaseIdx && phaseIdx < numPhases);
return phaseIdx != gPhaseIdx;
return phaseIdx != gasPhaseIdx;
}
//! \copydoc BaseFluidSystem::isCompressible
static bool isCompressible(int phaseIdx)
{
//assert(0 <= phaseIdx && phaseIdx < numPhases);
return (phaseIdx == gPhaseIdx)
return (phaseIdx == gasPhaseIdx)
// ideal gases are always compressible
? true
:
@@ -112,7 +112,7 @@ public:
static bool isIdealGas(int phaseIdx)
{
return
(phaseIdx == gPhaseIdx)
(phaseIdx == gasPhaseIdx)
? H2O::gasIsIdeal() && Air::gasIsIdeal()
: false;
}
@@ -272,7 +272,7 @@ public:
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
sumMoleFrac += fluidState.moleFraction(phaseIdx, compIdx);
if (phaseIdx == lPhaseIdx)
if (phaseIdx == liquidPhaseIdx)
{
if (!useComplexRelations)
// assume pure water
@@ -285,28 +285,28 @@ public:
return
clH2O
* (H2O::molarMass()*fluidState.moleFraction(lPhaseIdx, H2OIdx)
* (H2O::molarMass()*fluidState.moleFraction(liquidPhaseIdx, H2OIdx)
+
Air::molarMass()*fluidState.moleFraction(lPhaseIdx, AirIdx))
Air::molarMass()*fluidState.moleFraction(liquidPhaseIdx, AirIdx))
/ sumMoleFrac;
}
}
else if (phaseIdx == gPhaseIdx)
else if (phaseIdx == gasPhaseIdx)
{
if (!useComplexRelations)
// for the gas phase assume an ideal gas
return
IdealGas::molarDensity(T, p)
* fluidState.averageMolarMass(gPhaseIdx)
* fluidState.averageMolarMass(gasPhaseIdx)
/ std::max(1e-5, sumMoleFrac);
Scalar partialPressureH2O =
fluidState.moleFraction(gPhaseIdx, H2OIdx) *
fluidState.pressure(gPhaseIdx);
fluidState.moleFraction(gasPhaseIdx, H2OIdx) *
fluidState.pressure(gasPhaseIdx);
Scalar partialPressureAir =
fluidState.moleFraction(gPhaseIdx, AirIdx) *
fluidState.pressure(gPhaseIdx);
fluidState.moleFraction(gasPhaseIdx, AirIdx) *
fluidState.pressure(gasPhaseIdx);
return
H2O::gasDensity(T, partialPressureH2O) +
@@ -326,14 +326,14 @@ public:
Scalar T = fluidState.temperature(phaseIdx);
Scalar p = fluidState.pressure(phaseIdx);
if (phaseIdx == lPhaseIdx)
if (phaseIdx == liquidPhaseIdx)
{
// assume pure water for the liquid phase
// TODO: viscosity of mixture
// couldn't find a way to solve the mixture problem
return H2O::liquidViscosity(T, p);
}
else if (phaseIdx == gPhaseIdx)
else if (phaseIdx == gasPhaseIdx)
{
if(!useComplexRelations){
return Air::gasViscosity(T, p);
@@ -394,7 +394,7 @@ public:
Scalar T = fluidState.temperature(phaseIdx);
Scalar p = fluidState.pressure(phaseIdx);
if (phaseIdx == lPhaseIdx) {
if (phaseIdx == liquidPhaseIdx) {
if (compIdx == H2OIdx)
return H2O::vaporPressure(T)/p;
return Opm::BinaryCoeff::H2O_Air::henry(T)/p;
@@ -415,10 +415,10 @@ public:
Scalar T = fluidState.temperature(phaseIdx);
Scalar p = fluidState.pressure(phaseIdx);
if (phaseIdx == lPhaseIdx)
if (phaseIdx == liquidPhaseIdx)
return BinaryCoeff::H2O_Air::liquidDiffCoeff(T, p);
assert(phaseIdx == gPhaseIdx);
assert(phaseIdx == gasPhaseIdx);
return BinaryCoeff::H2O_Air::gasDiffCoeff(T, p);
}
@@ -433,22 +433,22 @@ public:
Valgrind::CheckDefined(T);
Valgrind::CheckDefined(p);
if (phaseIdx == lPhaseIdx)
if (phaseIdx == liquidPhaseIdx)
{
// TODO: correct way to deal with the solutes???
return H2O::liquidEnthalpy(T, p);
}
else if (phaseIdx == gPhaseIdx)
else if (phaseIdx == gasPhaseIdx)
{
Scalar result = 0.0;
result +=
H2O::gasEnthalpy(T, p) *
fluidState.massFraction(gPhaseIdx, H2OIdx);
fluidState.massFraction(gasPhaseIdx, H2OIdx);
result +=
Air::gasEnthalpy(T, p) *
fluidState.massFraction(gPhaseIdx, AirIdx);
fluidState.massFraction(gasPhaseIdx, AirIdx);
return result;
}
OPM_THROW(std::logic_error, "Invalid phase index " << phaseIdx);
@@ -465,7 +465,7 @@ public:
Scalar temperature = fluidState.temperature(phaseIdx) ;
Scalar pressure = fluidState.pressure(phaseIdx);
if (phaseIdx == lPhaseIdx){// liquid phase
if (phaseIdx == liquidPhaseIdx){// liquid phase
return H2O::liquidThermalConductivity(temperature, pressure);
}
else{// gas phase

114
opm/material/fluidsystems/H2OAirMesityleneFluidSystem.hpp
View File

@@ -80,11 +80,11 @@ public:
static const int numComponents = 3;
//! The index of the water phase
static const int wPhaseIdx = 0;
static const int waterPhaseIdx = 0;
//! The index of the NAPL phase
static const int nPhaseIdx = 1;
static const int naplPhaseIdx = 1;
//! The index of the gas phase
static const int gPhaseIdx = 2;
static const int gasPhaseIdx = 2;
//! The index of the water component
static const int H2OIdx = 0;
@@ -128,21 +128,21 @@ public:
static bool isLiquid(int phaseIdx)
{
//assert(0 <= phaseIdx && phaseIdx < numPhases);
return phaseIdx != gPhaseIdx;
return phaseIdx != gasPhaseIdx;
}
//! \copydoc BaseFluidSystem::isIdealGas
static bool isIdealGas(int phaseIdx)
{ return phaseIdx == gPhaseIdx && H2O::gasIsIdeal() && Air::gasIsIdeal() && NAPL::gasIsIdeal(); }
{ return phaseIdx == gasPhaseIdx && H2O::gasIsIdeal() && Air::gasIsIdeal() && NAPL::gasIsIdeal(); }
//! \copydoc BaseFluidSystem::isCompressible
static bool isCompressible(int phaseIdx)
{
//assert(0 <= phaseIdx && phaseIdx < numPhases);
// gases are always compressible
return (phaseIdx == gPhaseIdx)
return (phaseIdx == gasPhaseIdx)
? true
: (phaseIdx == wPhaseIdx)
: (phaseIdx == waterPhaseIdx)
? H2O::liquidIsCompressible()
: NAPL::liquidIsCompressible();
}
@@ -161,9 +161,9 @@ public:
static const char *phaseName(int phaseIdx)
{
switch (phaseIdx) {
case wPhaseIdx: return "w";
case nPhaseIdx: return "n";
case gPhaseIdx: return "g";;
case waterPhaseIdx: return "w";
case naplPhaseIdx: return "n";
case gasPhaseIdx: return "g";;
};
OPM_THROW(std::logic_error, "Invalid phase index " << phaseIdx);
}
@@ -201,7 +201,7 @@ public:
{
Scalar T = fluidState.temperature(phaseIdx) ;
if (phaseIdx == wPhaseIdx) {
if (phaseIdx == waterPhaseIdx) {
// See: Ochs 2008
Scalar p = H2O::liquidIsCompressible()?fluidState.pressure(phaseIdx):1e100;
Scalar rholH2O = H2O::liquidDensity(T, p);
@@ -210,28 +210,28 @@ public:
// this assumes each dissolved molecule displaces exactly one
// water molecule in the liquid
return
clH2O*(H2O::molarMass()*fluidState.moleFraction(wPhaseIdx, H2OIdx)
clH2O*(H2O::molarMass()*fluidState.moleFraction(waterPhaseIdx, H2OIdx)
+
Air::molarMass()*fluidState.moleFraction(wPhaseIdx, airIdx)
Air::molarMass()*fluidState.moleFraction(waterPhaseIdx, airIdx)
+
NAPL::molarMass()*fluidState.moleFraction(wPhaseIdx, NAPLIdx));
NAPL::molarMass()*fluidState.moleFraction(waterPhaseIdx, NAPLIdx));
}
else if (phaseIdx == nPhaseIdx) {
else if (phaseIdx == naplPhaseIdx) {
// assume pure NAPL for the NAPL phase
Scalar p = NAPL::liquidIsCompressible()?fluidState.pressure(phaseIdx):1e100;
return NAPL::liquidDensity(T, p);
}
assert (phaseIdx == gPhaseIdx);
assert (phaseIdx == gasPhaseIdx);
Scalar pH2O =
fluidState.moleFraction(gPhaseIdx, H2OIdx) *
fluidState.pressure(gPhaseIdx);
fluidState.moleFraction(gasPhaseIdx, H2OIdx) *
fluidState.pressure(gasPhaseIdx);
Scalar pAir =
fluidState.moleFraction(gPhaseIdx, airIdx) *
fluidState.pressure(gPhaseIdx);
fluidState.moleFraction(gasPhaseIdx, airIdx) *
fluidState.pressure(gasPhaseIdx);
Scalar pNAPL =
fluidState.moleFraction(gPhaseIdx, NAPLIdx) *
fluidState.pressure(gPhaseIdx);
fluidState.moleFraction(gasPhaseIdx, NAPLIdx) *
fluidState.pressure(gasPhaseIdx);
return
H2O::gasDensity(T, pH2O) +
Air::gasDensity(T, pAir) +
@@ -247,18 +247,18 @@ public:
Scalar T = fluidState.temperature(phaseIdx);
Scalar p = fluidState.pressure(phaseIdx);
if (phaseIdx == wPhaseIdx) {
if (phaseIdx == waterPhaseIdx) {
// assume pure water viscosity
return H2O::liquidViscosity(T,
p);
}
else if (phaseIdx == nPhaseIdx) {
else if (phaseIdx == naplPhaseIdx) {
// assume pure NAPL viscosity
return NAPL::liquidViscosity(T, p);
}
assert (phaseIdx == gPhaseIdx);
assert (phaseIdx == gasPhaseIdx);
/* Wilke method. See:
*
@@ -284,15 +284,15 @@ public:
NAPL::molarMass()
};
Scalar muAW = mu[airIdx]*fluidState.moleFraction(gPhaseIdx, airIdx)
+ mu[H2OIdx]*fluidState.moleFraction(gPhaseIdx, H2OIdx)
/ (fluidState.moleFraction(gPhaseIdx, airIdx)
+ fluidState.moleFraction(gPhaseIdx, H2OIdx));
Scalar xAW = fluidState.moleFraction(gPhaseIdx, airIdx)
+ fluidState.moleFraction(gPhaseIdx, H2OIdx);
Scalar muAW = mu[airIdx]*fluidState.moleFraction(gasPhaseIdx, airIdx)
+ mu[H2OIdx]*fluidState.moleFraction(gasPhaseIdx, H2OIdx)
/ (fluidState.moleFraction(gasPhaseIdx, airIdx)
+ fluidState.moleFraction(gasPhaseIdx, H2OIdx));
Scalar xAW = fluidState.moleFraction(gasPhaseIdx, airIdx)
+ fluidState.moleFraction(gasPhaseIdx, H2OIdx);
Scalar MAW = (fluidState.moleFraction(gPhaseIdx, airIdx)*Air::molarMass()
+ fluidState.moleFraction(gPhaseIdx, H2OIdx)*H2O::molarMass())
Scalar MAW = (fluidState.moleFraction(gasPhaseIdx, airIdx)*Air::molarMass()
+ fluidState.moleFraction(gasPhaseIdx, H2OIdx)*H2O::molarMass())
/ xAW;
Scalar phiCAW = 0.3; // simplification for this particular system
@@ -302,9 +302,9 @@ public:
*/
Scalar phiAWC = phiCAW * muAW*M[NAPLIdx]/(mu[NAPLIdx]*MAW);
muResult = (xAW*muAW)/(xAW+fluidState.moleFraction(gPhaseIdx, NAPLIdx)*phiAWC)
+ (fluidState.moleFraction(gPhaseIdx, NAPLIdx) * mu[NAPLIdx])
/ (fluidState.moleFraction(gPhaseIdx, NAPLIdx) + xAW*phiCAW);
muResult = (xAW*muAW)/(xAW+fluidState.moleFraction(gasPhaseIdx, NAPLIdx)*phiAWC)
+ (fluidState.moleFraction(gasPhaseIdx, NAPLIdx) * mu[NAPLIdx])
/ (fluidState.moleFraction(gasPhaseIdx, NAPLIdx) + xAW*phiCAW);
return muResult;
}
@@ -321,14 +321,14 @@ public:
Scalar p = fluidState.pressure(phaseIdx);
Scalar diffCont;
if (phaseIdx==gPhaseIdx) {
if (phaseIdx==gasPhaseIdx) {
Scalar diffAC = Opm::BinaryCoeff::Air_Mesitylene::gasDiffCoeff(T, p);
Scalar diffWC = Opm::BinaryCoeff::H2O_Mesitylene::gasDiffCoeff(T, p);
Scalar diffAW = Opm::BinaryCoeff::H2O_Air::gasDiffCoeff(T, p);
const Scalar xga = fluidState.moleFraction(gPhaseIdx, airIdx);
const Scalar xgw = fluidState.moleFraction(gPhaseIdx, H2OIdx);
const Scalar xgc = fluidState.moleFraction(gPhaseIdx, NAPLIdx);
const Scalar xga = fluidState.moleFraction(gasPhaseIdx, airIdx);
const Scalar xgw = fluidState.moleFraction(gasPhaseIdx, H2OIdx);
const Scalar xgc = fluidState.moleFraction(gasPhaseIdx, NAPLIdx);
if (compIdx==NAPLIdx) return (1 - xgw)/(xga/diffAW + xgc/diffWC);
else if (compIdx==H2OIdx) return (1 - xgc)/(xgw/diffWC + xga/diffAC);
@@ -336,14 +336,14 @@ public:
"Diffusivity of air in the gas phase "
"is constraint by sum of diffusive fluxes = 0 !\n");
}
else if (phaseIdx==wPhaseIdx){
else if (phaseIdx==waterPhaseIdx){
Scalar diffACl = 1.e-9; // BinaryCoeff::Air_Mesitylene::liquidDiffCoeff(temperature, pressure);
Scalar diffWCl = 1.e-9; // BinaryCoeff::H2O_Mesitylene::liquidDiffCoeff(temperature, pressure);
Scalar diffAWl = 1.e-9; // BinaryCoeff::H2O_Air::liquidDiffCoeff(temperature, pressure);
Scalar xwa = fluidState.moleFraction(wPhaseIdx, airIdx);
Scalar xww = fluidState.moleFraction(wPhaseIdx, H2OIdx);
Scalar xwc = fluidState.moleFraction(wPhaseIdx, NAPLIdx);
Scalar xwa = fluidState.moleFraction(waterPhaseIdx, airIdx);
Scalar xww = fluidState.moleFraction(waterPhaseIdx, H2OIdx);
Scalar xwc = fluidState.moleFraction(waterPhaseIdx, NAPLIdx);
switch (compIdx) {
case NAPLIdx:
@@ -358,7 +358,7 @@ public:
"is constraint by sum of diffusive fluxes = 0 !\n");
};
}
else if (phaseIdx==nPhaseIdx) {
else if (phaseIdx==naplPhaseIdx) {
OPM_THROW(std::logic_error,
"Diffusion coefficients of "
"substances in liquid phase are undefined!\n");
@@ -382,7 +382,7 @@ public:
Valgrind::CheckDefined(T);
Valgrind::CheckDefined(p);
if (phaseIdx == wPhaseIdx) {
if (phaseIdx == waterPhaseIdx) {
if (compIdx == H2OIdx)
return H2O::vaporPressure(T)/p;
else if (compIdx == airIdx)
@@ -396,7 +396,7 @@ public:
// component to the NAPL phase is much higher than for the
// other components, i.e. the fugacity cofficient is much
// smaller.
else if (phaseIdx == nPhaseIdx) {
else if (phaseIdx == naplPhaseIdx) {
Scalar phiNapl = NAPL::vaporPressure(T)/p;
if (compIdx == NAPLIdx)
return phiNapl;
@@ -409,7 +409,7 @@ public:
// for the gas phase, assume an ideal gas when it comes to
// fugacity (-> fugacity == partial pressure)
assert(phaseIdx == gPhaseIdx);
assert(phaseIdx == gasPhaseIdx);
return 1.0;
}
@@ -423,18 +423,18 @@ public:
Scalar T = fluidState.temperature(phaseIdx) ;
Scalar p = fluidState.pressure(phaseIdx);
if (phaseIdx == wPhaseIdx) {
if (phaseIdx == waterPhaseIdx) {
return H2O::liquidEnthalpy(T, p);
}
else if (phaseIdx == nPhaseIdx) {
else if (phaseIdx == naplPhaseIdx) {
return NAPL::liquidEnthalpy(T, p);
}
else if (phaseIdx == gPhaseIdx) {
else if (phaseIdx == gasPhaseIdx) {
// gas phase enthalpy depends strongly on composition
Scalar result = 0;
result += H2O::gasEnthalpy(T, p) * fluidState.massFraction(gPhaseIdx, H2OIdx);
result += NAPL::gasEnthalpy(T, p) * fluidState.massFraction(gPhaseIdx, airIdx);
result += Air::gasEnthalpy(T, p) * fluidState.massFraction(gPhaseIdx, NAPLIdx);
result += H2O::gasEnthalpy(T, p) * fluidState.massFraction(gasPhaseIdx, H2OIdx);
result += NAPL::gasEnthalpy(T, p) * fluidState.massFraction(gasPhaseIdx, airIdx);
result += Air::gasEnthalpy(T, p) * fluidState.massFraction(gasPhaseIdx, NAPLIdx);
return result;
}
@@ -452,15 +452,15 @@ public:
Scalar T = fluidState.temperature(phaseIdx) ;
Scalar p = fluidState.pressure(phaseIdx);
if (phaseIdx == wPhaseIdx){ // water phase
if (phaseIdx == waterPhaseIdx){ // water phase
return H2O::liquidThermalConductivity(T, p);
}
else if (phaseIdx == gPhaseIdx) { // gas phase
else if (phaseIdx == gasPhaseIdx) { // gas phase
Scalar lambdaDryAir = Air::gasThermalConductivity(T, p);
return lambdaDryAir;
}
assert(phaseIdx == nPhaseIdx);
assert(phaseIdx == naplPhaseIdx);
// Taken from:
//

108
opm/material/fluidsystems/H2OAirXyleneFluidSystem.hpp
View File

@@ -69,11 +69,11 @@ public:
static const int numComponents = 3;
//! The index of the water phase
static const int wPhaseIdx = 0;
static const int waterPhaseIdx = 0;
//! The index of the NAPL phase
static const int nPhaseIdx = 1;
static const int naplPhaseIdx = 1;
//! The index of the gas phase
static const int gPhaseIdx = 2;
static const int gasPhaseIdx = 2;
//! The index of the water component
static const int H2OIdx = 0;
@@ -90,12 +90,12 @@ public:
static bool isLiquid(int phaseIdx)
{
//assert(0 <= phaseIdx && phaseIdx < numPhases);
return phaseIdx != gPhaseIdx;
return phaseIdx != gasPhaseIdx;
}
//! \copydoc BaseFluidSystem::isIdealGas
static bool isIdealGas(int phaseIdx)
{ return phaseIdx == gPhaseIdx && H2O::gasIsIdeal() && Air::gasIsIdeal() && NAPL::gasIsIdeal(); }
{ return phaseIdx == gasPhaseIdx && H2O::gasIsIdeal() && Air::gasIsIdeal() && NAPL::gasIsIdeal(); }
//! \copydoc BaseFluidSystem::isIdealMixture
static bool isIdealMixture(int phaseIdx)
@@ -112,10 +112,10 @@ public:
static bool isCompressible(int phaseIdx)
{
return
(phaseIdx == gPhaseIdx)
(phaseIdx == gasPhaseIdx)
// gases are always compressible
? true
: (phaseIdx == wPhaseIdx)
: (phaseIdx == waterPhaseIdx)
// the water component decides for the water phase...
? H2O::liquidIsCompressible()
// the NAPL component decides for the napl phase...
@@ -126,9 +126,9 @@ public:
static const char *phaseName(int phaseIdx)
{
switch (phaseIdx) {
case wPhaseIdx: return "w";
case nPhaseIdx: return "n";
case gPhaseIdx: return "g";;
case waterPhaseIdx: return "w";
case naplPhaseIdx: return "n";
case gasPhaseIdx: return "g";;
};
OPM_THROW(std::logic_error, "Invalid phase index " << phaseIdx);
}
@@ -166,7 +166,7 @@ public:
const ParameterCache &paramCache,
int phaseIdx)
{
if (phaseIdx == wPhaseIdx) {
if (phaseIdx == waterPhaseIdx) {
// See: Ochs 2008
// \todo: proper citation
Scalar rholH2O = H2O::liquidDensity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
@@ -175,28 +175,28 @@ public:
// this assumes each dissolved molecule displaces exactly one
// water molecule in the liquid
return
clH2O*(H2O::molarMass()*fluidState.moleFraction(wPhaseIdx, H2OIdx)
clH2O*(H2O::molarMass()*fluidState.moleFraction(waterPhaseIdx, H2OIdx)
+
Air::molarMass()*fluidState.moleFraction(wPhaseIdx, airIdx)
Air::molarMass()*fluidState.moleFraction(waterPhaseIdx, airIdx)
+
NAPL::molarMass()*fluidState.moleFraction(wPhaseIdx, NAPLIdx));
NAPL::molarMass()*fluidState.moleFraction(waterPhaseIdx, NAPLIdx));
}
else if (phaseIdx == nPhaseIdx) {
else if (phaseIdx == naplPhaseIdx) {
// assume pure NAPL for the NAPL phase
Scalar pressure = NAPL::liquidIsCompressible()?fluidState.pressure(phaseIdx):1e100;
return NAPL::liquidDensity(fluidState.temperature(phaseIdx), pressure);
}
assert (phaseIdx == gPhaseIdx);
assert (phaseIdx == gasPhaseIdx);
Scalar pH2O =
fluidState.moleFraction(gPhaseIdx, H2OIdx) *
fluidState.pressure(gPhaseIdx);
fluidState.moleFraction(gasPhaseIdx, H2OIdx) *
fluidState.pressure(gasPhaseIdx);
Scalar pAir =
fluidState.moleFraction(gPhaseIdx, airIdx) *
fluidState.pressure(gPhaseIdx);
fluidState.moleFraction(gasPhaseIdx, airIdx) *
fluidState.pressure(gasPhaseIdx);
Scalar pNAPL =
fluidState.moleFraction(gPhaseIdx, NAPLIdx) *
fluidState.pressure(gPhaseIdx);
fluidState.moleFraction(gasPhaseIdx, NAPLIdx) *
fluidState.pressure(gasPhaseIdx);
return
H2O::gasDensity(fluidState.temperature(phaseIdx), pH2O) +
Air::gasDensity(fluidState.temperature(phaseIdx), pAir) +
@@ -209,17 +209,17 @@ public:
const ParameterCache &paramCache,
int phaseIdx)
{
if (phaseIdx == wPhaseIdx) {
if (phaseIdx == waterPhaseIdx) {
// assume pure water viscosity
return H2O::liquidViscosity(fluidState.temperature(phaseIdx),
fluidState.pressure(phaseIdx));
}
else if (phaseIdx == nPhaseIdx) {
else if (phaseIdx == naplPhaseIdx) {
// assume pure NAPL viscosity
return NAPL::liquidViscosity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
}
assert (phaseIdx == gPhaseIdx);
assert (phaseIdx == gasPhaseIdx);
/* Wilke method. See:
*
@@ -245,15 +245,15 @@ public:
NAPL::molarMass()
};
Scalar muAW = mu[airIdx]*fluidState.moleFraction(gPhaseIdx, airIdx)
+ mu[H2OIdx]*fluidState.moleFraction(gPhaseIdx, H2OIdx)
/ (fluidState.moleFraction(gPhaseIdx, airIdx)
+ fluidState.moleFraction(gPhaseIdx, H2OIdx));
Scalar xAW = fluidState.moleFraction(gPhaseIdx, airIdx)
+ fluidState.moleFraction(gPhaseIdx, H2OIdx);
Scalar muAW = mu[airIdx]*fluidState.moleFraction(gasPhaseIdx, airIdx)
+ mu[H2OIdx]*fluidState.moleFraction(gasPhaseIdx, H2OIdx)
/ (fluidState.moleFraction(gasPhaseIdx, airIdx)
+ fluidState.moleFraction(gasPhaseIdx, H2OIdx));
Scalar xAW = fluidState.moleFraction(gasPhaseIdx, airIdx)
+ fluidState.moleFraction(gasPhaseIdx, H2OIdx);
Scalar MAW = (fluidState.moleFraction(gPhaseIdx, airIdx)*Air::molarMass()
+ fluidState.moleFraction(gPhaseIdx, H2OIdx)*H2O::molarMass())
Scalar MAW = (fluidState.moleFraction(gasPhaseIdx, airIdx)*Air::molarMass()
+ fluidState.moleFraction(gasPhaseIdx, H2OIdx)*H2O::molarMass())
/ xAW;
/* TODO, please check phiCAW for the Xylene case here */
@@ -266,9 +266,9 @@ public:
*/
Scalar phiAWC = phiCAW * muAW*M[NAPLIdx]/(mu[NAPLIdx]*MAW);
muResult = (xAW*muAW)/(xAW+fluidState.moleFraction(gPhaseIdx, NAPLIdx)*phiAWC)
+ (fluidState.moleFraction(gPhaseIdx, NAPLIdx) * mu[NAPLIdx])
/ (fluidState.moleFraction(gPhaseIdx, NAPLIdx) + xAW*phiCAW);
muResult = (xAW*muAW)/(xAW+fluidState.moleFraction(gasPhaseIdx, NAPLIdx)*phiAWC)
+ (fluidState.moleFraction(gasPhaseIdx, NAPLIdx) * mu[NAPLIdx])
/ (fluidState.moleFraction(gasPhaseIdx, NAPLIdx) + xAW*phiCAW);
return muResult;
}
@@ -281,28 +281,28 @@ public:
{
Scalar diffCont;
if (phaseIdx==gPhaseIdx) {
if (phaseIdx==gasPhaseIdx) {
Scalar diffAC = Opm::BinaryCoeff::Air_Xylene::gasDiffCoeff(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
Scalar diffWC = Opm::BinaryCoeff::H2O_Xylene::gasDiffCoeff(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
Scalar diffAW = Opm::BinaryCoeff::H2O_Air::gasDiffCoeff(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
const Scalar xga = fluidState.moleFraction(gPhaseIdx, airIdx);
const Scalar xgw = fluidState.moleFraction(gPhaseIdx, H2OIdx);
const Scalar xgc = fluidState.moleFraction(gPhaseIdx, NAPLIdx);
const Scalar xga = fluidState.moleFraction(gasPhaseIdx, airIdx);
const Scalar xgw = fluidState.moleFraction(gasPhaseIdx, H2OIdx);
const Scalar xgc = fluidState.moleFraction(gasPhaseIdx, NAPLIdx);
if (compIdx==NAPLIdx) return (1.- xgw)/(xga/diffAW + xgc/diffWC);
else if (compIdx==H2OIdx) return (1.- xgc)/(xgw/diffWC + xga/diffAC);
else if (compIdx==airIdx) OPM_THROW(std::logic_error,
"Diffusivity of air in the gas phase "
"is constraint by sum of diffusive fluxes = 0 !\n");
} else if (phaseIdx==wPhaseIdx){
} else if (phaseIdx==waterPhaseIdx){
Scalar diffACl = 1.e-9; // BinaryCoeff::Air_Xylene::liquidDiffCoeff(temperature, pressure);
Scalar diffWCl = 1.e-9; // BinaryCoeff::H2O_Xylene::liquidDiffCoeff(temperature, pressure);
Scalar diffAWl = 1.e-9; // BinaryCoeff::H2O_Air::liquidDiffCoeff(temperature, pressure);
Scalar xwa = fluidState.moleFraction(wPhaseIdx, airIdx);
Scalar xww = fluidState.moleFraction(wPhaseIdx, H2OIdx);
Scalar xwc = fluidState.moleFraction(wPhaseIdx, NAPLIdx);
Scalar xwa = fluidState.moleFraction(waterPhaseIdx, airIdx);
Scalar xww = fluidState.moleFraction(waterPhaseIdx, H2OIdx);
Scalar xwc = fluidState.moleFraction(waterPhaseIdx, NAPLIdx);
switch (compIdx) {
case NAPLIdx:
@@ -316,7 +316,7 @@ public:
"Diffusivity of water in the water phase "
"is constraint by sum of diffusive fluxes = 0 !\n");
};
} else if (phaseIdx==nPhaseIdx) {
} else if (phaseIdx==naplPhaseIdx) {
OPM_THROW(std::logic_error,
"Diffusion coefficients of "
@@ -338,7 +338,7 @@ public:
Scalar T = fluidState.temperature(phaseIdx);
Scalar p = fluidState.pressure(phaseIdx);
if (phaseIdx == wPhaseIdx) {
if (phaseIdx == waterPhaseIdx) {
if (compIdx == H2OIdx)
return H2O::vaporPressure(T)/p;
else if (compIdx == airIdx)
@@ -352,7 +352,7 @@ public:
// component to the NAPL phase is much higher than for the
// other components, i.e. the fugacity cofficient is much
// smaller.
if (phaseIdx == nPhaseIdx) {
if (phaseIdx == naplPhaseIdx) {
Scalar phiNapl = NAPL::vaporPressure(T)/p;
if (compIdx == NAPLIdx)
return phiNapl;
@@ -364,7 +364,7 @@ public:
// for the gas phase, assume an ideal gas when it comes to
// fugacity (-> fugacity == partial pressure)
assert(phaseIdx == gPhaseIdx);
assert(phaseIdx == gasPhaseIdx);
return 1.0;
}
@@ -374,13 +374,13 @@ public:
const ParameterCache &paramCache,
int phaseIdx)
{
if (phaseIdx == wPhaseIdx) {
if (phaseIdx == waterPhaseIdx) {
return H2O::liquidEnthalpy(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
}
else if (phaseIdx == nPhaseIdx) {
else if (phaseIdx == naplPhaseIdx) {
return NAPL::liquidEnthalpy(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
}
else if (phaseIdx == gPhaseIdx) { // gas phase enthalpy depends strongly on composition
else if (phaseIdx == gasPhaseIdx) { // gas phase enthalpy depends strongly on composition
Scalar hgc = NAPL::gasEnthalpy(fluidState.temperature(phaseIdx),
fluidState.pressure(phaseIdx));
Scalar hgw = H2O::gasEnthalpy(fluidState.temperature(phaseIdx),
@@ -388,9 +388,9 @@ public:
Scalar hga = Air::gasEnthalpy(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)); // pressure is only a dummy here (not dependent on pressure, just temperature)
Scalar result = 0;
result += hgw * fluidState.massFraction(gPhaseIdx, H2OIdx);
result += hga * fluidState.massFraction(gPhaseIdx, airIdx);
result += hgc * fluidState.massFraction(gPhaseIdx, NAPLIdx);
result += hgw * fluidState.massFraction(gasPhaseIdx, H2OIdx);
result += hga * fluidState.massFraction(gasPhaseIdx, airIdx);
result += hgc * fluidState.massFraction(gasPhaseIdx, NAPLIdx);
return result;
}

44
opm/material/fluidsystems/H2ON2FluidSystem.hpp
View File

@@ -73,9 +73,9 @@ public:
static const int numPhases = 2;
//! Index of the liquid phase
static const int lPhaseIdx = 0;
static const int liquidPhaseIdx = 0;
//! Index of the gas phase
static const int gPhaseIdx = 1;
static const int gasPhaseIdx = 1;
//! \copydoc BaseFluidSystem::phaseName
static const char *phaseName(int phaseIdx)
@@ -93,7 +93,7 @@ public:
static bool isLiquid(int phaseIdx)
{
//assert(0 <= phaseIdx && phaseIdx < numPhases);
return phaseIdx != gPhaseIdx;
return phaseIdx != gasPhaseIdx;
}
//! \copydoc BaseFluidSystem::isCompressible
@@ -102,7 +102,7 @@ public:
//assert(0 <= phaseIdx && phaseIdx < numPhases);
// gases are always compressible
return
(phaseIdx == gPhaseIdx)
(phaseIdx == gasPhaseIdx)
? true
:H2O::liquidIsCompressible();// the water component decides for the liquid phase...
}
@@ -113,7 +113,7 @@ public:
//assert(0 <= phaseIdx && phaseIdx < numPhases);
return
(phaseIdx == gPhaseIdx)
(phaseIdx == gasPhaseIdx)
? H2O::gasIsIdeal() && N2::gasIsIdeal() // let the components decide
: false; // not a gas
}
@@ -278,7 +278,7 @@ public:
Valgrind::CheckDefined(sumMoleFrac);
// liquid phase
if (phaseIdx == lPhaseIdx) {
if (phaseIdx == liquidPhaseIdx) {
if (!useComplexRelations)
// assume pure water
return H2O::liquidDensity(T, p);
@@ -292,9 +292,9 @@ public:
// water molecule in the liquid
return
clH2O
* (H2O::molarMass()*fluidState.moleFraction(lPhaseIdx, H2OIdx)
* (H2O::molarMass()*fluidState.moleFraction(liquidPhaseIdx, H2OIdx)
+
N2::molarMass()*fluidState.moleFraction(lPhaseIdx, N2Idx))
N2::molarMass()*fluidState.moleFraction(liquidPhaseIdx, N2Idx))
/ sumMoleFrac;
}
}
@@ -304,13 +304,13 @@ public:
// for the gas phase assume an ideal gas
return
IdealGas::molarDensity(T, p)
* fluidState.averageMolarMass(gPhaseIdx)
* fluidState.averageMolarMass(gasPhaseIdx)
/ std::max(1e-5, sumMoleFrac);
// assume ideal mixture: steam and nitrogen don't "see" each
// other
Scalar rho_gH2O = H2O::gasDensity(T, p*fluidState.moleFraction(gPhaseIdx, H2OIdx));
Scalar rho_gN2 = N2::gasDensity(T, p*fluidState.moleFraction(gPhaseIdx, N2Idx));
Scalar rho_gH2O = H2O::gasDensity(T, p*fluidState.moleFraction(gasPhaseIdx, H2OIdx));
Scalar rho_gN2 = N2::gasDensity(T, p*fluidState.moleFraction(gasPhaseIdx, N2Idx));
return (rho_gH2O + rho_gN2) / std::max(1e-5, sumMoleFrac);
}
@@ -326,7 +326,7 @@ public:
Scalar p = fluidState.pressure(phaseIdx);
// liquid phase
if (phaseIdx == lPhaseIdx) {
if (phaseIdx == liquidPhaseIdx) {
// assume pure water for the liquid phase
return H2O::liquidViscosity(T, p);
}
@@ -384,7 +384,7 @@ public:
Scalar p = fluidState.pressure(phaseIdx);
// liquid phase
if (phaseIdx == lPhaseIdx) {
if (phaseIdx == liquidPhaseIdx) {
if (compIdx == H2OIdx)
return H2O::vaporPressure(T)/p;
return Opm::BinaryCoeff::H2O_N2::henry(T)/p;
@@ -407,11 +407,11 @@ public:
Scalar p = fluidState.pressure(phaseIdx);
// liquid phase
if (phaseIdx == lPhaseIdx)
if (phaseIdx == liquidPhaseIdx)
return BinaryCoeff::H2O_N2::liquidDiffCoeff(T, p);
// gas phase
assert(phaseIdx == gPhaseIdx);
assert(phaseIdx == gasPhaseIdx);
return BinaryCoeff::H2O_N2::gasDiffCoeff(T, p);
}
@@ -427,7 +427,7 @@ public:
Valgrind::CheckDefined(p);
// liquid phase
if (phaseIdx == lPhaseIdx) {
if (phaseIdx == liquidPhaseIdx) {
// TODO: correct way to deal with the solutes???
return H2O::liquidEnthalpy(T, p);
}
@@ -438,10 +438,10 @@ public:
// that the total specific enthalpy is the sum of the
// "partial specific enthalpies" of the components.
Scalar hH2O =
fluidState.massFraction(gPhaseIdx, H2OIdx)
fluidState.massFraction(gasPhaseIdx, H2OIdx)
* H2O::gasEnthalpy(T, p);
Scalar hN2 =
fluidState.massFraction(gPhaseIdx, N2Idx)
fluidState.massFraction(gasPhaseIdx, N2Idx)
* N2::gasEnthalpy(T, p);
return hH2O + hN2;
}
@@ -458,7 +458,7 @@ public:
Scalar temperature = fluidState.temperature(phaseIdx) ;
Scalar pressure = fluidState.pressure(phaseIdx);
if (phaseIdx == lPhaseIdx) { // liquid phase
if (phaseIdx == liquidPhaseIdx) { // liquid phase
return H2O::liquidThermalConductivity(temperature, pressure);
}
else { // gas phase
@@ -489,7 +489,7 @@ public:
const ParameterCache &paramCache,
int phaseIdx)
{
if (phaseIdx == lPhaseIdx) {
if (phaseIdx == liquidPhaseIdx) {
return H2O::liquidHeatCapacity(fluidState.temperature(phaseIdx),
fluidState.pressure(phaseIdx));
}
@@ -526,8 +526,8 @@ public:
// mangle both components together
return
c_pH2O*fluidState.massFraction(gPhaseIdx, H2OIdx)
+ c_pN2*fluidState.massFraction(gPhaseIdx, N2Idx);
c_pH2O*fluidState.massFraction(gasPhaseIdx, H2OIdx)
+ c_pN2*fluidState.massFraction(gasPhaseIdx, N2Idx);
}
};

24
opm/material/fluidsystems/H2ON2LiquidPhaseFluidSystem.hpp
View File

@@ -74,12 +74,12 @@ public:
static const int numPhases = 1;
//! Index of the liquid phase
static const int lPhaseIdx = 0;
static const int liquidPhaseIdx = 0;
//! \copydoc BaseFluidSystem::phaseName
static const char *phaseName(int phaseIdx)
{
assert(phaseIdx == lPhaseIdx);
assert(phaseIdx == liquidPhaseIdx);
return "l";
}
@@ -87,7 +87,7 @@ public:
//! \copydoc BaseFluidSystem::isLiquid
static bool isLiquid(int phaseIdx)
{
//assert(phaseIdx == lPhaseIdx);
//assert(phaseIdx == liquidPhaseIdx);
return true; //only water phase present
}
@@ -259,7 +259,7 @@ public:
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
sumMoleFrac += fluidState.moleFraction(phaseIdx, compIdx);
assert(phaseIdx == lPhaseIdx);
assert(phaseIdx == liquidPhaseIdx);
if (!useComplexRelations)
// assume pure water
@@ -274,9 +274,9 @@ public:
// water molecule in the liquid
return
clH2O
* (H2O::molarMass()*fluidState.moleFraction(lPhaseIdx, H2OIdx)
* (H2O::molarMass()*fluidState.moleFraction(liquidPhaseIdx, H2OIdx)
+
N2::molarMass()*fluidState.moleFraction(lPhaseIdx, N2Idx))
N2::molarMass()*fluidState.moleFraction(liquidPhaseIdx, N2Idx))
/ sumMoleFrac;
}
}
@@ -287,7 +287,7 @@ public:
const ParameterCache &paramCache,
int phaseIdx)
{
assert(phaseIdx == lPhaseIdx);
assert(phaseIdx == liquidPhaseIdx);
Scalar T = fluidState.temperature(phaseIdx);
Scalar p = fluidState.pressure(phaseIdx);
@@ -303,7 +303,7 @@ public:
int phaseIdx,
int compIdx)
{
assert(phaseIdx == lPhaseIdx);
assert(phaseIdx == liquidPhaseIdx);
assert(0 <= compIdx && compIdx < numComponents);
Scalar T = fluidState.temperature(phaseIdx);
@@ -322,7 +322,7 @@ public:
int compIdx)
{
assert(phaseIdx == lPhaseIdx);
assert(phaseIdx == liquidPhaseIdx);
Scalar T = fluidState.temperature(phaseIdx);
Scalar p = fluidState.pressure(phaseIdx);
@@ -336,7 +336,7 @@ public:
const ParameterCache &paramCache,
int phaseIdx)
{
assert (phaseIdx == lPhaseIdx);
assert (phaseIdx == liquidPhaseIdx);
Scalar T = fluidState.temperature(phaseIdx);
Scalar p = fluidState.pressure(phaseIdx);
@@ -353,7 +353,7 @@ public:
const ParameterCache &paramCache,
const int phaseIdx)
{
assert(phaseIdx == lPhaseIdx);
assert(phaseIdx == liquidPhaseIdx);
if(useComplexRelations){
Scalar temperature = fluidState.temperature(phaseIdx) ;
@@ -370,7 +370,7 @@ public:
const ParameterCache &paramCache,
int phaseIdx)
{
assert (phaseIdx == lPhaseIdx);
assert (phaseIdx == liquidPhaseIdx);
return H2O::liquidHeatCapacity(fluidState.temperature(phaseIdx),
fluidState.pressure(phaseIdx));

26
opm/material/fluidsystems/Spe5FluidSystem.hpp
View File

@@ -70,11 +70,11 @@ public:
static const int numPhases = 3;
//! Index of the gas phase
static const int gPhaseIdx = 0;
static const int gasPhaseIdx = 0;
//! Index of the water phase
static const int wPhaseIdx = 1;
static const int waterPhaseIdx = 1;
//! Index of the oil phase
static const int oPhaseIdx = 2;
static const int oilPhaseIdx = 2;
//! The component for pure water to be used
typedef Opm::H2O<Scalar> H2O;
@@ -96,7 +96,7 @@ public:
static bool isLiquid(int phaseIdx)
{
//assert(0 <= phaseIdx && phaseIdx < numPhases);
return phaseIdx != gPhaseIdx;
return phaseIdx != gasPhaseIdx;
}
/*!
@@ -123,7 +123,7 @@ public:
// always use the reference oil for the fugacity coefficents,
// so they cannot be dependent on composition and they the
// phases thus always an ideal mixture
return phaseIdx == wPhaseIdx;
return phaseIdx == waterPhaseIdx;
}
/****************************************
@@ -304,7 +304,7 @@ public:
Scalar minP = 1e4,
Scalar maxP = 100e6)
{
Opm::PengRobinsonParamsMixture<Scalar, ThisType, gPhaseIdx, /*useSpe5=*/true> prParams;
Opm::PengRobinsonParamsMixture<Scalar, ThisType, gasPhaseIdx, /*useSpe5=*/true> prParams;
// find envelopes of the 'a' and 'b' parameters for the range
// minT <= T <= maxT and minP <= p <= maxP. For
@@ -370,16 +370,16 @@ public:
{
assert(0 <= phaseIdx && phaseIdx <= numPhases);
if (phaseIdx == gPhaseIdx) {
if (phaseIdx == gasPhaseIdx) {
// given by SPE-5 in table on page 64. we use a constant
// viscosity, though...
return 0.0170e-2 * 0.1;
}
else if (phaseIdx == wPhaseIdx)
else if (phaseIdx == waterPhaseIdx)
// given by SPE-5: 0.7 centi-Poise = 0.0007 Pa s
return 0.7e-2 * 0.1;
else {
assert(phaseIdx == oPhaseIdx);
assert(phaseIdx == oilPhaseIdx);
// given by SPE-5 in table on page 64. we use a constant
// viscosity, though...
return 0.208e-2 * 0.1;
@@ -396,16 +396,16 @@ public:
assert(0 <= phaseIdx && phaseIdx <= numPhases);
assert(0 <= compIdx && compIdx <= numComponents);
if (phaseIdx == oPhaseIdx || phaseIdx == gPhaseIdx)
if (phaseIdx == oilPhaseIdx || phaseIdx == gasPhaseIdx)
return PengRobinsonMixture::computeFugacityCoefficient(fluidState,
paramCache,
phaseIdx,
compIdx);
else {
assert(phaseIdx == wPhaseIdx);
assert(phaseIdx == waterPhaseIdx);
return
henryCoeffWater_(compIdx, fluidState.temperature(wPhaseIdx))
/ fluidState.pressure(wPhaseIdx);
henryCoeffWater_(compIdx, fluidState.temperature(waterPhaseIdx))
/ fluidState.pressure(waterPhaseIdx);
}
}

58
opm/material/fluidsystems/Spe5ParameterCache.hpp
View File

@@ -13,7 +13,7 @@
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
You should have received a copy o2f the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
/*!
@@ -48,15 +48,15 @@ class Spe5ParameterCache
enum { numPhases = FluidSystem::numPhases };
enum { wPhaseIdx = FluidSystem::wPhaseIdx };
enum { oPhaseIdx = FluidSystem::oPhaseIdx };
enum { gPhaseIdx = FluidSystem::gPhaseIdx };
enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
public:
//! The cached parameters for the oil phase
typedef Opm::PengRobinsonParamsMixture<Scalar, FluidSystem, oPhaseIdx, /*useSpe5=*/true> OilPhaseParams;
typedef Opm::PengRobinsonParamsMixture<Scalar, FluidSystem, oilPhaseIdx, /*useSpe5=*/true> OilPhaseParams;
//! The cached parameters for the gas phase
typedef Opm::PengRobinsonParamsMixture<Scalar, FluidSystem, gPhaseIdx, /*useSpe5=*/true> GasPhaseParams;
typedef Opm::PengRobinsonParamsMixture<Scalar, FluidSystem, gasPhaseIdx, /*useSpe5=*/true> GasPhaseParams;
Spe5ParameterCache()
{
@@ -89,9 +89,9 @@ public:
int phaseIdx,
int compIdx)
{
if (phaseIdx == oPhaseIdx)
if (phaseIdx == oilPhaseIdx)
oilPhaseParams_.updateSingleMoleFraction(fluidState, compIdx);
else if (phaseIdx == gPhaseIdx)
else if (phaseIdx == gasPhaseIdx)
gasPhaseParams_.updateSingleMoleFraction(fluidState, compIdx);
// update the phase's molar volume
@@ -107,8 +107,8 @@ public:
{
switch (phaseIdx)
{
case oPhaseIdx: return oilPhaseParams_.a();
case gPhaseIdx: return gasPhaseParams_.a();
case oilPhaseIdx: return oilPhaseParams_.a();
case gasPhaseIdx: return gasPhaseParams_.a();
default:
OPM_THROW(std::logic_error,
"The a() parameter is only defined for "
@@ -125,8 +125,8 @@ public:
{
switch (phaseIdx)
{
case oPhaseIdx: return oilPhaseParams_.b();
case gPhaseIdx: return gasPhaseParams_.b();
case oilPhaseIdx: return oilPhaseParams_.b();
case gasPhaseIdx: return gasPhaseParams_.b();
default:
OPM_THROW(std::logic_error,
"The b() parameter is only defined for "
@@ -146,8 +146,8 @@ public:
{
switch (phaseIdx)
{
case oPhaseIdx: return oilPhaseParams_.pureParams(compIdx).a();
case gPhaseIdx: return gasPhaseParams_.pureParams(compIdx).a();
case oilPhaseIdx: return oilPhaseParams_.pureParams(compIdx).a();
case gasPhaseIdx: return gasPhaseParams_.pureParams(compIdx).a();
default:
OPM_THROW(std::logic_error,
"The a() parameter is only defined for "
@@ -166,8 +166,8 @@ public:
{
switch (phaseIdx)
{
case oPhaseIdx: return oilPhaseParams_.pureParams(compIdx).b();
case gPhaseIdx: return gasPhaseParams_.pureParams(compIdx).b();
case oilPhaseIdx: return oilPhaseParams_.pureParams(compIdx).b();
case gasPhaseIdx: return gasPhaseParams_.pureParams(compIdx).b();
default:
OPM_THROW(std::logic_error,
"The b() parameter is only defined for "
@@ -242,9 +242,9 @@ protected:
switch (phaseIdx)
{
case oPhaseIdx: oilPhaseParams_.updatePure(T, p); break;
case gPhaseIdx: gasPhaseParams_.updatePure(T, p); break;
//case wPhaseIdx: waterPhaseParams_.updatePure(phaseIdx, temperature, pressure);break;
case oilPhaseIdx: oilPhaseParams_.updatePure(T, p); break;
case gasPhaseIdx: gasPhaseParams_.updatePure(T, p); break;
//case waterPhaseIdx: waterPhaseParams_.updatePure(phaseIdx, temperature, pressure);break;
}
}
@@ -261,13 +261,13 @@ protected:
Valgrind::CheckDefined(fluidState.averageMolarMass(phaseIdx));
switch (phaseIdx)
{
case oPhaseIdx:
case oilPhaseIdx:
oilPhaseParams_.updateMix(fluidState);
break;
case gPhaseIdx:
case gasPhaseIdx:
gasPhaseParams_.updateMix(fluidState);
break;
case wPhaseIdx:
case waterPhaseIdx:
break;
}
}
@@ -281,44 +281,44 @@ protected:
// calculate molar volume of the phase (we will need this for the
// fugacity coefficients and the density anyway)
switch (phaseIdx) {
case gPhaseIdx: {
case gasPhaseIdx: {
// calculate molar volumes for the given composition. although
// this isn't a Peng-Robinson parameter strictly speaking, the
// molar volume appears in basically every quantity the fluid
// system can get queried, so it is okay to calculate it
// here...
Vm_[gPhaseIdx] =
Vm_[gasPhaseIdx] =
PengRobinson::computeMolarVolume(fluidState,
*this,
phaseIdx,
/*isGasPhase=*/true);
}
case oPhaseIdx: {
case oilPhaseIdx: {
// calculate molar volumes for the given composition. although
// this isn't a Peng-Robinson parameter strictly speaking, the
// molar volume appears in basically every quantity the fluid
// system can get queried, so it is okay to calculate it
// here...
Vm_[oPhaseIdx] =
Vm_[oilPhaseIdx] =
PengRobinson::computeMolarVolume(fluidState,
*this,
phaseIdx,
/*isGasPhase=*/false);
}
case wPhaseIdx: {
case waterPhaseIdx: {
// Density of water in the stock tank (i.e. atmospheric
// pressure) is specified as 62.4 lb/ft^3 by the SPE-5
// paper. Also 1 lb = 0.4535923 and 1 ft = 0.3048 m.
const Scalar stockTankWaterDensity = 62.4 * 0.45359237 / 0.028316847;
// Water compressibility is specified as 3.3e-6 per psi
// overpressure, where 1 psi = 6894.7573 Pa
Scalar overPressure = fluidState.pressure(wPhaseIdx) - 1.013e5; // [Pa]
Scalar overPressure = fluidState.pressure(waterPhaseIdx) - 1.013e5; // [Pa]
Scalar waterDensity =
stockTankWaterDensity * (1 + 3.3e-6*overPressure/6894.7573);
// convert water density [kg/m^3] to molar volume [m^3/mol]
Vm_[wPhaseIdx] = fluidState.averageMolarMass(wPhaseIdx)/waterDensity;
Vm_[waterPhaseIdx] = fluidState.averageMolarMass(waterPhaseIdx)/waterDensity;
};
};
}