releases/3.4/a18904_source.html

// -*- 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 3 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/>.   *

 *****************************************************************************/

#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 awn = fluidMatrixInteraction.wettingNonwettingInterface().area(Sw, pc);

        const auto ans = fluidMatrixInteraction.nonwettingSolidInterface().area(Sw, pc);

        interfacialArea_[phase0Idx][phase1Idx] = awn;

        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 - ans;

        }

        else

            interfacialArea_[phase0Idx][sPhaseIdx] = fluidMatrixInteraction.wettingSolidInterface().area(Sw, pc);


        interfacialArea_[sPhaseIdx][phase0Idx] = interfacialArea_[phase0Idx][sPhaseIdx];

        interfacialArea_[sPhaseIdx][sPhaseIdx] = 0.0;


        interfacialArea_[phase1Idx][sPhaseIdx] = ans;

        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 ans = 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 - ans;

        }

        else

            interfacialArea_[phase0Idx][sPhaseIdx] = fluidMatrixInteraction.wettingSolidInterface().area(Sw, pc);


        interfacialArea_[sPhaseIdx][phase0Idx] = interfacialArea_[phase0Idx][sPhaseIdx];

        interfacialArea_[sPhaseIdx][sPhaseIdx] = 0.;


        interfacialArea_[phase1Idx][sPhaseIdx] = ans;

        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

dimensionlessnumbers.hh
Collection of functions, calculating dimensionless numbers.

parameters.hh
The infrastructure to retrieve run-time parameters from Dune::ParameterTrees.

Dumux
Definition: adapt.hh:29

Dumux::IOName::density
std::string density(int phaseIdx) noexcept
I/O name of density for multiphase systems.
Definition: name.hh:65

Dumux::IOName::porosity
std::string porosity() noexcept
I/O name of porosity.
Definition: name.hh:139

Dumux::DimensionlessNumbers
Collection of functions which calculate dimensionless numbers. Each number has it's own function....
Definition: dimensionlessnumbers.hh:65

Dumux::NonEquilibriumVolumeVariablesImplementation
This class contains the volume variables required for the modules which require the specific interfac...
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:49

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, false, true, 2 >::updateInterfacialArea
void updateInterfacialArea(const ElemSol &elemSol, const FluidState &fluidState, const ParameterCache &paramCache, const Problem &problem, const Element &element, const Scv &scv)
Updates the volume specific interfacial area [m^2 / m^3] between the phases.
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:173

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, false, true, 2 >::factorEnergyTransfer
const Scalar factorEnergyTransfer() const
access function pre factor energy transfer
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:239

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, false, true, 2 >::prandtlNumber
const Scalar prandtlNumber(const unsigned int phaseIdx) const
access function Prandtl Number
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:227

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, false, true, 2 >::interfacialArea
const Scalar interfacialArea(const unsigned int phaseIIdx, const unsigned int phaseJIdx) const
The specific interfacial area between two fluid phases [m^2 / m^3].
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:215

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, false, true, 2 >::Indices
typename ModelTraits::Indices Indices
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:92

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, false, true, 2 >::update
void update(const ElemSol &elemSol, const Problem &problem, const Element &element, const Scv &scv)
Update all quantities for a given control volume.
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:102

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, false, true, 2 >::reynoldsNumber
const Scalar reynoldsNumber(const unsigned int phaseIdx) const
access function Reynolds Number
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:223

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, false, true, 2 >::updateDimLessNumbers
void updateDimLessNumbers(const ElemSol &elemSol, const FluidState &fluidState, const ParameterCache &paramCache, const Problem &problem, const Element &element, const Scv &scv)
Updates dimensionless numbers.
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:127

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, false, true, 2 >::nusseltNumber
const Scalar nusseltNumber(const unsigned int phaseIdx) const
access function Nusselt Number
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:231

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, false, true, 2 >::characteristicLength
const Scalar characteristicLength() const
access function characteristic length
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:235

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, false, true, 2 >::FluidState
typename Traits::FluidState FluidState
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:91

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, false, true, 1 >::Indices
typename ModelTraits::Indices Indices
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:278

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, false, true, 1 >::fluidSolidInterfacialArea
const Scalar fluidSolidInterfacialArea() const
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:395

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, false, true, 1 >::characteristicLength
const Scalar characteristicLength() const
access function characteristic length
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:388

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, false, true, 1 >::reynoldsNumber
const Scalar reynoldsNumber(const unsigned int phaseIdx) const
access function Reynolds Number
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:376

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, false, true, 1 >::updateDimLessNumbers
void updateDimLessNumbers(const ElemSol &elemSol, const FluidState &fluidState, const ParameterCache &paramCache, const Problem &problem, const Element &element, const Scv &scv)
Updates dimensionless numbers.
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:318

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, false, true, 1 >::updateInterfacialArea
void updateInterfacialArea(const ElemSol &elemSol, const FluidState &fluidState, const ParameterCache &paramCache, const Problem &problem, const Element &element, const Scv &scv)
Updates the volume specific interfacial area [m^2 / m^3] between the solid and the fluid phase.
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:364

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, false, true, 1 >::nusseltNumber
const Scalar nusseltNumber(const unsigned int phaseIdx) const
access function Nusselt Number
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:384

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, false, true, 1 >::prandtlNumber
const Scalar prandtlNumber(const unsigned int phaseIdx) const
access function Prandtl Number
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:380

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, false, true, 1 >::FluidState
typename Traits::FluidState FluidState
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:279

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, false, true, 1 >::factorEnergyTransfer
const Scalar factorEnergyTransfer() const
access function pre factor energy transfer
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:392

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, false, true, 1 >::update
void update(const ElemSol &elemSol, const Problem &problem, const Element &element, const Scv &scv)
Update all quantities for a given control volume.
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:291

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, true, false, 0 >::update
void update(const ElemSol &elemSol, const Problem &problem, const Element &element, const Scv &scv)
Update all quantities for a given control volume.
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:449

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, true, false, 0 >::reynoldsNumber
const Scalar reynoldsNumber(const unsigned int phaseIdx) const
access function Reynolds Number
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:549

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, true, false, 0 >::updateDimLessNumbers
void updateDimLessNumbers(const ElemSol &elemSol, const FluidState &fluidState, const ParameterCache &paramCache, const Problem &problem, const Element &element, const Scv &scv)
Updates the volume specific interfacial area [m^2 / m^3] between the phases.
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:474

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, true, false, 0 >::Indices
typename ModelTraits::Indices Indices
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:438

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, true, false, 0 >::characteristicLength
const Scalar characteristicLength() const
access function characteristic length
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:561

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, true, false, 0 >::sherwoodNumber
const Scalar sherwoodNumber(const unsigned int phaseIdx) const
access function Sherwood Number
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:557

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, true, false, 0 >::factorMassTransfer
const Scalar factorMassTransfer() const
access function pre factor mass transfer
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:565

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, true, false, 0 >::interfacialArea
const Scalar interfacialArea(const unsigned int phaseIIdx, const unsigned int phaseJIdx) const
The specific interfacial area between two fluid phases [m^2 / m^3].
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:540

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, true, false, 0 >::updateInterfacialArea
void updateInterfacialArea(const ElemSol &elemSol, const FluidState &fluidState, const ParameterCache &paramCache, const Problem &problem, const Element &element, const Scv &scv)
Updates the volume specific interfacial area [m^2 / m^3] between the phases.
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:521

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, true, false, 0 >::schmidtNumber
const Scalar schmidtNumber(const unsigned int phaseIdx) const
access function Schmidt Number
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:553

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, true, true, 2 >::characteristicLength
const Scalar characteristicLength() const
access function characteristic length
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:775

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, true, true, 2 >::sherwoodNumber
const Scalar sherwoodNumber(const unsigned int phaseIdx) const
access function Sherwood Number
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:771

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, true, true, 2 >::Indices
typename ModelTraits::Indices Indices
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:610

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, true, true, 2 >::update
void update(const ElemSol &elemSol, const Problem &problem, const Element &element, const Scv &scv)
Update all quantities for a given control volume.
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:622

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, true, true, 2 >::factorMassTransfer
const Scalar factorMassTransfer() const
access function pre factor mass transfer
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:783

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, true, true, 2 >::factorEnergyTransfer
const Scalar factorEnergyTransfer() const
access function pre factor energy transfer
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:779

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, true, true, 2 >::schmidtNumber
const Scalar schmidtNumber(const unsigned int phaseIdx) const
access function Schmidt Number
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:767

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, true, true, 2 >::nusseltNumber
const Scalar nusseltNumber(const unsigned int phaseIdx) const
access function Nusselt Number
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:763

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, true, true, 2 >::interfacialArea
const Scalar interfacialArea(const unsigned int phaseIIdx, const unsigned int phaseJIdx) const
The specific interfacial area between two fluid phases [m^2 / m^3].
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:747

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, true, true, 2 >::prandtlNumber
const Scalar prandtlNumber(const unsigned int phaseIdx) const
access function Prandtl Number
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:759

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, true, true, 2 >::updateDimLessNumbers
void updateDimLessNumbers(const ElemSol &elemSol, const FluidState &fluidState, const ParameterCache &paramCache, const Problem &problem, const Element &element, const Scv &scv)
Updates the volume specific interfacial area [m^2 / m^3] between the phases.
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:647

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, true, true, 2 >::reynoldsNumber
const Scalar reynoldsNumber(const unsigned int phaseIdx) const
access function Reynolds Number
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:755

Dumux::NonEquilibriumVolumeVariablesImplementation< Traits, EquilibriumVolumeVariables, true, true, 2 >::updateInterfacialArea
void updateInterfacialArea(const ElemSol &elemSol, const FluidState &fluidState, const ParameterCache &paramCache, const Problem &problem, const Element &element, const Scv &scv)
Updates the volume specific interfacial area [m^2 / m^3] between the phases.
Definition: porousmediumflow/nonequilibrium/volumevariables.hh:703