diff --git a/dumux/material/constraintsolvers/misciblemultiphasecomposition.hh b/dumux/material/constraintsolvers/misciblemultiphasecomposition.hh
index 5821a9610ab0199c6b72bf27bced5703a40ea17e..28bd6ca3fadd9f1d75629607e94de752e0f9e045 100644
--- a/dumux/material/constraintsolvers/misciblemultiphasecomposition.hh
+++ b/dumux/material/constraintsolvers/misciblemultiphasecomposition.hh
@@ -57,7 +57,7 @@ namespace Dumux {
* - if the setViscosity parameter is true, also dynamic viscosities of *all* phases
* - if the setEnthalpy parameter is true, also specific enthalpies of *all* phases
*/
-template <class Scalar, class FluidSystem>
+template <class Scalar, class FluidSystem, bool useKelvinEquation = false>
class MiscibleMultiPhaseComposition
{
static constexpr int numPhases = FluidSystem::numPhases;
@@ -72,6 +72,14 @@ public:
/*!
* \brief @copybrief Dumux::MiscibleMultiPhaseComposition
*
+ * This function additionally considers a lowering of the saturation vapor pressure
+ * of the wetting phase by the Kelvin equation:
+ * \f[
+ * p^\textrm{w}_\textrm{sat,Kelvin}
+ * = p^\textrm{w}_\textrm{sat}
+ * \exp \left( -\frac{p_\textrm{c}}{\varrho_\textrm{w} R_\textrm{w} T} \right)
+ * \f]
+ *
* \param fluidState A container with the current (physical) state of the fluid
* \param paramCache A container for iterative calculation of fluid composition
* \param setViscosity Should the viscosity be set in the fluidstate?
@@ -128,22 +136,43 @@ public:
// assemble the equations expressing the fact that the
// fugacities of each component are equal in all phases
- for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
- Scalar entryCol1 =
- fluidState.fugacityCoefficient(/*phaseIdx=*/0, compIdx)
- * fluidState.pressure(/*phaseIdx=*/0);
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ {
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;
+ Scalar entryPhase0 = 0.0;
+ for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ {
+ Scalar entry = fluidState.fugacityCoefficient(phaseIdx, compIdx)
+ * fluidState.pressure(phaseIdx);
+
+ // modify the saturation vapor pressure of the wetting component by the Kelvin equation
+ if (compIdx == FluidSystem::wCompIdx
+ && phaseIdx == FluidSystem::wPhaseIdx
+ && useKelvinEquation)
+ {
+ // a new fluidState is needed, because mole fractions are unknown
+ FluidState purePhaseFluidState;
+ // assign all phase pressures, needed for capillary pressure
+ for (unsigned int idx = 0; idx < numPhases; ++idx)
+ {
+ purePhaseFluidState.setPressure(idx, fluidState.pressure(idx));
+ }
+ purePhaseFluidState.setTemperature(fluidState.temperature());
+ purePhaseFluidState.setMoleFraction(phaseIdx, compIdx, 1.0);
+ entry = FluidSystem::kelvinVaporPressure(purePhaseFluidState, phaseIdx, compIdx);
+ }
+
+ if (phaseIdx == 0)
+ {
+ entryPhase0 = entry;
+ }
+ else
+ {
+ int rowIdx = (phaseIdx - 1)*numComponents + compIdx;
+ int col2Idx = phaseIdx*numComponents + compIdx;
+ M[rowIdx][col1Idx] = entryPhase0;
+ M[rowIdx][col2Idx] = -entry;
+ }
}
}
diff --git a/dumux/material/fluidsystems/h2oair.hh b/dumux/material/fluidsystems/h2oair.hh
index daeb4d0f4e7baaef11dcb4a1a1adff511fd61fe2..9d541e6cf06abaa0af4e64ec94c7ee9823d5a481 100644
--- a/dumux/material/fluidsystems/h2oair.hh
+++ b/dumux/material/fluidsystems/h2oair.hh
@@ -298,6 +298,30 @@ public:
DUNE_THROW(Dune::NotImplemented, "Invalid component index " << compIdx);
}
+ /*!
+ * \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);
+
+ return fugacityCoefficient(fluidState, phaseIdx, compIdx)
+ * fluidState.pressure(phaseIdx)
+ * std::exp(-(fluidState.pressure(nPhaseIdx)-fluidState.pressure(wPhaseIdx))
+ / density(fluidState, phaseIdx)
+ / (Dumux::Constants<Scalar>::R / molarMass(compIdx))
+ / fluidState.temperature());
+ }
+
/*!
* \brief Molar volume of a component at the critical point \f$\mathrm{[m^3/mol]}\f$.
*
diff --git a/dumux/material/fluidsystems/h2oairmesitylene.hh b/dumux/material/fluidsystems/h2oairmesitylene.hh
index 83082a156d97a7dc37da5c261e930deac6607c05..ca8763f8922ee57884286a47d10258fd655a1047 100644
--- a/dumux/material/fluidsystems/h2oairmesitylene.hh
+++ b/dumux/material/fluidsystems/h2oairmesitylene.hh
@@ -474,6 +474,14 @@ public:
return 1.0;
}
+ template <class FluidState>
+ static Scalar kelvinVaporPressure(const FluidState &fluidState,
+ const int phaseIdx,
+ const int compIdx)
+ {
+ DUNE_THROW(Dune::NotImplemented, "FluidSystems::H2OAirMesitylene::kelvinVaporPressure()");
+ }
+
using Base::enthalpy;
/*!
* \brief Given all mole fractions in a phase, return the specific
diff --git a/dumux/material/fluidsystems/h2oairxylene.hh b/dumux/material/fluidsystems/h2oairxylene.hh
index bfa947c48d78e9804d8d901e1f9606076bd295cf..41e652b32a54016b5548c31b48532ed3737d06dc 100644
--- a/dumux/material/fluidsystems/h2oairxylene.hh
+++ b/dumux/material/fluidsystems/h2oairxylene.hh
@@ -473,6 +473,14 @@ public:
return 1.0;
}
+ template <class FluidState>
+ static Scalar kelvinVaporPressure(const FluidState &fluidState,
+ const int phaseIdx,
+ const int compIdx)
+ {
+ DUNE_THROW(Dune::NotImplemented, "FluidSystems::H2OAirXylene::kelvinVaporPressure()");
+ }
+
using Base::enthalpy;
/*!
* \brief Given all mole fractions in a phase, return the specific
diff --git a/dumux/material/fluidsystems/h2on2.hh b/dumux/material/fluidsystems/h2on2.hh
index 50976ea8f953d134eb6acedb878ab81773547c98..27f4e2de1a7f10da6c66e55996180d8c8c482178 100644
--- a/dumux/material/fluidsystems/h2on2.hh
+++ b/dumux/material/fluidsystems/h2on2.hh
@@ -258,6 +258,30 @@ public:
return pcrit[compIdx];
}
+ /*!
+ * \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);
+
+ return fugacityCoefficient(fluidState, phaseIdx, compIdx)
+ * fluidState.pressure(phaseIdx)
+ * std::exp(-(fluidState.pressure(nPhaseIdx)-fluidState.pressure(wPhaseIdx))
+ / density(fluidState, phaseIdx)
+ / (Dumux::Constants<Scalar>::R / molarMass(compIdx))
+ / fluidState.temperature());
+ }
+
/*!
* \brief Molar volume of a component at the critical point \f$\mathrm{[m^3/mol]}\f$.
*
diff --git a/dumux/porousmediumflow/2p2c/implicit/properties.hh b/dumux/porousmediumflow/2p2c/implicit/properties.hh
index aa494f5df69d08fc616179888ae60665c5589f39..ddca6ae7db9409b1edac6448c313c28a317011f7 100644
--- a/dumux/porousmediumflow/2p2c/implicit/properties.hh
+++ b/dumux/porousmediumflow/2p2c/implicit/properties.hh
@@ -71,6 +71,7 @@ NEW_PROP_TAG(EffectiveDiffusivityModel); //!< The employed model for the computa
NEW_PROP_TAG(ProblemEnableGravity); //!< Returns whether gravity is considered in the problem
NEW_PROP_TAG(UseMoles); //!< Defines whether mole (true) or mass (false) fractions are used
NEW_PROP_TAG(UseConstraintSolver); //!< Determines whether the constraint solver should be used
+NEW_PROP_TAG(UseKelvinEquation); //!< Determines whether the Kelvin equation is used
NEW_PROP_TAG(ImplicitMassUpwindWeight); //!< The value of the upwind weight for the mass conservation equations
NEW_PROP_TAG(ImplicitMobilityUpwindWeight); //!< Weight for the upwind mobility in the velocity calculation
diff --git a/dumux/porousmediumflow/2p2c/implicit/propertydefaults.hh b/dumux/porousmediumflow/2p2c/implicit/propertydefaults.hh
index bedf66eca034fc3f4f0d1288ac84634e56dab378..c416c273f18836579b85f3d7834bb78fb4db8e91 100644
--- a/dumux/porousmediumflow/2p2c/implicit/propertydefaults.hh
+++ b/dumux/porousmediumflow/2p2c/implicit/propertydefaults.hh
@@ -172,6 +172,9 @@ SET_BOOL_PROP(TwoPTwoC, UseMoles, true);
//! Determines whether the constraint solver is used
SET_BOOL_PROP(TwoPTwoC, UseConstraintSolver, true);
+//! Determines whether the Kelvin equation is used to adapt the saturation vapor pressure
+SET_BOOL_PROP(TwoPTwoC, UseKelvinEquation, false);
+
//! Set default value for the Forchheimer coefficient
// Source: Ward, J.C. 1964 Turbulent flow in porous media. ASCE J. Hydraul. Div 90.
// Actually the Forchheimer coefficient is also a function of the dimensions of the
diff --git a/dumux/porousmediumflow/2p2c/implicit/volumevariables.hh b/dumux/porousmediumflow/2p2c/implicit/volumevariables.hh
index 1a0756ea3b344b43a0d48579bfd34306be72010a..6df37beb34f816ede77b74174c62f0bd626831d8 100644
--- a/dumux/porousmediumflow/2p2c/implicit/volumevariables.hh
+++ b/dumux/porousmediumflow/2p2c/implicit/volumevariables.hh
@@ -83,6 +83,9 @@ class TwoPTwoCVolumeVariables : public ImplicitVolumeVariables<TypeTag>
switchIdx = Indices::switchIdx,
pressureIdx = Indices::pressureIdx
};
+ static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
+ static const bool useConstraintSolver = GET_PROP_VALUE(TypeTag, UseConstraintSolver);
+ static const bool useKelvinEquation = GET_PROP_VALUE(TypeTag, UseKelvinEquation);
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GridView::template Codim<0>::Entity Element;
@@ -90,9 +93,7 @@ class TwoPTwoCVolumeVariables : public ImplicitVolumeVariables<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
- typedef Dumux::MiscibleMultiPhaseComposition<Scalar, FluidSystem> MiscibleMultiPhaseComposition;
- static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
- static const bool useConstraintSolver = GET_PROP_VALUE(TypeTag, UseConstraintSolver);
+ typedef Dumux::MiscibleMultiPhaseComposition<Scalar, FluidSystem, useKelvinEquation> MiscibleMultiPhaseComposition;
static_assert(useMoles || (!useMoles && useConstraintSolver),
"if UseMoles is set false, UseConstraintSolver has to be set to true");
typedef Dumux::ComputeFromReferencePhase<Scalar, FluidSystem> ComputeFromReferencePhase;
@@ -275,6 +276,9 @@ public:
wPhaseIdx,
wCompIdx) * pw;
+ if (useKelvinEquation)
+ partPressH2O = FluidSystem::kelvinVaporPressure(fluidState, wPhaseIdx, wCompIdx);
+
// get the partial pressure of the main component of the the nonwetting (gas) phase ("Air")
Scalar partPressAir = pn - partPressH2O;