From ae450cb3ccd9c57329db5f8a3bc20c847d56c7bb Mon Sep 17 00:00:00 2001
From: Katharina Heck <katharina.heck@iws.uni-stuttgart.de>
Date: Thu, 26 Jul 2018 16:45:28 +0200
Subject: [PATCH] [feature] add spatialparams with base stuff for
nonequilibrium
---
.../spatialparams/fvnonequilibrium.hh | 210 ++++++++++++++++++
.../nonequilibrium/volumevariables.hh | 107 +++++++--
.../mpnc/implicit/kinetic/spatialparams.hh | 71 ++----
.../thermalnonequilibrium/spatialparams.hh | 69 ++----
4 files changed, 338 insertions(+), 119 deletions(-)
create mode 100644 dumux/material/spatialparams/fvnonequilibrium.hh
diff --git a/dumux/material/spatialparams/fvnonequilibrium.hh b/dumux/material/spatialparams/fvnonequilibrium.hh
new file mode 100644
index 0000000000..98a0eef7f5
--- /dev/null
+++ b/dumux/material/spatialparams/fvnonequilibrium.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:
+/*****************************************************************************
+ * 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 SpatialParameters
+ * \brief base class for spatialparameters dealing with thermal and chemical nonequilibrium
+ *
+ */
+#ifndef DUMUX_FV_SPATIALPARAMS_NONEQUILIBRIUM_HH
+#define DUMUX_FV_SPATIALPARAMS_NONEQUILIBRIUM_HH
+
+#include <dumux/material/spatialparams/fv.hh>
+
+// material laws for interfacial area
+#include <dumux/material/fluidmatrixinteractions/2pia/efftoabslawia.hh>
+#include <dumux/material/fluidmatrixinteractions/2pia/awnsurfacepolynomial2ndorder.hh>
+#include <dumux/material/fluidmatrixinteractions/2pia/awnsurfacepcmaxfct.hh>
+#include <dumux/material/fluidmatrixinteractions/2pia/awnsurfaceexpswpcto3.hh>
+
+namespace Dumux {
+
+/**
+ * \brief Definition of the spatial parameters for nonequilibrium
+ */
+template<class FVGridGeometry, class Scalar, class Implementation>
+class FVNonEquilibriumSpatialParams
+: public FVSpatialParams<FVGridGeometry,
+ Scalar,
+ Implementation>
+{
+
+ using ParentType = FVSpatialParams<FVGridGeometry, Scalar, Implementation>;
+ using GridView = typename FVGridGeometry::GridView;
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolume = typename FVGridGeometry::SubControlVolume;
+ using Element = typename GridView::template Codim<0>::Entity;
+
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+public:
+ //! export the types used for interfacial area calculations
+
+ using AwnSurfaceParams = Scalar;
+ using AwsSurfaceParams = Scalar;
+ using AnsSurfaceParams = Scalar;
+
+ FVNonEquilibriumSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+ : ParentType(fvGridGeometry)
+ { }
+
+ /*!\brief Return a reference to the container object for the
+ * parametrization of the surface between wetting and non-Wetting phase.
+ *
+ * The position is determined based on the coordinate of
+ * the vertex belonging to the considered sub control volume.
+ * \param element The finite element
+ * \param fvGeometry The finite volume geometry
+ * \param scvIdx The local index of the sub-control volume */
+ template<class ElementSolution>
+ AwnSurfaceParams aWettingNonWettingSurfaceParams(const Element &element,
+ const SubControlVolume &scv,
+ const ElementSolution &elemSol) const
+ {
+ DUNE_THROW(Dune::InvalidStateException,
+ "The spatial parameters do not provide a aWettingNonWettingSurfaceParams() method.");
+ }
+
+ /*!\brief Return a reference to the container object for the
+ * parametrization of the surface between non-Wetting and solid phase.
+ *
+ * The position is determined based on the coordinate of
+ * the vertex belonging to the considered sub control volume.
+ * \param element The finite element
+ * \param fvGeometry The finite volume geometry
+ * \param scvIdx The local index of the sub-control volume */
+ template<class ElementSolution>
+ AnsSurfaceParams aNonWettingSolidSurfaceParams(const Element &element,
+ const SubControlVolume &scv,
+ const ElementSolution &elemSol) const
+ {
+ DUNE_THROW(Dune::InvalidStateException,
+ "The spatial parameters do not provide a aNonWettingSolidSurfaceParams() method.");
+ }
+
+ /*!\brief Return a reference to the container object for the
+ * parametrization of the surface between non-Wetting and solid phase.
+ *
+ * The position is determined based on the coordinate of
+ * the vertex belonging to the considered sub control volume.
+ * \param element The finite element
+ * \param fvGeometry The finite volume geometry
+ * \param scvIdx The local index of the sub-control volume */
+ template<class ElementSolution>
+ AwsSurfaceParams aWettingSolidSurfaceParams(const Element &element,
+ const SubControlVolume &scv,
+ const ElementSolution &elemSol) const
+ {
+ DUNE_THROW(Dune::InvalidStateException,
+ "The spatial parameters do not provide a aWettingSolidSurfaceParams() method.");
+ }
+
+ /*!\brief Return the maximum capillary pressure for the given pc-Sw curve
+ *
+ * Of course physically there is no such thing as a maximum capillary pressure.
+ * The parametrization (VG/BC) goes to infinity and physically there is only one pressure
+ * for single phase conditions.
+ * Here, this is used for fitting the interfacial area surface: the capillary pressure,
+ * where the interfacial area is zero.
+ * Technically this value is obtained as the capillary pressure of saturation zero.
+ * This value of course only exists for the case of a regularized pc-Sw relation.
+ * \param element The finite element
+ * \param fvGeometry The finite volume geometry
+ * \param scvIdx The local index of the sub-control volume */
+ template<class ElementSolution>
+ const Scalar pcMax(const Element &element,
+ const SubControlVolume &scv,
+ const ElementSolution &elemSol) const
+ { DUNE_THROW(Dune::InvalidStateException,
+ "The spatial parameters do not provide a aNonWettingSolidSurfaceParams() method.");
+ }
+
+
+ /*!\brief Return the characteristic length for the mass transfer.
+ *
+ * The position is determined based on the coordinate of
+ * the vertex belonging to the considered sub controle volume.
+ * \param element The finite element
+ * \param fvGeometry The finite volume geometry
+ * \param scvIdx The local index of the sub control volume */
+ template<class ElementSolution>
+ const Scalar characteristicLength(const Element & element,
+ const SubControlVolume &scv,
+ const ElementSolution &elemSol) const
+
+ { return this->asImp_().characteristicLengthAtPos(scv.dofPosition()); }
+
+ /*!\brief Return the characteristic length for the mass transfer.
+ * \param globalPos The position in global coordinates.*/
+ const Scalar characteristicLengthAtPos(const GlobalPosition & globalPos) const
+ {
+ DUNE_THROW(Dune::InvalidStateException,
+ "The spatial parameters do not provide "
+ "a characteristicLengthAtPos() method.");
+ }
+
+ /*!\brief Return the pre factor the the energy transfer
+ *
+ * The position is determined based on the coordinate of
+ * the vertex belonging to the considered sub controle volume.
+ * \param element The finite element
+ * \param fvGeometry The finite volume geometry
+ * \param scvIdx The local index of the sub control volume */
+ template<class ElementSolution>
+ const Scalar factorEnergyTransfer(const Element &element,
+ const SubControlVolume &scv,
+ const ElementSolution &elemSol) const
+ { return this->asImp_().factorEnergyTransferAtPos(scv.dofPosition()); }
+
+ /*!\brief Return the pre factor the the energy transfer
+ * \param globalPos The position in global coordinates.*/
+ const Scalar factorEnergyTransferAtPos(const GlobalPosition & globalPos) const
+ {
+ DUNE_THROW(Dune::InvalidStateException,
+ "The spatial parameters do not provide "
+ "a factorEnergyTransferAtPos() method.");
+ }
+
+ /*!\brief Return the pre factor the the mass transfer
+ *
+ * The position is determined based on the coordinate of
+ * the vertex belonging to the considered sub controle volume.
+ * \param element The finite element
+ * \param fvGeometry The finite volume geometry
+ * \param scvIdx The local index of the sub control volume */
+ template<class ElementSolution>
+ const Scalar factorMassTransfer(const Element &element,
+ const SubControlVolume &scv,
+ const ElementSolution &elemSol) const
+ { return this->asImp_().factorMassTransferAtPos(scv.dofPosition()); }
+
+
+ /*!\brief Return the pre factor the the mass transfer
+ * \param globalPos The position in global coordinates.*/
+ const Scalar factorMassTransferAtPos(const GlobalPosition & globalPos) const
+ {
+ DUNE_THROW(Dune::InvalidStateException,
+ "The spatial parameters do not provide "
+ "a factorMassTransferAtPos() method.");
+ }
+};
+
+}
+
+#endif // GUARDIAN
diff --git a/dumux/porousmediumflow/nonequilibrium/volumevariables.hh b/dumux/porousmediumflow/nonequilibrium/volumevariables.hh
index 09c73d50a4..af5d888243 100644
--- a/dumux/porousmediumflow/nonequilibrium/volumevariables.hh
+++ b/dumux/porousmediumflow/nonequilibrium/volumevariables.hh
@@ -68,10 +68,14 @@ class NonEquilibriumVolumeVariablesImplementation< Traits,
using ModelTraits = typename Traits::ModelTraits;
using Indices = typename ModelTraits::Indices;
+ using FS = typename Traits::FluidSystem;
static constexpr auto numEnergyEqFluid = ModelTraits::numEnergyEqFluid();
static constexpr auto numEnergyEqSolid = ModelTraits::numEnergyEqSolid();
- static_assert((numEnergyEqFluid < 2), "This model is a specialization for a energy transfer of a fluid mixture and a solid");
+ static constexpr auto phase0Idx = FS::phase0Idx;
+ static constexpr auto phase1Idx = FS::phase1Idx;
+ static constexpr auto sPhaseIdx = FS::numPhases;
+ static_assert((numEnergyEqFluid < 3), "this model only has interfacial area relationships for maximum 2 fluids with one solid");
using DimLessNum = DimensionlessNumbers<Scalar>;
@@ -97,7 +101,9 @@ public:
ParameterCache paramCache;
paramCache.updateAll(this->fluidState());
- updateInterfacialArea(elemSol, this->fluidState(), paramCache, problem, element, scv);
+ updateDimLessNumbers(elemSol, this->fluidState(), paramCache, problem, element, scv);
+ if (numEnergyEqFluid == 2)
+ updateInterfacialArea(elemSol, this->fluidState(), paramCache, problem, element, scv);
}
/*!
@@ -111,7 +117,7 @@ public:
* \param scv The sub-control volume
*/
template<class ElemSol, class Problem, class Element, class Scv>
- void updateInterfacialArea(const ElemSol& elemSol,
+ void updateDimLessNumbers(const ElemSol& elemSol,
const FluidState& fluidState,
const ParameterCache& paramCache,
const Problem& problem,
@@ -136,7 +142,7 @@ public:
const auto thermalConductivity = FluidSystem::thermalConductivity(fluidState, paramCache, phaseIdx);
const auto porosity = this->porosity();
- reynoldsNumber_[phaseIdx] = DimLessNum::reynoldsNumber(darcyMagVelocity, characteristicLength_, kinematicViscosity);
+ reynoldsNumber_[phaseIdx] = DimLessNum::reynoldsNumber(darcyMagVelocity, characteristicLength_,kinematicViscosity);
prandtlNumber_[phaseIdx] = DimLessNum::prandtlNumber(dynamicViscosity, heatCapacity, thermalConductivity);
nusseltNumber_[phaseIdx] = DimLessNum::nusseltNumberForced(reynoldsNumber_[phaseIdx],
prandtlNumber_[phaseIdx],
@@ -145,6 +151,80 @@ public:
}
}
+ /*!
+ * \brief Updates the volume specific interfacial area [m^2 / m^3] between the phases.
+ *
+ * \param elemSol A vector containing all primary variables connected to the element
+ * \param fluidState Container for all the secondary variables concerning the fluids
+ * \param paramCache The parameter cache corresponding to the fluid state
+ * \param problem The problem to be solved
+ * \param element An element which contains part of the control volume
+ * \param scv The sub-control volume
+ */
+ template<class ElemSol, class Problem, class Element, class Scv>
+ void updateInterfacialArea(const ElemSol& elemSol,
+ const FluidState& fluidState,
+ const ParameterCache& paramCache,
+ const Problem& problem,
+ const Element& element,
+ const Scv& scv)
+ {
+ // obtain (standard) material parameters (needed for the residual saturations)
+ const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
+
+ //obtain parameters for interfacial area constitutive relations
+ const auto& aWettingNonWettingSurfaceParams =problem.spatialParams().aWettingNonWettingSurfaceParams(element, scv, elemSol);
+ const auto& aNonWettingSolidSurfaceParams =problem.spatialParams().aNonWettingSolidSurfaceParams(element, scv, elemSol);
+
+ const Scalar pc = fluidState.pressure(phase1Idx) - fluidState.pressure(phase0Idx);
+ const Scalar Sw = fluidState.saturation(phase0Idx);
+
+ Scalar awn;
+
+ using AwnSurface = typename Problem::SpatialParams::AwnSurface;
+ awn = AwnSurface::interfacialArea(aWettingNonWettingSurfaceParams, materialParams, Sw, pc );
+ interfacialArea_[phase0Idx][phase1Idx] = awn;
+ interfacialArea_[phase1Idx][phase0Idx] = interfacialArea_[phase0Idx][phase1Idx];
+ interfacialArea_[phase0Idx][phase0Idx] = 0.;
+
+ using AnsSurface = typename Problem::SpatialParams::AnsSurface;
+ Scalar ans = AnsSurface::interfacialArea(aNonWettingSolidSurfaceParams, materialParams,Sw, pc);
+
+ // Switch for using a a_{wn} relations that has some "maximum capillary pressure" as parameter // That value is obtained by regularization of the pc(Sw) function.
+#if USE_PCMAX
+ const Scalar pcMax = problem.spatialParams().pcMax(element, scv, elemSol);
+ //I know the solid surface from the pore network. But it is more consistent to use the fitvalue.
+ using AnsSurface = typename Problem::SpatialParams::AnsSurface;
+ solidSurface_ = AnsSurface::interfacialArea(aNonWettingSolidSurfaceParams, materialParams, /*Sw=*/0., pcMax);
+
+ const Scalar aws = solidSurface_ - ans;
+ interfacialArea_[phase0Idx][sPhaseIdx] = aws;
+ interfacialArea_[sPhaseIdx][phase0Idx] = interfacialArea_[phase0Idx][sPhaseIdx];
+ interfacialArea_[sPhaseIdx][sPhaseIdx] = 0.;
+#else
+ using AwsSurface = typename Problem::SpatialParams::AwsSurface;
+ const auto& aWettingSolidSurfaceParams = problem.spatialParams().aWettingSolidSurfaceParams(element, scv, elemSol);
+ const auto aws = AwsSurface::interfacialArea(aWettingSolidSurfaceParams,materialParams, Sw, pc);
+ interfacialArea_[phase0Idx][sPhaseIdx] = aws ;
+ interfacialArea_[sPhaseIdx][phase0Idx] = interfacialArea_[phase0Idx][sPhaseIdx];
+ interfacialArea_[sPhaseIdx][sPhaseIdx] = 0.;
+#endif
+ interfacialArea_[phase1Idx][sPhaseIdx] = ans;
+ interfacialArea_[sPhaseIdx][phase1Idx] = interfacialArea_[phase1Idx][sPhaseIdx];
+ interfacialArea_[phase1Idx][phase1Idx] = 0.;
+ }
+
+ /*!
+ * \brief The specific interfacial area between two fluid phases [m^2 / m^3]
+ * \note This is _only_ required by the kinetic mass/energy modules
+ */
+ const Scalar interfacialArea(const unsigned int phaseIIdx, const unsigned int phaseJIdx) const
+ {
+ // there is no interfacial area between a phase and itself
+ assert(phaseIIdx not_eq phaseJIdx);
+ return interfacialArea_[phaseIIdx][phaseJIdx];
+ }
+
//! access function Reynolds Number
const Scalar reynoldsNumber(const unsigned int phaseIdx) const { return reynoldsNumber_[phaseIdx]; }
//! access function Prandtl Number
@@ -168,6 +248,7 @@ private:
Scalar factorEnergyTransfer_;
Scalar factorMassTransfer_;
Scalar solidSurface_ ;
+ Scalar interfacialArea_[ModelTraits::numPhases()+numEnergyEqSolid][ModelTraits::numPhases()+numEnergyEqSolid];
};
////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -390,24 +471,6 @@ public:
Scalar awn;
- // TODO can we delete this??? awn is overwritten anyway!
-#define AwnRegul 0
- // This regularizes the interfacial area between the fluid phases.
- // This makes sure, that
- // a) some saturation cannot be lost: Never leave two phase region.
- // b) We cannot leave the fit region: no crazy (e.g. negative) values possible
-
- // const Scalar Swr = aWettingNonWettingSurfaceParams.Swr() ;
- // const Scalar Snr = aWettingNonWettingSurfaceParams.Snr() ;
-
- // this just leads to a stalling newton error as soon as this kicks in.
- // May be a spline or sth like this would help, but I do not which derivatives
- // to specify.
-#if AwnRegul
- if(Sw < 5e-3 ) // or Sw > (1.-1e-5 )
- awn = 0. ; // 10.; //
-#endif
-
using AwnSurface = typename Problem::SpatialParams::AwnSurface;
awn = AwnSurface::interfacialArea(aWettingNonWettingSurfaceParams, materialParams, Sw, pc );
interfacialArea_[phase0Idx][phase1Idx] = awn;
diff --git a/test/porousmediumflow/mpnc/implicit/kinetic/spatialparams.hh b/test/porousmediumflow/mpnc/implicit/kinetic/spatialparams.hh
index 58f32f356b..b5088fa2d4 100644
--- a/test/porousmediumflow/mpnc/implicit/kinetic/spatialparams.hh
+++ b/test/porousmediumflow/mpnc/implicit/kinetic/spatialparams.hh
@@ -26,7 +26,7 @@
#define DUMUX_EVAPORATION_ATMOSPHERE_SPATIALPARAMS_HH
#include <dumux/porousmediumflow/properties.hh>
-#include <dumux/material/spatialparams/fv.hh>
+#include <dumux/material/spatialparams/fvnonequilibrium.hh>
#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
@@ -36,6 +36,12 @@
#include <dumux/material/fluidmatrixinteractions/mp/2poftadapter.hh>
#include <dumux/common/parameters.hh>
+// material laws for interfacial area
+#include <dumux/material/fluidmatrixinteractions/2pia/efftoabslawia.hh>
+#include <dumux/material/fluidmatrixinteractions/2pia/awnsurfacepolynomial2ndorder.hh>
+#include <dumux/material/fluidmatrixinteractions/2pia/awnsurfacepcmaxfct.hh>
+#include <dumux/material/fluidmatrixinteractions/2pia/awnsurfaceexpswpcto3.hh>
+
namespace Dumux {
/**
@@ -43,9 +49,9 @@ namespace Dumux {
*/
template<class TypeTag>
class EvaporationAtmosphereSpatialParams
-: public FVSpatialParams<GetPropType<TypeTag, Properties::FVGridGeometry>,
- GetPropType<TypeTag, Properties::Scalar>,
- EvaporationAtmosphereSpatialParams<TypeTag>>
+: public FVNonEquilibriumSpatialParams<GetPropType<TypeTag, Properties::FVGridGeometry>,
+ GetPropType<TypeTag, Properties::Scalar>,
+ EvaporationAtmosphereSpatialParams<TypeTag>>
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
@@ -53,7 +59,7 @@ class EvaporationAtmosphereSpatialParams
using FVElementGeometry = typename FVGridGeometry::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using Element = typename GridView::template Codim<0>::Entity;
- using ParentType = FVSpatialParams<FVGridGeometry, Scalar, EvaporationAtmosphereSpatialParams<TypeTag>>;
+ using ParentType = FVNonEquilibriumSpatialParams<FVGridGeometry, Scalar, EvaporationAtmosphereSpatialParams<TypeTag>>;
using GlobalPosition = Dune::FieldVector<Scalar, GridView::dimension>;
@@ -69,13 +75,17 @@ public:
using PermeabilityType = Scalar;
//! export the material law type used
using MaterialLaw = TwoPAdapter<liquidPhaseIdx, EffToAbsLaw<RegularizedBrooksCorey<Scalar>>>;
- //! export the types used for interfacial area calculations
- using AwnSurface = GetPropType<TypeTag, Properties::AwnSurface>;
- using AwsSurface = GetPropType<TypeTag, Properties::AwsSurface>;
- using AnsSurface = GetPropType<TypeTag, Properties::AnsSurface>;
-
//! convenience aliases of the law parameters
using MaterialLawParams = typename MaterialLaw::Params;
+
+ //! export the types used for interfacial area calculations
+ using EffectiveIALawAws = AwnSurfacePolynomial2ndOrder<Scalar>;
+ using EffectiveIALawAwn = AwnSurfacePcMaxFct<Scalar>;
+ using EffectiveIALawAns = AwnSurfaceExpSwPcTo3<Scalar>;
+ using AwnSurface = EffToAbsLawIA<EffectiveIALawAwn, MaterialLawParams>;
+ using AwsSurface = EffToAbsLawIA<EffectiveIALawAws, MaterialLawParams>;
+ using AnsSurface = EffToAbsLawIA<EffectiveIALawAns, MaterialLawParams>;
+
using AwnSurfaceParams = typename AwnSurface::Params;
using AwsSurfaceParams = typename AwsSurface::Params;
using AnsSurfaceParams = typename AnsSurface::Params;
@@ -284,21 +294,6 @@ public:
const ElementSolution &elemSol) const
{ return aWettingNonWettingSurfaceParams_.pcMax() ; }
-
- /*!\brief Return the characteristic length for the mass transfer.
- *
- * The position is determined based on the coordinate of
- * the vertex belonging to the considered sub controle volume.
- * \param element The finite element
- * \param fvGeometry The finite volume geometry
- * \param scvIdx The local index of the sub control volume */
- template<class ElementSolution>
- const Scalar characteristicLength(const Element & element,
- const SubControlVolume &scv,
- const ElementSolution &elemSol) const
-
- { return characteristicLengthAtPos(scv.dofPosition()); }
-
/*!\brief Return the characteristic length for the mass transfer.
* \param globalPos The position in global coordinates.*/
const Scalar characteristicLengthAtPos(const GlobalPosition & globalPos) const
@@ -310,19 +305,6 @@ public:
else DUNE_THROW(Dune::InvalidStateException, "You should not be here: x=" << globalPos[0] << " y= "<< globalPos[dimWorld-1]);
}
- /*!\brief Return the pre factor the the energy transfer
- *
- * The position is determined based on the coordinate of
- * the vertex belonging to the considered sub controle volume.
- * \param element The finite element
- * \param fvGeometry The finite volume geometry
- * \param scvIdx The local index of the sub control volume */
- template<class ElementSolution>
- const Scalar factorEnergyTransfer(const Element &element,
- const SubControlVolume &scv,
- const ElementSolution &elemSol) const
- { return factorEnergyTransferAtPos(scv.dofPosition()); }
-
/*!\brief Return the pre factor the the energy transfer
* \param globalPos The position in global coordinates.*/
const Scalar factorEnergyTransferAtPos(const GlobalPosition & globalPos) const
@@ -334,19 +316,6 @@ public:
else DUNE_THROW(Dune::InvalidStateException, "You should not be here: x=" << globalPos[0] << " y= "<< globalPos[dimWorld-1]);
}
- /*!\brief Return the pre factor the the mass transfer
- *
- * The position is determined based on the coordinate of
- * the vertex belonging to the considered sub controle volume.
- * \param element The finite element
- * \param fvGeometry The finite volume geometry
- * \param scvIdx The local index of the sub control volume */
- template<class ElementSolution>
- const Scalar factorMassTransfer(const Element &element,
- const SubControlVolume &scv,
- const ElementSolution &elemSol) const
- { return factorMassTransferAtPos(scv.dofPosition()); }
-
/*!\brief Return the pre factor the the mass transfer
* \param globalPos The position in global coordinates.*/
const Scalar factorMassTransferAtPos(const GlobalPosition & globalPos) const
diff --git a/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/spatialparams.hh b/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/spatialparams.hh
index edb715a669..8b80f1faaa 100644
--- a/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/spatialparams.hh
+++ b/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/spatialparams.hh
@@ -26,6 +26,7 @@
#include <dune/common/parametertreeparser.hh>
+#include <dumux/material/spatialparams/fvnonequilibrium.hh>
#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
#include <dumux/material/fluidmatrixinteractions/2p/heatpipelaw.hh>
#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
@@ -33,6 +34,12 @@
#include <dumux/porousmediumflow/properties.hh>
#include <dumux/material/spatialparams/fv.hh>
+// material laws for interfacial area
+#include <dumux/material/fluidmatrixinteractions/2pia/efftoabslawia.hh>
+#include <dumux/material/fluidmatrixinteractions/2pia/awnsurfacepolynomial2ndorder.hh>
+#include <dumux/material/fluidmatrixinteractions/2pia/awnsurfacepcmaxfct.hh>
+#include <dumux/material/fluidmatrixinteractions/2pia/awnsurfaceexpswpcto3.hh>
+
namespace Dumux {
/**
@@ -41,9 +48,9 @@ namespace Dumux {
*/
template<class TypeTag>
class CombustionSpatialParams
-: public FVSpatialParams<GetPropType<TypeTag, Properties::FVGridGeometry>,
- GetPropType<TypeTag, Properties::Scalar>,
- CombustionSpatialParams<TypeTag>>
+: public FVNonEquilibriumSpatialParams<GetPropType<TypeTag, Properties::FVGridGeometry>,
+ GetPropType<TypeTag, Properties::Scalar>,
+ CombustionSpatialParams<TypeTag>>
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
@@ -51,7 +58,7 @@ class CombustionSpatialParams
using FVElementGeometry = typename FVGridGeometry::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using Element = typename GridView::template Codim<0>::Entity;
- using ParentType = FVSpatialParams<FVGridGeometry, Scalar, CombustionSpatialParams<TypeTag>>;
+ using ParentType = FVNonEquilibriumSpatialParams<FVGridGeometry, Scalar, CombustionSpatialParams<TypeTag>>;
enum {dimWorld = GridView::dimensionworld};
using GlobalPosition = typename SubControlVolume::GlobalPosition;
@@ -68,6 +75,18 @@ public:
using MaterialLaw = TwoPAdapter<wPhaseIdx, EffToAbsLaw<EffectiveLaw>>;
using MaterialLawParams = typename MaterialLaw::Params;
+ //! export the types used for interfacial area calculations
+ using EffectiveIALawAws = AwnSurfacePolynomial2ndOrder<Scalar>;
+ using EffectiveIALawAwn = AwnSurfacePcMaxFct<Scalar>;
+ using EffectiveIALawAns = AwnSurfaceExpSwPcTo3<Scalar>;
+ using AwnSurface = EffToAbsLawIA<EffectiveIALawAwn, MaterialLawParams>;
+ using AwsSurface = EffToAbsLawIA<EffectiveIALawAws, MaterialLawParams>;
+ using AnsSurface = EffToAbsLawIA<EffectiveIALawAns, MaterialLawParams>;
+
+ using AwnSurfaceParams = typename AwnSurface::Params;
+ using AwsSurfaceParams = typename AwsSurface::Params;
+ using AnsSurfaceParams = typename AnsSurface::Params;
+
CombustionSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry) : ParentType(fvGridGeometry)
{
// this is the parameter value from file part
@@ -166,59 +185,17 @@ public:
const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition & globalPos) const
{ return materialParams_ ; }
- /*!\brief Return the characteristic length for the mass transfer.
- *
- * The position is determined based on the coordinate of
- * the vertex belonging to the considered sub controle volume.
- * \param element The finite element
- * \param fvGeometry The finite volume geometry
- * \param scvIdx The local index of the sub control volume */
- template<class ElementSolution>
- const Scalar characteristicLength(const Element & element,
- const SubControlVolume &scv,
- const ElementSolution &elemSol) const
-
- { return characteristicLengthAtPos(scv.center()); }
-
-
/*!\brief Return the characteristic length for the mass transfer.
* \param globalPos The position in global coordinates.*/
const Scalar characteristicLengthAtPos(const GlobalPosition & globalPos) const
{ return characteristicLength_ ; }
- /*!\brief Return the pre factor the the energy transfer
- *
- * The position is determined based on the coordinate of
- * the vertex belonging to the considered sub controle volume.
- * \param element The finite element
- * \param fvGeometry The finite volume geometry
- * \param scvIdx The local index of the sub control volume */
- template<class ElementSolution>
- const Scalar factorEnergyTransfer(const Element &element,
- const SubControlVolume &scv,
- const ElementSolution &elemSol) const
- { return factorEnergyTransferAtPos(scv.dofPosition()); }
-
-
/*!\brief Return the pre factor the the energy transfer
* \param globalPos The position in global coordinates.*/
const Scalar factorEnergyTransferAtPos(const GlobalPosition & globalPos) const
{ return factorEnergyTransfer_; }
- /*!\brief Return the pre factor the the mass transfer
- *
- * The position is determined based on the coordinate of
- * the vertex belonging to the considered sub controle volume.
- * \param element The finite element
- * \param fvGeometry The finite volume geometry
- * \param scvIdx The local index of the sub control volume */
- template<class ElementSolution>
- const Scalar factorMassTransfer(const Element &element,
- const SubControlVolume &scv,
- const ElementSolution &elemSol) const
- { return factorMassTransferAtPos(scv.dofPosition()); }
-
/*!\brief Return the pre factor the the mass transfer
* \param globalPos The position in global coordinates.*/
const Scalar factorMassTransferAtPos(const GlobalPosition & globalPos) const
--
GitLab