Skip to content
GitLab
Commit 217df792 authored by Timo Koch's avatar Timo Koch
Browse files

[co2] Port the 2p2c co2 model

parent 65bdbd31
Expand all Hide whitespace changes
Inline Side-by-side
dumux/porousmediumflow/2p2c/implicit/volumevariables.hh
View file @ 217df792
......@@ -110,9 +110,13 @@ class TwoPTwoCVolumeVariables : public ImplicitVolumeVariables<TypeTag>
enum { dofCodim = isBox ? dim : 0 };
public:
//! The type of the object returned by the fluidState() method
using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
//! Pull member functions of the parent for deriving classes
using ParentType::temperature;
using ParentType::enthalpy;
/*!
* \copydoc ImplicitVolumeVariables::update
*/
......@@ -123,7 +127,7 @@ public:
{
ParentType::update(elemSol, problem, element, scv);
completeFluidState(elemSol, problem, element, scv, fluidState_);
Implementation::completeFluidState(elemSol, problem, element, scv, fluidState_);
/////////////
// calculate the remaining quantities
......@@ -162,7 +166,6 @@ public:
/*!
* \copydoc ImplicitModel::completeFluidState
* \param isOldSol Specifies whether this is the previous solution or the current one
*/
static void completeFluidState(const ElementSolution& elemSol,
const Problem& problem,
......@@ -287,16 +290,6 @@ public:
fluidState.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
fluidState.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
fluidState.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
paramCache.updateComposition(fluidState, wPhaseIdx);
paramCache.updateComposition(fluidState, nPhaseIdx);
//set the phase densities
Scalar rhoW = FluidSystem::density(fluidState, paramCache, wPhaseIdx);
Scalar rhoN = FluidSystem::density(fluidState, paramCache, nPhaseIdx);
fluidState.setDensity(wPhaseIdx, rhoW);
fluidState.setDensity(nPhaseIdx, rhoN);
}
}
else if (phasePresence == nPhaseOnly)
......@@ -354,15 +347,6 @@ public:
fluidState.setMoleFraction(wPhaseIdx, wCompIdx, xww);
fluidState.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
paramCache.updateComposition(fluidState, wPhaseIdx);
paramCache.updateComposition(fluidState, nPhaseIdx);
Scalar rhoW = FluidSystem::density(fluidState, paramCache, wPhaseIdx);
Scalar rhoN = FluidSystem::density(fluidState, paramCache, nPhaseIdx);
fluidState.setDensity(wPhaseIdx, rhoW);
fluidState.setDensity(nPhaseIdx, rhoN);
}
}
else if (phasePresence == wPhaseOnly)
......@@ -416,23 +400,17 @@ public:
fluidState.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
fluidState.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
paramCache.updateComposition(fluidState, wPhaseIdx);
paramCache.updateComposition(fluidState, nPhaseIdx);
Scalar rhoW = FluidSystem::density(fluidState, paramCache, wPhaseIdx);
Scalar rhoN = FluidSystem::density(fluidState, paramCache, nPhaseIdx);
fluidState.setDensity(wPhaseIdx, rhoW);
fluidState.setDensity(nPhaseIdx, rhoN);
}
}
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
//set the viscosity here if constraintsolver is not used
// set the viscosity and desity here if constraintsolver is not used
if(!useConstraintSolver)
{
paramCache.updateComposition(fluidState, phaseIdx);
const Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
fluidState.setDensity(phaseIdx, rho);
const Scalar mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
fluidState.setViscosity(phaseIdx,mu);
}
......
dumux/porousmediumflow/co2/implicit/model.hh
View file @ 217df792
......@@ -21,17 +21,20 @@
*
* \brief Adaption of the fully implicit scheme to the CO2Model model.
*/
#ifndef DUMUX_CO2_MODEL_HH
#define DUMUX_CO2_MODEL_HH
#ifndef DUMUX_TWOP_TWOC_CO2_MODEL_HH
#define DUMUX_TWOP_TWOC_CO2_MODEL_HH
#include <dumux/common/properties.hh>
#include <dumux/porousmediumflow/2p2c/implicit/model.hh>
#include "primaryvariableswitch.hh"
#include "volumevariables.hh"
namespace Dumux
{
/*!
* \file
* \ingroup CO2Model
* \brief Adaption of the BOX or CC scheme to the non-isothermal two-phase two-component flow model.
* \brief Adaption of the non-isothermal two-phase two-component flow model to problems with CO2
*
* TODO: Put a doxgyen link refernce here
* See TwoPTwoCModel for reference to the equations used.
* The CO2 model is derived from the 2p2c model. In the CO2 model the phase switch criterion
* is different from the 2p2c model.
......@@ -45,253 +48,19 @@ namespace Dumux
* problem file. Make sure that the according units are used in the problem setup. useMoles is set to false by default.
*
*/
namespace Dumux {
namespace Properties {
template<class TypeTag>
class CO2Model: public TwoPTwoCModel<TypeTag>
{
typedef TwoPTwoCModel<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
enum {
switchIdx = Indices::switchIdx,
wPhaseIdx = Indices::wPhaseIdx,
nPhaseIdx = Indices::nPhaseIdx,
wCompIdx = Indices::wCompIdx,
nCompIdx = Indices::nCompIdx,
wPhaseOnly = Indices::wPhaseOnly,
nPhaseOnly = Indices::nPhaseOnly,
bothPhases = Indices::bothPhases,
pwsn = TwoPTwoCFormulation::pwsn,
pnsw = TwoPTwoCFormulation::pnsw,
formulation = GET_PROP_VALUE(TypeTag, Formulation)
};
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
enum {
dim = GridView::dimension,
dimWorld = GridView::dimensionworld
};
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 };
public:
/*!
* \brief Update the static data of all vertices in the grid.
*
* \param curGlobalSol The current global solution
* \param oldGlobalSol The previous global solution
*/
void updateStaticData(SolutionVector &curGlobalSol,
const SolutionVector &oldGlobalSol)
{
bool wasSwitched = false;
int succeeded;
try {
for (unsigned i = 0; i < ParentType::staticDat_.size(); ++i)
ParentType::staticDat_[i].visited = false;
FVElementGeometry fvGeometry;
static VolumeVariables volVars;
for (const auto& element : elements(this->gridView_()))
{
fvGeometry.update(this->gridView_(), element);
for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
{
int dofIdxGlobal = this->dofMapper().subIndex(element, scvIdx, dofCodim);
if (ParentType::staticDat_[dofIdxGlobal].visited)
continue;
ParentType::staticDat_[dofIdxGlobal].visited = true;
volVars.update(curGlobalSol[dofIdxGlobal],
this->problem_(),
element,
fvGeometry,
scvIdx,
false);
const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
if (primaryVarSwitch_(curGlobalSol,
volVars,
dofIdxGlobal,
globalPos))
{
this->jacobianAssembler().markDofRed(dofIdxGlobal);
wasSwitched = true;
}
}
}
succeeded = 1;
}
catch (NumericalProblem &e)
{
std::cout << "\n"
<< "Rank " << this->problem_().gridView().comm().rank()
<< " caught an exception while updating the static data." << e.what()
<< "\n";
succeeded = 0;
}
//make sure that all processes succeeded. If not throw a NumericalProblem to decrease the time step size.
if (this->gridView_().comm().size() > 1)
succeeded = this->gridView_().comm().min(succeeded);
if (!succeeded) {
DUNE_THROW(NumericalProblem,
"A process did not succeed in updating the static data.");
return;
}
// make sure that if there was a variable switch in an
// other partition we will also set the switch flag
// for our partition.
if (this->gridView_().comm().size() > 1)
wasSwitched = this->gridView_().comm().max(wasSwitched);
ParentType::setSwitched_(wasSwitched);
}
protected:
/*!
* \brief Performs variable switch at a vertex, returns true if a
* variable switch was performed.
*/
bool primaryVarSwitch_(SolutionVector &globalSol,
const VolumeVariables &volVars, int dofIdxGlobal,
const GlobalPosition &globalPos)
{
typename FluidSystem::ParameterCache paramCache;
// evaluate primary variable switch
bool wouldSwitch = false;
int phasePresence = ParentType::staticDat_[dofIdxGlobal].phasePresence;
int newPhasePresence = phasePresence;
// check if a primary var switch is necessary
if (phasePresence == nPhaseOnly)
{
Scalar xnw = volVars.moleFraction(nPhaseIdx, wCompIdx);
Scalar xnwMax = FluidSystem::equilibriumMoleFraction(volVars.fluidState(), paramCache, nPhaseIdx);
if(xnw > xnwMax)
wouldSwitch = true;
if (ParentType::staticDat_[dofIdxGlobal].wasSwitched)
xnwMax *= 1.02;
//If mole fraction is higher than the equilibrium mole fraction make a phase switch
if(xnw > xnwMax)
{
// wetting phase appears
std::cout << "wetting phase appears at vertex " << dofIdxGlobal
<< ", coordinates: " << globalPos << ", xnw > xnwMax: "
<< xnw << " > "<< xnwMax << std::endl;
newPhasePresence = bothPhases;
if (formulation == pnsw)
globalSol[dofIdxGlobal][switchIdx] = 0.0;
else if (formulation == pwsn)
globalSol[dofIdxGlobal][switchIdx] = 1.0;
}
}
else if (phasePresence == wPhaseOnly)
{
Scalar xwn = volVars.moleFraction(wPhaseIdx, nCompIdx);
Scalar xwnMax = FluidSystem::equilibriumMoleFraction(volVars.fluidState(), paramCache, wPhaseIdx);
//If mole fraction is higher than the equilibrium mole fraction make a phase switch
if(xwn > xwnMax)
wouldSwitch = true;
if (ParentType::staticDat_[dofIdxGlobal].wasSwitched)
xwnMax *= 1.02;
if(xwn > xwnMax)
{
// non-wetting phase appears
std::cout << "non-wetting phase appears at vertex " << dofIdxGlobal
<< ", coordinates: " << globalPos << ", xwn > xwnMax: "
<< xwn << " > "<< xwnMax << std::endl;
newPhasePresence = bothPhases;
if (formulation == pnsw)
globalSol[dofIdxGlobal][switchIdx] = 0.999;
else if (formulation == pwsn)
globalSol[dofIdxGlobal][switchIdx] = 0.001;
}
}
else if (phasePresence == bothPhases)
{
Scalar Smin = 0.0;
if (ParentType::staticDat_[dofIdxGlobal].wasSwitched)
Smin = -0.01;
if (volVars.saturation(nPhaseIdx) <= Smin)
{
wouldSwitch = true;
// nonwetting phase disappears
std::cout << "Nonwetting phase disappears at vertex " << dofIdxGlobal
<< ", coordinates: " << globalPos << ", sn: "
<< volVars.saturation(nPhaseIdx) << std::endl;
newPhasePresence = wPhaseOnly;
if(useMoles) // mole-fraction formulation
{
globalSol[dofIdxGlobal][switchIdx]
= volVars.moleFraction(wPhaseIdx, nCompIdx);
}
else // mass-fraction formulation
{
globalSol[dofIdxGlobal][switchIdx]
= volVars.massFraction(wPhaseIdx, nCompIdx);
}
}
else if (volVars.saturation(wPhaseIdx) <= Smin)
{
wouldSwitch = true;
// wetting phase disappears
std::cout << "Wetting phase disappears at vertex " << dofIdxGlobal
<< ", coordinates: " << globalPos << ", sw: "
<< volVars.saturation(wPhaseIdx) << std::endl;
newPhasePresence = nPhaseOnly;
if(useMoles) // mole-fraction formulation
{
globalSol[dofIdxGlobal][switchIdx]
= volVars.moleFraction(nPhaseIdx, wCompIdx);
}
else //mass-fraction formulation
{
globalSol[dofIdxGlobal][switchIdx]
= volVars.massFraction(nPhaseIdx, wCompIdx);
}
}
}
ParentType::staticDat_[dofIdxGlobal].phasePresence = newPhasePresence;
ParentType::staticDat_[dofIdxGlobal].wasSwitched = wouldSwitch;
return phasePresence != newPhasePresence;
}
NEW_TYPE_TAG(TwoPTwoCCO2, INHERITS_FROM(TwoPTwoC));
NEW_TYPE_TAG(TwoPTwoCCO2NI, INHERITS_FROM(TwoPTwoCNI));
};
//! the CO2 privarswitch and VolumeVariables properties
SET_TYPE_PROP(TwoPTwoCCO2, PrimaryVariableSwitch, TwoPTwoCCO2PrimaryVariableSwitch<TypeTag>);
SET_TYPE_PROP(TwoPTwoCCO2NI, PrimaryVariableSwitch, TwoPTwoCCO2PrimaryVariableSwitch<TypeTag>);
SET_TYPE_PROP(TwoPTwoCCO2, VolumeVariables, TwoPTwoCCO2VolumeVariables<TypeTag>);
SET_TYPE_PROP(TwoPTwoCCO2NI, VolumeVariables, TwoPTwoCCO2VolumeVariables<TypeTag>);
}
} // end namespace Properties
} // end namespace Dumux
#endif
dumux/porousmediumflow/co2/implicit/primaryvariableswitch.hh 0 → 100644
View file @ 217df792
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
*
* \brief The primary variable switch for the 2p2c-CO2 model
*/
#ifndef DUMUX_2P2C_CO2_PRIMARY_VARIABLE_SWITCH_HH
#define DUMUX_2P2C_CO2_PRIMARY_VARIABLE_SWITCH_HH
#include <dumux/porousmediumflow/compositional/primaryvariableswitch.hh>
namespace Dumux
{
/*!
* \ingroup TwoPTwoCModel
* \brief The primary variable switch for the 2p2c-CO2 model controlling the phase presence state variable
* The phase switch occurs when the equilibrium concentration
* of a component in a phase is exceeded, instead of the sum of the components in the virtual phase
* (the phase which is not present) being greater that unity as done in the 2p2c model.
*/
template<class TypeTag>
class TwoPTwoCCO2PrimaryVariableSwitch : public PrimaryVariableSwitch<TypeTag>
{
friend typename Dumux::PrimaryVariableSwitch<TypeTag>;
using ParentType = Dumux::PrimaryVariableSwitch<TypeTag>;
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using IndexType = typename GridView::IndexSet::IndexType;
using GlobalPosition = Dune::FieldVector<Scalar, GridView::dimensionworld>;
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
enum {
switchIdx = Indices::switchIdx,
wPhaseIdx = Indices::wPhaseIdx,
nPhaseIdx = Indices::nPhaseIdx,
wCompIdx = Indices::wCompIdx,
nCompIdx = Indices::nCompIdx,
wPhaseOnly = Indices::wPhaseOnly,
nPhaseOnly = Indices::nPhaseOnly,
bothPhases = Indices::bothPhases,
pwsn = TwoPTwoCFormulation::pwsn,
pnsw = TwoPTwoCFormulation::pnsw,
formulation = GET_PROP_VALUE(TypeTag, Formulation)
};
static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
public:
using ParentType::ParentType;
protected:
// perform variable switch at a degree of freedom location
bool update_(PrimaryVariables& priVars,
const VolumeVariables& volVars,
IndexType dofIdxGlobal,
const GlobalPosition& globalPos)
{
// evaluate primary variable switch
bool wouldSwitch = false;
int phasePresence = priVars.state();
int newPhasePresence = phasePresence;
// the param cache to evaluate the equilibrium mole fraction
typename FluidSystem::ParameterCache paramCache;
// check if a primary var switch is necessary
if (phasePresence == nPhaseOnly)
{
// calculate wetting component mole fraction in the non-wetting phase
Scalar xnw = volVars.moleFraction(nPhaseIdx, wCompIdx);
Scalar xnwMax = FluidSystem::equilibriumMoleFraction(volVars.fluidState(), paramCache, nPhaseIdx);
// if it is larger than the equilibirum mole fraction switch
if(xnw > xnwMax)
wouldSwitch = true;
if (this->wasSwitched_[dofIdxGlobal])
xnwMax *= 1.02;
// if it is larger than the equilibirum mole fraction switch wetting phase appears
if (xnw > xnwMax)
{
// wetting phase appears
std::cout << "wetting phase appears at vertex " << dofIdxGlobal
<< ", coordinates: " << globalPos << ", xnw > xnwMax: "
<< xnw << " > " << xnwMax << std::endl;
newPhasePresence = bothPhases;
if (formulation == pnsw)
priVars[switchIdx] = 0.0;
else if (formulation == pwsn)
priVars[switchIdx] = 1.0;
}
}
else if (phasePresence == wPhaseOnly)
{
// calculate non-wetting component mole fraction in the wetting phase
Scalar xwn = volVars.moleFraction(wPhaseIdx, nCompIdx);
Scalar xwnMax = FluidSystem::equilibriumMoleFraction(volVars.fluidState(), paramCache, wPhaseIdx);
// if it is larger than the equilibirum mole fraction switch
if(xwn > xwnMax)
wouldSwitch = true;
if (this->wasSwitched_[dofIdxGlobal])
xwnMax *= 1.02;
// if it is larger than the equilibirum mole fraction switch non-wetting phase appears
if(xwn > xwnMax)
{
// nonwetting phase appears
std::cout << "nonwetting phase appears at vertex " << dofIdxGlobal
<< ", coordinates: " << globalPos << ", xwn > xwnMax: "
<< xwn << " > " << xwnMax << std::endl;
newPhasePresence = bothPhases;
if (formulation == pnsw)
priVars[switchIdx] = 0.999;
else if (formulation == pwsn)
priVars[switchIdx] = 0.001;
}
}
// TODO: this is the same as for the 2p2c model maybe factor out
else if (phasePresence == bothPhases)
{
Scalar Smin = 0.0;
if (this->wasSwitched_[dofIdxGlobal])
Smin = -0.01;
if (volVars.saturation(nPhaseIdx) <= Smin)
{
wouldSwitch = true;
// nonwetting phase disappears
std::cout << "Nonwetting phase disappears at vertex " << dofIdxGlobal
<< ", coordinates: " << globalPos << ", sn: "
<< volVars.saturation(nPhaseIdx) << std::endl;
newPhasePresence = wPhaseOnly;
if(useMoles) // mole-fraction formulation
{
priVars[switchIdx]
= volVars.moleFraction(wPhaseIdx, nCompIdx);
}
else // mass-fraction formulation
{
priVars[switchIdx]
= volVars.massFraction(wPhaseIdx, nCompIdx);
}
}
else if (volVars.saturation(wPhaseIdx) <= Smin)
{
wouldSwitch = true;
// wetting phase disappears
std::cout << "Wetting phase disappears at vertex " << dofIdxGlobal
<< ", coordinates: " << globalPos << ", sw: "
<< volVars.saturation(wPhaseIdx) << std::endl;
newPhasePresence = nPhaseOnly;
if(useMoles) // mole-fraction formulation
{
priVars[switchIdx]
= volVars.moleFraction(nPhaseIdx, wCompIdx);
}
else // mass-fraction formulation
{
priVars[switchIdx]
= volVars.massFraction(nPhaseIdx, wCompIdx);
}
}
}
priVars.setState(newPhasePresence);
this->wasSwitched_[dofIdxGlobal] = wouldSwitch;
return phasePresence != newPhasePresence;
}
};
} // end namespace Dumux
#endif
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