From eee9935cfddf16e4f73cb2f455047e67ef69c358 Mon Sep 17 00:00:00 2001
From: Andreas Lauser <and@poware.org>
Date: Fri, 16 Dec 2011 13:23:13 +0000
Subject: [PATCH] 1p box model: make it use fluid systems
it used to call Fluid::foo() directly. now we do this by means of the
fluid system, so that the same problems for 1p2c can be reused by
1p. Also has the advantage that all box models now provide a fluid
state which makes writing generic code much simpler...
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@7053 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/boxmodels/1p/1pproperties.hh | 6 +-
dumux/boxmodels/1p/1ppropertydefaults.hh | 23 ++
dumux/boxmodels/1p/1pvolumevariables.hh | 81 +++-
dumux/boxmodels/2p/2pvolumevariables.hh | 39 +-
dumux/boxmodels/2p2c/2p2cvolumevariables.hh | 3 +-
dumux/boxmodels/2pni/2pnivolumevariables.hh | 126 +++---
.../MpNcfluidsystems/1pfluidsystem.hh | 371 ++++++++++++++++++
.../2pimmisciblefluidsystem.hh | 42 +-
.../MpNcfluidsystems/h2on2fluidsystem.hh | 53 +--
test/boxmodels/1p/1ptestproblem.hh | 3 +-
test/boxmodels/1p/grids/test_1p_2d.dgf | 2 +-
11 files changed, 588 insertions(+), 161 deletions(-)
create mode 100644 dumux/material/MpNcfluidsystems/1pfluidsystem.hh
diff --git a/dumux/boxmodels/1p/1pproperties.hh b/dumux/boxmodels/1p/1pproperties.hh
index 6040cf40f5..ada38886cf 100644
--- a/dumux/boxmodels/1p/1pproperties.hh
+++ b/dumux/boxmodels/1p/1pproperties.hh
@@ -53,13 +53,11 @@ NEW_TYPE_TAG(BoxOneP, INHERITS_FROM(BoxModel));
NEW_PROP_TAG(NumPhases); //!< Number of fluid phases in the system
NEW_PROP_TAG(OnePIndices); //!< Enumerations for the 1p models
NEW_PROP_TAG(SpatialParameters); //!< The type of the spatial parameters object
-NEW_PROP_TAG(Fluid); //!< The fluid for the single-phase problems
+NEW_PROP_TAG(FluidSystem); //!< The type of the fluid system to use
+NEW_PROP_TAG(Fluid); //!< The fluid used for the default fluid system
NEW_PROP_TAG(EnableGravity); //!< Returns whether gravity is considered in the problem
NEW_PROP_TAG(UpwindWeight); //!< Returns weight of the upwind cell when calculating fluxes
-// use the CG solver preconditioned by ILU-0
-SET_TYPE_PROP(BoxOneP, LinearSolver, Dumux::BoxCGILU0Solver<TypeTag> );
-
// \}
};
diff --git a/dumux/boxmodels/1p/1ppropertydefaults.hh b/dumux/boxmodels/1p/1ppropertydefaults.hh
index 106cc14691..2781aa9d9c 100644
--- a/dumux/boxmodels/1p/1ppropertydefaults.hh
+++ b/dumux/boxmodels/1p/1ppropertydefaults.hh
@@ -37,6 +37,12 @@
#include "1pfluxvariables.hh"
#include "1pindices.hh"
+#include <dumux/material/fluidsystems/gasphase.hh>
+#include <dumux/material/fluidsystems/liquidphase.hh>
+#include <dumux/material/components/nullcomponent.hh>
+
+#include <dumux/material/MpNcfluidsystems/1pfluidsystem.hh>
+
namespace Dumux
{
// \{
@@ -69,6 +75,23 @@ SET_TYPE_PROP(BoxOneP, OnePIndices, OnePIndices);
//! fluxes. Use central differences by default.
SET_SCALAR_PROP(BoxOneP, UpwindWeight, 0.5);
+//! The fluid system to use by default
+SET_PROP(BoxOneP, FluidSystem)
+{ private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Fluid)) Fluid;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+
+public:
+ typedef OnePFluidSystem<Scalar, Fluid> type;
+};
+
+SET_PROP(BoxOneP, Fluid)
+{ private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+public:
+ typedef Dumux::LiquidPhase<Scalar, Dumux::NullComponent<Scalar> > type;
+};
+
// \}
};
diff --git a/dumux/boxmodels/1p/1pvolumevariables.hh b/dumux/boxmodels/1p/1pvolumevariables.hh
index 8cb85d474b..356d29d5dc 100644
--- a/dumux/boxmodels/1p/1pvolumevariables.hh
+++ b/dumux/boxmodels/1p/1pvolumevariables.hh
@@ -29,6 +29,8 @@
#include "1pproperties.hh"
#include <dumux/boxmodels/common/boxvolumevariables.hh>
+#include <dumux/material/MpNcfluidstates/immisciblefluidstate.hh>
+
namespace Dumux
{
@@ -42,13 +44,13 @@ template <class TypeTag>
class OnePVolumeVariables : public BoxVolumeVariables<TypeTag>
{
typedef BoxVolumeVariables<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) Implementation;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) Implementation;
typedef typename GridView::template Codim<0>::Entity Element;
@@ -59,12 +61,15 @@ class OnePVolumeVariables : public BoxVolumeVariables<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(OnePIndices)) Indices;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Fluid)) Fluid;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef Dune::FieldVector<Scalar, dim> LocalPosition;
public:
+ //! Type of the fluid state
+ typedef Dumux::ImmiscibleFluidState<Scalar, FluidSystem> FluidState;
+
/*!
* \brief Update all quantities for a given control volume.
*
@@ -84,20 +89,35 @@ public:
{
ParentType::update(priVars, problem, element, elemGeom, scvIdx, isOldSol);
- this->asImp_().updateTemperature_(priVars,
- element,
- elemGeom,
- scvIdx,
- problem);
+ // set the phase temperature and pressure
+ asImp_().updateTemperature_(priVars,
+ element,
+ elemGeom,
+ scvIdx,
+ problem);
+ fluidState_.setPressure(/*phaseIdx=*/0, priVars[Indices::pressureIdx]);
- pressure_ = priVars[Indices::pressureIdx];
- density_ = Fluid::density(temperature(), pressure());
- viscosity_ = Fluid::viscosity(temperature(), pressure());
+ // saturation in a single phase is always 1 and thus redundant
+ // to set. But since we use the fluid state shared by the
+ // immiscible multi-phase models, so we have to set it here...
+ fluidState_.setSaturation(/*phaseIdx=*/0, 1.0);
+
+ typename FluidSystem::ParameterCache paramCache;
+ paramCache.updatePhase(fluidState_, /*phaseIdx=*/0);
+
+ Scalar value = FluidSystem::density(fluidState_, paramCache, /*phaseIdx=*/0);
+ fluidState_.setDensity(/*phaseIdx=*/0, value);
+
+ value = FluidSystem::viscosity(fluidState_, paramCache, /*phaseIdx=*/0);
+ fluidState_.setViscosity(/*phaseIdx=*/0, value);
// porosity
porosity_ = problem.spatialParameters().porosity(element,
elemGeom,
scvIdx);
+
+ // energy related quantities
+ asImp_().updateEnergy_(paramCache, priVars, problem, element, elemGeom, scvIdx, isOldSol);
};
/*!
@@ -108,28 +128,28 @@ public:
* identical.
*/
Scalar temperature() const
- { return temperature_; }
+ { return fluidState_.temperature(); }
/*!
* \brief Returns the effective pressure of a given phase within
* the control volume.
*/
Scalar pressure() const
- { return pressure_; }
+ { return fluidState_.pressure(/*phaseIdx=*/0); }
/*!
* \brief Returns the mass density of a given phase within the
* control volume.
*/
Scalar density() const
- { return density_; }
+ { return fluidState_.density(/*phaseIdx=*/0); }
/*!
* \brief Returns the dynamic viscosity of the fluid within the
* control volume.
*/
Scalar viscosity() const
- { return viscosity_; }
+ { return fluidState_.viscosity(/*phaseIdx=*/0); }
/*!
* \brief Returns the average porosity within the control volume.
@@ -137,19 +157,42 @@ public:
Scalar porosity() const
{ return porosity_; }
+ /*!
+ * \brief Returns the fluid state of the control volume.
+ */
+ const FluidState &fluidState() const
+ { return fluidState_; }
+
protected:
void updateTemperature_(const PrimaryVariables &priVars,
const Element &element,
const FVElementGeometry &elemGeom,
int scvIdx,
const Problem &problem)
- { temperature_ = problem.boxTemperature(element, elemGeom, scvIdx); }
+ { fluidState_.setTemperature(problem.boxTemperature(element, elemGeom, scvIdx)); }
- Scalar temperature_;
- Scalar pressure_;
- Scalar density_;
- Scalar viscosity_;
+ /*!
+ * \brief Called by update() to compute the energy related quantities
+ */
+ template <class ParameterCache>
+ void updateEnergy_(ParameterCache ¶mCache,
+ const PrimaryVariables &sol,
+ const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int vertIdx,
+ bool isOldSol)
+ { }
+
+ FluidState fluidState_;
Scalar porosity_;
+
+private:
+ Implementation &asImp_()
+ { return *static_cast<Implementation*>(this); }
+
+ const Implementation &asImp_() const
+ { return *static_cast<const Implementation*>(this); }
};
}
diff --git a/dumux/boxmodels/2p/2pvolumevariables.hh b/dumux/boxmodels/2p/2pvolumevariables.hh
index 97f9216239..50ec393383 100644
--- a/dumux/boxmodels/2p/2pvolumevariables.hh
+++ b/dumux/boxmodels/2p/2pvolumevariables.hh
@@ -102,12 +102,12 @@ public:
scvIdx,
isOldSol);
- asImp().updateTemperature_(priVars,
- element,
- elemGeom,
- scvIdx,
- problem);
-
+ asImp_().updateTemperature_(priVars,
+ element,
+ elemGeom,
+ scvIdx,
+ problem);
+
// material law parameters
const MaterialLawParams &materialParams =
problem.spatialParameters().materialLawParams(element, elemGeom, scvIdx);
@@ -126,6 +126,16 @@ public:
porosity_ = problem.spatialParameters().porosity(element,
elemGeom,
scvIdx);
+
+
+#warning "TODO: it would be nice to avoid creating the parameter cache a second time here," \
+ " but that requires to abandon/modify model.completeFluidState()." \
+ " (alternatively a completeFluidState() method can be introduced in the 2pni model.)"
+ typename FluidSystem::ParameterCache paramCache;
+ paramCache.updateAll(fluidState_);
+
+ // energy related quantities
+ asImp_().updateEnergy_(paramCache, priVars, problem, element, elemGeom, scvIdx, isOldSol);
}
/*!
@@ -202,15 +212,28 @@ protected:
fluidState_.setTemperature(problem.boxTemperature(element, elemGeom, scvIdx));
}
+ /*!
+ * \brief Called by update() to compute the energy related quantities
+ */
+ template <class ParameterCache>
+ void updateEnergy_(ParameterCache ¶mCache,
+ const PrimaryVariables &sol,
+ const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int vertIdx,
+ bool isOldSol)
+ { }
+
FluidState fluidState_;
Scalar porosity_;
Scalar mobility_[numPhases];
private:
- Implementation &asImp()
+ Implementation &asImp_()
{ return *static_cast<Implementation*>(this); }
- const Implementation &asImp() const
+ const Implementation &asImp_() const
{ return *static_cast<const Implementation*>(this); }
};
diff --git a/dumux/boxmodels/2p2c/2p2cvolumevariables.hh b/dumux/boxmodels/2p2c/2p2cvolumevariables.hh
index 01beda45f6..6ee0c2461c 100644
--- a/dumux/boxmodels/2p2c/2p2cvolumevariables.hh
+++ b/dumux/boxmodels/2p2c/2p2cvolumevariables.hh
@@ -266,8 +266,7 @@ protected:
const FVElementGeometry &elemGeom,
int vertIdx,
bool isOldSol)
- {
- }
+ { }
Scalar porosity_; //!< Effective porosity within the control volume
Scalar relativePermeability_[numPhases]; //!< Relative permeability within the control volume
diff --git a/dumux/boxmodels/2pni/2pnivolumevariables.hh b/dumux/boxmodels/2pni/2pnivolumevariables.hh
index f5d636b67f..eb04350fab 100644
--- a/dumux/boxmodels/2pni/2pnivolumevariables.hh
+++ b/dumux/boxmodels/2pni/2pnivolumevariables.hh
@@ -64,76 +64,6 @@ class TwoPNIVolumeVariables : public TwoPVolumeVariables<TypeTag>
typedef Dune::FieldVector<Scalar, numPhases> PhasesVector;
public:
- /*!
- * \brief Update all quantities for a given control volume.
- *
- * \param priVars The local primary variable vector
- * \param problem The problem object
- * \param element The current element
- * \param elemGeom The finite-volume geometry in the box scheme
- * \param scvIdx 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 &priVars,
- const Problem &problem,
- const Element &element,
- const FVElementGeometry &elemGeom,
- int scvIdx,
- bool isOldSol)
- {
- typedef Indices I;
-
- // vertex update data for the mass balance
- ParentType::update(priVars,
- problem,
- element,
- elemGeom,
- scvIdx,
- isOldSol);
-
- typename FluidSystem::ParameterCache paramCache;
- // TODO: this calculates the cached parameters a second time
- // and is thus inefficient!
- paramCache.updateAll(this->fluidState());
-
- // the internal energies
- for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
- internalEnergy_[phaseIdx] =
- FluidSystem::internalEnergy(this->fluidState(),
- paramCache,
- phaseIdx);
- }
- Valgrind::CheckDefined(internalEnergy_);
- }
-
- /*!
- * \brief Update the temperature for a given control volume.
- *
- * \param priVars The local primary variable vector
- * \param element The current element
- * \param elemGeom The finite-volume geometry in the box scheme
- * \param scvIdx The local index of the SCV (sub-control volume)
- * \param problem The problem object
- *
- */
-
- // this method gets called by the parent class
- void updateTemperature_(const PrimaryVariables &priVars,
- const Element &element,
- const FVElementGeometry &elemGeom,
- int scvIdx,
- const Problem &problem)
- {
- // retrieve temperature from primary variables
- this->fluidState_.setTemperature(priVars[temperatureIdx]);
-
- heatCapacity_ =
- problem.spatialParameters().heatCapacity(element, elemGeom, scvIdx);
-
- Valgrind::CheckDefined(heatCapacity_);
- }
-
/*!
* \brief Returns the total internal energy of a phase in the
* sub-control volume.
@@ -142,7 +72,7 @@ public:
*
*/
Scalar internalEnergy(int phaseIdx) const
- { return internalEnergy_[phaseIdx]; };
+ { return this->fluidState_.internalEnergy(phaseIdx); };
/*!
* \brief Returns the total enthalpy of a phase in the sub-control
@@ -151,10 +81,7 @@ public:
* \param phaseIdx The phase index
*/
Scalar enthalpy(int phaseIdx) const
- {
- return internalEnergy_[phaseIdx]
- + this->fluidState().pressure(phaseIdx)
- / this->fluidState().density(phaseIdx); };
+ { return this->fluidState_.enthalpy(phaseIdx); };
/*!
* \brief Returns the total heat capacity \f$\mathrm{[J/K*m^3]}\f$ of the rock matrix in
@@ -164,7 +91,54 @@ public:
{ return heatCapacity_; };
protected:
- Scalar internalEnergy_[numPhases];
+ // this method gets called by the parent class. since this method
+ // is protected, we are friends with our parent..
+ friend class TwoPVolumeVariables<TypeTag>;
+
+ /*!
+ * \brief Update the temperature for a given control volume.
+ *
+ * \param priVars The local primary variable vector
+ * \param element The current element
+ * \param elemGeom The finite-volume geometry in the box scheme
+ * \param scvIdx The local index of the SCV (sub-control volume)
+ * \param problem The problem object
+ *
+ */
+ void updateTemperature_(const PrimaryVariables &priVars,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int scvIdx,
+ const Problem &problem)
+ {
+ // retrieve temperature from primary variables
+ this->fluidState_.setTemperature(priVars[temperatureIdx]);
+ }
+
+ /*!
+ * \brief Called by update() to compute the energy related quantities
+ */
+ template <class ParameterCache>
+ void updateEnergy_(ParameterCache ¶mCache,
+ const PrimaryVariables &sol,
+ const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int vertIdx,
+ bool isOldSol)
+ {
+ // copmute and set the internal energies of the fluid phases
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ Scalar u = FluidSystem::internalEnergy(this->fluidState_, paramCache, phaseIdx);
+
+ this->fluidState_.setInternalEnergy(phaseIdx, u);
+ }
+
+ // copmute and set the heat capacity of the solid phase
+ heatCapacity_ = problem.spatialParameters().heatCapacity(element, elemGeom, scvIdx);
+ Valgrind::CheckDefined(heatCapacity_);
+ }
+
Scalar heatCapacity_;
};
diff --git a/dumux/material/MpNcfluidsystems/1pfluidsystem.hh b/dumux/material/MpNcfluidsystems/1pfluidsystem.hh
new file mode 100644
index 0000000000..bac653903f
--- /dev/null
+++ b/dumux/material/MpNcfluidsystems/1pfluidsystem.hh
@@ -0,0 +1,371 @@
+/*****************************************************************************
+ * Copyright (C) 2011 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 for single phase models.
+ */
+#ifndef DUMUX_1P_FLUIDSYSTEM_HH
+#define DUMUX_1P_FLUIDSYSTEM_HH
+
+#include <dumux/material/fluidsystems/liquidphase.hh>
+#include <dumux/material/fluidsystems/gasphase.hh>
+
+#include <dune/common/exceptions.hh>
+
+#include "nullparametercache.hh"
+
+namespace Dumux {
+/*!
+ * \ingroup Fluidsystems
+ *
+ * \brief A fluid system for single phase models.
+ *
+ * The fluid is defined as a template parameter. For existing
+ * components the Dumux::LiquidPhase<Component> and
+ * Dumux::GasPhase<Component> may be used.
+ */
+template <class Scalar, class Fluid>
+class OnePFluidSystem
+{
+ typedef OnePFluidSystem<Scalar, Fluid> ThisType;
+
+public:
+ //! The type of parameter cache objects
+ typedef Dumux::NullParameterCache ParameterCache;
+
+ /****************************************
+ * Fluid phase related static parameters
+ ****************************************/
+
+ //! Number of phases in the fluid system
+ static constexpr int numPhases = 1;
+
+ /*!
+ * \brief Return the human readable name of a fluid phase
+ *
+ * \param phaseIdx The index of the fluid phase to consider
+ */
+ static const char *phaseName(int phaseIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+
+ return Fluid::name();
+ }
+
+ /*!
+ * \brief Return whether a phase is liquid
+ *
+ * \param phaseIdx The index of the fluid phase to consider
+ */
+ static bool isLiquid(int phaseIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+
+ return Fluid::isLiquid();
+ }
+
+ /*!
+ * \brief Returns true if and only if a fluid phase is assumed to
+ * be an ideal mixture.
+ *
+ * We define an ideal mixture as a fluid phase where the fugacity
+ * coefficients of all components times the pressure of the phase
+ * are indepent on the fluid composition. This assumtion is true
+ * if only a single component is involved. If you are unsure what
+ * this function should return, it is safe to return false. The
+ * only damage done will be (slightly) increased computation times
+ * in some cases.
+ *
+ * \param phaseIdx The index of the fluid phase to consider
+ */
+ static bool isIdealMixture(int phaseIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+
+ // we assume immisibility
+ return true;
+ }
+
+ /****************************************
+ * Component related static parameters
+ ****************************************/
+
+ //! Number of components in the fluid system
+ static constexpr int numComponents = 1;
+
+ /*!
+ * \brief Return the human readable name of a component
+ *
+ * \param compIdx The index of the component to consider
+ */
+ static const char *componentName(int compIdx)
+ {
+ assert(0 <= compIdx && compIdx < numComponents);
+
+ return Fluid::name();
+ }
+
+ /*!
+ * \brief Return the molar mass of a component in [kg/mol].
+ *
+ * \param compIdx index of the component
+ */
+ static Scalar molarMass(int compIdx)
+ {
+ assert(0 <= compIdx && compIdx < numComponents);
+
+ return Fluid::molarMass();
+ }
+
+ /*!
+ * \brief Critical temperature of a component [K].
+ *
+ * \param compIdx The index of the component to consider
+ */
+ static Scalar criticalTemperature(int compIdx)
+ {
+ assert(0 <= compIdx && compIdx < numComponents);
+
+ return Fluid::criticalTemperature();
+ };
+
+ /*!
+ * \brief Critical pressure of a component [Pa].
+ *
+ * \param compIdx The index of the component to consider
+ */
+ static Scalar criticalPressure(int compIdx)
+ {
+ assert(0 <= compIdx && compIdx < numComponents);
+
+ return Fluid::criticalPressure();
+ };
+
+ /*!
+ * \brief The acentric factor of a component [].
+ *
+ * \param compIdx The index of the component to consider
+ */
+ static Scalar acentricFactor(int compIdx)
+ {
+ assert(0 <= compIdx && compIdx < numComponents);
+
+ return Fluid::acentricFactor();
+ };
+
+ /****************************************
+ * thermodynamic relations
+ ****************************************/
+
+ /*!
+ * \brief Initialize the fluid system's static parameters
+ */
+ static void init()
+ { }
+
+ /*!
+ * \brief Calculate the density [kg/m^3] of a fluid phase
+ *
+ * \param fluidState An abitrary fluid state
+ * \param paramCache The fluid system's parameter cache
+ * \param phaseIdx The index of the fluid phase to consider
+ */
+ template <class FluidState>
+ static Scalar density(const FluidState &fluidState,
+ const ParameterCache ¶mCache,
+ int phaseIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+
+ Scalar temperature = fluidState.temperature(phaseIdx);
+ Scalar pressure = fluidState.pressure(phaseIdx);
+ return Fluid::density(temperature, pressure);
+ }
+
+ /*!
+ * \brief Calculate the dynamic viscosity of a fluid phase [Pa*s]
+ *
+ * \param fluidState An abitrary fluid state
+ * \param paramCache The fluid system's parameter cache
+ * \param phaseIdx The index of the fluid phase to consider
+ */
+ template <class FluidState>
+ static Scalar viscosity(const FluidState &fluidState,
+ const ParameterCache ¶mCache,
+ int phaseIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+
+ Scalar temperature = fluidState.temperature(phaseIdx);
+ Scalar pressure = fluidState.pressure(phaseIdx);
+ return Fluid::viscosity(temperature, pressure);
+ }
+
+ /*!
+ * \brief Calculate the fugacity coefficient [Pa] of an individual
+ * component in a fluid phase
+ *
+ * The fugacity coefficient \f$\phi_\kappa\f$ is connected to the
+ * fugacity \f$f_\kappa\f$ and the component's molarity
+ * \f$x_\kappa\f$ by means of the relation
+ *
+ * \f[ f_\kappa = \phi_\kappa * x_{\kappa} \f]
+ */
+ template <class FluidState>
+ static Scalar fugacityCoefficient(const FluidState &fluidState,
+ int phaseIdx,
+ int compIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+ assert(0 <= compIdx && compIdx < numComponents);
+
+ if (phaseIdx == compIdx)
+ // TODO (?): calculate the real fugacity coefficient of
+ // the component in the fluid. Probably that's not worth
+ // the effort, since the fugacity coefficient of the other
+ // component is infinite anyway...
+ return 1.0;
+ return std::numeric_limits<Scalar>::infinity();
+ }
+
+ /*!
+ * \brief Calculate the molecular diffusion coefficient for a
+ * component in a fluid phase [mol^2 * s / (kg*m^3)]
+ *
+ * Molecular diffusion of a compoent \f$\kappa\f$ is caused by a
+ * gradient of the chemical potential and follows the law
+ *
+ * \f[ J = - D \grad mu_\kappa \f]
+ *
+ * where \f$\mu_\kappa\f$ is the component's chemical potential,
+ * \f$D\f$ is the diffusion coefficient and \f$J\f$ is the
+ * diffusive flux. \f$mu_\kappa\f$ is connected to the component's
+ * fugacity \f$f_\kappa\f$ by the relation
+ *
+ * \f[ \mu_\kappa = R T_\alpha \mathrm{ln} \frac{f_\kappa}{p_\alpha} \f]
+ *
+ * where \f$p_\alpha\f$ and \f$T_\alpha\f$ are the fluid phase'
+ * pressure and temperature.
+ *
+ * \param fluidState An abitrary fluid state
+ * \param paramCache The fluid system's parameter cache
+ * \param phaseIdx The index of the fluid phase to consider
+ * \param compIdx The index of the component to consider
+ */
+ template <class FluidState>
+ static Scalar diffusionCoefficient(const FluidState &fluidState,
+ const ParameterCache ¶mCache,
+ int phaseIdx,
+ int compIdx)
+ {
+ DUNE_THROW(Dune::InvalidStateException, "Not applicable: Diffusion coefficients");
+ };
+
+ /*!
+ * \brief Given a phase's composition, temperature and pressure,
+ * return the binary diffusion coefficient for components
+ * \f$i\f$ and \f$j\f$ in this phase.
+ *
+ * \param fluidState An abitrary fluid state
+ * \param paramCache The fluid system's parameter cache
+ * \param phaseIdx The index of the fluid phase to consider
+ * \param compIIdx The index of the first component to consider
+ * \param compJIdx The index of the second component to consider
+ */
+ template <class FluidState>
+ static Scalar binaryDiffusionCoefficient(const FluidState &fluidState,
+ const ParameterCache ¶mCache,
+ int phaseIdx,
+ int compIIdx,
+ int compJIdx)
+
+ {
+ DUNE_THROW(Dune::InvalidStateException, "Not applicable: Binary diffusion coefficients");
+ };
+
+ /*!
+ * \brief Given a phase's composition, temperature, pressure and
+ * density, calculate its specific internal energy [J/kg].
+ *
+ * \param fluidState An abitrary fluid state
+ * \param paramCache The fluid system's parameter cache
+ * \param phaseIdx The index of the fluid phase to consider
+ */
+ template <class FluidState>
+ static Scalar internalEnergy(const FluidState &fluidState,
+ const ParameterCache ¶mCache,
+ int phaseIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+
+ Scalar temperature = fluidState.temperature(phaseIdx);
+ Scalar pressure = fluidState.pressure(phaseIdx);
+ return Fluid::internalEnergy(temperature, pressure);
+ }
+
+ /*!
+ * \brief Thermal conductivity of a fluid phase [W/(m^2 K/m)].
+ *
+ * Use the conductivity of air and water as a first approximation.
+ * Source:
+ * http://en.wikipedia.org/wiki/List_of_thermal_conductivities
+ *
+ * \param fluidState An abitrary fluid state
+ * \param paramCache The fluid system's parameter cache
+ * \param phaseIdx The index of the fluid phase to consider
+ */
+ template <class FluidState>
+ static Scalar thermalConductivity(const FluidState &fluidState,
+ const ParameterCache ¶mCache,
+ int phaseIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+
+ Scalar temperature = fluidState.temperature(phaseIdx);
+ Scalar pressure = fluidState.pressure(phaseIdx);
+ return Fluid::thermalConductivity(temperature, pressure);
+ }
+
+ /*!
+ * \brief Specific isobaric heat capacity of a fluid phase.
+ * \f$\mathrm{[J/kg]}\f$.
+ *
+ * \param fluidState An abitrary fluid state
+ * \param paramCache The fluid system's parameter cache
+ * \param phaseIdx The index of the fluid phase to consider
+ */
+ template <class FluidState>
+ static Scalar heatCapacity(const FluidState &fluidState,
+ const ParameterCache ¶mCache,
+ int phaseIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+
+ Scalar temperature = fluidState.temperature(phaseIdx);
+ Scalar pressure = fluidState.pressure(phaseIdx);
+ return Fluid::heatCapacity(temperature, pressure);
+ }
+
+};
+
+} // end namepace
+
+#endif
diff --git a/dumux/material/MpNcfluidsystems/2pimmisciblefluidsystem.hh b/dumux/material/MpNcfluidsystems/2pimmisciblefluidsystem.hh
index 85b5ddfc08..dc36667a27 100644
--- a/dumux/material/MpNcfluidsystems/2pimmisciblefluidsystem.hh
+++ b/dumux/material/MpNcfluidsystems/2pimmisciblefluidsystem.hh
@@ -234,6 +234,24 @@ public:
return NonWettingPhase::density(temperature, pressure);
}
+ /*!
+ * \brief Return the viscosity of a phase [Pa*s].
+ *
+ */
+ using Base::viscosity;
+ template <class FluidState>
+ static Scalar viscosity(const FluidState &fluidState,
+ int phaseIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+
+ Scalar temperature = fluidState.temperature(phaseIdx);
+ Scalar pressure = fluidState.pressure(phaseIdx);
+ if (phaseIdx == wPhaseIdx)
+ return WettingPhase::viscosity(temperature, pressure);
+ return NonWettingPhase::viscosity(temperature, pressure);
+ }
+
/*!
* \brief Calculate the fugacity coefficient [Pa] of an individual
* component in a fluid phase
@@ -262,30 +280,6 @@ public:
return std::numeric_limits<Scalar>::infinity();
}
- /*!
- * \brief Return the viscosity of a phase [Pa*s].
- *
- * \param phaseIdx index of the phase
- * \param temperature phase temperature in [K]
- * \param pressure phase pressure in [Pa]
- * \param fluidState The fluid state of the two-phase model
- * \tparam FluidState the fluid state class of the two-phase model
- * \return returns the viscosity of the phase [Pa*s]
- */
- using Base::viscosity;
- template <class FluidState>
- static Scalar viscosity(const FluidState &fluidState,
- int phaseIdx)
- {
- assert(0 <= phaseIdx && phaseIdx < numPhases);
-
- Scalar temperature = fluidState.temperature(phaseIdx);
- Scalar pressure = fluidState.pressure(phaseIdx);
- if (phaseIdx == wPhaseIdx)
- return WettingPhase::viscosity(temperature, pressure);
- return NonWettingPhase::viscosity(temperature, pressure);
- }
-
/*!
* \brief Calculate the binary molecular diffusion coefficient for
* a component in a fluid phase [mol^2 * s / (kg*m^3)]
diff --git a/dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh b/dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh
index fcece55f1e..f643546ed4 100644
--- a/dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh
+++ b/dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh
@@ -279,7 +279,7 @@ public:
}
/*!
- * \brief Calculate the molar volume [m^3/mol] of a fluid phase
+ * \brief Calculate the density [kg/m^3] of a fluid phase
*
* \param fluidState An abitrary fluid state
* \param paramCache The fluid system's parameter cache
@@ -306,6 +306,31 @@ public:
* fluidState.averageMolarMass(gPhaseIdx);
};
+ /*!
+ * \brief Calculate the dynamic viscosity of a fluid phase [Pa*s]
+ *
+ * \param fluidState An abitrary fluid state
+ * \param paramCache The fluid system's parameter cache
+ * \param phaseIdx The index of the fluid phase to consider
+ */
+ template <class FluidState>
+ static Scalar viscosity(const FluidState &fluidState,
+ const ParameterCache ¶mCache,
+ int phaseIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+
+ Scalar T = fluidState.temperature(phaseIdx);
+ Scalar p = fluidState.pressure(phaseIdx);
+ if (phaseIdx == lPhaseIdx) {
+ // assume pure water for the liquid phase
+ return H2O::liquidViscosity(T, p);
+ }
+
+ // assume pure nitrogen for the gas phase
+ return N2::gasViscosity(T, p);
+ };
+
/*!
* \brief Calculate the fugacity coefficient [Pa] of an individual
* component in a fluid phase
@@ -342,30 +367,6 @@ public:
return 1.0; // ideal gas
}
- /*!
- * \brief Calculate the dynamic viscosity of a fluid phase [Pa*s]
- *
- * \param fluidState An abitrary fluid state
- * \param paramCache The fluid system's parameter cache
- * \param phaseIdx The index of the fluid phase to consider
- */
- template <class FluidState>
- static Scalar viscosity(const FluidState &fluidState,
- const ParameterCache ¶mCache,
- int phaseIdx)
- {
- assert(0 <= phaseIdx && phaseIdx < numPhases);
-
- Scalar T = fluidState.temperature(phaseIdx);
- Scalar p = fluidState.pressure(phaseIdx);
- if (phaseIdx == lPhaseIdx) {
- // assume pure water for the liquid phase
- return H2O::liquidViscosity(T, p);
- }
-
- // assume pure nitrogen for the gas phase
- return N2::gasViscosity(T, p);
- };
/*!
* \brief Calculate the molecular diffusion coefficient for a
@@ -455,7 +456,7 @@ public:
/*!
* \brief Given a phase's composition, temperature, pressure and
- * density, calculate its specific enthalpy [J/kg].
+ * density, calculate its specific internal energy [J/kg].
*
* \todo This fluid system neglects the contribution of
* gas-molecules in the liquid phase. This contribution is
diff --git a/test/boxmodels/1p/1ptestproblem.hh b/test/boxmodels/1p/1ptestproblem.hh
index 9e67cf863c..2a5aad85cd 100644
--- a/test/boxmodels/1p/1ptestproblem.hh
+++ b/test/boxmodels/1p/1ptestproblem.hh
@@ -78,7 +78,8 @@ SET_PROP(OnePTestProblem, SpatialParameters)
// Linear solver settings
SET_TYPE_PROP(OnePTestProblem, LinearSolver, Dumux::BoxCGILU0Solver<TypeTag> );
-SET_INT_PROP(OnePTestProblem, LinearSolverVerbosity, 1);
+SET_INT_PROP(OnePTestProblem, LinearSolverVerbosity, 0);
+SET_SCALAR_PROP(OnePTestProblem, LinearSolverResidualReduction, 1e-12);
SET_INT_PROP(OnePTestProblem, PreconditionerIterations, 1);
SET_SCALAR_PROP(OnePTestProblem, PreconditionerRelaxation, 1.0);
diff --git a/test/boxmodels/1p/grids/test_1p_2d.dgf b/test/boxmodels/1p/grids/test_1p_2d.dgf
index 536c2946ca..74b5ecc050 100644
--- a/test/boxmodels/1p/grids/test_1p_2d.dgf
+++ b/test/boxmodels/1p/grids/test_1p_2d.dgf
@@ -2,7 +2,7 @@ DGF
Interval
0 0 % first corner
1 1 % second corner
-10 10 % cells in x and y direction
+100 100 % cells in x and y direction
#
BOUNDARYDOMAIN
--
GitLab