From 10237266818f9438213a7c411e68dd76b921f26c Mon Sep 17 00:00:00 2001
From: Dominik Riesterer <thedom.89@googlemail.com>
Date: Fri, 20 Sep 2013 14:25:14 +0000
Subject: [PATCH] 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
---
dumux/implicit/2p2c/2p2clocalresidual.hh | 367 ++++++++++++++------
dumux/implicit/2p2c/2p2cmodel.hh | 55 ++-
dumux/implicit/2p2c/2p2cproperties.hh | 1 +
dumux/implicit/2p2c/2p2cpropertydefaults.hh | 3 +
dumux/implicit/2p2c/2p2cvolumevariables.hh | 83 +++--
test/implicit/2p2c/injectionproblem.hh | 38 +-
6 files changed, 393 insertions(+), 154 deletions(-)
diff --git a/dumux/implicit/2p2c/2p2clocalresidual.hh b/dumux/implicit/2p2c/2p2clocalresidual.hh
index 08f1101126..9e58d4ae6e 100644
--- a/dumux/implicit/2p2c/2p2clocalresidual.hh
+++ b/dumux/implicit/2p2c/2p2clocalresidual.hh
@@ -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;
+ }
+ }
}
/*!
diff --git a/dumux/implicit/2p2c/2p2cmodel.hh b/dumux/implicit/2p2c/2p2cmodel.hh
index 677b1ce17f..6cf1fa1e65 100644
--- a/dumux/implicit/2p2c/2p2cmodel.hh
+++ b/dumux/implicit/2p2c/2p2cmodel.hh
@@ -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);
+ }
}
}
diff --git a/dumux/implicit/2p2c/2p2cproperties.hh b/dumux/implicit/2p2c/2p2cproperties.hh
index 1c87bd9ea8..98b26ad907 100644
--- a/dumux/implicit/2p2c/2p2cproperties.hh
+++ b/dumux/implicit/2p2c/2p2cproperties.hh
@@ -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)
diff --git a/dumux/implicit/2p2c/2p2cpropertydefaults.hh b/dumux/implicit/2p2c/2p2cpropertydefaults.hh
index 15814127ab..a5b3e5509f 100644
--- a/dumux/implicit/2p2c/2p2cpropertydefaults.hh
+++ b/dumux/implicit/2p2c/2p2cpropertydefaults.hh
@@ -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
diff --git a/dumux/implicit/2p2c/2p2cvolumevariables.hh b/dumux/implicit/2p2c/2p2cvolumevariables.hh
index 537090802e..dbb34041ae 100644
--- a/dumux/implicit/2p2c/2p2cvolumevariables.hh
+++ b/dumux/implicit/2p2c/2p2cvolumevariables.hh
@@ -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[wCompIdx];
+
+ // set the fluid state
+ fluidState.setMoleFraction(nPhaseIdx, wCompIdx, moleFractionN[wCompIdx]);
+ fluidState.setMoleFraction(nPhaseIdx, nCompIdx, moleFractionN[nCompIdx]);
+ }
// calculate the composition of the remaining phases (as
// well as the densities of all phases). this is the job
// of the "ComputeFromReferencePhase" constraint solver
@@ -277,21 +291,34 @@ public:
// only the wetting phase is present, i.e. wetting phase
// composition is stored explicitly.
- // extract _mass_ fractions in the nonwetting phase
- Scalar massFractionW[numComponents];
- massFractionW[nCompIdx] = priVars[switchIdx];
- massFractionW[wCompIdx] = 1 - massFractionW[nCompIdx];
-
- // calculate average molar mass of the nonwetting phase
- Scalar M1 = FluidSystem::molarMass(wCompIdx);
- Scalar M2 = FluidSystem::molarMass(nCompIdx);
- Scalar X2 = massFractionW[nCompIdx];
- Scalar avgMolarMass = M1*M2/(M2 + X2*(M1 - M2));
-
- // convert mass to mole fractions and set the fluid state
- fluidState.setMoleFraction(wPhaseIdx, wCompIdx, massFractionW[wCompIdx]*avgMolarMass/M1);
- fluidState.setMoleFraction(wPhaseIdx, nCompIdx, massFractionW[nCompIdx]*avgMolarMass/M2);
+ if(!useMoles) //mass-fraction formulation
+ {
+ // extract _mass_ fractions in the nonwetting phase
+ Scalar massFractionW[numComponents];
+ massFractionW[nCompIdx] = priVars[switchIdx];
+ massFractionW[wCompIdx] = 1 - massFractionW[nCompIdx];
+
+ // calculate average molar mass of the nonwetting phase
+ Scalar M1 = FluidSystem::molarMass(wCompIdx);
+ Scalar M2 = FluidSystem::molarMass(nCompIdx);
+ Scalar X2 = massFractionW[nCompIdx];
+ Scalar avgMolarMass = M1*M2/(M2 + X2*(M1 - M2));
+
+ // convert mass to mole fractions and set the fluid state
+ fluidState.setMoleFraction(wPhaseIdx, wCompIdx, massFractionW[wCompIdx]*avgMolarMass/M1);
+ fluidState.setMoleFraction(wPhaseIdx, nCompIdx, massFractionW[nCompIdx]*avgMolarMass/M2);
+ }
+ else //mole-fraction formulation
+ {
+ Scalar moleFractionW[numComponents];
+ moleFractionW[nCompIdx] = priVars[switchIdx];
+ moleFractionW[wCompIdx] = 1 - moleFractionW[nCompIdx];
+
+ // set the fluid state
+ fluidState.setMoleFraction(wPhaseIdx, wCompIdx, moleFractionW[wCompIdx]);
+ fluidState.setMoleFraction(wPhaseIdx, nCompIdx, moleFractionW[nCompIdx]);
+ }
// calculate the composition of the remaining phases (as
// well as the densities of all phases). this is the job
// of the "ComputeFromReferencePhase" constraint solver
diff --git a/test/implicit/2p2c/injectionproblem.hh b/test/implicit/2p2c/injectionproblem.hh
index c12bc1934a..75b43f816f 100644
--- a/test/implicit/2p2c/injectionproblem.hh
+++ b/test/implicit/2p2c/injectionproblem.hh
@@ -68,6 +68,9 @@ SET_PROP(InjectionProblem, FluidSystem)
// Enable gravity
SET_BOOL_PROP(InjectionProblem, ProblemEnableGravity, true);
+// Define whether mole(true) or mass (false) fractions are used
+SET_BOOL_PROP(InjectionProblem, UseMoles, true);
+
SET_BOOL_PROP(InjectionProblem, ImplicitEnableJacobianRecycling, true);
SET_BOOL_PROP(InjectionProblem, VtkAddVelocity, false);
}
@@ -86,6 +89,9 @@ SET_BOOL_PROP(InjectionProblem, VtkAddVelocity, false);
* (\f$ 5m<y<15m\f$) and migrates upwards due to buoyancy. It accumulates and
* partially enters the lower permeable aquitard.
*
+ * 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. The default setting for useMoles is true.
+ *
* This problem uses the \ref TwoPTwoCModel.
*
* To run the simulation execute the following line in shell:
@@ -137,6 +143,9 @@ class InjectionProblem : public ImplicitPorousMediaProblem<TypeTag>
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ //! property that defines whether mole or mass fractions are used
+ static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
+
public:
/*!
* \brief The constructor
@@ -190,6 +199,16 @@ public:
/*pmin=*/pressureLow_,
/*pmax=*/pressureHigh_,
/*np=*/nPressure_);
+
+ //stateing in the console whether mole or mass fractions are used
+ if(!useMoles)
+ {
+ std::cout<<"problem uses mass-fractions"<<std::endl;
+ }
+ else
+ {
+ std::cout<<"problem uses mole-fractions"<<std::endl;
+ }
}
/*!
@@ -287,6 +306,8 @@ public:
*
* For this method, the \a values parameter stores the mass flux
* in normal direction of each phase. Negative values mean influx.
+ *
+ * The units must be according to either using mole or mass fractions. (mole/(m^2*s) or kg/(m^2*s))
*/
void neumann(PrimaryVariables &values,
const Element &element,
@@ -304,7 +325,7 @@ public:
globalPos = intersection.geometry().center();
if (globalPos[1] < 15 && globalPos[1] > 7) {
- values[contiN2EqIdx] = -1e-3; // kg/(s*m^2)
+ values[contiN2EqIdx] = -1e-3/FluidSystem::molarMass(nCompIdx); //mole/(m^2*s) //-1e-3; // kg/(s*m^2)
}
}
@@ -357,11 +378,16 @@ private:
Scalar meanM =
FluidSystem::molarMass(wCompIdx)*moleFracLiquidH2O +
FluidSystem::molarMass(nCompIdx)*moleFracLiquidN2;
-
- Scalar massFracLiquidN2 = moleFracLiquidN2*FluidSystem::molarMass(nCompIdx)/meanM;
-
+ if(!useMoles) //mass fraction formulation
+ {
+ Scalar massFracLiquidN2 = moleFracLiquidN2*FluidSystem::molarMass(nCompIdx)/meanM;
+ values[Indices::switchIdx] = massFracLiquidN2;
+ }
+ else //mole-fraction formulation
+ {
+ values[Indices::switchIdx] = moleFracLiquidN2;
+ }
values[Indices::pressureIdx] = pl;
- values[Indices::switchIdx] = massFracLiquidN2;
}
Scalar temperature_;
@@ -377,8 +403,6 @@ private:
Scalar temperatureLow_, temperatureHigh_;
-
-
};
} //end namespace
--
GitLab