Commit ad8a5257 authored by Katharina Heck's avatar Katharina Heck Committed by Bernd Flemisch
Browse files

[test][nonequilibrium][2p2c] make test nonisothermal and set better boundary...

[test][nonequilibrium][2p2c] make test nonisothermal and set better boundary conditions. add reference solution for box and tpfa
parent b63eb573
......@@ -87,7 +87,11 @@
#include <dumux/porousmediumflow/2p/formulation.hh>
#include <dumux/porousmediumflow/nonisothermal/model.hh>
#include <dumux/porousmediumflow/nonisothermal/iofields.hh>
#include <dumux/porousmediumflow/nonequilibrium/model.hh>
#include <dumux/porousmediumflow/nonequilibrium/volumevariables.hh>
#include <dumux/material/fluidmatrixinteractions/2p/thermalconductivitysomerton.hh>
#include <dumux/material/fluidmatrixinteractions/2p/thermalconductivitysimplefluidlumping.hh>
#include "volumevariables.hh"
......@@ -178,6 +182,133 @@ struct IOFields<TypeTag, TTag::TwoPTwoCNI> { using type = EnergyIOFields<TwoPNCI
template<class TypeTag>
struct ThermalConductivityModel<TypeTag, TTag::TwoPTwoCNI> { using type = ThermalConductivitySomerton<GetPropType<TypeTag, Properties::Scalar>>; };
} // end namespace Properties
template<class TwoPTwoCModelTraits>
struct TwoPTwoCUnconstrainedModelTraits : public TwoPTwoCModelTraits
{
static constexpr int numConstraintEq() { return 0; }
};
namespace Properties {
//////////////////////////////////////////////////////////////////
// Type tags
//////////////////////////////////////////////////////////////////
namespace TTag {
struct TwoPTwoCNonEquil { using InheritsFrom = std::tuple<NonEquilibrium, TwoPTwoC>; };
} // end namespace TTag
/////////////////////////////////////////////////
// Properties for the non-equilibrium TwoPTwoC model
/////////////////////////////////////////////////
template<class TypeTag>
struct EquilibriumLocalResidual<TypeTag, TTag::TwoPTwoCNonEquil> { using type = CompositionalLocalResidual<TypeTag>; };
//! Set the vtk output fields specific to this model
template<class TypeTag>
struct EquilibriumIOFields<TypeTag, TTag::TwoPTwoCNonEquil> { using type = TwoPNCIOFields; };
template<class TypeTag>
struct ModelTraits<TypeTag, TTag::TwoPTwoCNonEquil>
{
private:
using EquiTraits = GetPropType<TypeTag, Properties::EquilibriumModelTraits>;
static constexpr bool enableTNE = getPropValue<TypeTag, Properties::EnableThermalNonEquilibrium>();
static constexpr bool enableCNE = getPropValue<TypeTag, Properties::EnableChemicalNonEquilibrium>();
static constexpr int numEF = getPropValue<TypeTag, Properties::NumEnergyEqFluid>();
static constexpr int numES = getPropValue<TypeTag, Properties::NumEnergyEqSolid>();
static constexpr auto nf = getPropValue<TypeTag, Properties::NusseltFormulation>();
static constexpr auto ns = getPropValue<TypeTag, Properties::SherwoodFormulation>();
using NonEquilTraits = NonEquilibriumModelTraits<EquiTraits, enableCNE, enableTNE, numEF, numES, nf, ns>;
public:
using type = NonEquilTraits;
};
//! Set equilibrium model traits
template<class TypeTag>
struct EquilibriumModelTraits<TypeTag, TTag::TwoPTwoCNonEquil>
{
private:
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using EquilibriumTraits = GetPropType<TypeTag, Properties::BaseModelTraits>;
public:
using type = TwoPTwoCUnconstrainedModelTraits<EquilibriumTraits>;
};
//! In case we do not assume full thermal non-equilibrium (e.g. only an energy balance for the solid phase and a fluid mixture) one needs a law for calculating the thermal conductivity of the fluid mixture
template<class TypeTag>
struct ThermalConductivityModel<TypeTag, TTag::TwoPTwoCNonEquil>
{
private:
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
public:
using type = ThermalConductivitySimpleFluidLumping<Scalar, getPropValue<TypeTag, Properties::NumEnergyEqFluid>()>;
};
//! Use the nonequilibrium volume variables together with the 2p2c vol vars
template<class TypeTag>
struct VolumeVariables<TypeTag, TTag::TwoPTwoCNonEquil>
{
private:
using PV = GetPropType<TypeTag, Properties::PrimaryVariables>;
using FSY = GetPropType<TypeTag, Properties::FluidSystem>;
using FST = GetPropType<TypeTag, Properties::FluidState>;
using SSY = GetPropType<TypeTag, Properties::SolidSystem>;
using SST = GetPropType<TypeTag, Properties::SolidState>;
using MT = GetPropType<TypeTag, Properties::ModelTraits>;
using PT = typename GetPropType<TypeTag, Properties::SpatialParams>::PermeabilityType;
static constexpr bool useConstraintSolver = getPropValue<TypeTag, Properties::UseConstraintSolver>();
using Traits = TwoPNCVolumeVariablesTraits<PV, FSY, FST, SSY, SST, PT, MT>;
using EquilibriumVolVars = TwoPTwoCVolumeVariables<Traits, useConstraintSolver>;
public:
using type = NonEquilibriumVolumeVariables<Traits, EquilibriumVolVars>;
};
/////////////////////////////////////////////////
// Properties for the non-equilibrium nonisothermal TwoPTwoC model which assumes thermal equilibrium but chemical nonequilibrium
/////////////////////////////////////////////////
namespace TTag {
struct TwoPTwoCNINonEquil { using InheritsFrom = std::tuple<TwoPTwoCNonEquil>; };
} // end namespace TTag
//! Set the non-isothermal model traits with the nonequilibrium model traits as isothermal traits
template<class TypeTag>
struct ModelTraits<TypeTag, TTag::TwoPTwoCNINonEquil>
{
private:
private:
using EquiTraits = GetPropType<TypeTag, Properties::EquilibriumModelTraits>;
static constexpr bool enableTNE = getPropValue<TypeTag, Properties::EnableThermalNonEquilibrium>();
static constexpr bool enableCNE = getPropValue<TypeTag, Properties::EnableChemicalNonEquilibrium>();
static constexpr int numEF = getPropValue<TypeTag, Properties::NumEnergyEqFluid>();
static constexpr int numES = getPropValue<TypeTag, Properties::NumEnergyEqSolid>();
static constexpr auto nf = getPropValue<TypeTag, Properties::NusseltFormulation>();
static constexpr auto ns = getPropValue<TypeTag, Properties::SherwoodFormulation>();
using IsothermalTraits = NonEquilibriumModelTraits<EquiTraits, enableCNE, enableTNE, numEF, numES, nf, ns>;
public:
using type = PorousMediumFlowNIModelTraits<IsothermalTraits>;
};
//! Set the equilibrium IO fields which are in that case the nonisothermal io fields
template<class TypeTag>
struct EquilibriumIOFields<TypeTag, TTag::TwoPTwoCNINonEquil>
{
private:
using NonisothermalIOFields = EnergyIOFields<TwoPNCIOFields>;
public:
using type = NonisothermalIOFields;
};
//! Somerton is used as default model to compute the effective thermal heat conductivity
template<class TypeTag>
struct ThermalConductivityModel<TypeTag, TTag::TwoPTwoCNINonEquil> { using type = ThermalConductivitySomerton<GetPropType<TypeTag, Properties::Scalar>>; };
} // end namespace Properties
} // end namespace Dumux
......
......@@ -128,6 +128,8 @@ struct TwoPNCModelTraits
static constexpr bool enableAdvection() { return true; }
static constexpr bool enableMolecularDiffusion() { return true; }
static constexpr bool enableEnergyBalance() { return false; }
static constexpr bool enableThermalNonEquilibrium() { return false; }
static constexpr bool enableChemicalNonEquilibrium() { return false; }
static constexpr bool useMoles() { return useMol; }
static constexpr bool setMoleFractionsForFirstPhase() { return setMoleFractionForFP; }
......
......@@ -19,7 +19,7 @@
/*!
* \file
* \ingroup NonEquilibriumModel
* \brief Adds I/O fields specific to non-isothermal models
* \brief Adds I/O fields specific to non-equilibrium models
*/
#ifndef DUMUX_NONEQUILBRIUM_OUTPUT_FIELDS_HH
......@@ -29,12 +29,17 @@
namespace Dumux {
template<class ModelTraits, class EquilibriumIOFields, bool enableThermalNonEquilibrium>
class NonEquilibriumIOFieldsImplementation;
template<class ModelTraits, class EquilibriumIOFields>
using NonEquilibriumIOFields = NonEquilibriumIOFieldsImplementation<ModelTraits, EquilibriumIOFields, ModelTraits::enableThermalNonEquilibrium()>;
/*!
* \ingroup NonEquilibriumModel
* \brief Adds I/O fields specific to non-isothermal models
* \brief Adds I/O fields specific to non-equilibrium models with chemical and thermal nonequilbirum or thermal non-equilibrium only
*/
template<class ModelTraits, class EquilibriumIOFields>
class NonEquilibriumIOFields
class NonEquilibriumIOFieldsImplementation<ModelTraits, EquilibriumIOFields, true>
{
public:
template <class OutputModule>
......@@ -61,6 +66,31 @@ public:
}
};
template<class ModelTraits, class EquilibriumIOFields>
class NonEquilibriumIOFieldsImplementation<ModelTraits, EquilibriumIOFields, false>
{
public:
template <class OutputModule>
static void initOutputModule(OutputModule& out)
{
using FluidSystem = typename OutputModule::VolumeVariables::FluidSystem;
EquilibriumIOFields::initOutputModule(out);
for (int i = 0; i < ModelTraits::numFluidPhases(); ++i)
{
out.addVolumeVariable( [i](const auto& v){ return v.reynoldsNumber(i); }, "reynoldsNumber_" + FluidSystem::phaseName(i) );
}
}
template <class OutputModule>
DUNE_DEPRECATED_MSG("use initOutputModule instead")
static void init(OutputModule& out)
{
initOutputModule(out);
}
};
} // end namespace Dumux
#endif
......@@ -32,18 +32,18 @@
namespace Dumux {
// forward declaration
template<class TypeTag, bool enableThermalNonEquilibrium, bool enableChemicalNonEquilibrium>
template<class TypeTag, bool enableChemicalNonEquilibrium>
class NonEquilibriumLocalResidualImplementation;
template <class TypeTag>
using NonEquilibriumLocalResidual = NonEquilibriumLocalResidualImplementation<TypeTag, GetPropType<TypeTag, Properties::ModelTraits>::enableThermalNonEquilibrium(), GetPropType<TypeTag, Properties::ModelTraits>::enableChemicalNonEquilibrium()>;
using NonEquilibriumLocalResidual = NonEquilibriumLocalResidualImplementation<TypeTag, GetPropType<TypeTag, Properties::ModelTraits>::enableChemicalNonEquilibrium()>;
/*!
* \ingroup NonEquilibriumModel
* \brief The mass conservation part of the nonequilibrium model for a model without chemical non-equilibrium
*/
template<class TypeTag>
class NonEquilibriumLocalResidualImplementation<TypeTag, true, false>: public GetPropType<TypeTag, Properties::EquilibriumLocalResidual>
class NonEquilibriumLocalResidualImplementation<TypeTag, false>: public GetPropType<TypeTag, Properties::EquilibriumLocalResidual>
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using ParentType = GetPropType<TypeTag, Properties::EquilibriumLocalResidual>;
......@@ -120,7 +120,7 @@ public:
}
//! Add advective and diffusive phase energy fluxes
EnergyLocalResidual::heatConvectionFlux(flux, fluxVars, elemVolVars,scvf, phaseIdx);
EnergyLocalResidual::heatConvectionFlux(flux, fluxVars, phaseIdx);
}
//! Add diffusive energy fluxes. For isothermal model the contribution is zero.
EnergyLocalResidual::heatConductionFlux(flux, fluxVars);
......@@ -155,12 +155,11 @@ public:
};
/*!
* \brief The mass conservation part of the nonequilibrium model for a model assuming chemical non-equilibrium and two phases
* \brief The mass conservation part of the nonequilibrium model for a model assuming chemical non-equilibrium and two phases but assuming thermal equilibrium
*/
template<class TypeTag>
class NonEquilibriumLocalResidualImplementation<TypeTag, true, true>: public GetPropType<TypeTag, Properties::EquilibriumLocalResidual>
class NonEquilibriumLocalResidualImplementation<TypeTag, true>: public GetPropType<TypeTag, Properties::EquilibriumLocalResidual>
{
using ParentType = GetPropType<TypeTag, Properties::EquilibriumLocalResidual>;
using Implementation = GetPropType<TypeTag, Properties::LocalResidual>;
......@@ -186,7 +185,6 @@ class NonEquilibriumLocalResidualImplementation<TypeTag, true, true>: public Get
static constexpr int numPhases = ModelTraits::numFluidPhases();
static constexpr int numComponents = ModelTraits::numFluidComponents();
static constexpr bool useMoles = getPropValue<TypeTag, Properties::UseMoles>();
using ComponentVector = Dune::FieldVector<Scalar, numComponents>;
static constexpr auto conti0EqIdx = Indices::conti0EqIdx;
static constexpr auto comp1Idx = FluidSystem::comp1Idx;
......@@ -194,6 +192,11 @@ class NonEquilibriumLocalResidualImplementation<TypeTag, true, true>: public Get
static constexpr auto phase0Idx = FluidSystem::phase0Idx;
static constexpr auto phase1Idx = FluidSystem::phase1Idx;
static_assert(numPhases == numComponents,
"currently chemical non-equilibrium is only available when numPhases equals numComponents");
static_assert(useMoles == true,
"chemical nonequilibrium can only be calculated based on mole fractions not mass fractions");
public:
using ParentType::ParentType;
/*!
......@@ -209,16 +212,16 @@ public:
{
NumEqVector storage(0.0);
// compute storage term of all components within all phases
// compute storage term of all components within all phases
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
{
auto eqIdx = conti0EqIdx + phaseIdx*numComponents + compIdx;
storage[eqIdx] += volVars.porosity()
* volVars.saturation(phaseIdx)
* volVars.molarDensity(phaseIdx)
* volVars.moleFraction(phaseIdx, compIdx);
* volVars.saturation(phaseIdx)
* volVars.molarDensity(phaseIdx)
* volVars.moleFraction(phaseIdx, compIdx);
}
//! The energy storage in the fluid phase with index phaseIdx
EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, phaseIdx);
......@@ -257,20 +260,21 @@ public:
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
{
// get equation index
const auto eqIdx = conti0EqIdx + phaseIdx*numComponents + compIdx;
const auto eqIdx = conti0EqIdx + phaseIdx*numComponents + compIdx;
// the physical quantities for which we perform upwinding
const auto upwindTerm = [phaseIdx, compIdx] (const auto& volVars)
{ return volVars.molarDensity(phaseIdx)*volVars.moleFraction(phaseIdx, compIdx)*volVars.mobility(phaseIdx); };
flux[eqIdx] += fluxVars.advectiveFlux(phaseIdx, upwindTerm);
flux[eqIdx] += fluxVars.advectiveFlux(phaseIdx, upwindTerm);
// diffusive fluxes (only for the component balances)
if (compIdx == phaseIdx) //do not add diffusive flux of main component, as that is not done in master as well
// do not add diffusive flux of main component, as that is not done in master as well
if (compIdx == phaseIdx)
continue;
flux[eqIdx] += diffusiveFluxes[compIdx];
flux[eqIdx] += diffusiveFluxes[compIdx];
}
//! Add advective phase energy fluxes. For isothermal model the contribution is zero.
EnergyLocalResidual::heatConvectionFlux(flux, fluxVars, elemVolVars,scvf, phaseIdx);
EnergyLocalResidual::heatConvectionFlux(flux, fluxVars, phaseIdx);
}
//! Add diffusive energy fluxes. For isothermal model the contribution is zero.
......@@ -294,104 +298,72 @@ public:
const ElementVolumeVariables& elemVolVars,
const SubControlVolume &scv) const
{
NumEqVector source(0.0);
NumEqVector source(0.0);
// In the case of a kinetic consideration, mass transfer
// between phases is realized via source terms there is a
// balance equation for each component in each phase
const auto& volVars = elemVolVars[scv];
ComponentVector componentIntoPhaseMassTransfer[numPhases];
#define FUNKYMASSTRANSFER 0
#if FUNKYMASSTRANSFER
const Scalar mu_nPhaseNCompEquil = volVars.chemicalPotentialEquil(phase1Idx, comp1Idx) ; // very 2p2c
const Scalar mu_wPhaseWCompEquil = volVars.chemicalPotentialEquil(phase0Idx, comp0Idx); // very 2p2c
const Scalar mu_wPhaseNComp = volVars.chemicalPotential(phase0Idx, comp1Idx) ; // very 2p2c
const Scalar mu_nPhaseWComp = volVars.chemicalPotential(phase1Idx, comp0Idx); // very 2p2c
Valgrind::CheckDefined(mu_nPhaseNCompEquil);
Valgrind::CheckDefined(mu_wPhaseWCompEquil);
Valgrind::CheckDefined(mu_wPhaseNComp);
Valgrind::CheckDefined(mu_nPhaseWComp);
std::array<std::array<Scalar, numComponents>, numPhases> componentIntoPhaseMassTransfer = {{{0.0},{0.0}}};
//get characteristic length
const Scalar characteristicLength = volVars.characteristicLength() ;
const Scalar temperature = volVars.temperatureFluid(phase0Idx);
const Scalar pn = volVars.pressure(phase1Idx);
const Scalar henry = FluidSystem::henry(temperature) ;
const Scalar gradNinWApprox = ( mu_wPhaseNComp - mu_nPhaseNCompEquil) / characteristicLength; // very 2p2c // 1. / henry *
const Scalar gradWinNApprox = ( mu_nPhaseWComp - mu_wPhaseWCompEquil) / characteristicLength; // very 2p2c // 1. / pn *
#else // FUNKYMASSTRANSFER
Scalar x[numPhases][numComponents]; // mass fractions in wetting phase
for(int phaseIdx=0; phaseIdx<numPhases; ++phaseIdx){
for (int compIdx=0; compIdx< numComponents; ++ compIdx){
x[phaseIdx][compIdx] = volVars.moleFraction(phaseIdx, compIdx);
}
}
Valgrind::CheckDefined(x);
// "equilibrium" values: calculated in volume variables
const Scalar x_wPhaseNCompEquil = volVars.xEquil(phase0Idx, comp1Idx) ; // very 2p2c
const Scalar x_nPhaseWCompEquil = volVars.xEquil(phase1Idx, comp0Idx); // very 2p2c
Valgrind::CheckDefined(x_wPhaseNCompEquil);
Valgrind::CheckDefined(x_nPhaseWCompEquil);
const Scalar characteristicLength = volVars.characteristicLength() ;
const Scalar gradNinWApprox = (x[phase0Idx][comp1Idx] - x_wPhaseNCompEquil) / characteristicLength; // very 2p2c
const Scalar gradWinNApprox = (x[phase1Idx][comp0Idx] - x_nPhaseWCompEquil) / characteristicLength; // very 2p2c
#endif
Scalar phaseDensity[numPhases];
for(int phaseIdx=0; phaseIdx<numPhases; ++phaseIdx){
phaseDensity[phaseIdx] = volVars.molarDensity(phaseIdx);
}
// diffusion coefficients in wetting phase
const Scalar diffCoeffNinW = volVars.diffusionCoefficient(phase0Idx, comp1Idx) ;
Valgrind::CheckDefined(diffCoeffNinW);
// diffusion coefficients in non-wetting phase
const Scalar diffCoeffWinN = volVars.diffusionCoefficient(phase1Idx, comp0Idx) ;
Valgrind::CheckDefined(diffCoeffWinN);
const Scalar factorMassTransfer = volVars.factorMassTransfer() ;
const Scalar awn = volVars.interfacialArea(phase0Idx, phase1Idx);
const Scalar sherwoodWPhase = volVars.sherwoodNumber(phase0Idx);
const Scalar sherwoodNPhase = volVars.sherwoodNumber(phase1Idx);
for(int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
const Scalar sherwoodNumber = volVars.sherwoodNumber(phaseIdx);
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
{
if (compIdx <= numPhases)
{
//we only have to do the source calculation for one component going into the other phase as for each component: n_wn -> n_n = - n_n -> n_wn
if (phaseIdx == compIdx)
continue;
// actual diffusion is always calculated for eq's 2,3
// Eq's 1,4 have to be the same with different sign, because no mass is accumulated in the interface
// i.e. automatically conserving mass that mvoes across the interface
const Scalar xNonEquil = volVars.moleFraction(phaseIdx, compIdx);
const Scalar nCompIntoWPhase = - factorMassTransfer * gradNinWApprox * awn * phaseDensity[phase0Idx] * diffCoeffNinW * sherwoodWPhase;
const Scalar nCompIntoNPhase = - nCompIntoWPhase ;
const Scalar wCompIntoNPhase = - factorMassTransfer * gradWinNApprox * awn * phaseDensity[phase1Idx] * diffCoeffWinN * sherwoodNPhase;
const Scalar wCompIntoWPhase = - wCompIntoNPhase ;
//additionally get equilibrium values from volume variables
const Scalar xEquil = volVars.xEquil(phaseIdx, compIdx);
//get the diffusion coefficient
const Scalar diffCoeff = volVars.diffusionCoefficient(phaseIdx, compIdx);
componentIntoPhaseMassTransfer[phase0Idx][comp1Idx] = nCompIntoWPhase;
componentIntoPhaseMassTransfer[phase1Idx][comp1Idx] = nCompIntoNPhase;
componentIntoPhaseMassTransfer[phase1Idx][comp0Idx] = wCompIntoNPhase;
componentIntoPhaseMassTransfer[phase0Idx][comp0Idx] = wCompIntoWPhase;
//now compute the flux
const Scalar compFluxIntoOtherPhase = factorMassTransfer * (xEquil-xNonEquil)*characteristicLength * awn * volVars.molarDensity(phaseIdx) * diffCoeff * sherwoodNumber;
Valgrind::CheckDefined(componentIntoPhaseMassTransfer);
componentIntoPhaseMassTransfer[phaseIdx][compIdx] += compFluxIntoOtherPhase;
componentIntoPhaseMassTransfer[compIdx][compIdx] += -compFluxIntoOtherPhase;
}
}
}
// Actually add the mass transfer to the sources which might
// exist externally
for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
for (int compIdx = 0; compIdx < numComponents; ++ compIdx) {
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
{
const unsigned int eqIdx = conti0EqIdx + compIdx + phaseIdx*numComponents;
source[eqIdx] += componentIntoPhaseMassTransfer[phaseIdx][compIdx] ;
source[eqIdx] += componentIntoPhaseMassTransfer[phaseIdx][compIdx];
using std::isfinite;
if (!isfinite(source[eqIdx]))
DUNE_THROW(NumericalProblem, "Calculated non-finite source");
using std::isfinite;
if (!isfinite(source[eqIdx]))
DUNE_THROW(NumericalProblem, "Calculated non-finite source");
}
}
// Call the (kinetic) Energy module, for the source term.
// it has to be called from here, because the mass transfered has to be known.
EnergyLocalResidual::computeSourceEnergy(source,
element,
fvGeometry,
elemVolVars,
scv);
if (ModelTraits::enableThermalNonEquilibrium())
{
// Call the (kinetic) Energy module, for the source term.
// it has to be called from here, because the mass transfered has to be known.
EnergyLocalResidual::computeSourceEnergy(source,
element,
fvGeometry,
elemVolVars,
scv);
}
// add contributions from volume flux sources
source += problem.source(element, fvGeometry, elemVolVars, scv);
......
......@@ -71,6 +71,8 @@ struct NonEquilibriumModelTraits : public ET
static constexpr NusseltFormulation nusseltFormulation() { return nf; }
static constexpr SherwoodFormulation sherwoodFormulation() { return sf; }
static_assert(!(ET::enableEnergyBalance() && therm), "It is not possible to use a nonisothermal model assuming local thermal equilibrium in combination with a model using thermal non-equilibrium");
using Indices = NonEquilbriumIndices<typename ET::Indices, numEnergyEqFluid(), numEnergyEqSolid(), numEq()>;
};
......
......@@ -64,7 +64,7 @@ class EnergyLocalResidualNonEquilibrium<TypeTag, 1/*numEnergyEqFluid*/>
static constexpr auto energyEq0Idx = Indices::energyEq0Idx;
static constexpr auto energyEqSolidIdx = Indices::energyEqSolidIdx;
static constexpr auto numPhases = ModelTraits::numFluidPhases();
static constexpr auto numPhases = ModelTraits::numFluidPhases();
static constexpr auto numComponents = ModelTraits::numFluidComponents();
public:
......@@ -76,9 +76,9 @@ public:
{
//in case we have one energy equation for more than one fluid phase, add up parts on the one energy equation
storage[energyEq0Idx] += volVars.porosity()
* volVars.density(phaseIdx)
* volVars.internalEnergy(phaseIdx)
* volVars.saturation(phaseIdx);
* volVars.density(phaseIdx)
* volVars.internalEnergy(phaseIdx)
* volVars.saturation(phaseIdx);
}
......@@ -88,28 +88,20 @@ public:
const SubControlVolume& scv,
const VolumeVariables& volVars)
{
// heat conduction for the fluid phases
for(int sPhaseIdx=0; sPhaseIdx<numEnergyEqSolid; ++sPhaseIdx)
{
// heat conduction for the fluid phases
for(int sPhaseIdx = 0; sPhaseIdx<numEnergyEqSolid; ++sPhaseIdx)
{
storage[energyEqSolidIdx+sPhaseIdx] += volVars.temperatureSolid()
* volVars.solidHeatCapacity()
* volVars.solidDensity()
* (1.0 - volVars.porosity());
}
* volVars.solidHeatCapacity()
* volVars.solidDensity()
* (1.0 - volVars.porosity());
}
}
// this is to make nonequilibrium work with compositional local residual, compositional calls that for non-isothermal models
static void heatConvectionFlux(NumEqVector& flux,
FluxVariables& fluxVars,
int phaseIdx)
{}
//! The advective phase energy fluxes
//! The advective phase energy fluxes
static void heatConvectionFlux(NumEqVector& flux,
FluxVariables& fluxVars,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf,
int phaseIdx)
{
auto upwindTerm = [phaseIdx](const auto& volVars)
......@@ -120,8 +112,12 @@ public:
//now add the diffusive part
const auto diffusiveFluxes = fluxVars.molecularDiffusionFlux(phaseIdx);
const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
const auto& elemVolVars = fluxVars.elemVolVars();
const auto& scvf = fluxVars.scvFace();
const auto& insideScv = fluxVars.fvGeometry().scv(scvf.insideScvIdx());
const auto& outsideScv = fluxVars.fvGeometry().scv(scvf.outsideScvIdx());
const auto& insideVolVars = elemVolVars[insideScv];
const auto& outsideVolVars = elemVolVars[outsideScv];
auto insideEnthalpy = insideVolVars.enthalpy(phaseIdx);
auto outsideEnthalpy = outsideVolVars.enthalpy(phaseIdx);
......@@ -130,7 +126,7 @@ public:
//no diffusion of the main component, this is a hack to use normal fick's law which computes both diffusions (main and component). We only add the part from the component here
if (phaseIdx == compIdx)
continue;
//we need the upwind enthapy. Even better would be the componentEnthalpy
//we need the upwind enthalpy. Even better would be the componentEnthalpy
auto enthalpy = 0.0;
if (diffusiveFluxes[compIdx] > 0)
enthalpy += insideEnthalpy;
......@@ -147,12 +143,12 @@ public:
//in case we have one energy equation for more than one fluid phase we use an effective law in the nonequilibrium fourierslaw
flux[energyEq0Idx] += fluxVars.heatConductionFlux(0);
//heat conduction for the solid phases
for(int sPhaseIdx=0; sPhaseIdx