Commit 75a46dcc 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.

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
parent 4908e266
...@@ -41,6 +41,8 @@ namespace Dumux ...@@ -41,6 +41,8 @@ namespace Dumux
* The CO2VolumeVariables do not use a constraint solver for calculating the mole fractions as is the * 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 * case in the 2p2c model. Instead mole fractions are calculated in the FluidSystem with a given
* temperature, pressurem and salinity. * 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> ...@@ -81,6 +83,7 @@ class CO2Model: public TwoPTwoCModel<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition; typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) }; enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
enum { dofCodim = isBox ? dim : 0 }; enum { dofCodim = isBox ? dim : 0 };
...@@ -228,8 +231,16 @@ public: ...@@ -228,8 +231,16 @@ public:
<< volVars.saturation(nPhaseIdx) << std::endl; << volVars.saturation(nPhaseIdx) << std::endl;
newPhasePresence = wPhaseOnly; newPhasePresence = wPhaseOnly;
globalSol[globalIdx][switchIdx] if(!useMoles) //mass-fraction formulation
= volVars.fluidState().massFraction(wPhaseIdx, nCompIdx); {
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) else if (volVars.saturation(wPhaseIdx) <= Smin)
{ {
...@@ -240,8 +251,16 @@ public: ...@@ -240,8 +251,16 @@ public:
<< volVars.saturation(wPhaseIdx) << std::endl; << volVars.saturation(wPhaseIdx) << std::endl;
newPhasePresence = nPhaseOnly; newPhasePresence = nPhaseOnly;
globalSol[globalIdx][switchIdx] if(!useMoles) //mass-fraction formulation
= volVars.fluidState().massFraction(nPhaseIdx, wCompIdx); {
globalSol[globalIdx][switchIdx]
= volVars.fluidState().massFraction(nPhaseIdx, wCompIdx);
}
else //mole-fraction formulation
{
globalSol[globalIdx][switchIdx]
= volVars.fluidState().moleFraction(nPhaseIdx, wCompIdx);
}
} }
} }
......
...@@ -95,6 +95,7 @@ class CO2VolumeVariables: public TwoPTwoCVolumeVariables<TypeTag> ...@@ -95,6 +95,7 @@ class CO2VolumeVariables: public TwoPTwoCVolumeVariables<TypeTag>
enum { dofCodim = isBox ? dim : 0 }; enum { dofCodim = isBox ? dim : 0 };
static const Scalar R; // universial nonwetting constant static const Scalar R; // universial nonwetting constant
static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
public: public:
//! The type of the object returned by the fluidState() method //! The type of the object returned by the fluidState() method
...@@ -217,50 +218,74 @@ public: ...@@ -217,50 +218,74 @@ public:
// only the nonwetting phase is present, i.e. nonwetting phase // only the nonwetting phase is present, i.e. nonwetting phase
// composition is stored explicitly. // 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 if(!useMoles) //mass-fraction formulation
Scalar M1 = FluidSystem::molarMass(wCompIdx); {
Scalar M2 = FluidSystem::molarMass(nCompIdx); // extract _mass_ fractions in the nonwetting phase
Scalar X2 = massFractionN[nCompIdx]; Scalar massFractionN[numComponents];
Scalar avgMolarMass = M1*M2/(M2 + X2*(M1 - M2)); massFractionN[wCompIdx] = priVars[switchIdx];
massFractionN[nCompIdx] = 1 - massFractionN[wCompIdx];
// convert mass to mole fractions and set the fluid state
ParentType::fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, massFractionN[wCompIdx]*avgMolarMass/M1); // calculate average molar mass of the nonwetting phase
ParentType::fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, massFractionN[nCompIdx]*avgMolarMass/M2); Scalar M1 = FluidSystem::molarMass(wCompIdx);
Scalar M2 = FluidSystem::molarMass(nCompIdx);
// TODO give values for non-existing wetting phase Scalar X2 = massFractionN[nCompIdx];
Scalar xwCO2 = FluidSystem::equilibriumMoleFraction(ParentType::fluidState_, paramCache, wPhaseIdx); Scalar avgMolarMass = M1*M2/(M2 + X2*(M1 - M2));
Scalar xwH2O = 1 - xwCO2;
// Scalar xwCO2 = FluidSystem::equilibriumMoleFraction(ParentType::fluidState_, paramCache, wPhaseIdx, nPhaseOnly); // convert mass to mole fractions and set the fluid state
// Scalar xwH2O = 1 - xwCO2; ParentType::fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, massFractionN[wCompIdx]*avgMolarMass/M1);
ParentType::fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwCO2); ParentType::fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, massFractionN[nCompIdx]*avgMolarMass/M2);
ParentType::fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xwH2O);
// TODO give values for non-existing wetting phase
Scalar xwCO2 = FluidSystem::equilibriumMoleFraction(ParentType::fluidState_, paramCache, wPhaseIdx);
//Get the phase densities from the FluidSystem and set them in the fluidState Scalar xwH2O = 1 - xwCO2;
// Scalar xwCO2 = FluidSystem::equilibriumMoleFraction(ParentType::fluidState_, paramCache, wPhaseIdx, nPhaseOnly);
Scalar rhoW = FluidSystem::density(ParentType::fluidState_, paramCache, wPhaseIdx); // Scalar xwH2O = 1 - xwCO2;
Scalar rhoN = FluidSystem::density(ParentType::fluidState_, paramCache, nPhaseIdx); ParentType::fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwCO2);
ParentType::fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xwH2O);
ParentType::fluidState_.setDensity(wPhaseIdx, rhoW); }
ParentType::fluidState_.setDensity(nPhaseIdx, rhoN); else //mole-fraction formulation
{
//Get the phase viscosities from the FluidSystem and set them in the fluidState // extract _mole_ fractions in the nonwetting phase
Scalar moleFractionN[numComponents];
Scalar muW = FluidSystem::viscosity(ParentType::fluidState_, paramCache, wPhaseIdx); moleFractionN[wCompIdx] = priVars[switchIdx];
Scalar muN = FluidSystem::viscosity(ParentType::fluidState_, paramCache, nPhaseIdx); moleFractionN[nCompIdx] = 1 - moleFractionN[wCompIdx];
ParentType::fluidState_.setViscosity(wPhaseIdx, muW); // set the fluid state
ParentType::fluidState_.setViscosity(nPhaseIdx, muN); 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) { else if (phasePresence == wPhaseOnly) {
// only the wetting phase is present, i.e. wetting phase // only the wetting phase is present, i.e. wetting phase
// composition is stored explicitly. // composition is stored explicitly.
if(!useMoles) //mass-fraction formulation
{
// extract _mass_ fractions in the nonwetting phase // extract _mass_ fractions in the nonwetting phase
Scalar massFractionW[numComponents]; Scalar massFractionW[numComponents];
massFractionW[nCompIdx] = priVars[switchIdx]; massFractionW[nCompIdx] = priVars[switchIdx];
...@@ -283,7 +308,26 @@ public: ...@@ -283,7 +308,26 @@ public:
// Scalar xnCO2 = 1 - xnH2O; // Scalar xnCO2 = 1 - xnH2O;
ParentType::fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnCO2); ParentType::fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnCO2);
ParentType::fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnH2O); 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 rhoW = FluidSystem::density(ParentType::fluidState_, paramCache, wPhaseIdx);
Scalar rhoN = FluidSystem::density(ParentType::fluidState_, paramCache, nPhaseIdx); Scalar rhoN = FluidSystem::density(ParentType::fluidState_, paramCache, nPhaseIdx);
......
...@@ -101,6 +101,8 @@ SET_BOOL_PROP(HeterogeneousProblem, ProblemEnableGravity, true); ...@@ -101,6 +101,8 @@ SET_BOOL_PROP(HeterogeneousProblem, ProblemEnableGravity, true);
SET_BOOL_PROP(HeterogeneousProblem, ImplicitEnableJacobianRecycling, false); SET_BOOL_PROP(HeterogeneousProblem, ImplicitEnableJacobianRecycling, false);
SET_BOOL_PROP(HeterogeneousProblem, VtkAddVelocity, false); SET_BOOL_PROP(HeterogeneousProblem, VtkAddVelocity, false);
SET_BOOL_PROP(HeterogeneousProblem, UseMoles, false);
} }
...@@ -121,6 +123,9 @@ SET_BOOL_PROP(HeterogeneousProblem, VtkAddVelocity, false); ...@@ -121,6 +123,9 @@ SET_BOOL_PROP(HeterogeneousProblem, VtkAddVelocity, false);
* between different parts of the boundary. * between different parts of the boundary.
* These boundary ids can be imported into the problem where the boundary conditions can then be assigned accordingly. * 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. * 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): * 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> ...@@ -150,6 +155,8 @@ class HeterogeneousProblem : public ImplicitPorousMediaProblem<TypeTag>
lPhaseIdx = Indices::wPhaseIdx, lPhaseIdx = Indices::wPhaseIdx,
gPhaseIdx = Indices::nPhaseIdx, gPhaseIdx = Indices::nPhaseIdx,
wCompIdx = FluidSystem::wCompIdx,
nCompIdx = FluidSystem::nCompIdx,
BrineIdx = FluidSystem::BrineIdx, BrineIdx = FluidSystem::BrineIdx,
CO2Idx = FluidSystem::CO2Idx, CO2Idx = FluidSystem::CO2Idx,
...@@ -175,6 +182,7 @@ class HeterogeneousProblem : public ImplicitPorousMediaProblem<TypeTag> ...@@ -175,6 +182,7 @@ class HeterogeneousProblem : public ImplicitPorousMediaProblem<TypeTag>
typedef Dumux::CO2<Scalar, CO2Table> CO2; typedef Dumux::CO2<Scalar, CO2Table> CO2;
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) }; enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
enum { dofCodim = isBox ? dim : 0 }; enum { dofCodim = isBox ? dim : 0 };
static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
public: public:
/*! /*!
...@@ -236,6 +244,16 @@ public: ...@@ -236,6 +244,16 @@ public:
/*pmin=*/pressureLow_, /*pmin=*/pressureLow_,
/*pmax=*/pressureHigh_, /*pmax=*/pressureHigh_,
/*np=*/nPressure_); /*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: ...@@ -335,6 +353,9 @@ public:
* *
* \param values Stores the source values, acts as return value * \param values Stores the source values, acts as return value
* \param globalPos The global position * \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, void sourceAtPos(PrimaryVariables &values,
const GlobalPosition &globalPos) const const GlobalPosition &globalPos) const
...@@ -410,6 +431,9 @@ public: ...@@ -410,6 +431,9 @@ public:
* *
* For this method, the \a values parameter stores the mass flux * For this method, the \a values parameter stores the mass flux
* in normal direction of each phase. Negative values mean influx. * 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, void neumann(PrimaryVariables &values,
const Element &element, const Element &element,
...@@ -423,7 +447,7 @@ public: ...@@ -423,7 +447,7 @@ public:
values = 0; values = 0;
if (boundaryId == injectionBottom_) 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: ...@@ -478,11 +502,16 @@ private:
Scalar meanM = Scalar meanM =
FluidSystem::molarMass(BrineIdx)*moleFracLiquidBrine + FluidSystem::molarMass(BrineIdx)*moleFracLiquidBrine +
FluidSystem::molarMass(CO2Idx)*moleFracLiquidCO2; FluidSystem::molarMass(CO2Idx)*moleFracLiquidCO2;
if(!useMoles) //mass-fraction formulation
Scalar massFracLiquidCO2 = moleFracLiquidCO2*FluidSystem::molarMass(CO2Idx)/meanM; {
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::pressureIdx] = pl;
values[Indices::switchIdx] = massFracLiquidCO2;
} }
Scalar temperature_(const GlobalPosition globalPos) const Scalar temperature_(const GlobalPosition globalPos) const
......
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