diff --git a/configure.ac b/configure.ac
index 27105dd96b78206ceb5ae928790b32293a812b18..cbaa65656a47ec21a750d38f4fff2b2c801f1c06 100644
--- a/configure.ac
+++ b/configure.ac
@@ -21,6 +21,8 @@ AC_CONFIG_FILES([dumux.pc
dumux/boxmodels/2p2c/Makefile
dumux/boxmodels/2p2cni/Makefile
dumux/boxmodels/2pni/Makefile
+ dumux/boxmodels/3p3c/Makefile
+ dumux/boxmodels/3p3cni/Makefile
dumux/boxmodels/common/Makefile
dumux/boxmodels/MpNc/Makefile
dumux/boxmodels/MpNc/diffusion/Makefile
@@ -55,6 +57,7 @@ AC_CONFIG_FILES([dumux.pc
dumux/material/components/iapws/Makefile
dumux/material/fluidmatrixinteractions/Makefile
dumux/material/fluidmatrixinteractions/2p/Makefile
+ dumux/material/fluidmatrixinteractions/3p/Makefile
dumux/material/fluidmatrixinteractions/Mp/Makefile
dumux/material/spatialparameters/Makefile
dumux/material/old_fluidsystems/Makefile
@@ -73,6 +76,8 @@ AC_CONFIG_FILES([dumux.pc
test/boxmodels/1p2c/Makefile
test/boxmodels/2p2c/Makefile
test/boxmodels/2p2cni/Makefile
+ test/boxmodels/3p3c/Makefile
+ test/boxmodels/3p3cni/Makefile
test/boxmodels/MpNc/Makefile
test/boxmodels/richards/Makefile
test/common/Makefile
diff --git a/dumux/boxmodels/3p3cni/3p3cnifluxvariables.hh b/dumux/boxmodels/3p3cni/3p3cnifluxvariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..b7227e3e11f6039e0af742e81cae9370ad96b17c
--- /dev/null
+++ b/dumux/boxmodels/3p3cni/3p3cnifluxvariables.hh
@@ -0,0 +1,133 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * Copyright (C) 2011- by Holger Class *
+ * Copyright (C) 2008-2009 by Andreas Lauser *
+ * Copyright (C) 2008-2009 by Melanie Darcis *
+ * Copyright (C) 2008-2009 by Klaus Mosthaf *
+ * Copyright (C) 2008-2009 by Bernd Flemisch *
+ * 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 This file contains the data which is required to calculate
+ * all fluxes (mass of components and energy) over a face of a finite volume.
+ *
+ * This means pressure, concentration and temperature gradients, phase
+ * densities at the integration point, etc.
+ */
+#ifndef DUMUX_3P3CNI_FLUX_VARIABLES_HH
+#define DUMUX_3P3CNI_FLUX_VARIABLES_HH
+
+#include <dumux/common/math.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup ThreePThreeCNIModel
+ * \brief This template class contains the data which is required to
+ * calculate all fluxes (mass of components and energy) over a face of a finite
+ * volume for the non-isothermal three-phase, three-component model.
+ *
+ * This means pressure and concentration gradients, phase densities at
+ * the integration point, etc.
+ */
+template <class TypeTag>
+class ThreePThreeCNIFluxVariables : public ThreePThreeCFluxVariables<TypeTag>
+{
+ typedef ThreePThreeCFluxVariables<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
+
+ typedef typename GridView::ctype CoordScalar;
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+
+ enum {
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld,
+
+ numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
+ };
+
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename FVElementGeometry::SubControlVolume SCV;
+ typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
+
+ typedef Dune::FieldVector<CoordScalar, dimWorld> Vector;
+
+public:
+ /*
+ * \brief The constructor
+ *
+ * \param problem The problem
+ * \param element The finite element
+ * \param elemGeom The finite-volume geometry in the box scheme
+ * \param faceIdx The local index of the SCV (sub-control-volume) face
+ * \param elemDat The volume variables of the current element
+ */
+ ThreePThreeCNIFluxVariables(const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int scvfIdx,
+ const ElementVolumeVariables &elemDat)
+ : ParentType(problem, element, elemGeom, scvfIdx, elemDat)
+ {
+ // calculate temperature gradient using finite element
+ // gradients
+ Vector temperatureGrad(0);
+ Vector tmp(0.0);
+ for (int vertIdx = 0; vertIdx < elemGeom.numVertices; vertIdx++)
+ {
+ tmp = elemGeom.subContVolFace[scvfIdx].grad[vertIdx];
+ tmp *= elemDat[vertIdx].temperature();
+ temperatureGrad += tmp;
+ }
+
+ // The spatial parameters calculates the actual heat flux vector
+ problem.spatialParameters().matrixHeatFlux(tmp,
+ *this,
+ elemDat,
+ temperatureGrad,
+ element,
+ elemGeom,
+ scvfIdx);
+ // project the heat flux vector on the face's normal vector
+ normalMatrixHeatFlux_ = tmp*elemGeom.subContVolFace[scvfIdx].normal;
+ }
+
+ /*!
+ * \brief The total heat flux \f$\mathrm{[J/s]}\f$ due to heat conduction
+ * of the rock matrix over the sub-control volume's face in
+ * direction of the face normal.
+ */
+ Scalar normalMatrixHeatFlux() const
+ { return normalMatrixHeatFlux_; }
+
+private:
+ Scalar normalMatrixHeatFlux_;
+};
+
+} // end namepace
+
+#endif
diff --git a/dumux/boxmodels/3p3cni/3p3cniindices.hh b/dumux/boxmodels/3p3cni/3p3cniindices.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d7e566dd168f8d3a54c32cb9cfa1715970fb32fe
--- /dev/null
+++ b/dumux/boxmodels/3p3cni/3p3cniindices.hh
@@ -0,0 +1,58 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * Copyright (C) 2011 by Holger Class *
+ * Copyright (C) 2008-2010 by Andreas Lauser *
+ * Copyright (C) 2008-2009 by Melanie Darcis *
+ * Copyright (C) 2008-2009 by Klaus Mosthaf *
+ * Copyright (C) 2008-2009 by Bernd Flemisch *
+ * 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 Defines the indices used by the 3p3cni box model
+ */
+#ifndef DUMUX_3P3CNI_INDICES_HH
+#define DUMUX_3P3CNI_INDICES_HH
+
+#include <dumux/boxmodels/3p3c/3p3cindices.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup ThreePThreeCNIModel
+ */
+// \{
+
+/*!
+ * \brief Enumerations for the non-isothermal 3-phase 3-component model
+ *
+ * \tparam PVOffset The first index in a primary variable vector.
+ */
+template <class TypeTag, int formulation, int PVOffset>
+class ThreePThreeCNIIndices : public ThreePThreeCIndices<TypeTag, formulation, PVOffset>
+{
+public:
+ static const int temperatureIdx = PVOffset + 3; //! The index for temperature in primary variable vectors.
+ static const int energyEqIdx = PVOffset + 3; //! The index for energy in equation vectors.
+};
+
+}
+#endif
diff --git a/dumux/boxmodels/3p3cni/3p3cnilocalresidual.hh b/dumux/boxmodels/3p3cni/3p3cnilocalresidual.hh
new file mode 100644
index 0000000000000000000000000000000000000000..69a857e59c6bcead47052300477dc842b7b165af
--- /dev/null
+++ b/dumux/boxmodels/3p3cni/3p3cnilocalresidual.hh
@@ -0,0 +1,210 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * Copyright (C) 2011- by Holger Class *
+ * Copyright (C) 2008-2009 by Andreas Lauser *
+ * Copyright (C) 2008-2009 by Melanie Darcis *
+ * Copyright (C) 2008-2009 by Klaus Mosthaf *
+ * Copyright (C) 2008 by Bernd Flemisch *
+ * 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 Element-wise calculation of the Jacobian matrix for problems
+ * using the non-isothermal three-phase three-component box model.
+ *
+ */
+#ifndef DUMUX_NEW_3P3CNI_LOCAL_RESIDUAL_HH
+#define DUMUX_NEW_3P3CNI_LOCAL_RESIDUAL_HH
+
+#include <dumux/boxmodels/3p3c/3p3clocalresidual.hh>
+
+
+#include "3p3cnivolumevariables.hh"
+#include "3p3cnifluxvariables.hh"
+
+#include "3p3cniproperties.hh"
+
+namespace Dumux
+{
+/*!
+ * \ingroup ThreePThreeCNIModel
+ * \brief Element-wise calculation of the Jacobian matrix for problems
+ * using the three-phase three-component box model.
+ */
+template<class TypeTag>
+class ThreePThreeCNILocalResidual : public ThreePThreeCLocalResidual<TypeTag>
+{
+ typedef ThreePThreeCNILocalResidual<TypeTag> ThisType;
+ typedef ThreePThreeCLocalResidual<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+
+
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+
+ typedef typename GET_PROP_TYPE(TypeTag, ThreePThreeCIndices) Indices;
+
+ enum {
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld,
+
+ numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
+
+ energyEqIdx = Indices::energyEqIdx,
+ temperatureIdx = Indices::temperatureIdx,
+
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx,
+ gPhaseIdx = Indices::gPhaseIdx
+ };
+
+
+ typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+
+ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+ static constexpr Scalar massUpwindWeight = GET_PROP_VALUE(TypeTag, MassUpwindWeight);
+
+public:
+ /*!
+ * \brief Evaluate the amount all conservation quantities
+ * (e.g. phase mass) within a sub-control volume.
+ *
+ * The result should be averaged over the volume (e.g. phase mass
+ * inside a sub control volume divided by the volume)
+ *
+ * \param result The storage of the conservation quantitiy (mass or energy) within the sub-control volume
+ * \param scvIdx The SCV (sub-control-volume) index
+ * \param usePrevSol Evaluate function with solution of current or previous time step
+ */
+ void computeStorage(PrimaryVariables &result, int scvIdx, bool usePrevSol) const
+ {
+ // compute the storage term for phase mass
+ ParentType::computeStorage(result, scvIdx, usePrevSol);
+
+ // if flag usePrevSol is set, the solution from the previous
+ // time step is used, otherwise the current solution is
+ // used. The secondary variables are used accordingly. This
+ // is required to compute the derivative of the storage term
+ // using the implicit euler method.
+ const ElementVolumeVariables &elemDat = usePrevSol ? this->prevVolVars_() : this->curVolVars_();
+ const VolumeVariables &vertDat = elemDat[scvIdx];
+
+ // compute the energy storage
+ Scalar wdens = vertDat.density(wPhaseIdx);
+ Scalar ndens = vertDat.density(nPhaseIdx);
+ Scalar gdens = vertDat.density(gPhaseIdx);
+ Scalar wInEnerg = vertDat.internalEnergy(wPhaseIdx);
+ Scalar nInEnerg = vertDat.internalEnergy(nPhaseIdx);
+ Scalar gInEnerg = vertDat.internalEnergy(gPhaseIdx);
+ Scalar wsat = vertDat.saturation(wPhaseIdx);
+ Scalar nsat = vertDat.saturation(nPhaseIdx);
+ Scalar gsat = vertDat.saturation(gPhaseIdx);
+ Scalar temp = vertDat.temperature();
+ Scalar heatCap = vertDat.heatCapacity();
+ Scalar poro = vertDat.porosity();
+
+ result[energyEqIdx] = temp*heatCap
+ + poro * (gdens*gInEnerg*gsat
+ + wdens*wInEnerg*wsat
+ + ndens*nInEnerg*nsat);
+ /*
+ vertDat.porosity()*(vertDat.density(wPhaseIdx) *
+ vertDat.internalEnergy(wPhaseIdx) *
+ vertDat.saturation(wPhaseIdx)
+ +
+ vertDat.density(nPhaseIdx) *
+ vertDat.internalEnergy(nPhaseIdx) *
+ vertDat.saturation(nPhaseIdx)
+ +
+ vertDat.density(gPhaseIdx) *
+ vertDat.internalEnergy(gPhaseIdx) *
+ vertDat.saturation(gPhaseIdx))
+ +
+ vertDat.temperature()*vertDat.heatCapacity();
+ */
+ }
+
+ /*!
+ * \brief Evaluates the advective mass flux and the heat flux
+ * over a face of a subcontrol volume and writes the result in
+ * the flux vector.
+ *
+ * \param flux The advective flux over the SCV (sub-control-volume) face for each component
+ * \param fluxData The flux variables at the current SCV face
+ *
+ * This method is called by compute flux (base class)
+ */
+ void computeAdvectiveFlux(PrimaryVariables &flux,
+ const FluxVariables &fluxData) const
+ {
+ // advective mass flux
+ ParentType::computeAdvectiveFlux(flux, fluxData);
+
+ // advective heat flux in all phases
+ flux[energyEqIdx] = 0;
+ for (int phase = 0; phase < numPhases; ++phase) {
+ // vertex data of the upstream and the downstream vertices
+ const VolumeVariables &up = this->curVolVars_(fluxData.upstreamIdx(phase));
+ const VolumeVariables &dn = this->curVolVars_(fluxData.downstreamIdx(phase));
+
+ flux[energyEqIdx] +=
+ fluxData.KmvpNormal(phase) * (
+ massUpwindWeight * // upstream vertex
+ ( up.density(phase) *
+ up.mobility(phase) *
+ up.enthalpy(phase))
+ +
+ (1-massUpwindWeight) * // downstream vertex
+ ( dn.density(phase) *
+ dn.mobility(phase) *
+ dn.enthalpy(phase)) );
+ }
+ }
+
+ /*!
+ * \brief Adds the diffusive heat flux to the flux vector over
+ * the face of a sub-control volume.
+ *
+ * \param flux The diffusive flux over the SCV (sub-control-volume) face for each conservation quantity (mass, energy)
+ * \param fluxData The flux variables at the current SCV face
+ *
+ * This method is called by compute flux (base class)
+ */
+ void computeDiffusiveFlux(PrimaryVariables &flux,
+ const FluxVariables &fluxData) const
+ {
+ // diffusive mass flux
+ ParentType::computeDiffusiveFlux(flux, fluxData);
+
+ // diffusive heat flux
+ flux[temperatureIdx] +=
+ fluxData.normalMatrixHeatFlux();
+ }
+};
+
+}
+
+#endif
diff --git a/dumux/boxmodels/3p3cni/3p3cnimodel.hh b/dumux/boxmodels/3p3cni/3p3cnimodel.hh
new file mode 100644
index 0000000000000000000000000000000000000000..c3d57c228bc450c77ce69a7289dd5e0113074188
--- /dev/null
+++ b/dumux/boxmodels/3p3cni/3p3cnimodel.hh
@@ -0,0 +1,120 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * Copyright (C) 2011- by Holger Class *
+ * Copyright (C) 2008-2009 by Andreas Lauser *
+ * Copyright (C) 2008-2009 by Melanie Darcis *
+ * Copyright (C) 2008-2009 by Klaus Mosthaf *
+ * Copyright (C) 2008 by Bernd Flemisch *
+ * 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 Adaption of the BOX scheme to the non-isothermal three-phase three-component flow model.
+ */
+#ifndef DUMUX_NEW_3P3CNI_MODEL_HH
+#define DUMUX_NEW_3P3CNI_MODEL_HH
+
+#include <dumux/boxmodels/3p3c/3p3cmodel.hh>
+
+namespace Dumux {
+/*!
+ * \ingroup BoxModels
+ * \defgroup ThreePThreeCNIModel Non-isothermal three-phase three-component box model
+ */
+
+/*!
+ * \ingroup ThreePThreeCNIModel
+ * \brief Adaption of the BOX scheme to the non-isothermal three-phase three-component flow model.
+ *
+ * This model implements three-phase three-component flow of three fluid phases
+ * \f$\alpha \in \{ water, gas, NAPL \}\f$ each composed of up to three components
+ * \f$\kappa \in \{ water, air, contaminant \}\f$. The standard multiphase Darcy
+ * approach is used as the equation for the conservation of momentum:
+ * \f[
+ v_\alpha = - \frac{k_{r\alpha}}{\mu_\alpha} \mbox{\bf K}
+ \left(\text{grad}\, p_\alpha - \varrho_{\alpha} \mbox{\bf g} \right)
+ * \f]
+ *
+ * By inserting this into the equations for the conservation of the
+ * components, one transport equation for each component is obtained as
+ * \f{eqnarray*}
+ && \phi \frac{\partial (\sum_\alpha \varrho_{\text{mol}, \alpha} x_\alpha^\kappa
+ S_\alpha )}{\partial t}
+ - \sum\limits_\alpha \nabla \cdot \left\{ \frac{k_{r\alpha}}{\mu_\alpha}
+ \varrho_{\text{mol}, \alpha} x_\alpha^\kappa \mbox{\bf K}
+ (\nabla p_\alpha - \varrho_{\text{mass}, \alpha} \mbox{\bf g}) \right\}
+ \nonumber \\
+ \nonumber \\
+ && - \sum\limits_\alpha \nabla \cdot \left\{ D_{pm}^\kappa \varrho_{\text{mol},
+ \alpha } \nabla x_\alpha^\kappa \right\}
+ - q^\kappa = 0 \qquad \forall \kappa , \; \forall \alpha
+ \f}
+ *
+ * Note that these balance equations above are molar.
+ * In addition to that, a single balance of thermal energy is formulated
+ * for the fluid-filled porous medium under the assumption of local thermal
+ * equilibrium
+ * \f{eqnarray*}
+ && \phi \frac{\partial \left( \sum_\alpha \varrho_\alpha u_\alpha S_\alpha \right)}{\partial t}
+ + \left( 1 - \phi \right) \frac{\partial (\varrho_s c_s T)}{\partial t}
+ - \sum_\alpha \text{div} \left\{ \varrho_\alpha h_\alpha
+ \frac{k_{r\alpha}}{\mu_\alpha} \mathbf{K} \left( \text{grad}\,
+ p_\alpha
+ - \varrho_\alpha \mathbf{g} \right) \right\} \\
+ &-& \text{div} \left( \lambda_{pm} \text{grad} \, T \right)
+ - q^h = 0 \qquad \alpha \in \{w, n, g\}
+ \f}
+ *
+
+ *
+ * The equations are discretized using a fully-coupled vertex
+ * centered finite volume (BOX) scheme as spatial scheme and
+ * the implicit Euler method as temporal discretization.
+ *
+ * The model uses commonly applied auxiliary conditions like
+ * \f$S_w + S_n + S_g = 1\f$ for the saturations and
+ * \f$x^w_\alpha + x^a_\alpha + x^c_\alpha = 1\f$ for the mole fractions.
+ * Furthermore, the phase pressures are related to each other via
+ * capillary pressures between the fluid phases, which are functions of
+ * the saturation, e.g. according to the approach of Parker et al.
+ *
+ * The used primary variables are dependent on the locally present fluid phases
+ * An adaptive primary variable switch is included. The phase state is stored for all nodes
+ * of the system. The following cases can be distinguished:
+ * <ul>
+ * <li> All three phases are present: Primary variables are two saturations (\f$S_w\f$ and \f$S_n\f$, a pressure, in this case \f$p_g\f$, and the temperature \f$T\f$. </li>
+ * <li> Only the water phase is present: Primary variables are now the mole fractions of air and contaminant in the water phase (\f$x_w^a\f$ and \f$x_w^c\f$), as well as temperature and the gas pressure, which is, of course, in a case where only the water phase is present, just the same as the water pressure. </li>
+ * <li> Gas and NAPL phases are present: Primary variables (\f$S_n\f$, \f$x_g^w\f$, \f$p_g\f$, \f$T\f$). </li>
+ * <li> Water and NAPL phases are present: Primary variables (\f$S_n\f$, \f$x_w^a\f$, \f$p_g\f$, \f$T\f$). </li>
+ * <li> Only gas phase is present: Primary variables (\f$x_g^w\f$, \f$x_g^c\f$, \f$p_g\f$, \f$T\f$). </li>
+ * <li> Water and gas phases are present: Primary variables (\f$S_w\f$, \f$x_w^g\f$, \f$p_g\f$, \f$T\f$). </li>
+ * </ul>
+ *
+ */
+template<class TypeTag>
+class ThreePThreeCNIModel : public ThreePThreeCModel<TypeTag>
+{
+};
+
+}
+
+#include "3p3cnipropertydefaults.hh"
+
+#endif
diff --git a/dumux/boxmodels/3p3cni/3p3cniproblem.hh b/dumux/boxmodels/3p3cni/3p3cniproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..0dac26925c32f2c9d4f1dffaa07b151233426a60
--- /dev/null
+++ b/dumux/boxmodels/3p3cni/3p3cniproblem.hh
@@ -0,0 +1,81 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * Copyright (C) 2009 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 Base class for all problems which use the non-isothermal three-phase,
+ * three-component box model
+ */
+#ifndef DUMUX_3P3CNI_PROBLEM_HH
+#define DUMUX_3P3CNI_PROBLEM_HH
+
+#include <dumux/boxmodels/3p3c/3p3cproblem.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup ThreePThreeCNIModel
+ * \brief Base class for all problems which use the non-isothermal
+ * three-phase, three-component box model.
+ */
+template<class TypeTag>
+class ThreePThreeCNIProblem : public ThreePThreeCProblem<TypeTag>
+{
+ typedef ThreePThreeCProblem<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+
+public:
+ /*!
+ * \brief The constructor
+ *
+ * \param timeManager The time manager
+ * \param gridView The grid view
+ */
+ ThreePThreeCNIProblem(TimeManager &timeManager, const GridView &gridView)
+ : ParentType(timeManager, gridView)
+ {
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief Returns the temperature within the domain.
+ *
+ * This is a non-isothermal model, so this method just throws an
+ * exception. This method MUST NOT be overwritten by the actual
+ * problem.
+ */
+ Scalar temperature() const
+ { DUNE_THROW(Dune::Exception, "temperature() method called for a 3p3cni problem"); };
+
+ // \}
+};
+
+}
+
+#endif
diff --git a/dumux/boxmodels/3p3cni/3p3cniproperties.hh b/dumux/boxmodels/3p3cni/3p3cniproperties.hh
new file mode 100644
index 0000000000000000000000000000000000000000..b3c259fe17838380b69b19ddbfed5930f02110d0
--- /dev/null
+++ b/dumux/boxmodels/3p3cni/3p3cniproperties.hh
@@ -0,0 +1,56 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * Copyright (C) 2011 by Holger Class *
+ * Copyright (C) 2008 by Klaus Mosthaf, Andreas Lauser, Bernd Flemisch *
+ * 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/>. *
+ *****************************************************************************/
+/*!
+ * \ingroup ThreePThreeCNIModel
+ */
+/*!
+ * \file
+ *
+ * \brief Defines the properties required for the non-isothermal three-phase,
+ * three-component BOX model.
+ */
+#ifndef DUMUX_3P3CNI_PROPERTIES_HH
+#define DUMUX_3P3CNI_PROPERTIES_HH
+
+#include <dumux/boxmodels/3p3c/3p3cproperties.hh>
+
+namespace Dumux
+{
+
+namespace Properties
+{
+//////////////////////////////////////////////////////////////////
+// Type tags
+//////////////////////////////////////////////////////////////////
+
+//! The type tag for the non-isothermal three-phase, three-component problems
+NEW_TYPE_TAG(BoxThreePThreeCNI, INHERITS_FROM(BoxThreePThreeC));
+
+//////////////////////////////////////////////////////////////////
+// Property tags
+//////////////////////////////////////////////////////////////////
+NEW_PROP_TAG(ThreePThreeCNIIndices); //!< Enumerations for the 3p3cni models
+}
+}
+
+#endif
diff --git a/dumux/boxmodels/3p3cni/3p3cnipropertydefaults.hh b/dumux/boxmodels/3p3cni/3p3cnipropertydefaults.hh
new file mode 100644
index 0000000000000000000000000000000000000000..b4cc81f5fe594402e159b3371e091da97eb5c303
--- /dev/null
+++ b/dumux/boxmodels/3p3cni/3p3cnipropertydefaults.hh
@@ -0,0 +1,87 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * Copyright (C) 2011- by Holger Class *
+ * Copyright (C) 2008-2010 by Andreas Lauser *
+ * Copyright (C) 2008-2009 by Melanie Darcis *
+ * Copyright (C) 2008-2009 by Klaus Mosthaf *
+ * Copyright (C) 2008-2009 by Bernd Flemisch *
+ * 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/>. *
+ *****************************************************************************/
+/*!
+ * \ingroup ThreePThreeCNIModel
+ */
+/*!
+ * \file
+ *
+ * \brief Defines default values for most properties required by the 3p3cni
+ * box model.
+ */
+#ifndef DUMUX_3P3CNI_PROPERTY_DEFAULTS_HH
+#define DUMUX_3P3CNI_PROPERTY_DEFAULTS_HH
+
+#include <dumux/boxmodels/3p3c/3p3cpropertydefaults.hh>
+
+#include "3p3cnimodel.hh"
+#include "3p3cniproblem.hh"
+#include "3p3cniindices.hh"
+#include "3p3cnilocalresidual.hh"
+#include "3p3cnivolumevariables.hh"
+#include "3p3cnifluxvariables.hh"
+
+namespace Dumux
+{
+
+namespace Properties
+{
+//////////////////////////////////////////////////////////////////
+// Property values
+//////////////////////////////////////////////////////////////////
+
+SET_INT_PROP(BoxThreePThreeCNI, NumEq, 4); //!< set the number of equations to 4
+
+//! Use the 3p3cni local jacobian operator for the 3p3cni model
+SET_TYPE_PROP(BoxThreePThreeCNI,
+ LocalResidual,
+ ThreePThreeCNILocalResidual<TypeTag>);
+
+//! the Model property
+SET_TYPE_PROP(BoxThreePThreeCNI, Model, ThreePThreeCNIModel<TypeTag>);
+
+//! the VolumeVariables property
+SET_TYPE_PROP(BoxThreePThreeCNI, VolumeVariables, ThreePThreeCNIVolumeVariables<TypeTag>);
+
+
+//! the FluxVariables property
+SET_TYPE_PROP(BoxThreePThreeCNI, FluxVariables, ThreePThreeCNIFluxVariables<TypeTag>);
+
+//! The indices required by the non-isothermal 3p3c model
+SET_PROP(BoxThreePThreeCNI, ThreePThreeCIndices)
+{ typedef typename GET_PROP_TYPE(TypeTag, ThreePThreeCNIIndices) type; };
+
+SET_PROP(BoxThreePThreeCNI, ThreePThreeCNIIndices)
+{ private:
+ enum { formulation = GET_PROP_VALUE(TypeTag, Formulation) };
+ public:
+ typedef ThreePThreeCNIIndices<TypeTag, formulation, 0> type;
+};
+
+}
+
+}
+#endif
diff --git a/dumux/boxmodels/3p3cni/3p3cnivolumevariables.hh b/dumux/boxmodels/3p3cni/3p3cnivolumevariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..673e23b338a6c47145ea49ad163cebb7bf0c1718
--- /dev/null
+++ b/dumux/boxmodels/3p3cni/3p3cnivolumevariables.hh
@@ -0,0 +1,181 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * Copyright (C) 2011- by Holger Class *
+ * Copyright (C) 2008-2009 by Melanie Darcis *
+ * Copyright (C) 2008-2010 by Andreas Lauser *
+ * Copyright (C) 2008-2009 by Klaus Mosthaf *
+ * Copyright (C) 2008-2009 by Bernd Flemisch *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * *
+ * 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 Contains the quantities which are constant within a
+ * finite volume in the non-isothermal three-phase, three-component
+ * model.
+ */
+#ifndef DUMUX_3P3CNI_VOLUME_VARIABLES_HH
+#define DUMUX_3P3CNI_VOLUME_VARIABLES_HH
+
+#include <dumux/boxmodels/3p3c/3p3cvolumevariables.hh>
+
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup ThreePThreeCNIModel
+ * \brief Contains the quantities which are are constant within a
+ * finite volume in the non-isothermal three-phase, three-component
+ * model.
+ */
+template <class TypeTag>
+class ThreePThreeCNIVolumeVariables : public ThreePThreeCVolumeVariables<TypeTag>
+{
+ //! \cond 0
+ typedef ThreePThreeCVolumeVariables<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+
+ enum {
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld,
+ };
+
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, ThreePThreeCIndices) Indices;
+ enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
+ enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
+ enum { temperatureIdx = Indices::temperatureIdx };
+
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+ typedef Dune::FieldVector<Scalar, numPhases> PhasesVector;
+ //! \endcond
+
+public:
+ /*!
+ * \brief Update all quantities for a given control volume.
+ *
+ * \param sol The solution primary variables
+ * \param problem The problem
+ * \param element The element
+ * \param elemGeom Evaluate function with solution of current or previous time step
+ * \param vertIdx The local index of the SCV (sub-control volume)
+ * \param isOldSol Evaluate function with solution of current or previous time step
+ */
+ void update(const PrimaryVariables &sol,
+ const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int vertIdx,
+ bool isOldSol)
+ {
+ // vertex update data for the mass balance
+ ParentType::update(sol,
+ problem,
+ element,
+ elemGeom,
+ vertIdx,
+ isOldSol);
+
+ // the internal energies and the enthalpies
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ enthalpy_[phaseIdx] =
+ FluidSystem::enthalpy(phaseIdx, this->fluidState());
+ internalEnergy_[phaseIdx] =
+ FluidSystem::internalEnergy(phaseIdx, this->fluidState());
+ }
+ Valgrind::CheckDefined(internalEnergy_);
+ Valgrind::CheckDefined(enthalpy_);
+ };
+
+ /*!
+ * \brief Returns the total internal energy of a phase in the
+ * sub-control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar internalEnergy(int phaseIdx) const
+ { return internalEnergy_[phaseIdx]; };
+
+ /*!
+ * \brief Returns the total enthalpy of a phase in the sub-control
+ * volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar enthalpy(int phaseIdx) const
+ { return enthalpy_[phaseIdx]; };
+
+ /*!
+ * \brief Returns the total heat capacity \f$\mathrm{[J/(K*m^3]}\f$ of the rock matrix in
+ * the sub-control volume.
+ */
+ Scalar heatCapacity() const
+ { return heatCapacity_; };
+
+protected:
+ // this method gets called by the parent class. since this method
+ // is protected, we are friends with our parent..
+ friend class ThreePThreeCVolumeVariables<TypeTag>;
+
+ static Scalar temperature_(const PrimaryVariables &primaryVars,
+ const Problem& problem,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int scvIdx)
+ {
+ return primaryVars[temperatureIdx];
+ }
+
+ /*!
+ * \brief Update all quantities for a given control volume.
+ *
+ * \param sol The solution primary variables
+ * \param problem The problem
+ * \param element The element
+ * \param elemGeom Evaluate function with solution of current or previous time step
+ * \param scvIdx The local index of the SCV (sub-control volume)
+ * \param isOldSol Evaluate function with solution of current or previous time step
+ */
+ void updateEnergy_(const PrimaryVariables &sol,
+ const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int scvIdx,
+ bool isOldSol)
+ {
+ // copmute and set the heat capacity of the solid phase
+ heatCapacity_ = problem.spatialParameters().heatCapacity(element, elemGeom, scvIdx);
+ Valgrind::CheckDefined(heatCapacity_);
+ };
+
+
+ Scalar internalEnergy_[numPhases];
+ Scalar enthalpy_[numPhases];
+ Scalar heatCapacity_;
+};
+
+} // end namepace
+
+#endif
diff --git a/dumux/boxmodels/3p3cni/Makefile.am b/dumux/boxmodels/3p3cni/Makefile.am
new file mode 100644
index 0000000000000000000000000000000000000000..9045919660e66804c1a28d087e4dccfdc067b712
--- /dev/null
+++ b/dumux/boxmodels/3p3cni/Makefile.am
@@ -0,0 +1,4 @@
+3p3cnidir = $(includedir)/dumux/boxmodels/3p3cni
+3p3cni_HEADERS = *.hh
+
+include $(top_srcdir)/am/global-rules
diff --git a/dumux/boxmodels/Makefile.am b/dumux/boxmodels/Makefile.am
index 5f795772a5c222ffceeef8e8af13849a0e8a712c..b1c3d1c898c924b28f8f9a3fb7485d1b368423c2 100644
--- a/dumux/boxmodels/Makefile.am
+++ b/dumux/boxmodels/Makefile.am
@@ -1,4 +1,4 @@
-SUBDIRS = common 1p 1p2c 2p 2p2c 2p2cni 2pni MpNc richards
+SUBDIRS = common 1p 1p2c 2p 2p2c 2p2cni 2pni 3p3c 3p3cni MpNc richards
boxmodelsdir = $(includedir)/dumux/boxmodels
diff --git a/dumux/material/binarycoefficients/air_mesitylene.hh b/dumux/material/binarycoefficients/air_mesitylene.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7be22a9a0c1c26eb8a23a53ebb4b75a663a1b22d
--- /dev/null
+++ b/dumux/material/binarycoefficients/air_mesitylene.hh
@@ -0,0 +1,111 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * Copyright (C) 2011 by Holger Class *
+ * 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 Binary coefficients for air and mesitylene.
+ */
+#ifndef DUMUX_BINARY_COEFF_AIR_MESITYLENE_HH
+#define DUMUX_BINARY_COEFF_AIR_MESITYLENE_HH
+
+#include <dumux/material/components/air.hh>
+#include <dumux/material/components/mesitylene.hh>
+
+namespace Dumux
+{
+namespace BinaryCoeff
+{
+
+/*!
+ * \brief Binary coefficients for water and mesitylene.
+ */
+class Air_Mesitylene
+{
+public:
+ /*!
+ *
+ */
+ template <class Scalar>
+ static Scalar henry(Scalar temperature)
+ { DUNE_THROW(Dune::NotImplemented,
+ "Henry coefficient of air in mesitylene");
+ };
+
+ /*!
+ * \brief Binary diffusion coefficent [m^2/s] for air and mesitylene.
+ * I used the method according to Wilke and Lee
+ * see Handbook of chem. property's Estimation Methods
+ * W.J. Lyman, W.F. Reehl, D.H. Rosenblatt
+ *
+ */
+ template <class Scalar>
+ static Scalar gasDiffCoeff(Scalar temperature, Scalar pressure)
+ {
+ typedef Dumux::Air<Scalar> Air;
+ typedef Dumux::Mesitylene<Scalar> Mesitylene;
+
+ temperature = std::max(temperature, 1e-9); // regularization
+ temperature = std::min(temperature, 500.0); // regularization
+ pressure = std::max(pressure, 0.0); // regularization
+ pressure = std::min(pressure, 1e8); // regularization
+
+ const Scalar M_m = 1e3*Mesitylene::molarMass(); // [g/mol] molecular weight of mesitylene
+ const Scalar M_a = 1e3*Air::molarMass(); // [g/mol] molecular weight of air
+ const Scalar Tb_m = 437.9; // [K] boiling temperature of mesitylene
+ const Scalar sigma_a = 3.711; // charact. length of air
+ const Scalar T_scal_a = 78.6; // [K] (molec. energy of attraction/Boltzmann constant)
+ const Scalar V_B_m = 162.6; // [cm^3/mol] LeBas molal volume of mesitylene
+ const Scalar sigma_m = 1.18*std::pow(V_B_m, 0.333); // charact. length of mesitylene
+ const Scalar sigma_am = 0.5*(sigma_a + sigma_m);
+ const Scalar T_scal_m = 1.15*Tb_m;
+ const Scalar T_scal_am = std::sqrt(T_scal_a*T_scal_m);
+
+ Scalar T_star = temperature/T_scal_am;
+ T_star = std::max(T_star, 1e-5); // regularization
+
+ const Scalar Omega = 1.06036/std::pow(T_star, 0.1561) + 0.193/std::exp(T_star*0.47635)
+ + 1.03587/std::exp(T_star*1.52996) + 1.76474/std::exp(T_star*3.89411);
+ const Scalar B_ = 0.00217 - 0.0005*std::sqrt(1.0/M_a + 1.0/M_m);
+ const Scalar Mr = (M_a + M_m)/(M_a*M_m);
+ const Scalar D_am = (B_*std::pow(temperature, 1.5) * std::sqrt(Mr))
+ /(1e-5*pressure*std::pow(sigma_am, 2.0) * Omega); // [cm^2/s]
+
+ return 1e-4*D_am; // [m^2/s]
+ };
+
+ /*!
+ * \brief Diffusion coefficent [m^2/s] for molecular mesitylene in liquid water.
+ *
+ * \todo
+ */
+ template <class Scalar>
+ static Scalar liquidDiffCoeff(Scalar temperature, Scalar pressure)
+ { DUNE_THROW(Dune::NotImplemented,
+ "Binary liquid diffusion coefficients of air and mesitylene");
+ };
+};
+
+}
+} // end namepace
+
+#endif
diff --git a/dumux/material/binarycoefficients/air_xylene.hh b/dumux/material/binarycoefficients/air_xylene.hh
new file mode 100644
index 0000000000000000000000000000000000000000..e27ed3e66397df53922905d98d26b1b74ea5aac4
--- /dev/null
+++ b/dumux/material/binarycoefficients/air_xylene.hh
@@ -0,0 +1,111 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * Copyright (C) 2011 by Holger Class *
+ * 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 Binary coefficients for air and xylene.
+ */
+#ifndef DUMUX_BINARY_COEFF_AIR_XYLENE_HH
+#define DUMUX_BINARY_COEFF_AIR_XYLENE_HH
+
+#include <dumux/material/components/air.hh>
+#include <dumux/material/components/xylene.hh>
+
+namespace Dumux
+{
+namespace BinaryCoeff
+{
+
+/*!
+ * \brief Binary coefficients for water and xylene.
+ */
+class Air_Xylene
+{
+public:
+ /*!
+ *
+ */
+ template <class Scalar>
+ static Scalar henry(Scalar temperature)
+ { DUNE_THROW(Dune::NotImplemented,
+ "Henry coefficient of air in xylene");
+ };
+
+ /*!
+ * \brief Binary diffusion coefficent [m^2/s] for air and xylene.
+ * method according to Wilke and Lee
+ * see Handbook of chem. property's Estimation Methods
+ * W.J. Lyman, W.F. Reehl, D.H. Rosenblatt
+ *
+ */
+ template <class Scalar>
+ static Scalar gasDiffCoeff(Scalar temperature, Scalar pressure)
+ {
+ typedef Dumux::Air<Scalar> Air;
+ typedef Dumux::Xylene<Scalar> Xylene;
+
+ temperature = std::max(temperature, 1e-9); // regularization
+ temperature = std::min(temperature, 500.0); // regularization
+ pressure = std::max(pressure, 0.0); // regularization
+ pressure = std::min(pressure, 1e8); // regularization
+
+ const Scalar M_x = 1e3*Xylene::molarMass(); // [g/mol] molecular weight of xylene
+ const Scalar M_a = 1e3*Air::molarMass(); // [g/mol] molecular weight of air
+ const Scalar Tb_x = 412.0; // [K] boiling temperature of xylene
+ const Scalar sigma_a = 3.711; // charact. length of air
+ const Scalar T_scal_a = 78.6; // [K] (molec. energy of attraction/Boltzmann constant)
+ const Scalar V_B_x = 140.4; // [cm^3/mol] LeBas molal volume of xylene
+ const Scalar sigma_x = 1.18*std::pow(V_B_x, 0.333); // charact. length of xylene
+ const Scalar sigma_ax = 0.5*(sigma_a + sigma_x);
+ const Scalar T_scal_x = 1.15*Tb_x;
+ const Scalar T_scal_ax = std::sqrt(T_scal_a*T_scal_x);
+
+ Scalar T_star = temperature/T_scal_ax;
+ T_star = std::max(T_star, 1e-5); // regularization
+
+ const Scalar Omega = 1.06036/std::pow(T_star, 0.1561) + 0.193/std::exp(T_star*0.47635)
+ + 1.03587/std::exp(T_star*1.52996) + 1.76474/std::exp(T_star*3.89411);
+ const Scalar B_ = 0.00217 - 0.0005*std::sqrt(1.0/M_a + 1.0/M_x);
+ const Scalar Mr = (M_a + M_x)/(M_a*M_x);
+ const Scalar D_ax = (B_*std::pow(temperature,1.5)*std::sqrt(Mr))
+ /(1e-5*pressure*std::pow(sigma_ax, 2.0)*Omega); // [cm^2/s]
+
+ return D_ax*1e-4; // [m^2/s]
+ };
+
+ /*!
+ * \brief Diffusion coefficent [m^2/s] for molecular xylene in liquid water.
+ *
+ * \todo
+ */
+ template <class Scalar>
+ static Scalar liquidDiffCoeff(Scalar temperature, Scalar pressure)
+ { DUNE_THROW(Dune::NotImplemented,
+ "Binary liquid diffusion coefficients of air and xylene");
+ };
+};
+
+}
+} // end namepace
+
+#endif
diff --git a/dumux/material/binarycoefficients/h2o_mesitylene.hh b/dumux/material/binarycoefficients/h2o_mesitylene.hh
new file mode 100644
index 0000000000000000000000000000000000000000..a123cde5e1de0623ee5c61b553669f54bfa80876
--- /dev/null
+++ b/dumux/material/binarycoefficients/h2o_mesitylene.hh
@@ -0,0 +1,119 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * Copyright (C) 2011 by Holger Class *
+ * 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 Binary coefficients for water and mesitylene.
+ */
+#ifndef DUMUX_BINARY_COEFF_H2O_MESITYLENE_HH
+#define DUMUX_BINARY_COEFF_H2O_MESITYLENE_HH
+
+#include <dumux/material/components/h2o.hh>
+#include <dumux/material/components/mesitylene.hh>
+
+namespace Dumux
+{
+namespace BinaryCoeff
+{
+
+/*!
+ * \brief Binary coefficients for water and mesitylene.
+ */
+class H2O_Mesitylene
+{
+public:
+ /*!
+ * \brief Henry coefficent \f$[N/m^2]\f$ for mesitylene in liquid water.
+ *
+ * See:
+ *
+ * Sanders1999 Henry collection
+ */
+ template <class Scalar>
+ static Scalar henry(Scalar temperature)
+ {
+ // after Sanders
+ Scalar sanderH = 1.7e-1; // [M/atm]
+ //conversion to our Henry definition
+ Scalar dumuxH = sanderH / 101.325; // has now [(mol/m^3)/Pa]
+ dumuxH *= 18.02e-6; // multiplied by molar volume of reference phase = water
+ return 1.0/dumuxH; // [Pa]
+ };
+
+ /*!
+ * \brief Binary diffusion coefficent [m^2/s] for molecular water and mesitylene.
+ *
+ */
+ template <class Scalar>
+ static Scalar gasDiffCoeff(Scalar temperature, Scalar pressure)
+ {
+ typedef Dumux::H2O<Scalar> H2O;
+ typedef Dumux::Mesitylene<Scalar> Mesitylene;
+
+ temperature = std::max(temperature, 1e-9); // regularization
+ temperature = std::min(temperature, 500.0); // regularization
+ pressure = std::max(pressure, 0.0); // regularization
+ pressure = std::min(pressure, 1e8); // regularization
+
+ const Scalar M_m = 1e3*Mesitylene::molarMass(); // [g/mol] molecular weight of mesitylene
+ const Scalar M_w = 1e3*H2O::molarMass(); // [g/mol] molecular weight of water
+ const Scalar Tb_m = 437.9; // [K] boiling temperature of mesitylen
+ const Scalar Tb_w = 373.15; // [K] boiling temperature of water (at p_atm)
+ const Scalar V_B_w = 18.0; // [cm^3/mol] LeBas molal volume of water
+ const Scalar sigma_w = 1.18*std::pow(V_B_w, 0.333); // charact. length of air
+ const Scalar T_scal_w = 1.15*Tb_w; // [K] (molec. energy of attraction/Boltzmann constant)
+ const Scalar V_B_m = 162.6; // [cm^3/mol] LeBas molal volume of mesitylen
+ const Scalar sigma_m = 1.18*std::pow(V_B_m, 0.333); // charact. length of mesitylen
+ const Scalar sigma_wm = 0.5*(sigma_w + sigma_m);
+ const Scalar T_scal_m = 1.15*Tb_m;
+ const Scalar T_scal_wm = std::sqrt(T_scal_w*T_scal_m);
+
+ Scalar T_star = temperature/T_scal_wm;
+ T_star = std::max(T_star, 1e-5); // regularization
+
+ const Scalar Omega = 1.06036/std::pow(T_star,0.1561) + 0.193/std::exp(T_star*0.47635)
+ + 1.03587/std::exp(T_star*1.52996) + 1.76474/std::exp(T_star*3.89411);
+ const Scalar B_ = 0.00217 - 0.0005*std::sqrt(1.0/M_w + 1.0/M_m);
+ const Scalar Mr = (M_w + M_m)/(M_w*M_m);
+ const Scalar D_wm = (B_*std::pow(temperature,1.6)*std::sqrt(Mr))
+ /(1e-5*pressure*std::pow(sigma_wm, 2.0)*Omega); // [cm^2/s]
+
+ return D_wm*1e-4; // [m^2/s]
+ };
+
+ /*!
+ * \brief Diffusion coefficent [m^2/s] for mesitylene in liquid water.
+ *
+ * \todo
+ */
+ template <class Scalar>
+ static Scalar liquidDiffCoeff(Scalar temperature, Scalar pressure)
+ {
+ return 1.e-9; // This is just an order of magnitude. Please improve it!
+ };
+};
+
+}
+} // end namepace
+
+#endif
diff --git a/dumux/material/binarycoefficients/h2o_xylene.hh b/dumux/material/binarycoefficients/h2o_xylene.hh
new file mode 100644
index 0000000000000000000000000000000000000000..a6e94b2eaa38caa7168583bec7a791d310ecbd82
--- /dev/null
+++ b/dumux/material/binarycoefficients/h2o_xylene.hh
@@ -0,0 +1,120 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * Copyright (C) 2011 by Holger Class *
+ * 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 Binary coefficients for water and xylene.
+ */
+#ifndef DUMUX_BINARY_COEFF_H2O_XYLENE_HH
+#define DUMUX_BINARY_COEFF_H2O_XYLENE_HH
+
+#include <dumux/material/components/h2o.hh>
+#include <dumux/material/components/xylene.hh>
+
+namespace Dumux
+{
+namespace BinaryCoeff
+{
+
+/*!
+ * \brief Binary coefficients for water and xylene.
+ */
+class H2O_Xylene
+{
+public:
+ /*!
+ * \brief Henry coefficent \f$[N/m^2]\f$ for xylene in liquid water.
+ *
+ * See:
+ *
+ * Sanders1999 Henry collection
+ */
+
+ template <class Scalar>
+ static Scalar henry(Scalar temperature)
+ {
+ // after Sanders
+ Scalar sanderH = 1.5e-1; //[M/atm]
+ //conversion to our Henry definition
+ Scalar dumuxH = sanderH / 101.325; // has now [(mol/m^3)/Pa]
+ dumuxH *= 18.02e-6; //multiplied by molar volume of reference phase = water
+ return 1.0/dumuxH; // [Pa]
+ };
+
+ /*!
+ * \brief Binary diffusion coefficent [m^2/s] for molecular water and xylene.
+ *
+ */
+ template <class Scalar>
+ static Scalar gasDiffCoeff(Scalar temperature, Scalar pressure)
+ {
+ typedef Dumux::H2O<Scalar> H2O;
+ typedef Dumux::Xylene<Scalar> Xylene;
+
+ temperature = std::max(temperature, 1e-9); // regularization
+ temperature = std::min(temperature, 500.0); // regularization
+ pressure = std::max(pressure, 0.0); // regularization
+ pressure = std::min(pressure, 1e8); // regularization
+
+ const Scalar M_x = 1e3*Xylene::molarMass(); // [g/mol] molecular weight of xylene
+ const Scalar M_w = 1e3*H2O::molarMass(); // [g/mol] molecular weight of water
+ const Scalar Tb_x = 412.9; // [K] boiling temperature of xylene
+ const Scalar Tb_w = 373.15; // [K] boiling temperature of water (at p_atm)
+ const Scalar V_B_w = 18.0; // [cm^3/mol] LeBas molal volume of water
+ const Scalar sigma_w = 1.18*std::pow(V_B_w, 0.333); // charact. length of air
+ const Scalar T_scal_w = 1.15*Tb_w; // [K] (molec. energy of attraction/Boltzmann constant)
+ const Scalar V_B_x = 140.4; // [cm^3/mol] LeBas molal volume of xylene
+ const Scalar sigma_x = 1.18*std::pow(V_B_x, 0.333); // charact. length of xylene
+ const Scalar sigma_wx = 0.5*(sigma_w + sigma_x);
+ const Scalar T_scal_x = 1.15*Tb_x;
+ const Scalar T_scal_wx = std::sqrt(T_scal_w*T_scal_x);
+
+ Scalar T_star = temperature/T_scal_wx;
+ T_star = std::max(T_star, 1e-5); // regularization
+
+ const Scalar Omega = 1.06036/std::pow(T_star, 0.1561) + 0.193/std::exp(T_star*0.47635)
+ + 1.03587/std::exp(T_star*1.52996) + 1.76474/std::exp(T_star*3.89411);
+ const Scalar B_ = 0.00217 - 0.0005*std::sqrt(1.0/M_w + 1.0/M_x);
+ const Scalar Mr = (M_w + M_x)/(M_w*M_x);
+ const Scalar D_wx = (B_*std::pow(temperature,1.6)*std::sqrt(Mr))
+ /(1e-5*pressure*std::pow(sigma_wx, 2.0)*Omega); // [cm^2/s]
+
+ return D_wx*1e-4; // [m^2/s]
+ };
+
+ /*!
+ * \brief Diffusion coefficent [m^2/s] for xylene in liquid water.
+ *
+ * \todo
+ */
+ template <class Scalar>
+ static Scalar liquidDiffCoeff(Scalar temperature, Scalar pressure)
+ {
+ return 1.e-9; // This is just an order of magnitude. Please improve it!
+ };
+};
+
+}
+} // end namepace
+
+#endif
diff --git a/dumux/material/components/mesitylene.hh b/dumux/material/components/mesitylene.hh
new file mode 100644
index 0000000000000000000000000000000000000000..848a219828ec4db70a174afb424f819a81a1c15a
--- /dev/null
+++ b/dumux/material/components/mesitylene.hh
@@ -0,0 +1,280 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * Copyright (C) 2011 by Holger Class *
+ * Copyright (C) 2009 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/>. *
+ *****************************************************************************/
+#ifndef DUMUX_MESITYLENE_HH
+#define DUMUX_MESITYLENE_HH
+
+#include <dumux/material/components/component.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup Components
+ * \brief mesitylene
+ *
+ * \tparam Scalar The type used for scalar values
+ */
+template <class Scalar>
+class Mesitylene : public Component<Scalar, Mesitylene<Scalar> >
+{
+ typedef Component<Scalar, Mesitylene<Scalar> > ParentType;
+
+public:
+ /*!
+ * \brief A human readable name for the mesitylene
+ */
+ static const char *name()
+ { return "mesitylene"; }
+
+ /*!
+ * \brief The molar mass in \f$\mathrm{[kg/mol]}\f$ of mesitylene
+ */
+ static Scalar molarMass()
+ { return 0.120; }
+
+ /*!
+ * \brief Returns the critical temperature \f$\mathrm{[K]}\f$ of mesitylene
+ */
+ static Scalar criticalTemperature()
+ { return 637.3; }
+
+ /*!
+ * \brief Returns the critical pressure \f$\mathrm{[Pa]}\f$ of mesitylene
+ */
+ static Scalar criticalPressure()
+ { return 31.3e5; }
+
+ /*!
+ * \brief Returns the temperature \f$\mathrm{[K]}\f$ at mesitylene's triple point.
+ */
+ static Scalar tripleTemperature()
+ {
+ DUNE_THROW(Dune::NotImplemented, "tripleTemperature for mesitylene");
+ };
+
+ /*!
+ * \brief Returns the pressure \f$\mathrm{[Pa]}\f$ at mesitylene's triple point.
+ */
+ static Scalar triplePressure()
+ {
+ DUNE_THROW(Dune::NotImplemented, "triplePressure for mesitylene");
+ };
+
+ /*!
+ * \brief The saturation vapor pressure in \f$\mathrm{[Pa]}\f$ of pure mesitylene
+ * at a given temperature according to Antoine after Betz 1997 -> Gmehling et al 1980
+ *
+ * \param T temperature of component in \f$\mathrm{[K]}\f$
+ */
+ static Scalar vaporPressure(Scalar temperature)
+ {
+ const Scalar A = 7.07638;
+ const Scalar B = 1571.005;
+ const Scalar C = 209.728;
+
+ const Scalar T = temperature - 273.15;
+
+ return 100 * 1.334 * std::pow(10.0, (A - (B / (T + C))));
+ }
+
+
+ /*!
+ * \brief Specific enthalpy of liquid mesitylene \f$\mathrm{[J/kg]}\f$.
+ *
+ * \param temp temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar liquidEnthalpy(Scalar temp,
+ Scalar pressure)
+ {
+ const Scalar DTemp = temp - 273.0; // K -> degC
+ return 0.5*DTemp*(spHeatCapLiquidPhase_(0.2113*DTemp,pressure) + spHeatCapLiquidPhase_(0.7887*DTemp,pressure));
+ }
+
+ /*!
+ * \brief Specific enthalpy of mesitylene vapor \f$\mathrm{[J/kg]}\f$.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar gasEnthalpy(Scalar temperature, Scalar pressure)
+ {
+ return liquidEnthalpy(temperature,pressure) + vaporizationHeat_(temperature,pressure);
+ }
+
+ /*!
+ * \brief
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar gasDensity(Scalar temperature, Scalar pressure)
+ {
+ return IdealGas<Scalar>::density(molarMass(),
+ temperature,
+ pressure);
+ }
+
+ static Scalar molarGasDensity(Scalar temperature, Scalar pressure)
+ {
+ return (gasDensity(temperature, pressure)*molarMass());
+ }
+
+
+ /*!
+ * \brief The density of pure mesitylene at a given pressure and temperature \f$\mathrm{[kg/m^3]}\f$.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar liquidDensity(Scalar temperature, Scalar pressure)
+ {
+ return molarLiquidDensity_(temperature, pressure)*molarMass(); // [kg/m^3]
+ }
+
+ /*!
+ * \brief The dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of mesitylene vapor
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ * \param regularize defines, if the functions is regularized or not, set to true by default
+ */
+ static Scalar gasViscosity(Scalar temp, Scalar pressure, bool regularize=true)
+ {
+ temp = std::min(temp, 500.0); // regularization
+ temp = std::max(temp, 250.0);
+
+ // reduced temperature
+ Scalar Tr = temp/criticalTemperature();
+
+ Scalar Fp0 = 1.0;
+ Scalar xi = 0.00474;
+ Scalar eta_xi =
+ Fp0*(0.807*std::pow(Tr,0.618)
+ - 0.357*std::exp(-0.449*Tr)
+ + 0.34*std::exp(-4.058*Tr)
+ + 0.018);
+
+ return eta_xi/xi/1e7; // [Pa s]
+ }
+
+ /*!
+ * \brief The dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of pure mesitylene.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar liquidViscosity(Scalar temp, Scalar pressure)
+ {
+ temp = std::min(temp, 500.0); // regularization
+ temp = std::max(temp, 250.0);
+
+ const Scalar A = -6.749;
+ const Scalar B = 2010.0;
+
+ return std::exp(A + B/temp)*1e-3; // [Pa s]
+ }
+
+protected:
+ /*!
+ * \brief The molar density of pure mesitylene at a given pressure and temperature
+ * \f$\mathrm{[mol/m^3]}\f$.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar molarLiquidDensity_(Scalar temp, Scalar pressure)
+ {
+ temp = std::min(temp, 500.0); // regularization
+ temp = std::max(temp, 250.0);
+
+ const Scalar Z_RA = 0.2556; // from equation
+ const Scalar expo = 1.0 + std::pow(1.0 - temp/criticalTemperature(), 2.0/7.0);
+ Scalar V = (10.0*8.314*criticalTemperature()/criticalPressure()/1e5)*std::pow(Z_RA, expo); // liquid molar volume [cm^3/mol]
+ V *= 1e-6;
+
+ return 1.0/V; // molar density [mol/m^3]
+ }
+
+ /*!
+ * \brief latent heat of vaporization for mesitylene \f$\mathrm{[J/kg]}\f$.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar vaporizationHeat_(Scalar temp, Scalar pressure)
+ {
+ temp = std::min(temp, criticalTemperature()); // regularization
+ temp = std::max(temp, 0.0);
+
+ const Scalar DH_v_b = 39086.0; // [J/mol] Chen method (at boiling point 437.9 K */
+ // Variation with Temp according to Watson relation
+ const Scalar Tr1 = 0.687;
+ const Scalar Tr2 = temp/criticalTemperature();
+ const Scalar n = 0.375;
+ const Scalar DH_vap = DH_v_b * std::pow(((1.0 - Tr2)/(1.0 - Tr1)), n);
+
+ return (DH_vap/0.120); // we need [J/kg]
+ }
+
+ /*!
+ * \brief Specific heat cap of liquid mesitylene \f$\mathrm{[J/kg]}\f$.
+ *
+ * \param temp temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar spHeatCapLiquidPhase_(Scalar temp, Scalar pressure)
+ {
+ // according Reid et al. : Missenard group contrib. method (s. example 5-8)
+ // Mesitylen: C9H12 : 3* CH3 ; 1* C6H5 (phenyl-ring) ; -2* H (this was to much!)
+ // linear interpolation between table values [J/(mol K)]
+
+ Scalar H, CH3, C6H5;
+ if(temp < 298.0) {
+ // extrapolation for Temp<273
+ H = 13.4 + 1.2*(temp - 273.0)/25.0;
+ CH3 = 40.0 + 1.6*(temp - 273.0)/25.0;
+ C6H5 = 113.0 + 4.2*(temp - 273.0)/25.0;
+ }
+ else if(temp < 323.0){
+ H = 14.6 + 0.9*(temp - 298.0)/25.0;
+ CH3 = 41.6 + 1.9*(temp - 298.0)/25.0;
+ C6H5 = 117.2 + 6.2*(temp - 298.0)/25.0;
+ }
+ else if(temp < 348.0){
+ H = 15.5 + 1.2*(temp - 323.0)/25.0;
+ CH3 = 43.5 + 2.3*(temp - 323.0)/25.0;
+ C6H5 = 123.4 + 6.3*(temp - 323.0)/25.0;
+ }
+ else { // extrapolation for Temp>373
+ H = 16.7 + 2.1*(temp - 348.0)/25.0; // leads probably to underestimates
+ CH3 = 45.8 + 2.5*(temp - 348.0)/25.0;
+ C6H5 = 129.7 + 6.3*(temp - 348.0)/25.0;
+ }
+
+ return (C6H5 + 3*CH3 - 2*H)/molarMass();
+ }
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/material/components/xylene.hh b/dumux/material/components/xylene.hh
new file mode 100644
index 0000000000000000000000000000000000000000..adb248a77b25d5214e7bf9228fdf8a3387d2b2df
--- /dev/null
+++ b/dumux/material/components/xylene.hh
@@ -0,0 +1,293 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * Copyright (C) 2011 by Holger Class *
+ * Copyright (C) 2009 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/>. *
+ *****************************************************************************/
+#ifndef DUMUX_XYLENE_HH
+#define DUMUX_XYLENE_HH
+
+
+#include <dumux/material/components/component.hh>
+
+
+namespace Dumux
+{
+/*!
+ * \ingroup Components
+ * \brief xylene
+ *
+ * \tparam Scalar The type used for scalar values
+ */
+template <class Scalar>
+class Xylene : public Component<Scalar, Xylene<Scalar> >
+{
+ typedef Component<Scalar, Xylene<Scalar> > ParentType;
+
+public:
+ /*!
+ * \brief A human readable name for the xylene
+ */
+ static const char *name()
+ { return "xylene"; }
+
+ /*!
+ * \brief The molar mass in \f$\mathrm{[kg/mol]}\f$ of xylene
+ */
+ static Scalar molarMass()
+ { return 0.106; }
+
+ /*!
+ * \brief Returns the critical temperature \f$\mathrm{[K]}\f$ of xylene
+ */
+ static Scalar criticalTemperature()
+ { return 617.1; }
+
+ /*!
+ * \brief Returns the critical pressure \f$\mathrm{[Pa]}\f$ of xylene
+ */
+ static Scalar criticalPressure()
+ { return 35.4e5; }
+
+ /*!
+ * \brief Returns the temperature \f$\mathrm{[K]}\f$ at xylene's triple point.
+ */
+ static Scalar tripleTemperature()
+ {
+ DUNE_THROW(Dune::NotImplemented, "tripleTemperature for xylene");
+ };
+
+ /*!
+ * \brief Returns the pressure \f$\mathrm{[Pa]}\f$ at xylene's triple point.
+ */
+ static Scalar triplePressure()
+ {
+ DUNE_THROW(Dune::NotImplemented, "triplePressure for xylene");
+ };
+
+ /*!
+ * \brief The saturation vapor pressure in \f$\mathrm{[Pa]}\f$ of pure xylene
+ * at a given temperature according to Antoine after Betz 1997 -> Gmehling et al 1980
+ *
+ * \param T temperature of component in \f$\mathrm{[K]}\f$
+ */
+
+ static Scalar vaporPressure(Scalar temperature)
+ {
+ const Scalar A = 7.00909;;
+ const Scalar B = 1462.266;;
+ const Scalar C = 215.110;;
+
+ Scalar T = temperature - 273.15;
+ Scalar psat = 1.334*std::pow(10.0, (A - (B/(T + C)))); // in [mbar]
+ psat *= 100.0; // in [Pa] (0.001*1.E5)
+
+ return psat;
+ }
+
+ /*!
+ * \brief Specific heat cap of liquid xylene \f$\mathrm{[J/kg]}\f$.
+ *
+ * \param temp temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar spHeatCapLiquidPhase(Scalar temp, Scalar pressure)
+ {
+ Scalar CH3,C6H5,H;
+ // after Reid et al. : Missenard group contrib. method (s. example 5-8)
+ // Xylene: C9H12 : 3* CH3 ; 1* C6H5 (phenyl-ring) ; -2* H (this was too much!)
+ // linear interpolation between table values [J/(mol K)]
+
+ if(temp < 298.0){ // take care: extrapolation for Temp<273
+ H = 13.4 + 1.2*(temp - 273.0)/25.0;
+ CH3 = 40.0 + 1.6*(temp - 273.0)/25.0;
+ C6H5 = 113.0 + 4.2*(temp - 273.0)/25.0;
+ }
+ else if(temp < 323.0){
+ H = 14.6 + 0.9*(temp - 298.0)/25.0;
+ CH3 = 41.6 + 1.9*(temp - 298.0)/25.0;
+ C6H5 = 117.2 + 6.2*(temp - 298.0)/25.0;
+ }
+ else if(temp < 348.0){
+ H = 15.5 + 1.2*(temp - 323.0)/25.0;
+ CH3 = 43.5 + 2.3*(temp - 323.0)/25.0;
+ C6H5 = 123.4 + 6.3*(temp - 323.0)/25.0;
+ }
+ else { // take care: extrapolation for Temp>373
+ H = 16.7 + 2.1*(temp - 348.0)/25.0; // most likely leads to underestimation
+ CH3 = 45.8 + 2.5*(temp - 348.0)/25.0;
+ C6H5 = 129.7 + 6.3*(temp - 348.0)/25.0;
+ }
+
+ return (C6H5 + 2*CH3 - H)/molarMass();
+ }
+
+
+ /*!
+ * \brief Specific enthalpy of liquid xylene \f$\mathrm{[J/kg]}\f$.
+ *
+ * \param temp temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar liquidEnthalpy(Scalar temp,
+ Scalar pressure)
+ {
+ const Scalar T_ref = 273.0; // typically T_ref=0 C, but not necessarily
+ const Scalar DTemp = temp - T_ref;
+
+ // numerical integration using 2-point Gauss */
+ return 0.5*DTemp*(spHeatCapLiquidPhase(0.2113*DTemp,pressure)
+ + spHeatCapLiquidPhase(0.7887*DTemp,pressure));
+ }
+
+ /*!
+ * \brief latent heat of vaporization for xylene \f$\mathrm{[J/kg]}\f$.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar heatVap(Scalar temp, Scalar pressure)
+ {
+ // Variation with Temp according to Watson relation
+ temp = std::min(temp, criticalTemperature()); // regularization
+ temp = std::max(temp, 0.0); // regularization
+
+ const Scalar DH_v_b = 36195.0; // [J/mol] Chen method (at boiling point 437.9 K
+ const Scalar Tr1 = 0.668; // Tb/T_crit
+ const Scalar Tr2 = temp/criticalTemperature();
+ const Scalar n = 0.375;
+ const Scalar DH_vap = DH_v_b * std::pow(((1.0 - Tr2)/(1.0 - Tr1)), n);
+
+ return (DH_vap/molarMass()); // we need [J/kg]
+ }
+
+ /*!
+ * \brief Specific enthalpy of xylene vapor \f$\mathrm{[J/kg]}\f$.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar gasEnthalpy(Scalar temperature, Scalar pressure)
+ {
+ return liquidEnthalpy(temperature, pressure) + heatVap(temperature, pressure);
+ }
+
+ /*!
+ * \brief
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar gasDensity(Scalar temperature, Scalar pressure)
+ {
+ return IdealGas<Scalar>::density(molarMass(),
+ temperature,
+ pressure);
+ }
+
+ static Scalar molarGasDensity(Scalar temperature, Scalar pressure)
+ {
+ return (gasDensity(temperature, pressure)*molarMass());
+ }
+
+ /*!
+ * \brief The molar density of pure xylene at a given pressure and temperature
+ * \f$\mathrm{[mol/m^3]}\f$.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar molarLiquidDensity(Scalar temp, Scalar pressure)
+ {
+ // saturated molar volume according to Lide, CRC Handbook of
+ // Thermophysical and Thermochemical Data, CRC Press, 1994
+ // valid for 245 < Temp < 600
+ temp = std::min(temp, 500.0); // regularization
+ temp = std::max(temp, 250.0); // regularization
+
+ const Scalar A1 = 0.25919; // from table
+ const Scalar A2 = 0.0014569; // from table
+ const Scalar expo = 1.0 + std::pow((1.0 - temp/criticalTemperature()), (2.0/7.0));
+ const Scalar V = A2*std::pow(A1, expo); // liquid molar volume [m^3/mol]
+
+ return 1.0/V; // molar density [mol/m^3]
+ }
+
+ /*!
+ * \brief The density of pure xylene at a given pressure and temperature \f$\mathrm{[kg/m^3]}\f$.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar liquidDensity(Scalar temperature, Scalar pressure)
+ {
+ return (molarLiquidDensity(temperature, pressure)*molarMass()); // [kg/m^3]
+ }
+
+ /*!
+ * \brief The dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of xylene vapor
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ * \param regularize defines, if the functions is regularized or not, set to true by default
+ */
+ static Scalar gasViscosity(Scalar temp, Scalar pressure, bool regularize=true)
+ {
+ temp = std::min(temp, 500.0); // regularization
+ temp = std::max(temp, 250.0); // regularization
+
+ const Scalar Tr = std::max(temp/criticalTemperature(), 1e-10);
+ const Scalar Fp0 = 1.0;
+ const Scalar xi = 0.004623;
+ const Scalar eta_xi = Fp0*(0.807*std::pow(Tr, 0.618)
+ - 0.357*std::exp(-0.449*Tr)
+ + 0.34*std::exp(-4.058*Tr)
+ + 0.018);
+ Scalar r = eta_xi/xi; // [1e-6 P]
+ r /= 1.0e7; // [Pa s]
+
+ return r;
+ }
+
+ /*!
+ * \brief The dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of pure xylene.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar liquidViscosity(Scalar temp, Scalar pressure)
+ {
+ temp = std::min(temp, 500.0); // regularization
+ temp = std::max(temp, 250.0); // regularization
+
+ const Scalar A = -3.82;
+ const Scalar B = 1027.0;
+ const Scalar C = -6.38e-4;
+ const Scalar D = 4.52e-7;
+
+ Scalar r = std::exp(A + B/temp + C*temp + D*temp*temp); // in [cP]
+ r *= 1.0e-3; // in [Pa s]
+
+ return r; // [Pa s]
+ }
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/material/fluidmatrixinteractions/3p/parkerVanGen3p.hh b/dumux/material/fluidmatrixinteractions/3p/parkerVanGen3p.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d64ed5d41f1f43ae271d1320a5d2e114834b4284
--- /dev/null
+++ b/dumux/material/fluidmatrixinteractions/3p/parkerVanGen3p.hh
@@ -0,0 +1,381 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * Copyright (C) 2011 by Holger Class *
+ * Copyright (C) 2010 by Jochen Fritz, Benjamin Faigle *
+ * 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 VanGenuchten.hh Implementation of van Genuchten's capillary
+ * pressure <-> saturation relation
+ */
+#ifndef PARKERVANGEN_3P_HH
+#define PARKERVANGEN_3P_HH
+
+
+#include "parkerVanGen3pparams.hh"
+
+#include <algorithm>
+
+
+namespace Dumux
+{
+/*!
+ * \ingroup material
+ *
+ * \brief Implementation of van Genuchten's capillary pressure <->
+ * saturation relation. This class bundles the "raw" curves
+ * as static members and doesn't concern itself converting
+ * absolute to effective saturations and vince versa.
+ *
+ * \sa VanGenuchten, VanGenuchtenThreephase
+ */
+template <class ScalarT, class ParamsT = ParkerVanGen3PParams<ScalarT> >
+class ParkerVanGen3P
+{
+ static const int wPhaseIdx = 0;
+ static const int nPhaseIdx = 1;
+ static const int gPhaseIdx = 2;
+
+public:
+ typedef ParamsT Params;
+ typedef typename Params::Scalar Scalar;
+
+ /*!
+ * \brief The capillary pressure-saturation curve.
+ *
+ */
+ static Scalar pC(const Params ¶ms, Scalar Sw)
+ {
+ DUNE_THROW(Dune::NotImplemented, "Capillary pressures for three phases is not so simple! Use pCGN, pCNW, and pcGW");
+ }
+
+ static Scalar pCGW(const Params ¶ms, Scalar Sw)
+ {
+ /*
+ Sw = wetting phase saturation, or,
+ sum of wetting phase saturations
+ alpha : VanGenuchten-alpha
+ this function is just copied from MUFTE/pml/constrel3p3cni.c
+ that is why variable names do not yet fulfill Dumux rules, TODO Change */
+
+ Scalar r,Se,x,vg_m;
+ Scalar pc,pc_prime,Se_regu;
+ Scalar PC_VG_REG = 0.01;
+
+ Se = (Sw-params.Swr())/(1.-params.Sgr());
+
+ /* Snr = 0.0; test version */
+
+ /* regularization */
+ if (Se<0.0) Se=0.0;
+ if (Se>1.0) Se=1.0;
+ vg_m = 1.-1./params.vgN();
+
+ if (Se>PC_VG_REG && Se<1-PC_VG_REG)
+ {
+ r = std::pow(Se,-1/vg_m);
+ x = r-1;
+ vg_m = 1-vg_m;
+ x = std::pow(x,vg_m);
+ r = x/params.vgAlpha();
+ return(r);
+ }
+ else
+ {
+ /* value and derivative at regularization point */
+ if (Se<=PC_VG_REG) Se_regu = PC_VG_REG; else Se_regu = 1-PC_VG_REG;
+ pc = std::pow(std::pow(Se_regu,-1/vg_m)-1,1/params.vgN())/params.vgAlpha();
+ pc_prime = std::pow(std::pow(Se_regu,-1/vg_m)-1,1/params.vgN()-1)*std::pow(Se_regu,-1/vg_m-1)*(-1/vg_m)/params.vgAlpha()/(1-params.Sgr()-params.Swr())/params.vgN();
+
+ /* evaluate tangential */
+ r = (Se-Se_regu)*pc_prime+pc;
+ return(r);
+ }
+ }
+
+ static Scalar pCNW(const Params ¶ms, Scalar Sw)
+ {
+ /*
+ Sw = wetting phase saturation, or,
+ sum of wetting phase saturations
+ alpha : VanGenuchten-alpha
+ this function is just copied from MUFTE/pml/constrel3p3cni.c
+ that is why variable names do not yet fulfill Dumux rules, TODO Change */
+
+ Scalar r,Se,x,vg_m;
+ Scalar pc,pc_prime,Se_regu;
+ Scalar PC_VG_REG = 0.01;
+
+ Se = (Sw-params.Swr())/(1.-params.Snr());
+
+ /* Snr = 0.0; test version */
+
+ /* regularization */
+ if (Se<0.0) Se=0.0;
+ if (Se>1.0) Se=1.0;
+ vg_m = 1.-1./params.vgN();
+
+ if (Se>PC_VG_REG && Se<1-PC_VG_REG)
+ {
+ r = std::pow(Se,-1/vg_m);
+ x = r-1;
+ vg_m = 1-vg_m;
+ x = std::pow(x,vg_m);
+ r = x/params.vgAlpha();
+ return(r);
+ }
+ else
+ {
+ /* value and derivative at regularization point */
+ if (Se<=PC_VG_REG) Se_regu = PC_VG_REG; else Se_regu = 1-PC_VG_REG;
+ pc = std::pow(std::pow(Se_regu,-1/vg_m)-1,1/params.vgN())/params.vgAlpha();
+ pc_prime = std::pow(std::pow(Se_regu,-1/vg_m)-1,1/params.vgN()-1)*std::pow(Se_regu,-1/vg_m-1)*(-1/vg_m)/params.vgAlpha()/(1-params.Snr()-params.Swr())/params.vgN();
+
+ /* evaluate tangential */
+ r = (Se-Se_regu)*pc_prime+pc;
+ return(r);
+ }
+ }
+
+ static Scalar pCGN(const Params ¶ms, Scalar St)
+ {
+ /*
+ St = sum of wetting (liquid) phase saturations
+ alpha : VanGenuchten-alpha
+ this function is just copied from MUFTE/pml/constrel3p3cni.c
+ that is why variable names do not yet fulfill Dumux rules, TODO Change */
+
+ Scalar r,Se,x,vg_m;
+ Scalar pc,pc_prime,Se_regu;
+ Scalar PC_VG_REG = 0.01;
+
+ Se = (St-params.Swrx())/(1.-params.Swrx());
+
+ /* Snr = 0.0; test version */
+
+ /* regularization */
+ if (Se<0.0) Se=0.0;
+ if (Se>1.0) Se=1.0;
+ vg_m = 1.-1./params.vgN();
+
+ if (Se>PC_VG_REG && Se<1-PC_VG_REG)
+ {
+ r = std::pow(Se,-1/vg_m);
+ x = r-1;
+ vg_m = 1-vg_m;
+ x = std::pow(x,vg_m);
+ r = x/params.vgAlpha();
+ return(r);
+ }
+ else
+ {
+ /* value and derivative at regularization point */
+ if (Se<=PC_VG_REG) Se_regu = PC_VG_REG; else Se_regu = 1-PC_VG_REG;
+ pc = std::pow(std::pow(Se_regu,-1/vg_m)-1,1/params.vgN())/params.vgAlpha();
+ pc_prime = std::pow(std::pow(Se_regu,-1/vg_m)-1,1/params.vgN()-1)*std::pow(Se_regu,-1/vg_m-1)*(-1/vg_m)/params.vgAlpha()/(1-params.Sgr()-params.Swrx())/params.vgN();
+
+ /* evaluate tangential */
+ r = (Se-Se_regu)*pc_prime+pc;
+ return(r);
+ }
+ }
+
+ static Scalar pCAlpha(const Params ¶ms, Scalar Sn)
+ {
+ /* continuous transition to zero */
+ Scalar alpha,Sne;
+
+ Sne=Sn;
+ /* regularization */
+ if (Sne<=0.001) Sne=0.0;
+ if (Sne>=1.0) Sne=1.0;
+
+ if (Sne>params.Snr()) alpha = 1.0;
+ else
+ {
+ if (params.Snr()>=0.001) alpha = Sne/params.Snr();
+ else alpha = 0.0;
+ }
+ return(alpha);
+ }
+
+ /*!
+ * \brief The saturation-capillary pressure curve.
+ *
+ */
+ static Scalar Sw(const Params ¶ms, Scalar pC)
+ {
+ DUNE_THROW(Dune::NotImplemented, "Sw(pc) for three phases not implemented! Do it yourself!");
+ }
+
+ /*!
+ * \brief Returns the partial derivative of the capillary
+ * pressure to the effective saturation.
+ *
+ */
+ static Scalar dpC_dSw(const Params ¶ms, Scalar Sw)
+ {
+ DUNE_THROW(Dune::NotImplemented, "dpC/dSw for three phases not implemented! Do it yourself!");
+ }
+
+ /*!
+ * \brief Returns the partial derivative of the effective
+ * saturation to the capillary pressure.
+ */
+ static Scalar dSw_dpC(const Params ¶ms, Scalar pC)
+ {
+ DUNE_THROW(Dune::NotImplemented, "dSw/dpC for three phases not implemented! Do it yourself!");
+ }
+
+ /*!
+ * \brief The relative permeability for the wetting phase of
+ * the medium implied by van Genuchten's
+ * parameterization.
+ *
+ * The permeability of water in a 3p system equals the standard 2p description.
+ * (see p61. in "Comparison of the Three-Phase Oil Relative Permeability Models"
+ * MOJDEH DELSHAD and GARY A. POPE, Transport in Porous Media 4 (1989), 59-83.)
+ *
+ * \param The mobile saturation of all phases. (Sw used)
+ */
+ static Scalar krw(const Params ¶ms, Scalar saturation, Scalar Sn, Scalar Sg)
+ {
+
+ //transformation to effective saturation
+ Scalar Se = (saturation - params.Swr()) / (1-params.Swr());
+
+ /* regularization */
+ if(Se > 1.0) return 1.;
+ if(Se < 0.0) return 0.;
+
+ Scalar r = 1. - std::pow(1 - std::pow(Se, 1/params.vgM()), params.vgM());
+ return std::sqrt(Se)*r*r;
+ };
+
+ /*!
+ * \brief The relative permeability for the non-wetting phase
+ * after the Model of Parker et al. (1987).
+ *
+ * See model 7 in "Comparison of the Three-Phase Oil Relative Permeability Models"
+ * MOJDEH DELSHAD and GARY A. POPE, Transport in Porous Media 4 (1989), 59-83.
+ * or more comprehensive in
+ * "Estimation of primary drainage three-phase relative permeability for organic
+ * liquid transport in the vadose zone", Leonardo I. Oliveira, Avery H. Demond,
+ * Journal of Contaminant Hydrology 66 (2003), 261-285
+ *
+ *
+ * \param The mobile saturation of all phases. (Sw and Sn used)
+ */
+ static Scalar krn(const Params ¶ms, Scalar Sw, Scalar saturation, Scalar Sg)
+ {
+
+ Scalar Swe = std::min((Sw - params.Swr()) / (1 - params.Swr()), 1.);
+ Scalar Ste = std::min((Sw + saturation - params.Swr()) / (1 - params.Swr()), 1.);
+
+ // regularization
+ if(Swe <= 0.0) Swe = 0.;
+ if(Ste <= 0.0) Ste = 0.;
+ if(Ste - Swe <= 0.0) return 0.;
+
+ Scalar krn_;
+ krn_ = std::pow(1 - std::pow(Swe, 1/params.vgM()), params.vgM());
+ krn_ -= std::pow(1 - std::pow(Ste, 1/params.vgM()), params.vgM());
+ krn_ *= krn_;
+
+ if (params.krRegardsSnr())
+ {
+ // regard Snr in the permeability of the n-phase, see Helmig1997
+ Scalar resIncluded = std::max(std::min((saturation - params.Snr()/ (1-params.Swr())), 1.), 0.);
+ krn_ *= std::sqrt(resIncluded );
+ }
+ else
+ krn_ *= std::sqrt(saturation / (1 - params.Swr())); // Hint: (Ste - Swe) = Sn / (1-Srw)
+
+
+ return krn_;
+ };
+
+
+ /*!
+ * \brief The relative permeability for the non-wetting phase
+ * of the medium implied by van Genuchten's
+ * parameterization.
+ *
+ * The permeability of gas in a 3p system equals the standard 2p description.
+ * (see p61. in "Comparison of the Three-Phase Oil Relative Permeability Models"
+ * MOJDEH DELSHAD and GARY A. POPE, Transport in Porous Media 4 (1989), 59-83.)
+ *
+ * \param The mobile saturation of all phases. (Sg used)
+ */
+ static Scalar krg(const Params ¶ms, Scalar Sw, Scalar Sn, Scalar saturation)
+ {
+
+ // Se = (Sw+Sn - Sgr)/(1-Sgr)
+ Scalar Se = std::min(((1-saturation) - params.Sgr()) / (1 - params.Sgr()), 1.);
+
+
+ /* regularization */
+ if(Se > 1.0) return 0.0;
+ if(Se < 0.0) return 1.0;
+ Scalar scalFact = 1.;
+ if (saturation<=0.1)
+ {
+ scalFact = (saturation - params.Sgr())/(0.1 - params.Sgr());
+ if (scalFact < 0.) scalFact = 0.;
+ }
+
+ Scalar result = scalFact * std::pow(1 - Se, 1.0/3.) * std::pow(1 - std::pow(Se, 1/params.vgM()), 2*params.vgM());
+
+ return result;
+ };
+
+ /*!
+ * \brief The relative permeability for a phase.
+ *
+ * \param Phase indicator, The saturation of all phases.
+ */
+ static Scalar kr(const Params ¶ms, const int phase, const Scalar Sw, const Scalar Sn, const Scalar Sg)
+ {
+ switch (phase)
+ {
+ case 0:
+ return krw(params, Sw, Sn, Sg);
+ break;
+ case 1:
+ return krn(params, Sw, Sn, Sg);
+ break;
+ case 2:
+ return krg(params, Sw, Sn, Sg);
+ break;
+ }
+ return 0;
+ };
+
+ /*
+ * \brief the basis for calculating adsorbed NAPL in storage term
+ * \param bulk density of porous medium, adsorption coefficient
+ */
+ static Scalar bulkDensTimesAdsorpCoeff (const Params ¶ms)
+ {
+ return params.rhoBulk() * params.KdNAPL();
+ }
+};
+}
+
+#endif
diff --git a/dumux/material/fluidmatrixinteractions/3p/parkerVanGen3pparams.hh b/dumux/material/fluidmatrixinteractions/3p/parkerVanGen3pparams.hh
new file mode 100644
index 0000000000000000000000000000000000000000..af25bad9453d7371aaaa8cc11a4805e30b5b3929
--- /dev/null
+++ b/dumux/material/fluidmatrixinteractions/3p/parkerVanGen3pparams.hh
@@ -0,0 +1,233 @@
+/*****************************************************************************
+ * Copyright (C) 2011 by Holger Class *
+ * Copyright (C) 2010 by Jochen Fritz, Benjamin Faigle *
+ * 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 vangenuchtenparams.hh Specification of the material params
+ * for the van Genuchten capillary pressure model.
+ *
+ * In comparison to the 2p version, this parameter container also includes
+ * the residual saturations, as their inclusion is very model-specific.
+ */
+#ifndef PARKERVANGEN_PARAMS_3P_HH
+#define PARKERVANGEN_PARAMS_3P_HH
+
+namespace Dumux
+{
+/*!
+ * \brief Reference implementation of a van Genuchten params
+ */
+template<class ScalarT>
+class ParkerVanGen3PParams
+{
+public:
+ typedef ScalarT Scalar;
+
+ ParkerVanGen3PParams()
+ {}
+
+ ParkerVanGen3PParams(Scalar vgAlpha, Scalar vgN, Scalar KdNAPL, Scalar rhoBulk, Dune::FieldVector<Scalar, 4> residualSaturation, bool regardSnr=false)
+ {
+ setVgAlpha(vgAlpha);
+ setVgN(vgN);
+ setSwr(residualSaturation[0]);
+ setSnr(residualSaturation[1]);
+ setSgr(residualSaturation[2]);
+ setSwrx(residualSaturation[3]);
+ setkrRegardsSnr(regardSnr);
+ setKdNAPL(KdNAPL);
+ setRhoBulk(rhoBulk);
+ };
+
+ /*!
+ * \brief Return the \f$\alpha\f$ shape parameter of van Genuchten's
+ * curve.
+ */
+ Scalar vgAlpha() const
+ { return vgAlpha_; }
+
+ /*!
+ * \brief Set the \f$\alpha\f$ shape parameter of van Genuchten's
+ * curve.
+ */
+ void setVgAlpha(Scalar v)
+ { vgAlpha_ = v; }
+
+ /*!
+ * \brief Return the \f$m\f$ shape parameter of van Genuchten's
+ * curve.
+ */
+ Scalar vgM() const
+ { return vgM_; }
+
+ /*!
+ * \brief Set the \f$m\f$ shape parameter of van Genuchten's
+ * curve.
+ *
+ * The \f$n\f$ shape parameter is set to \f$n = \frac{1}{1 - m}\f$
+ */
+ void setVgM(Scalar m)
+ { vgM_ = m; vgN_ = 1/(1 - vgM_); }
+
+ /*!
+ * \brief Return the \f$n\f$ shape parameter of van Genuchten's
+ * curve.
+ */
+ Scalar vgN() const
+ { return vgN_; }
+
+ /*!
+ * \brief Set the \f$n\f$ shape parameter of van Genuchten's
+ * curve.
+ *
+ * The \f$n\f$ shape parameter is set to \f$m = 1 - \frac{1}{n}\f$
+ */
+ void setVgN(Scalar n)
+ { vgN_ = n; vgM_ = 1 - 1/vgN_; }
+
+ /*!
+ * \brief Return the residual saturation.
+ */
+ Scalar satResidual(int phaseIdx) const
+ {
+ switch (phaseIdx)
+ {
+ case 0:
+ return Swr_;
+ break;
+ case 1:
+ return Snr_;
+ break;
+ case 2:
+ return Sgr_;
+ break;
+ };
+ DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
+ }
+
+ /*!
+ * \brief Set all residual saturations.
+ */
+ void setResiduals(Dune::FieldVector<Scalar, 3> residualSaturation)
+ {
+ setSwr(residualSaturation[0]);
+ setSnr(residualSaturation[1]);
+ setSgr(residualSaturation[2]);
+ }
+
+
+ /*!
+ * \brief Return the residual wetting saturation.
+ */
+ Scalar Swr() const
+ { return Swr_; }
+
+ /*!
+ * \brief Set the residual wetting saturation.
+ */
+ void setSwr(Scalar input)
+ { Swr_ = input; }
+
+ /*!
+ * \brief Return the residual non-wetting saturation.
+ */
+ Scalar Snr() const
+ { return Snr_; }
+
+ /*!
+ * \brief Set the residual non-wetting saturation.
+ */
+ void setSnr(Scalar input)
+ { Snr_ = input; }
+
+ /*!
+ * \brief Return the residual gas saturation.
+ */
+ Scalar Sgr() const
+ { return Sgr_; }
+
+ /*!
+ * \brief Set the residual gas saturation.
+ */
+ void setSgr(Scalar input)
+ { Sgr_ = input; }
+
+ Scalar Swrx() const
+ { return Swrx_; }
+
+ /*!
+ * \brief Set the residual gas saturation.
+ */
+ void setSwrx(Scalar input)
+ { Swrx_ = input; }
+
+ /*!
+ * \brief defines if residual n-phase saturation should be regarded in its relative permeability.
+ */
+ void setkrRegardsSnr(bool input)
+ { krRegardsSnr_ = input; }
+ /*!
+ * \brief Calls if residual n-phase saturation should be regarded in its relative permeability.
+ */
+ bool krRegardsSnr() const
+ { return krRegardsSnr_; }
+
+
+ /*!
+ * \brief Return the bulk density of the porous medium
+ */
+ Scalar rhoBulk() const
+ { return rhoBulk_; }
+
+ /*!
+ * \brief Set the bulk density of the porous medium
+ */
+ void setRhoBulk(Scalar input)
+ { rhoBulk_ = input; }
+
+ /*!
+ * \brief Return the adsorption coefficient
+ */
+ Scalar KdNAPL() const
+ { return KdNAPL_; }
+
+ /*!
+ * \brief Set the adsorption coefficient
+ */
+ void setKdNAPL(Scalar input)
+ { KdNAPL_ = input; }
+
+
+private:
+ Scalar vgAlpha_;
+ Scalar vgM_;
+ Scalar vgN_;
+ Scalar Swr_;
+ Scalar Snr_;
+ Scalar Sgr_;
+ Scalar Swrx_; /* (Sw+Sn)_r */
+
+ Scalar KdNAPL_;
+ Scalar rhoBulk_;
+
+ bool krRegardsSnr_ ;
+};
+} // namespace Dumux
+
+#endif
diff --git a/dumux/material/fluidmatrixinteractions/Makefile.am b/dumux/material/fluidmatrixinteractions/Makefile.am
index 1e1fef223fd8ee5915237290d11ba55b31b6615d..355e5bf5c62f6295cb0d550d4eb52393e38ba24c 100644
--- a/dumux/material/fluidmatrixinteractions/Makefile.am
+++ b/dumux/material/fluidmatrixinteractions/Makefile.am
@@ -1,4 +1,4 @@
-SUBDIRS = 2p Mp
+SUBDIRS = 2p 3p Mp
fluidmatrixinteractionsdir = $(includedir)/dumux/material/fluidmatrixinteractions
diff --git a/dumux/material/fluidsystems/h2o_air_xylene_system.hh b/dumux/material/fluidsystems/h2o_air_xylene_system.hh
new file mode 100644
index 0000000000000000000000000000000000000000..69969c535424df66ed9935d46852e10d6c07e544
--- /dev/null
+++ b/dumux/material/fluidsystems/h2o_air_xylene_system.hh
@@ -0,0 +1,451 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * Copyright (C) 2011 by Holger Class *
+ * 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 A fluid system with water, gas and NAPL as phases and
+ * \f$H_2O\f$ and \f$Air\f$ and \f$NAPL (contaminant)\f$ as components.
+ */
+#ifndef DUMUX_H2O_AIR_XYLENE_SYSTEM_HH
+#define DUMUX_H2O_AIR_XYLENE_SYSTEM_HH
+
+#include <dumux/material/idealgas.hh>
+#include <dumux/material/components/air.hh>
+#include <dumux/material/components/h2o.hh>
+#include <dumux/material/components/xylene.hh>
+#include <dumux/material/components/tabulatedcomponent.hh>
+#include <dumux/material/binarycoefficients/h2o_air.hh>
+#include <dumux/material/binarycoefficients/h2o_xylene.hh>
+#include <dumux/material/binarycoefficients/air_xylene.hh>
+#include <dumux/common/propertysystem.hh>
+
+namespace Dumux
+{
+
+// forward declarations of property tags
+namespace Properties
+{
+NEW_PROP_TAG(Scalar);
+}
+
+/*!
+ * \brief A compositional fluid with water and molecular nitrogen as
+ * components in both, the liquid and the gas phase.
+ */
+template <class TypeTag>
+class H2O_Air_Xylene_System
+{
+ typedef H2O_Air_Xylene_System<TypeTag> ThisType;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+
+ typedef Dumux::IdealGas<Scalar> IdealGas;
+
+public:
+ typedef Dumux::H2O<Scalar> H2O;
+ typedef Dumux::Xylene<Scalar> NAPL;
+ typedef Dumux::Air<Scalar> Air;
+
+ static const int numComponents = 3;
+ static const int numPhases = 3;
+
+ static const int wPhaseIdx = 0; // index of the water phase
+ static const int nPhaseIdx = 1; // index of the NAPL phase
+ static const int gPhaseIdx = 2; // index of the gas phase
+
+ static const int H2OIdx = 0;
+ static const int NAPLIdx = 1;
+ static const int airIdx = 2;
+
+ // TODO is it possible to import the phase state definitions directly from 3p3cindices.hh ??
+ static const int threePhases = 1;
+ static const int wPhaseOnly = 2;
+ static const int gnPhaseOnly = 3;
+ static const int wnPhaseOnly = 4;
+ static const int gPhaseOnly = 5;
+ static const int wgPhaseOnly = 6;
+
+ static void init()
+ { }
+
+ /*!
+ * \brief Return the human readable name of a phase (used in indices)
+ */
+ static const char *phaseName(int phaseIdx)
+ {
+ switch (phaseIdx) {
+ case wPhaseIdx: return "w";
+ case nPhaseIdx: return "n";
+ case gPhaseIdx: return "g";;
+ };
+ DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
+ }
+
+ /*!
+ * \brief Return the human readable name of a component (used in indices)
+ */
+ static const char *componentName(int compIdx)
+ {
+ switch (compIdx) {
+ case H2OIdx: return H2O::name();
+ case airIdx: return Air::name();
+ case NAPLIdx: return NAPL::name();
+ };
+ DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
+ }
+
+ /*!
+ * \brief Return the molar mass of a component in [kg/mol].
+ */
+ static Scalar molarMass(int compIdx)
+ {
+ switch (compIdx) {
+ case H2OIdx: return H2O::molarMass();
+ case airIdx: return Air::molarMass();
+ case NAPLIdx: return NAPL::molarMass();
+ };
+ DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
+ }
+
+
+ /*!
+ * \brief Given all mole fractions in a phase, return the phase
+ * density [kg/m^3].
+ */
+ template <class FluidState>
+ static Scalar density(int phaseIdx,
+ const FluidState &fluidState)
+ {
+ if (phaseIdx == wPhaseIdx) {
+ // See: Ochs 2008
+ // \todo: proper citation
+ Scalar rholH2O = H2O::liquidDensity(fluidState.temperature(), fluidState.pressure(phaseIdx));
+ Scalar clH2O = rholH2O/H2O::molarMass();
+
+ // this assumes each dissolved molecule displaces exactly one
+ // water molecule in the liquid
+ return
+ clH2O*(H2O::molarMass()*fluidState.moleFraction(wPhaseIdx, H2OIdx)
+ +
+ Air::molarMass()*fluidState.moleFraction(wPhaseIdx, airIdx)
+ +
+ NAPL::molarMass()*fluidState.moleFraction(wPhaseIdx, NAPLIdx));
+ }
+ else if (phaseIdx == nPhaseIdx) {
+ // assume pure NAPL for the NAPL phase
+ return NAPL::liquidDensity(fluidState.temperature(), fluidState.pressure(phaseIdx));
+ }
+ else if (phaseIdx == gPhaseIdx) {
+ Scalar fugH2O =
+ fluidState.moleFraction(gPhaseIdx, H2OIdx) *
+ fluidState.pressure(gPhaseIdx);
+ Scalar fugAir =
+ fluidState.moleFraction(gPhaseIdx, airIdx) *
+ fluidState.pressure(gPhaseIdx);
+ Scalar fugNAPL =
+ fluidState.moleFraction(gPhaseIdx, NAPLIdx) *
+ fluidState.pressure(gPhaseIdx);
+ return
+ H2O::gasDensity(fluidState.temperature(), fugH2O) +
+ Air::gasDensity(fluidState.temperature(), fugAir) +
+ NAPL::gasDensity(fluidState.temperature(), fugNAPL);
+
+ }
+ DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
+ }
+
+ /*!
+ * \brief Return the viscosity of a phase.
+ */
+ template <class FluidState>
+ static Scalar viscosity(int phaseIdx,
+ const FluidState &fluidState)
+ {
+ if (phaseIdx == wPhaseIdx) {
+ // assume pure water viscosity
+ return H2O::liquidViscosity(fluidState.temperature(),
+ fluidState.pressure(phaseIdx));
+ }
+ else if (phaseIdx == nPhaseIdx) {
+ // assume pure NAPL viscosity
+ return NAPL::liquidViscosity(fluidState.temperature(), fluidState.pressure(phaseIdx));
+ }
+ else if (phaseIdx == gPhaseIdx) {
+ /* Wilke method. See:
+ *
+ * See: R. Reid, et al.: The Properties of Gases and Liquids,
+ * 4th edition, McGraw-Hill, 1987, 407-410
+ * 5th edition, McGraw-Hill, 20001, p. 9.21/22
+ *
+ * in this case, we use a simplified version in order to avoid
+ * computationally costly evaluation of sqrt and pow functions and
+ * divisions
+ * -- compare e.g. with Promo Class p. 32/33
+ */
+ Scalar muResult;
+ const Scalar mu[numComponents] = {
+ H2O::gasViscosity(fluidState.temperature(), H2O::vaporPressure(fluidState.temperature())),
+ Air::simpleGasViscosity(fluidState.temperature(), fluidState.pressure(phaseIdx)),
+ NAPL::gasViscosity(fluidState.temperature(), NAPL::vaporPressure(fluidState.temperature()))
+ };
+ // molar masses
+ const Scalar M[numComponents] = {
+ H2O::molarMass(),
+ Air::molarMass(),
+ NAPL::molarMass()
+ };
+
+ Scalar muAW = mu[airIdx]*fluidState.moleFraction(gPhaseIdx, airIdx)
+ + mu[H2OIdx]*fluidState.moleFraction(gPhaseIdx, H2OIdx)
+ / (fluidState.moleFraction(gPhaseIdx, airIdx)
+ + fluidState.moleFraction(gPhaseIdx, H2OIdx));
+ Scalar xAW = fluidState.moleFraction(gPhaseIdx, airIdx)
+ + fluidState.moleFraction(gPhaseIdx, H2OIdx);
+
+ Scalar MAW = (fluidState.moleFraction(gPhaseIdx, airIdx)*Air::molarMass()
+ + fluidState.moleFraction(gPhaseIdx, H2OIdx)*H2O::molarMass());
+
+ /* TODO, please check phiCAW for the Xylene case here */
+
+ Scalar phiCAW = 0.3; // simplification for this particular system
+ /* actually like this
+ * Scalar phiCAW = std::pow(1.+std::sqrt(mu[NAPLIdx]/muAW)*std::pow(MAW/M[NAPLIdx],0.25),2)
+ * / std::sqrt(8.*(1.+M[NAPLIdx]/MAW));
+ */
+ Scalar phiAWC = phiCAW * muAW*M[NAPLIdx]/(mu[NAPLIdx]*MAW);
+
+ muResult = (xAW*muAW)/(xAW+fluidState.moleFraction(gPhaseIdx, NAPLIdx)*phiAWC)
+ + (fluidState.moleFraction(gPhaseIdx, NAPLIdx) * mu[NAPLIdx])
+ / (fluidState.moleFraction(gPhaseIdx, NAPLIdx) + xAW*phiCAW);
+ return muResult;
+ }
+ DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
+ }
+
+
+ /*!
+ * \brief Given all mole fractions, return the diffusion
+ * coefficent of a component in a phase.
+ */
+ template <class FluidState>
+ static Scalar diffusionCoefficient(int phaseIdx,
+ int compIdx,
+ const FluidState &fluidState)
+ {
+ Scalar diffCont;
+
+ if (phaseIdx==gPhaseIdx) {
+ Scalar diffAC = Dumux::BinaryCoeff::Air_Xylene::gasDiffCoeff(fluidState.temperature(), fluidState.pressure(phaseIdx));
+ Scalar diffWC = Dumux::BinaryCoeff::H2O_Xylene::gasDiffCoeff(fluidState.temperature(), fluidState.pressure(phaseIdx));
+ Scalar diffAW = Dumux::BinaryCoeff::H2O_Air::gasDiffCoeff(fluidState.temperature(), fluidState.pressure(phaseIdx));
+
+ const Scalar xga = fluidState.moleFraction(gPhaseIdx, airIdx);
+ const Scalar xgw = fluidState.moleFraction(gPhaseIdx, H2OIdx);
+ const Scalar xgc = fluidState.moleFraction(gPhaseIdx, NAPLIdx);
+
+ if (compIdx==NAPLIdx) return (1.- xgw)/(xga/diffAW + xgc/diffWC);
+ else if (compIdx==H2OIdx) return (1.- xgc)/(xgw/diffWC + xga/diffAC);
+ else if (compIdx==airIdx) DUNE_THROW(Dune::InvalidStateException,
+ "Diffusivity of air in the gas phase "
+ "is constraint by sum of diffusive fluxes = 0 !\n");
+ } else if (phaseIdx==wPhaseIdx){
+ Scalar diffACl = 1.e-9; // BinaryCoeff::Air_Xylene::liquidDiffCoeff(temperature, pressure);
+ Scalar diffWCl = 1.e-9; // BinaryCoeff::H2O_Xylene::liquidDiffCoeff(temperature, pressure);
+ Scalar diffAWl = 1.e-9; // BinaryCoeff::H2O_Air::liquidDiffCoeff(temperature, pressure);
+
+ Scalar xwa = fluidState.moleFraction(wPhaseIdx, airIdx);
+ Scalar xww = fluidState.moleFraction(wPhaseIdx, H2OIdx);
+ Scalar xwc = fluidState.moleFraction(wPhaseIdx, NAPLIdx);
+
+ switch (compIdx) {
+ case NAPLIdx:
+ diffCont = (1.- xww)/(xwa/diffAWl + xwc/diffWCl);
+ return diffCont;
+ case airIdx:
+ diffCont = (1.- xwc)/(xww/diffWCl + xwa/diffACl);
+ return diffCont;
+ case H2OIdx:
+ DUNE_THROW(Dune::InvalidStateException,
+ "Diffusivity of water in the water phase "
+ "is constraint by sum of diffusive fluxes = 0 !\n");
+ }
+ } else if (phaseIdx==nPhaseIdx) {
+
+ DUNE_THROW(Dune::InvalidStateException,
+ "Diffusion coefficients of "
+ "substances in liquid phase are undefined!\n");
+ }
+ return 0;
+ }
+
+ /*!
+ * \brief Returns the vapor pressures of components in the gas phase which
+ * are normally present as liquid: water and NAPL
+ */
+ template <class FluidState>
+ static Scalar vaporPressure(const FluidState &fluidState, int phaseIdx, int compIdx)
+ {
+ if (phaseIdx == gPhaseIdx) {
+ Scalar T = fluidState.temperature(phaseIdx);
+ Scalar p = fluidState.pressure(phaseIdx);
+
+ if (compIdx == H2OIdx)
+ return H2O::vaporPressure(T);
+ if (compIdx == NAPLIdx)
+ return NAPL::vaporPressure(T);
+
+ // should not be required, since it is no physically transparent implementation
+ if (compIdx == airIdx)
+ return p - H2O::vaporPressure(T) - NAPL::vaporPressure(T);
+ }
+ else
+ DUNE_THROW(Dune::NotImplemented, "vapor pressure method implemented only for gas");
+
+ return 0;
+ }
+
+
+ /*!
+ * \brief Returns the fugacity coefficient [-] of a component in a
+ * phase.
+ *
+ * In this case, things are actually pretty simple. We have an ideal
+ * solution. Thus, the fugacity coefficient is 1 in the gas phase
+ * (fugacity equals the partial pressure of the component in the gas phase
+ * respectively in the liquid phases it is the inverse of the
+ * Henry coefficients scaled by pressure
+ */
+ template <class FluidState>
+ static Scalar fugacityCoefficient(const FluidState &fluidState,
+ int phaseIdx,
+ int compIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+ assert(0 <= compIdx && compIdx < numComponents);
+
+ Scalar T = fluidState.temperature(phaseIdx);
+ Scalar p = fluidState.pressure(phaseIdx);
+
+ if (phaseIdx == wPhaseIdx) {
+ if (compIdx == H2OIdx)
+ return H2O::vaporPressure(T)/p;
+ if (compIdx == airIdx)
+ return Dumux::BinaryCoeff::H2O_Air::henry(T)/p;
+ if (compIdx == NAPLIdx)
+ return Dumux::BinaryCoeff::H2O_Xylene::henry(T)/p;
+ }
+ // for the NAPL phase, we assume currently that nothing is dissolved
+ else if (phaseIdx == nPhaseIdx) {
+ if (compIdx == NAPLIdx)
+ return 1;
+ if (compIdx == airIdx)
+ return 0;
+ if (compIdx == H2OIdx)
+ return 0;
+ }
+ // for the gas phase, assume an ideal gas when it comes to
+ // fugacity (-> fugacity == partial pressure)
+ else if (phaseIdx == gPhaseIdx)
+ {
+ return 1.0;
+ }
+ else
+ DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
+
+ return 0;
+ }
+
+
+ /*!
+ * \brief Given all mole fractions in a phase, return the specific
+ * phase enthalpy [J/kg].
+ */
+ /*!
+ * \todo This system neglects the contribution of gas-molecules in the liquid phase.
+ * This contribution is probably not big. Somebody would have to find out the enthalpy of solution for this system. ...
+ */
+ template <class FluidState>
+ static Scalar enthalpy(int phaseIdx,
+ const FluidState &fluidState)
+ {
+ if (phaseIdx == wPhaseIdx) {
+ return H2O::liquidEnthalpy(fluidState.temperature(), fluidState.pressure(phaseIdx));
+ }
+ else if (phaseIdx == nPhaseIdx) {
+ return NAPL::liquidEnthalpy(fluidState.temperature(), fluidState.pressure(phaseIdx));
+ }
+ else if (phaseIdx == gPhaseIdx) { // gas phase enthalpy depends strongly on composition
+ Scalar hgc = NAPL::gasEnthalpy(fluidState.temperature(),
+ fluidState.pressure(phaseIdx));
+ Scalar hgw = H2O::gasEnthalpy(fluidState.temperature(), fluidState.pressure(phaseIdx));
+ Scalar hga = Air::gasEnthalpy(fluidState.temperature(), fluidState.pressure(phaseIdx)); // pressure is only a dummy here
+
+ Scalar result = 0;
+ result += hgw * fluidState.massFraction(gPhaseIdx, H2OIdx);
+ result += hga * fluidState.massFraction(gPhaseIdx, airIdx);
+ result += hgc * fluidState.massFraction(gPhaseIdx, NAPLIdx);
+
+ return result;
+ }
+ DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
+ }
+
+ /*!
+ * \brief Given all mole fractions in a phase, return the phase's
+ * internal energy [J/kg].
+ */
+ template <class FluidState>
+ static Scalar internalEnergy(int phaseIdx,
+ const FluidState &fluidState)
+ {
+ return
+ enthalpy(phaseIdx, fluidState) -
+ fluidState.pressure(phaseIdx)/density(phaseIdx, fluidState);
+ }
+
+private:
+ static Scalar waterPhaseDensity_(Scalar T, Scalar pw, Scalar xww, Scalar xwa, Scalar xwc)
+ {
+ Scalar rholH2O = H2O::liquidDensity(T, pw);
+ Scalar clH2O = rholH2O/H2O::molarMass();
+
+ // this assumes each dissolved molecule displaces exactly one
+ // water molecule in the liquid
+ return
+ clH2O*(xww*H2O::molarMass() + xwa*Air::molarMass() + xwc*NAPL::molarMass());
+ }
+
+ static Scalar gasPhaseDensity_(Scalar T, Scalar pg, Scalar xgw, Scalar xga, Scalar xgc)
+ {
+ return H2O::gasDensity(T, pg*xgw) + Air::gasDensity(T, pg*xga) + NAPL::gasDensity(T, pg*xgc);
+ };
+
+ static Scalar NAPLPhaseDensity_(Scalar T, Scalar pn)
+ {
+ return NAPL::liquidDensity(T, pn);
+ }
+
+};
+
+} // end namespace
+
+#endif
diff --git a/test/boxmodels/2p/test_2p.cc b/test/boxmodels/2p/test_2p.cc
index 49b19114b65a45ca1db9c36283f2788231e31203..4bd2fae28e30e7ff66a1adecb581145c28094099 100644
--- a/test/boxmodels/2p/test_2p.cc
+++ b/test/boxmodels/2p/test_2p.cc
@@ -35,6 +35,7 @@
#include <dune/common/exceptions.hh>
#include <dune/common/mpihelper.hh>
+#include <dune/grid/uggrid.hh>
#include <iostream>
diff --git a/test/boxmodels/Makefile.am b/test/boxmodels/Makefile.am
index bcebc872c57a4ace2ab7ce7452f4161db679cfd5..3c69177f01fc3558ea794c9b56988d568f075a36 100644
--- a/test/boxmodels/Makefile.am
+++ b/test/boxmodels/Makefile.am
@@ -1,4 +1,4 @@
-SUBDIRS = 1p 1p2c 2p 2p2c 2p2cni 2pni MpNc richards
+SUBDIRS = 1p 1p2c 2p 2p2c 2p2cni 2pni 3p3c 3p3cni MpNc richards
EXTRA_DIST = CMakeLists.txt