diff --git a/dumux/discretization/cellcentered/tpfa/fourierslawnonequilibrium.hh b/dumux/discretization/cellcentered/tpfa/fourierslawnonequilibrium.hh
new file mode 100644
index 0000000000000000000000000000000000000000..f4ace0d089b05c1f4c6548354a84929ba92efec5
--- /dev/null
+++ b/dumux/discretization/cellcentered/tpfa/fourierslawnonequilibrium.hh
@@ -0,0 +1,190 @@
+// -*- 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 CCTpfaDiscretization
+* \brief Fourier's law for cell-centered finite volume schemes with two-point flux approximation
+*/
+#ifndef DUMUX_DISCRETIZATION_CC_TPFA_FOURIERS_LAW_NONEQUILIBRIUM_HH
+#define DUMUX_DISCRETIZATION_CC_TPFA_FOURIERS_LAW_NONEQUILIBRIUM_HH
+
+#include <dumux/common/parameters.hh>
+#include <dumux/common/properties.hh>
+
+#include <dumux/discretization/methods.hh>
+#include <dumux/discretization/cellcentered/tpfa/computetransmissibility.hh>
+
+namespace Dumux
+{
+// forward declaration
+template<class TypeTag, DiscretizationMethod discMethod>
+class FouriersLawNonEquilibriumImplementation;
+
+/*!
+* \ingroup CCTpfaDiscretization
+* \brief Fourier's law for cell-centered finite volume schemes with two-point flux approximation
+*/
+template <class TypeTag>
+class FouriersLawNonEquilibriumImplementation<TypeTag, DiscretizationMethod::cctpfa>
+{
+ using Implementation = FouriersLawNonEquilibriumImplementation<TypeTag, DiscretizationMethod::cctpfa>;
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using IndexType = typename GridView::IndexSet::IndexType;
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using ElementFluxVarsCache = typename GET_PROP_TYPE(TypeTag, GridFluxVariablesCache)::LocalView;
+ using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
+
+ static const int dim = GridView::dimension;
+ static const int dimWorld = GridView::dimensionworld;
+
+ using DimWorldMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
+
+ using ThermalConductivityModel = typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel);
+
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ enum { numEnergyEqFluid = GET_PROP_VALUE(TypeTag, NumEnergyEqFluid) };
+ enum {sPhaseIdx = FluidSystem::numPhases};
+
+public:
+ //! state the discretization method this implementation belongs to
+ static const DiscretizationMethod discMethod = DiscretizationMethod::cctpfa;
+
+ using Cache = FluxVariablesCaching::EmptyHeatConductionCache;
+
+ //! Compute the heat condution flux assuming thermal equilibrium
+ static Scalar flux(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const int phaseIdx,
+ const ElementFluxVarsCache& elemFluxVarsCache)
+ {
+ Scalar tInside = 0.0;
+ Scalar tOutside = 0.0;
+ // get the inside/outside temperatures
+ if (phaseIdx < numEnergyEqFluid)
+ {
+ tInside += elemVolVars[scvf.insideScvIdx()].temperatureFluid(phaseIdx);
+ tOutside += elemVolVars[scvf.outsideScvIdx()].temperatureFluid(phaseIdx);
+ }
+ else //temp solid
+ {
+ tInside += elemVolVars[scvf.insideScvIdx()].temperatureSolid();
+ tOutside += elemVolVars[scvf.outsideScvIdx()].temperatureSolid();
+ }
+
+ Scalar tij = calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx);
+ return tij*(tInside - tOutside);
+ }
+
+ //! Compute transmissibilities
+ static Scalar calculateTransmissibility(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const int phaseIdx)
+ {
+ Scalar tij;
+
+ const auto insideScvIdx = scvf.insideScvIdx();
+ const auto& insideScv = fvGeometry.scv(insideScvIdx);
+ const auto& insideVolVars = elemVolVars[insideScvIdx];
+ Scalar insideLambda = 0.0;
+ Scalar outsideLambda = 0.0;
+
+ // effective diffusion tensors
+ if (phaseIdx != sPhaseIdx)
+ {
+ if (numEnergyEqFluid == 1)
+ { //when only one energy equation for fluids is used, we need an effective law for that
+ insideLambda += ThermalConductivityModel::effectiveThermalConductivity(insideVolVars, problem.spatialParams(), element, fvGeometry, insideScv);
+ }
+ else //numEnergyEqFluid >1
+ {
+ insideLambda += insideVolVars.fluidThermalConductivity(phaseIdx)*insideVolVars.saturation(phaseIdx)*insideVolVars.porosity();
+ }
+ }
+ //solid phase
+ else
+ {
+ insideLambda += insideVolVars.solidThermalConductivity()*(1-insideVolVars.porosity());
+ }
+
+ const Scalar ti = computeTpfaTransmissibility(scvf, insideScv, insideLambda, insideVolVars.extrusionFactor());
+
+ // for the boundary (dirichlet) or at branching points we only need ti
+ if (scvf.boundary() || scvf.numOutsideScvs() > 1)
+ {
+ tij = scvf.area()*ti;
+ }
+ // otherwise we compute a tpfa harmonic mean
+ else
+ {
+ const auto outsideScvIdx = scvf.outsideScvIdx();
+ const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
+ const auto& outsideVolVars = elemVolVars[outsideScvIdx];
+
+ // effective diffusion tensors
+ if (phaseIdx != sPhaseIdx)
+ {
+ if (numEnergyEqFluid == 1)
+ { //when only one energy equation for fluids is used, we need an effective law for that
+ outsideLambda += ThermalConductivityModel::effectiveThermalConductivity(outsideVolVars, problem.spatialParams(), element, fvGeometry, outsideScv);
+ }
+ else
+ {
+ outsideLambda += outsideVolVars.fluidThermalConductivity(phaseIdx)*outsideVolVars.saturation(phaseIdx)*outsideVolVars.porosity();
+ }
+ }
+ //solid phase
+ else
+ {
+ outsideLambda +=outsideVolVars.solidThermalConductivity()*(1-outsideVolVars.porosity());
+ }
+ Scalar tj;
+ if (dim == dimWorld)
+ // assume the normal vector from outside is anti parallel so we save flipping a vector
+ tj = -1.0*computeTpfaTransmissibility(scvf, outsideScv, outsideLambda, outsideVolVars.extrusionFactor());
+ else
+ tj = computeTpfaTransmissibility(fvGeometry.flipScvf(scvf.index()), outsideScv, outsideLambda, outsideVolVars.extrusionFactor());
+
+ // check for division by zero!
+ if (ti*tj <= 0.0)
+ tij = 0;
+ else
+ tij = scvf.area()*(ti * tj)/(ti + tj);
+ }
+ return tij;
+ }
+
+private:
+
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/flux/fourierslawnonequilibrium.hh b/dumux/flux/fourierslawnonequilibrium.hh
index cbd9c38946d9d29decd08e89b7f6cc2795383d20..ad6822283a5826727ca3837d65dd69c2d2a1f91f 100644
--- a/dumux/flux/fourierslawnonequilibrium.hh
+++ b/dumux/flux/fourierslawnonequilibrium.hh
@@ -43,6 +43,7 @@ using FouriersLawNonEquilibrium = FouriersLawNonEquilibriumImplementation<TypeTa
} // end namespace Dumux
+#include <dumux/discretization/cellcentered/tpfa/fourierslawnonequilibrium.hh>
#include <dumux/discretization/box/fourierslawnonequilibrium.hh>
#endif
diff --git a/dumux/material/spatialparams/fvnonequilibrium.hh b/dumux/material/spatialparams/fvnonequilibrium.hh
index 98a0eef7f58d6cee494ee9603b254ad95b89f5e1..73ed50d919d37f82e38692b2b902ed7129fabf15 100644
--- a/dumux/material/spatialparams/fvnonequilibrium.hh
+++ b/dumux/material/spatialparams/fvnonequilibrium.hh
@@ -27,12 +27,6 @@
#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 {
/**
@@ -64,78 +58,6 @@ public:
: 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
diff --git a/dumux/porousmediumflow/nonequilibrium/volumevariables.hh b/dumux/porousmediumflow/nonequilibrium/volumevariables.hh
index af5d888243edcc574da785cc9bd3803d506825ca..4ca1cc3b8d7faaeb46540b808ce5cc15278a722c 100644
--- a/dumux/porousmediumflow/nonequilibrium/volumevariables.hh
+++ b/dumux/porousmediumflow/nonequilibrium/volumevariables.hh
@@ -42,7 +42,7 @@ namespace Dumux {
* modules which require the specific interfacial area between
* fluid phases.
*/
-template<class Traits, class EquilibriumVolumeVariables, bool enableChemicalNonEquilibrium ,bool enableThermalNonEquilibrium>
+template<class Traits, class EquilibriumVolumeVariables, bool enableChemicalNonEquilibrium ,bool enableThermalNonEquilibrium, int numEnergyEqFluid>
class NonEquilibriumVolumeVariablesImplementation;
template<class Traits, class EquilibriumVolumeVariables>
@@ -50,16 +50,18 @@ using NonEquilibriumVolumeVariables =
NonEquilibriumVolumeVariablesImplementation< Traits,
EquilibriumVolumeVariables,
Traits::ModelTraits::enableChemicalNonEquilibrium(),
- Traits::ModelTraits::enableThermalNonEquilibrium() >;
+ Traits::ModelTraits::enableThermalNonEquilibrium(),
+ Traits::ModelTraits::numEnergyEqFluid()>;
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
-// specialization for the case of NO kinetic mass but kinetic energy transfer of a fluid mixture and solid
+// specialization for the case of NO kinetic mass but kinetic energy transfer of two fluids and a solid
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
template<class Traits, class EquilibriumVolumeVariables>
class NonEquilibriumVolumeVariablesImplementation< Traits,
EquilibriumVolumeVariables,
false/*chemicalNonEquilibrium?*/,
- true/*thermalNonEquilibrium?*/ >
+ true/*thermalNonEquilibrium?*/,
+ 2>
:public EquilibriumVolumeVariables
{
using ParentType = EquilibriumVolumeVariables;
@@ -67,7 +69,7 @@ class NonEquilibriumVolumeVariablesImplementation< Traits,
using Scalar = typename Traits::PrimaryVariables::value_type;
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();
@@ -75,12 +77,12 @@ class NonEquilibriumVolumeVariablesImplementation< Traits,
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>;
public:
using FluidState = typename Traits::FluidState;
+ using Indices = typename ModelTraits::Indices;
/*!
* \brief Update all quantities for a given control volume
*
@@ -102,12 +104,11 @@ public:
ParameterCache paramCache;
paramCache.updateAll(this->fluidState());
updateDimLessNumbers(elemSol, this->fluidState(), paramCache, problem, element, scv);
- if (numEnergyEqFluid == 2)
- updateInterfacialArea(elemSol, this->fluidState(), paramCache, problem, element, scv);
+ updateInterfacialArea(elemSol, this->fluidState(), paramCache, problem, element, scv);
}
/*!
- * \brief Updates the volume specific interfacial area [m^2 / m^3] between the phases.
+ * \brief Updates dimensionless numbers
*
* \param elemSol A vector containing all primary variables connected to the element
* \param fluidState Container for all the secondary variables concerning the fluids
@@ -148,10 +149,11 @@ public:
prandtlNumber_[phaseIdx],
porosity,
ModelTraits::nusseltFormulation());
+
}
}
- /*!
+ /*!
* \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
@@ -251,6 +253,131 @@ private:
Scalar interfacialArea_[ModelTraits::numPhases()+numEnergyEqSolid][ModelTraits::numPhases()+numEnergyEqSolid];
};
+/////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// specialization for the case of NO kinetic mass but kinetic energy transfer of a fluid mixture and a solid
+/////////////////////////////////////////////////////////////////////////////////////////////////////////////
+template<class Traits, class EquilibriumVolumeVariables>
+class NonEquilibriumVolumeVariablesImplementation< Traits,
+ EquilibriumVolumeVariables,
+ false/*chemicalNonEquilibrium?*/,
+ true/*thermalNonEquilibrium?*/,
+ 1>
+ :public EquilibriumVolumeVariables
+{
+ using ParentType = EquilibriumVolumeVariables;
+ using ParameterCache = typename Traits::FluidSystem::ParameterCache;
+ using Scalar = typename Traits::PrimaryVariables::value_type;
+
+ 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 constexpr auto phase0Idx = FS::phase0Idx;
+ static constexpr auto phase1Idx = FS::phase1Idx;
+ static constexpr auto sPhaseIdx = FS::numPhases;
+
+ using DimLessNum = DimensionlessNumbers<Scalar>;
+
+public:
+ using FluidState = typename Traits::FluidState;
+ /*!
+ * \brief Update all quantities for a given control volume
+ *
+ * \param elemSol A vector containing all primary variables connected to the element
+ * \param problem The object specifying the problem which ought to
+ * be simulated
+ * \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 update(const ElemSol &elemSol,
+ const Problem &problem,
+ const Element &element,
+ const Scv& scv)
+ {
+ // Update parent type (also completes the fluid state)
+ ParentType::update(elemSol, problem, element, scv);
+
+ ParameterCache paramCache;
+ paramCache.updateAll(this->fluidState());
+ //only update of DimLessNumbers is necessary here, as interfacial area is easy due to only one fluid with a solid and is directly computed in localresidual
+ updateDimLessNumbers(elemSol, this->fluidState(), paramCache, problem, element, scv);
+ }
+
+ /*!
+ * \brief Updates dimensionless numbers
+ *
+ * \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 updateDimLessNumbers(const ElemSol& elemSol,
+ const FluidState& fluidState,
+ const ParameterCache& paramCache,
+ const Problem& problem,
+ const Element& element,
+ const Scv& scv)
+ {
+ factorMassTransfer_ = problem.spatialParams().factorMassTransfer(element, scv, elemSol);
+ factorEnergyTransfer_ = problem.spatialParams().factorEnergyTransfer(element, scv, elemSol);
+ characteristicLength_ = problem.spatialParams().characteristicLength(element, scv, elemSol);
+
+ // set the dimensionless numbers and obtain respective quantities
+ const unsigned int vIdxGlobal = scv.dofIndex();
+ for (int phaseIdx = 0; phaseIdx < ModelTraits::numPhases(); ++phaseIdx)
+ {
+ const auto darcyMagVelocity = problem.gridVariables().volumeDarcyMagVelocity(phaseIdx, vIdxGlobal);
+ const auto dynamicViscosity = fluidState.viscosity(phaseIdx);
+ const auto density = fluidState.density(phaseIdx);
+ const auto kinematicViscosity = dynamicViscosity/density;
+
+ using FluidSystem = typename Traits::FluidSystem;
+ const auto heatCapacity = FluidSystem::heatCapacity(fluidState, paramCache, phaseIdx);
+ const auto thermalConductivity = FluidSystem::thermalConductivity(fluidState, paramCache, phaseIdx);
+ const auto porosity = this->porosity();
+
+ reynoldsNumber_[phaseIdx] = DimLessNum::reynoldsNumber(darcyMagVelocity, characteristicLength_,kinematicViscosity);
+ prandtlNumber_[phaseIdx] = DimLessNum::prandtlNumber(dynamicViscosity, heatCapacity, thermalConductivity);
+ nusseltNumber_[phaseIdx] = DimLessNum::nusseltNumberForced(reynoldsNumber_[phaseIdx],
+ prandtlNumber_[phaseIdx],
+ porosity,
+ ModelTraits::nusseltFormulation());
+
+ }
+ }
+
+ //! access function Reynolds Number
+ const Scalar reynoldsNumber(const unsigned int phaseIdx) const { return reynoldsNumber_[phaseIdx]; }
+ //! access function Prandtl Number
+ const Scalar prandtlNumber(const unsigned int phaseIdx) const { return prandtlNumber_[phaseIdx]; }
+ //! access function Nusselt Number
+ const Scalar nusseltNumber(const unsigned int phaseIdx) const { return nusseltNumber_[phaseIdx]; }
+ //! access function characteristic length
+ const Scalar characteristicLength() const { return characteristicLength_; }
+ //! access function pre factor energy transfer
+ const Scalar factorEnergyTransfer() const { return factorEnergyTransfer_; }
+ //! access function pre factor mass transfer
+ const Scalar factorMassTransfer() const { return factorMassTransfer_; }
+
+private:
+ //! dimensionless numbers
+ Scalar reynoldsNumber_[ModelTraits::numPhases()];
+ Scalar prandtlNumber_[ModelTraits::numPhases()];
+ Scalar nusseltNumber_[ModelTraits::numPhases()];
+
+ Scalar characteristicLength_;
+ Scalar factorEnergyTransfer_;
+ Scalar factorMassTransfer_;
+ Scalar solidSurface_ ;
+ Scalar interfacialArea_[ModelTraits::numPhases()+numEnergyEqSolid][ModelTraits::numPhases()+numEnergyEqSolid];
+};
+
////////////////////////////////////////////////////////////////////////////////////////////////////
// specialization for the case of (only) kinetic mass transfer. Be careful, we do not test this!
////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -258,7 +385,8 @@ template<class Traits, class EquilibriumVolumeVariables>
class NonEquilibriumVolumeVariablesImplementation< Traits,
EquilibriumVolumeVariables,
true/*chemicalNonEquilibrium?*/,
- false/*thermalNonEquilibrium?*/>
+ false/*thermalNonEquilibrium?*/,
+ 0>
:public EquilibriumVolumeVariables
{
using ParentType = EquilibriumVolumeVariables;
@@ -297,6 +425,7 @@ public:
ParameterCache paramCache;
paramCache.updateAll(this->fluidState());
+ updateDimLessNumbers(elemSol, this->fluidState(), paramCache, problem, element, scv);
updateInterfacialArea(elemSol, this->fluidState(), paramCache, problem, element, scv);
}
@@ -311,24 +440,13 @@ public:
* \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)
+ void updateDimLessNumbers(const ElemSol& elemSol,
+ const FluidState& fluidState,
+ const ParameterCache& paramCache,
+ const Problem& problem,
+ const Element& element,
+ const Scv& scv)
{
- // obtain parameters for awnsurface and material law
- const auto& awnSurfaceParams = problem.spatialParams().aWettingNonWettingSurfaceParams(element, scv, elemSol) ;
- const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol) ;
-
- const auto Sw = fluidState.saturation(phase0Idx) ;
- const auto pc = fluidState.pressure(phase1Idx) - fluidState.pressure(phase0Idx);
-
- // so far there is only a model for kinetic mass transfer between fluid phases
- using AwnSurface = typename Problem::SpatialParams::AwnSurface;
- interfacialArea_ = AwnSurface::interfacialArea(awnSurfaceParams, materialParams, Sw, pc);
-
factorMassTransfer_ = problem.spatialParams().factorMassTransfer(element, scv, elemSol);
characteristicLength_ = problem.spatialParams().characteristicLength(element, scv, elemSol);
@@ -357,6 +475,36 @@ 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 parameters for awnsurface and material law
+ const auto& awnSurfaceParams = problem.spatialParams().aWettingNonWettingSurfaceParams(element, scv, elemSol) ;
+ const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol) ;
+
+ const auto Sw = fluidState.saturation(phase0Idx) ;
+ const auto pc = fluidState.pressure(phase1Idx) - fluidState.pressure(phase0Idx);
+
+ // so far there is only a model for kinetic mass transfer between fluid phases
+ using AwnSurface = typename Problem::SpatialParams::AwnSurface;
+ interfacialArea_ = AwnSurface::interfacialArea(awnSurfaceParams, materialParams, Sw, pc);
+ }
+
/*!
* \brief The specific interfacial area between two fluid phases [m^2 / m^3]
*/
@@ -394,7 +542,8 @@ template<class Traits, class EquilibriumVolumeVariables>
class NonEquilibriumVolumeVariablesImplementation< Traits,
EquilibriumVolumeVariables,
true/*chemicalNonEquilibrium?*/,
- true/*thermalNonEquilibrium?*/>
+ true/*thermalNonEquilibrium?*/,
+ 2>
:public EquilibriumVolumeVariables
{
using ParentType = EquilibriumVolumeVariables;
@@ -438,6 +587,7 @@ public:
ParameterCache paramCache;
paramCache.updateAll(this->fluidState());
+ updateDimLessNumbers(elemSol, this->fluidState(), paramCache, problem, element, scv);
updateInterfacialArea(elemSol, this->fluidState(), paramCache, problem, element, scv);
}
@@ -452,6 +602,61 @@ public:
* \param scv The sub-control volume
*/
template<class ElemSol, class Problem, class Element, class Scv>
+ void updateDimLessNumbers(const ElemSol& elemSol,
+ const FluidState& fluidState,
+ const ParameterCache& paramCache,
+ const Problem& problem,
+ const Element& element,
+ const Scv& scv)
+ {
+ factorMassTransfer_ = problem.spatialParams().factorMassTransfer(element, scv, elemSol);
+ factorEnergyTransfer_ = problem.spatialParams().factorEnergyTransfer(element, scv, elemSol);
+ characteristicLength_ = problem.spatialParams().characteristicLength(element, scv, elemSol);
+
+ const auto vIdxGlobal = scv.dofIndex();
+ using FluidSystem = typename Traits::FluidSystem;
+ for (int phaseIdx = 0; phaseIdx < ModelTraits::numPhases(); ++phaseIdx)
+ {
+ const auto darcyMagVelocity = problem.gridVariables().volumeDarcyMagVelocity(phaseIdx, vIdxGlobal);
+ const auto dynamicViscosity = fluidState.viscosity(phaseIdx);
+ const auto density = fluidState.density(phaseIdx);
+ const auto kinematicViscosity = dynamicViscosity/density;
+ const auto heatCapacity = FluidSystem::heatCapacity(fluidState, paramCache, phaseIdx);
+ const auto thermalConductivity = FluidSystem::thermalConductivity(fluidState, paramCache, phaseIdx);
+
+ // diffusion coefficient of non-wetting component in wetting phase
+ const auto porosity = this->porosity();
+ const auto diffCoeff = FluidSystem::binaryDiffusionCoefficient(fluidState,
+ paramCache,
+ phaseIdx,
+ wCompIdx,
+ nCompIdx);
+
+ reynoldsNumber_[phaseIdx] = DimLessNum::reynoldsNumber(darcyMagVelocity, characteristicLength_, kinematicViscosity);
+ prandtlNumber_[phaseIdx] = DimLessNum::prandtlNumber(dynamicViscosity, heatCapacity, thermalConductivity);
+ schmidtNumber_[phaseIdx] = DimLessNum::schmidtNumber(dynamicViscosity, density, diffCoeff);
+ nusseltNumber_[phaseIdx] = DimLessNum::nusseltNumberForced(reynoldsNumber_[phaseIdx],
+ prandtlNumber_[phaseIdx],
+ porosity,
+ ModelTraits::nusseltFormulation());
+ // If Diffusion is not enabled, Sherwood is divided by zero
+ sherwoodNumber_[phaseIdx] = DimLessNum::sherwoodNumber(reynoldsNumber_[phaseIdx],
+ schmidtNumber_[phaseIdx],
+ ModelTraits::sherwoodFormulation());
+ }
+ }
+
+ /*!
+ * \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,
@@ -504,42 +709,6 @@ public:
interfacialArea_[phase1Idx][sPhaseIdx] = ans;
interfacialArea_[sPhaseIdx][phase1Idx] = interfacialArea_[phase1Idx][sPhaseIdx];
interfacialArea_[phase1Idx][phase1Idx] = 0.;
-
- factorMassTransfer_ = problem.spatialParams().factorMassTransfer(element, scv, elemSol);
- factorEnergyTransfer_ = problem.spatialParams().factorEnergyTransfer(element, scv, elemSol);
- characteristicLength_ = problem.spatialParams().characteristicLength(element, scv, elemSol);
-
- const auto vIdxGlobal = scv.dofIndex();
- using FluidSystem = typename Traits::FluidSystem;
- for (int phaseIdx = 0; phaseIdx < ModelTraits::numPhases(); ++phaseIdx)
- {
- const auto darcyMagVelocity = problem.gridVariables().volumeDarcyMagVelocity(phaseIdx, vIdxGlobal);
- const auto dynamicViscosity = fluidState.viscosity(phaseIdx);
- const auto density = fluidState.density(phaseIdx);
- const auto kinematicViscosity = dynamicViscosity/density;
- const auto heatCapacity = FluidSystem::heatCapacity(fluidState, paramCache, phaseIdx);
- const auto thermalConductivity = FluidSystem::thermalConductivity(fluidState, paramCache, phaseIdx);
-
- // diffusion coefficient of non-wetting component in wetting phase
- const auto porosity = this->porosity();
- const auto diffCoeff = FluidSystem::binaryDiffusionCoefficient(fluidState,
- paramCache,
- phaseIdx,
- wCompIdx,
- nCompIdx);
-
- reynoldsNumber_[phaseIdx] = DimLessNum::reynoldsNumber(darcyMagVelocity, characteristicLength_, kinematicViscosity);
- prandtlNumber_[phaseIdx] = DimLessNum::prandtlNumber(dynamicViscosity, heatCapacity, thermalConductivity);
- schmidtNumber_[phaseIdx] = DimLessNum::schmidtNumber(dynamicViscosity, density, diffCoeff);
- nusseltNumber_[phaseIdx] = DimLessNum::nusseltNumberForced(reynoldsNumber_[phaseIdx],
- prandtlNumber_[phaseIdx],
- porosity,
- ModelTraits::nusseltFormulation());
- // If Diffusion is not enabled, Sherwood is divided by zero
- sherwoodNumber_[phaseIdx] = DimLessNum::sherwoodNumber(reynoldsNumber_[phaseIdx],
- schmidtNumber_[phaseIdx],
- ModelTraits::sherwoodFormulation());
- }
}
/*!
diff --git a/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/spatialparams.hh b/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/spatialparams.hh
index 8b80f1faaaab6a38085d5e665df76a12e9702a74..7b96235b0f767b6c50134fa280b704c49ebe0a83 100644
--- a/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/spatialparams.hh
+++ b/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/spatialparams.hh
@@ -34,12 +34,6 @@
#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 {
/**
@@ -75,18 +69,6 @@ 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