Commit 916f27c9 authored by Andreas Lauser's avatar Andreas Lauser
Browse files

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
parent fb0b3a75
......@@ -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
])
......
......@@ -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
}
/*!
......
......@@ -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
......
......@@ -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 {
/*!
......
......@@ -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);
......
This diff is collapsed.
......@@ -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);
......
......@@ -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
......
......@@ -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:
......
SUBDIRS = tabulation .
SUBDIRS = tabulation ncpflash
include $(top_srcdir)/am/global-rules
check_PROGRAMS = test_ncpflash
test_ncpflash_SOURCES = test_ncpflash.cc
include $(top_srcdir)/am/global-rules
/*****************************************************************************
* 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;
}
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