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

change the interface of the h2o_n2 fluid system to allow arbitrary equations of state

Browse Source
This commit is contained in:
Andreas Lauser
2010-03-01 16:42:39 +00:00
committed by Andreas Lauser
parent 29350f1f3a
commit 3c29ccfef1
7 changed files with 527 additions and 97 deletions
Download Patch File Download Diff File Expand all files Collapse all files

24
dumux/new_material/components/h2.hh
View File

@@ -109,9 +109,18 @@ public:
static Scalar density(Scalar temperature, Scalar pressure)
{
// Assume an ideal gas
return IdealGas::density(molarMass(), pressure, temperature);
return IdealGas::density(molarMass(), temperature, pressure);
}
/*
* \brief The pressure of gaseous N2 at a given density and temperature [Pa].
*/
static Scalar gasPressure(Scalar temperature, Scalar density)
{
// Assume an ideal gas
return IdealGas::pressure(temperature, density/molarMass());
}
/*!
* \brief Specific enthalpy [J/kg] of pure hydrogen gas.
*/
@@ -122,6 +131,19 @@ public:
return temperature * (cvHat + 1) * IdealGas::R / molarMass();
}
/*!
* \brief The density [kg/m^3] of liquid hydrogen at a given pressure and temperature.
*/
static Scalar liquidDensity(Scalar temperature, Scalar pressure)
{ DUNE_THROW(NotImplemented, "liquidDensity for H2"); }
/*
* \brief The pressure of liquid hydrogen at a given density and
* temperature [Pa].
*/
static Scalar liquidPressure(Scalar temperature, Scalar density)
{ DUNE_THROW(NotImplemented, "liquidPressure for H2"); }
/*!
* \brief Specific enthalpy [J/kg] of pure liquid H2 .
*/

53
dumux/new_material/components/h2o.hh
View File

@@ -248,7 +248,33 @@ public:
*/
static Scalar gasPressure(Scalar temperature, Scalar density)
{
DUNE_THROW(NotImplemented, "H2O::gasPressure");
Valgrind::CheckDefined(temperature);
Valgrind::CheckDefined(density);
// We use the newton method for this. For the initial value we
// assume steam to be an ideal gas
Scalar pressure = IdealGas<Scalar>::pressure(temperature, density/molarMass());
Scalar eps = pressure*1e-7;
Scalar deltaP = pressure*2;
Valgrind::CheckDefined(pressure);
Valgrind::CheckDefined(deltaP);
for (int i = 0; i < 5 && std::abs(pressure*1e-9) < std::abs(deltaP); ++i) {
Scalar f = gasDensity(temperature, pressure) - density;
Scalar df_dp;
df_dp = gasDensity(temperature, pressure + eps);
df_dp -= gasDensity(temperature, pressure - eps);
df_dp /= 2*eps;
deltaP = - f/df_dp;
pressure += deltaP;
Valgrind::CheckDefined(pressure);
Valgrind::CheckDefined(deltaP);
}
return pressure;
}
/*!
@@ -277,7 +303,8 @@ public:
}
/*!
* \brief The pressure of liquid at a given density and temperature [Pa].
* \brief The pressure of liquid water at a given density and
* temperature [Pa].
*
* See:
*
@@ -287,7 +314,27 @@ public:
*/
static Scalar liquidPressure(Scalar temperature, Scalar density)
{
DUNE_THROW(NotImplemented, "H2O::liquidPressure");
// We use the newton method for this. For the initial value we
// assume the pressure to be 10% higher than the vapor
// pressure
Scalar pressure = 1.1*vaporPressure(temperature);
Scalar eps = pressure*1e-7;
Scalar deltaP = pressure*2;
for (int i = 0; i < 5 && std::abs(pressure*1e-9) < std::abs(deltaP); ++i) {
Scalar f = liquidDensity(temperature, pressure) - density;
Scalar df_dp;
df_dp = liquidDensity(temperature, pressure + eps);
df_dp -= liquidDensity(temperature, pressure - eps);
df_dp /= 2*eps;
deltaP = - f/df_dp;
pressure += deltaP;
}
return pressure;
}
/*!

18
dumux/new_material/components/n2.hh
View File

@@ -125,7 +125,16 @@ public:
static Scalar gasDensity(Scalar temperature, Scalar pressure)
{
// Assume an ideal gas
return IdealGas::density(molarMass(), pressure, temperature);
return IdealGas::density(molarMass(), temperature, pressure);
}
/*
* \brief The pressure of gaseous N2 at a given density and temperature [Pa].
*/
static Scalar gasPressure(Scalar temperature, Scalar density)
{
// Assume an ideal gas
return IdealGas::pressure(temperature, density/molarMass());
}
/*!
@@ -134,6 +143,13 @@ public:
static Scalar liquidDensity(Scalar temperature, Scalar pressure)
{ DUNE_THROW(NotImplemented, "liquidDensity for N2"); }
/*
* \brief The pressure of liquid nitrogen at a given density and
* temperature [Pa].
*/
static Scalar liquidPressure(Scalar temperature, Scalar density)
{ DUNE_THROW(NotImplemented, "liquidPressure for N2"); }
/*!
* \brief Specific enthalpy [J/kg] of pure nitrogen gas.
*/

24
dumux/new_material/components/o2.hh
View File

@@ -119,9 +119,18 @@ public:
static Scalar density(Scalar temperature, Scalar pressure)
{
// Assume an ideal gas
return IdealGas::density(molarMass(), pressure, temperature);
return IdealGas::density(molarMass(), temperature, pressure);
}
/*
* \brief The pressure of gaseous N2 at a given density and temperature [Pa].
*/
static Scalar gasPressure(Scalar temperature, Scalar density)
{
// Assume an ideal gas
return IdealGas::pressure(temperature, density/molarMass());
}
/*!
* \brief Specific enthalpy [J/kg] of pure nitrogen gas.
*/
@@ -132,6 +141,19 @@ public:
return temperature * (cvHat + 1) * IdealGas::R / molarMass();
}
/*!
* \brief The density [kg/m^3] of gaseous O2 at a given pressure and temperature.
*/
static Scalar liquidDensity(Scalar temperature, Scalar pressure)
{ DUNE_THROW(NotImplemented, "liquidDensity for O2"); }
/*
* \brief The pressure of liquid oxygen at a given density and
* temperature [Pa].
*/
static Scalar liquidPressure(Scalar temperature, Scalar density)
{ DUNE_THROW(NotImplemented, "liquidPressure for O2"); }
/*!
* \brief Specific enthalpy [J/kg] of pure liquid O2.
*/

25
dumux/new_material/components/simpleh2o.hh
View File

@@ -158,11 +158,21 @@ public:
}
/*!
* \brief The density of steam at a given pressure and temperature [kg/m^3].
*/
* \brief The density [kg/m^3] of steam at a given pressure and temperature.
*/
static Scalar gasDensity(Scalar temperature, Scalar pressure)
{
return 1/IdealGas::density(molarMass(), temperature, pressure);
// Assume an ideal gas
return IdealGas::density(molarMass(), temperature, pressure);
}
/*
* \brief The pressure of steam at a given density and temperature [Pa].
*/
static Scalar gasPressure(Scalar temperature, Scalar density)
{
// Assume an ideal gas
return IdealGas::pressure(temperature, density/molarMass());
}
/*!
@@ -173,6 +183,15 @@ public:
return 1000;
}
/*
* \brief The pressure of water at a given density and temperature [Pa].
*/
static Scalar liquidPressure(Scalar temperature, Scalar density)
{
DUNE_THROW(InvalidStateException,
"The liquid pressure is undefined for incompressible fluids");
}
/*!
* \brief The dynamic viscosity [N/m^3*s] of steam.
*/

216
dumux/new_material/components/tabulatedcomponent.hh
View File

@@ -45,9 +45,14 @@ public:
pressMin_ = pressMin;
pressMax_ = pressMax;
nPress_ = nPress;
nDensity_ = nPress_;
// allocate the arrays
vaporPressure_ = new Scalar[nTemp_];
minGasDensity__ = new Scalar[nTemp_];
maxGasDensity__ = new Scalar[nTemp_];
minLiquidDensity__ = new Scalar[nTemp_];
maxLiquidDensity__ = new Scalar[nTemp_];
gasEnthalpy_ = new Scalar[nTemp_*nPress_];
liquidEnthalpy_ = new Scalar[nTemp_*nPress_];
@@ -57,6 +62,8 @@ public:
liquidDensity_ = new Scalar[nTemp_*nPress_];
gasViscosity_ = new Scalar[nTemp_*nPress_];
liquidViscosity_ = new Scalar[nTemp_*nPress_];
gasPressure_ = new Scalar[nTemp_*nDensity_];
liquidPressure_ = new Scalar[nTemp_*nDensity_];
assert(std::numeric_limits<Scalar>::has_quiet_NaN);
Scalar NaN = std::numeric_limits<Scalar>::quiet_NaN();
@@ -70,6 +77,8 @@ public:
Scalar pgMax = maxGasPressure_(iT);
Scalar pgMin = minGasPressure_(iT);
// fill the temperature, pressure gas arrays
for (unsigned iP = 0; iP < nPress_; ++ iP) {
Scalar pressure = iP * (pgMax - pgMin)/(nPress_ - 1) + pgMin;
@@ -88,6 +97,26 @@ public:
catch (NumericalProblem) { gasViscosity_[i] = NaN; };
};
// calculate the minimum and maximum values for the gas
// densities
minGasDensity__[iT] = RawComponent::gasDensity(temperature, pgMin);
maxGasDensity__[iT] = RawComponent::gasDensity(temperature, pgMax);
// fill the temperature, density gas arrays
for (unsigned iRho = 0; iRho < nDensity_; ++ iRho) {
Scalar density =
iRho * (maxGasDensity__[iT] - minGasDensity__[iT])
/
(nDensity_ - 1)
+
minGasDensity__[iT];
unsigned i = iT + iRho*nTemp_;
try { gasPressure_[i] = RawComponent::gasPressure(temperature, density); }
catch (NumericalProblem) { gasPressure_[i] = NaN; };
};
Scalar plMin = minLiquidPressure_(iT);
Scalar plMax = maxLiquidPressure_(iT);
for (unsigned iP = 0; iP < nPress_; ++ iP) {
@@ -107,6 +136,26 @@ public:
try { liquidViscosity_[i] = RawComponent::liquidViscosity(temperature, pressure); }
catch (NumericalProblem) { liquidViscosity_[i] = NaN; };
}
// calculate the minimum and maximum values for the liquid
// densities
minLiquidDensity__[iT] = RawComponent::liquidDensity(temperature, plMin);
maxLiquidDensity__[iT] = RawComponent::liquidDensity(temperature, plMax);
// fill the temperature, density liquid arrays
for (unsigned iRho = 0; iRho < nDensity_; ++ iRho) {
Scalar density =
iRho * (maxLiquidDensity__[iT] - minLiquidDensity__[iT])
/
(nDensity_ - 1)
+
minLiquidDensity__[iT];
unsigned i = iT + iRho*nTemp_;
try { liquidPressure_[i] = RawComponent::liquidPressure(temperature, density); }
catch (NumericalProblem) { liquidPressure_[i] = NaN; };
};
}
}
@@ -219,6 +268,39 @@ public:
return result;
}
/*!
* \brief The pressure of gas at a given density and temperature [Pa].
*/
static Scalar gasPressure(Scalar temperature, Scalar density)
{
Scalar result = interpolateGasTRho_(liquidPressure_,
temperature,
density);
if (std::isnan(result)) {
std::cerr << "forward gasPressure("<<temperature<<", "<<density<<")\n";
return RawComponent::gasPressure(temperature,
density);
}
return result;
};
/*!
* \brief The pressure of liquid at a given density and temperature [Pa].
*/
static Scalar liquidPressure(Scalar temperature, Scalar density)
{
Scalar result = interpolateLiquidTRho_(liquidPressure_,
temperature,
density);
if (std::isnan(result)) {
std::cerr << "forward liquidPressure("<<temperature<<", "<<density<<")\n";
return RawComponent::liquidPressure(temperature,
density);
}
return result;
};
/*!
* \brief The density of gas at a given pressure and temperature
* [kg/m^3].
@@ -353,6 +435,63 @@ private:
values[(iT + 1) + (iP2 + 1)*nTemp_]*( alphaT)*( alphaP2);
}
// returns an interpolated value for gas depending on
// temperature and density
static Scalar interpolateGasTRho_(const Scalar *values, Scalar T, Scalar rho)
{
Scalar alphaT = tempIdx_(T);
if (alphaT < 0 || alphaT >= nTemp_ - 1) {
// std::cerr << __LINE__ << " T: " << T << "\n";
return std::numeric_limits<Scalar>::quiet_NaN();
}
unsigned iT = (unsigned) alphaT;
alphaT -= iT;
Scalar alphaP1 = densityGasIdx_(rho, iT);
Scalar alphaP2 = densityGasIdx_(rho, iT + 1);
unsigned iP1 = std::min(nDensity_ - 2, (unsigned) alphaP1);
alphaP1 -= iP1;
unsigned iP2 = std::min(nDensity_ - 2, (unsigned) alphaP2);
alphaP2 -= iP2;
return
values[(iT ) + (iP1 )*nTemp_]*(1 - alphaT)*(1 - alphaP1) +
values[(iT ) + (iP1 + 1)*nTemp_]*(1 - alphaT)*( alphaP1) +
values[(iT + 1) + (iP2 )*nTemp_]*( alphaT)*(1 - alphaP2) +
values[(iT + 1) + (iP2 + 1)*nTemp_]*( alphaT)*( alphaP2);
}
// returns an interpolated value for liquid depending on
// temperature and density
static Scalar interpolateLiquidTRho_(const Scalar *values, Scalar T, Scalar rho)
{
Scalar alphaT = tempIdx_(T);
if (alphaT < 0 || alphaT >= nTemp_ - 1) {
// std::cerr << __LINE__ << " T: " << T << "\n";
return std::numeric_limits<Scalar>::quiet_NaN();
}
unsigned iT = (unsigned) alphaT;
alphaT -= iT;
Scalar alphaP1 = densityLiquidIdx_(rho, iT);
Scalar alphaP2 = densityLiquidIdx_(rho, iT + 1);
unsigned iP1 = std::min(nDensity_ - 2, (unsigned) alphaP1);
alphaP1 -= iP1;
unsigned iP2 = std::min(nDensity_ - 2, (unsigned) alphaP2);
alphaP2 -= iP2;
return
values[(iT ) + (iP1 )*nTemp_]*(1 - alphaT)*(1 - alphaP1) +
values[(iT ) + (iP1 + 1)*nTemp_]*(1 - alphaT)*( alphaP1) +
values[(iT + 1) + (iP2 )*nTemp_]*( alphaT)*(1 - alphaP2) +
values[(iT + 1) + (iP2 + 1)*nTemp_]*( alphaT)*( alphaP2);
}
// returns the index of an entry in a temperature field
static Scalar tempIdx_(Scalar temperature)
{
@@ -373,6 +512,22 @@ private:
Scalar pressMax = maxGasPressure_(tempIdx);
return (nPress_ - 1)*(pressure - pressMin_)/(pressMax - pressMin_);
}
// returns the index of an entry in a density field
static Scalar densityLiquidIdx_(Scalar density, unsigned tempIdx)
{
Scalar densityMin = minLiquidDensity_(tempIdx);
Scalar densityMax = maxLiquidDensity_(tempIdx);
return (nDensity_ - 1)*(density - densityMin)/(densityMax - densityMin);
}
// returns the index of an entry in a density field
static Scalar densityGasIdx_(Scalar density, unsigned tempIdx)
{
Scalar densityMin = minGasDensity_(tempIdx);
Scalar densityMax = maxGasDensity_(tempIdx);
return (nDensity_ - 1)*(density - densityMin)/(densityMax - densityMin);
}
// returns the minimum tabulized liquid pressure at a given
// temperature index
@@ -394,8 +549,39 @@ private:
static Scalar maxGasPressure_(int tempIdx)
{ return std::min<Scalar>(pressMax_, vaporPressure_[tempIdx] * 1.1); }
// returns the minimum tabulized liquid density at a given
// temperature index
static Scalar minLiquidDensity_(int tempIdx)
{ return minLiquidDensity__[tempIdx]; }
// returns the maximum tabulized liquid density at a given
// temperature index
static Scalar maxLiquidDensity_(int tempIdx)
{ return maxLiquidDensity__[tempIdx]; }
// returns the minumum tabulized gas density at a given
// temperature index
static Scalar minGasDensity_(int tempIdx)
{ return minGasDensity__[tempIdx]; }
// returns the maximum tabulized gas density at a given
// temperature index
static Scalar maxGasDensity_(int tempIdx)
{ return maxGasDensity__[tempIdx]; }
// 1D fields with the temperature as degree of freedom
static Scalar *vaporPressure_;
static Scalar *minLiquidDensity__;
static Scalar *maxLiquidDensity__;
static Scalar *minGasDensity__;
static Scalar *maxGasDensity__;
// 2D fields with the temperature and pressure as degrees of
// freedom
static Scalar *gasEnthalpy_;
static Scalar *liquidEnthalpy_;
@@ -408,6 +594,12 @@ private:
static Scalar *gasViscosity_;
static Scalar *liquidViscosity_;
// 2D fields with the temperature and density as degrees of
// freedom
static Scalar *gasPressure_;
static Scalar *liquidPressure_;
// temperature, pressure and density ranges
static Scalar tempMin_;
static Scalar tempMax_;
static unsigned nTemp_;
@@ -415,11 +607,23 @@ private:
static Scalar pressMin_;
static Scalar pressMax_;
static unsigned nPress_;
static Scalar densityMin_;
static Scalar densityMax_;
static unsigned nDensity_;
};
template <class Scalar, class RawComponent>
Scalar* TabulatedComponent<Scalar, RawComponent>::vaporPressure_;
template <class Scalar, class RawComponent>
Scalar* TabulatedComponent<Scalar, RawComponent>::minLiquidDensity__;
template <class Scalar, class RawComponent>
Scalar* TabulatedComponent<Scalar, RawComponent>::maxLiquidDensity__;
template <class Scalar, class RawComponent>
Scalar* TabulatedComponent<Scalar, RawComponent>::minGasDensity__;
template <class Scalar, class RawComponent>
Scalar* TabulatedComponent<Scalar, RawComponent>::maxGasDensity__;
template <class Scalar, class RawComponent>
Scalar* TabulatedComponent<Scalar, RawComponent>::gasEnthalpy_;
template <class Scalar, class RawComponent>
Scalar* TabulatedComponent<Scalar, RawComponent>::liquidEnthalpy_;
@@ -436,6 +640,10 @@ Scalar* TabulatedComponent<Scalar, RawComponent>::gasViscosity_;
template <class Scalar, class RawComponent>
Scalar* TabulatedComponent<Scalar, RawComponent>::liquidViscosity_;
template <class Scalar, class RawComponent>
Scalar* TabulatedComponent<Scalar, RawComponent>::gasPressure_;
template <class Scalar, class RawComponent>
Scalar* TabulatedComponent<Scalar, RawComponent>::liquidPressure_;
template <class Scalar, class RawComponent>
Scalar TabulatedComponent<Scalar, RawComponent>::tempMin_;
template <class Scalar, class RawComponent>
Scalar TabulatedComponent<Scalar, RawComponent>::tempMax_;
@@ -447,6 +655,12 @@ template <class Scalar, class RawComponent>
Scalar TabulatedComponent<Scalar, RawComponent>::pressMax_;
template <class Scalar, class RawComponent>
unsigned TabulatedComponent<Scalar, RawComponent>::nPress_;
template <class Scalar, class RawComponent>
Scalar TabulatedComponent<Scalar, RawComponent>::densityMin_;
template <class Scalar, class RawComponent>
Scalar TabulatedComponent<Scalar, RawComponent>::densityMax_;
template <class Scalar, class RawComponent>
unsigned TabulatedComponent<Scalar, RawComponent>::nDensity_;
} // end namepace

264
dumux/new_material/fluidsystems/h2o_n2_system.hh
View File

@@ -69,10 +69,6 @@ public:
static const int H2OIdx = 0;
static const int N2Idx = 1;
H2O_N2_System()
{
}
static void init()
{
std::cout << "Initializing tables for the H2O fluid properties.\n";
@@ -93,7 +89,7 @@ public:
}
/*!
* \brief Return the molar mass of a component in [kg/mol].
* \brief Return the molar mass of a component [kg/mol].
*/
static Scalar molarMass(int compIdx)
{
@@ -105,37 +101,78 @@ public:
}
/*!
* \brief Return the vapor pressure of a component in [Pa].
* \brief Given the gas phase's composition, temperature and
* pressure, compute the partial presures of all components
* in [Pa] and assign it to the PhaseState.
*
* This is required for models which cannot calculate the the
* partial pressures of the components in the gas phase from the
* degasPressure(). To use this method, the PhaseState must have a
* setPartialPressure(componentIndex, pressure) method.
*/
static Scalar vaporPressure(int compIdx,
Scalar temperature)
{
switch (compIdx) {
case H2OIdx: return H2O::vaporPressure(temperature);
case N2Idx: return N2::vaporPressure(temperature);
};
DUNE_THROW(InvalidStateException, "Invalid component index " << compIdx);
}
template <class PhaseState>
static void computePartialPressures(Scalar temperature,
Scalar pg,
PhaseState &phaseState)
{
Scalar X1 = phaseState.massFrac(gPhaseIdx, H2OIdx);
// We use the newton method for this. For the initial value we
// assume all components to be an ideal gas
Scalar pH2O =
phaseState.moleFrac(gPhaseIdx, H2OIdx) * pg;
Scalar eps = pg*1e-9;
Scalar deltaP = pH2O;
Valgrind::CheckDefined(pH2O);
Valgrind::CheckDefined(deltaP);
for (int i = 0; i < 5 && std::abs(deltaP/pg) > 1e-9; ++i) {
Scalar f =
H2O::gasDensity(temperature, pH2O)*(1 - 1/X1) +
N2::gasDensity(temperature, pg - pH2O);
Scalar df_dp;
df_dp =
H2O::gasDensity(temperature, pH2O + eps)*(1 - 1/X1) +
N2::gasDensity(temperature, pg - (pH2O + eps));
df_dp -=
H2O::gasDensity(temperature, pH2O - eps)*(1 - 1/X1) +
N2::gasDensity(temperature, pg - (pH2O - eps));
df_dp /=
2*eps;
deltaP = - f/df_dp;
pH2O += deltaP;
Valgrind::CheckDefined(pH2O);
Valgrind::CheckDefined(deltaP);
}
phaseState.setPartialPressure(H2OIdx, pH2O);
phaseState.setPartialPressure(N2Idx, pg - pH2O);
};
/*!
* \brief Given all mole fractions in a phase, return the phase
* \brief Given a phase's composition, temperature, pressure, and
* the partial pressures of all components, return its
* density [kg/m^3].
*/
template <class PhaseState>
static Scalar phaseDensity(int phaseIdx,
int temperature,
int pressure,
const PhaseState &phaseState)
{
switch (phaseIdx) {
case lPhaseIdx:
{
Scalar pVap = 0.9*H2O::vaporPressure(phaseState.temperature());
Scalar pressure = phaseState.phasePressure(lPhaseIdx);
Scalar pVap = 0.9*H2O::vaporPressure(temperature);
if (pressure < pVap)
pressure = pVap;
// See: Ochs 2008
// \todo: proper citation
Scalar rhoWater = H2O::liquidDensity(phaseState.temperature(),
Scalar rhoWater = H2O::liquidDensity(temperature,
pressure);
Scalar cWater = rhoWater/H2O::molarMass();
return
@@ -145,41 +182,37 @@ public:
}
case gPhaseIdx:
{
// assume ideal gas
Scalar avgMolarMass =
phaseState.moleFrac(gPhaseIdx, N2Idx)*N2::molarMass();
Scalar pWater = phaseState.partialPressure(H2OIdx);
Scalar pH2O = phaseState.partialPressure(H2OIdx);
Scalar pN2 = phaseState.partialPressure(N2Idx);
Scalar pVap = 1.1*H2O::vaporPressure(phaseState.temperature());
if (pWater > pVap)
pWater = pVap;
Scalar rhoWater = H2O::gasDensity(phaseState.temperature(), pWater);
// using the two partial pressures we can calculate the
// density of the phase. This assumes that Dalton's law is
// valid
return
rhoWater +
IdealGas::density(avgMolarMass,
phaseState.temperature(),
phaseState.phasePressure(gPhaseIdx) - pWater);
H2O::gasDensity(temperature, pH2O) +
N2::gasDensity(temperature, pN2);
};
}
DUNE_THROW(InvalidStateException, "Invalid phase index " << phaseIdx);
}
/*!
* \brief Return the viscosity of a phase.
* \brief Given a phase's composition, temperature and pressure,
* return its viscosity.
*/
template <class PhaseState>
template <class PhaseCompo>
static Scalar phaseViscosity(int phaseIdx,
const PhaseState &phaseState)
Scalar temperature,
Scalar pressure,
const PhaseCompo &phaseCompo)
{
if (phaseIdx == lPhaseIdx) {
Scalar pVap = 0.9*H2O::vaporPressure(phaseState.temperature());
Scalar pressure = phaseState.phasePressure(lPhaseIdx);
Scalar pVap = 0.9*H2O::vaporPressure(temperature);
if (pressure < pVap)
pressure = pVap;
// assume pure water for the liquid phase
// TODO: viscosity of mixture
return H2O::liquidViscosity(phaseState.temperature(),
return H2O::liquidViscosity(phaseCompo.temperature(),
pressure);
}
else {
@@ -196,10 +229,10 @@ public:
*/
Scalar muResult = 0;
const Scalar mu[numComponents] = {
H2O::gasViscosity(phaseState.temperature(),
H2O::vaporPressure(phaseState.temperature())),
N2::gasViscosity(phaseState.temperature(),
phaseState.phasePressure(gPhaseIdx))
H2O::gasViscosity(temperature,
H2O::vaporPressure(temperature)),
N2::gasViscosity(temperature,
pressure)
};
// molar masses
const Scalar M[numComponents] = {
@@ -214,9 +247,9 @@ public:
pow(M[i]/M[j], 1/4.0));
phiIJ *= phiIJ;
phiIJ /= sqrt(8*(1 + M[i]/M[j]));
divisor += phaseState.moleFrac(phaseIdx, j)*phiIJ;
divisor += phaseCompo.moleFrac(phaseIdx, j)*phiIJ;
}
muResult += phaseState.moleFrac(phaseIdx, i)*mu[i] / divisor;
muResult += phaseCompo.moleFrac(phaseIdx, i)*mu[i] / divisor;
}
return muResult;
@@ -224,34 +257,89 @@ public:
}
/*!
* \brief Returns the derivative of the equilibrium partial
* pressure \f$\partial p^\kappa_g / \partial x^\kappa_l\$
* to the mole fraction of a component in the liquid phase.
* \brief Return the pressure which a component degases from the
* liquid scaled to \f$100\%\f$ of the component.
*
* For solutions with only traces in a solvent this boils down to
* the inverse Henry constant for the solutes and the partial
* pressure for the solvent.
*/
template <class PhaseState>
static Scalar dPg_dxl(int compIdx,
const PhaseState &phaseState)
static Scalar degasPressure(int compIdx,
Scalar temperature,
Scalar pg)
{
switch (compIdx) {
case H2OIdx: return H2O::vaporPressure(phaseState.temperature());
case N2Idx: return BinaryCoeff::H2O_N2::henry(phaseState.temperature());
case H2OIdx: return H2O::vaporPressure(temperature);
case N2Idx: return BinaryCoeff::H2O_N2::henry(temperature);
};
DUNE_THROW(InvalidStateException, "Invalid component index " << compIdx);
}
/*!
* \brief Given a component's pressure and temperature, return its
* density in a phase [kg/m^3].
*/
static Scalar componentDensity(int phaseIdx,
int compIdx,
Scalar temperature,
Scalar pressure)
{
if (phaseIdx == lPhaseIdx) {
if (compIdx == H2OIdx)
return H2O::liquidDensity(temperature, pressure);
else if (compIdx == N2Idx)
return N2::liquidDensity(temperature, pressure);
DUNE_THROW(InvalidStateException, "Invalid component index " << compIdx);
}
else if (phaseIdx == gPhaseIdx) {
if (compIdx == H2OIdx) {
return H2O::gasDensity(temperature, pressure);
}
else if (compIdx == N2Idx)
return N2::gasDensity(temperature, pressure);
DUNE_THROW(InvalidStateException, "Invalid component index " << compIdx);
}
DUNE_THROW(InvalidStateException, "Invalid phase index " << phaseIdx);
};
/*!
* \brief Given all mole fractions, return the diffusion
* coefficent of a component in a phase.
* \brief Given a component's density and temperature, return the
* corresponding pressure in a phase [Pa].
*/
template <class PhaseState>
static Scalar componentPressure(int phaseIdx,
int compIdx,
Scalar temperature,
Scalar density)
{
if (phaseIdx == lPhaseIdx) {
if (compIdx == H2OIdx)
return H2O::liquidPressure(temperature, density);
else if (compIdx == N2Idx)
return N2::liquidPressure(temperature, density);
DUNE_THROW(InvalidStateException, "Invalid component index " << compIdx);
}
else if (phaseIdx == gPhaseIdx) {
if (compIdx == H2OIdx)
return H2O::gasPressure(temperature, density);
else if (compIdx == N2Idx)
return N2::gasPressure(temperature, density);
DUNE_THROW(InvalidStateException, "Invalid component index " << compIdx);
}
DUNE_THROW(InvalidStateException, "Invalid phase index " << phaseIdx);
};
/*!
* \brief Given a phase's composition, temperature and pressure,
* return the binary diffusion coefficent for components
* \f$i\f$ and \f$j\f$ in this phase.
*/
template <class PhaseCompo>
static Scalar diffCoeff(int phaseIdx,
int compIIdx,
int compJIdx,
const PhaseState &phaseState)
Scalar temperature,
Scalar pressure,
const PhaseCompo &phaseCompo)
{
if (compIIdx > compJIdx)
std::swap(compIIdx, compJIdx);
@@ -274,8 +362,8 @@ public:
switch (compIIdx) {
case H2OIdx:
switch (compJIdx) {
case N2Idx: return BinaryCoeff::H2O_N2::liquidDiffCoeff(phaseState.temperature(),
phaseState.phasePressure(phaseIdx));
case N2Idx: return BinaryCoeff::H2O_N2::liquidDiffCoeff(temperature,
pressure);
}
default:
DUNE_THROW(InvalidStateException,
@@ -286,8 +374,8 @@ public:
switch (compIIdx) {
case H2OIdx:
switch (compJIdx) {
case N2Idx: return BinaryCoeff::H2O_N2::gasDiffCoeff(phaseState.temperature(),
phaseState.phasePressure(phaseIdx));
case N2Idx: return BinaryCoeff::H2O_N2::gasDiffCoeff(temperature,
pressure);
}
}
}
@@ -299,92 +387,94 @@ public:
};
/*!
* \brief Given all mole fractions in a phase, return the specific
* phase enthalpy [J/kg].
* \brief Given a phase's composition, temperature and pressure,
* return its specific enthalpy [J/kg].
*/
template <class PhaseState>
template <class PhaseCompo>
static Scalar enthalpy(int phaseIdx,
const PhaseState &phaseState)
Scalar temperature,
Scalar pressure,
const PhaseCompo &phaseCompo)
{
Scalar temperature = phaseState.temperature();
if (phaseIdx == lPhaseIdx) {
Scalar pVap = 1.1*H2O::vaporPressure(temperature);
Scalar pWater = phaseState.phasePressure(lPhaseIdx);
Scalar pWater = pressure;
if (pWater < pVap)
pWater = pVap;
Scalar cN2 = phaseState.concentration(lPhaseIdx, N2Idx);
Scalar cN2 = phaseCompo.concentration(lPhaseIdx, N2Idx);
Scalar pN2 = IdealGas::pressure(temperature, cN2);
// TODO: correct way to deal with the solutes???
return
phaseState.massFrac(lPhaseIdx, H2OIdx)*
phaseCompo.massFrac(lPhaseIdx, H2OIdx)*
H2O::liquidEnthalpy(temperature, pWater)
+
phaseState.massFrac(lPhaseIdx, N2Idx)*
phaseCompo.massFrac(lPhaseIdx, N2Idx)*
N2::gasEnthalpy(temperature, pN2);
}
else {
Scalar pVap = 0.9*H2O::vaporPressure(temperature);
Scalar pWater = phaseState.partialPressure(H2OIdx);
Scalar pWater = phaseCompo.partialPressure(H2OIdx);
if (pWater > pVap)
pWater = pVap;
Scalar pN2 = phaseState.partialPressure(N2Idx);
Scalar pN2 = phaseCompo.partialPressure(N2Idx);
Scalar result = 0;
result +=
H2O::gasEnthalpy(temperature, pWater) *
phaseState.massFrac(gPhaseIdx, H2OIdx);
phaseCompo.massFrac(gPhaseIdx, H2OIdx);
result +=
N2::gasEnthalpy(temperature, pN2) *
phaseState.massFrac(gPhaseIdx, N2Idx);
phaseCompo.massFrac(gPhaseIdx, N2Idx);
return result;
}
}
/*!
* \brief Given all mole fractions in a phase, return the phase's
* internal energy [J/kg].
* \brief Given a phase's composition, temperature and pressure,
* return its specific internal energy [J/kg].
*/
template <class PhaseState>
template <class PhaseCompo>
static Scalar internalEnergy(int phaseIdx,
const PhaseState &phaseState)
Scalar temperature,
Scalar pressure,
const PhaseCompo &phaseCompo)
{
Scalar temperature = phaseState.temperature();
if (phaseIdx == lPhaseIdx) {
Scalar pVap = 1.1*H2O::vaporPressure(temperature);
Scalar pWater = phaseState.phasePressure(lPhaseIdx);
Scalar pWater = pressure;
if (pWater < pVap)
pWater = pVap;
Scalar cN2 = phaseState.concentration(lPhaseIdx, N2Idx);
Scalar cN2 = phaseCompo.concentration(lPhaseIdx, N2Idx);
Scalar pN2 = IdealGas::pressure(temperature, cN2);
// TODO: correct way to deal with the solutes???
return
phaseState.massFrac(lPhaseIdx, H2OIdx)*
phaseCompo.massFrac(lPhaseIdx, H2OIdx)*
H2O::liquidInternalEnergy(temperature, pWater)
+
phaseState.massFrac(lPhaseIdx, N2Idx)*
phaseCompo.massFrac(lPhaseIdx, N2Idx)*
N2::gasInternalEnergy(temperature, pN2);
}
else {
Scalar pVap = 0.9*H2O::vaporPressure(temperature);
Scalar pWater = phaseState.partialPressure(H2OIdx);
Scalar pWater = phaseCompo.partialPressure(H2OIdx);
if (pWater > pVap)
pWater = pVap;
Scalar pN2 = phaseState.partialPressure(N2Idx);
Scalar pN2 = phaseCompo.partialPressure(N2Idx);
Scalar result = 0;
result +=
H2O::gasInternalEnergy(temperature, pWater)*
phaseState.massFrac(gPhaseIdx, H2OIdx);
phaseCompo.massFrac(gPhaseIdx, H2OIdx);
result +=
N2::gasInternalEnergy(temperature, pN2)*
phaseState.massFrac(gPhaseIdx, N2Idx);
phaseCompo.massFrac(gPhaseIdx, N2Idx);
return result;
}