From 92d7f5a17919f5a3305ce34e72cd5ab033717510 Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Fri, 17 Feb 2017 18:36:59 +0100
Subject: [PATCH] [staggeredGrid][NavierStokesNC] Improve model
---
dumux/freeflow/staggered/localresidual.hh | 134 +-----------------
dumux/freeflow/staggerednc/fluxvariables.hh | 11 +-
dumux/freeflow/staggerednc/indices.hh | 16 ++-
dumux/freeflow/staggerednc/localresidual.hh | 94 +++---------
dumux/freeflow/staggerednc/model.hh | 14 +-
dumux/freeflow/staggerednc/properties.hh | 1 +
.../freeflow/staggerednc/propertydefaults.hh | 78 ++--------
dumux/freeflow/staggerednc/volumevariables.hh | 80 ++++++-----
dumux/freeflow/stokesnc/localresidual.hh | 61 +++++---
.../staggerednc/channeltestproblem.hh | 51 ++++---
10 files changed, 186 insertions(+), 354 deletions(-)
diff --git a/dumux/freeflow/staggered/localresidual.hh b/dumux/freeflow/staggered/localresidual.hh
index 7764160d5f..e6f655970f 100644
--- a/dumux/freeflow/staggered/localresidual.hh
+++ b/dumux/freeflow/staggered/localresidual.hh
@@ -160,7 +160,7 @@ public:
const VolumeVariables& volVars)
{
CellCenterPrimaryVariables storage;
- storage[0] = volVars.density(0);
+ storage[0] = volVars.density();
return storage;
}
@@ -180,7 +180,7 @@ public:
{
FacePrimaryVariables storage(0.0);
const Scalar velocity = globalFaceVars.faceVars(scvf.dofIndex()).velocity();
- storage[0] = volVars.density(0) * velocity;
+ storage[0] = volVars.density() * velocity;
return storage;
}
@@ -345,132 +345,6 @@ private:
}
};
-// // specialization for miscible, isothermal flow
-// template<class TypeTag>
-// class StaggeredNavierStokesResidualImpl<TypeTag, true, false> : public StaggeredNavierStokesResidualImpl<TypeTag, false, false>
-// {
-// using ParentType = StaggeredNavierStokesResidualImpl<TypeTag, false, false>;
-// friend class StaggeredLocalResidual<TypeTag>;
-// using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
-//
-// enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
-//
-// using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-// using Implementation = typename GET_PROP_TYPE(TypeTag, LocalResidual);
-// using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
-// using Element = typename GridView::template Codim<0>::Entity;
-// using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
-// using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
-// using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
-// using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
-// using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
-// using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
-// using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
-// using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
-// using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
-// using CellCenterSolutionVector = typename GET_PROP_TYPE(TypeTag, CellCenterSolutionVector);
-// using FaceSolutionVector = typename GET_PROP_TYPE(TypeTag, FaceSolutionVector);
-// using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
-// using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
-// using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
-// using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
-//
-//
-// using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
-// typename DofTypeIndices::CellCenterIdx cellCenterIdx;
-// typename DofTypeIndices::FaceIdx faceIdx;
-//
-// enum {
-// // grid and world dimension
-// dim = GridView::dimension,
-// dimWorld = GridView::dimensionworld,
-//
-// pressureIdx = Indices::pressureIdx,
-// velocityIdx = Indices::velocityIdx,
-//
-// massBalanceIdx = Indices::massBalanceIdx,
-// momentumBalanceIdx = Indices::momentumBalanceIdx,
-//
-// conti0EqIdx = Indices::conti0EqIdx
-// };
-//
-// using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
-// using GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
-//
-// static constexpr bool navierStokes = GET_PROP_VALUE(TypeTag, EnableInertiaTerms);
-// static constexpr int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
-//
-// static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
-//
-// //! The index of the component balance equation that gets replaced with the total mass balance
-// static constexpr int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
-// /*!
-// * \brief Evaluate the rate of change of all conservation
-// * quantites (e.g. phase mass) within a sub-control
-// * volume of a finite volume element for the immiscible models.
-// * \param scv The sub control volume
-// * \param volVars The current or previous volVars
-// * \note This function should not include the source and sink terms.
-// * \note The volVars can be different to allow computing
-// * the implicit euler time derivative here
-// */
-// CellCenterPrimaryVariables computeStorageForCellCenter(const SubControlVolume& scv,
-// const VolumeVariables& volVars,
-// bool useMoles = true) const
-// {
-// PrimaryVariables storage(0.0);
-//
-// // formulation with mole balances
-// if (useMoles)
-// {
-// // compute storage term of all components
-// for (int compIdx = 0; compIdx < numComponents; ++compIdx)
-// {
-// auto eqIdx = conti0EqIdx + compIdx;
-// auto s = volVars.molarDensity(0)
-// * volVars.moleFraction(0, compIdx);
-//
-// if (eqIdx != replaceCompEqIdx)
-// storage[eqIdx] += s;
-//
-// // in case one balance is substituted by the total mass balance
-// if (replaceCompEqIdx < numComponents)
-// storage[replaceCompEqIdx] += s;
-// }
-//
-// //! The energy storage in the fluid phase with index phaseIdx
-// // EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, phaseIdx);
-//
-// //! The energy storage in the solid matrix
-// // EnergyLocalResidual::solidPhaseStorage(storage, scv, volVars);
-// }
-// // formulation with mass balances
-// else
-// {
-// for (int compIdx = 0; compIdx < numComponents; ++compIdx)
-// {
-// auto eqIdx = conti0EqIdx + compIdx;
-// auto s = volVars.density(0)
-// * volVars.massFraction(0, compIdx);
-//
-// if (eqIdx != replaceCompEqIdx)
-// storage[eqIdx] += s;
-//
-// // in case one balance is substituted by the total mass balance
-// if (replaceCompEqIdx < numComponents)
-// storage[replaceCompEqIdx] += s;
-// }
-//
-// //! The energy storage in the fluid phase with index phaseIdx
-// // EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, phaseIdx);
-//
-// //! The energy storage in the solid matrix
-// // EnergyLocalResidual::solidPhaseStorage(storage, scv, volVars);
-// }
-//
-// return storage;
-// }
-// };
// specialization for immiscible, non-isothermal flow
template<class TypeTag>
@@ -531,8 +405,8 @@ class StaggeredNavierStokesResidualImpl<TypeTag, false, true> : public Staggered
const VolumeVariables& volVars)
{
CellCenterPrimaryVariables storage;
- storage[massBalanceIdx] = volVars.density(0);
- storage[energyBalanceIdx] = volVars.density(0) * volVars.internalEnergy(0);
+ storage[massBalanceIdx] = volVars.density();
+ storage[energyBalanceIdx] = volVars.density() * volVars.internalEnergy();
return storage;
}
};
diff --git a/dumux/freeflow/staggerednc/fluxvariables.hh b/dumux/freeflow/staggerednc/fluxvariables.hh
index 44598f6c60..da1f309a1b 100644
--- a/dumux/freeflow/staggerednc/fluxvariables.hh
+++ b/dumux/freeflow/staggerednc/fluxvariables.hh
@@ -93,7 +93,8 @@ class FreeFlowFluxVariablesImpl<TypeTag, true, false> : public FreeFlowFluxVaria
massBalanceIdx = Indices::massBalanceIdx,
momentumBalanceIdx = Indices::momentumBalanceIdx,
- conti0EqIdx = Indices::conti0EqIdx
+ conti0EqIdx = Indices::conti0EqIdx,
+ phaseIdx = Indices::phaseIdx
};
public:
@@ -130,11 +131,13 @@ private:
const Scalar velocity = globalFaceVars.faceVars(scvf.dofIndex()).velocity();
- const bool insideIsUpstream = sign(scvf.outerNormalScalar()) == sign(velocity) ? true : false;
+ const bool insideIsUpstream = sign(scvf.outerNormalScalar()) == sign(velocity);
const auto& upstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars;
const auto& downstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars;
const Scalar upWindWeight = GET_PROP_VALUE(TypeTag, ImplicitUpwindWeight);
+ const Scalar upstreamDensity = useMoles ? upstreamVolVars.molarDensity() : upstreamVolVars.density();
+ const Scalar downstreamDensity = useMoles ? downstreamVolVars.molarDensity() : downstreamVolVars.density();
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
{
@@ -143,8 +146,8 @@ private:
const Scalar upstreamDensity = useMoles ? upstreamVolVars.molarDensity() : upstreamVolVars.density();
const Scalar downstreamDensity = useMoles ? downstreamVolVars.molarDensity() : downstreamVolVars.density();
- const Scalar upstreamFraction = useMoles ? upstreamVolVars.moleFraction(0, compIdx) : upstreamVolVars.massFraction(0, compIdx);
- const Scalar downstreamFraction = useMoles ? downstreamVolVars.moleFraction(0, compIdx) : downstreamVolVars.massFraction(0, compIdx);
+ const Scalar upstreamFraction = useMoles ? upstreamVolVars.moleFraction(phaseIdx, compIdx) : upstreamVolVars.massFraction(phaseIdx, compIdx);
+ const Scalar downstreamFraction = useMoles ? downstreamVolVars.moleFraction(phaseIdx, compIdx) : downstreamVolVars.massFraction(phaseIdx, compIdx);
Scalar advFlux = 0.0;
diff --git a/dumux/freeflow/staggerednc/indices.hh b/dumux/freeflow/staggerednc/indices.hh
index 3d2383f480..b68c4d4ea0 100644
--- a/dumux/freeflow/staggerednc/indices.hh
+++ b/dumux/freeflow/staggerednc/indices.hh
@@ -18,7 +18,7 @@
*****************************************************************************/
/*!
* \file
- * \brief Defines the indices for the one-phase fully implicit model.
+ * \brief Defines the indices for the staggered Navier-Stokes NC model.
*/
#ifndef DUMUX_STAGGERED_NAVIERSTOKES_NC_INDICES_HH
#define DUMUX_STAGGERED_NAVIERSTOKES_NC_INDICES_HH
@@ -29,15 +29,25 @@ namespace Dumux
{
// \{
/*!
- * \ingroup NavierStokesModel
+ * \ingroup NavierStokesNCModel
* \ingroup ImplicitIndices
- * \brief The common indices for the isothermal stokes model.
+ * \brief Indices for the staggered Navier-Stokes NC model model.
*
* \tparam PVOffset The first index in a primary variable vector.
*/
template <class TypeTag, int PVOffset = 0>
struct NavierStokesNCIndices : public NavierStokesCommonIndices<TypeTag, PVOffset>
{
+private:
+ using ParentType = NavierStokesCommonIndices<TypeTag, PVOffset>;
+
+public:
+
+ static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
+ static constexpr int mainCompIdx = phaseIdx;
+ static constexpr int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
+ // TODO: componentIdx?
+ //TODO: what about the offset?
};
diff --git a/dumux/freeflow/staggerednc/localresidual.hh b/dumux/freeflow/staggerednc/localresidual.hh
index 9496f9ca7a..b96ff69a91 100644
--- a/dumux/freeflow/staggerednc/localresidual.hh
+++ b/dumux/freeflow/staggerednc/localresidual.hh
@@ -60,29 +60,11 @@ class StaggeredNavierStokesResidualImpl;
template<class TypeTag>
class StaggeredNavierStokesResidualImpl<TypeTag, true, false> : public StaggeredNavierStokesResidualImpl<TypeTag, false, false>
{
- using ParentType = StaggeredNavierStokesResidualImpl<TypeTag, false, false>;
friend class StaggeredLocalResidual<TypeTag>;
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
-
- enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using Implementation = typename GET_PROP_TYPE(TypeTag, LocalResidual);
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using Element = typename GridView::template Codim<0>::Entity;
- using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
- using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
- using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
- using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
- using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
- using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
- using CellCenterSolutionVector = typename GET_PROP_TYPE(TypeTag, CellCenterSolutionVector);
- using FaceSolutionVector = typename GET_PROP_TYPE(TypeTag, FaceSolutionVector);
using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
- using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
@@ -92,31 +74,21 @@ class StaggeredNavierStokesResidualImpl<TypeTag, true, false> : public Staggered
typename DofTypeIndices::FaceIdx faceIdx;
enum {
- // grid and world dimension
- dim = GridView::dimension,
- dimWorld = GridView::dimensionworld,
-
- pressureIdx = Indices::pressureIdx,
- velocityIdx = Indices::velocityIdx,
-
- massBalanceIdx = Indices::massBalanceIdx,
- momentumBalanceIdx = Indices::momentumBalanceIdx,
+ conti0EqIdx = Indices::conti0EqIdx,
+ phaseIdx = Indices::phaseIdx,
- conti0EqIdx = Indices::conti0EqIdx
+ // The index of the component balance equation
+ // that gets replaced with the total mass balance
+ replaceCompEqIdx = Indices::replaceCompEqIdx
};
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
- static constexpr bool navierStokes = GET_PROP_VALUE(TypeTag, EnableInertiaTerms);
static constexpr int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
-
static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
- static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
- //! The index of the component balance equation that gets replaced with the total mass balance
- static constexpr int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
/*!
* \brief Evaluate the rate of change of all conservation
* quantites (e.g. phase mass) within a sub-control
@@ -132,54 +104,24 @@ class StaggeredNavierStokesResidualImpl<TypeTag, true, false> : public Staggered
{
CellCenterPrimaryVariables storage(0.0);
- // formulation with mole balances
- if (useMoles)
- {
- // compute storage term of all components
- for (int compIdx = 0; compIdx < numComponents; ++compIdx)
- {
- auto eqIdx = conti0EqIdx + compIdx;
- auto s = volVars.molarDensity(0)
- * volVars.moleFraction(0, compIdx);
-
- if (eqIdx != replaceCompEqIdx)
- storage[eqIdx] += s;
-
- // in case one balance is substituted by the total mass balance
- if (replaceCompEqIdx < numComponents)
- storage[replaceCompEqIdx] += s;
- }
-
- //! The energy storage in the fluid phase with index phaseIdx
-// EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, phaseIdx);
-
- //! The energy storage in the solid matrix
-// EnergyLocalResidual::solidPhaseStorage(storage, scv, volVars);
- }
- // formulation with mass balances
- else
- {
- for (int compIdx = 0; compIdx < numComponents; ++compIdx)
- {
- auto eqIdx = conti0EqIdx + compIdx;
- auto s = volVars.density(0)
- * volVars.massFraction(0, compIdx);
+ const Scalar density = useMoles ? volVars.molarDensity() : volVars.density();
- if (eqIdx != replaceCompEqIdx)
- storage[eqIdx] += s;
-
- // in case one balance is substituted by the total mass balance
- if (replaceCompEqIdx < numComponents)
- storage[replaceCompEqIdx] += s;
- }
+ // compute storage term of all components
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ {
+ const int eqIdx = conti0EqIdx + compIdx;
- //! The energy storage in the fluid phase with index phaseIdx
-// EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, phaseIdx);
+ const Scalar massOrMoleFraction = useMoles ? volVars.moleFraction(phaseIdx, compIdx) : volVars.massFraction(phaseIdx, compIdx);
+ const Scalar s = density * massOrMoleFraction;
- //! The energy storage in the solid matrix
-// EnergyLocalResidual::solidPhaseStorage(storage, scv, volVars);
+ if (eqIdx != replaceCompEqIdx)
+ storage[eqIdx] += s;
}
+ // in case one balance is substituted by the total mass balance
+ if(replaceCompEqIdx < numComponents)
+ storage[replaceCompEqIdx] += density;
+ //TODO: energy balance
return storage;
}
};
diff --git a/dumux/freeflow/staggerednc/model.hh b/dumux/freeflow/staggerednc/model.hh
index 01198e50f7..040319ebe8 100644
--- a/dumux/freeflow/staggerednc/model.hh
+++ b/dumux/freeflow/staggerednc/model.hh
@@ -85,6 +85,8 @@ class NavierStokesNCModel : public NavierStokesModel<TypeTag>
typename DofTypeIndices::CellCenterIdx cellCenterIdx;
typename DofTypeIndices::FaceIdx faceIdx;
+ enum { phaseIdx = Indices::phaseIdx };
+
public:
void init(Problem& problem)
@@ -93,10 +95,16 @@ public:
// register standardized vtk output fields
auto& vtkOutputModule = problem.vtkOutputModule();
- vtkOutputModule.addSecondaryVariable("rho",[](const VolumeVariables& v){ return v.molarDensity(); });
+ vtkOutputModule.addSecondaryVariable("rhoMolar",[](const VolumeVariables& v){ return v.molarDensity(); });
+ vtkOutputModule.addSecondaryVariable("rho",[](const VolumeVariables& v){ return v.density(); });
for (int j = 0; j < numComponents; ++j)
- vtkOutputModule.addSecondaryVariable("x" + FluidSystem::componentName(j) + FluidSystem::phaseName(0),
- [j](const VolumeVariables& v){ return v.massFraction(0,j); });
+ {
+ vtkOutputModule.addSecondaryVariable("X^" + FluidSystem::componentName(j) + "_" + FluidSystem::phaseName(phaseIdx),
+ [j](const VolumeVariables& v){ return v.massFraction(phaseIdx,j); });
+
+ vtkOutputModule.addSecondaryVariable("x^" + FluidSystem::componentName(j) + "_" + FluidSystem::phaseName(phaseIdx),
+ [j](const VolumeVariables& v){ return v.moleFraction(phaseIdx,j); });
+ }
// NonIsothermalModel::maybeAddTemperature(vtkOutputModule);
}
diff --git a/dumux/freeflow/staggerednc/properties.hh b/dumux/freeflow/staggerednc/properties.hh
index ef5c28a23f..b50ad8f482 100644
--- a/dumux/freeflow/staggerednc/properties.hh
+++ b/dumux/freeflow/staggerednc/properties.hh
@@ -68,6 +68,7 @@ NEW_PROP_TAG(EnableEnergyTransport); //!< Returns whether to consider energy tr
NEW_PROP_TAG(FaceVariables); //!< Returns whether to consider energy transport or not
NEW_PROP_TAG(ReplaceCompEqIdx); //!< Returns whether to consider energy transport or not
NEW_PROP_TAG(UseMoles); //!< Defines whether molar (true) or mass (false) density is used
+NEW_PROP_TAG(PhaseIdx); //!< Defines the phaseIdx
// \}
}
diff --git a/dumux/freeflow/staggerednc/propertydefaults.hh b/dumux/freeflow/staggerednc/propertydefaults.hh
index 31d7d5b87a..329c732640 100644
--- a/dumux/freeflow/staggerednc/propertydefaults.hh
+++ b/dumux/freeflow/staggerednc/propertydefaults.hh
@@ -63,12 +63,16 @@ NEW_PROP_TAG(FluxVariablesCache);
// default property values for the isothermal single phase model
///////////////////////////////////////////////////////////////////////////
namespace Properties {
-// SET_INT_PROP(NavierStokes, NumEqCellCenter, 1); //!< set the number of equations to 1
-// SET_INT_PROP(NavierStokes, NumEqFace, 1); //!< set the number of equations to 1
-// SET_INT_PROP(NavierStokes, NumPhases, 1); //!< The number of phases in the 1p model is 1
-//
-//
-SET_INT_PROP(NavierStokesNC, NumEqCellCenter, 2);
+
+SET_PROP(NavierStokesNC, NumEqCellCenter)
+{
+private:
+ static constexpr int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
+public:
+ static constexpr int value = numComponents;
+};
+
+
SET_INT_PROP(NavierStokesNC, ReplaceCompEqIdx, 0);
@@ -92,6 +96,7 @@ public:
//! the VolumeVariables property
SET_TYPE_PROP(NavierStokesNC, VolumeVariables, NavierStokesNCVolumeVariables<TypeTag>);
SET_TYPE_PROP(NavierStokesNC, Model, NavierStokesNCModel<TypeTag>);
+SET_TYPE_PROP(NavierStokesNC, Indices, NavierStokesNCIndices<TypeTag>);
/*!
@@ -109,72 +114,13 @@ SET_PROP(NavierStokesNC, FluidState)
typedef CompositionalFluidState<Scalar, FluidSystem> type;
};
-// //! Enable advection
-// SET_BOOL_PROP(NavierStokes, EnableAdvection, true);
-//
-// //! The one-phase model has no molecular diffusion
-// SET_BOOL_PROP(NavierStokes, EnableMolecularDiffusion, false);
-//
-// //! Isothermal model by default
-// SET_BOOL_PROP(NavierStokes, EnableEnergyBalance, false);
-//
-// //! The indices required by the isothermal single-phase model
-// SET_TYPE_PROP(NavierStokes, Indices, NavierStokesCommonIndices<TypeTag>);
-//
-// //! The weight of the upwind control volume when calculating
-// //! fluxes. Use central differences by default.
-// SET_SCALAR_PROP(NavierStokes, ImplicitMassUpwindWeight, 0.5);
-//
-// //! weight for the upwind mobility in the velocity calculation
-// //! fluxes. Use central differences by default.
-// SET_SCALAR_PROP(NavierStokes, ImplicitMobilityUpwindWeight, 0.5);
-//
-// //! The fluid system to use by default
-// SET_TYPE_PROP(NavierStokes, FluidSystem, Dumux::FluidSystems::OneP<typename GET_PROP_TYPE(TypeTag, Scalar), typename GET_PROP_TYPE(TypeTag, Fluid)>);
-//
-// SET_PROP(NavierStokes, Fluid)
-// { private:
-// typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-// public:
-// typedef FluidSystems::LiquidPhase<Scalar, Dumux::NullComponent<Scalar> > type;
-// };
-//
-// /*!
-// * \brief The fluid state which is used by the volume variables to
-// * store the thermodynamic state. This should be chosen
-// * appropriately for the model ((non-)isothermal, equilibrium, ...).
-// * This can be done in the problem.
-// */
-// SET_PROP(NavierStokes, FluidState){
-// private:
-// typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-// typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
-// public:
-// typedef Dumux::ImmiscibleFluidState<Scalar, FluidSystem> type;
-// };
-//
-// // disable velocity output by default
-// SET_BOOL_PROP(NavierStokes, VtkAddVelocity, true);
-//
-// // enable gravity by default
-// SET_BOOL_PROP(NavierStokes, ProblemEnableGravity, true);
-//
-// SET_BOOL_PROP(NavierStokes, EnableInertiaTerms, true);
-//
-// SET_BOOL_PROP(NavierStokes, EnableEnergyTransport, false);
//
SET_BOOL_PROP(NavierStokesNC, EnableComponentTransport, true);
SET_BOOL_PROP(NavierStokesNC, UseMoles, false); //!< Defines whether molar (true) or mass (false) density is used
+SET_INT_PROP(NavierStokesNC, PhaseIdx, 0); //!< Defines the phaseIdx
-//! average is used as default model to compute the effective thermal heat conductivity
-// SET_PROP(NavierStokesNI, ThermalConductivityModel)
-// { private :
-// typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-// public:
-// typedef ThermalConductivityAverage<Scalar> type;
-// };
//////////////////////////////////////////////////////////////////
// Property values for isothermal model required for the general non-isothermal model
diff --git a/dumux/freeflow/staggerednc/volumevariables.hh b/dumux/freeflow/staggerednc/volumevariables.hh
index 9d8515285c..673cce0e4e 100644
--- a/dumux/freeflow/staggerednc/volumevariables.hh
+++ b/dumux/freeflow/staggerednc/volumevariables.hh
@@ -54,7 +54,14 @@ class NavierStokesNCVolumeVariables : public ImplicitVolumeVariables<TypeTag>
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Element = typename GridView::template Codim<0>::Entity;
- enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents)};
+ enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents),
+ numPhases = FluidSystem::numPhases,
+ phaseIdx = Indices::phaseIdx,
+ mainCompIdx = Indices::mainCompIdx,
+ pressureIdx = Indices::pressureIdx
+ };
+
+ static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
public:
@@ -75,25 +82,20 @@ public:
typename FluidSystem::ParameterCache paramCache;
paramCache.updateAll(fluidState_);
- int compIIdx = 0;
+ int compIIdx = phaseIdx;
for (unsigned int compJIdx = 0; compJIdx < numComponents; ++compJIdx)
{
// binary diffusion coefficents
if(compIIdx!= compJIdx)
{
- setDiffusionCoefficient_(0, compJIdx,
+ setDiffusionCoefficient_(phaseIdx, compJIdx,
FluidSystem::binaryDiffusionCoefficient(fluidState_,
paramCache,
- 0,
+ phaseIdx,
compIIdx,
compJIdx));
}
}
-
-// std::cout << "mole frac 0 is : " << moleFraction(0,0) << std::endl;
-// std::cout << "mole frac 1 is : " << moleFraction(0,1) << std::endl;
- if(moleFraction(0,0) + moleFraction(0,1) > 1.000001)
- std::cout << "sum: " << moleFraction(0,0) + moleFraction(0,1) << std::endl;
};
/*!
@@ -107,32 +109,46 @@ public:
{
Scalar t = ParentType::temperature(elemSol, problem, element, scv);
fluidState.setTemperature(t);
- fluidState.setSaturation(/*phaseIdx=*/0, 1.);
-
- fluidState.setPressure(/*phaseIdx=*/0, elemSol[0][Indices::pressureIdx]);
-
- fluidState.setMassFraction(0, 1, elemSol[0][1]);
- fluidState.setMassFraction(0, 0, 1.0 - elemSol[0][1]);
+ fluidState.setPressure(phaseIdx, elemSol[0][Indices::pressureIdx]);
+ Scalar fracMinor = 0.0;
+ int transportEqIdx = 1;
+
+ for(int compIdx = 0; compIdx < numComponents; ++compIdx)
+ {
+ if(compIdx == mainCompIdx)
+ continue;
+
+ const Scalar moleOrMassFraction = elemSol[0][transportEqIdx++] + 1.0;
+ if(useMoles)
+ fluidState.setMoleFraction(phaseIdx, compIdx, moleOrMassFraction -1.0);
+ else
+ fluidState.setMassFraction(phaseIdx, compIdx, moleOrMassFraction -1.0);
+ fracMinor += moleOrMassFraction - 1.0;
+ }
+ if(useMoles)
+ fluidState.setMoleFraction(phaseIdx, mainCompIdx, 1.0 - fracMinor);
+ else
+ fluidState.setMassFraction(phaseIdx, mainCompIdx, 1.0 - fracMinor);
// saturation in a single phase is always 1 and thus redundant
// to set. But since we use the fluid state shared by the
// immiscible multi-phase models, so we have to set it here...
- fluidState.setSaturation(/*phaseIdx=*/0, 1.0);
+ fluidState.setSaturation(phaseIdx, 1.0);
typename FluidSystem::ParameterCache paramCache;
- paramCache.updatePhase(fluidState, /*phaseIdx=*/0);
+ paramCache.updatePhase(fluidState, phaseIdx);
- Scalar value = FluidSystem::density(fluidState, paramCache, /*phaseIdx=*/0);
- fluidState.setDensity(/*phaseIdx=*/0, value);
+ Scalar value = FluidSystem::density(fluidState, paramCache, phaseIdx);
+ fluidState.setDensity(phaseIdx, value);
- value = FluidSystem::viscosity(fluidState, paramCache, /*phaseIdx=*/0);
- fluidState.setViscosity(/*phaseIdx=*/0, value);
+ value = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
+ fluidState.setViscosity(phaseIdx, value);
// compute and set the enthalpy
- value = ParentType::enthalpy(fluidState, paramCache, /*phaseIdx=*/0);
- fluidState.setEnthalpy(/*phaseIdx=*/0, value);
+ value = ParentType::enthalpy(fluidState, paramCache, phaseIdx);
+ fluidState.setEnthalpy(phaseIdx, value);
}
/*!
@@ -149,27 +165,27 @@ public:
* \brief Return the effective pressure \f$\mathrm{[Pa]}\f$ of a given phase within
* the control volume.
*/
- Scalar pressure(int phaseIdx = 0) const
+ Scalar pressure() const
{ return fluidState_.pressure(phaseIdx); }
/*!
* \brief Return the saturation
*/
- Scalar saturation(int phaseIdx = 0) const
+ Scalar saturation() const
{ return 1.0; }
/*!
* \brief Return the mass density \f$\mathrm{[kg/m^3]}\f$ of a given phase within the
* control volume.
*/
- Scalar density(int phaseIdx = 0) const
+ Scalar density() const
{ return fluidState_.density(phaseIdx); }
/*!
* \brief Return the dynamic viscosity \f$\mathrm{[Pa s]}\f$ of the fluid within the
* control volume.
*/
- Scalar viscosity(int phaseIdx = 0) const
+ Scalar viscosity() const
{ return fluidState_.viscosity(phaseIdx); }
/*!
@@ -196,13 +212,9 @@ public:
*
* \param phaseIdx The phase index
*/
- Scalar molarDensity(int phaseIdx = 0) const
+ Scalar molarDensity() const
{
- if (phaseIdx < 1/*numPhases*/)
- return fluidState_.molarDensity(phaseIdx);
-
- else
- DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
+ return fluidState_.molarDensity(phaseIdx);
}
/*!
@@ -244,7 +256,7 @@ protected:
}
- std::array<std::array<Scalar, numComponents-1>, 1> diffCoefficient_;
+ std::array<std::array<Scalar, numComponents-1>, numPhases> diffCoefficient_;
};
}
diff --git a/dumux/freeflow/stokesnc/localresidual.hh b/dumux/freeflow/stokesnc/localresidual.hh
index b23474c4a3..86bda5cfa6 100644
--- a/dumux/freeflow/stokesnc/localresidual.hh
+++ b/dumux/freeflow/stokesnc/localresidual.hh
@@ -111,36 +111,55 @@ public:
this->prevVolVars_() : this->curVolVars_();
const VolumeVariables &volVars = elemVolVars[scvIdx];
- if (useMoles)
+ const Scalar density = useMoles ? volVars.molarDensity() : volVars.density();
+ // mass balance and transport equations
+ for (int compIdx=0; compIdx<numComponents; compIdx++)
{
- // mass balance and transport equations
- for (int compIdx=0; compIdx<numComponents; compIdx++)
- {
- if (conti0EqIdx+compIdx != massBalanceIdx)
+ const Scalar moleOrMassFraction = useMoles ? volVars.moleFraction(compIdx) : volVars.massFraction(transportCompIdx);
+
+ if (conti0EqIdx+compIdx != massBalanceIdx)
//else // transport equations
{
- storage[conti0EqIdx+compIdx] = volVars.molarDensity()
- * volVars.moleFraction(compIdx);
+ storage[conti0EqIdx+compIdx] = density * moleOrMassFraction;
- Valgrind::CheckDefined(volVars.molarDensity());
- Valgrind::CheckDefined(volVars.moleFraction(compIdx));
+ Valgrind::CheckDefined(density);
+ Valgrind::CheckDefined(moleOrMassFraction);
}
}
- }
- else
- {
- /* works for a maximum of two components, for more components
- mole fractions must be used (set property UseMoles to true) */
- //storage of transported component
- storage[transportEqIdx] = volVars.density()
- * volVars.massFraction(transportCompIdx);
+ }
- Valgrind::CheckDefined(volVars.density());
- Valgrind::CheckDefined(volVars.massFraction(transportCompIdx));
- }
- }
+// if (useMoles)
+// {
+// // mass balance and transport equations
+// for (int compIdx=0; compIdx<numComponents; compIdx++)
+// {
+// if (conti0EqIdx+compIdx != massBalanceIdx)
+// //else // transport equations
+// {
+// storage[conti0EqIdx+compIdx] = volVars.molarDensity()
+// * volVars.moleFraction(compIdx);
+//
+// Valgrind::CheckDefined(volVars.molarDensity());
+// Valgrind::CheckDefined(volVars.moleFraction(compIdx));
+// }
+// }
+// }
+// else
+// {
+// /* works for a maximum of two components, for more components
+// mole fractions must be used (set property UseMoles to true) */
+//
+// //storage of transported component
+// storage[transportEqIdx] = volVars.density()
+// * volVars.massFraction(transportCompIdx);
+//
+// Valgrind::CheckDefined(volVars.density());
+// Valgrind::CheckDefined(volVars.massFraction(transportCompIdx));
+//
+// }
+// }
/*!
* \brief Evaluates the advective component (mass) flux
diff --git a/test/freeflow/staggerednc/channeltestproblem.hh b/test/freeflow/staggerednc/channeltestproblem.hh
index cf22592c81..ec0f0b7162 100644
--- a/test/freeflow/staggerednc/channeltestproblem.hh
+++ b/test/freeflow/staggerednc/channeltestproblem.hh
@@ -49,17 +49,21 @@ namespace Properties
{
NEW_TYPE_TAG(ChannelTestProblem, INHERITS_FROM(StaggeredModel, NavierStokesNC));
-// SET_PROP(ChannelTestProblem, Fluid)
-// {
-// private:
-// typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-// public:
-// typedef FluidSystems::LiquidPhase<Scalar, Dumux::Constant<TypeTag, Scalar> > type;
-// };
+NEW_PROP_TAG(FluidSystem);
+
// Select the fluid system
SET_TYPE_PROP(ChannelTestProblem, FluidSystem,
- FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>);
+ FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)/*, SimpleH2O<typename GET_PROP_TYPE(TypeTag, Scalar)>, false*/>);
+SET_PROP(ChannelTestProblem, PhaseIdx)
+{
+private:
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+public:
+ static constexpr int value = 0;//FluidSystem::wPhaseIdx;
+};
+
+SET_INT_PROP(ChannelTestProblem, ReplaceCompEqIdx, 0);
// Set the grid type
SET_TYPE_PROP(ChannelTestProblem, Grid, Dune::YaspGrid<2>);
@@ -74,6 +78,7 @@ SET_BOOL_PROP(ChannelTestProblem, EnableGlobalVolumeVariablesCache, true);
// Enable gravity
SET_BOOL_PROP(ChannelTestProblem, ProblemEnableGravity, true);
+SET_BOOL_PROP(ChannelTestProblem, UseMoles, true);
#if ENABLE_NAVIERSTOKES
SET_BOOL_PROP(ChannelTestProblem, EnableInertiaTerms, true);
@@ -93,6 +98,7 @@ class ChannelTestProblem : public NavierStokesProblem<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
// copy some indices for convenience
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
@@ -103,12 +109,14 @@ class ChannelTestProblem : public NavierStokesProblem<TypeTag>
};
enum {
massBalanceIdx = Indices::massBalanceIdx,
+ transportEqIdx = 1,
momentumBalanceIdx = Indices::momentumBalanceIdx,
momentumXBalanceIdx = Indices::momentumXBalanceIdx,
momentumYBalanceIdx = Indices::momentumYBalanceIdx,
pressureIdx = Indices::pressureIdx,
velocityXIdx = Indices::velocityXIdx,
- velocityYIdx = Indices::velocityYIdx
+ velocityYIdx = Indices::velocityYIdx,
+ transportCompIdx = 1/*FluidSystem::wCompIdx*/
};
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
@@ -122,6 +130,7 @@ class ChannelTestProblem : public NavierStokesProblem<TypeTag>
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
@@ -142,6 +151,13 @@ public:
Scalar,
Problem,
InletVelocity);
+
+ std::cout << "Indices: " << std::endl;
+ std::cout << "massBalanceIdx: " << massBalanceIdx << std::endl;
+ std::cout << "pressureIdx: " << pressureIdx << std::endl;
+ std::cout << "transportCompIdx: " << transportCompIdx << std::endl;
+
+ FluidSystem::init();
}
/*!
@@ -199,21 +215,22 @@ public:
// set Dirichlet values for the velocity everywhere
values.setDirichlet(momentumBalanceIdx);
- values.setNeumann(1);
+ values.setNeumann(transportEqIdx);
+ values.setNeumann(massBalanceIdx);
// values.setDirichlet(1);
- // set a fixed pressure in one cell
+// set a fixed pressure in one cell
if (isOutlet(globalPos))
{
values.setDirichlet(massBalanceIdx);
- values.setOutflow(momentumBalanceIdx);
- values.setOutflow(1);
+// values.setOutflow(momentumBalanceIdx);
+// values.setOutflow(transportEqIdx);
}
- else
- values.setNeumann(massBalanceIdx);
+// // else
+// values.setNeumann(massBalanceIdx);
if(isInlet(globalPos))
- values.setDirichlet(1);
+ values.setDirichlet(transportEqIdx);
return values;
}
@@ -234,7 +251,7 @@ public:
{
values[velocityXIdx] = inletVelocity_;
values[velocityYIdx] = 0.0;
- values[pressureIdx+1] =1e-6;
+ values[transportCompIdx] =1e-3;
}
else if(isWall(globalPos))
{
--
GitLab