Commit 71325638 authored by Kilian Weishaupt's avatar Kilian Weishaupt Committed by Ned Coltman
Browse files

[2p2c][2pnc] Change volVars for effective diffusivities

* add in volVars of 2p2c and 2pnc.
* add effective diffusivity model to volumevariables traits
* add in wettingPhaseIdx()
parent f4d8a7c6
......@@ -145,13 +145,15 @@ private:
using SST = GetPropType<TypeTag, Properties::SolidState>;
using MT = GetPropType<TypeTag, Properties::ModelTraits>;
using PT = typename GetPropType<TypeTag, Properties::SpatialParams>::PermeabilityType;
using DT = GetPropType<TypeTag, Properties::MolecularDiffusionType>;
using EDM = GetPropType<TypeTag, Properties::EffectiveDiffusivityModel>;
static_assert(FSY::numComponents == 2, "Only fluid systems with 2 components are supported by the 2p2c model!");
static_assert(FSY::numPhases == 2, "Only fluid systems with 2 phases are supported by the 2p2c model!");
static constexpr bool useConstraintSolver = getPropValue<TypeTag, Properties::UseConstraintSolver>();
using Traits = TwoPNCVolumeVariablesTraits<PV, FSY, FST, SSY, SST, PT, MT>;
using Traits = TwoPNCVolumeVariablesTraits<PV, FSY, FST, SSY, SST, PT, MT, DT, EDM>;
public:
using type = TwoPTwoCVolumeVariables<Traits, useConstraintSolver>;
};
......@@ -174,6 +176,28 @@ public:
using type = PorousMediumFlowNIModelTraits<IsothermalTraits>;
};
//! Set the volume variables property
template<class TypeTag>
struct VolumeVariables<TypeTag, TTag::TwoPTwoCNI>
{
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;
using DT = GetPropType<TypeTag, Properties::MolecularDiffusionType>;
using EDM = GetPropType<TypeTag, Properties::EffectiveDiffusivityModel>;
using ETCM = GetPropType< TypeTag, Properties:: ThermalConductivityModel>;
using Traits = TwoPNCNIVolumeVariablesTraits<PV, FSY, FST, SSY, SST, PT, MT, DT, EDM, ETCM>;
public:
using type = TwoPTwoCVolumeVariables<Traits>;
};
//! Set non-isothermal output fields
template<class TypeTag>
struct IOFields<TypeTag, TTag::TwoPTwoCNI> { using type = EnergyIOFields<TwoPNCIOFields>; };
......@@ -259,10 +283,12 @@ private:
using SST = GetPropType<TypeTag, Properties::SolidState>;
using MT = GetPropType<TypeTag, Properties::ModelTraits>;
using PT = typename GetPropType<TypeTag, Properties::SpatialParams>::PermeabilityType;
using DT = GetPropType<TypeTag, Properties::MolecularDiffusionType>;
using EDM = GetPropType<TypeTag, Properties::EffectiveDiffusivityModel>;
static constexpr bool useConstraintSolver = getPropValue<TypeTag, Properties::UseConstraintSolver>();
using Traits = TwoPNCVolumeVariablesTraits<PV, FSY, FST, SSY, SST, PT, MT>;
using Traits = TwoPNCVolumeVariablesTraits<PV, FSY, FST, SSY, SST, PT, MT, DT, EDM>;
using EquilibriumVolVars = TwoPTwoCVolumeVariables<Traits, useConstraintSolver>;
public:
using type = NonEquilibriumVolumeVariables<Traits, EquilibriumVolVars>;
......@@ -305,6 +331,30 @@ public:
using type = NonisothermalIOFields;
};
//! Use the nonequilibrium volume variables together with the 2p2c vol vars
template<class TypeTag>
struct VolumeVariables<TypeTag, TTag::TwoPTwoCNINonEquil>
{
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;
using DT = GetPropType<TypeTag, Properties::MolecularDiffusionType>;
using EDM = GetPropType<TypeTag, Properties::EffectiveDiffusivityModel>;
using ETCM = GetPropType< TypeTag, Properties:: ThermalConductivityModel>;
static constexpr bool useConstraintSolver = getPropValue<TypeTag, Properties::UseConstraintSolver>();
using Traits = TwoPNCNIVolumeVariablesTraits<PV, FSY, FST, SSY, SST, PT, MT, DT, EDM, ETCM>;
using EquilibriumVolVars = TwoPTwoCVolumeVariables<Traits, useConstraintSolver>;
public:
using type = NonEquilibriumVolumeVariables<Traits, EquilibriumVolVars>;
};
//! 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>>; };
......
......@@ -100,6 +100,8 @@ class TwoPTwoCVolumeVariablesBase
using PermeabilityType = typename Traits::PermeabilityType;
using ComputeFromReferencePhase = Dumux::ComputeFromReferencePhase<Scalar, typename Traits::FluidSystem>;
using MiscibleMultiPhaseComposition = Dumux::MiscibleMultiPhaseComposition< Scalar, typename Traits::FluidSystem >;
using EffDiffModel = typename Traits::EffectiveDiffusivityModel;
using DiffusionCoefficients = typename Traits::DiffusionType::template DiffusionCoefficientsContainer<2, 2>;
public:
//! The type of the object returned by the fluidState() method
using FluidState = typename Traits::FluidState;
......@@ -145,21 +147,25 @@ public:
using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
const auto& matParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
const int wPhaseIdx = problem.spatialParams().template wettingPhase<FluidSystem>(element, scv, elemSol);
const int nPhaseIdx = 1 - wPhaseIdx;
const int nPhaseIdx = 1 - wPhaseIdx_;
// relative permeabilities -> require wetting phase saturation as parameter!
relativePermeability_[wPhaseIdx] = MaterialLaw::krw(matParams, saturation(wPhaseIdx));
relativePermeability_[nPhaseIdx] = MaterialLaw::krn(matParams, saturation(wPhaseIdx));
relativePermeability_[wPhaseIdx_] = MaterialLaw::krw(matParams, saturation(wPhaseIdx_));
relativePermeability_[nPhaseIdx] = MaterialLaw::krn(matParams, saturation(wPhaseIdx_));
// binary diffusion coefficients
diffCoeff_[phase0Idx] = FluidSystem::binaryDiffusionCoefficient(fluidState_, paramCache, phase0Idx, comp0Idx, comp1Idx);
diffCoeff_[phase1Idx] = FluidSystem::binaryDiffusionCoefficient(fluidState_, paramCache, phase1Idx, comp0Idx, comp1Idx);
diffCoeff_(phase0Idx, comp0Idx, comp1Idx) = FluidSystem::binaryDiffusionCoefficient(fluidState_, paramCache, phase0Idx, comp0Idx, comp1Idx);
diffCoeff_(phase1Idx, comp1Idx, comp0Idx) = FluidSystem::binaryDiffusionCoefficient(fluidState_, paramCache, phase1Idx, comp0Idx, comp1Idx);
// porosity & permeabilty
updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, numFluidComps);
EnergyVolVars::updateSolidEnergyParams(elemSol, problem, element, scv, solidState_);
permeability_ = problem.spatialParams().permeability(element, scv, elemSol);
effectiveDiffCoeff_(phase0Idx, comp0Idx, comp1Idx) = EffDiffModel::effectiveDiffusivity(*this, phase0Idx, comp0Idx, comp1Idx);
effectiveDiffCoeff_(phase1Idx, comp1Idx, comp0Idx) = EffDiffModel::effectiveDiffusivity(*this, phase1Idx, comp1Idx, comp0Idx);
EnergyVolVars::updateEffectiveThermalConductivity();
}
/*!
......@@ -190,8 +196,8 @@ public:
using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
const int wPhaseIdx = problem.spatialParams().template wettingPhase<FluidSystem>(element, scv, elemSol);
fluidState.setWettingPhase(wPhaseIdx);
wPhaseIdx_ = problem.spatialParams().template wettingPhase<FluidSystem>(element, scv, elemSol);
fluidState.setWettingPhase(wPhaseIdx_);
// set the saturations
if (phasePresence == firstPhaseOnly)
......@@ -221,17 +227,17 @@ public:
DUNE_THROW(Dune::InvalidStateException, "Invalid phase presence.");
// set pressures of the fluid phases
pc_ = MaterialLaw::pc(materialParams, fluidState.saturation(wPhaseIdx));
pc_ = MaterialLaw::pc(materialParams, fluidState.saturation(wPhaseIdx_));
if (formulation == TwoPFormulation::p0s1)
{
fluidState.setPressure(phase0Idx, priVars[pressureIdx]);
fluidState.setPressure(phase1Idx, (wPhaseIdx == phase0Idx) ? priVars[pressureIdx] + pc_
fluidState.setPressure(phase1Idx, (wPhaseIdx_ == phase0Idx) ? priVars[pressureIdx] + pc_
: priVars[pressureIdx] - pc_);
}
else
{
fluidState.setPressure(phase1Idx, priVars[pressureIdx]);
fluidState.setPressure(phase0Idx, (wPhaseIdx == phase0Idx) ? priVars[pressureIdx] - pc_
fluidState.setPressure(phase0Idx, (wPhaseIdx_ == phase0Idx) ? priVars[pressureIdx] - pc_
: priVars[pressureIdx] + pc_);
}
}
......@@ -371,26 +377,52 @@ public:
/*!
* \brief Returns the binary diffusion coefficients for a phase in \f$[m^2/s]\f$.
*/
[[deprecated("Signature deprecated. Use diffusionCoefficient(phaseIdx, compIIdx, compJIdx)!")]]
Scalar diffusionCoefficient(int phaseIdx, int compIdx) const
{
if(phaseIdx == compIdx)
DUNE_THROW(Dune::InvalidStateException, "Diffusion coefficient called for phaseIdx = compIdx");
else
return diffCoeff_[phaseIdx];
}
{ return diffCoeff_(phaseIdx, 0, 0); }
/*!
* \brief Returns the binary diffusion coefficients for a phase in \f$[m^2/s]\f$.
*/
Scalar diffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
{ return diffCoeff_(phaseIdx, compIIdx, compJIdx); }
/*!
* \brief Returns the effective diffusion coefficients for a phase in \f$[m^2/s]\f$.
*/
[[deprecated("Signature deprecated. Use effectiveDiffusivity(phaseIdx, compIIdx, compJIdx)!")]]
Scalar effectiveDiffusivity(int phaseIdx, int compI) const
{ return effectiveDiffCoeff_(phaseIdx, 0, 0); }
/*!
* \brief Returns the effective diffusion coefficients for a phase in \f$[m^2/s]\f$.
*/
Scalar effectiveDiffusivity(int phaseIdx, int compIIdx, int compJIdx) const
{ return effectiveDiffCoeff_(phaseIdx, compIIdx, compJIdx); }
/*!
* \brief Returns the wetting phase index
*/
const int wettingPhaseIdx() const
{ return wPhaseIdx_; }
private:
FluidState fluidState_;
SolidState solidState_;
int wPhaseIdx_;
Scalar pc_; // The capillary pressure
PermeabilityType permeability_; // Effective permeability within the control volume
// Relative permeability within the control volume
std::array<Scalar, ModelTraits::numFluidPhases()> relativePermeability_;
// Binary diffusion coefficients for the phases
std::array<Scalar, ModelTraits::numFluidPhases()> diffCoeff_;
// Binary diffusion coefficient
DiffusionCoefficients diffCoeff_;
// Effective diffusion coefficients for the phases
DiffusionCoefficients effectiveDiffCoeff_;
protected:
const Impl &asImp_() const { return *static_cast<const Impl*>(this); }
......
......@@ -146,8 +146,10 @@ struct TwoPNCModelTraits
* \tparam FST The fluid state type
* \tparam PT The type used for permeabilities
* \tparam MT The model traits
*/
template<class PV, class FSY, class FST, class SSY, class SST, class PT, class MT>
* \tparam DT The diffusion type
* \tparam EDM The effective diffusivity model
*/
template<class PV, class FSY, class FST, class SSY, class SST, class PT, class MT, class DT, class EDM>
struct TwoPNCVolumeVariablesTraits
{
using PrimaryVariables = PV;
......@@ -157,6 +159,27 @@ struct TwoPNCVolumeVariablesTraits
using SolidState = SST;
using PermeabilityType = PT;
using ModelTraits = MT;
using DiffusionType = DT;
using EffectiveDiffusivityModel = EDM;
};
/*!
* \ingroup TwoPNCModel
* \brief Traits class for the volume variables of the single-phase model.
*
* \tparam PV The type used for primary variables
* \tparam FSY The fluid system type
* \tparam FST The fluid state type
* \tparam PT The type used for permeabilities
* \tparam MT The model traits
* \tparam DT The diffusion type
* \tparam EDM The effective diffusivity model
* \tparam ETCM The effective thermal conductivity model
*/
template<class PV, class FSY, class FST, class SSY, class SST, class PT, class MT, class DT, class EDM, class ETCM>
struct TwoPNCNIVolumeVariablesTraits : public TwoPNCVolumeVariablesTraits<PV, FSY, FST, SSY, SST, PT, MT, DT, EDM>
{
using EffectiveThermalConductivityModel = ETCM;
};
namespace Properties {
......@@ -195,8 +218,10 @@ private:
using SST = GetPropType<TypeTag, Properties::SolidState>;
using MT = GetPropType<TypeTag, Properties::ModelTraits>;
using PT = typename GetPropType<TypeTag, Properties::SpatialParams>::PermeabilityType;
using DT = GetPropType<TypeTag, Properties::MolecularDiffusionType>;
using EDM = GetPropType<TypeTag, Properties::EffectiveDiffusivityModel>;
using Traits = TwoPNCVolumeVariablesTraits<PV, FSY, FST, SSY, SST, PT, MT>;
using Traits = TwoPNCVolumeVariablesTraits<PV, FSY, FST, SSY, SST, PT, MT, DT, EDM>;
public:
using type = TwoPNCVolumeVariables<Traits>;
};
......@@ -267,6 +292,27 @@ public:
using type = PorousMediumFlowNIModelTraits<IsothermalTraits>;
};
//! Set the volume variables property
template<class TypeTag>
struct VolumeVariables<TypeTag, TTag::TwoPNCNI>
{
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;
using DT = GetPropType<TypeTag, Properties::MolecularDiffusionType>;
using EDM = GetPropType<TypeTag, Properties::EffectiveDiffusivityModel>;
using ETCM = GetPropType< TypeTag, Properties:: ThermalConductivityModel>;
using Traits = TwoPNCNIVolumeVariablesTraits<PV, FSY, FST, SSY, SST, PT, MT, DT, EDM, ETCM>;
public:
using type = TwoPNCVolumeVariables<Traits>;
};
//! Set non-isothermal output fields
template<class TypeTag>
struct IOFields<TypeTag, TTag::TwoPNCNI> { using type = EnergyIOFields<TwoPNCIOFields>; };
......
......@@ -90,6 +90,7 @@ class TwoPNCVolumeVariables
using MiscibleMultiPhaseComposition = Dumux::MiscibleMultiPhaseComposition<Scalar, FS>;
using ComputeFromReferencePhase = Dumux::ComputeFromReferencePhase<Scalar, FS>;
using EffDiffModel = typename Traits::EffectiveDiffusivityModel;
public:
......@@ -144,37 +145,44 @@ public:
using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
const auto& matParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
const int wPhaseIdx = problem.spatialParams().template wettingPhase<FluidSystem>(element, scv, elemSol);
const int nPhaseIdx = 1 - wPhaseIdx;
const int nPhaseIdx = 1 - wPhaseIdx_;
// mobilities -> require wetting phase saturation as parameter!
mobility_[wPhaseIdx] = MaterialLaw::krw(matParams, saturation(wPhaseIdx))/fluidState_.viscosity(wPhaseIdx);
mobility_[nPhaseIdx] = MaterialLaw::krn(matParams, saturation(wPhaseIdx))/fluidState_.viscosity(nPhaseIdx);
mobility_[wPhaseIdx_] = MaterialLaw::krw(matParams, saturation(wPhaseIdx_))/fluidState_.viscosity(wPhaseIdx_);
mobility_[nPhaseIdx] = MaterialLaw::krn(matParams, saturation(wPhaseIdx_))/fluidState_.viscosity(nPhaseIdx);
//update porosity before calculating the effective properties depending on it
updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, numFluidComps);
// binary diffusion coefficients
for (unsigned int compJIdx = 0; compJIdx < ModelTraits::numFluidComponents(); ++compJIdx)
{
if(compJIdx != comp0Idx)
{
setDiffusionCoefficient_( phase0Idx, compJIdx,
FluidSystem::binaryDiffusionCoefficient(fluidState_,
paramCache,
phase0Idx,
comp0Idx,
compJIdx) );
setEffectiveDiffusionCoefficient_(phase0Idx, compJIdx);
}
if(compJIdx != comp1Idx)
{
setDiffusionCoefficient_( phase1Idx, compJIdx,
FluidSystem::binaryDiffusionCoefficient(fluidState_,
paramCache,
phase1Idx,
comp1Idx,
compJIdx) );
setEffectiveDiffusionCoefficient_(phase1Idx, compJIdx);
}
}
// calculate the remaining quantities
updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, numFluidComps);
EnergyVolVars::updateSolidEnergyParams(elemSol, problem, element, scv, solidState_);
permeability_ = problem.spatialParams().permeability(element, scv, elemSol);
EnergyVolVars::updateEffectiveThermalConductivity();
}
/*!
......@@ -204,8 +212,8 @@ public:
using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
const int wPhaseIdx = problem.spatialParams().template wettingPhase<FluidSystem>(element, scv, elemSol);
fluidState.setWettingPhase(wPhaseIdx);
wPhaseIdx_ = problem.spatialParams().template wettingPhase<FluidSystem>(element, scv, elemSol);
fluidState.setWettingPhase(wPhaseIdx_);
// set the saturations
if (phasePresence == secondPhaseOnly)
......@@ -235,17 +243,17 @@ public:
DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
// set pressures of the fluid phases
pc_ = MaterialLaw::pc(materialParams, fluidState.saturation(wPhaseIdx));
pc_ = MaterialLaw::pc(materialParams, fluidState.saturation(wPhaseIdx_));
if (formulation == TwoPFormulation::p0s1)
{
fluidState.setPressure(phase0Idx, priVars[pressureIdx]);
fluidState.setPressure(phase1Idx, (wPhaseIdx == phase0Idx) ? priVars[pressureIdx] + pc_
fluidState.setPressure(phase1Idx, (wPhaseIdx_ == phase0Idx) ? priVars[pressureIdx] + pc_
: priVars[pressureIdx] - pc_);
}
else
{
fluidState.setPressure(phase1Idx, priVars[pressureIdx]);
fluidState.setPressure(phase0Idx, (wPhaseIdx == phase0Idx) ? priVars[pressureIdx] - pc_
fluidState.setPressure(phase0Idx, (wPhaseIdx_ == phase0Idx) ? priVars[pressureIdx] - pc_
: priVars[pressureIdx] + pc_);
}
......@@ -463,6 +471,19 @@ public:
DUNE_THROW(Dune::InvalidStateException, "Diffusion coefficient called for phaseIdx = compIdx");
}
/*!
* \brief Returns the effective diffusion coefficients for a phase in \f$[m^2/s]\f$.
*/
Scalar effectiveDiffusivity(int phaseIdx, int compIdx) const
{
if (compIdx < phaseIdx)
return effectiveDiffCoeff_[phaseIdx][compIdx];
else if (compIdx > phaseIdx)
return effectiveDiffCoeff_[phaseIdx][compIdx-1];
else
DUNE_THROW(Dune::InvalidStateException, "Diffusion coefficient called for phaseIdx = compIdx");
}
/*!
* \brief Returns the mass fraction of a component in the phase
*
......@@ -481,6 +502,12 @@ public:
Scalar moleFraction(int phaseIdx, int compIdx) const
{ return fluidState_.moleFraction(phaseIdx, compIdx); }
/*!
* \brief Returns the wetting phase index
*/
const int wettingPhaseIdx() const
{ return wPhaseIdx_; }
protected:
FluidState fluidState_;
SolidState solidState_;
......@@ -496,11 +523,23 @@ private:
DUNE_THROW(Dune::InvalidStateException, "Diffusion coefficient for phaseIdx = compIdx doesn't exist");
}
void setEffectiveDiffusionCoefficient_(int phaseIdx, int compIdx)
{
if (compIdx < phaseIdx)
effectiveDiffCoeff_[phaseIdx][compIdx] = EffDiffModel::effectiveDiffusivity(*this, diffCoefficient_[phaseIdx][compIdx], phaseIdx);
else if (compIdx > phaseIdx)
effectiveDiffCoeff_[phaseIdx][compIdx-1] = EffDiffModel::effectiveDiffusivity(*this, diffCoefficient_[phaseIdx][compIdx-1], phaseIdx);
else
DUNE_THROW(Dune::InvalidStateException, "Effective diffusion coefficient for phaseIdx = compIdx doesn't exist");
}
Scalar pc_; // The capillary pressure
Scalar porosity_; // Effective porosity within the control volume
PermeabilityType permeability_; // Effective permeability within the control volume
Scalar mobility_[ModelTraits::numFluidPhases()]; // Effective mobility within the control volume
std::array<std::array<Scalar, ModelTraits::numFluidComponents()-1>, ModelTraits::numFluidPhases()> diffCoefficient_;
std::array<std::array<Scalar, ModelTraits::numFluidComponents()-1>, ModelTraits::numFluidPhases()> effectiveDiffCoeff_;
int wPhaseIdx_;
};
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment