Commit ebd395cb authored by Katharina Heck's avatar Katharina Heck Committed by Timo Koch
Browse files

[material][fluidsystems] introduce molar density in all fluid systems

    Co-authored-by: Simon Scholz simon.scholz@iws.uni-stuttgart.de
    Co-authored-by: Beatrix Becker beatrix.becker@iws.uni-stuttgart.de
parent a53dd95e
......@@ -188,6 +188,36 @@ public:
fluidState.pressure(phaseIdx));
}
/*!
* \brief The molar density \f$\rho_{mol,\alpha}\f$
* of a fluid phase \f$\alpha\f$ in \f$\mathrm{[mol/m^3]}\f$
*
* The molar density is defined by the
* mass density \f$\rho_\alpha\f$ and the component molar mass \f$M_\alpha\f$:
*
* \f[\rho_{mol,\alpha} = \frac{\rho_\alpha}{M_\alpha} \;.\f]
*/
static Scalar molarDensity(Scalar temperature, Scalar pressure)
{ return Component::gasMolarDensity(temperature, pressure); }
/*!
* \brief The molar density \f$\rho_{mol,\alpha}\f$
* of a fluid phase \f$\alpha\f$ in \f$\mathrm{[mol/m^3]}\f$
*
* The molar density is defined by the
* mass density \f$\rho_\alpha\f$ and the component molar mass \f$M_\alpha\f$:
*
* \f[\rho_{mol,\alpha} = \frac{\rho_\alpha}{M_\alpha} \;.\f]
*/
using Base::molarDensity;
template <class FluidState>
static Scalar molarDensity(const FluidState &fluidState,
const int phaseIdx)
{
return molarDensity(fluidState.temperature(phaseIdx),
fluidState.pressure(phaseIdx));
}
/*!
* \brief The pressure \f$\mathrm{[Pa]}\f$ of the component at a given density and temperature.
* \param temperature The given temperature \f$\mathrm{[K]}\f$
......
......@@ -191,6 +191,29 @@ public:
fluidState.pressure(phaseIdx));
}
using Base::molarDensity;
/*!
* \brief The molar density \f$\rho_{mol,\alpha}\f$
* of a fluid phase \f$\alpha\f$ in \f$\mathrm{[mol/m^3]}\f$
*
* The molar density is defined by the
* mass density \f$\rho_\alpha\f$ and the main component molar mass \f$M_\alpha\f$:
*
* \f[\rho_{mol,\alpha} = \frac{\rho_\alpha}{M_\alpha} \;.\f]
*/
template <class FluidState>
static Scalar molarDensity(const FluidState &fluidState, const int phaseIdx)
{
return molarDensity(fluidState.temperature(phaseIdx),
fluidState.pressure(phaseIdx));
}
/*!
* \brief The density \f$\mathrm{[kg/m^3]}\f$ of the component at a given pressure and temperature.
*/
static Scalar molarDensity(Scalar temperature, Scalar pressure)
{ return Component::liquidMolarDensity(temperature, pressure); }
/*!
* \brief The pressure \f$\mathrm{[Pa]}\f$ of the component at a given density and temperature.
*/
......
......@@ -279,21 +279,43 @@ public:
{
assert(0 <= phaseIdx && phaseIdx < numPhases);
Scalar T = fluidState.temperature(phaseIdx);
Scalar p = fluidState.pressure(phaseIdx);
Scalar temperature = fluidState.temperature(phaseIdx);
Scalar pressure = fluidState.pressure(phaseIdx);
// liquid phase
if (phaseIdx == liquidPhaseIdx) {
return Component::liquidDensity(T, p);
return Component::liquidDensity(temperature, pressure);
}
else if (phaseIdx == gasPhaseIdx)// gas phase
{
return Component::gasDensity(T, p);
return Component::gasDensity(temperature, pressure);
}
else DUNE_THROW(Dune::NotImplemented,
"wrong index");
}
using Base::molarDensity;
/*!
* \brief The molar density \f$\rho_{mol,\alpha}\f$
* of a fluid phase \f$\alpha\f$ in \f$\mathrm{[mol/m^3]}\f$
*
* The molar density is defined by the
* mass density \f$\rho_\alpha\f$ and the molar mass \f$M_\alpha\f$:
*
* \f[\rho_{mol,\alpha} = \frac{\rho_\alpha}{M_\alpha} \;.\f]
*/
template <class FluidState>
static Scalar molarDensity(const FluidState &fluidState, int phaseIdx)
{
Scalar temperature = fluidState.temperature(phaseIdx);
Scalar pressure = fluidState.temperature(phaseIdx);
if (phaseIdx == liquidPhaseIdx)
return Component::liquidMolarDensity(temperature, pressure);
return Component::gasMolarDensity(temperature, pressure);
}
/*!
* \brief Calculate the dynamic viscosity of a fluid phase [Pa*s]
*
......@@ -307,16 +329,16 @@ public:
{
assert(0 <= phaseIdx && phaseIdx < numPhases);
Scalar T = fluidState.temperature(phaseIdx);
Scalar p = fluidState.pressure(phaseIdx);
Scalar temperature = fluidState.temperature(phaseIdx);
Scalar pressure = fluidState.pressure(phaseIdx);
// liquid phase
if (phaseIdx == liquidPhaseIdx) {
return Component::liquidViscosity(T, p);
return Component::liquidViscosity(temperature, pressure);
}
else if (phaseIdx == gasPhaseIdx) // gas phase
{
return Component::gasViscosity(T, p) ;
return Component::gasViscosity(temperature, pressure) ;
}
else DUNE_THROW(Dune::NotImplemented,
"wrong index");
......@@ -374,12 +396,14 @@ public:
assert(0 <= phaseIdx && phaseIdx < numPhases);
assert(0 <= compIdx && compIdx < numComponents);
Scalar T = fluidState.temperature(phaseIdx);
Scalar p = fluidState.pressure(phaseIdx);
Scalar temperature = fluidState.temperature(phaseIdx);
Scalar pressure = fluidState.pressure(phaseIdx);
// liquid phase
if (phaseIdx == liquidPhaseIdx)
return Component::vaporPressure(T)/p;
{
return Component::vaporPressure(temperature)/pressure;
}
// for the gas phase, assume an ideal gas when it comes to
// fugacity (-> fugacity == partial pressure)
......@@ -443,12 +467,12 @@ public:
// liquid phase
if (phaseIdx == liquidPhaseIdx) {
return Component::liquidEnthalpy(fluidState.temperature(phaseIdx),
fluidState.pressure(phaseIdx));
fluidState.pressure(phaseIdx));
}
else if (phaseIdx == gasPhaseIdx) // gas phase
{
return Component::gasEnthalpy(fluidState.temperature(phaseIdx),
fluidState.pressure(phaseIdx));
fluidState.pressure(phaseIdx));
}
else DUNE_THROW(Dune::NotImplemented,
"wrong index");
......@@ -471,12 +495,12 @@ public:
// liquid phase
if (phaseIdx == liquidPhaseIdx) {
return Component::liquidThermalConductivity(fluidState.temperature(phaseIdx),
fluidState.pressure(phaseIdx)); //0.68 ;
fluidState.pressure(phaseIdx)); //0.68 ;
}
else if (phaseIdx == gasPhaseIdx) // gas phase
{
return Component::gasThermalConductivity(fluidState.temperature(phaseIdx),
fluidState.pressure(phaseIdx)); //0.0248;
fluidState.pressure(phaseIdx)); //0.0248;
}
else DUNE_THROW(Dune::NotImplemented,
"wrong index");
......@@ -497,13 +521,13 @@ public:
assert(0 <= phaseIdx && phaseIdx < numPhases);
// liquid phase
if (phaseIdx == liquidPhaseIdx) {
return Component::liquidHeatCapacity(fluidState.temperature(phaseIdx),
fluidState.pressure(phaseIdx));//4.217e3 ;
return Component::liquidHeatCapacity(fluidState.temperature(phaseIdx),
fluidState.pressure(phaseIdx));//4.217e3 ;
}
else if (phaseIdx == gasPhaseIdx) // gas phase
{
return Component::gasHeatCapacity(fluidState.temperature(phaseIdx),
fluidState.pressure(phaseIdx));//2.029e3;
fluidState.pressure(phaseIdx));//2.029e3;
}
else DUNE_THROW(Dune::NotImplemented, "wrong index");
}
......
......@@ -294,6 +294,26 @@ public:
return Fluid1::density(temperature, pressure);
}
using Base::molarDensity;
/*!
* \brief The molar density \f$\rho_{mol,\alpha}\f$
* of a fluid phase \f$\alpha\f$ in \f$\mathrm{[mol/m^3]}\f$
*
* The molar density is defined by the
* mass density \f$\rho_\alpha\f$ and the component molar mass \f$M_\alpha\f$:
*
* \f[\rho_{mol,\alpha} = \frac{\rho_\alpha}{M_\alpha} \;.\f]
*/
template <class FluidState>
static Scalar molarDensity(const FluidState &fluidState, int phaseIdx)
{
Scalar temperature = fluidState.temperature(phaseIdx);
Scalar pressure = fluidState.pressure(phaseIdx);
if (phaseIdx == phase0Idx)
return Fluid0::molarDensity(temperature, pressure);
return Fluid1::molarDensity(temperature, pressure);
}
using Base::viscosity;
/*!
* \brief Return the viscosity of a phase \f$\mathrm{[Pa*s]}\f$.
......
......@@ -365,6 +365,34 @@ public:
}
}
using Base::molarDensity;
/*!
* \brief The molar density \f$\rho_{mol,\alpha}\f$
* of a fluid phase \f$\alpha\f$ in \f$\mathrm{[mol/m^3]}\f$
*
* The molar density is defined by the
* mass density \f$\rho_\alpha\f$ and the component molar mass \f$M_\alpha\f$:
*
* \f[\rho_{mol,\alpha} = \frac{\rho_\alpha}{M_\alpha} \;.\f]
*/
template <class FluidState>
static Scalar molarDensity(const FluidState &fluidState,
int phaseIdx)
{
assert(0 <= phaseIdx && phaseIdx < numPhases);
Scalar temperature = fluidState.temperature(phaseIdx);
Scalar pressure = fluidState.pressure(phaseIdx);
switch(phaseIdx)
{
case wPhaseIdx: return WettingFluid::molarDensity(temperature, pressure);
case nPhaseIdx: return NonwettingFluid::molarDensity(temperature, pressure);
case gPhaseIdx: return Gas::molarDensity(temperature, pressure);
default: DUNE_THROW(Dune::InvalidStateException, "Invalid phase index");
}
}
using Base::viscosity;
/*!
* \brief Return the viscosity of a phase \f$\mathrm{[Pa*s]}\f$.
......
......@@ -127,7 +127,21 @@ public:
}
/*!
* \brief Calculate the fugacity coefficient \f$\mathrm{[-]}\f$ of an individual
* \brief Calculate the molar density \f$\mathrm{[mol/m^3]}\f$ of a fluid phase
* \param fluidState The fluid state
* \param paramCache mutable parameters
* \param phaseIdx Index of the fluid phase
*/
template <class FluidState>
static Scalar molarDensity(const FluidState &fluidState,
const ParameterCache &paramCache,
int phaseIdx)
{
return Implementation::molarDensity(fluidState, phaseIdx);
}
/*!
* \brief Calculate the fugacity coefficient \f$\mathrm{[Pa]}\f$ of an individual
* component in a fluid phase
*
* The fugacity coefficient \f$\mathrm{\phi^\kappa_\alpha}\f$ is connected to the
......
......@@ -124,6 +124,17 @@ public:
return phaseIdx != gasPhaseIdx;
}
/*!
* \brief Return whether a phase is gaseous
*
* \param phaseIdx The index of the fluid phase to consider
*/
static constexpr bool isGas(int phaseIdx)
{
assert(0 <= phaseIdx && phaseIdx < numPhases);
return phaseIdx == phase0Idx;
}
/*!
* \brief Returns true if and only if a fluid phase is assumed to
* be an ideal mixture.
......@@ -224,9 +235,9 @@ public:
* \param Temperature temperature of the liquid phase
* \param salinity salinity (salt mass fraction) of the liquid phase
*/
static Scalar vaporPressure(Scalar Temperature, Scalar salinity) //Vapor pressure dependence on Osmosis
static Scalar vaporPressure(Scalar temperature, Scalar salinity) //Vapor pressure dependence on Osmosis
{
return vaporPressure_(Temperature,salinity);
return Brine::vaporPressure(temperature, salinity);
}
......@@ -279,7 +290,7 @@ public:
}
using Base::density;
/*!
/*!
* \brief Given a phase's composition, temperature, pressure, and
* the partial pressures of all components, return its
* density \f$\mathrm{[kg/m^3]}\f$.
......@@ -301,20 +312,71 @@ public:
Scalar temperature = fluidState.temperature(phaseIdx);
Scalar pressure = fluidState.pressure(phaseIdx);
switch (phaseIdx) {
case liquidPhaseIdx:
if (phaseIdx == phase0Idx){
if (!useComplexRelations)
// assume pure brine
return Brine::liquidDensity(temperature,
pressure,
fluidState.massFraction(liquidPhaseIdx, NaClIdx));
case gasPhaseIdx:
return gasDensity_(temperature,
pressure,
fluidState.moleFraction(gasPhaseIdx, H2OIdx));
default:
DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
pressure,
fluidState.massFraction(phase0Idx, NaClIdx));
else
{
return Brine::liquidMolarDensity(temperature,
pressure,
fluidState.massFraction(phase0Idx, NaClIdx))
*(Brine::molarMass()*fluidState.moleFraction(liquidPhaseIdx, H2OIdx)
+ Brine::molarMass()*fluidState.moleFraction(liquidPhaseIdx, NaClIdx)
+ Air::molarMass()*fluidState.moleFraction(liquidPhaseIdx, AirIdx));
}
}
}
else if (phaseIdx == phase1Idx){
if (!useComplexRelations)
// for the gas phase assume an ideal gas
{
const Scalar averageMolarMass = fluidState.averageMolarMass(phase1Idx);
return IdealGas::density(averageMolarMass, temperature, pressure);
}
return
Brine::gasDensity(temperature, fluidState.partialPressure(phase1Idx, H2OIdx)) +
Air::gasDensity(temperature, fluidState.partialPressure(phase1Idx, AirIdx));
}
else
DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
}
using Base::molarDensity;
/*!
* \brief The molar density \f$\rho_{mol,\alpha}\f$
* of a fluid phase \f$\alpha\f$ in \f$\mathrm{[mol/m^3]}\f$
*
* The molar density for the complex relation is defined by the
* mass density \f$\rho_\alpha\f$ and the mean molar mass \f$\overline M_\alpha\f$:
*
* \f[\rho_{mol,\alpha} = \frac{\rho_\alpha}{\overline M_\alpha} \;.\f]
*/
template <class FluidState>
static Scalar molarDensity(const FluidState &fluidState, int phaseIdx)
{
const Scalar temperature = fluidState.temperature(phaseIdx);
const Scalar pressure = fluidState.pressure(phaseIdx);
if (phaseIdx == liquidPhaseIdx)
{
// assume pure water or that each gas molecule displaces exactly one
// molecule in the liquid.
Scalar salinity = fluidState.massFraction(liquidPhaseIdx, NaClIdx);
return Brine::liquidMolarDensity(temperature, pressure, salinity);
}
else if (phaseIdx == phase1Idx)
{
if (!useComplexRelations)
// for the gas phase assume an ideal gas
{ return IdealGas::molarDensity(temperature, pressure); }
return
Brine::gasMolarDensity(temperature, fluidState.partialPressure(phase1Idx, H2OIdx)) +
Air::gasMolarDensity(temperature, fluidState.partialPressure(phase1Idx, AirIdx));
}
else
DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
}
using Base::viscosity;
/*!
......@@ -343,7 +405,7 @@ public:
Scalar XNaCl = fluidState.massFraction(liquidPhaseIdx, NaClIdx);
result = Brine::liquidViscosity(temperature, pressure, XNaCl);
}
else if (phaseIdx == gasPhaseIdx)
else if (phaseIdx == phase1Idx)
{
result = Air::gasViscosity(temperature, pressure);
}
......@@ -392,12 +454,13 @@ public:
if (phaseIdx == gasPhaseIdx)
return 1.0;
else if (phaseIdx == liquidPhaseIdx)
else if (phaseIdx == phase0Idx)
{
if (compIdx == H2OIdx)
return Brine::vaporPressure(T)/p;
else if (compIdx == AirIdx)
return BinaryCoeff::H2O_Air::henry(T)/p;
if (compIdx == H2OIdx)
return H2O::vaporPressure(T)/p;
else if (compIdx == AirIdx)
return BinaryCoeff::H2O_Air::henry(T)/p;
else
return 1/p;
}
......@@ -625,17 +688,6 @@ public:
return Brine_Air::molalityNaCl(salinity);
}
private:
static Scalar gasDensity_(Scalar T, Scalar pg, Scalar xgH2O)
{
//Dalton' Law
const Scalar pH2O = xgH2O*pg;
const Scalar pAir = pg - pH2O;
const Scalar gasDensityAir = Air::gasDensity(T, pAir);
const Scalar gasDensityH2O = H2O::gasDensity(T, pH2O);
const Scalar gasDensity = gasDensityAir + gasDensityH2O;
return gasDensity;
}
};
} // end namespace FluidSystems
......
......@@ -287,24 +287,33 @@ public:
else {
assert(phaseIdx == gasPhaseIdx);
// use normalized composition for to calculate the density
// (the relations don't seem to take non-normalized
// compositions too well...)
using std::min;
using std::max;
Scalar xgBrine = min(1.0, max(0.0, fluidState.moleFraction(gasPhaseIdx, BrineIdx)));
Scalar xgCO2 = min(1.0, max(0.0, fluidState.moleFraction(gasPhaseIdx, CO2Idx)));
Scalar sumx = xgBrine + xgCO2;
xgBrine /= sumx;
xgCO2 /= sumx;
Scalar result = gasDensity_(temperature,
pressure,
xgBrine,
xgCO2);
Valgrind::CheckDefined(result);
return result;
//in the end it comes down to a simplification of just CO2
return CO2::gasDensity(temperature, pressure);
}
}
using Base::molarDensity;
/*!
* \brief The molar density \f$\rho_{mol,\alpha}\f$
* of a fluid phase \f$\alpha\f$ in \f$\mathrm{[mol/m^3]}\f$
*
* The molar density is defined by the
* mass density \f$\rho_\alpha\f$ and the mean molar mass \f$\overline M_\alpha\f$:
*
* \f[\rho_{mol,\alpha} = \frac{\rho_\alpha}{\overline M_\alpha} \;.\f]
*/
template <class FluidState>
static Scalar molarDensity(const FluidState &fluidState, int phaseIdx)
{
Scalar temperature = fluidState.temperature(phaseIdx);
Scalar pressure = fluidState.pressure(phaseIdx);
if (phaseIdx == liquidPhaseIdx)
{
return density(fluidState, phaseIdx)/fluidState.averageMolarMass(phaseIdx);
}
else
return CO2::gasMolarDensity(temperature,pressure);
}
using Base::viscosity;
......@@ -623,19 +632,6 @@ public:
}
private:
static Scalar gasDensity_(Scalar T,
Scalar pg,
Scalar xgH2O,
Scalar xgCO2)
{
Valgrind::CheckDefined(T);
Valgrind::CheckDefined(pg);
Valgrind::CheckDefined(xgH2O);
Valgrind::CheckDefined(xgCO2);
Scalar gasDensity = CO2::gasDensity(T, pg);
return gasDensity;
}
/***********************************************************************/
/* */
......
......@@ -89,9 +89,10 @@ public:
*/
static std::string phaseName(int phaseIdx)
{
switch (phaseIdx) {
case liquidPhaseIdx: return "liquid";
case gasPhaseIdx: return "gas";
switch (phaseIdx)
{
case liquidPhaseIdx: return "liquid";
case gasPhaseIdx: return "gas";
}
DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
}
......@@ -194,8 +195,8 @@ public:
{
switch (compIdx)
{
case H2OIdx: return H2O::name();
case AirIdx: return Air::name();
case H2OIdx: return H2O::name();
case AirIdx: return Air::name();
}
DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
}
......@@ -209,13 +210,12 @@ public:
{
switch (compIdx)
{
case H2OIdx: return H2O::molarMass();
case AirIdx: return Air::molarMass();
case H2OIdx: return H2O::molarMass();
case AirIdx: return Air::molarMass();
}
DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
}
/*!
* \brief Critical temperature of a component \f$\mathrm{[K]}\f$.
*
......@@ -351,7 +351,7 @@ public:
<< " entries).\n";
H2O::init(tempMin, tempMax, nTemp,
pressMin, pressMax, nPress);
pressMin, pressMax, nPress);
}
}
......@@ -378,12 +378,7 @@ public:
const Scalar T = fluidState.temperature(phaseIdx);
const Scalar p = fluidState.pressure(phaseIdx);
Scalar sumMoleFrac = 0;
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
sumMoleFrac += fluidState.moleFraction(phaseIdx, compIdx);
if (phaseIdx == liquidPhaseIdx)
if (phaseIdx == phase0Idx)
{
if (!useComplexRelations)
// assume pure water
......@@ -391,34 +386,63 @@ public:
else
{
// See: Eq. (7) in Class et al. (2002a)
const Scalar rholH2O = H2O::liquidDensity(T, p);
const Scalar clH2O = rholH2O/H2O::molarMass();
return
clH2O
* (H2O::molarMass()*fluidState.moleFraction(liquidPhaseIdx, H2OIdx)
+
Air::molarMass()*fluidState.moleFraction(liquidPhaseIdx, AirIdx))
/ sumMoleFrac;
// This assumes each gas molecule displaces exactly one