Commit 10237266 authored by Dominik Riesterer's avatar Dominik Riesterer
Browse files

The model can now be used with either mole or mass fractions. The property

useMoles has to be set in the problem file and the boundary conditions 
have to be choosen accordingly.

In order to be consistent with the multicomponent models (3p3c,
3p3cni,...) the model uses MOLE fractions by default now.

The documentation is changed and extended accordingly.

reviewed by kathinka


git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@11526 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 6f2f0c47
......@@ -74,6 +74,9 @@ class TwoPTwoCLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
static constexpr unsigned int replaceCompEqIdx =
GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
//! property that defines whether mole or mass fractions are used
static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
public:
/*!
* \brief Constructor. Sets the upwind weight.
......@@ -137,24 +140,46 @@ class TwoPTwoCLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
// compute storage term of all components within all phases
storage = 0;
for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
if(!useMoles) //mass fraction formulation
{
for (unsigned int compIdx = contiCompIdx1_(); compIdx <= contiCompIdx2_(); ++compIdx)
{
unsigned int eqIdx = (compIdx == wCompIdx) ? contiWEqIdx : contiNEqIdx;
storage[eqIdx] += volVars.density(phaseIdx)
* volVars.saturation(phaseIdx)
* volVars.fluidState().massFraction(phaseIdx, compIdx);
}
// this is only processed, if one component mass balance equation
// is replaced by the total mass balance equation
if (replaceCompEqIdx < numComponents)
storage[replaceCompEqIdx] +=
volVars.density(phaseIdx)
* volVars.saturation(phaseIdx);
for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
for (unsigned int compIdx = contiCompIdx1_(); compIdx <= contiCompIdx2_(); ++compIdx)
{
unsigned int eqIdx = (compIdx == wCompIdx) ? contiWEqIdx : contiNEqIdx;
storage[eqIdx] += volVars.density(phaseIdx)
* volVars.saturation(phaseIdx)
* volVars.fluidState().massFraction(phaseIdx, compIdx);
}
// this is only processed, if one component mass balance equation
// is replaced by the total mass balance equation
if (replaceCompEqIdx < numComponents)
storage[replaceCompEqIdx] +=
volVars.density(phaseIdx)
* volVars.saturation(phaseIdx);
}
storage *= volVars.porosity();
}
else //mole-fraction formulation
{
for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
for (unsigned int compIdx = contiCompIdx1_(); compIdx <= contiCompIdx2_(); ++compIdx)
{
unsigned int eqIdx = (compIdx == wCompIdx) ? contiWEqIdx : contiNEqIdx;
storage[eqIdx] += volVars.molarDensity(phaseIdx)
* volVars.saturation(phaseIdx)
* volVars.fluidState().moleFraction(phaseIdx, compIdx);
}
// this is only processed, if one component mass balance equation
// is replaced by the total mass balance equation
if (replaceCompEqIdx < numComponents)
storage[replaceCompEqIdx] +=
volVars.molarDensity(phaseIdx)
* volVars.saturation(phaseIdx);
}
storage *= volVars.porosity();
}
storage *= volVars.porosity();
}
/*!
......@@ -193,65 +218,131 @@ class TwoPTwoCLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
////////
// advective fluxes of all components in all phases
////////
for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
// data attached to upstream and the downstream vertices
// of the current phase
const VolumeVariables &up =
this->curVolVars_(fluxVars.upstreamIdx(phaseIdx));
const VolumeVariables &dn =
this->curVolVars_(fluxVars.downstreamIdx(phaseIdx));
for (unsigned int compIdx = contiCompIdx1_(); compIdx <= contiCompIdx2_(); ++compIdx)
if(!useMoles) //mass-fraction formulation
{
for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
// data attached to upstream and the downstream vertices
// of the current phase
const VolumeVariables &up =
this->curVolVars_(fluxVars.upstreamIdx(phaseIdx));
const VolumeVariables &dn =
this->curVolVars_(fluxVars.downstreamIdx(phaseIdx));
for (unsigned int compIdx = contiCompIdx1_(); compIdx <= contiCompIdx2_(); ++compIdx)
{
unsigned int eqIdx = (compIdx == wCompIdx) ? contiWEqIdx : contiNEqIdx;
// add advective flux of current component in current
// phase
if (massUpwindWeight_ > 0.0)
// upstream vertex
flux[eqIdx] +=
fluxVars.volumeFlux(phaseIdx)
* massUpwindWeight_
* up.density(phaseIdx)
* up.fluidState().massFraction(phaseIdx, compIdx);
if (massUpwindWeight_ < 1.0)
// downstream vertex
flux[eqIdx] +=
fluxVars.volumeFlux(phaseIdx)
* (1 - massUpwindWeight_)
* dn.density(phaseIdx)
* dn.fluidState().massFraction(phaseIdx, compIdx);
Valgrind::CheckDefined(fluxVars.volumeFlux(phaseIdx));
Valgrind::CheckDefined(up.density(phaseIdx));
Valgrind::CheckDefined(up.fluidState().massFraction(phaseIdx, compIdx));
Valgrind::CheckDefined(dn.density(phaseIdx));
Valgrind::CheckDefined(dn.fluidState().massFraction(phaseIdx, compIdx));
}
// flux of the total mass balance;
// this is only processed, if one component mass balance equation
// is replaced by a total mass balance equation
if (replaceCompEqIdx < numComponents)
{
// upstream vertex
if (massUpwindWeight_ > 0.0)
flux[replaceCompEqIdx] +=
fluxVars.volumeFlux(phaseIdx)
* massUpwindWeight_
* up.density(phaseIdx);
// downstream vertex
if (massUpwindWeight_ < 1.0)
flux[replaceCompEqIdx] +=
fluxVars.volumeFlux(phaseIdx)
* (1 - massUpwindWeight_)
* dn.density(phaseIdx);
Valgrind::CheckDefined(fluxVars.volumeFlux(phaseIdx));
Valgrind::CheckDefined(up.density(phaseIdx));
Valgrind::CheckDefined(dn.density(phaseIdx));
}
}
}
else //mole-fraction formulation
{
for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
unsigned int eqIdx = (compIdx == wCompIdx) ? contiWEqIdx : contiNEqIdx;
// add advective flux of current component in current
// phase
if (massUpwindWeight_ > 0.0)
// data attached to upstream and the downstream vertices
// of the current phase
const VolumeVariables &up =
this->curVolVars_(fluxVars.upstreamIdx(phaseIdx));
const VolumeVariables &dn =
this->curVolVars_(fluxVars.downstreamIdx(phaseIdx));
for (unsigned int compIdx = contiCompIdx1_(); compIdx <= contiCompIdx2_(); ++compIdx)
{
unsigned int eqIdx = (compIdx == wCompIdx) ? contiWEqIdx : contiNEqIdx;
// add advective flux of current component in current
// phase
if (massUpwindWeight_ > 0.0)
// upstream vertex
flux[eqIdx] +=
fluxVars.volumeFlux(phaseIdx)
* massUpwindWeight_
* up.molarDensity(phaseIdx)
* up.fluidState().moleFraction(phaseIdx, compIdx);
if (massUpwindWeight_ < 1.0)
// downstream vertex
flux[eqIdx] +=
fluxVars.volumeFlux(phaseIdx)
* (1 - massUpwindWeight_)
* dn.molarDensity(phaseIdx)
* dn.fluidState().moleFraction(phaseIdx, compIdx);
Valgrind::CheckDefined(fluxVars.volumeFlux(phaseIdx));
Valgrind::CheckDefined(up.molarDensity(phaseIdx));
Valgrind::CheckDefined(up.fluidState().moleFraction(phaseIdx, compIdx));
Valgrind::CheckDefined(dn.molarDensity(phaseIdx));
Valgrind::CheckDefined(dn.fluidState().moleFraction(phaseIdx, compIdx));
}
// flux of the total mass balance;
// this is only processed, if one component mass balance equation
// is replaced by a total mass balance equation
if (replaceCompEqIdx < numComponents)
{
// upstream vertex
flux[eqIdx] +=
fluxVars.volumeFlux(phaseIdx)
* massUpwindWeight_
* up.density(phaseIdx)
* up.fluidState().massFraction(phaseIdx, compIdx);
if (massUpwindWeight_ < 1.0)
if (massUpwindWeight_ > 0.0)
flux[replaceCompEqIdx] +=
fluxVars.volumeFlux(phaseIdx)
* massUpwindWeight_
* up.molarDensity(phaseIdx);
// downstream vertex
flux[eqIdx] +=
fluxVars.volumeFlux(phaseIdx)
* (1 - massUpwindWeight_)
* dn.density(phaseIdx)
* dn.fluidState().massFraction(phaseIdx, compIdx);
Valgrind::CheckDefined(fluxVars.volumeFlux(phaseIdx));
Valgrind::CheckDefined(up.density(phaseIdx));
Valgrind::CheckDefined(up.fluidState().massFraction(phaseIdx, compIdx));
Valgrind::CheckDefined(dn.density(phaseIdx));
Valgrind::CheckDefined(dn.fluidState().massFraction(phaseIdx, compIdx));
}
// flux of the total mass balance;
// this is only processed, if one component mass balance equation
// is replaced by a total mass balance equation
if (replaceCompEqIdx < numComponents)
{
// upstream vertex
if (massUpwindWeight_ > 0.0)
flux[replaceCompEqIdx] +=
fluxVars.volumeFlux(phaseIdx)
* massUpwindWeight_
* up.density(phaseIdx);
// downstream vertex
if (massUpwindWeight_ < 1.0)
flux[replaceCompEqIdx] +=
fluxVars.volumeFlux(phaseIdx)
* (1 - massUpwindWeight_)
* dn.density(phaseIdx);
Valgrind::CheckDefined(fluxVars.volumeFlux(phaseIdx));
Valgrind::CheckDefined(up.density(phaseIdx));
Valgrind::CheckDefined(dn.density(phaseIdx));
if (massUpwindWeight_ < 1.0)
flux[replaceCompEqIdx] +=
fluxVars.volumeFlux(phaseIdx)
* (1 - massUpwindWeight_)
* dn.molarDensity(phaseIdx);
Valgrind::CheckDefined(fluxVars.volumeFlux(phaseIdx));
Valgrind::CheckDefined(up.molarDensity(phaseIdx));
Valgrind::CheckDefined(dn.molarDensity(phaseIdx));
}
}
}
}
}
}
/*!
......@@ -262,30 +353,56 @@ class TwoPTwoCLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
* \param fluxVars The flux variables at the current sub control volume face
*/
void computeDiffusiveFlux(PrimaryVariables &flux, const FluxVariables &fluxVars) const
{
// add diffusive flux of gas component in liquid phase
Scalar tmp = fluxVars.moleFractionGrad(wPhaseIdx)*fluxVars.face().normal;
tmp *= -1;
tmp *=
fluxVars.porousDiffCoeff(wPhaseIdx) *
fluxVars.molarDensity(wPhaseIdx);
// add the diffusive fluxes only to the component mass balance
if (replaceCompEqIdx != contiNEqIdx)
flux[contiNEqIdx] += tmp * FluidSystem::molarMass(nCompIdx);
if (replaceCompEqIdx != contiWEqIdx)
flux[contiWEqIdx] -= tmp * FluidSystem::molarMass(wCompIdx);
// add diffusive flux of liquid component in non-wetting phase
tmp = fluxVars.moleFractionGrad(nPhaseIdx)*fluxVars.face().normal;
tmp *= -1;
tmp *=
fluxVars.porousDiffCoeff(nPhaseIdx) *
fluxVars.molarDensity(nPhaseIdx);
// add the diffusive fluxes only to the component mass balance
if (replaceCompEqIdx != contiWEqIdx)
flux[contiWEqIdx] += tmp * FluidSystem::molarMass(wCompIdx);
if (replaceCompEqIdx != contiNEqIdx)
flux[contiNEqIdx] -= tmp * FluidSystem::molarMass(nCompIdx);
if(!useMoles) //mass-fraction formulation
{
// add diffusive flux of gas component in liquid phase
Scalar tmp = - (fluxVars.moleFractionGrad(wPhaseIdx)*fluxVars.face().normal);
tmp *=
fluxVars.porousDiffCoeff(wPhaseIdx) *
fluxVars.molarDensity(wPhaseIdx);
// add the diffusive fluxes only to the component mass balance
if (replaceCompEqIdx != contiNEqIdx)
flux[contiNEqIdx] += tmp * FluidSystem::molarMass(nCompIdx);
if (replaceCompEqIdx != contiWEqIdx)
flux[contiWEqIdx] -= tmp * FluidSystem::molarMass(wCompIdx);
// add diffusive flux of liquid component in non-wetting phase
tmp = -(fluxVars.moleFractionGrad(nPhaseIdx)*fluxVars.face().normal);
tmp *=
fluxVars.porousDiffCoeff(nPhaseIdx) *
fluxVars.molarDensity(nPhaseIdx);
// add the diffusive fluxes only to the component mass balance
if (replaceCompEqIdx != contiWEqIdx)
flux[contiWEqIdx] += tmp * FluidSystem::molarMass(wCompIdx);
if (replaceCompEqIdx != contiNEqIdx)
flux[contiNEqIdx] -= tmp * FluidSystem::molarMass(nCompIdx);
}
else //mole-fraction formulation
{
// add diffusive flux of gas component in liquid phase
Scalar tmp = - (fluxVars.moleFractionGrad(wPhaseIdx)*fluxVars.face().normal);
tmp *=
fluxVars.porousDiffCoeff(wPhaseIdx) *
fluxVars.molarDensity(wPhaseIdx);
// add the diffusive fluxes only to the component mass balance
if (replaceCompEqIdx != contiNEqIdx)
flux[contiNEqIdx] += tmp;
if (replaceCompEqIdx != contiWEqIdx)
flux[contiWEqIdx] -= tmp;
// add diffusive flux of liquid component in non-wetting phase
tmp = -(fluxVars.moleFractionGrad(nPhaseIdx)*fluxVars.face().normal);
tmp *=
fluxVars.porousDiffCoeff(nPhaseIdx) *
fluxVars.molarDensity(nPhaseIdx);
// add the diffusive fluxes only to the component mass balance
if (replaceCompEqIdx != contiWEqIdx)
flux[contiWEqIdx] += tmp;
if (replaceCompEqIdx != contiNEqIdx)
flux[contiNEqIdx] -= tmp;
}
}
/*!
......@@ -293,6 +410,7 @@ class TwoPTwoCLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
*
* \param source The source/sink in the sub-control volume for each component
* \param scvIdx The index of the sub-control volume
* \be careful what you use! (mole or mass Fraction!) Think of the units!
*/
void computeSource(PrimaryVariables& source, const int scvIdx)
{
......@@ -306,25 +424,50 @@ class TwoPTwoCLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
protected:
void evalPhaseStorage_(const int phaseIdx)
{
// evaluate the storage terms of a single phase
for (int i=0; i < this->fvGeometry_().numScv; i++) {
PrimaryVariables &storage = this->storageTerm_[i];
const ElementVolumeVariables &elemVolVars = this->curVolVars_();
const VolumeVariables &volVars = elemVolVars[i];
// compute storage term of all components within all phases
storage = 0;
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
{
int eqIdx = (compIdx == wCompIdx) ? contiWEqIdx : contiNEqIdx;
storage[eqIdx] += volVars.density(phaseIdx)
* volVars.saturation(phaseIdx)
* volVars.fluidState().massFraction(phaseIdx, compIdx);
}
storage *= volVars.porosity();
storage *= this->fvGeometry_().subContVol[i].volume;
}
if(!useMoles) //mass-fraction formulation
{
// evaluate the storage terms of a single phase
for (int i=0; i < this->fvGeometry_().numScv; i++) {
PrimaryVariables &storage = this->storageTerm_[i];
const ElementVolumeVariables &elemVolVars = this->curVolVars_();
const VolumeVariables &volVars = elemVolVars[i];
// compute storage term of all components within all phases
storage = 0;
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
{
int eqIdx = (compIdx == wCompIdx) ? contiWEqIdx : contiNEqIdx;
storage[eqIdx] += volVars.density(phaseIdx)
* volVars.saturation(phaseIdx)
* volVars.fluidState().massFraction(phaseIdx, compIdx);
}
storage *= volVars.porosity();
storage *= this->fvGeometry_().subContVol[i].volume;
}
}
else //mole-fraction formulation
{
// evaluate the storage terms of a single phase
for (int i=0; i < this->fvGeometry_().numScv; i++) {
PrimaryVariables &storage = this->storageTerm_[i];
const ElementVolumeVariables &elemVolVars = this->curVolVars_();
const VolumeVariables &volVars = elemVolVars[i];
// compute storage term of all components within all phases
storage = 0;
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
{
int eqIdx = (compIdx == wCompIdx) ? contiWEqIdx : contiNEqIdx;
storage[eqIdx] += volVars.molarDensity(phaseIdx)
* volVars.saturation(phaseIdx)
* volVars.fluidState().moleFraction(phaseIdx, compIdx);
}
storage *= volVars.porosity();
storage *= this->fvGeometry_().subContVol[i].volume;
}
}
}
/*!
......
......@@ -73,14 +73,17 @@ namespace Dumux
* default, the model uses \f$p_w\f$ and \f$S_n\f$.
* Moreover, the second primary variable depends on the phase state, since a
* primary variable switch is included. The phase state is stored for all nodes
* of the system. Following cases can be distinguished:
* of the system.
* The model is able to use either mole or mass fractions. The property useMoles can be set to either true or false in the
* problem file. Make sure that the according units are used in the problem setup. useMoles is set to true by default.
* Following cases can be distinguished:
* <ul>
* <li> Both phases are present: The saturation is used (either \f$S_n\f$ or \f$S_w\f$, dependent on the chosen <tt>Formulation</tt>),
* as long as \f$ 0 < S_\alpha < 1\f$</li>.
* <li> Only wetting phase is present: The mass fraction of, e.g., air in the wetting phase \f$X^a_w\f$ is used,
* as long as the maximum mass fraction is not exceeded \f$(X^a_w<X^a_{w,max})\f$</li>
* as long as the maximum mass/mole fraction is not exceeded \f$(X^a_w<X^a_{w,max})\f$</li>
* <li> Only non-wetting phase is present: The mass fraction of, e.g., water in the non-wetting phase, \f$X^w_n\f$, is used,
* as long as the maximum mass fraction is not exceeded \f$(X^w_n<X^w_{n,max})\f$</li>
* as long as the maximum mass/mole fraction is not exceeded \f$(X^w_n<X^w_{n,max})\f$</li>
* </ul>
*/
......@@ -129,6 +132,7 @@ class TwoPTwoCModel: public GET_PROP_TYPE(TypeTag, BaseModel)
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
enum { dofCodim = isBox ? dim : 0 };
......@@ -301,6 +305,10 @@ public:
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
massFrac[phaseIdx][compIdx] = writer.allocateManagedBuffer(numDofs);
ScalarField *moleFrac[numPhases][numComponents];
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
moleFrac[phaseIdx][compIdx] = writer.allocateManagedBuffer(numDofs);
ScalarField *temperature = writer.allocateManagedBuffer(numDofs);
ScalarField *poro = writer.allocateManagedBuffer(numDofs);
VectorField *velocityN = writer.template allocateManagedBuffer<double, dim>(numDofs);
......@@ -357,6 +365,14 @@ public:
Valgrind::CheckDefined((*massFrac[phaseIdx][compIdx])[globalIdx][0]);
}
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
{
(*moleFrac[phaseIdx][compIdx])[globalIdx]
= elemVolVars[scvIdx].fluidState().moleFraction(phaseIdx, compIdx);
Valgrind::CheckDefined((*moleFrac[phaseIdx][compIdx])[globalIdx][0]);
}
(*poro)[globalIdx] = elemVolVars[scvIdx].porosity();
(*temperature)[globalIdx] = elemVolVars[scvIdx].temperature();
(*phasePresence)[globalIdx]
......@@ -387,6 +403,15 @@ public:
writer.attachDofData(*massFrac[phaseIdx][compIdx], oss.str(), isBox);
}
}
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
{
std::ostringstream oss;
oss << "x_" << FluidSystem::phaseName(phaseIdx) << "^" << FluidSystem::componentName(compIdx);
writer.attachDofData(*moleFrac[phaseIdx][compIdx], oss.str(), isBox);
}
}
writer.attachDofData(*poro, "porosity", isBox);
writer.attachDofData(*temperature, "temperature", isBox);
writer.attachDofData(*phasePresence, "phase presence", isBox);
......@@ -632,8 +657,16 @@ public:
<< volVars.saturation(nPhaseIdx) << std::endl;
newPhasePresence = wPhaseOnly;
globalSol[globalIdx][switchIdx]
= volVars.fluidState().massFraction(wPhaseIdx, nCompIdx);
if(!useMoles) //mass-fraction formulation
{
globalSol[globalIdx][switchIdx]
= volVars.fluidState().massFraction(wPhaseIdx, nCompIdx);
}
else //mole-fraction formulation
{
globalSol[globalIdx][switchIdx]
= volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
}
}
else if (volVars.saturation(wPhaseIdx) <= Smin)
{
......@@ -644,8 +677,16 @@ public:
<< volVars.saturation(wPhaseIdx) << std::endl;
newPhasePresence = nPhaseOnly;
globalSol[globalIdx][switchIdx]
= volVars.fluidState().massFraction(nPhaseIdx, wCompIdx);
if(!useMoles) //mass-fraction formulation
{
globalSol[globalIdx][switchIdx]
= volVars.fluidState().massFraction(nPhaseIdx, wCompIdx);
}
else //mole-fraction formulation
{
globalSol[globalIdx][switchIdx]
= volVars.fluidState().moleFraction(nPhaseIdx, wCompIdx);
}
}
}
......
......@@ -62,6 +62,7 @@ NEW_PROP_TAG(MaterialLawParams); //!< The parameters of the material law (extrac
NEW_PROP_TAG(EffectiveDiffusivityModel); //!< The employed model for the computation of the effective diffusivity
NEW_PROP_TAG(ProblemEnableGravity); //!< Returns whether gravity is considered in the problem
NEW_PROP_TAG(UseMoles); //!Defines whether mole (true) or mass (false) fractions are used
NEW_PROP_TAG(ImplicitMassUpwindWeight); //!< The value of the upwind weight for the mass conservation equations
NEW_PROP_TAG(ImplicitMobilityUpwindWeight); //!< Weight for the upwind mobility in the velocity calculation
NEW_PROP_TAG(ReplaceCompEqIdx); //!< The index of the total mass balance equation, if one component balance is replaced (ReplaceCompEqIdx < NumComponents)
......
......@@ -158,6 +158,9 @@ SET_BOOL_PROP(TwoPTwoC, VtkAddVelocity, false);
// enable gravity by default
SET_BOOL_PROP(TwoPTwoC, ProblemEnableGravity, true);
SET_BOOL_PROP(TwoPTwoC, UseMoles, true); //!< Define that mole fractions are used in the balance equations per default
//! default value for the forchheimer coefficient
// Source: Ward, J.C. 1964 Turbulent flow in porous media. ASCE J. Hydraul. Div 90.
// Actually the Forchheimer coefficient is also a function of the dimensions of the
......
......@@ -99,6 +99,7 @@ class TwoPTwoCVolumeVariables : public ImplicitVolumeVariables<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef Dumux::MiscibleMultiPhaseComposition<Scalar, FluidSystem> MiscibleMultiPhaseComposition;
static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
typedef Dumux::ComputeFromReferencePhase<Scalar, FluidSystem> ComputeFromReferencePhase;
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
......@@ -249,21 +250,34 @@ public:
// only the nonwetting phase is present, i.e. nonwetting phase
// composition is stored explicitly.
// extract _mass_ fractions in the nonwetting phase
Scalar massFractionN[numComponents];
massFractionN[wCompIdx] = priVars[switchIdx];
massFractionN[nCompIdx] = 1 - massFractionN[wCompIdx];
// calculate average molar mass of the nonwetting phase
Scalar M1 = FluidSystem::molarMass(wCompIdx);
Scalar M2 = FluidSystem::molarMass(nCompIdx);
Scalar X2 = massFractionN[nCompIdx];
Scalar avgMolarMass = M1*M2/(M2 + X2*(M1 - M2));
// convert mass to mole fractions and set the fluid state
fluidState.setMoleFraction(nPhaseIdx, wCompIdx, massFractionN[wCompIdx]*avgMolarMass/M1);
fluidState.setMoleFraction(nPhaseIdx, nCompIdx, massFractionN[nCompIdx]*avgMolarMass/M2);
if(!useMoles) //mass-fraction formulation
{
// extract _mass_ fractions in the nonwetting phase
Scalar massFractionN[numComponents];
massFractionN[wCompIdx] = priVars[switchIdx];
massFractionN[nCompIdx] = 1 - massFractionN[wCompIdx];
// calculate average molar mass of the nonwetting phase
Scalar M1 = FluidSystem::molarMass(wCompIdx);
Scalar M2 = FluidSystem::molarMass(nCompIdx);
Scalar X2 = massFractionN[nCompIdx];
Scalar avgMolarMass = M1*M2/(M2 + X2*(M1 - M2));
// convert mass to mole fractions and set the fluid state
fluidState.setMoleFraction(nPhaseIdx, wCompIdx, massFractionN[wCompIdx]*avgMolarMass/M1);
fluidState.setMoleFraction(nPhaseIdx, nCompIdx, massFractionN[nCompIdx]*avgMolarMass/M2);
}
else //mole-fraction formulation
{
Scalar moleFractionN[numComponents];
moleFractionN[wCompIdx] = priVars[switchIdx];
moleFractionN[nCompIdx] = 1 - moleFractionN