Commit c2b86a8f authored by Timo Koch's avatar Timo Koch
Browse files

Actual stokesnc files replace dummies from last post. Now the...

Actual stokesnc files replace dummies from last post. Now the freeflow/stokesnc model is fully functioning. Reviewed by katherina and klaus.

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@12078 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 544ae323
......@@ -18,14 +18,15 @@
*****************************************************************************/
/*!
* \file
*
* \brief This file contains the data which is required to calculate the
* component fluxes over a face of a finite volume.
*
* This means concentration gradients, diffusion coefficients, mass fractions, etc.
* at the integration point.
*/
#ifndef DUMUX_STOKES2C_FLUX_VARIABLES_HH
#define DUMUX_STOKES2C_FLUX_VARIABLES_HH
#ifndef DUMUX_STOKESNC_FLUX_VARIABLES_HH
#define DUMUX_STOKESNC_FLUX_VARIABLES_HH
#include <dumux/common/math.hh>
#include <dumux/freeflow/stokes/stokesfluxvariables.hh>
......@@ -33,23 +34,18 @@
namespace Dumux
{
namespace Properties
{
NEW_PROP_TAG(Stokes2cIndices); //!< Enumerations for the compositional stokes models
}
/*!
* \ingroup BoxStokes2cModel
* \ingroup ImplicitFluxVariables
* \ingroup BoxStokesncModel
* \ingroup BoxFluxVariables
* \brief This template class contains data which is required to
* calculate the component fluxes over a face of a finite
* volume for the compositional Stokes model.
* volume for the compositional n component Stokes model.
*
* This means concentration gradients, diffusion coefficient, mass fractions, etc.
* at the integration point of a SCV or boundary face.
*/
template <class TypeTag>
class Stokes2cFluxVariables : public StokesFluxVariables<TypeTag>
class StokesncFluxVariables : public StokesFluxVariables<TypeTag>
{
typedef StokesFluxVariables<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
......@@ -59,97 +55,113 @@ class Stokes2cFluxVariables : public StokesFluxVariables<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
enum { dim = GridView::dimension };
enum { phaseIdx = Indices::phaseIdx };
enum { transportCompIdx = Indices::transportCompIdx };
typedef typename GridView::template Codim<0>::Entity Element;
//dimensions
enum { dim = GridView::dimension };
//phase indices
enum { phaseIdx = Indices::phaseIdx };
//component indices
enum { phaseCompIdx = Indices::phaseCompIdx,
transportCompIdx = Indices::transportCompIdx };
//number of components
enum { numComponents = Indices::numComponents };
typedef typename GridView::template Codim<0>::Entity Element;
typedef Dune::FieldVector<Scalar, dim> DimVector;
typedef Dune::FieldMatrix<Scalar, dim, dim> DimMatrix;
public:
Stokes2cFluxVariables(const Problem &problem,
//Constrcutor calls ParentType function
StokesncFluxVariables(const Problem &problem,
const Element &element,
const FVElementGeometry &fvGeometry,
const int faceIdx,
const ElementVolumeVariables &elemVolVars,
const bool onBoundary = false)
: ParentType(problem, element, fvGeometry, faceIdx, elemVolVars, onBoundary)
: ParentType(problem, element, fvGeometry, faceIdx, elemVolVars, onBoundary)
{
calculateValues_(problem, element, elemVolVars);
}
/*!
/*!
* \brief Return the molar density \f$ \mathrm{[mol/m^3]} \f$ at the integration point.
*/
const Scalar molarDensity() const
{ return this->molarDensity_; }
/*!
* \brief Return the mass fraction of a transported component at the integration point.
*/
const Scalar massFraction(int compIdx) const
{ return massFraction_[compIdx]; }
/*!
* \brief Return the mass fraction of the transported component at the integration point.
* \brief Return the molar diffusion coefficient of a transported component at the integration point.
*/
const Scalar massFraction() const
{ return massFraction_; }
const Scalar diffusionCoeff(int compIdx) const
{ return diffusionCoeff_[compIdx]; }
/*!
* \brief Return the molar diffusion coefficient at the integration point.
* \brief Return the gradient of the mole fraction at the integration point.
*/
const Scalar diffusionCoeff() const
{ return diffusionCoeff_; }
const DimVector &moleFractionGrad(int compIdx) const
{ return moleFractionGrad_[compIdx]; }
/*!
/*!
* \brief Return the eddy diffusivity (if implemented).
*/
const Scalar eddyDiffusivity() const
{ return 0; }
/*!
* \brief Return the gradient of the mole fraction at the integration point.
*/
const DimVector &moleFractionGrad() const
{ return moleFractionGrad_; }
protected:
void calculateValues_(const Problem &problem,
const Element &element,
const ElementVolumeVariables &elemVolVars)
{
molarDensity_ = Scalar(0);
massFraction_ = Scalar(0);
diffusionCoeff_ = Scalar(0);
moleFractionGrad_ = Scalar(0);
// calculate gradients and secondary variables at IPs
for (int idx = 0;
idx < this->fvGeometry_.numScv;
idx++) // loop over vertices of the element
{
molarDensity_ += elemVolVars[idx].molarDensity()*
this->face().shapeValue[idx];
massFraction_ += elemVolVars[idx].massFraction(transportCompIdx) *
this->face().shapeValue[idx];
diffusionCoeff_ += elemVolVars[idx].diffusionCoeff() *
this->face().shapeValue[idx];
// the gradient of the mass fraction at the IP
for (int dimIdx=0; dimIdx<dim; ++dimIdx)
{
moleFractionGrad_ +=
this->face().grad[idx][dimIdx] *
elemVolVars[idx].moleFraction(transportCompIdx);
}
}
Valgrind::CheckDefined(massFraction_);
Valgrind::CheckDefined(diffusionCoeff_);
Valgrind::CheckDefined(moleFractionGrad_);
// loop over all components
for (int compIdx=0; compIdx<numComponents; compIdx++){
if (phaseCompIdx!=compIdx) //no transport equationen parameters needed for the mass balance
{
this->molarDensity_ = Scalar(0.0);
massFraction_[compIdx] = Scalar(0.0);
diffusionCoeff_[compIdx] = Scalar(0.0);
moleFractionGrad_[compIdx] = Scalar(0.0);
// calculate gradients and secondary variables at IPs
for (int scvIdx = 0;
scvIdx < this->fvGeometry_.numScv;
scvIdx++) // loop over vertices of the element
{
this->molarDensity_ += elemVolVars[scvIdx].molarDensity()*
this->face().shapeValue[scvIdx];
massFraction_[compIdx] += elemVolVars[scvIdx].fluidState().massFraction(phaseIdx, compIdx) *
this->face().shapeValue[scvIdx];
diffusionCoeff_[compIdx] += elemVolVars[scvIdx].diffusionCoeff(compIdx) *
this->face().shapeValue[scvIdx];
// the gradient of the mole fraction at the IP
for (int dimIdx=0; dimIdx<dim; ++dimIdx)
{
moleFractionGrad_[compIdx] +=
this->face().grad[scvIdx][dimIdx] *
elemVolVars[scvIdx].fluidState().moleFraction(phaseIdx, compIdx);
}
}
Valgrind::CheckDefined(massFraction_[compIdx]);
Valgrind::CheckDefined(diffusionCoeff_[compIdx]);
Valgrind::CheckDefined(moleFractionGrad_[compIdx]);
}
}
}
Scalar molarDensity_;
Scalar massFraction_;
Scalar diffusionCoeff_;
DimVector moleFractionGrad_;
Scalar molarDensity_;
Scalar massFraction_[numComponents];
Scalar diffusionCoeff_[numComponents];
DimVector moleFractionGrad_[numComponents];
};
} // end namespace
......
......@@ -19,11 +19,13 @@
/*!
* \file
* \brief Defines the indices required for the compositional Stokes box model.
*
* \brief Defines the indices required for the compositional n component Stokes box model.
*/
#ifndef DUMUX_STOKES2C_INDICES_HH
#define DUMUX_STOKES2C_INDICES_HH
#ifndef DUMUX_STOKESNC_INDICES_HH
#define DUMUX_STOKESNC_INDICES_HH
#include "stokesncproperties.hh"
#include <dumux/freeflow/stokes/stokesindices.hh>
namespace Dumux
......@@ -31,27 +33,41 @@ namespace Dumux
// \{
/*!
* \ingroup BoxStokes2cModel
* \ingroup ImplicitIndices
* \brief The common indices for the compositional Stokes box model.
* \ingroup BoxStokesncModel
* \ingroup BoxIndices
* \brief The common indices for the compositional n component Stokes box model.
*
* \tparam PVOffset The first index in a primary variable vector.
*/
template <class TypeTag, int PVOffset = 0>
struct Stokes2cCommonIndices : public StokesCommonIndices<TypeTag>
struct StokesncCommonIndices : public StokesCommonIndices<TypeTag>
{
// Phase index
static const int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx); //!< Index of the employed phase in case of a two-phase fluidsystem (set by default to nPhase)
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
// Component indices
static const int phaseCompIdx = phaseIdx; //!< The index of the main component of the considered phase
static const int transportCompIdx = (unsigned int)(1-phaseIdx); //!< The index of the transported (minor) component; ASSUMES phase indices of 0 and 1
public:
// Dimension (copied for convenience)
static const int dim = StokesCommonIndices<TypeTag>::dim; //!< Number of dimensions
// equation and primary variable indices
static const int dim = StokesCommonIndices<TypeTag>::dim;
static const int transportEqIdx = PVOffset + dim+1; //!< The index for the transport equation
// Phase Index
static const int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx); //!< Index of the employed phase in case of a two-phase fluidsystem (set by default to nPhase)
static const int massOrMoleFracIdx = transportEqIdx; //!< The index of the mass or mole fraction of the transported component in primary variable vectors
// Number of Components
static const int numComponents = FluidSystem::numComponents; //!< Number of components in employed fluidsystem
// Component indices
static const int phaseCompIdx = phaseIdx; //!< The index of the main component of the considered phase
static const int transportCompIdx = (unsigned int)(1-phaseIdx); //!< The index of the first transported component; ASSUMES phase indices of 0 and 1
// Transport equation indices
static const int conti0EqIdx = PVOffset + dim; //!< The index of the mass conservation equation of the first component. In analogy to porous media models "conti" is used here to describe mass conservation equations, i.e total mass balance and transport equations.
static const int massBalanceIdx = conti0EqIdx + phaseCompIdx; //!< The index of the mass balance equation sits on the slot of the employed phase
static const int transportEqIdx = conti0EqIdx + transportCompIdx; //!< The index of the transport equation for a two component model. For n>2 please don't use this index, because it looses its actual meaning.
// Primary variables
static const int massOrMoleFracIdx = transportEqIdx; //!< The index of the first mass or mole fraction of the transported component in primary variable vectors
static const int pressureIdx = massBalanceIdx; //!< The index of the pressure in primary variable vectors
};
} // end namespace
......
......@@ -23,49 +23,73 @@
* using the compositional Stokes box model.
*
*/
#ifndef DUMUX_STOKES2C_LOCAL_RESIDUAL_HH
#define DUMUX_STOKES2C_LOCAL_RESIDUAL_HH
#ifndef DUMUX_STOKESNC_LOCAL_RESIDUAL_HH
#define DUMUX_STOKESNC_LOCAL_RESIDUAL_HH
#include <dumux/freeflow/stokes/stokeslocalresidual.hh>
#include <dumux/freeflow/stokes2c/stokes2cvolumevariables.hh>
#include <dumux/freeflow/stokes2c/stokes2cfluxvariables.hh>
#include "stokesncvolumevariables.hh"
#include "stokesncfluxvariables.hh"
namespace Dumux
{
/*!
* \ingroup BoxStokes2cModel
* \ingroup ImplicitLocalResidual
* \ingroup BoxStokesncModel
* \ingroup BoxLocalResidual
* \brief Element-wise calculation of the Jacobian matrix for problems
* using the compositional Stokes box model. This is derived
* from the Stokes box model.
*/
template<class TypeTag>
class Stokes2cLocalResidual : public StokesLocalResidual<TypeTag>
class StokesncLocalResidual : public StokesLocalResidual<TypeTag>
{
typedef StokesLocalResidual<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
enum { dim = GridView::dimension };
enum { transportEqIdx = Indices::transportEqIdx }; //!< Index of the transport equation
enum { phaseIdx = Indices::phaseIdx }; //!< Index of the considered phase (only of interest when using two-phase fluidsystems)
// component indices
enum { phaseCompIdx = Indices::phaseCompIdx }; //!< Index of the main component of the fluid phase
enum { transportCompIdx = Indices::transportCompIdx }; //!< Index of the minor component of the fluid phase
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
//dimensions
enum { dim = GridView::dimension };
//number of equations
enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
//number of components
enum { numComponents = Indices::numComponents };
//equation indices
enum { massBalanceIdx = Indices::massBalanceIdx,
momentumXIdx = Indices::momentumXIdx,
lastMomentumIdx = Indices::lastMomentumIdx,
transportEqIdx = Indices::transportEqIdx,
conti0EqIdx = Indices::conti0EqIdx };
//primary variable indices
enum { pressureIdx = Indices::pressureIdx };
//phase employed
enum { phaseIdx = Indices::phaseIdx };
//component indices
enum { phaseCompIdx = Indices::phaseCompIdx,
transportCompIdx = Indices::transportCompIdx };
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef Dune::FieldVector<Scalar, dim> DimVector;
typedef typename GridView::Intersection Intersection;
typedef typename GridView::IntersectionIterator IntersectionIterator;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
static const bool calculateNavierStokes = GET_PROP_VALUE(TypeTag, EnableNavierStokes);
//! property that defines whether mole or mass fractions are used
static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
public:
/*!
* \brief Evaluate the stored amount of quantities additional to the Stokes model
* (transport equation).
* (transport equations). For using mole fraction also momentum balances and mass balance
* have to be calculated using molar quantities
*
* The result should be averaged over the volume (e.g. phase mass
* inside a sub control volume divided by the volume)
......@@ -76,24 +100,60 @@ public:
*/
void computeStorage(PrimaryVariables &storage, const int scvIdx, const bool usePrevSol) const
{
// compute the storage term for the transport equation
ParentType::computeStorage(storage, scvIdx, usePrevSol);
// if flag usePrevSol is set, the solution from the previous
// if flag usePrevSol is set, the solution from the previous
// time step is used, otherwise the current solution is
// used. The secondary variables are used accordingly. This
// is required to compute the derivative of the storage term
// using the implicit euler method.
const ElementVolumeVariables &elemVolVars = usePrevSol ? this->prevVolVars_() : this->curVolVars_();
// using the implicit Euler method.
const ElementVolumeVariables &elemVolVars = usePrevSol ? this->prevVolVars_()
: this->curVolVars_();
const VolumeVariables &volVars = elemVolVars[scvIdx];
storage = 0.0;
if (!useMoles)
{
/* works for a maximum of two components, for more components
mole fractions must be used (property useMoles => true) */
// mass and momentum balance
ParentType::computeStorage(storage, scvIdx, usePrevSol);
//storage of transported component
storage[transportEqIdx] = volVars.density()
* volVars.fluidState().massFraction(phaseIdx, transportCompIdx);
Valgrind::CheckDefined(volVars.density());
Valgrind::CheckDefined(volVars.fluidState().massFraction(phaseIdx, transportCompIdx));
}
else
{
/*//TODO call parent type function
// momentum balance
for (int momentumIdx = momentumXIdx; momentumIdx <= lastMomentumIdx; ++momentumIdx)
storage[momentumIdx] = volVars.molarDensity()
* volVars.velocity()[momentumIdx-momentumXIdx];*/
// mass and momentum balance
ParentType::computeStorage(storage, scvIdx, usePrevSol);
// compute the storage of the component
storage[transportEqIdx] =
volVars.density() *
volVars.massFraction(transportCompIdx);
Valgrind::CheckDefined(volVars.density());
Valgrind::CheckDefined(volVars.massFraction(transportCompIdx));
// mass balance and transport equations
for (int compIdx=0; compIdx<numComponents; compIdx++)
{
if (conti0EqIdx+compIdx != massBalanceIdx)
//storage[massBalanceIdx] = volVars.molarDensity();
//else // transport equations
{
storage[conti0EqIdx+compIdx] = volVars.molarDensity()
* volVars.fluidState().moleFraction(phaseIdx, compIdx);
Valgrind::CheckDefined(volVars.molarDensity());
Valgrind::CheckDefined(volVars.fluidState().moleFraction(phaseIdx, compIdx));
}
}
}
}
/*!
......@@ -109,26 +169,61 @@ public:
void computeAdvectiveFlux(PrimaryVariables &flux,
const FluxVariables &fluxVars) const
{
// call computation of the advective fluxes of the stokes model
// (momentum and mass fluxes)
ParentType::computeAdvectiveFlux(flux, fluxVars);
// vertex data of the upstream and the downstream vertices
const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx());
const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx());
Scalar tmp = fluxVars.normalVelocity();
if (this->massUpwindWeight_ > 0.0)
tmp *= this->massUpwindWeight_ * // upwind data
up.density() * up.massFraction(transportCompIdx);
if (this->massUpwindWeight_ < 1.0)
tmp += (1.0 - this->massUpwindWeight_) * // rest
dn.density() * dn.massFraction(transportCompIdx);
flux[transportEqIdx] += tmp;
Valgrind::CheckDefined(flux[transportEqIdx]);
}
// data attached to upstream and the downstream vertices
const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx());
const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx());
Scalar tmp = 0.0;
if(!useMoles)
{
// call ParentType function
ParentType::computeAdvectiveFlux(flux,fluxVars);
// for transport equations
tmp = fluxVars.normalVelocity();
if (this->massUpwindWeight_ > 0.0)
tmp *= this->massUpwindWeight_ // upwind data
* up.density()
* up.fluidState().massFraction(phaseIdx, transportCompIdx);
if (this->massUpwindWeight_ < 1.0)
tmp += (1.0 - this->massUpwindWeight_) // rest
* dn.density()
* dn.fluidState().massFraction(phaseIdx, transportCompIdx);
flux[transportEqIdx] += tmp;
Valgrind::CheckDefined(flux[transportEqIdx]);
}
else
{
// call ParentType function
ParentType::computeAdvectiveFlux(flux,fluxVars);
//transport equations
for (int compIdx=0; compIdx<numComponents; compIdx++)
{
if (conti0EqIdx+compIdx != massBalanceIdx) //mass balance is calculated above
{
tmp = fluxVars.normalVelocity();
if (this->massUpwindWeight_ > 0.0)
tmp *= this->massUpwindWeight_ // upwind data
* up.molarDensity()
* up.fluidState().moleFraction(phaseIdx, compIdx);
if (this->massUpwindWeight_ < 1.0)
tmp += (1.0 - this->massUpwindWeight_) // rest
* dn.molarDensity()
* dn.fluidState().moleFraction(phaseIdx, compIdx);
flux[conti0EqIdx+compIdx] += tmp;
Valgrind::CheckDefined(flux[conti0EqIdx+compIdx]);
}
}
}
}
/*!
* \brief Adds the diffusive component flux to the flux vector over
......@@ -140,19 +235,36 @@ public:
void computeDiffusiveFlux(PrimaryVariables &flux,
const FluxVariables &fluxVars) const
{
// diffusive mass flux
ParentType::computeDiffusiveFlux(flux, fluxVars);
// diffusive component flux
for (int dimIdx = 0; dimIdx < dim; ++dimIdx){
// diffusive component flux
for (int dimIdx = 0; dimIdx < dim; ++dimIdx)
flux[transportEqIdx] -=
fluxVars.moleFractionGrad()[dimIdx] *
fluxVars.face().normal[dimIdx] *
(fluxVars.diffusionCoeff() + fluxVars.eddyDiffusivity()) *
fluxVars.molarDensity() *
FluidSystem::molarMass(transportCompIdx);
if(!useMoles)
{
flux[transportEqIdx] -= fluxVars.moleFractionGrad(transportCompIdx)[dimIdx]
* fluxVars.face().normal[dimIdx]
*(fluxVars.diffusionCoeff(transportCompIdx) + fluxVars.eddyDiffusivity())
* fluxVars.molarDensity()
* FluidSystem::molarMass(transportCompIdx);// Multipled by molarMass [kg/mol] to convert form [mol/m^3 s] to [kg/m^3 s]
Valgrind::CheckDefined(flux[transportEqIdx]);
}
else
{
//loop over secondary components
for (int compIdx=0; compIdx<numComponents; compIdx++)
{
if (conti0EqIdx+compIdx != massBalanceIdx)
{
flux[conti0EqIdx+compIdx] -= fluxVars.moleFractionGrad(compIdx)[dimIdx]
* fluxVars.face().normal[dimIdx]
*(fluxVars.diffusionCoeff(compIdx) + fluxVars.eddyDiffusivity())
* fluxVars.molarDensity();
Valgrind::CheckDefined(flux[conti0EqIdx+compIdx]);
}
}
}
Valgrind::CheckDefined(flux[transportEqIdx]);
}
}
};
......
......@@ -19,23 +19,25 @@
/*!
* \file
*
* \brief Adaptation of the box scheme to the two-component Stokes model.
* \brief Adaptation of the box scheme to the n-component Stokes model.
*/
#ifndef DUMUX_STOKES2C_MODEL_HH
#define DUMUX_STOKES2C_MODEL_HH
#ifndef DUMUX_STOKESNC_MODEL_HH
#define DUMUX_STOKESNC_MODEL_HH