From 916f27c94860ef1efa4e1b4f5fe67bf6ec62a4c9 Mon Sep 17 00:00:00 2001
From: Andreas Lauser <and@poware.org>
Date: Fri, 16 Dec 2011 13:23:35 +0000
Subject: [PATCH] add a flash constraint solver based on non-linear
complementarity functions
also add a test
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@7068 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
configure.ac | 1 +
dumux/boxmodels/MpNc/MpNcvolumevariables.hh | 43 +-
.../compositionfromfugacities.hh | 10 +-
.../computefromreferencephase.hh | 6 +
.../misciblemultiphasecomposition.hh | 8 +-
.../MpNcconstraintsolvers/ncpflash.hh | 662 ++++++++++++++++++
.../MpNcfluidsystems/h2on2fluidsystem.hh | 15 +-
.../MpNcfluidsystems/nullparametercache.hh | 12 -
.../MpNcfluidsystems/parametercachebase.hh | 36 +-
test/material/Makefile.am | 2 +-
test/material/ncpflash/Makefile.am | 5 +
test/material/ncpflash/test_ncpflash.cc | 294 ++++++++
12 files changed, 1058 insertions(+), 36 deletions(-)
create mode 100644 dumux/material/MpNcconstraintsolvers/ncpflash.hh
create mode 100644 test/material/ncpflash/Makefile.am
create mode 100644 test/material/ncpflash/test_ncpflash.cc
diff --git a/configure.ac b/configure.ac
index 3f08f975a4..e7decf83d7 100644
--- a/configure.ac
+++ b/configure.ac
@@ -84,6 +84,7 @@ AC_CONFIG_FILES([dumux.pc
test/decoupled/2p2c/Makefile
test/material/Makefile
test/material/tabulation/Makefile
+ test/material/ncpflash/Makefile
tutorial/Makefile
])
diff --git a/dumux/boxmodels/MpNc/MpNcvolumevariables.hh b/dumux/boxmodels/MpNc/MpNcvolumevariables.hh
index 163225029b..4afe21e577 100644
--- a/dumux/boxmodels/MpNc/MpNcvolumevariables.hh
+++ b/dumux/boxmodels/MpNc/MpNcvolumevariables.hh
@@ -35,6 +35,7 @@
#include "MpNcvolumevariablesia.hh"
#include <dumux/boxmodels/common/boxvolumevariables.hh>
+#include <dumux/material/MpNcconstraintsolvers/ncpflash.hh>
namespace Dumux
{
@@ -240,8 +241,48 @@ public:
elemGeom,
scvIdx);
IAVolumeVariables::checkDefined();
-
checkDefined();
+
+
+#warning "HACK for testing the NCP flash calculation!"
+#if 0
+ FluidState foo;
+ foo.assign(fluidState_);
+
+ typedef Dumux::NcpFlash<Scalar, FluidSystem> NcpFlash;
+ typedef typename NcpFlash::ComponentVector ComponentVector;
+ ComponentVector globalMolarities(0.0);
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ globalMolarities[compIdx] +=
+ fluidState_.saturation(phaseIdx)
+ * fluidState_.molarity(phaseIdx, compIdx);
+
+ NcpFlash::guessInitial(foo, paramCache, globalMolarities);
+ NcpFlash::template solve<MaterialLaw>(foo, paramCache, materialParams, globalMolarities);
+
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ if (foo.pressure(phaseIdx) == fluidState_.pressure(phaseIdx))
+ continue;
+ Scalar error = 1 - fluidState_.pressure(phaseIdx)/foo.pressure(phaseIdx);
+ if (std::abs(error) > 1e-6) {
+ std::cout << "pressure error phase " << phaseIdx << ": " << foo.pressure(phaseIdx) << " vs " << fluidState_.pressure(phaseIdx) << " error=" << 1 - fluidState_.pressure(phaseIdx)/foo.pressure(phaseIdx) << "\n";
+ std::cout << "x_0: " << fluidState_.moleFraction(0, 0) << " vs " << foo.moleFraction(0, 0) << "\n";
+ std::cout << "x_1: " << fluidState_.moleFraction(0, 1) << " vs " << foo.moleFraction(0, 1) << "\n";
+ }
+ };
+
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ if (foo.saturation(phaseIdx) == fluidState_.saturation(phaseIdx))
+ continue;
+
+ Scalar error = 1 - (1 + fluidState_.saturation(phaseIdx))/(1 + foo.saturation(phaseIdx));
+ if (std::abs(error) > 1e-6)
+ std::cout << "saturation error phase " << phaseIdx << ": " << foo.saturation(phaseIdx) << " vs " << fluidState_.saturation(phaseIdx) << " error=" << 1 - fluidState_.saturation(phaseIdx)/foo.saturation(phaseIdx) << "\n";
+ };
+
+ //exit(1);
+#endif
}
/*!
diff --git a/dumux/material/MpNcconstraintsolvers/compositionfromfugacities.hh b/dumux/material/MpNcconstraintsolvers/compositionfromfugacities.hh
index 79cf8f60ed..e236eacba5 100644
--- a/dumux/material/MpNcconstraintsolvers/compositionfromfugacities.hh
+++ b/dumux/material/MpNcconstraintsolvers/compositionfromfugacities.hh
@@ -26,7 +26,11 @@
#ifndef DUMUX_COMPOSITION_FROM_FUGACITIES_HH
#define DUMUX_COMPOSITION_FROM_FUGACITIES_HH
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+
#include <dumux/common/exceptions.hh>
+#include <dumux/common/valgrind.hh>
namespace Dumux {
@@ -237,7 +241,7 @@ protected:
// deviate the mole fraction of the i-th component
Scalar x_i = fluidState.moleFraction(phaseIdx, i);
fluidState.setMoleFraction(phaseIdx, i, x_i + eps);
- paramCache.updatePhaseSingleMoleFraction(fluidState, phaseIdx, i);
+ paramCache.updateSingleMoleFraction(fluidState, phaseIdx, i);
// compute new defect and derivative for all component
// fugacities
@@ -262,7 +266,7 @@ protected:
// reset composition to original value
fluidState.setMoleFraction(phaseIdx, i, x_i);
- paramCache.updatePhaseSingleMoleFraction(fluidState, phaseIdx, i);
+ paramCache.updateSingleMoleFraction(fluidState, phaseIdx, i);
// end forward differences
////////
@@ -320,7 +324,7 @@ protected:
//sumMoleFrac += std::abs(newx);
}
- paramCache.updatePhaseComposition(fluidState, phaseIdx);
+ paramCache.updateComposition(fluidState, phaseIdx);
/*
// if the sum of the mole fractions gets 0, we take the
diff --git a/dumux/material/MpNcconstraintsolvers/computefromreferencephase.hh b/dumux/material/MpNcconstraintsolvers/computefromreferencephase.hh
index f251dc06c0..1848ea40f2 100644
--- a/dumux/material/MpNcconstraintsolvers/computefromreferencephase.hh
+++ b/dumux/material/MpNcconstraintsolvers/computefromreferencephase.hh
@@ -32,6 +32,12 @@
#include <dumux/material/MpNcconstraintsolvers/compositionfromfugacities.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+
+#include <dumux/common/exceptions.hh>
+#include <dumux/common/valgrind.hh>
+
namespace Dumux {
/*!
diff --git a/dumux/material/MpNcconstraintsolvers/misciblemultiphasecomposition.hh b/dumux/material/MpNcconstraintsolvers/misciblemultiphasecomposition.hh
index d7db5a5512..a75ec2fb92 100644
--- a/dumux/material/MpNcconstraintsolvers/misciblemultiphasecomposition.hh
+++ b/dumux/material/MpNcconstraintsolvers/misciblemultiphasecomposition.hh
@@ -27,6 +27,12 @@
#ifndef DUMUX_MISCIBLE_MULTIPHASE_COMPOSITION_HH
#define DUMUX_MISCIBLE_MULTIPHASE_COMPOSITION_HH
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+
+#include <dumux/common/exceptions.hh>
+#include <dumux/common/valgrind.hh>
+
namespace Dumux {
/*!
* \brief Computes the composition of all phases of a N-phase,
@@ -167,7 +173,7 @@ public:
int rowIdx = phaseIdx*numComponents + compIdx;
fluidState.setMoleFraction(phaseIdx, compIdx, x[rowIdx]);
}
- paramCache.updatePhaseComposition(fluidState, phaseIdx);
+ paramCache.updateComposition(fluidState, phaseIdx);
Scalar value = FluidSystem::density(fluidState, paramCache, phaseIdx);
fluidState.setDensity(phaseIdx, value);
diff --git a/dumux/material/MpNcconstraintsolvers/ncpflash.hh b/dumux/material/MpNcconstraintsolvers/ncpflash.hh
new file mode 100644
index 0000000000..89827025c1
--- /dev/null
+++ b/dumux/material/MpNcconstraintsolvers/ncpflash.hh
@@ -0,0 +1,662 @@
+/*****************************************************************************
+ * Copyright (C) 2011 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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 Determines the phase compositions, pressures and saturations
+ * given the total mass of all components.
+ */
+#ifndef DUMUX_NCP_FLASH_HH
+#define DUMUX_NCP_FLASH_HH
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+
+#include <dumux/common/exceptions.hh>
+#include <dumux/common/valgrind.hh>
+
+namespace Dumux {
+
+/*!
+ * \brief Determines the phase compositions, pressures and saturations
+ * given the total mass of all components.
+ */
+template <class Scalar, class FluidSystem>
+class NcpFlash
+{
+ enum { numPhases = FluidSystem::numPhases };
+ enum { numComponents = FluidSystem::numComponents };
+
+ typedef typename FluidSystem::ParameterCache ParameterCache;
+
+ // In a M-phase, N-component context, we have the following
+ // unknowns:
+ //
+ // - M pressures
+ // - M saturations
+ // - M*N mole fractions
+ //
+ // This sums up to M*(N + 2). On the equations side of things,
+ // we have:
+ //
+ // - (M - 1)*N equation stemming from the fact that the
+ // fugacity of any component is the same in all phases
+ // - 1 equation from the closure condition of all saturations
+ // (they sum up to 1)
+ // - M - 1 constraints from the capillary pressures
+ // (-> p_\beta = p_\alpha + p_c\alpha,\beta)
+ // - N constraints from the fact that the total mass of each
+ // component is given (-> sum_\alpha rhoMolar_\alpha *
+ // x_\alpha^\kappa = const)
+ // - M model constraints. Here we use the NCP constraints
+ // (-> 0 = min{S_\alpha, 1 - \sum_\kappa x_\alpha^\kappa})
+ //
+ // this also sums up to M*(N + 2).
+ //
+ // We use the following catches: Capillary pressures are taken
+ // into accout expicitly, so that only the pressure of the first
+ // phase is solved implicitly, also the closure condition for the
+ // saturations is taken into account explicitly, which means, that
+ // we don't need to implicitly solve for the last
+ // saturation. These two measures reduce the number of unknowns to
+ // M*(N + 1), namely:
+ //
+ // - 1 pressure
+ // - M - 1 saturations
+ // - M*N mole fractions
+ static constexpr int numEq = numPhases*(numComponents + 1);
+
+ typedef Dune::FieldMatrix<Scalar, numEq, numEq> Matrix;
+ typedef Dune::FieldVector<Scalar, numEq> Vector;
+
+public:
+ typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
+
+ /*!
+ * \brief Guess initial values for all quantities.
+ */
+ template <class FluidState>
+ static void guessInitial(FluidState &fluidState,
+ ParameterCache ¶mCache,
+ const ComponentVector &globalMolarities)
+ {
+ // the sum of all molarities
+ Scalar sumMoles = 0;
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ sumMoles += globalMolarities[compIdx];
+
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+ // composition
+ for (int compIdx = 0; compIdx < numComponents; ++ compIdx)
+ fluidState.setMoleFraction(phaseIdx,
+ compIdx,
+ globalMolarities[compIdx]/sumMoles);
+
+ // pressure. use atmospheric pressure as initial guess
+ fluidState.setPressure(phaseIdx, 2e5);
+
+ // saturation. assume all fluids to be equally distributed
+ fluidState.setSaturation(phaseIdx, 1.0/numPhases);
+ }
+
+ // set the fugacity coefficients of all components in all phases
+ paramCache.updateAll(fluidState);
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+ for (int compIdx = 0; compIdx < numComponents; ++ compIdx) {
+ Scalar phi = FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
+ fluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
+ }
+ }
+ };
+
+ /*!
+ * \brief Calculates the chemical equilibrium from the component
+ * fugacities in a phase.
+ *
+ * The phase's fugacities must already be set.
+ */
+ template <class MaterialLaw, class FluidState>
+ static void solve(FluidState &fluidState,
+ ParameterCache ¶mCache,
+ const typename MaterialLaw::Params &matParams,
+ const ComponentVector &globalMolarities)
+ {
+ Dune::FMatrixPrecision<Scalar>::set_singular_limit(1e-25);
+
+ /////////////////////////
+ // Newton method
+ /////////////////////////
+
+ // Jacobian matrix
+ Matrix J;
+ // solution, i.e. phase composition
+ Vector deltaX;
+ // right hand side
+ Vector b;
+
+ Valgrind::SetUndefined(J);
+ Valgrind::SetUndefined(deltaX);
+ Valgrind::SetUndefined(b);
+
+ // make the fluid state consistent with the fluid system.
+ completeFluidState_<MaterialLaw>(fluidState,
+ paramCache,
+ matParams);
+
+ /*
+ std::cout << "--------------------\n";
+ std::cout << "globalMolarities: ";
+ for (int compIdx = 0; compIdx < numComponents; ++ compIdx)
+ std::cout << globalMolarities[compIdx] << " ";
+ std::cout << "\n";
+ */
+
+ const int nMax = 50; // <- maximum number of newton iterations
+ for (int nIdx = 0; nIdx < nMax; ++nIdx) {
+ // calculate Jacobian matrix and right hand side
+ linearize_<MaterialLaw>(J, b, fluidState, paramCache, matParams, globalMolarities);
+ Valgrind::CheckDefined(J);
+ Valgrind::CheckDefined(b);
+
+ // Solve J*x = b
+ deltaX = 0;
+
+ try { J.solve(deltaX, b); }
+ catch (Dune::FMatrixError e)
+ {
+ /*
+ std::cout << "error: " << e << "\n";
+ std::cout << "b: " << b << "\n";
+ */
+
+ throw Dumux::NumericalProblem(e.what());
+ }
+ Valgrind::CheckDefined(deltaX);
+
+ /*
+ printFluidState_(fluidState);
+ std::cout << "J:\n";
+ for (int i = 0; i < numEq; ++i) {
+ for (int j = 0; j < numEq; ++j) {
+ std::ostringstream os;
+ os << J[i][j];
+
+ std::string s(os.str());
+ do {
+ s += " ";
+ } while (s.size() < 20);
+ std::cout << s;
+ }
+ std::cout << "\n";
+ };
+
+ std::cout << "deltaX: " << deltaX << "\n";
+ std::cout << "---------------\n";
+ */
+
+ // update the fluid quantities.
+ Scalar relError = update_<MaterialLaw>(fluidState, paramCache, matParams, deltaX);
+
+ if (relError < 1e-9)
+ return;
+ }
+
+ /*
+ printFluidState_(fluidState);
+ std::cout << "globalMolarities: ";
+ for (int compIdx = 0; compIdx < numComponents; ++ compIdx)
+ std::cout << globalMolarities[compIdx] << " ";
+ std::cout << "\n";
+ */
+
+ DUNE_THROW(NumericalProblem,
+ "Flash calculation failed."
+ " {c_alpha^kappa} = {" << globalMolarities << "}, T = "
+ << fluidState.temperature(/*phaseIdx=*/0));
+ };
+
+
+protected:
+ template <class FluidState>
+ static void printFluidState_(const FluidState &fs)
+ {
+ std::cout << "saturations: ";
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ std::cout << fs.saturation(phaseIdx) << " ";
+ std::cout << "\n";
+
+ std::cout << "pressures: ";
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ std::cout << fs.pressure(phaseIdx) << " ";
+ std::cout << "\n";
+
+ std::cout << "densities: ";
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ std::cout << fs.density(phaseIdx) << " ";
+ std::cout << "\n";
+
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ std::cout << "composition " << FluidSystem::phaseName(phaseIdx) << "Phase: ";
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
+ std::cout << fs.moleFraction(phaseIdx, compIdx) << " ";
+ }
+ std::cout << "\n";
+ }
+
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ std::cout << "fugacities " << FluidSystem::phaseName(phaseIdx) << "Phase: ";
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
+ std::cout << fs.fugacity(phaseIdx, compIdx) << " ";
+ }
+ std::cout << "\n";
+ }
+
+ std::cout << "global component molarities: ";
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
+ Scalar sum = 0;
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ sum += fs.saturation(phaseIdx)*fs.molarity(phaseIdx, compIdx);
+ }
+ std::cout << sum << " ";
+ }
+ std::cout << "\n";
+ }
+
+ template <class MaterialLaw, class FluidState>
+ static void linearize_(Matrix &J,
+ Vector &b,
+ FluidState &fluidState,
+ ParameterCache ¶mCache,
+ const typename MaterialLaw::Params &matParams,
+ const ComponentVector &globalMolarities)
+ {
+ FluidState origFluidState(fluidState);
+ ParameterCache origParamCache(paramCache);
+
+ Vector tmp;
+
+ // reset jacobian
+ J = 0;
+
+ Valgrind::SetUndefined(b);
+ calculateDefect_(b, fluidState, fluidState, globalMolarities);
+ Valgrind::CheckDefined(b);
+
+ // assemble jacobian matrix
+ for (int pvIdx = 0; pvIdx < numEq; ++ pvIdx) {
+ ////////
+ // approximately calculate partial derivatives of the
+ // fugacity defect of all components in regard to the mole
+ // fraction of the i-th component. This is done via
+ // forward differences
+
+ // deviate the mole fraction of the i-th component
+ Scalar x_i = getQuantity_(fluidState, pvIdx);
+ const Scalar eps = 1e-8/quantityWeight_(fluidState, pvIdx);
+ setQuantity_<MaterialLaw>(fluidState, paramCache, matParams, pvIdx, x_i + eps);
+ assert(getQuantity_(fluidState, pvIdx) == x_i + eps);
+
+ // compute derivative of the defect
+ calculateDefect_(tmp, origFluidState, fluidState, globalMolarities);
+ tmp -= b;
+ tmp /= eps;
+
+ // store derivative in jacobian matrix
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ J[eqIdx][pvIdx] = tmp[eqIdx];
+
+ // fluid state and parameter cache to their original values
+ fluidState = origFluidState;
+ paramCache = origParamCache;
+
+ // end forward differences
+ ////////
+ }
+ }
+
+ template <class FluidState>
+ static void calculateDefect_(Vector &b,
+ const FluidState &fluidStateEval,
+ const FluidState &fluidState,
+ const ComponentVector &globalMolarities)
+ {
+ int eqIdx = 0;
+
+ // fugacity of any component must be equal in all phases
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
+ for (int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx) {
+ b[eqIdx] =
+ fluidState.fugacity(/*phaseIdx=*/0, compIdx) -
+ fluidState.fugacity(phaseIdx, compIdx);
+ ++eqIdx;
+ }
+ }
+
+ assert(eqIdx == numComponents*(numPhases - 1));
+
+ // the fact saturations must sum up to 1 is included explicitly!
+
+ // capillary pressures are explicitly included!
+
+ // global molarities are given
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
+ b[eqIdx] = 0.0;
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ b[eqIdx] +=
+ fluidState.saturation(phaseIdx)
+ * fluidState.molarity(phaseIdx, compIdx);
+ }
+
+ b[eqIdx] -= globalMolarities[compIdx];
+ ++eqIdx;
+ }
+
+ // model assumptions (-> non-linear complementarity functions)
+ // must be adhered
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ Scalar sumMoleFracEval = 0.0;
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ sumMoleFracEval += fluidStateEval.moleFraction(phaseIdx, compIdx);
+
+ if (1.0 - sumMoleFracEval > fluidStateEval.saturation(phaseIdx)) {
+ b[eqIdx] = fluidState.saturation(phaseIdx);
+ }
+ else {
+ Scalar sumMoleFrac = 0.0;
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ sumMoleFrac += fluidState.moleFraction(phaseIdx, compIdx);
+ b[eqIdx] = 1.0 - sumMoleFrac;
+ }
+
+ ++eqIdx;
+ }
+ }
+
+ template <class MaterialLaw, class FluidState>
+ static Scalar update_(FluidState &fluidState,
+ ParameterCache ¶mCache,
+ const typename MaterialLaw::Params &matParams,
+ const Vector &deltaX)
+ {
+ Scalar relError = 0;
+ for (int pvIdx = 0; pvIdx < numEq; ++ pvIdx) {
+ Scalar tmp = getQuantity_(fluidState, pvIdx);
+ Scalar delta = deltaX[pvIdx];
+
+ relError = std::max(relError, std::abs(delta)*quantityWeight_(fluidState, pvIdx));
+
+ if (isSaturationIdx_(pvIdx)) {
+ // dampen to at most 20% change in saturation per
+ // iteration
+ delta = std::min(0.2, std::max(-0.2, delta));
+ }
+ else if (isMoleFracIdx_(pvIdx)) {
+ // dampen to at most 15% change in mole fraction per
+ // iteration
+ delta = std::min(0.15, std::max(-0.15, delta));
+ }
+ else if (isPressureIdx_(pvIdx)) {
+ // dampen to at most 15% change in pressure per
+ // iteration
+ delta = std::min(0.15*fluidState.pressure(0), std::max(-0.15*fluidState.pressure(0), delta));
+ };
+
+ setQuantityRaw_(fluidState, pvIdx, tmp - delta);
+ }
+
+ // make sure all saturations, pressures and mole fractions are non-negative
+ Scalar sumSat = 0;
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ Scalar value = fluidState.saturation(phaseIdx);
+ if (value < -0.05) {
+ value = -0.05;
+ fluidState.setSaturation(phaseIdx, value);
+ }
+ sumSat += value;
+
+ value = fluidState.pressure(phaseIdx);
+ if (value < 0)
+ fluidState.setPressure(phaseIdx, 0.0);
+
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
+ value = fluidState.moleFraction(phaseIdx, compIdx);
+ if (value < 0)
+ fluidState.setMoleFraction(phaseIdx, compIdx, 0.0);
+ }
+ };
+
+ // last saturation
+ if (sumSat > 1.05) {
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ Scalar value = fluidState.saturation(phaseIdx)/(0.95*sumSat);
+ fluidState.setSaturation(phaseIdx, value);
+ };
+ };
+
+ completeFluidState_<MaterialLaw>(fluidState, paramCache, matParams);
+
+ return relError;
+ }
+
+ template <class MaterialLaw, class FluidState>
+ static void completeFluidState_(FluidState &fluidState,
+ ParameterCache ¶mCache,
+ const typename MaterialLaw::Params &matParams)
+ {
+ // calculate the saturation of the last phase as a function of
+ // the other saturations
+ Scalar sumSat = 0.0;
+ for (int phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
+ sumSat += fluidState.saturation(phaseIdx);
+ fluidState.setSaturation(/*phaseIdx=*/numPhases - 1, 1.0 - sumSat);
+
+ // update the pressures using the material law (saturations
+ // and first pressure are already set because it is implicitly
+ // solved for.)
+ ComponentVector pC;
+ MaterialLaw::capillaryPressures(pC, matParams, fluidState);
+ for (int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
+ fluidState.setPressure(phaseIdx,
+ fluidState.pressure(0)
+ + (pC[phaseIdx] - pC[0]));
+
+ // update the parameter cache
+ paramCache.updateAll(fluidState, /*except=*/ParameterCache::Temperature);
+
+ // update all densities and fugacity coefficients
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+ Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
+ fluidState.setDensity(phaseIdx, rho);
+
+ for (int compIdx = 0; compIdx < numComponents; ++ compIdx) {
+ Scalar phi = FluidSystem::fugacityCoefficient( fluidState, paramCache, phaseIdx, compIdx);
+ fluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
+ }
+ }
+ }
+
+ static bool isPressureIdx_(int pvIdx)
+ { return pvIdx == 0; }
+
+ static bool isSaturationIdx_(int pvIdx)
+ { return 1 <= pvIdx && pvIdx < numPhases; }
+
+ static bool isMoleFracIdx_(int pvIdx)
+ { return numPhases <= pvIdx && pvIdx < numPhases + numPhases*numComponents; }
+
+ // retrieves a quantity from the fluid state
+ template <class FluidState>
+ static Scalar getQuantity_(const FluidState &fs, int pvIdx)
+ {
+ assert(pvIdx < numEq);
+
+ // first pressure
+ if (pvIdx < 1) {
+ int phaseIdx = 0;
+ return fs.pressure(phaseIdx);
+ }
+ // first M - 1 saturations
+ else if (pvIdx < numPhases) {
+ int phaseIdx = pvIdx - 1;
+ return fs.saturation(phaseIdx);
+ }
+ // mole fractions
+ else // if (pvIdx < numPhases + numPhases*numComponents)
+ {
+ int phaseIdx = (pvIdx - numPhases)/numComponents;
+ int compIdx = (pvIdx - numPhases)%numComponents;
+ return fs.moleFraction(phaseIdx, compIdx);
+ }
+ };
+
+ // set a quantity in the fluid state
+ template <class MaterialLaw, class FluidState>
+ static void setQuantity_(FluidState &fs,
+ ParameterCache ¶mCache,
+ const typename MaterialLaw::Params &matParams,
+ int pvIdx,
+ Scalar value)
+ {
+ assert(0 <= pvIdx && pvIdx < numEq);
+
+ if (pvIdx < 1) { // <- first pressure
+ Scalar delta = value - fs.pressure(0);
+
+ // set all pressures. here we assume that the capillary
+ // pressure does not depend on absolute pressure.
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ fs.setPressure(phaseIdx, fs.pressure(phaseIdx) + delta);
+ paramCache.updateAllPressures(fs);
+
+ // update all densities and fugacity coefficients
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
+ fs.setDensity(phaseIdx, rho);
+
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
+ Scalar phi = FluidSystem::fugacityCoefficient(fs, paramCache, phaseIdx, compIdx);
+ fs.setFugacityCoefficient(phaseIdx, compIdx, phi);
+ }
+ }
+ }
+ else if (pvIdx < numPhases) { // <- first M - 1 saturations
+ Scalar delta = value - fs.saturation(/*phaseIdx=*/pvIdx - 1);
+ fs.setSaturation(/*phaseIdx=*/pvIdx - 1, value);
+
+ // set last saturation (-> minus the change of the
+ // satuation of the other phase)
+ fs.setSaturation(/*phaseIdx=*/numPhases - 1,
+ fs.saturation(numPhases - 1) - delta);
+
+ // update all fluid pressures using the capillary pressure
+ // law
+ ComponentVector pC;
+ MaterialLaw::capillaryPressures(pC, matParams, fs);
+ for (int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
+ fs.setPressure(phaseIdx,
+ fs.pressure(0)
+ + (pC[phaseIdx] - pC[0]));
+ paramCache.updateAllPressures(fs);
+
+ // update all densities and fugacity coefficients
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
+ fs.setDensity(phaseIdx, rho);
+
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
+ Scalar phi = FluidSystem::fugacityCoefficient(fs, paramCache, phaseIdx, compIdx);
+ fs.setFugacityCoefficient(phaseIdx, compIdx, phi);
+ }
+ }
+ }
+ else if (pvIdx < numPhases + numPhases*numComponents) // <- mole fractions
+ {
+ int phaseIdx = (pvIdx - numPhases)/numComponents;
+ int compIdx = (pvIdx - numPhases)%numComponents;
+
+ fs.setMoleFraction(phaseIdx, compIdx, value);
+ paramCache.updateSingleMoleFraction(fs, phaseIdx, compIdx);
+
+ // update the density of the phase
+ Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
+ fs.setDensity(phaseIdx, rho);
+
+ // if the phase's fugacity coefficients are composition
+ // dependent, update them as well.
+ if (!FluidSystem::isIdealMixture(phaseIdx)) {
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
+ Scalar phi = FluidSystem::fugacityCoefficient(fs, paramCache, phaseIdx, compIdx);
+ fs.setFugacityCoefficient(phaseIdx, compIdx, phi);
+ }
+ }
+ }
+ else {
+ assert(false);
+ }
+
+ // make the fluid state consistent with the fluid system.
+ completeFluidState_<MaterialLaw>(fs,
+ paramCache,
+ matParams);
+ };
+
+ // set a quantity in the fluid state
+ template <class FluidState>
+ static void setQuantityRaw_(FluidState &fs, int pvIdx, Scalar value)
+ {
+ assert(pvIdx < numEq);
+
+ // first pressure
+ if (pvIdx < 1) {
+ int phaseIdx = 0;
+ fs.setPressure(phaseIdx, value);
+ }
+ // first M - 1 saturations
+ else if (pvIdx < numPhases) {
+ int phaseIdx = pvIdx - 1;
+ fs.setSaturation(phaseIdx, value);
+ }
+ // mole fractions
+ else // if (pvIdx < numPhases + numPhases*numComponents)
+ {
+ int phaseIdx = (pvIdx - numPhases)/numComponents;
+ int compIdx = (pvIdx - numPhases)%numComponents;
+ fs.setMoleFraction(phaseIdx, compIdx, value);
+ }
+ };
+
+ template <class FluidState>
+ static Scalar quantityWeight_(const FluidState &fs, int pvIdx)
+ {
+ // first pressure
+ if (pvIdx < 1)
+ return 1.0/1e5;
+ // first M - 1 saturations
+ else if (pvIdx < numPhases)
+ return 1.0;
+ // mole fractions
+ else // if (pvIdx < numPhases + numPhases*numComponents)
+ return 1.0;
+ };
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh b/dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh
index d06a8abe36..a7ae7d4439 100644
--- a/dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh
+++ b/dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh
@@ -291,11 +291,11 @@ public:
Scalar T = fluidState.temperature(phaseIdx);
Scalar p = fluidState.pressure(phaseIdx);
-
- if (phaseIdx == lPhaseIdx)
- // assume pure water where one water molecule gets
- // replaced by one nitrogen molecule
+
+ if (phaseIdx == lPhaseIdx) {
+ // assume pure water
return H2O::liquidDensity(T, p);
+ }
// for the gas phase assume an ideal gas
return
@@ -470,8 +470,7 @@ public:
Valgrind::CheckDefined(p);
if (phaseIdx == lPhaseIdx) {
// TODO: correct way to deal with the solutes???
- return
- H2O::liquidEnthalpy(T, p) ;
+ return H2O::liquidEnthalpy(T, p);
}
else {
// assume ideal gas
@@ -503,7 +502,7 @@ public:
// Scalar p = fluidState.pressure(phaseIdx);
// Scalar T = fluidState.temperature(phaseIdx);
// Scalar x = fluidState.moleFrac(phaseIdx,compIdx);
-#warning: so far rough estimates from wikipedia
+//#warning: so far rough estimates from wikipedia
if (phaseIdx == lPhaseIdx)
return 0.6; // conductivity of water[W / (m K ) ]
@@ -524,7 +523,7 @@ public:
int phaseIdx)
{
// http://en.wikipedia.org/wiki/Heat_capacity
-#warning: so far rough estimates from wikipedia
+//#warning: so far rough estimates from wikipedia
// TODO heatCapacity is a function of composition.
// Scalar p = fluidState.pressure(phaseIdx);
// Scalar T = fluidState.temperature(phaseIdx);
diff --git a/dumux/material/MpNcfluidsystems/nullparametercache.hh b/dumux/material/MpNcfluidsystems/nullparametercache.hh
index 7efe8a584c..7952aa5b9a 100644
--- a/dumux/material/MpNcfluidsystems/nullparametercache.hh
+++ b/dumux/material/MpNcfluidsystems/nullparametercache.hh
@@ -37,18 +37,6 @@ class NullParameterCache : public ParameterCacheBase<NullParameterCache>
public:
NullParameterCache()
{};
-
- template <class FluidState>
- void updateAll(const FluidState &fs)
- {
- };
-
- /*!
- * \brief Update all cached parameters of a specific fluid phase
- */
- template <class FluidState>
- void updatePhase(const FluidState &fs, int phaseIdx)
- {};
};
} // end namepace
diff --git a/dumux/material/MpNcfluidsystems/parametercachebase.hh b/dumux/material/MpNcfluidsystems/parametercachebase.hh
index 1c381b4f41..2914ade4e1 100644
--- a/dumux/material/MpNcfluidsystems/parametercachebase.hh
+++ b/dumux/material/MpNcfluidsystems/parametercachebase.hh
@@ -34,21 +34,37 @@ template <class Implementation>
class ParameterCacheBase
{
public:
+ enum ExceptQuantities {
+ None = 0,
+ Temperature = 1,
+ Pressure = 2,
+ Composition = 2,
+ };
+
ParameterCacheBase()
{};
template <class FluidState>
- void updateAll(const FluidState &fs)
+ void updateAll(const FluidState &fs, int exceptQuantities = None)
{
for (int phaseIdx = 0; phaseIdx < FluidState::numPhases; ++phaseIdx)
updatePhase(fs, phaseIdx);
};
+
+ template <class FluidState>
+ void updateAllPressures(const FluidState &fs)
+ {
+ for (int phaseIdx = 0; phaseIdx < FluidState::numPhases; ++phaseIdx)
+ updatePhase(fs, phaseIdx);
+ };
+
+
/*!
* \brief Update all cached parameters of a specific fluid phase
*/
template <class FluidState>
- void updatePhase(const FluidState &fs, int phaseIdx)
+ void updatePhase(const FluidState &fs, int phaseIdx, int exceptQuantities = None)
{};
/*!
@@ -59,8 +75,8 @@ public:
* changed between two update*() calls. If more changed, call
* updatePhase()!
*/
- template <class FluidState> void
- updatePhaseTemperature(const FluidState &fs, int phaseIdx)
+ template <class FluidState>
+ void updateTemperature(const FluidState &fs, int phaseIdx)
{
asImp_().updatePhase(fs, phaseIdx);
};
@@ -74,7 +90,7 @@ public:
* updatePhase()!
*/
template <class FluidState>
- void updatePhasePressure(const FluidState &fs, int phaseIdx)
+ void updateSinglePressure(const FluidState &fs, int phaseIdx)
{
asImp_().updatePhase(fs, phaseIdx);
};
@@ -88,7 +104,7 @@ public:
* calls. If more changed, call updatePhase()!
*/
template <class FluidState>
- void updatePhaseComposition(const FluidState &fs, int phaseIdx)
+ void updateComposition(const FluidState &fs, int phaseIdx)
{
asImp_().updatePhase(fs, phaseIdx);
};
@@ -103,11 +119,11 @@ public:
* updatePhase()!
*/
template <class FluidState>
- void updatePhaseSingleMoleFraction(const FluidState &fs,
- int phaseIdx,
- int compIdx)
+ void updateSingleMoleFraction(const FluidState &fs,
+ int phaseIdx,
+ int compIdx)
{
- asImp_().updatePhaseComposition(fs, phaseIdx);
+ asImp_().updateComposition(fs, phaseIdx);
};
private:
diff --git a/test/material/Makefile.am b/test/material/Makefile.am
index e01d18a4d7..17f3316dcc 100644
--- a/test/material/Makefile.am
+++ b/test/material/Makefile.am
@@ -1,4 +1,4 @@
-SUBDIRS = tabulation .
+SUBDIRS = tabulation ncpflash
include $(top_srcdir)/am/global-rules
diff --git a/test/material/ncpflash/Makefile.am b/test/material/ncpflash/Makefile.am
new file mode 100644
index 0000000000..3c18c7f283
--- /dev/null
+++ b/test/material/ncpflash/Makefile.am
@@ -0,0 +1,5 @@
+check_PROGRAMS = test_ncpflash
+
+test_ncpflash_SOURCES = test_ncpflash.cc
+
+include $(top_srcdir)/am/global-rules
diff --git a/test/material/ncpflash/test_ncpflash.cc b/test/material/ncpflash/test_ncpflash.cc
new file mode 100644
index 0000000000..20fb52f37b
--- /dev/null
+++ b/test/material/ncpflash/test_ncpflash.cc
@@ -0,0 +1,294 @@
+/*****************************************************************************
+ * Copyright (C) 2011 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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 This is a program to test the flash calculation which uses
+ * non-linear complementarity problems (NCP)
+ *
+ * A flash calculation determines the pressures, saturations and
+ * composition of all phases given the total mass (or, as in this case
+ * the total number of moles) in a given amount of pore space.
+ */
+#include "config.h"
+
+#include <dumux/material/MpNcconstraintsolvers/misciblemultiphasecomposition.hh>
+#include <dumux/material/MpNcconstraintsolvers/computefromreferencephase.hh>
+#include <dumux/material/MpNcconstraintsolvers/ncpflash.hh>
+
+#include <dumux/material/MpNcfluidstates/compositionalfluidstate.hh>
+
+#include <dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh>
+
+#include <dumux/material/fluidmatrixinteractions/Mp/Mplinearmaterial.hh>
+#include <dumux/material/fluidmatrixinteractions/Mp/2padapter.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedlinearmaterial.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+
+template <class Scalar, class FluidState>
+void checkSame(const FluidState &fsRef, const FluidState &fsFlash)
+{
+ enum { numPhases = FluidState::numPhases };
+ enum { numComponents = FluidState::numComponents };
+
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ Scalar error;
+
+ // check the pressures
+ error = 1 - fsRef.pressure(phaseIdx)/fsFlash.pressure(phaseIdx);
+ if (std::abs(error) > 1e-6) {
+ std::cout << "pressure error phase " << phaseIdx << ": "
+ << fsFlash.pressure(phaseIdx) << " flash vs "
+ << fsRef.pressure(phaseIdx) << " reference"
+ << " error=" << error << "\n";
+ }
+
+ // check the saturations
+ error = fsRef.saturation(phaseIdx) - fsFlash.saturation(phaseIdx);
+ if (std::abs(error) > 1e-6)
+ std::cout << "saturation error phase " << phaseIdx << ": "
+ << fsFlash.saturation(phaseIdx) << " flash vs "
+ << fsRef.saturation(phaseIdx) << " reference"
+ << " error=" << error << "\n";
+
+ // check the compositions
+ for (int compIdx = 0; compIdx < numComponents; ++ compIdx) {
+ error = fsRef.moleFraction(phaseIdx, compIdx) - fsFlash.moleFraction(phaseIdx, compIdx);
+ if (std::abs(error) > 1e-6)
+ std::cout << "composition error phase " << phaseIdx << ", component " << compIdx << ": "
+ << fsFlash.moleFraction(phaseIdx, compIdx) << " flash vs "
+ << fsRef.moleFraction(phaseIdx, compIdx) << " reference"
+ << " error=" << error << "\n";
+ };
+ }
+}
+
+template <class Scalar, class FluidSystem, class MaterialLaw, class FluidState>
+void checkNcpFlash(const FluidState &fsRef,
+ typename MaterialLaw::Params &matParams)
+{
+ enum { numPhases = FluidSystem::numPhases };
+ enum { numComponents = FluidSystem::numComponents };
+ typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
+
+ // calculate the total amount of stuff in the reference fluid
+ // phase
+ ComponentVector globalMolarities(0.0);
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ globalMolarities[compIdx] +=
+ fsRef.saturation(phaseIdx)*fsRef.molarity(phaseIdx, compIdx);
+ }
+ }
+
+ // initialize the fluid state for the flash calculation
+ typedef Dumux::NcpFlash<Scalar, FluidSystem> NcpFlash;
+ FluidState fsFlash;
+
+ fsFlash.setTemperature(fsRef.temperature(/*phaseIdx=*/0));
+
+ // run the flash calculation
+ typename FluidSystem::ParameterCache paramCache;
+ NcpFlash::guessInitial(fsFlash, paramCache, globalMolarities);
+ NcpFlash::template solve<MaterialLaw>(fsFlash, paramCache, matParams, globalMolarities);
+
+ // compare the "flashed" fluid state with the reference one
+ checkSame<Scalar>(fsRef, fsFlash);
+};
+
+
+template <class Scalar, class FluidSystem, class MaterialLaw, class FluidState>
+void completeReferenceFluidState(FluidState &fs,
+ typename MaterialLaw::Params &matParams,
+ int refPhaseIdx)
+{
+ enum { numPhases = FluidSystem::numPhases };
+
+ typedef Dumux::ComputeFromReferencePhase<Scalar, FluidSystem> ComputeFromReferencePhase;
+ typedef Dune::FieldVector<Scalar, numPhases> PhaseVector;
+
+ int otherPhaseIdx = 1 - refPhaseIdx;
+
+ // calculate the other saturation
+ fs.setSaturation(otherPhaseIdx, 1.0 - fs.saturation(refPhaseIdx));
+
+ // calulate the capillary pressure
+ PhaseVector pC;
+ MaterialLaw::capillaryPressures(pC, matParams, fs);
+ fs.setPressure(otherPhaseIdx,
+ fs.pressure(refPhaseIdx)
+ + (pC[otherPhaseIdx] - pC[refPhaseIdx]));
+
+ // make the fluid state consistent with local thermodynamic
+ // equilibrium
+ typename FluidSystem::ParameterCache paramCache;
+ ComputeFromReferencePhase::solve(fs,
+ paramCache,
+ refPhaseIdx,
+ /*setViscosity=*/false,
+ /*setEnthalpy=*/false);
+}
+
+
+int main()
+{
+ typedef double Scalar;
+ typedef Dumux::H2ON2FluidSystem<Scalar> FluidSystem;
+ typedef Dumux::CompositionalFluidState<Scalar, FluidSystem> CompositionalFluidState;
+
+ enum { numPhases = FluidSystem::numPhases };
+ enum { numComponents = FluidSystem::numComponents };
+ enum { lPhaseIdx = FluidSystem::lPhaseIdx };
+ enum { gPhaseIdx = FluidSystem::gPhaseIdx };
+
+ enum { H2OIdx = FluidSystem::H2OIdx };
+ enum { N2Idx = FluidSystem::N2Idx };
+
+ typedef Dumux::RegularizedBrooksCorey<Scalar> EffMaterialLaw;
+ typedef Dumux::EffToAbsLaw<EffMaterialLaw> TwoPMaterialLaw;
+ typedef Dumux::TwoPAdapter<lPhaseIdx, TwoPMaterialLaw> MaterialLaw;
+ typedef MaterialLaw::Params MaterialLawParams;
+
+ Scalar T = 273.15 + 25;
+
+ // initialize the tables of the fluid system
+ Scalar Tmin = T - 1.0;
+ Scalar Tmax = T + 1.0;
+ int nT = 3;
+
+ Scalar pmin = 0.75 * 1e5;
+ Scalar pmax = 1.25 * 2e5;
+ int np = 100;
+
+ FluidSystem::init(Tmin, Tmax, nT, pmin, pmax, np);
+
+ // set the parameters for the capillary pressure law
+ MaterialLawParams matParams;
+ matParams.setSwr(0.0);
+ matParams.setSnr(0.0);
+ matParams.setPe(0);
+ matParams.setLambda(2.0);
+
+ CompositionalFluidState fsRef;
+
+ // create an fluid state which is consistent
+
+ // set the fluid temperatures
+ fsRef.setTemperature(T);
+
+ ////////////////
+ // only liquid
+ ////////////////
+
+ // set liquid saturation
+ fsRef.setSaturation(lPhaseIdx, 1.0);
+
+ // set pressure of the liquid phase
+ fsRef.setPressure(lPhaseIdx, 2e5);
+
+ // set the liquid composition to pure water
+ fsRef.setMoleFraction(lPhaseIdx, N2Idx, 0.0);
+ fsRef.setMoleFraction(lPhaseIdx, H2OIdx, 1.0);
+
+ // "complete" the fluid state
+ completeReferenceFluidState<Scalar, FluidSystem, MaterialLaw>(fsRef, matParams, lPhaseIdx);
+
+ // check the flash calculation
+ checkNcpFlash<Scalar, FluidSystem, MaterialLaw>(fsRef, matParams);
+
+ ////////////////
+ // only gas
+ ////////////////
+
+ // set gas saturation
+ fsRef.setSaturation(gPhaseIdx, 1.0);
+
+ // set pressure of the gas phase
+ fsRef.setPressure(gPhaseIdx, 1e6);
+
+ // set the gas composition to 99.9% nitrogen and 0.1% water
+ fsRef.setMoleFraction(gPhaseIdx, N2Idx, 0.999);
+ fsRef.setMoleFraction(gPhaseIdx, H2OIdx, 0.001);
+
+ // "complete" the fluid state
+ completeReferenceFluidState<Scalar, FluidSystem, MaterialLaw>(fsRef, matParams, gPhaseIdx);
+
+ // check the flash calculation
+ checkNcpFlash<Scalar, FluidSystem, MaterialLaw>(fsRef, matParams);
+
+ ////////////////
+ // both pgases gas
+ ////////////////
+
+ // set saturations
+ fsRef.setSaturation(lPhaseIdx, 0.5);
+ fsRef.setSaturation(gPhaseIdx, 0.5);
+
+ // set pressures
+ fsRef.setPressure(lPhaseIdx, 1e6);
+ fsRef.setPressure(gPhaseIdx, 1e6);
+
+ FluidSystem::ParameterCache paramCache;
+ typedef Dumux::MiscibleMultiPhaseComposition<Scalar, FluidSystem> MiscibleMultiPhaseComposition;
+ MiscibleMultiPhaseComposition::solve(fsRef, paramCache,
+ /*setViscosity=*/false,
+ /*setEnthalpy=*/false);
+
+ // check the flash calculation
+ checkNcpFlash<Scalar, FluidSystem, MaterialLaw>(fsRef, matParams);
+
+ ////////////////
+ // with capillary pressure
+ ////////////////
+
+ MaterialLawParams matParams2;
+ matParams2.setSwr(0.0);
+ matParams2.setSnr(0.0);
+ matParams2.setPe(1e3);
+ matParams2.setLambda(2.0);
+
+ // set gas saturation
+ fsRef.setSaturation(gPhaseIdx, 0.5);
+ fsRef.setSaturation(lPhaseIdx, 0.5);
+
+ // set pressure of the liquid phase
+ fsRef.setPressure(lPhaseIdx, 1e6);
+
+ // calulate the capillary pressure
+ typedef Dune::FieldVector<Scalar, numPhases> PhaseVector;
+ PhaseVector pC;
+ MaterialLaw::capillaryPressures(pC, matParams2, fsRef);
+ fsRef.setPressure(gPhaseIdx,
+ fsRef.pressure(lPhaseIdx)
+ + (pC[gPhaseIdx] - pC[lPhaseIdx]));
+
+ typedef Dumux::MiscibleMultiPhaseComposition<Scalar, FluidSystem> MiscibleMultiPhaseComposition;
+ MiscibleMultiPhaseComposition::solve(fsRef, paramCache,
+ /*setViscosity=*/false,
+ /*setEnthalpy=*/false);
+
+
+ // check the flash calculation
+ checkNcpFlash<Scalar, FluidSystem, MaterialLaw>(fsRef, matParams2);
+
+ return 0;
+}
--
GitLab