diff --git a/dumux/implicit/co2/co2model.hh b/dumux/implicit/co2/co2model.hh
index 9000e68b0672ab7e030ec4a4fe975e6818a4b3c9..70ab5cb74c26214c87ce121a7b8b6602d94a50b3 100644
--- a/dumux/implicit/co2/co2model.hh
+++ b/dumux/implicit/co2/co2model.hh
@@ -41,6 +41,8 @@ namespace Dumux
* The CO2VolumeVariables do not use a constraint solver for calculating the mole fractions as is the
* case in the 2p2c model. Instead mole fractions are calculated in the FluidSystem with a given
* temperature, pressurem and salinity.
+ * 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 false by default.
*
*/
@@ -81,6 +83,7 @@ class CO2Model: public TwoPTwoCModel<TypeTag>
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 };
@@ -228,8 +231,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)
{
@@ -240,8 +251,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/co2/co2volumevariables.hh b/dumux/implicit/co2/co2volumevariables.hh
index 6fc918cd4446a0454f2f9571409b68afad801c2c..ed6cfa7ac8cbd57c0f0782c005d8963aafedf162 100644
--- a/dumux/implicit/co2/co2volumevariables.hh
+++ b/dumux/implicit/co2/co2volumevariables.hh
@@ -95,6 +95,7 @@ class CO2VolumeVariables: public TwoPTwoCVolumeVariables<TypeTag>
enum { dofCodim = isBox ? dim : 0 };
static const Scalar R; // universial nonwetting constant
+ static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
public:
//! The type of the object returned by the fluidState() method
@@ -217,50 +218,74 @@ 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
- ParentType::fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, massFractionN[wCompIdx]*avgMolarMass/M1);
- ParentType::fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, massFractionN[nCompIdx]*avgMolarMass/M2);
-
- // TODO give values for non-existing wetting phase
- Scalar xwCO2 = FluidSystem::equilibriumMoleFraction(ParentType::fluidState_, paramCache, wPhaseIdx);
- Scalar xwH2O = 1 - xwCO2;
-// Scalar xwCO2 = FluidSystem::equilibriumMoleFraction(ParentType::fluidState_, paramCache, wPhaseIdx, nPhaseOnly);
-// Scalar xwH2O = 1 - xwCO2;
- ParentType::fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwCO2);
- ParentType::fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xwH2O);
-
-
- //Get the phase densities from the FluidSystem and set them in the fluidState
-
- Scalar rhoW = FluidSystem::density(ParentType::fluidState_, paramCache, wPhaseIdx);
- Scalar rhoN = FluidSystem::density(ParentType::fluidState_, paramCache, nPhaseIdx);
-
- ParentType::fluidState_.setDensity(wPhaseIdx, rhoW);
- ParentType::fluidState_.setDensity(nPhaseIdx, rhoN);
-
- //Get the phase viscosities from the FluidSystem and set them in the fluidState
-
- Scalar muW = FluidSystem::viscosity(ParentType::fluidState_, paramCache, wPhaseIdx);
- Scalar muN = FluidSystem::viscosity(ParentType::fluidState_, paramCache, nPhaseIdx);
-
- ParentType::fluidState_.setViscosity(wPhaseIdx, muW);
- ParentType::fluidState_.setViscosity(nPhaseIdx, muN);
+ 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
+ ParentType::fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, massFractionN[wCompIdx]*avgMolarMass/M1);
+ ParentType::fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, massFractionN[nCompIdx]*avgMolarMass/M2);
+
+ // TODO give values for non-existing wetting phase
+ Scalar xwCO2 = FluidSystem::equilibriumMoleFraction(ParentType::fluidState_, paramCache, wPhaseIdx);
+ Scalar xwH2O = 1 - xwCO2;
+ // Scalar xwCO2 = FluidSystem::equilibriumMoleFraction(ParentType::fluidState_, paramCache, wPhaseIdx, nPhaseOnly);
+ // Scalar xwH2O = 1 - xwCO2;
+ ParentType::fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwCO2);
+ ParentType::fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xwH2O);
+ }
+ else //mole-fraction formulation
+ {
+ // extract _mole_ fractions in the nonwetting phase
+ Scalar moleFractionN[numComponents];
+ moleFractionN[wCompIdx] = priVars[switchIdx];
+ moleFractionN[nCompIdx] = 1 - moleFractionN[wCompIdx];
+
+ // set the fluid state
+ ParentType::fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, moleFractionN[wCompIdx]);
+ ParentType::fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, moleFractionN[nCompIdx]);
+
+ // TODO give values for non-existing wetting phase
+ Scalar xwCO2 = FluidSystem::equilibriumMoleFraction(ParentType::fluidState_, paramCache, wPhaseIdx);
+ Scalar xwH2O = 1 - xwCO2;
+ // Scalar xwCO2 = FluidSystem::equilibriumMoleFraction(ParentType::fluidState_, paramCache, wPhaseIdx, nPhaseOnly);
+ // Scalar xwH2O = 1 - xwCO2;
+ ParentType::fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwCO2);
+ ParentType::fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xwH2O);
+ }
+
+ //Get the phase densities from the FluidSystem and set them in the fluidState
+
+ Scalar rhoW = FluidSystem::density(ParentType::fluidState_, paramCache, wPhaseIdx);
+ Scalar rhoN = FluidSystem::density(ParentType::fluidState_, paramCache, nPhaseIdx);
+
+ ParentType::fluidState_.setDensity(wPhaseIdx, rhoW);
+ ParentType::fluidState_.setDensity(nPhaseIdx, rhoN);
+
+ //Get the phase viscosities from the FluidSystem and set them in the fluidState
+
+ Scalar muW = FluidSystem::viscosity(ParentType::fluidState_, paramCache, wPhaseIdx);
+ Scalar muN = FluidSystem::viscosity(ParentType::fluidState_, paramCache, nPhaseIdx);
+
+ ParentType::fluidState_.setViscosity(wPhaseIdx, muW);
+ ParentType::fluidState_.setViscosity(nPhaseIdx, muN);
}
else if (phasePresence == wPhaseOnly) {
// only the wetting phase is present, i.e. wetting phase
// composition is stored explicitly.
+ if(!useMoles) //mass-fraction formulation
+ {
// extract _mass_ fractions in the nonwetting phase
Scalar massFractionW[numComponents];
massFractionW[nCompIdx] = priVars[switchIdx];
@@ -283,7 +308,26 @@ public:
// Scalar xnCO2 = 1 - xnH2O;
ParentType::fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnCO2);
ParentType::fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnH2O);
-
+ }
+ else //mole-fraction formulation
+ {
+ // extract _mole_ fractions in the nonwetting phase
+ Scalar moleFractionW[numComponents];
+ moleFractionW[nCompIdx] = priVars[switchIdx];
+ moleFractionW[wCompIdx] = 1 - moleFractionW[nCompIdx];
+
+ // convert mass to mole fractions and set the fluid state
+ ParentType::fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, moleFractionW[wCompIdx]);
+ ParentType::fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, moleFractionW[nCompIdx]);
+
+ // TODO give values for non-existing nonwetting phase
+ Scalar xnH2O = FluidSystem::equilibriumMoleFraction(ParentType::fluidState_, paramCache, nPhaseIdx);
+ Scalar xnCO2 = 1 - xnH2O; //FluidSystem::equilibriumMoleFraction(ParentType::fluidState_, paramCache, nPhaseIdx);
+ // Scalar xnH2O = FluidSystem::equilibriumMoleFraction(ParentType::fluidState_, paramCache, nPhaseIdx, wPhaseOnly);
+ // Scalar xnCO2 = 1 - xnH2O;
+ ParentType::fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnCO2);
+ ParentType::fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnH2O);
+ }
Scalar rhoW = FluidSystem::density(ParentType::fluidState_, paramCache, wPhaseIdx);
Scalar rhoN = FluidSystem::density(ParentType::fluidState_, paramCache, nPhaseIdx);
diff --git a/test/implicit/co2/heterogeneousproblem.hh b/test/implicit/co2/heterogeneousproblem.hh
index cac5b849201470e9ebb24b61113248fb463c2c45..31c58b9ad4ed3b02fbd43220d620a19efb96888b 100644
--- a/test/implicit/co2/heterogeneousproblem.hh
+++ b/test/implicit/co2/heterogeneousproblem.hh
@@ -101,6 +101,8 @@ SET_BOOL_PROP(HeterogeneousProblem, ProblemEnableGravity, true);
SET_BOOL_PROP(HeterogeneousProblem, ImplicitEnableJacobianRecycling, false);
SET_BOOL_PROP(HeterogeneousProblem, VtkAddVelocity, false);
+
+SET_BOOL_PROP(HeterogeneousProblem, UseMoles, false);
}
@@ -121,6 +123,9 @@ SET_BOOL_PROP(HeterogeneousProblem, VtkAddVelocity, false);
* between different parts of the boundary.
* These boundary ids can be imported into the problem where the boundary conditions can then be assigned accordingly.
*
+ * 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 false.
+ *
* This problem uses the \ref OnePTwoCBoxModel model.
*
* To run the simulation execute the following line in shell (works with the box and cell centered spatial discretization method):
@@ -150,6 +155,8 @@ class HeterogeneousProblem : public ImplicitPorousMediaProblem<TypeTag>
lPhaseIdx = Indices::wPhaseIdx,
gPhaseIdx = Indices::nPhaseIdx,
+ wCompIdx = FluidSystem::wCompIdx,
+ nCompIdx = FluidSystem::nCompIdx,
BrineIdx = FluidSystem::BrineIdx,
CO2Idx = FluidSystem::CO2Idx,
@@ -175,6 +182,7 @@ class HeterogeneousProblem : public ImplicitPorousMediaProblem<TypeTag>
typedef Dumux::CO2<Scalar, CO2Table> CO2;
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
enum { dofCodim = isBox ? dim : 0 };
+ static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
public:
/*!
@@ -236,6 +244,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;
+ }
}
/*!
@@ -335,6 +353,9 @@ public:
*
* \param values Stores the source values, acts as return value
* \param globalPos The global position
+ *
+ * Depending on whether useMoles is set on true or false, the flux has to be given either in
+ * kg/(m^3*s) or mole/(m^3*s) in the input file!!
*/
void sourceAtPos(PrimaryVariables &values,
const GlobalPosition &globalPos) const
@@ -410,6 +431,9 @@ public:
*
* For this method, the \a values parameter stores the mass flux
* in normal direction of each phase. Negative values mean influx.
+ *
+ * Depending on whether useMoles is set on true or false, the flux has to be given either in
+ * kg/(m^2*s) or mole/(m^2*s) in the input file!!
*/
void neumann(PrimaryVariables &values,
const Element &element,
@@ -423,7 +447,7 @@ public:
values = 0;
if (boundaryId == injectionBottom_)
{
- values[contiCO2EqIdx] = -injectionRate_; // kg/(s*m^2)
+ values[contiCO2EqIdx] = -injectionRate_; ///FluidSystem::molarMass(nCompIdx); // kg/(s*m^2) or mole/(m^2*s) !!
}
}
@@ -478,11 +502,16 @@ private:
Scalar meanM =
FluidSystem::molarMass(BrineIdx)*moleFracLiquidBrine +
FluidSystem::molarMass(CO2Idx)*moleFracLiquidCO2;
-
- Scalar massFracLiquidCO2 = moleFracLiquidCO2*FluidSystem::molarMass(CO2Idx)/meanM;
-
+ if(!useMoles) //mass-fraction formulation
+ {
+ Scalar massFracLiquidCO2 = moleFracLiquidCO2*FluidSystem::molarMass(CO2Idx)/meanM;
+ values[Indices::switchIdx] = massFracLiquidCO2;
+ }
+ else //mole-fraction formulation
+ {
+ values[Indices::switchIdx] = moleFracLiquidCO2;
+ }
values[Indices::pressureIdx] = pl;
- values[Indices::switchIdx] = massFracLiquidCO2;
}
Scalar temperature_(const GlobalPosition globalPos) const