Commit a916147c authored by Katherina Baber's avatar Katherina Baber
Browse files

all indices are now set in 1p2cindices.hh

changed x1Idx for second component to xIdx to avoid confusion
added molar concentration gradient to 1p2cfluxvariables.hh
usage of comp2Idx instead of 1 bei 2-point gradients,
calculate concentrationgradient mit concentrations in
1p2cfluxvariables.hh
splitted computeflux in advective and diffusice component
1p2c/1p2clocalresidual.hh
temperature from problem-file in 1p2cvolumevariables.hh



git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@5304 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 69b006f3
......@@ -48,14 +48,14 @@ class OnePTwoCFluidState : public FluidState<typename GET_PROP_TYPE(TypeTag, PTA
enum {
pressureIdx = Indices::pressureIdx,
x1Idx = Indices::x1Idx,
xIdx = Indices::xIdx,
contiEqIdx = Indices::contiEqIdx,
transEqIdx = Indices::transEqIdx,
phaseIndex = GET_PROP_VALUE(TypeTag, PTAG(PhaseIndex)),
comp1Index = GET_PROP_VALUE(TypeTag, PTAG(Comp1Index)),
comp2Index = GET_PROP_VALUE(TypeTag, PTAG(Comp2Index)),
phaseIdx = Indices::phaseIdx,
comp1Idx = Indices::comp1Idx,
comp2Idx = Indices::comp2Idx,
};
public:
......@@ -76,15 +76,15 @@ public:
temperature_ = temperature;
phasePressure_ = primaryVars[pressureIdx];
x1_ = primaryVars[x1Idx];
x_ = primaryVars[xIdx]; //mole fraction of component 2
meanMolarMass_ =
(1 - x1_)*FluidSystem::molarMass(comp1Index) +
(x1_ )*FluidSystem::molarMass(comp2Index);
(1 - x_)*FluidSystem::molarMass(comp1Idx) +
(x_ )*FluidSystem::molarMass(comp2Idx);
density_ = FluidSystem::phaseDensity(phaseIndex, temperature_, phasePressure_, *this);
density_ = FluidSystem::phaseDensity(phaseIdx, temperature_, phasePressure_, *this);
molarDensity_ = density_ / meanMolarMass_;
Valgrind::CheckDefined(x1_);
Valgrind::CheckDefined(x_);
Valgrind::CheckDefined(phasePressure_);
Valgrind::CheckDefined(density_);
Valgrind::CheckDefined(meanMolarMass_);
......@@ -97,8 +97,8 @@ public:
*
* \param phaseIdx The index of the considered phase
*/
Scalar saturation(int phaseIdx) const
{ return (phaseIndex == phaseIdx)?1.0:0.0; };
Scalar saturation(int phaseIndex) const
{ return (phaseIdx == phaseIndex)?1.0:0.0; };
/*!
* \brief Returns the molar fraction of a component in a fluid phase.
......@@ -106,15 +106,15 @@ public:
* \param phaseIdx The index of the considered phase
* \param compIdx The index of the considered component
*/
Scalar moleFrac(int phaseIdx, int compIdx) const
Scalar moleFrac(int phaseIndex, int compIdx) const
{
// we are a single phase model!
if (phaseIdx != phaseIndex) return 0.0;
if (phaseIndex != phaseIdx) return 0.0;
if (compIdx==comp1Index)
return 1-x1_;
else if (compIdx==comp2Index)
return x1_;
if (compIdx==comp1Idx)
return 1-x_;
else if (compIdx==comp2Idx)
return x_;
return 0.0;
}
......@@ -124,9 +124,9 @@ public:
* This is equivalent to the sum of all component concentrations.
* \param phaseIdx The index of the considered phase
*/
Scalar phaseConcentration(int phaseIdx) const
Scalar phaseConcentration(int phaseIndex) const
{
if (phaseIdx != phaseIndex)
if (phaseIndex != phaseIdx)
return 0;
return density_/meanMolarMass_;
};
......@@ -137,8 +137,8 @@ public:
* \param phaseIdx The index of the considered phase
* \param compIdx The index of the considered component
*/
Scalar concentration(int phaseIdx, int compIdx) const
{ return phaseConcentration(phaseIdx)*moleFrac(phaseIdx, compIdx); };
Scalar concentration(int phaseIndex, int compIdx) const
{ return phaseConcentration(phaseIndex)*moleFrac(phaseIndex, compIdx); };
/*!
......@@ -147,12 +147,12 @@ public:
* \param phaseIdx The index of the considered phase
* \param compIdx The index of the considered component
*/
Scalar massFrac(int phaseIdx, int compIdx) const
Scalar massFrac(int phaseIndex, int compIdx) const
{
if (phaseIdx != phaseIndex)
if (phaseIndex != phaseIdx)
return 0;
return
moleFrac(phaseIdx, compIdx)*
moleFrac(phaseIndex, compIdx)*
FluidSystem::molarMass(compIdx)
/
meanMolarMass_;
......@@ -164,16 +164,16 @@ public:
* \param phaseIdx The index of the considered phase
*
*/
Scalar density(int phaseIdx) const
Scalar density(int phaseIndex) const
{
if (phaseIdx != phaseIndex)
if (phaseIndex != phaseIdx)
return 0;
return density_;
}
Scalar molarDensity(int phaseIdx) const
Scalar molarDensity(int phaseIndex) const
{
if (phaseIdx != phaseIndex)
if (phaseIndex != phaseIdx)
return 0;
return molarDensity_;
}
......@@ -186,9 +186,9 @@ public:
*
* \param phaseIdx The index of the considered phase
*/
Scalar meanMolarMass(int phaseIdx) const
Scalar meanMolarMass(int phaseIndex) const
{
if (phaseIdx != phaseIndex)
if (phaseIndex != phaseIdx)
return 0;
return meanMolarMass_;
};
......@@ -198,7 +198,7 @@ public:
*
* \param phaseIdx The index of the considered phase
*/
Scalar phasePressure(int phaseIdx) const
Scalar phasePressure(int phaseIndex) const
{
return phasePressure_;
}
......@@ -213,7 +213,7 @@ public:
{ return temperature_; };
public:
Scalar x1_;
Scalar x_;
Scalar phasePressure_;
Scalar density_;
Scalar molarDensity_;
......
......@@ -59,6 +59,7 @@ class OnePTwoCFluxVariables
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters)) SpatialParameters;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(OnePTwoCIndices)) Indices;
enum { dim = GridView::dimension };
enum { dimWorld = GridView::dimensionworld };
......@@ -70,7 +71,11 @@ class OnePTwoCFluxVariables
typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
typedef typename GridView::template Codim<0>::Entity Element;
enum {
phaseIdx = Indices::phaseIdx,
comp1Idx = Indices::comp1Idx,
comp2Idx = Indices::comp2Idx,
};
public:
/*
* \brief The constructor
......@@ -136,6 +141,19 @@ public:
return concentrationGrad_;
};
/*!
* \brief The molar concentration gradient of a component in a phase.
*/
const Vector &molarConcGrad(int compIdx) const
{
if (compIdx != 1)
{ DUNE_THROW(Dune::InvalidStateException,
"The 1p2c model is supposed to need "
"only the concentration gradient of "
"the second component!"); }
return molarConcGrad_;
};
Scalar porousDiffCoeff() const
{
// TODO: tensorial diffusion coefficients
......@@ -189,6 +207,7 @@ protected:
potentialGrad_ = 0.0;
concentrationGrad_ = 0.0;
molarConcGrad_ = 0.0;
Vector tmp;
//The decision of the if-statement depends on the function useTwoPointGradient(const Element &elem,
......@@ -209,9 +228,13 @@ protected:
// the concentration gradient
tmp = feGrad;
tmp *= elemDat[idx].concentration(1);
tmp *= elemDat[idx].concentration(comp2Idx);
concentrationGrad_ += tmp;
tmp = feGrad;
tmp *= elemDat[idx].moleFrac(comp2Idx);
molarConcGrad_ += tmp;
// phase viscosity
viscosityAtIP_ += elemDat[idx].viscosity()*face().shapeValue[idx];
......@@ -231,7 +254,9 @@ protected:
potentialGrad_ = tmp;
potentialGrad_ *= vVars_j.pressure() - vVars_i.pressure();
concentrationGrad_ = tmp;
concentrationGrad_ *= vVars_j.moleFrac(1) - vVars_i.moleFrac(1);
concentrationGrad_ *= vVars_j.concentration(comp2Idx) - vVars_i.concentration(comp2Idx);
molarConcGrad_ = tmp;
molarConcGrad_ *= vVars_j.moleFrac(comp2Idx) - vVars_i.moleFrac(comp2Idx);
}
// correct the pressure by the hydrostatic pressure due to
......@@ -340,6 +365,8 @@ protected:
Vector potentialGrad_;
//! concentratrion gradient
Vector concentrationGrad_;
//! molar concentratrion gradient
Vector molarConcGrad_;
//! the effective diffusion coefficent in the porous medium
Scalar diffCoeffPM_;
......
......@@ -43,13 +43,22 @@ namespace Dumux
template <int PVOffset = 0>
struct OnePTwoCIndices
{
//! Set the default phase used by the fluid system to the first one
static const int phaseIdx = 0;
//! Set the default for the first component the fluid system's first component
static const int comp1Idx = 0;
//! Set the default for the second component the fluid system's second component
static const int comp2Idx = 1;
// Equation indices
static const int contiEqIdx = PVOffset + 0; //!< continuity equation index
static const int transEqIdx = PVOffset + 1; //!< transport equation index
// primary variable indices
static const int pressureIdx = PVOffset + 0; //!< pressure
static const int x1Idx = PVOffset + 1; //!< mole fraction of the second component
static const int xIdx = PVOffset + 1; //!< mole fraction of the second component
//in my case the therapeutic agent
};
......
......@@ -56,6 +56,7 @@ class OnePTwoCLocalResidual : public BoxLocalResidual<TypeTag>
{
protected:
typedef OnePTwoCLocalResidual<TypeTag> ThisType;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalResidual)) Implementation;
typedef BoxLocalResidual<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
......@@ -73,7 +74,11 @@ protected:
// indices of the primary variables
pressureIdx = Indices::pressureIdx,
x1Idx = Indices::x1Idx,
xIdx = Indices::xIdx,
phaseIdx = Indices::phaseIdx,
comp1Idx = Indices::comp1Idx,
comp2Idx = Indices::comp2Idx,
// indices of the equations
contiEqIdx = Indices::contiEqIdx,
......@@ -115,7 +120,7 @@ public:
// storage term of the transport equation
result[transEqIdx] =
volVars.concentration(1) *
volVars.concentration(comp2Idx) *
volVars.porosity();
}
......@@ -135,7 +140,23 @@ public:
faceId,
this->curVolVars_());
Vector tmpVec;
asImp_()->computeAdvectiveFlux(flux, fluxVars);
asImp_()->computeDiffusiveFlux(flux, fluxVars);
}
/*!
* \brief Evaluates the advective mass flux of all components over
* a face of a subcontrol volume.
*
* \param flux The advective flux over the sub-control-volume face for each component
* \param vars The flux variables at the current SCV
*/
void computeAdvectiveFlux(PrimaryVariables &flux, const FluxVariables &fluxVars) const
{
////////
// advective fluxes of all components in all phases
////////
Vector tmpVec(0);
fluxVars.intrinsicPermeability().mv(fluxVars.potentialGrad(), tmpVec);
......@@ -154,22 +175,32 @@ public:
// advective flux of the second component
flux[transEqIdx] +=
normalFlux *
(( upwindAlpha)*up.concentration(1)/up.viscosity()
(( upwindAlpha)*up.concentration(comp2Idx)/up.viscosity()
+
(1 - upwindAlpha)*dn.concentration(1)/dn.viscosity());
(1 - upwindAlpha)*dn.concentration(comp2Idx)/dn.viscosity());
// diffusive flux of second component
flux[transEqIdx] -=
fluxVars.porousDiffCoeff() *
(fluxVars.concentrationGrad(1) * fluxVars.face().normal);
// dispersive flux of second component
Vector normalDisp;
fluxVars.dispersionTensor().mv(fluxVars.face().normal, normalDisp);
flux[transEqIdx] -=
(normalDisp * fluxVars.concentrationGrad(1));
}
/*!
* \brief Adds the diffusive mass flux of all components over
* a face of a subcontrol volume.
*
* \param flux The diffusive flux over the sub-control-volume face for each component
* \param vars The flux variables at the current SCV
*/
void computeDiffusiveFlux(PrimaryVariables &flux, const FluxVariables &fluxVars) const
{
// diffusive flux of second component
flux[transEqIdx] -=
fluxVars.porousDiffCoeff() *
(fluxVars.concentrationGrad(1) * fluxVars.face().normal);
// dispersive flux of second component
Vector normalDisp;
fluxVars.dispersionTensor().mv(fluxVars.face().normal, normalDisp);
flux[transEqIdx] -=
(normalDisp * fluxVars.concentrationGrad(comp2Idx));
}
/*!
* \brief Calculate the source term of the equation
* \param q The source/sink in the SCV for each component
......@@ -183,6 +214,11 @@ public:
this->fvElemGeom_(),
localVertexIdx);
}
Implementation *asImp_()
{ return static_cast<Implementation *> (this); }
const Implementation *asImp_() const
{ return static_cast<const Implementation *> (this); }
};
}
......
......@@ -59,9 +59,6 @@ NEW_PROP_TAG(SpatialParameters); //!< The type of the spatial parameters
NEW_PROP_TAG(FluidSystem); //!< Type of the multi-component relations
NEW_PROP_TAG(UpwindAlpha); //!< The default value of the upwind parameter
NEW_PROP_TAG(EnableGravity); //!< Returns whether gravity is considered in the problem
NEW_PROP_TAG(PhaseIndex); //!< The index of the phase in the fluid system
NEW_PROP_TAG(Comp1Index); //!< The index of the first component in the fluid system
NEW_PROP_TAG(Comp2Index); //!< The index of the second component in the fluid system
}
// \}
}
......
......@@ -66,20 +66,22 @@ SET_TYPE_PROP(BoxOnePTwoC, VolumeVariables, OnePTwoCVolumeVariables<TypeTag>);
//! define the FluxVariables
SET_TYPE_PROP(BoxOnePTwoC, FluxVariables, OnePTwoCFluxVariables<TypeTag>);
// the BoundaryVariables property
SET_TYPE_PROP(BoxOnePTwoC, BoundaryVariables, OnePTwoCBoundaryVariables<TypeTag>);
//! set default upwind factor to 1.0, i.e. fully upwind
SET_SCALAR_PROP(BoxOnePTwoC, UpwindAlpha, 1.0);
//! Set the indices used by the 1p2c model
SET_TYPE_PROP(BoxOnePTwoC, OnePTwoCIndices, Dumux::OnePTwoCIndices<0>);
//! Set the default phase used by the fluid system to the first one
SET_INT_PROP(BoxOnePTwoC, PhaseIndex, 0);
//! Set the default for the first component the fluid system's first component
SET_INT_PROP(BoxOnePTwoC, Comp1Index, 0);
// enable jacobian matrix recycling by default
SET_BOOL_PROP(BoxOnePTwoC, EnableJacobianRecycling, true);
// enable partial reassembling by default
SET_BOOL_PROP(BoxOnePTwoC, EnablePartialReassemble, true);
// enable time-step ramp up by default
SET_BOOL_PROP(BoxOnePTwoC, EnableTimeStepRampUp, false);
//! Set the default for the second component the fluid system's second component
SET_INT_PROP(BoxOnePTwoC, Comp2Index, 1);
}
// \}
}
......
......@@ -62,9 +62,9 @@ class OnePTwoCVolumeVariables : public BoxVolumeVariables<TypeTag>
dimWorld = GridView::dimensionworld,
phaseIndex = GET_PROP_VALUE(TypeTag, PTAG(PhaseIndex)),
comp1Index = GET_PROP_VALUE(TypeTag, PTAG(Comp1Index)),
comp2Index = GET_PROP_VALUE(TypeTag, PTAG(Comp2Index)),
phaseIdx = Indices::phaseIdx,
comp1Idx = Indices::comp1Idx,
comp2Idx = Indices::comp2Idx,
contiEqIdx = Indices::contiEqIdx,
transEqIdx = Indices::transEqIdx
......@@ -93,22 +93,28 @@ public:
bool isOldSol)
{
ParentType::update(priVars, problem, element, elemGeom, scvIdx, isOldSol);
asImp().updateTemperature_(priVars,
element,
elemGeom,
scvIdx,
problem);
fluidState_.update(priVars, temperature_);
porosity_ = problem.spatialParameters().porosity(element, elemGeom, scvIdx);
tortuosity_ = problem.spatialParameters().tortuosity(element, elemGeom, scvIdx);
dispersivity_ = problem.spatialParameters().dispersivity(element, elemGeom, scvIdx);
Scalar temperature = problem.temperature(element, elemGeom, scvIdx);
fluidState_.update(priVars, temperature);
viscosity_ = FluidSystem::phaseViscosity(phaseIndex,
temperature,
viscosity_ = FluidSystem::phaseViscosity(phaseIdx,
temperature_,
pressure(),
fluidState_);
diffCoeff_ = FluidSystem::diffCoeff(phaseIndex,
comp1Index,
comp2Index,
temperature,
diffCoeff_ = FluidSystem::diffCoeff(phaseIdx,
comp1Idx,
comp2Idx,
temperature_,
pressure(),
*this);
......@@ -132,10 +138,10 @@ public:
* \brief Returns density the of the fluid phase.
*/
Scalar density() const
{ return fluidState_.density(phaseIndex); }
{ return fluidState_.density(phaseIdx); }
Scalar molarDensity() const
{ return fluidState_.molarDensity(phaseIndex);}
{ return fluidState_.molarDensity(phaseIdx);}
/*!
* \brief Returns mole fraction of a component in the phase
......@@ -143,28 +149,28 @@ public:
* \param compIdx The index of the component
*/
Scalar moleFrac(int compIdx) const
{ return fluidState_.moleFrac(phaseIndex, (compIdx==0)?comp1Index:comp2Index); }
{ return fluidState_.moleFrac(phaseIdx, (compIdx==0)?comp1Idx:comp2Idx); }
/*!
* \brief Returns mass fraction of a component in the phase
* \param compIdx The index of the component
*/
Scalar massFrac(int compIdx) const
{ return fluidState_.massFrac(phaseIndex, (compIdx==0)?comp1Index:comp2Index); }
{ return fluidState_.massFrac(phaseIdx, (compIdx==0)?comp1Idx:comp2Idx); }
/*!
* \brief Returns concentration of a component in the phase
* \param compIdx The index of the component
*/
Scalar concentration(int compIdx) const
{ return fluidState_.concentration(phaseIndex, (compIdx==0)?comp1Index:comp2Index); }
{ return fluidState_.concentration(phaseIdx, (compIdx==0)?comp1Idx:comp2Idx); }
/*!
* \brief Returns the effective pressure of a given phase within
* the control volume.
*/
Scalar pressure() const
{ return fluidState_.phasePressure(phaseIndex); }
{ return fluidState_.phasePressure(phaseIdx); }
/*!
* \brief Returns the binary diffusion coefficient in the fluid
......@@ -192,7 +198,7 @@ public:
* identical.
*/
Scalar temperature() const
{ return fluidState_.temperature(); }
{ return temperature_; }
/*!
* \brief Returns the dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of a given phase
......@@ -208,14 +214,32 @@ public:
{ return porosity_; }
protected:
Scalar porosity_;
void updateTemperature_(const PrimaryVariables &priVars,
const Element &element,
const FVElementGeometry &elemGeom,
int scvIdx,
const Problem &problem)
{
temperature_ = problem.temperature(element, elemGeom, scvIdx);
}
Scalar temperature_; //!< Temperature within the control volume
Scalar porosity_; //!< Effective porosity within the control volume
Scalar viscosity_;
Scalar tortuosity_;
Vector dispersivity_;
Scalar diffCoeff_;
FluidState fluidState_;
private:
Implementation &asImp()
{ return *static_cast<Implementation*>(this); }
const Implementation &asImp() const
{ return *static_cast<const Implementation*>(this); }
};
}
}// end namepace
#endif
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