diff --git a/dumux/material/fluidmatrixinteractions/CMakeLists.txt b/dumux/material/fluidmatrixinteractions/CMakeLists.txt
index 07a9487406b622961e4a1c71ccf3bde1e2615e9f..3b4e7d43340bcb73cea04253db035a45bb7aeae0 100644
--- a/dumux/material/fluidmatrixinteractions/CMakeLists.txt
+++ b/dumux/material/fluidmatrixinteractions/CMakeLists.txt
@@ -12,5 +12,4 @@ diffusivitymillingtonquirk.hh
permeabilitykozenycarman.hh
permeabilityrutqvisttsang.hh
porosityprecipitation.hh
-porosityreactivebed.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/material/fluidmatrixinteractions)
diff --git a/dumux/material/fluidmatrixinteractions/permeabilitykozenycarman.hh b/dumux/material/fluidmatrixinteractions/permeabilitykozenycarman.hh
index 7bc259d984993c2bea9671934cb353050dc859ca..5864e210abfb3dc312f1df878b073df11b437105 100644
--- a/dumux/material/fluidmatrixinteractions/permeabilitykozenycarman.hh
+++ b/dumux/material/fluidmatrixinteractions/permeabilitykozenycarman.hh
@@ -35,71 +35,27 @@ namespace Dumux {
* \brief The Kozeny-Carman relationship for the calculation of a porosity-dependent permeability.
* When the porosity is implemented as solution-independent, using this relationship for the
* permeability leads to unnecessary overhead.
+ *
+ * \tparam PermeabilityType The type used for the intrinsic permeability
*/
-template<class TypeTag>
+template<class PermeabilityType>
class PermeabilityKozenyCarman
{
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using SpatialParams = typename GET_PROP_TYPE(TypeTag, SpatialParams);
- using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
- using SubControlVolume = typename FVElementGeometry::SubControlVolume;
-
- static constexpr int dim = GridView::dimension;
- static constexpr int dimWorld = GridView::dimensionworld;
- using Element = typename GridView::template Codim<0>::Entity;
- using Tensor = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
-
- using PermType = typename SpatialParams::PermeabilityType;
-
public:
-
- /*!
- * \brief The initial parameter distribution.
- * \param spatialParams the spatial parameters
- */
- void init(const SpatialParams& spatialParams)
- {
- spatialParamsPtr_ = &spatialParams;
- }
-
/*!
* \brief calculates the permeability for a given sub-control volume
- * \param element element
- * \param elemSol the element solution
- * \param scv sub control volume
+ * \param refPerm Reference permeability before porosity changes
+ * \param refPoro The poro corresponding to the reference permeability
+ * \param poro The porosity for which permeability is to be evaluated
*/
- template<class ElementSolution>
- PermType evaluatePermeability(const Element& element,
- const SubControlVolume& scv,
- const ElementSolution& elemSol) const
+ template<class Scalar>
+ PermeabilityType evaluatePermeability(PermeabilityType refPerm, Scalar refPoro, Scalar poro) const
{
- auto refPoro = spatialParams().referencePorosity(element, scv);
- auto poro = spatialParams().porosity(element, scv, elemSol);
-
using std::pow;
auto factor = pow((1.0 - refPoro)/(1.0 - poro), 2) * pow(poro/refPoro, 3);
- return applyFactorToPermeability_(spatialParams().referencePermeability(element, scv), factor);
- }
-
-private:
- const SpatialParams& spatialParams() const
- { return *spatialParamsPtr_; }
-
- Scalar applyFactorToPermeability_(Scalar k, Scalar factor) const
- { return k*factor; }
-
- Tensor applyFactorToPermeability_(Tensor K, Scalar factor) const
- {
- Tensor result(K);
- for (int i = 0; i < dim; ++i)
- result[i][i] *= factor;
- return result;
+ refPerm *= factor;
+ return refPerm;
}
-
- const SpatialParams* spatialParamsPtr_;
};
} // namespace Dumux
diff --git a/dumux/material/fluidmatrixinteractions/porosityprecipitation.hh b/dumux/material/fluidmatrixinteractions/porosityprecipitation.hh
index 68e3491e61d81c30effa9b7cc0c58ea1be9325a8..cdf9af827f3a265dc6893497a7c07b9f0d5a1a2e 100644
--- a/dumux/material/fluidmatrixinteractions/porosityprecipitation.hh
+++ b/dumux/material/fluidmatrixinteractions/porosityprecipitation.hh
@@ -24,7 +24,6 @@
#ifndef DUMUX_POROSITY_PRECIPITATION_HH
#define DUMUX_POROSITY_PRECIPITATION_HH
-#include <dumux/common/properties.hh>
#include <dumux/discretization/evalsolution.hh>
namespace Dumux {
@@ -32,58 +31,39 @@ namespace Dumux {
/*!
* \ingroup Fluidmatrixinteractions
* \brief Calculates the porosity depending on the volume fractions of precipitated minerals.
+ *
+ * \tparam Scalar The type used for scalar values
+ * \numComp The number of components in the fluid phases
+ * \numSolidPhases The number of precipitating solid phases
*/
-template<class TypeTag>
+template<class Scalar, int numComp, int numSolidPhases>
class PorosityPrecipitation
{
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using SpatialParams = typename GET_PROP_TYPE(TypeTag, SpatialParams);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
- using SubControlVolume = typename FVElementGeometry::SubControlVolume;
-
- static const int dim = GridView::dimension;
- static const int dimWorld = GridView::dimensionworld;
- static const int numComponents = GET_PROP_TYPE(TypeTag, ModelTraits)::numComponents();
- static const int numSolidPhases = GET_PROP_TYPE(TypeTag, ModelTraits)::numSPhases();
-
- using Element = typename GridView::template Codim<0>:: Entity;
-
public:
- void init(const SpatialParams& spatialParams)
- {
- spatialParamsPtr_ = &spatialParams;
- }
-
/*!
* \brief calculates the porosity in a sub-control volume
* \param element element
* \param elemSol the element solution
* \param scv sub control volume
+ * \param refPoro The solid matrix porosity without precipitates
+ * \param minPoro A minimum porosity value
*/
- template<class ElementSolution>
+ template<class Element, class SubControlVolume, class ElemSol>
Scalar evaluatePorosity(const Element& element,
const SubControlVolume& scv,
- const ElementSolution& elemSol) const
+ const ElemSol& elemSol,
+ Scalar refPoro,
+ Scalar minPoro = 0.0) const
{
auto priVars = evalSolution(element, element.geometry(), elemSol, scv.center());
Scalar sumPrecipitates = 0.0;
for (unsigned int solidPhaseIdx = 0; solidPhaseIdx < numSolidPhases; ++solidPhaseIdx)
- sumPrecipitates += priVars[numComponents + solidPhaseIdx];
-
- auto minPoro = spatialParams_().minPorosity(element, scv);
+ sumPrecipitates += priVars[numComp + solidPhaseIdx];
using std::max;
- return max(minPoro, spatialParams_().referencePorosity(element, scv) - sumPrecipitates);
+ return max(minPoro, refPoro - sumPrecipitates);
}
-
-private:
- const SpatialParams& spatialParams_() const
- { return *spatialParamsPtr_; }
-
- const SpatialParams* spatialParamsPtr_;
};
} // namespace Dumux
diff --git a/dumux/material/fluidmatrixinteractions/porosityreactivebed.hh b/dumux/material/fluidmatrixinteractions/porosityreactivebed.hh
deleted file mode 100644
index df985ae6e48fd213c7f87dcd54ef9d2c32603223..0000000000000000000000000000000000000000
--- a/dumux/material/fluidmatrixinteractions/porosityreactivebed.hh
+++ /dev/null
@@ -1,89 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * See the file COPYING for full copying permissions. *
- * *
- * 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
- * \ingroup Fluidmatrixinteractions
- * \brief Class for the evaluation of the porosity subject to precipitation.
- */
-#ifndef DUMUX_POROSITY_REACTIVE_BED_HH
-#define DUMUX_POROSITY_REACTIVE_BED_HH
-
-#include <dumux/common/properties.hh>
-#include <dumux/discretization/evalsolution.hh>
-
-namespace Dumux {
-
-/*!
- * \ingroup Fluidmatrixinteractions
- * \brief Calculates the porosity depeding on the volume fractions of different solid species.
- */
-template<class TypeTag>
-class PorosityReactiveBed
-{
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using SpatialParams = typename GET_PROP_TYPE(TypeTag, SpatialParams);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
- using SubControlVolume = typename FVElementGeometry::SubControlVolume;
-
- static const int dim = GridView::dimension;
- static const int dimWorld = GridView::dimensionworld;
- static const int numComponents = GET_PROP_TYPE(TypeTag, ModelTraits)::numComponents();
- static const int numSolidPhases = GET_PROP_TYPE(TypeTag, ModelTraits)::numSPhases();
-
- using Element = typename GridView::template Codim<0>:: Entity;
-
-public:
- void init(const SpatialParams& spatialParams)
- {
- spatialParamsPtr_ = &spatialParams;
- }
-
- /*!
- * \brief Calculates the porosity in a sub-control volume.
- *
- * \param element element
- * \param elemSol the element solution
- * \param scv sub control volume
- */
- template<class ElementSolution>
- Scalar evaluatePorosity(const Element& element,
- const SubControlVolume& scv,
- const ElementSolution& elemSol) const
- {
- auto priVars = evalSolution(element, element.geometry(), elemSol, scv.center());
-
- Scalar sumPrecipitates = 0.0;
- for (unsigned int solidPhaseIdx = 0; solidPhaseIdx < numSolidPhases; ++solidPhaseIdx)
- sumPrecipitates += priVars[numComponents + solidPhaseIdx];
-
- return (1 - sumPrecipitates);
- }
-
-private:
- const SpatialParams& spatialParams_() const
- { return *spatialParamsPtr_; }
-
- const SpatialParams* spatialParamsPtr_;
-};
-
-} // namespace Dumux
-
-#endif
diff --git a/test/porousmediumflow/1pncmin/implicit/thermochemspatialparams.hh b/test/porousmediumflow/1pncmin/implicit/thermochemspatialparams.hh
index 561a578d20e925d14b440eff2dcf4302953ced86..f5d23d4f3992797c52e34ff6e3d18839eb9558f8 100644
--- a/test/porousmediumflow/1pncmin/implicit/thermochemspatialparams.hh
+++ b/test/porousmediumflow/1pncmin/implicit/thermochemspatialparams.hh
@@ -28,8 +28,6 @@
#include <dumux/material/spatialparams/fv1p.hh>
-#include <dumux/material/fluidmatrixinteractions/porosityreactivebed.hh>
-#include <dumux/material/fluidmatrixinteractions/permeabilitykozenycarman.hh>
#include <dumux/material/fluidmatrixinteractions/mineralization/effectivesoliddensity.hh>
#include <dumux/material/fluidmatrixinteractions/mineralization/effectivesolidheatcapacity.hh>
@@ -79,8 +77,6 @@ class ThermoChemSpatialParams : public FVSpatialParamsOneP<TypeTag>
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using Element = typename GridView::template Codim<0>::Entity;
- using PorosityLaw = PorosityReactiveBed<TypeTag>;
- using PermeabilityLaw = PermeabilityKozenyCarman<TypeTag>;
using EffectiveSolidRho = EffectiveSolidDensity<TypeTag>;
using EffectiveSolidCp = EffectiveSolidHeatCapacity<TypeTag>;
@@ -102,32 +98,10 @@ public:
cp_[cPhaseIdx-numPhases] = 934; //[J/kgK] heat capacity of CaO (see Shao et al. 2014)
cp_[hPhaseIdx-numPhases] = 1530; //[J/kgK] heat capacity of Ca(OH)_2 (see Shao et al. 2014)
eps_ = 1e-6;
- poroLaw_.init(*this);
- permLaw_.init(*this);
effSolRho_.init(*this);
effSolCp_.init(*this);
}
- /*!
- * \brief Define the reference permeability \f$[m^2]\f$ distribution
- *
- * \param element The finite element
- * \param scv The sub-control volume
- */
- Scalar referencePermeability(const Element& element, const SubControlVolume &scv) const
- { return 8.53e-12; }
-
- /*!
- * \brief Define the reference porosity \f$[-]\f$ distribution
- *
- * \param element The finite element
- * \param scv The sub-control volume
- */
- Scalar referencePorosity(const Element& element, const SubControlVolume &scv) const
- {
- return 0.8;
- }
-
/*! Intrinsic permeability tensor K \f$[m^2]\f$ depending
* on the position in the domain
*
@@ -141,19 +115,7 @@ public:
Scalar permeability(const Element& element,
const SubControlVolume& scv,
const ElementSolution& elemSol) const
- { return permLaw_.evaluatePermeability(element, scv, elemSol); }
-
- /*!
- * \brief Define the minimum porosity \f$[-]\f$ distribution
- *
- * \param element The finite element
- * \param scv The sub-control volume
- */
- Scalar minPorosity(const Element& element, const SubControlVolume &scv) const
- {
-// return 0.604; //intrinsic porosity of CaO2H2 (see Nagel et al. 2014)
- return 0.8;
- }
+ { return 8.53e-12; }
/*!
* \brief Define the minimum porosity \f$[-]\f$ after clogging caused by mineralization
@@ -166,10 +128,7 @@ public:
Scalar porosity(const Element& element,
const SubControlVolume& scv,
const ElementSolution& elemSol) const
- {
-// return poroLaw_.evaluatePorosity(element, scv, elemSol);
- return 0.8;
- }
+ { return 0.8; }
/*!
* \brief Returns the average heat capacity \f$[J / (kg K)]\f$ of solid phases.
@@ -248,8 +207,6 @@ private:
Scalar lambdaSolid_;
std::array<Scalar, numSPhases> rho_;
std::array<Scalar, numSPhases> cp_;
- PorosityLaw poroLaw_;
- PermeabilityLaw permLaw_;
EffectiveSolidRho effSolRho_;
EffectiveSolidCp effSolCp_;
};
diff --git a/test/porousmediumflow/2pncmin/implicit/dissolutionspatialparams.hh b/test/porousmediumflow/2pncmin/implicit/dissolutionspatialparams.hh
index 5ad810c65607c1cd52433f5faee5460dcbb53328..0cd4526dcddbc812e43cd09a43a57be926628907 100644
--- a/test/porousmediumflow/2pncmin/implicit/dissolutionspatialparams.hh
+++ b/test/porousmediumflow/2pncmin/implicit/dissolutionspatialparams.hh
@@ -83,9 +83,6 @@ class DissolutionSpatialparams : public FVSpatialParams<TypeTag>
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using Element = typename GridView::template Codim<0>::Entity;
- using PorosityLaw = PorosityPrecipitation<TypeTag>;
- using PermeabilityLaw = PermeabilityKozenyCarman<TypeTag>;
-
public:
// type used for the permeability (i.e. tensor or scalar)
using PermeabilityType = Scalar;
@@ -108,10 +105,6 @@ public:
// parameters of Brooks & Corey Law
materialParams_.setPe(pEntry1_);
materialParams_.setLambda(bcLambda1_);
-
- //! Initialize the parameter laws
- poroLaw_.init(*this);
- permLaw_.init(*this);
}
/*!
@@ -145,17 +138,7 @@ public:
Scalar porosity(const Element& element,
const SubControlVolume& scv,
const ElementSolution& elemSol) const
- { return poroLaw_.evaluatePorosity(element, scv, elemSol); }
-
- /*!
- * \brief Define the reference permeability \f$[m^2]\f$ distribution
- *
- * \param element The finite element
- * \param scv The sub-control volume
- */
- PermeabilityType referencePermeability(const Element& element, const SubControlVolume &scv) const
- { return referencePermeability_; }
-
+ { return poroLaw_.evaluatePorosity(element, scv, elemSol, referencePorosity_, 1e-5); }
/*! Intrinsic permeability tensor K \f$[m^2]\f$ depending
* on the position in the domain
@@ -169,7 +152,10 @@ public:
PermeabilityType permeability(const Element& element,
const SubControlVolume& scv,
const ElementSolution& elemSol) const
- { return permLaw_.evaluatePermeability(element, scv, elemSol); }
+ {
+ const auto poro = porosity(element, scv, elemSol);
+ return permLaw_.evaluatePermeability(referencePermeability_, referencePorosity_, poro);
+ }
Scalar solidity(const SubControlVolume &scv) const
{ return 1.0 - porosityAtPos(scv.center()); }
@@ -188,8 +174,11 @@ private:
MaterialLawParams materialParams_;
+ using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
+ using PorosityLaw = PorosityPrecipitation<Scalar, ModelTraits::numComponents(), ModelTraits::numSPhases()>;
PorosityLaw poroLaw_;
- PermeabilityLaw permLaw_;
+ PermeabilityKozenyCarman<PermeabilityType> permLaw_;
+
Scalar solubilityLimit_;
Scalar referencePorosity_;
PermeabilityType referencePermeability_ = 0.0;