From 75a46dccb3e17037436e49c000465ebcc758f766 Mon Sep 17 00:00:00 2001
From: Dominik Riesterer <thedom.89@googlemail.com>
Date: Fri, 20 Sep 2013 14:27:30 +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.
The model still uses MASS fractions by default.
The documentation is changed and extended accordingly.
reviewed by kathinka
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@11528 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/implicit/co2/co2model.hh | 27 ++++-
dumux/implicit/co2/co2volumevariables.hh | 122 +++++++++++++++-------
test/implicit/co2/heterogeneousproblem.hh | 39 ++++++-
3 files changed, 140 insertions(+), 48 deletions(-)
diff --git a/dumux/implicit/co2/co2model.hh b/dumux/implicit/co2/co2model.hh
index 9000e68b06..70ab5cb74c 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 6fc918cd44..ed6cfa7ac8 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 cac5b84920..31c58b9ad4 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
--
GitLab