porousmediumflow/nonequilibrium/volumevariables.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
//
// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
#ifndef DUMUX_NONEQUILIBRIUM_VOLUME_VARIABLES_HH
#define DUMUX_NONEQUILIBRIUM_VOLUME_VARIABLES_HH
 
#include <cassert>
#include <array>
 
#include <dumux/common/dimensionlessnumbers.hh>
#include <dumux/common/parameters.hh>
 
namespace Dumux {
 
template<class Traits, class EquilibriumVolumeVariables, bool enableChemicalNonEquilibrium,
         bool enableThermalNonEquilibrium, int numEnergyEqFluid>
class NonEquilibriumVolumeVariablesImplementation;
 
template<class Traits, class EquilibriumVolumeVariables>
using NonEquilibriumVolumeVariables =
      NonEquilibriumVolumeVariablesImplementation<Traits,
                                                  EquilibriumVolumeVariables,
                                                  Traits::ModelTraits::enableChemicalNonEquilibrium(),
                                                  Traits::ModelTraits::enableThermalNonEquilibrium(),
                                                  Traits::ModelTraits::numEnergyEqFluid()>;
 
// 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?*/,
                                                  2>
: public EquilibriumVolumeVariables
{
    using ParentType = EquilibriumVolumeVariables;
    using ParameterCache = typename Traits::FluidSystem::ParameterCache;
    using Scalar = typename Traits::PrimaryVariables::value_type;
 
    using ModelTraits = typename Traits::ModelTraits;
 
    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>;
 
    using NumFluidPhasesArray = std::array<Scalar, ModelTraits::numFluidPhases()>;
    using InterfacialAreasArray =  std::array<std::array<Scalar, ModelTraits::numFluidPhases()+numEnergyEqSolid>,
                                              ModelTraits::numFluidPhases()+numEnergyEqSolid>;
 
public:
     using FluidState = typename Traits::FluidState;
     using Indices = typename ModelTraits::Indices;
    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());
        updateDimLessNumbers(elemSol, this->fluidState(), paramCache, problem, element, scv);
        updateInterfacialArea(elemSol, this->fluidState(), paramCache, problem, element, scv);
    }
 
    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)
    {
        const auto& spatialParams = problem.spatialParams();
        factorEnergyTransfer_ = spatialParams.factorEnergyTransfer(element, scv, elemSol);
        characteristicLength_ = 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::numFluidPhases(); ++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());
 
        }
    }
 
    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)
    {
        const Scalar pc = fluidState.pressure(phase1Idx) - fluidState.pressure(phase0Idx);
        const Scalar Sw = fluidState.saturation(phase0Idx);
 
        const auto fluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, scv, elemSol);
        const auto areaWN = fluidMatrixInteraction.wettingNonwettingInterface().area(Sw, pc);
        const auto areaNS = fluidMatrixInteraction.nonwettingSolidInterface().area(Sw, pc);
        interfacialArea_[phase0Idx][phase1Idx] = areaWN;
        interfacialArea_[phase1Idx][phase0Idx] = interfacialArea_[phase0Idx][phase1Idx];
        interfacialArea_[phase0Idx][phase0Idx] = 0.0;
 
        // 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.
        static const bool computeAwsFromAnsAndPcMax = getParam<bool>("SpatialParams.ComputeAwsFromAnsAndPcMax", true);
        if (computeAwsFromAnsAndPcMax)
        {
            // I know the solid surface from the pore network. But it is more consistent to use the fit value.
            const Scalar pcMax = fluidMatrixInteraction.wettingNonwettingInterface().basicParams().pcMax();
            const auto solidSurface = fluidMatrixInteraction.nonwettingSolidInterface().area(/*Sw=*/0., pcMax);
            interfacialArea_[phase0Idx][sPhaseIdx] = solidSurface - areaNS;
        }
        else
            interfacialArea_[phase0Idx][sPhaseIdx] = fluidMatrixInteraction.wettingSolidInterface().area(Sw, pc);
 
        interfacialArea_[sPhaseIdx][phase0Idx] = interfacialArea_[phase0Idx][sPhaseIdx];
        interfacialArea_[sPhaseIdx][sPhaseIdx] = 0.0;
 
        interfacialArea_[phase1Idx][sPhaseIdx] = areaNS;
        interfacialArea_[sPhaseIdx][phase1Idx] = interfacialArea_[phase1Idx][sPhaseIdx];
        interfacialArea_[phase1Idx][phase1Idx] = 0.0;
    }
 
    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];
    }
 
    const Scalar reynoldsNumber(const unsigned int phaseIdx) const
    { return reynoldsNumber_[phaseIdx]; }
 
    const Scalar prandtlNumber(const unsigned int phaseIdx) const
    { return prandtlNumber_[phaseIdx]; }
 
    const Scalar nusseltNumber(const unsigned int phaseIdx) const
    { return nusseltNumber_[phaseIdx]; }
 
    const Scalar characteristicLength() const
    { return characteristicLength_; }
 
    const Scalar factorEnergyTransfer() const
    { return factorEnergyTransfer_; }
 
private:
    NumFluidPhasesArray reynoldsNumber_;
    NumFluidPhasesArray prandtlNumber_;
    NumFluidPhasesArray nusseltNumber_;
 
    Scalar characteristicLength_;
    Scalar factorEnergyTransfer_;
    InterfacialAreasArray interfacialArea_;
};
 
// 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 FS = typename Traits::FluidSystem;
    static constexpr auto numEnergyEqFluid = ModelTraits::numEnergyEqFluid();
    static constexpr auto numEnergyEqSolid = ModelTraits::numEnergyEqSolid();
 
    using DimLessNum = DimensionlessNumbers<Scalar>;
 
    using NumFluidPhasesArray = std::array<Scalar, ModelTraits::numFluidPhases()>;
 
public:
    using Indices = typename ModelTraits::Indices;
    using FluidState = typename Traits::FluidState;
 
    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);
        updateInterfacialArea(elemSol, this->fluidState(), paramCache, problem, element, scv);
    }
 
    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)
    {
        const auto& spatialParams = problem.spatialParams();
        factorEnergyTransfer_ = spatialParams.factorEnergyTransfer(element, scv, elemSol);
        characteristicLength_ = 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::numFluidPhases(); ++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());
        }
    }
 
 
    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)
    {
       using FluidSolidInterfacialAreaFormulation = typename Problem::SpatialParams::FluidSolidInterfacialAreaFormulation;
       interfacialArea_ = FluidSolidInterfacialAreaFormulation::fluidSolidInterfacialArea(this->porosity(), characteristicLength());
    }
 
    const Scalar reynoldsNumber(const unsigned int phaseIdx) const
    { return reynoldsNumber_[phaseIdx]; }
 
    const Scalar prandtlNumber(const unsigned int phaseIdx) const
    { return prandtlNumber_[phaseIdx]; }
 
    const Scalar nusseltNumber(const unsigned int phaseIdx) const
    { return nusseltNumber_[phaseIdx]; }
 
    const Scalar characteristicLength() const
    { return characteristicLength_; }
 
    const Scalar factorEnergyTransfer() const
    { return factorEnergyTransfer_; }
 
    const Scalar fluidSolidInterfacialArea() const
    {return interfacialArea_;}
 
private:
    NumFluidPhasesArray reynoldsNumber_;
    NumFluidPhasesArray prandtlNumber_;
    NumFluidPhasesArray nusseltNumber_;
 
    Scalar characteristicLength_;
    Scalar factorEnergyTransfer_;
    Scalar interfacialArea_ ;
};
 
// specialization for the case of (only) kinetic mass transfer. Be careful, we do not test this!
template<class Traits, class EquilibriumVolumeVariables>
class NonEquilibriumVolumeVariablesImplementation<Traits,
                                                  EquilibriumVolumeVariables,
                                                  true/*chemicalNonEquilibrium?*/,
                                                  false/*thermalNonEquilibrium?*/,
                                                  0>
: public EquilibriumVolumeVariables
{
    using ParentType = EquilibriumVolumeVariables;
    using FluidState = typename Traits::FluidState;
    using FS = typename Traits::FluidSystem;
    using ParameterCache = typename Traits::FluidSystem::ParameterCache;
    using Scalar = typename Traits::PrimaryVariables::value_type;
 
    using ModelTraits = typename Traits::ModelTraits;
 
    static constexpr auto phase0Idx = FS::phase0Idx;
    static constexpr auto phase1Idx = FS::phase1Idx;
    static constexpr auto wCompIdx = FS::comp0Idx;
    static constexpr auto nCompIdx = FS::comp1Idx;
 
    using DimLessNum = DimensionlessNumbers<Scalar>;
 
    using NumFluidPhasesArray = std::array<Scalar, ModelTraits::numFluidPhases()>;
 
public:
    using Indices = typename ModelTraits::Indices;
    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());
        updateDimLessNumbers(elemSol, this->fluidState(), paramCache, problem, element, scv);
        updateInterfacialArea(elemSol, this->fluidState(), paramCache, problem, element, scv);
    }
 
    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)
    {
        const auto& spatialParams = problem.spatialParams();
        factorMassTransfer_ = spatialParams.factorMassTransfer(element, scv, elemSol);
        characteristicLength_ = 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::numFluidPhases(); ++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;
 
            // diffusion coefficient of nonwetting component in wetting phase
            using FluidSystem = typename Traits::FluidSystem;
            const auto diffCoeff = FluidSystem::binaryDiffusionCoefficient(fluidState,
                                                                           paramCache,
                                                                           phaseIdx,
                                                                           wCompIdx,
                                                                           nCompIdx);
 
            reynoldsNumber_[phaseIdx] = DimLessNum::reynoldsNumber(darcyMagVelocity, characteristicLength_, kinematicViscosity);
            schmidtNumber_[phaseIdx] = DimLessNum::schmidtNumber(dynamicViscosity, density, diffCoeff);
            sherwoodNumber_[phaseIdx] = DimLessNum::sherwoodNumber(reynoldsNumber_[phaseIdx],
                                                                   schmidtNumber_[phaseIdx],
                                                                   ModelTraits::sherwoodFormulation());
        }
    }
 
    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)
    {
        const auto Sw = fluidState.saturation(phase0Idx) ;
        const auto pc = fluidState.pressure(phase1Idx) - fluidState.pressure(phase0Idx);
 
        // when we only consider chemical non-equilibrium there is only mass transfer between
        // the fluid phases, so in 2p only interfacial area between wetting and non-wetting
        const auto fluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, scv, elemSol);
        interfacialArea_ = fluidMatrixInteraction.wettingNonwettingInterface().area(Sw, pc);
    }
 
    const Scalar interfacialArea(const unsigned int phaseIIdx, const unsigned int phaseJIdx) const
    {
        // so far there is only a model for kinetic mass transfer between fluid phases
        assert( (phaseIIdx == phase1Idx && phaseJIdx == phase0Idx)
                || (phaseIIdx == phase0Idx && phaseJIdx == phase1Idx) );
        return interfacialArea_;
    }
 
    const Scalar reynoldsNumber(const unsigned int phaseIdx) const
    { return reynoldsNumber_[phaseIdx]; }
 
    const Scalar schmidtNumber(const unsigned int phaseIdx) const
    { return schmidtNumber_[phaseIdx]; }
 
    const Scalar sherwoodNumber(const unsigned int phaseIdx) const
    { return sherwoodNumber_[phaseIdx]; }
 
    const Scalar characteristicLength() const
    { return characteristicLength_; }
 
    const Scalar factorMassTransfer() const
    { return factorMassTransfer_; }
 
private:
    Scalar characteristicLength_;
    Scalar factorMassTransfer_;
    Scalar interfacialArea_ ;
    NumFluidPhasesArray sherwoodNumber_;
    NumFluidPhasesArray schmidtNumber_;
    NumFluidPhasesArray reynoldsNumber_;
};
 
// Specialization where everything is in non-equilibrium.
template<class Traits, class EquilibriumVolumeVariables>
class NonEquilibriumVolumeVariablesImplementation<Traits,
                                                  EquilibriumVolumeVariables,
                                                  true/*chemicalNonEquilibrium?*/,
                                                  true/*thermalNonEquilibrium?*/,
                                                  2>
: public EquilibriumVolumeVariables
{
    using ParentType = EquilibriumVolumeVariables;
    using FluidState = typename Traits::FluidState;
    using FS = typename Traits::FluidSystem;
    using ParameterCache = typename Traits::FluidSystem::ParameterCache;
    using Scalar = typename Traits::PrimaryVariables::value_type;
 
    using ModelTraits = typename Traits::ModelTraits;
    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;
    static constexpr auto wCompIdx = FS::comp0Idx;
    static constexpr auto nCompIdx = FS::comp1Idx;
 
    using DimLessNum = DimensionlessNumbers<Scalar>;
    static_assert((numEnergyEqFluid > 1), "This model only deals with energy transfer between two fluids and one solid phase");
 
    using NumFluidPhasesArray = std::array<Scalar, ModelTraits::numFluidPhases()>;
    using InterfacialAreasArray =  std::array<std::array<Scalar, ModelTraits::numFluidPhases()+numEnergyEqSolid>,
                                              ModelTraits::numFluidPhases()+numEnergyEqSolid>;
 
public:
    using Indices = typename ModelTraits::Indices;
 
    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());
        updateDimLessNumbers(elemSol, this->fluidState(), paramCache, problem, element, scv);
        updateInterfacialArea(elemSol, this->fluidState(), paramCache, problem, element, scv);
    }
 
    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)
    {
        const auto& spatialParams = problem.spatialParams();
        factorMassTransfer_ = spatialParams.factorMassTransfer(element, scv, elemSol);
        factorEnergyTransfer_ = spatialParams.factorEnergyTransfer(element, scv, elemSol);
        characteristicLength_ = spatialParams.characteristicLength(element, scv, elemSol);
 
        const auto vIdxGlobal = scv.dofIndex();
        using FluidSystem = typename Traits::FluidSystem;
        for (int phaseIdx = 0; phaseIdx < ModelTraits::numFluidPhases(); ++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 nonwetting 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());
        }
    }
 
    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)
    {
        const Scalar pc = fluidState.pressure(phase1Idx) - fluidState.pressure(phase0Idx);
        const Scalar Sw = fluidState.saturation(phase0Idx);
 
        const auto fluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, scv, elemSol);
 
        const auto awn = fluidMatrixInteraction.wettingNonwettingInterface().area(Sw, pc);
        interfacialArea_[phase0Idx][phase1Idx] = awn;
        interfacialArea_[phase1Idx][phase0Idx] = interfacialArea_[phase0Idx][phase1Idx];
        interfacialArea_[phase0Idx][phase0Idx] = 0.;
 
        const auto areaNS = fluidMatrixInteraction.nonwettingSolidInterface().area(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.
        static const bool computeAwsFromAnsAndPcMax = getParam<bool>("SpatialParams.ComputeAwsFromAnsAndPcMax", true);
        if (computeAwsFromAnsAndPcMax)
        {
            // I know the solid surface from the pore network. But it is more consistent to use the fit value.
            const Scalar pcMax = fluidMatrixInteraction.wettingNonwettingInterface().basicParams().pcMax();
            const auto solidSurface = fluidMatrixInteraction.nonwettingSolidInterface().area(/*Sw=*/0., pcMax);
            interfacialArea_[phase0Idx][sPhaseIdx] = solidSurface - areaNS;
        }
        else
            interfacialArea_[phase0Idx][sPhaseIdx] = fluidMatrixInteraction.wettingSolidInterface().area(Sw, pc);
 
        interfacialArea_[sPhaseIdx][phase0Idx] = interfacialArea_[phase0Idx][sPhaseIdx];
        interfacialArea_[sPhaseIdx][sPhaseIdx] = 0.;
 
        interfacialArea_[phase1Idx][sPhaseIdx] = areaNS;
        interfacialArea_[sPhaseIdx][phase1Idx] = interfacialArea_[phase1Idx][sPhaseIdx];
        interfacialArea_[phase1Idx][phase1Idx] = 0.;
    }
 
    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];
    }
 
    const Scalar reynoldsNumber(const unsigned int phaseIdx) const
    { return reynoldsNumber_[phaseIdx]; }
 
    const Scalar prandtlNumber(const unsigned int phaseIdx) const
    { return prandtlNumber_[phaseIdx]; }
 
    const Scalar nusseltNumber(const unsigned int phaseIdx) const
    { return nusseltNumber_[phaseIdx]; }
 
    const Scalar schmidtNumber(const unsigned int phaseIdx) const
    { return schmidtNumber_[phaseIdx]; }
 
    const Scalar sherwoodNumber(const unsigned int phaseIdx) const
    { return sherwoodNumber_[phaseIdx]; }
 
    const Scalar characteristicLength() const
    { return characteristicLength_; }
 
    const Scalar factorEnergyTransfer() const
    { return factorEnergyTransfer_; }
 
    const Scalar factorMassTransfer() const
    { return factorMassTransfer_; }
 
private:
    NumFluidPhasesArray reynoldsNumber_;
    NumFluidPhasesArray prandtlNumber_;
    NumFluidPhasesArray nusseltNumber_;
    NumFluidPhasesArray schmidtNumber_;
    NumFluidPhasesArray sherwoodNumber_;
    Scalar characteristicLength_;
    Scalar factorEnergyTransfer_;
    Scalar factorMassTransfer_;
    InterfacialAreasArray interfacialArea_;
};
 
} // namespace Dumux
 
#endif