From 39434007a9d9e10e650fec120badbdd2a3bc964a Mon Sep 17 00:00:00 2001
From: Andreas Lauser <and@poware.org>
Date: Thu, 26 Jan 2012 18:04:04 +0000
Subject: [PATCH] add a test for the peng robinson EOS
This implied to move the SPE5 fluid system to the main
repository. alas, the test itself does not work yet because the flash
solver blows up.
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@7530 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
configure.ac | 17 +-
dumux/material/eos/pengrobinson.hh | 6 +-
dumux/material/eos/pengrobinsonmixture.hh | 2 +-
.../material/eos/pengrobinsonparamsmixture.hh | 171 +++++--
dumux/material/eos/pengrobinsonparamspure.hh | 158 ------
.../material/fluidsystems/spe5fluidsystem.hh | 470 ++++++++++++++++++
.../fluidsystems/spe5parametercache.hh | 351 +++++++++++++
test/material/Makefile.am | 2 +-
.../material/fluidsystems/checkfluidsystem.hh | 1 +
9 files changed, 959 insertions(+), 219 deletions(-)
delete mode 100644 dumux/material/eos/pengrobinsonparamspure.hh
create mode 100644 dumux/material/fluidsystems/spe5fluidsystem.hh
create mode 100644 dumux/material/fluidsystems/spe5parametercache.hh
diff --git a/configure.ac b/configure.ac
index dae11aec8a..27105dd96b 100644
--- a/configure.ac
+++ b/configure.ac
@@ -43,10 +43,10 @@ AC_CONFIG_FILES([dumux.pc
dumux/decoupled/2p2c/Makefile
dumux/decoupled/common/Makefile
dumux/decoupled/common/fv/Makefile
- dumux/freeflow/Makefile
- dumux/freeflow/stokes/Makefile
- dumux/freeflow/stokes2c/Makefile
- dumux/freeflow/stokes2cni/Makefile
+ dumux/freeflow/Makefile
+ dumux/freeflow/stokes/Makefile
+ dumux/freeflow/stokes2c/Makefile
+ dumux/freeflow/stokes2cni/Makefile
dumux/io/Makefile
dumux/linear/Makefile
dumux/material/Makefile
@@ -85,12 +85,13 @@ AC_CONFIG_FILES([dumux.pc
test/decoupled/2padaptive/Makefile
test/decoupled/2p2c/Makefile
test/freeflow/Makefile
- test/freeflow/stokes/Makefile
- test/freeflow/stokes2c/Makefile
- test/freeflow/stokes2cni/Makefile
- test/material/Makefile
+ test/freeflow/stokes/Makefile
+ test/freeflow/stokes2c/Makefile
+ test/freeflow/stokes2cni/Makefile
+ test/material/Makefile
test/material/tabulation/Makefile
test/material/ncpflash/Makefile
+ test/material/pengrobinson/Makefile
test/material/fluidsystems/Makefile
tutorial/Makefile
])
diff --git a/dumux/material/eos/pengrobinson.hh b/dumux/material/eos/pengrobinson.hh
index ac0afef8da..a7883d8c49 100644
--- a/dumux/material/eos/pengrobinson.hh
+++ b/dumux/material/eos/pengrobinson.hh
@@ -247,7 +247,7 @@ public:
std::pow(tmp, expo);
return fugCoeff;
- };
+ }
/*!
* \brief Returns the fugacity coefficient for a given pressure
@@ -463,7 +463,7 @@ protected:
Scalar f2 = (tau*(-0.64771 + std::sqrt(tau)*(2.41539 - tau*4.26979)) + 3.25259*std::pow(tau, 5))/Tr;
return Component::criticalPressure()*std::exp(f0 + omega * (f1 + omega*f2));
- };
+ }
/*!
* \brief Returns the difference between the liquid and the gas phase
@@ -481,7 +481,7 @@ protected:
Scalar p,
Scalar VmLiquid,
Scalar VmGas)
- { return fugacity(params, T, p, VmLiquid) - fugacity(params, T, p, VmGas); };
+ { return fugacity(params, T, p, VmLiquid) - fugacity(params, T, p, VmGas); }
static Tabulated2DFunction<Scalar> criticalTemperature_;
static Tabulated2DFunction<Scalar> criticalPressure_;
diff --git a/dumux/material/eos/pengrobinsonmixture.hh b/dumux/material/eos/pengrobinsonmixture.hh
index 9624803aa9..86e86c2b78 100644
--- a/dumux/material/eos/pengrobinsonmixture.hh
+++ b/dumux/material/eos/pengrobinsonmixture.hh
@@ -154,7 +154,7 @@ public:
}
return fugCoeff;
- };
+ }
};
} // end namepace Dumux
diff --git a/dumux/material/eos/pengrobinsonparamsmixture.hh b/dumux/material/eos/pengrobinsonparamsmixture.hh
index 1cfabdc620..5b1f87c581 100644
--- a/dumux/material/eos/pengrobinsonparamsmixture.hh
+++ b/dumux/material/eos/pengrobinsonparamsmixture.hh
@@ -1,7 +1,5 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
- * Copyright (C) 2010 by Andreas Lauser *
+ * Copyright (C) 2010-2012 by Andreas Lauser *
* Institute of Hydraulic Engineering *
* University of Stuttgart, Germany *
* email: <givenname>.<name>@iws.uni-stuttgart.de *
@@ -22,49 +20,62 @@
/*!
* \file
*
- * \brief The Peng-Robinson parameters for a mixture
+ * \brief The mixing rule for the oil and the gas phases of the SPE5 problem.
+ *
+ * This problem comprises \f$H_2O\f$, \f$C_1\f$, \f$C_3\f$, \f$C_6\f$,
+ * \f$C_10\f$, \f$C_15\f$ and \f$C_20\f$ as components.
*
* See:
*
* R. Reid, et al.: The Properties of Gases and Liquids, 4th edition,
* McGraw-Hill, 1987, pp. 43-44
+ *
+ * and
+ *
+ * J.E. Killough, et al.: Fifth Comparative Solution Project:
+ * Evaluation of Miscible Flood Simulators, Ninth SPE Symposium on
+ * Reservoir Simulation, 1987
*/
#ifndef DUMUX_PENG_ROBINSON_PARAMS_MIXTURE_HH
#define DUMUX_PENG_ROBINSON_PARAMS_MIXTURE_HH
#include "pengrobinsonparams.hh"
-#include "pengrobinson.hh"
#include <dumux/material/constants.hh>
namespace Dumux
{
+
/*!
- * \brief The mixing rule for the Peng-Robinson equation of state as given in Reid, p. 82
+ * \brief The mixing rule for the oil and the gas phases of the SPE5 problem.
+ *
+ * This problem comprises \f$H_2O\f$, \f$C_1\f$, \f$C_3\f$, \f$C_6\f$,
+ * \f$C_10\f$, \f$C_15\f$ and \f$C_20\f$ as components.
*
* See:
*
* R. Reid, et al.: The Properties of Gases and Liquids, 4th edition,
- * McGraw-Hill, 1987, p. 82
+ * McGraw-Hill, 1987, pp. 43-44
+ *
+ * and
+ *
+ * J.E. Killough, et al.: Fifth Comparative Solution Project:
+ * Evaluation of Miscible Flood Simulators, Ninth SPE Symposium on
+ * Reservoir Simulation, 1987
*/
-template <class Scalar, class StaticParams, int phaseIdx>
-class PengRobinsonParamsMixture : public PengRobinsonParams<Scalar>
+template <class Scalar, class FluidSystem, int phaseIdx, bool useSpe5Relations=false>
+class PengRobinsonParamsMixture
+ : public PengRobinsonParams<Scalar>
{
- typedef Dumux::PengRobinsonParams<Scalar> ParentType;
+ enum { numComponents = FluidSystem::numComponents };
// Peng-Robinson parameters for pure substances
typedef Dumux::PengRobinsonParams<Scalar> PureParams;
- // The Peng-Robinson EOS for this mixture
-
- // number of components of which the fluid is composed
- enum { numComponents = StaticParams::numComponents };
-
- // ideal gas constant
+ // the ideal gas constant
static constexpr Scalar R = Dumux::Constants<Scalar>::R;
public:
-
/*!
* \brief Update Peng-Robinson parameters for the pure components.
*/
@@ -76,7 +87,7 @@ public:
}
/*!
- * \brief Update Peng-Robinson parameters for the pure components.
+ * \brief Peng-Robinson parameters for the pure components.
*
* This method is given by the SPE5 paper.
*/
@@ -88,18 +99,37 @@ public:
// See: R. Reid, et al.: The Properties of Gases and Liquids,
// 4th edition, McGraw-Hill, 1987, p. 43
for (int i = 0; i < numComponents; ++i) {
- Scalar pc = StaticParams::criticalPressure(i);
- Scalar omega = StaticParams::acentricFactor(i);
- Scalar Tr = temperature/StaticParams::criticalTemperature(i);
- Scalar RTc = R*StaticParams::criticalTemperature(i);
- Scalar f_omega = 0.37464 + omega*(1.54226 - omega*0.26992);
+ Scalar pc = FluidSystem::criticalPressure(i);
+ Scalar omega = FluidSystem::acentricFactor(i);
+ Scalar Tr = temperature/FluidSystem::criticalTemperature(i);
+ Scalar RTc = R*FluidSystem::criticalTemperature(i);
+
+ Scalar f_omega;
+
+ if (useSpe5Relations) {
+ if (omega < 0.49) f_omega = 0.37464 + omega*(1.54226 + omega*(-0.26992));
+ else f_omega = 0.379642 + omega*(1.48503 + omega*(-0.164423 + omega*0.016666));
+ }
+ else
+ f_omega = 0.37464 + omega*(1.54226 - omega*0.26992);
+
+ Valgrind::CheckDefined(f_omega);
Scalar tmp = 1 + f_omega*(1 - std::sqrt(Tr));
- this->pureParams_[i].setA(0.4572355*RTc*RTc/pc
- *
- tmp*tmp);
- this->pureParams_[i].setB(0.0777961 * RTc / pc);
+ tmp = tmp*tmp;
+
+ Scalar a = 0.4572355*RTc*RTc/pc * tmp;
+ Scalar b = 0.0777961 * RTc / pc;
+ assert(std::isfinite(a));
+ assert(std::isfinite(b));
+
+ this->pureParams_[i].setA(a);
+ this->pureParams_[i].setB(b);
+ Valgrind::CheckDefined(this->pureParams_[i].a());
+ Valgrind::CheckDefined(this->pureParams_[i].b());
}
+
+ updateACache_();
}
/*!
@@ -110,34 +140,67 @@ public:
* this method!
*/
template <class FluidState>
- void updateMix(const FluidState &fluidState)
+ void updateMix(const FluidState &fs)
{
+ Scalar sumx = 0.0;
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ sumx += fs.moleFraction(phaseIdx, compIdx);
+ sumx = std::max(1e-10, sumx);
+
// Calculate the Peng-Robinson parameters of the mixture
//
// See: R. Reid, et al.: The Properties of Gases and Liquids,
// 4th edition, McGraw-Hill, 1987, p. 82
Scalar a = 0;
Scalar b = 0;
- for (int i = 0; i < numComponents; ++i) {
- Scalar xi = fluidState.moleFrac(phaseIdx, i);
- for (int j = i; j < numComponents; ++j) {
- Scalar xj = fluidState.moleFrac(phaseIdx, j);
-
- // interaction coefficient as given in SPE5
- Scalar Psi = StaticParams::interactionCoefficient(i, j);
-
+ for (int compIIdx = 0; compIIdx < numComponents; ++compIIdx) {
+ Scalar xi = fs.moleFraction(phaseIdx, compIIdx) / sumx;
+ Valgrind::CheckDefined(xi);
+
+ for (int compJIdx = 0; compJIdx < numComponents; ++compJIdx) {
+ Scalar xj = fs.moleFraction(phaseIdx, compJIdx) / sumx;
+ Valgrind::CheckDefined(xj);
+
// mixing rule from Reid, page 82
- a += xi*xj*std::sqrt(pureParams_[i].a()*pureParams_[j].a())*(1 - Psi);
+ a += xi * xj * aCache_[compIIdx][compJIdx];
+
+ assert(std::isfinite(a));
}
// mixing rule from Reid, page 82
- b += xi * pureParams_[i].b();
+ b += xi * this->pureParams_[compIIdx].b();
+ assert(std::isfinite(b));
}
-
+
this->setA(a);
this->setB(b);
+ Valgrind::CheckDefined(this->a());
+ Valgrind::CheckDefined(this->b());
+
+ }
+
+ /*!
+ * \brief Calculates the "a" and "b" Peng-Robinson parameters for
+ * the mixture provided that only a single mole fraction
+ * was changed.
+ *
+ * The updatePure() method needs to be called _before_ calling
+ * this method!
+ */
+ template <class FluidState>
+ void updateSingleMoleFraction(const FluidState &fs,
+ int compIdx,
+ Scalar delta)
+ {
+ updateMix(fs);
}
+ /*!
+ * \brief Return the Peng-Robinson parameters of a pure substance,
+ */
+ const PureParams &pureParams(int compIdx) const
+ { return pureParams_[compIdx]; }
+
/*!
* \brief Returns the Peng-Robinson parameters for a pure component.
*/
@@ -154,23 +217,35 @@ public:
void checkDefined() const
{
#ifndef NDEBUG
- ParentType::checkDefined();
for (int i = 0; i < numComponents; ++i)
pureParams_[i].checkDefined();
+
+ Valgrind::CheckDefined(this->a());
+ Valgrind::CheckDefined(this->b());
#endif
};
- /*!
- * \brief Return the Peng-Robinson parameters of a pure substance,
- */
- const PureParams &pureParams(int compIdx) const
- { return pureParams_[compIdx]; }
-
-
protected:
PureParams pureParams_[numComponents];
-};
+private:
+ void updateACache_()
+ {
+ for (int compIIdx = 0; compIIdx < numComponents; ++ compIIdx) {
+ for (int compJIdx = 0; compJIdx < numComponents; ++ compJIdx) {
+ // interaction coefficient as given in SPE5
+ Scalar Psi = FluidSystem::interactionCoefficient(compIIdx, compJIdx);
+
+ aCache_[compIIdx][compJIdx] =
+ std::sqrt(this->pureParams_[compIIdx].a()
+ * this->pureParams_[compJIdx].a())
+ * (1 - Psi);
+ }
+ }
+ }
+
+ Scalar aCache_[numComponents][numComponents];
+};
} // end namepace
diff --git a/dumux/material/eos/pengrobinsonparamspure.hh b/dumux/material/eos/pengrobinsonparamspure.hh
deleted file mode 100644
index 6fe2fa3196..0000000000
--- a/dumux/material/eos/pengrobinsonparamspure.hh
+++ /dev/null
@@ -1,158 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * Copyright (C) 2010 by Andreas Lauser *
- * Institute of Hydraulic Engineering *
- * University of Stuttgart, Germany *
- * email: <givenname>.<name>@iws.uni-stuttgart.de *
- * *
- * This program is free software: you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation, either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * 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 *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief The Peng-Robinson parameters for a pure component
- *
- * See:
- *
- * R. Reid, et al.: The Properties of Gases and Liquids, 4th edition,
- * McGraw-Hill, 1987, pp. 43-44
- */
-#ifndef DUMUX_PENG_ROBINSON_PARAMS_PURE_HH
-#define DUMUX_PENG_ROBINSON_PARAMS_PURE_HH
-
-#include "pengrobinsonparams.hh"
-
-#include <dumux/material/constants.hh>
-
-namespace Dumux
-{
-/*!
- * \brief Stores, provides access to and calculates the Peng-Robinson
- * parameters of a pure component.
- *
- * See:
- *
- * R. Reid, et al.: The Properties of Gases and Liquids, 4th edition,
- * McGraw-Hill, 1987, pp. 43-44
- */
-template <class Scalar, class ComponentT>
-class PengRobinsonParamsPure : public PengRobinsonParams<Scalar>
-{
- typedef PengRobinsonParams<Scalar> ParentType;
-
- // the ideal gas constant
- static const Scalar R = Dumux::Constants<Scalar>::R;
-
-public:
- typedef ComponentT Component;
-
- /*!
- * \brief Calculates the "a" and "b" Peng-Robinson parameters for
- * the component.
- *
- * See:
- *
- * R. Reid, et al.: The Properties of Gases and Liquids, 4th edition,
- * McGraw-Hill, 1987, pp. 43-44
- */
- void update(Scalar T, Scalar p)
- {
- temperature_ = T;
- pressure_ = p;
-
- Scalar pc = Component::criticalPressure();
- Scalar omega = Component::acentricFactor();
- Scalar Tr = T/Component::criticalTemperature();
- Scalar RTc = R*Component::criticalTemperature();
- Scalar f_omega = 0.37464 + omega*(1.54226 - omega*0.26992);
- Scalar tmp = 1.0 + f_omega*(1.0 - std::sqrt(Tr));
- tmp = tmp*tmp;
-
- Scalar a = 0.4572355 * RTc*RTc/pc * tmp;
- Scalar b = 0.0777961 * RTc/pc;
-
- this->setA(a);
- this->setB(b);
- }
-
- /*!
- * \brief Sets the molar volume [m^3/mol] of the substance.
- *
- * The phaseIdx parameter is there to adhere to the common
- * interface with the multi-phase stuff and is just ignored.
- */
- void setMolarVolume(int phaseIdx, Scalar Vm)
- { setMolarVolume(Vm); }
-
- /*!
- * \brief Sets the molar volume [m^3/mol] of the substance.
- */
- void setMolarVolume(Scalar Vm)
- { molarVolume_ = Vm; }
-
- /*!
- * \brief Returns the temperature [K] of the system.
- *
- * The phaseIdx parameter is there to adhere to the common
- * interface with the multi-phase stuff and is just ignored.
- */
- Scalar temperature(int phaseIdx = 0) const
- { return temperature_; }
-
- /*!
- * \brief Returns the pressure [Pa] of the system.
- *
- * The phaseIdx parameter is there to adhere to the common
- * interface with the multi-phase stuff and is just ignored.
- */
- Scalar pressure(int phaseIdx = 0) const
- { return pressure_; }
-
- /*!
- * \brief Returns the molar volume [m^3/mol] of the substance.
- *
- * The phaseIdx parameter is there to adhere to the common
- * interface with the multi-phase stuff and is just ignored.
- */
- Scalar molarVolume(int phaseIdx = 0) const
- { return molarVolume_; }
-
- /*!
- * \brief Returns the attractive parameter "a" [Pa (m^3/mol)^2] for the cubic EOS.
- *
- * The phaseIdx parameter is there to adhere to the common
- * interface with the multi-phase stuff and is just ignored.
- */
- Scalar a(int phaseIdx = 0) const
- { return ParentType::a(); }
-
- /*!
- * \brief Returns the covolume of the substance [m^3/mol]
- *
- * The phaseIdx parameter is there to adhere to the common
- * interface with the multi-phase stuff and is just ignored.
- */
- Scalar b(int phaseIdx = 0) const
- { return ParentType::b(); }
-
-
-private:
- Scalar temperature_;
- Scalar pressure_;
- Scalar molarVolume_;
-};
-} // end namepace
-
-#endif
diff --git a/dumux/material/fluidsystems/spe5fluidsystem.hh b/dumux/material/fluidsystems/spe5fluidsystem.hh
new file mode 100644
index 0000000000..b7f80609b5
--- /dev/null
+++ b/dumux/material/fluidsystems/spe5fluidsystem.hh
@@ -0,0 +1,470 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2010 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * 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 *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief The mixing rule for the oil and the gas phases of the SPE5 problem.
+ *
+ * This problem comprises \f$H_2O\f$, \f$C_1\f$, \f$C_3\f$, \f$C_6\f$,
+ * \f$C_10\f$, \f$C_15\f$ and \f$C_20\f$ as components.
+ *
+ * See:
+ *
+ * J.E. Killough, et al.: Fifth Comparative Solution Project:
+ * Evaluation of Miscible Flood Simulators, Ninth SPE Symposium on
+ * Reservoir Simulation, 1987
+ */
+#ifndef DUMUX_SPE5_FLUID_SYSTEM_HH
+#define DUMUX_SPE5_FLUID_SYSTEM_HH
+
+#include "spe5parametercache.hh"
+
+#include <dumux/common/spline.hh>
+#include <dumux/material/eos/pengrobinsonmixture.hh>
+
+namespace Dumux
+{
+namespace FluidSystems
+{
+/*!
+ * \brief The fluid system for the SPE-5 benchmark problem.
+ *
+ * This problem comprises \f$H_2O\f$, \f$C_1\f$, \f$C_3\f$, \f$C_6\f$,
+ * \f$C_10\f$, \f$C_15\f$ and \f$C_20\f$ as components.
+ *
+ * See:
+ *
+ * J.E. Killough, et al.: Fifth Comparative Solution Project:
+ * Evaluation of Miscible Flood Simulators, Ninth SPE Symposium on
+ * Reservoir Simulation, 1987
+ */
+template <class Scalar>
+class Spe5
+{
+ typedef Dumux::FluidSystems::Spe5<Scalar> ThisType;
+
+ typedef typename Dumux::PengRobinsonMixture<Scalar, ThisType> PengRobinsonMixture;
+ typedef typename Dumux::PengRobinson<Scalar> PengRobinson;
+
+public:
+ typedef Dumux::Spe5ParameterCache<Scalar, ThisType> ParameterCache;
+
+ /****************************************
+ * Fluid phase parameters
+ ****************************************/
+
+ //! Number of phases in the fluid system
+ static const int numPhases = 3;
+
+ //! Index of the gas phase
+ static const int gPhaseIdx = 0;
+ //! Index of the water phase
+ static const int wPhaseIdx = 1;
+ //! Index of the oil phase
+ static const int oPhaseIdx = 2;
+
+ //! The component for pure water to be used
+ typedef Dumux::H2O<Scalar> H2O;
+
+ /*!
+ * \brief Return the human readable name of a fluid phase
+ */
+ static const char *phaseName(int phaseIdx)
+ {
+ static const char *name[] = {
+ "g",
+ "w",
+ "o",
+ };
+
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+ return name[phaseIdx];
+ }
+
+ /*!
+ * \brief Return whether a phase is liquid
+ */
+ static bool isLiquid(int phaseIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+ return phaseIdx != gPhaseIdx;
+ }
+
+ /*!
+ * \brief Return whether a phase is compressible
+ *
+ * In the SPE-5 problems all fluids are compressible...
+ */
+ static bool isCompressible(int phaseIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+ return true;
+ }
+
+ /*!
+ * \brief Returns true if and only if a fluid phase is assumed to
+ * be an ideal mixture.
+ *
+ * We define an ideal mixture as a fluid phase where the fugacity
+ * coefficients of all components times the pressure of the phase
+ * are indepent on the fluid composition. This assumtion is true
+ * if Henry's law and Rault's law apply. If you are unsure what
+ * this function should return, it is safe to return false. The
+ * only damage done will be (slightly) increased computation times
+ * in some cases.
+ */
+ static bool isIdealMixture(int phaseIdx)
+ {
+ // 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;
+ }
+
+ /****************************************
+ * Component related parameters
+ ****************************************/
+
+ //! Number of components in the fluid system
+ static const int numComponents = 7;
+
+ static const int H2OIdx = 0;
+ static const int C1Idx = 1;
+ static const int C3Idx = 2;
+ static const int C6Idx = 3;
+ static const int C10Idx = 4;
+ static const int C15Idx = 5;
+ static const int C20Idx = 6;
+
+ /*!
+ * \brief Return the human readable name of a component
+ */
+ static const char *componentName(int compIdx)
+ {
+ static const char *name[] = {
+ H2O::name(),
+ "C1",
+ "C3",
+ "C6",
+ "C10",
+ "C15",
+ "C20"
+ };
+
+ assert(0 <= compIdx && compIdx < numComponents);
+ return name[compIdx];
+ }
+
+ /*!
+ * \brief Return the molar mass of a component in [kg/mol].
+ */
+ static Scalar molarMass(int compIdx)
+ {
+ static const Scalar M[] = {
+ H2O::molarMass(),
+ 16.04e-3, // C1
+ 44.10e-3, // C3
+ 86.18e-3, // C6
+ 142.29e-3, // C10
+ 206.00e-3, // C15
+ 282.00e-3 // C20
+ };
+
+ assert(0 <= compIdx && compIdx < numComponents);
+ return M[compIdx];
+ }
+
+ /*!
+ * \brief Critical temperature of a component [K].
+ */
+ static Scalar criticalTemperature(int compIdx)
+ {
+ static const Scalar Tcrit[] = {
+ H2O::criticalTemperature(), // H2O
+ 343.0*5/9, // C1
+ 665.7*5/9, // C3
+ 913.4*5/9, // C6
+ 1111.8*5/9, // C10
+ 1270.0*5/9, // C15
+ 1380.0*5/9 // C20
+ };
+
+ assert(0 <= compIdx && compIdx < numComponents);
+ return Tcrit[compIdx];
+ };
+
+ /*!
+ * \brief Critical pressure of a component [Pa].
+ */
+ static Scalar criticalPressure(int compIdx)
+ {
+ static const Scalar pcrit[] = {
+ H2O::criticalPressure(), // H2O
+ 667.8 * 6894.7573, // C1
+ 616.3 * 6894.7573, // C3
+ 436.9 * 6894.7573, // C6
+ 304.0 * 6894.7573, // C10
+ 200.0 * 6894.7573, // C15
+ 162.0 * 6894.7573 // C20
+ };
+
+ assert(0 <= compIdx && compIdx < numComponents);
+ return pcrit[compIdx];
+ };
+
+ /*!
+ * \brief Molar volume of a component at the critical point [m^3/mol].
+ */
+ static Scalar criticalMolarVolume(int compIdx)
+ {
+ static const Scalar R = 8.314472;
+ static const Scalar vcrit[] = {
+ H2O::criticalMolarVolume(), // H2O
+ 0.290*R*criticalTemperature(1)/criticalPressure(1), // C1
+ 0.277*R*criticalTemperature(2)/criticalPressure(2), // C3
+ 0.264*R*criticalTemperature(3)/criticalPressure(3), // C6
+ 0.257*R*criticalTemperature(4)/criticalPressure(4), // C10
+ 0.245*R*criticalTemperature(5)/criticalPressure(5), // C15
+ 0.235*R*criticalTemperature(6)/criticalPressure(6) // C20
+ };
+
+ assert(0 <= compIdx && compIdx < numComponents);
+ return vcrit[compIdx];
+ };
+
+ /*!
+ * \brief The acentric factor of a component [].
+ */
+ static Scalar acentricFactor(int compIdx)
+ {
+ static const Scalar accFac[] = {
+ 0.344, // H2O (from Reid, et al.)
+ 0.0130, // C1
+ 0.1524, // C3
+ 0.3007, // C6
+ 0.4885, // C10
+ 0.6500, // C15
+ 0.8500 // C20
+ };
+
+ assert(0 <= compIdx && compIdx < numComponents);
+ return accFac[compIdx];
+ };
+
+ /*!
+ * \brief Returns the interaction coefficient for two components.
+ *
+ * The values are given by the SPE5 paper.
+ */
+ static Scalar interactionCoefficient(int comp1Idx, int comp2Idx)
+ {
+ int i = std::min(comp1Idx, comp2Idx);
+ int j = std::max(comp1Idx, comp2Idx);
+ if (i == C1Idx && (j == C15Idx || j == C20Idx))
+ return 0.05;
+ else if (i == C3Idx && (j == C15Idx || j == C20Idx))
+ return 0.005;
+ return 0;
+ }
+
+ /****************************************
+ * Methods which compute stuff
+ ****************************************/
+
+ /*!
+ * \brief Initialize the fluid system's static parameters
+ */
+ static void init()
+ {
+#warning "TODO: find the envelopes of 'a' and 'b'"
+ PengRobinson::init(/*aMin=*/1, /*aMax=*/100, /*na=*/100,
+ /*bMin=*/1e-6, /*bMax=*/1e-3, /*nb=*/100);
+ }
+
+ /*!
+ * \brief Calculate the density [kg/m^3] of a fluid phase
+ *
+ *
+ * \param fluidState An abitrary fluid state
+ * \param phaseIdx The index of the fluid phase to consider
+ */
+ template <class FluidState>
+ static Scalar density(const FluidState &fluidState,
+ const ParameterCache ¶mCache,
+ int phaseIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+ return fluidState.averageMolarMass(phaseIdx)/paramCache.molarVolume(phaseIdx);
+ }
+
+ /*!
+ * \brief Calculate the dynamic viscosity of a fluid phase [Pa*s]
+ */
+ template <class FluidState>
+ static Scalar viscosity(const FluidState &fs,
+ const ParameterCache ¶mCache,
+ int phaseIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx <= numPhases);
+
+ if (phaseIdx == gPhaseIdx) {
+ // 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)
+ // given by SPE-5: 0.7 centi-Poise = 0.0007 Pa s
+ return 0.7e-2 * 0.1;
+ else {
+ assert(phaseIdx == oPhaseIdx);
+ // given by SPE-5 in table on page 64. we use a constant
+ // viscosity, though...
+ return 0.208e-2 * 0.1;
+ }
+ }
+
+ /*!
+ * \brief Calculate the fugacity coefficient [Pa] of an individual
+ * component in a fluid phase
+ *
+ * The fugacity coefficient \f$\phi^\kappa_\alpha\f$ is connected
+ * to the fugacity \f$f^\kappa_\alpha\f$ and the component's mole
+ * fraction in a phase \f$x^\kappa_\alpha\f$ by means of the
+ * relation
+ *
+ * \f[ f^\kappa_\alpha = \phi^\kappa_\alpha \cdot x^\kappa_\alpha \cdot p_alpha \f]
+ */
+ template <class FluidState>
+ static Scalar fugacityCoefficient(const FluidState &fs,
+ const ParameterCache ¶mCache,
+ int phaseIdx,
+ int compIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx <= numPhases);
+ assert(0 <= compIdx && compIdx <= numComponents);
+
+ if (phaseIdx == oPhaseIdx || phaseIdx == gPhaseIdx)
+ return PengRobinsonMixture::computeFugacityCoefficient(fs,
+ paramCache,
+ phaseIdx,
+ compIdx);
+ else {
+ assert(phaseIdx == wPhaseIdx);
+ return
+ henryCoeffWater_(compIdx, fs.temperature(wPhaseIdx))
+ / fs.pressure(wPhaseIdx);
+ }
+ }
+
+
+ /*!
+ * \brief Calculate the binary molecular diffusion coefficient for
+ * a component in a fluid phase [mol^2 * s / (kg*m^3)]
+ *
+ * Molecular diffusion of a compoent $\kappa$ is caused by a
+ * gradient of the chemical potential and follows the law
+ *
+ * \f[ J = - D \grad mu_\kappa \f]
+ *
+ * where \f$\mu_\kappa\$ is the component's chemical potential,
+ * \f$D\f$ is the diffusion coefficient and \f$J\f$ is the
+ * diffusive flux. \f$mu_\kappa\f$ is connected to the component's
+ * fugacity \f$f_\kappa\f$ by the relation
+ *
+ * \f[ \mu_\kappa = R T_\alpha \mathrm{ln} \frac{f_\kappa}{p_\alpha} \f]
+ *
+ * where \f$p_\alpha\f$ and \f$T_\alpha\f$ are the fluid phase'
+ * pressure and temperature.
+ */
+ template <class FluidState>
+ static Scalar diffusionCoefficient(const FluidState &fs,
+ const ParameterCache ¶mCache,
+ int phaseIdx,
+ int compIdx)
+ { DUNE_THROW(Dune::NotImplemented, "Diffusion coefficients"); }
+
+ /*!
+ * \brief Given a phase's composition, temperature and pressure,
+ * return the binary diffusion coefficient for components
+ * \f$i\f$ and \f$j\f$ in this phase.
+ *
+ * \param fluidState An abitrary fluid state
+ * \param phaseIdx The index of the fluid phase to consider
+ * \param compIIdx The index of the first component to consider
+ * \param compJIdx The index of the second component to consider
+ */
+ template <class FluidState>
+ static Scalar binaryDiffusionCoefficient(const FluidState &fluidState,
+ const ParameterCache ¶mCache,
+ int phaseIdx,
+ int compIIdx,
+ int compJIdx)
+ { DUNE_THROW(Dune::NotImplemented, "Binary diffusion coefficients"); }
+
+ /*!
+ * \brief Given a phase's composition, temperature and pressure,
+ * calculate its specific enthalpy [J/kg].
+ */
+ template <class FluidState>
+ static Scalar enthalpy(const FluidState &fs,
+ const ParameterCache ¶mCache,
+ int phaseIdx)
+ { DUNE_THROW(Dune::NotImplemented, "Enthalpies"); }
+
+ template <class FluidState>
+ static Scalar thermalConductivity(const FluidState &fluidState,
+ const ParameterCache ¶mCache,
+ int phaseIdx)
+ { DUNE_THROW(Dune::NotImplemented, "Thermal conductivities"); }
+
+ template <class FluidState>
+ static Scalar heatCapacity(const FluidState &fluidState,
+ const ParameterCache ¶mCache,
+ int phaseIdx)
+ { DUNE_THROW(Dune::NotImplemented, "Heat capacities"); }
+
+
+private:
+ static Scalar henryCoeffWater_(int compIdx, Scalar temperature)
+ {
+ // use henry's law for the solutes and the vapor pressure for
+ // the solvent.
+ switch (compIdx) {
+ case H2OIdx: return H2O::vaporPressure(temperature);
+
+ // the values of the Henry constant for the solutes have
+ // been computed using the Peng-Robinson equation of state
+ // (-> slope of the component's fugacity function at
+ // almost 100% water content)
+ case C1Idx: return 5.57601e+09;
+ case C3Idx: return 1.89465e+10;
+ case C6Idx: return 5.58969e+12;
+ case C10Idx: return 4.31947e+17;
+ case C15Idx: return 4.27283e+28;
+ case C20Idx: return 3.39438e+36;
+ default: DUNE_THROW(Dune::InvalidStateException, "Unknown component index " << compIdx);
+ }
+ };
+};
+
+} // end namepace
+} // end namepace
+
+#endif
diff --git a/dumux/material/fluidsystems/spe5parametercache.hh b/dumux/material/fluidsystems/spe5parametercache.hh
new file mode 100644
index 0000000000..f04d00ba23
--- /dev/null
+++ b/dumux/material/fluidsystems/spe5parametercache.hh
@@ -0,0 +1,351 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2010 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * 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 *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Specifies the parameters required by the SPE5 problem which
+ * are dependend on the thermodynamic state.
+ */
+#ifndef SPE5_PARAMETER_CACHE_HH
+#define SPE5_PARAMETER_CACHE_HH
+
+#include <dumux/material/components/h2o.hh>
+#include <dumux/material/fluidsystems/parametercachebase.hh>
+
+#include <dumux/material/eos/pengrobinson.hh>
+#include <dumux/material/eos/pengrobinsonparamsmixture.hh>
+
+#include <assert.h>
+
+namespace Dumux
+{
+/*!
+ * \brief Specifies the parameters required by the SPE5 problem which
+ * are dependend on the thermodynamic state.
+ */
+template <class Scalar, class FluidSystem>
+class Spe5ParameterCache
+ : public Dumux::ParameterCacheBase<Spe5ParameterCache<Scalar, FluidSystem> >
+{
+ typedef Spe5ParameterCache<Scalar, FluidSystem> ThisType;
+ typedef Dumux::ParameterCacheBase<ThisType> ParentType;
+
+ typedef Dumux::PengRobinson<Scalar> PengRobinson;
+
+ enum { numPhases = FluidSystem::numPhases };
+ enum { numComponents = FluidSystem::numComponents };
+
+ enum { wPhaseIdx = FluidSystem::wPhaseIdx };
+ enum { oPhaseIdx = FluidSystem::oPhaseIdx };
+ enum { gPhaseIdx = FluidSystem::gPhaseIdx };
+
+public:
+ // types of the parameter objects for each phase
+ typedef Dumux::PengRobinsonParamsMixture<Scalar, FluidSystem, oPhaseIdx, /*useSpe5=*/true> OilPhaseParams;
+ typedef Dumux::PengRobinsonParamsMixture<Scalar, FluidSystem, gPhaseIdx, /*useSpe5=*/true> GasPhaseParams;
+
+ /*!
+ * \brief The constructor
+ */
+ Spe5ParameterCache()
+ {
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ VmUpToDate_[phaseIdx] = false;
+ Valgrind::SetUndefined(Vm_[phaseIdx]);
+ }
+ };
+
+ /*!
+ * \brief Update all parameters required by the fluid system to
+ * calculate some quantities for the phase.
+ */
+ template <class FluidState>
+ void updatePhase(const FluidState &fs,
+ int phaseIdx,
+ int except = ParentType::None)
+ {
+ updateEosParams(fs, phaseIdx, except);
+
+ // if we don't need to recalculate the molar volume, we exit
+ // here
+ if (VmUpToDate_[phaseIdx])
+ return;
+
+ // update the phase's molar volume
+ updateMolarVolume_(fs, phaseIdx);
+ }
+
+ /*!
+ * \brief Update all cached parameters of a specific fluid phase
+ * which depend on the mole fraction of a single component
+ *
+ * *Only* use this method if just a single component's
+ * concentration changed between two update*() calls. If more than
+ * one concentration changed, call updatePhaseComposition() of
+ * updatePhase()!
+ */
+ template <class FluidState>
+ void updateSingleMoleFraction(const FluidState &fs,
+ int phaseIdx,
+ int compIdx)
+ {
+ if (phaseIdx == oPhaseIdx)
+ oilPhaseParams_.updateSingleMoleFraction(fs,
+ compIdx,
+ fs.moleFraction(phaseIdx, compIdx)
+ - moleFrac_[phaseIdx][compIdx]);
+ else if (phaseIdx == gPhaseIdx)
+ gasPhaseParams_.updateSingleMoleFraction(fs,
+ compIdx,
+ fs.moleFraction(phaseIdx, compIdx)
+ - moleFrac_[phaseIdx][compIdx]);
+
+ // update the mole fraction which the parameters are
+ // calculated for
+ moleFrac_[phaseIdx][compIdx] = fs.moleFraction(phaseIdx, compIdx);
+
+ // update the phase's molar volume
+ updateMolarVolume_(fs, phaseIdx);
+ }
+
+ /*!
+ * \brief The Peng-Robinson attractive parameter for a phase.
+ */
+ Scalar a(int phaseIdx) const
+ {
+ switch (phaseIdx)
+ {
+ case oPhaseIdx: return oilPhaseParams_.a();
+ case gPhaseIdx: return gasPhaseParams_.a();
+ default:
+ DUNE_THROW(Dune::InvalidStateException,
+ "The a() parameter is only defined for "
+ "oil and gas phases");
+ };
+ }
+
+ /*!
+ * \brief The Peng-Robinson covolume for a phase.
+ */
+ Scalar b(int phaseIdx) const
+ {
+ switch (phaseIdx)
+ {
+ case oPhaseIdx: return oilPhaseParams_.b();
+ case gPhaseIdx: return gasPhaseParams_.b();
+ default:
+ DUNE_THROW(Dune::InvalidStateException,
+ "The b() parameter is only defined for "
+ "oil and gas phases");
+ };
+ }
+
+ /*!
+ * \brief The Peng-Robinson attractive parameter for a pure
+ * component given the same temperature and pressure of the
+ * phase.
+ */
+ Scalar aPure(int phaseIdx, int compIdx) const
+ {
+ switch (phaseIdx)
+ {
+ case oPhaseIdx: return oilPhaseParams_.pureParams(compIdx).a();
+ case gPhaseIdx: return gasPhaseParams_.pureParams(compIdx).a();
+ default:
+ DUNE_THROW(Dune::InvalidStateException,
+ "The a() parameter is only defined for "
+ "oil and gas phases");
+ };
+ }
+
+ /*!
+ * \brief The Peng-Robinson covolume for a pure component given
+ * the same temperature and pressure of the phase.
+ */
+ Scalar bPure(int phaseIdx, int compIdx) const
+ {
+ switch (phaseIdx)
+ {
+ case oPhaseIdx: return oilPhaseParams_.pureParams(compIdx).b();
+ case gPhaseIdx: return gasPhaseParams_.pureParams(compIdx).b();
+ default:
+ DUNE_THROW(Dune::InvalidStateException,
+ "The b() parameter is only defined for "
+ "oil and gas phases");
+ };
+ }
+
+ /*!
+ * \brief Returns the molar volume of a phase [m^3/mol]
+ */
+ Scalar molarVolume(int phaseIdx) const
+ { assert(VmUpToDate_[phaseIdx]); return Vm_[phaseIdx]; }
+
+
+ /*!
+ * \brief Returns the Peng-Robinson mixture parameters for the oil
+ * phase.
+ */
+ const OilPhaseParams &oilPhaseParams() const
+ { return oilPhaseParams_; }
+
+ /*!
+ * \brief Returns the Peng-Robinson mixture parameters for the gas
+ * phase.
+ */
+ const GasPhaseParams &gasPhaseParams() const
+ { return gasPhaseParams_; }
+
+ /*!
+ * \brief Update all parameters required by the equation of state to
+ * calculate some quantities for the phase.
+ */
+ template <class FluidState>
+ void updateEosParams(const FluidState &fs,
+ int phaseIdx,
+ int exceptQuantities = ParentType::None)
+ {
+ if (!(exceptQuantities & ParentType::Temperature))
+ {
+ updatePure_(fs, phaseIdx);
+ updateMix_(fs, phaseIdx);
+
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ moleFrac_[phaseIdx][compIdx] = fs.moleFraction(phaseIdx, compIdx);
+ VmUpToDate_[phaseIdx] = false;
+ }
+ else if (!(exceptQuantities & ParentType::Composition))
+ {
+ updateMix_(fs, phaseIdx);
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ moleFrac_[phaseIdx][compIdx] = fs.moleFraction(phaseIdx, compIdx);
+ VmUpToDate_[phaseIdx] = false;
+ }
+ else if (!(exceptQuantities & ParentType::Pressure))
+ VmUpToDate_[phaseIdx] = false;
+ }
+
+protected:
+ /*!
+ * \brief Update all parameters of a phase which only depend on
+ * temperature and/or pressure.
+ *
+ * This usually means the parameters for the pure components.
+ */
+ template <class FluidState>
+ void updatePure_(const FluidState &fs, int phaseIdx)
+ {
+ Scalar T = fs.temperature(phaseIdx);
+ Scalar p = fs.pressure(phaseIdx);
+
+ switch (phaseIdx)
+ {
+ case oPhaseIdx: oilPhaseParams_.updatePure(T, p); break;
+ case gPhaseIdx: gasPhaseParams_.updatePure(T, p); break;
+ //case wPhaseIdx: waterPhaseParams_.updatePure(phaseIdx, temperature, pressure);break;
+ }
+ }
+
+ /*!
+ * \brief Update all parameters of a phase which depend on the
+ * fluid composition. It is assumed that updatePure() has
+ * been called before this method.
+ *
+ * Here, the mixing rule kicks in.
+ */
+ template <class FluidState>
+ void updateMix_(const FluidState &fs, int phaseIdx)
+ {
+ Valgrind::CheckDefined(fs.averageMolarMass(phaseIdx));
+ switch (phaseIdx)
+ {
+ case oPhaseIdx:
+ oilPhaseParams_.updateMix(fs);
+ break;
+ case gPhaseIdx:
+ gasPhaseParams_.updateMix(fs);
+ break;
+ case wPhaseIdx:
+ break;
+ }
+ }
+
+ template <class FluidState>
+ void updateMolarVolume_(const FluidState &fs,
+ int phaseIdx)
+ {
+ VmUpToDate_[phaseIdx] = true;
+
+ // calculate molar volume of the phase (we will need this for the
+ // fugacity coefficients and the density anyway)
+ switch (phaseIdx) {
+ case gPhaseIdx: {
+ // 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] =
+ PengRobinson::computeMolarVolume(fs,
+ *this,
+ phaseIdx,
+ /*isGasPhase=*/true);
+ }
+ case oPhaseIdx: {
+ // 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] =
+ PengRobinson::computeMolarVolume(fs,
+ *this,
+ phaseIdx,
+ /*isGasPhase=*/false);
+
+ }
+ case wPhaseIdx: {
+ // 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 = fs.pressure(wPhaseIdx) - 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] = fs.averageMolarMass(wPhaseIdx)/waterDensity;
+ };
+ };
+ }
+
+ bool VmUpToDate_[numPhases];
+ Scalar Vm_[numPhases];
+ Scalar moleFrac_[numPhases][numComponents];
+
+ OilPhaseParams oilPhaseParams_;
+ GasPhaseParams gasPhaseParams_;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/material/Makefile.am b/test/material/Makefile.am
index 9e4e2db118..1e361cae2c 100644
--- a/test/material/Makefile.am
+++ b/test/material/Makefile.am
@@ -1,4 +1,4 @@
-SUBDIRS = fluidsystems tabulation ncpflash
+SUBDIRS = fluidsystems ncpflash pengrobinson tabulation
include $(top_srcdir)/am/global-rules
diff --git a/test/material/fluidsystems/checkfluidsystem.hh b/test/material/fluidsystems/checkfluidsystem.hh
index 0afe5519f4..15712c1c9b 100644
--- a/test/material/fluidsystems/checkfluidsystem.hh
+++ b/test/material/fluidsystems/checkfluidsystem.hh
@@ -36,6 +36,7 @@
#include <dumux/material/fluidsystems/1pfluidsystem.hh>
#include <dumux/material/fluidsystems/2pimmisciblefluidsystem.hh>
#include <dumux/material/fluidsystems/h2on2fluidsystem.hh>
+#include <dumux/material/fluidsystems/spe5fluidsystem.hh>
// include all fluid states
#include <dumux/material/fluidstates/pressureoverlayfluidstate.hh>
--
GitLab