Commit 8713e368 authored by Katharina Heck's avatar Katharina Heck
Browse files

unite misciblemultiphaseconstraintsolver and miscible2pncconstraintsolver and...

unite misciblemultiphaseconstraintsolver and miscible2pncconstraintsolver and deprecate unnecessary constraintsolvers. use misciblemultiphase now in all models
parent e84f0785
......@@ -62,94 +62,8 @@ namespace Dumux {
* - if the setEnthalpy parameter is true, also specific enthalpies and internal energies of *the reference* phase
*/
template <class Scalar, class FluidSystem>
class FluidSystemComputeFromReferencePhase
{
static constexpr int numPhases = FluidSystem::numPhases;
static constexpr int numComponents = FluidSystem::numComponents;
using ComponentVector = Dune::FieldVector<Scalar, numComponents>;
public:
/*!
* \brief Computes the composition of all phases of from a function in the fluidsystem.
*
* This constraint solver assumes that there is a function
* calculateEquilibriumMoleFractionOtherPhase
* in the fluidsystem. I.e. Either this function only has lookup tables or
* short-circuits the solution of a linear system of equations
* Therefore, this approach cannot be used if the solubility is a function
* of composition. For this nonlinear case the "fugacity coefficient"
* approach has to be used.
*
* The constraint solver assumes the following quantities to be set:
*
* - temperatures of *all* phases
* - saturations of *all* phases
* - pressures of *all* phases
*
* - the composition of the *reference* phase
*
* - the Fluidsystem has to be able to calculate the other mole fraction(s) from that information.
*
* After calling the solve() method the following quantities
* are calculated in addition:
*
* - density, molar density, molar volume of *the reference* phase
* - composition in mole and mass fractions and molarities of *all* phases
* - mean molar masses of *all* phases
* - if the setViscosity parameter is true, also dynamic viscosities of *the reference* phase
* - if the setEnthalpy parameter is true, also specific enthalpies and internal energies of *the reference* phase
*
* \param fluidState A container with the current (physical) state of the fluid.
* \param paramCache A container for iterative calculation of fluid composition
* \param refPhaseIdx The index of the phase whose's composition is known.
* \param setViscosity Should the viscosity be set in the fluidstate?
* \param setEnthalpy Should the enthalpy be set in the fluidstate?
*/
template <class FluidState, class ParameterCache>
static void solve(FluidState & fluidState,
ParameterCache & paramCache,
const unsigned int refPhaseIdx,
const bool setViscosity,
const bool setEnthalpy)
{
for(int compIdx=0; compIdx<numComponents; compIdx++){
// In this function the actual work is done.
// Either tables for solubility or functional relations are used therein.
// I give you phase and component and you set me the component in the other phase
FluidSystem::calculateEquilibriumMoleFractionOtherPhase(fluidState,
paramCache,
refPhaseIdx,
compIdx);
}
for(int phaseIdx=0; phaseIdx < numPhases; phaseIdx++){
Scalar value = FluidSystem::density(fluidState, paramCache, phaseIdx);
fluidState.setDensity(phaseIdx, value);
if (setViscosity) {
value = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
fluidState.setViscosity(phaseIdx, value);
}
if (setEnthalpy) {
value = FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
fluidState.setEnthalpy(phaseIdx, value);
}
checkDefinedMoleFractions(fluidState);
}
}
/*!
* \brief checks whether all the mole fractions which are stored in the fluidstate are founded on defined values.
*/
template <class FluidState>
static void checkDefinedMoleFractions(const FluidState & fluidState){
for(int phaseIdx=0; phaseIdx < numPhases; phaseIdx++)
for(int compIdx=0; compIdx<numComponents; compIdx++)
Valgrind::CheckDefined(fluidState.moleFraction(phaseIdx, compIdx));
// fluidState.checkDefined(); // cannot be fully defined for the case of viscosity and enthalpy not set
}
};
using FluidSystemComputeFromReferencePhase DUNE_DEPRECATED_MSG("Use MiscibleMultiPhaseComposition from dumux/material/constraintsolvers/misciblemultiphasecomposition.hh")
= MiscibleMultiPhaseComposition<Scalar, FluidSystem>;
} // end namespace Dumux
#endif
......@@ -31,6 +31,7 @@
#include <dumux/common/exceptions.hh>
#include <dumux/common/valgrind.hh>
#include <dumux/material/constraintsolvers/misciblemultiphasecomposition.hh>
namespace Dumux {
/*!
......@@ -58,152 +59,8 @@ namespace Dumux {
* - if the setEnthalpy parameter is true, also specific enthalpies of *all* phases
*/
template <class Scalar, class FluidSystem>
class Miscible2pNCComposition
{
static const int numPhases = FluidSystem::numPhases;
static const int numComponents = FluidSystem::numComponents;
static const int numMajorComponents = FluidSystem::numPhases;
public:
/*!
* \brief @copybrief Dumux::Miscible2pNCComposition
*/
template <class FluidState, class ParameterCache>
static void solve(FluidState &fluidState,
ParameterCache &paramCache,
const int knownPhaseIdx,
bool setViscosity,
bool setEnthalpy)
{
#ifndef NDEBUG
// currently this solver can only handle fluid systems which
// assume ideal mixtures of all fluids. TODO: relax this
// (requires solving a non-linear system of equations, i.e. using
// newton method.)
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
assert(FluidSystem::isIdealMixture(phaseIdx));
}
#endif
//get the known mole fractions from the fluidState
Dune::FieldVector<Scalar, numComponents-numMajorComponents> xKnown(0.0);
for (int knownCompIdx = 0; knownCompIdx < numComponents-numMajorComponents; ++knownCompIdx)
{
xKnown[knownCompIdx] = fluidState.moleFraction(knownPhaseIdx, knownCompIdx + numMajorComponents);
}
// compute all fugacity coefficients
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
paramCache.updatePhase(fluidState, phaseIdx);
// since we assume ideal mixtures, the fugacity
// coefficients of the components cannot depend on
// composition, i.e. the parameters in the cache are valid
for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
Scalar fugCoeff = FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
fluidState.setFugacityCoefficient(phaseIdx, compIdx, fugCoeff);
}
}
// create the linear system of equations which defines the
// mole fractions
Dune::FieldMatrix<Scalar, numComponents*numPhases, numComponents*numPhases> M(0.0);
Dune::FieldVector<Scalar, numComponents*numPhases> x(0.0);
Dune::FieldVector<Scalar, numComponents*numPhases> b(0.0);
// assemble the equations expressing the assumption that the
// sum of all mole fractions in each phase must be 1
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
int rowIdx = numComponents*(numPhases - 1) + phaseIdx;
b[rowIdx] = 1.0;
for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
int colIdx = phaseIdx*numComponents + compIdx;
M[rowIdx][colIdx] = 1.0;
}
}
//set the additional equations for the numComponents-numMajorComponents
//Components, of which the molefractions are known, set to molefraction(knownCompIdx)=xKnown
for(int knownCompIdx = 0; knownCompIdx < numComponents-numMajorComponents; ++knownCompIdx)
{
int rowIdx = numComponents + numPhases + knownCompIdx;
int colIdx = knownPhaseIdx*numComponents + knownCompIdx + numMajorComponents;
M[rowIdx][colIdx] = 1.0;
b[rowIdx] = xKnown[knownCompIdx];
}
// assemble the equations expressing the fact that the
// fugacities of each component is equal in all phases
for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
Scalar entryCol1 =
fluidState.fugacityCoefficient(/*phaseIdx=*/0, compIdx)
* fluidState.pressure(/*phaseIdx=*/0);
int col1Idx = compIdx;
for (int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx) {
int rowIdx = (phaseIdx - 1)*numComponents + compIdx;
int col2Idx = phaseIdx*numComponents + compIdx;
Scalar entryCol2 =
fluidState.fugacityCoefficient(phaseIdx, compIdx)
* fluidState.pressure(phaseIdx);
M[rowIdx][col1Idx] = entryCol1;
M[rowIdx][col2Idx] = -entryCol2;
}
}
//preconditioning of M to reduce condition number
//prevents matrix meeting dune's singularity criteria
for (int compIdx = 0; compIdx < numComponents; compIdx++)
{
if (compIdx < numMajorComponents)
{
for (int colIdx = 0; colIdx < numPhases*numComponents; colIdx++)
{
//Multiply row of main component (Raoult's Law) with 10e-5 (order of magn. of pressure)
M[compIdx][colIdx] *= 10e-5;
}
} else {
for (int colIdx = 0; colIdx < numPhases*numComponents; colIdx++)
{
//Multiply row of sec. components (Henry's Law) with 10e-9 (order of magn. of Henry constant)
M[compIdx][colIdx] *= 10e-9;
}
}
}
// solve for all mole fractions
M.solve(x, b);
// set all mole fractions and the the additional quantities in
// the fluid state
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
int rowIdx = phaseIdx*numComponents + compIdx;
fluidState.setMoleFraction(phaseIdx, compIdx, x[rowIdx]);
}
paramCache.updateComposition(fluidState, phaseIdx);
Scalar value = FluidSystem::density(fluidState, paramCache, phaseIdx);
fluidState.setDensity(phaseIdx, value);
if (setViscosity) {
value = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
fluidState.setViscosity(phaseIdx, value);
}
if (setEnthalpy) {
value = FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
fluidState.setEnthalpy(phaseIdx, value);
}
}
}
};
using Miscible2pNCComposition DUNE_DEPRECATED_MSG("Use MiscibleMultiPhaseComposition from dumux/material/constraintsolvers/misciblemultiphasecomposition.hh")
= MiscibleMultiPhaseComposition<Scalar, FluidSystem>;
} // end namespace Dumux
#endif
......@@ -62,11 +62,7 @@ class MiscibleMultiPhaseComposition
{
static constexpr int numPhases = FluidSystem::numPhases;
static constexpr int numComponents = FluidSystem::numComponents;
static_assert(numComponents == numPhases,
"This solver requires that the number fluid phases is equal "
"to the number of components");
static const int numMajorComponents = FluidSystem::numPhases;
public:
/*!
......@@ -89,7 +85,8 @@ public:
static void solve(FluidState &fluidState,
ParameterCache &paramCache,
bool setViscosity,
bool setEnthalpy)
bool setEnthalpy,
int knownPhaseIdx = 0)
{
#ifndef NDEBUG
// currently this solver can only handle fluid systems which
......@@ -98,9 +95,18 @@ public:
// newton method.)
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
assert(FluidSystem::isIdealMixture(phaseIdx));
}
#endif
//get the known mole fractions from the fluidState
//in a 2pnc system the n>2 mole fractions are primary variables and are already set in the fluidstate
Dune::FieldVector<Scalar, numComponents-numMajorComponents> xKnown(0.0);
for (int knownCompIdx = 0; knownCompIdx < numComponents-numMajorComponents; ++knownCompIdx)
{
xKnown[knownCompIdx] = fluidState.moleFraction(knownPhaseIdx, knownCompIdx + numMajorComponents);
}
// compute all fugacity coefficients
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
paramCache.updatePhase(fluidState, phaseIdx);
......@@ -134,6 +140,17 @@ public:
}
}
//set the additional equations for the numComponents-numMajorComponents
// this is only relevant if numphases != numcomponents e.g. in a 2pnc system
//Components, of which the molefractions are known, set to molefraction(knownCompIdx)=xKnown
for(int knownCompIdx = 0; knownCompIdx < numComponents-numMajorComponents; ++knownCompIdx)
{
int rowIdx = numComponents + numPhases + knownCompIdx;
int colIdx = knownPhaseIdx*numComponents + knownCompIdx + numMajorComponents;
M[rowIdx][colIdx] = 1.0;
b[rowIdx] = xKnown[knownCompIdx];
}
// assemble the equations expressing the fact that the
// fugacities of each component are equal in all phases
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
......@@ -176,6 +193,28 @@ public:
}
}
//preconditioning of M to reduce condition number
//prevents matrix meeting dune's singularity criteria
for (int compIdx = 0; compIdx < numComponents; compIdx++)
{
if (compIdx < numMajorComponents)
{
for (int colIdx = 0; colIdx < numPhases*numComponents; colIdx++)
{
//Multiply row of main component (Raoult's Law) with 10e-5 (order of magn. of pressure)
M[compIdx][colIdx] *= 10e-5;
}
} else {
for (int colIdx = 0; colIdx < numPhases*numComponents; colIdx++)
{
//Multiply row of sec. components (Henry's Law) with 10e-9 (order of magn. of Henry constant)
M[compIdx][colIdx] *= 10e-9;
}
}
}
// solve for all mole fractions
try
{
......
......@@ -435,6 +435,14 @@ public:
DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
}
template <class FluidState>
static Scalar kelvinVaporPressure(const FluidState &fluidState,
const int phaseIdx,
const int compIdx)
{
DUNE_THROW(Dune::NotImplemented, "FluidSystems::BrineAir::kelvinVaporPressure()");
}
using Base::diffusionCoefficient;
template <class FluidState>
......
......@@ -321,6 +321,31 @@ public:
return kelvinVaporPressure;
}
/*!
* \brief Vapor pressure including the Kelvin equation in \f$\mathrm{[Pa]}\f$
*
* Calculate the decreased vapor pressure due to capillarity
*
* \param fluidState An abitrary fluid state
* \param phaseIdx The index of the fluid phase to consider
* \param compIdx The index of the component to consider
*/
template <class FluidState>
static Scalar kelvinVaporPressure(const FluidState &fluidState,
const int phaseIdx,
const int compIdx)
{
assert(compIdx == wCompIdx && phaseIdx == wPhaseIdx);
using std::exp;
return fugacityCoefficient(fluidState, phaseIdx, compIdx)
* fluidState.pressure(phaseIdx)
* exp(-(fluidState.pressure(nPhaseIdx)-fluidState.pressure(wPhaseIdx))
/ density(fluidState, phaseIdx)
/ (Dumux::Constants<Scalar>::R / molarMass(compIdx))
/ fluidState.temperature());
}
/*!
* \brief Calculate the surface tension between water and air in \f$\mathrm{[\frac{N}{m}]}\f$,
* according to IAPWS Release on Surface Tension from September 1994.
......
......@@ -35,7 +35,7 @@
#include <dumux/material/fluidstates/compositional.hh>
#include <dumux/porousmediumflow/volumevariables.hh>
#include <dumux/material/constraintsolvers/computefromreferencephase.hh>
#include <dumux/material/constraintsolvers/miscible2pnccomposition.hh>
#include <dumux/material/constraintsolvers/misciblemultiphasecomposition.hh>
#include "indices.hh" // for formulation
......@@ -94,7 +94,7 @@ class TwoPNCVolumeVariables : public PorousMediumFlowVolumeVariables<TypeTag>
};
using Element = typename GridView::template Codim<0>::Entity;
using Miscible2pNCComposition = Dumux::Miscible2pNCComposition<Scalar, FluidSystem>;
using MiscibleMultiPhaseComposition = Dumux::MiscibleMultiPhaseComposition<Scalar, FluidSystem>;
using ComputeFromReferencePhase = Dumux::ComputeFromReferencePhase<Scalar, FluidSystem>;
static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
static_assert(useMoles, "use moles has to be set true in the 2pnc model");
......@@ -263,7 +263,7 @@ public:
// constraint solver
// set the known mole fractions in the fluidState so that they
// can be used by the Miscible2pNCComposition constraint solver
// can be used by the MiscibleMultiPhaseComposition constraint solver
unsigned int knownPhaseIdx = nPhaseIdx;
if (GET_PROP_VALUE(TypeTag, SetMoleFractionsForWettingPhase))
{
......@@ -275,11 +275,11 @@ public:
fluidState.setMoleFraction(knownPhaseIdx, compIdx, priVars[compIdx]);
}
Miscible2pNCComposition::solve(fluidState,
MiscibleMultiPhaseComposition::solve(fluidState,
paramCache,
knownPhaseIdx,
/*setViscosity=*/true,
/*setEnthalpy=*/false);
/*setEnthalpy=*/false,
knownPhaseIdx);
}
else if (phasePresence == nPhaseOnly)
{
......
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