diff --git a/configure.ac b/configure.ac
index dae11aec8a402fcd7dfc06c3574c458fcbfd6cc2..27105dd96b78206ceb5ae928790b32293a812b18 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 ac0afef8da56fc2a5217d18f2aafa37a289fa058..a7883d8c498d46f87a0016899d90bc44bc00209e 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 9624803aa90d60960651925cb7c2589385d40136..86e86c2b787cb5a67f10657b2770bbaab21a8936 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 1cfabdc6206809f78f0135fc523a5ff06cc17139..5b1f87c58132e3442233f56da71f735dc676da64 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 6fe2fa3196ba7bf50dbc1c883840b69746ceb1d8..0000000000000000000000000000000000000000
--- 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 0000000000000000000000000000000000000000..b7f80609b5b108f236a3b3b3cd65c450d9f78e4e
--- /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 0000000000000000000000000000000000000000..f04d00ba23ca1eb9f817bd2f15f8a75367260442
--- /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 9e4e2db118421116cbae50609290a70760720a9d..1e361cae2c1b7a7e4b61f6f74ccee6675dcdadb9 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 0afe5519f472de05150278ec57a2f5274fc423e0..15712c1c9ba1348e6f7b38edc36b9ac498d4fb31 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>