From 7c890d1cbc64f5c8605e26ceb857e05a7874d3d0 Mon Sep 17 00:00:00 2001
From: Johannes Hommel <johannes.hommel@iws.uni-stuttgart.de>
Date: Mon, 19 Nov 2018 18:23:38 +0100
Subject: [PATCH] Feature/user specified solid params in spatial params
---
.../nonisothermal/volumevariables.hh | 229 +++++++++++++++++-
.../3p3c/implicit/columnxylol/params.input | 10 +-
.../3p3c/implicit/columnxylol/problem.hh | 13 +-
.../implicit/columnxylol/spatialparams.hh | 62 ++---
4 files changed, 262 insertions(+), 52 deletions(-)
diff --git a/dumux/porousmediumflow/nonisothermal/volumevariables.hh b/dumux/porousmediumflow/nonisothermal/volumevariables.hh
index 3c7d08247f..107dec550e 100644
--- a/dumux/porousmediumflow/nonisothermal/volumevariables.hh
+++ b/dumux/porousmediumflow/nonisothermal/volumevariables.hh
@@ -25,10 +25,55 @@
#ifndef DUMUX_ENERGY_VOLUME_VARIABLES_HH
#define DUMUX_ENERGY_VOLUME_VARIABLES_HH
+#include <type_traits>
+
+#include <dumux/material/solidsystems/inertsolidphase.hh>
#include <dumux/porousmediumflow/volumevariables.hh>
namespace Dumux {
+#ifndef DOXYGEN
+namespace Detail {
+// helper struct detecting if the user-defined spatial params class has user-specified functions
+// for solidHeatCapacity, solidDensity, and solidThermalConductivity.
+// for g++ > 5.3, this can be replaced by a lambda
+template<class E, class SCV, class ES, class SS>
+struct hasSolidHeatCapacity
+{
+ template<class SpatialParams>
+ auto operator()(const SpatialParams& a)
+ -> decltype(a.solidHeatCapacity(std::declval<E>(), std::declval<SCV>(), std::declval<ES>(), std::declval<SS>()))
+ {}
+};
+
+template<class E, class SCV, class ES, class SS>
+struct hasSolidDensity
+{
+ template<class SpatialParams>
+ auto operator()(const SpatialParams& a)
+ -> decltype(a.solidDensity(std::declval<E>(), std::declval<SCV>(), std::declval<ES>(), std::declval<SS>()))
+ {}
+};
+
+template<class E, class SCV, class ES, class SS>
+struct hasSolidThermalConductivity
+{
+ template<class SpatialParams>
+ auto operator()(const SpatialParams& a)
+ -> decltype(a.solidThermalConductivity(std::declval<E>(), std::declval<SCV>(), std::declval<ES>(), std::declval<SS>()))
+ {}
+};
+
+template<class SolidSystem>
+struct isInertSolidPhase : public std::false_type {};
+
+template<class Scalar, class Component>
+struct isInertSolidPhase<SolidSystems::InertSolidPhase<Scalar, Component>> : public std::true_type {};
+
+} // end namespace Detail
+#endif
+
+
// forward declaration
template <class IsothermalTraits, class Impl, bool enableEnergyBalance>
class EnergyVolumeVariablesImplementation;
@@ -55,6 +100,7 @@ public:
using FluidState = typename IsothermalTraits::FluidState;
using SolidState = typename IsothermalTraits::SolidState;
using FluidSystem = typename IsothermalTraits::FluidSystem;
+ using SpatialParams = typename IsothermalTraits::SpatialParams;
//! The temperature is obtained from the problem as a constant for isothermal models
template<class ElemSol, class Problem, class Element, class Scv>
@@ -163,15 +209,14 @@ public:
const Scv &scv,
SolidState & solidState)
{
+ Scalar cs = solidHeatCapacity_(elemSol, problem, element, scv, solidState);
+ solidState.setHeatCapacity(cs);
- Scalar solidHeatCapacity = SolidSystem::heatCapacity(solidState);
- solidState.setHeatCapacity(solidHeatCapacity);
-
- Scalar solidDensity = SolidSystem::density(solidState);
- solidState.setDensity(solidDensity);
+ Scalar rhos = solidDensity_(elemSol, problem, element, scv, solidState);
+ solidState.setDensity(rhos);
- Scalar solidThermalConductivity = SolidSystem::thermalConductivity(solidState);
- solidState.setThermalConductivity(solidThermalConductivity);
+ Scalar lambdas = solidThermalConductivity_(elemSol, problem, element, scv, solidState);
+ solidState.setThermalConductivity(lambdas);
}
/*!
@@ -250,6 +295,176 @@ protected:
const Impl &asImp_() const { return *static_cast<const Impl*>(this); }
Impl &asImp_() { return *static_cast<Impl*>(this); }
+private:
+ /*!
+ * It has to be decided if the full solid system / solid state interface is used (general option, but more complicated),
+ * or the simple nonisothermal spatial params interface (simpler but less general).
+ * In the simple nonisothermal spatial params interface the functions solidHeatCapacity, solidDensity, and solidThermalConductivity
+ * in the spatial params overwrite the parameters given in the solid system. This only makes sense in combination
+ * with the simplest solid system InertSolidPhase, and can be used to quickly change parameters in certain domain regions.
+ * For setups with more general solids with several components these functions should not exist. Instead, the solid system
+ * determines the values for solidHeatCapacity, solidDensity, and solidThermalConductivity depending on the given composition.
+ */
+
+ /*!
+ * \name Access functions for the solidsystem / solidstate interface
+ */
+ // \{
+
+ /*!
+ * \brief get the solid heat capacity in an scv
+ * \param elemSol the element solution vector
+ * \param problem the problem to solve
+ * \param element the element (codim-0-entity) the scv belongs to
+ * \param scv the sub control volume
+ * \param solidState the solid state
+ * \note this gets selected if the user uses the solidsystem / solidstate interface
+ */
+ template<class ElemSol, class Problem, class Element, class Scv>
+ auto solidHeatCapacity_(const ElemSol& elemSol,
+ const Problem& problem,
+ const Element& element,
+ const Scv& scv,
+ const SolidState& solidState)
+ -> typename std::enable_if_t<!decltype(
+ isValid(Detail::hasSolidHeatCapacity<Element, Scv, ElemSol, SolidState>())(problem.spatialParams())
+ )::value, Scalar>
+ {
+ return SolidSystem::heatCapacity(solidState);
+ }
+
+ /*!
+ * \brief get the solid density in an scv
+ * \param elemSol the element solution vector
+ * \param problem the problem to solve
+ * \param element the element (codim-0-entity) the scv belongs to
+ * \param scv the sub control volume
+ * \param solidState the solid state
+ * \note this gets selected if the user uses the solidsystem / solidstate interface
+ */
+ template<class ElemSol, class Problem, class Element, class Scv>
+ auto solidDensity_(const ElemSol& elemSol,
+ const Problem& problem,
+ const Element& element,
+ const Scv& scv,
+ const SolidState& solidState)
+ -> typename std::enable_if_t<!decltype(
+ isValid(Detail::hasSolidDensity<Element, Scv, ElemSol, SolidState>())(problem.spatialParams())
+ )::value, Scalar>
+ {
+ return SolidSystem::density(solidState);
+ }
+
+ /*!
+ * \brief get the solid's thermal conductivity in an scv
+ * \param elemSol the element solution vector
+ * \param problem the problem to solve
+ * \param element the element (codim-0-entity) the scv belongs to
+ * \param scv the sub control volume
+ * \param solidState the solid state
+ * \note this gets selected if the user uses the solidsystem / solidstate interface
+ */
+ template<class ElemSol, class Problem, class Element, class Scv>
+ auto solidThermalConductivity_(const ElemSol& elemSol,
+ const Problem& problem,
+ const Element& element,
+ const Scv& scv,
+ const SolidState& solidState)
+ -> typename std::enable_if_t<!decltype(
+ isValid(Detail::hasSolidThermalConductivity<Element, Scv, ElemSol, SolidState>())(problem.spatialParams())
+ )::value, Scalar>
+ {
+ return SolidSystem::thermalConductivity(solidState);
+ }
+
+ // \}
+
+ /*!
+ * \name Access functions for the simple nonisothermal spatial params interface in
+ * combination with an InertSolidPhase as solid system
+ */
+ // \{
+
+ /*!
+ * \brief get the solid heat capacity in an scv
+ * \param elemSol the element solution vector
+ * \param problem the problem to solve
+ * \param element the element (codim-0-entity) the scv belongs to
+ * \param scv the sub control volume
+ * \param solidState the solid state
+ * \note this gets selected if the user uses the simple spatial params interface in
+ * combination with an InertSolidPhase as solid system
+ */
+ template<class ElemSol, class Problem, class Element, class Scv>
+ auto solidHeatCapacity_(const ElemSol& elemSol,
+ const Problem& problem,
+ const Element& element,
+ const Scv& scv,
+ const SolidState& solidState)
+ -> typename std::enable_if_t<decltype(
+ isValid(Detail::hasSolidHeatCapacity<Element, Scv, ElemSol, SolidState>())(problem.spatialParams())
+ )::value, Scalar>
+ {
+ static_assert(Detail::isInertSolidPhase<SolidSystem>::value,
+ "solidHeatCapacity can only be overwritten in the spatial params when the solid system is a simple InertSolidPhase\n"
+ "If you select a proper solid system, the solid heat capacity will be computed as stated in the solid system!");
+ return problem.spatialParams().solidHeatCapacity(element, scv, elemSol, solidState);
+ }
+
+ /*!
+ * \brief get the solid density in an scv
+ * \param elemSol the element solution vector
+ * \param problem the problem to solve
+ * \param element the element (codim-0-entity) the scv belongs to
+ * \param scv the sub control volume
+ * \param solidState the solid state
+ * \note this gets selected if the user uses the simple spatial params interface in
+ * combination with an InertSolidPhase as solid system
+ */
+ template<class ElemSol, class Problem, class Element, class Scv>
+ auto solidDensity_(const ElemSol& elemSol,
+ const Problem& problem,
+ const Element& element,
+ const Scv& scv,
+ const SolidState& solidState)
+ -> typename std::enable_if_t<decltype(
+ isValid(Detail::hasSolidDensity<Element, Scv, ElemSol, SolidState>())(problem.spatialParams())
+ )::value, Scalar>
+ {
+ static_assert(Detail::isInertSolidPhase<SolidSystem>::value,
+ "solidDensity can only be overwritten in the spatial params when the solid system is a simple InertSolidPhase\n"
+ "If you select a proper solid system, the solid density will be computed as stated in the solid system!");
+ return problem.spatialParams().solidDensity(element, scv, elemSol, solidState);
+ }
+
+ /*!
+ * \brief get the solid's heat capacity in an scv
+ * \param elemSol the element solution vector
+ * \param problem the problem to solve
+ * \param element the element (codim-0-entity) the scv belongs to
+ * \param scv the sub control volume
+ * \param solidState the solid state
+ * \note this gets selected if the user uses the simple spatial params interface in
+ * combination with an InertSolidPhase as solid system
+ */
+ template<class ElemSol, class Problem, class Element, class Scv>
+ auto solidThermalConductivity_(const ElemSol& elemSol,
+ const Problem& problem,
+ const Element& element,
+ const Scv& scv,
+ const SolidState& solidState)
+ -> typename std::enable_if_t<decltype(
+ isValid(Detail::hasSolidThermalConductivity<Element, Scv, ElemSol, SolidState>())(problem.spatialParams())
+ )::value, Scalar>
+ {
+ static_assert(Detail::isInertSolidPhase<SolidSystem>::value,
+ "solidThermalConductivity can only be overwritten in the spatial params when the solid system is a simple InertSolidPhase\n"
+ "If you select a proper solid system, the solid thermal conductivity will be computed as stated in the solid system!");
+ return problem.spatialParams().solidThermalConductivity(element, scv, elemSol, solidState);
+ }
+
+ // \}
+
};
} // end namespace Dumux
diff --git a/test/porousmediumflow/3p3c/implicit/columnxylol/params.input b/test/porousmediumflow/3p3c/implicit/columnxylol/params.input
index f383e20058..f502ce89f7 100644
--- a/test/porousmediumflow/3p3c/implicit/columnxylol/params.input
+++ b/test/porousmediumflow/3p3c/implicit/columnxylol/params.input
@@ -13,12 +13,8 @@ Name = columnxylol # name passed to the output routines
[Assembly]
NumericDifferenceMethod = 0 # -1 backward differences, 0: central differences, +1: forward differences
-[1.Component]
+[Component]
SolidDensity = 2650
SolidThermalConductivity = 2.8
-SolidHeatCapacity = 850
-
-[2.Component]
-SolidDensity = 2650
-SolidThermalConductivity = 2.8
-SolidHeatCapacity = 84000
+SolidHeatCapacityFine = 850
+SolidHeatCapacityCoarse = 84000
diff --git a/test/porousmediumflow/3p3c/implicit/columnxylol/problem.hh b/test/porousmediumflow/3p3c/implicit/columnxylol/problem.hh
index 9769d9117d..22fbccaa6f 100644
--- a/test/porousmediumflow/3p3c/implicit/columnxylol/problem.hh
+++ b/test/porousmediumflow/3p3c/implicit/columnxylol/problem.hh
@@ -28,8 +28,8 @@
#include <dune/grid/yaspgrid.hh>
#include <dumux/material/fluidsystems/h2oairxylene.hh>
-#include <dumux/material/solidstates/compositionalsolidstate.hh>
-#include <dumux/material/solidsystems/compositionalsolidphase.hh>
+#include <dumux/material/solidstates/inertsolidstate.hh>
+#include <dumux/material/solidsystems/inertsolidphase.hh>
#include <dumux/material/components/constant.hh>
#include <dumux/discretization/cellcentered/tpfa/properties.hh>
@@ -75,13 +75,10 @@ template<class TypeTag>
struct SolidSystem<TypeTag, TTag::Column>
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
- using ComponentOne = Dumux::Components::Constant<1, Scalar>;
- using ComponentTwo = Dumux::Components::Constant<2, Scalar>;
- static constexpr int numInertComponents = 2;
- using type = SolidSystems::CompositionalSolidPhase<Scalar, ComponentOne, ComponentTwo, numInertComponents>;
+ using Component = Dumux::Components::Constant<1, Scalar>;
+ using type = SolidSystems::InertSolidPhase<Scalar, Component>;
};
-
//! The two-phase model uses the immiscible fluid state
template<class TypeTag>
struct SolidState<TypeTag, TTag::Column>
@@ -90,7 +87,7 @@ private:
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using SolidSystem = GetPropType<TypeTag, Properties::SolidSystem>;
public:
- using type = CompositionalSolidState<Scalar, SolidSystem>;
+ using type = InertSolidState<Scalar, SolidSystem>;
};
// Set the spatial parameters
diff --git a/test/porousmediumflow/3p3c/implicit/columnxylol/spatialparams.hh b/test/porousmediumflow/3p3c/implicit/columnxylol/spatialparams.hh
index 2187574c32..c942f2a062 100644
--- a/test/porousmediumflow/3p3c/implicit/columnxylol/spatialparams.hh
+++ b/test/porousmediumflow/3p3c/implicit/columnxylol/spatialparams.hh
@@ -70,12 +70,11 @@ public:
coarseK_ = 1.4e-8;
// porosities
- finePorosity_ = 0.46;
- coarsePorosity_ = 0.46;
+ porosity_ = 0.46;
// specific heat capacities
- fineHeatCap_ = 850.;
- coarseHeatCap_ = 84000.;
+ fineHeatCap_ = getParam<Scalar>("Component.SolidHeatCapacityFine", 850.0);
+ coarseHeatCap_ = getParam<Scalar>("Component.SolidHeatCapacityCoarse", 84000.0);
// residual saturations
fineMaterialParams_.setSwr(0.12);
@@ -121,31 +120,13 @@ public:
return coarseK_;
}
- /*!
- * \brief Define the porosity \f$[-]\f$ of the spatial parameters
- *
- * \param element The current element
- * \param scv The sub-control volume inside the element.
- * \param elemSol The solution at the dofs connected to the element.
- */
- template<class SolidSystem>
- Scalar inertVolumeFractionAtPos(const GlobalPosition& globalPos,
- int compIdx) const
+ /*! \brief Define the porosity in [-].
+ *
+ * \param globalPos The global position where we evaluate
+ */
+ Scalar porosityAtPos(const GlobalPosition& globalPos) const
{
- if (compIdx == SolidSystem::comp0Idx)
- {
- if (isFineMaterial_(globalPos))
- return 1-finePorosity_;
- else
- return 0;
- }
- else
- {
- if (isFineMaterial_(globalPos))
- return 0;
- else
- return 1-coarsePorosity_;
- }
+ return porosity_;
}
/*!
@@ -168,6 +149,28 @@ public:
return coarseMaterialParams_;
}
+ /*!
+ * \brief User-defined solid heat capacity.
+ *
+ * \param element The current element
+ * \param scv The sub-control volume inside the element.
+ * \param elemSol The solution at the dofs connected to the element.
+ * \param solidState The solid state
+ * \return the solid heat capacity
+ */
+ template <class ElementSolution, class SolidState>
+ Scalar solidHeatCapacity(const Element& element,
+ const SubControlVolume& scv,
+ const ElementSolution& elemSol,
+ const SolidState& solidState) const
+ {
+ const auto& globalPos = scv.dofPosition();
+ if (isFineMaterial_(globalPos))
+ return fineHeatCap_;
+ else
+ return coarseHeatCap_;
+ }
+
private:
bool isFineMaterial_(const GlobalPosition &globalPos) const
{
@@ -177,8 +180,7 @@ private:
Scalar fineK_;
Scalar coarseK_;
- Scalar finePorosity_;
- Scalar coarsePorosity_;
+ Scalar porosity_;
Scalar fineHeatCap_;
Scalar coarseHeatCap_;
--
GitLab