Commit 22c4daa9 authored by Andreas Lauser's avatar Andreas Lauser
Browse files

Revert "convert the 2p2c(ni) models to the new fluid systems"

This reverts commit f7086e6326a239654762b525ac3f0692ad0cc83b.

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@4955 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 18109017
......@@ -26,7 +26,6 @@
* \brief Calculates the phase state from the primary variables in the
* 2p2c model.
*/
#if 0
#ifndef DUMUX_2P2C_PHASE_STATE_HH
#define DUMUX_2P2C_PHASE_STATE_HH
......@@ -368,4 +367,3 @@ public:
} // end namepace
#endif
#endif
......@@ -244,7 +244,6 @@ private:
const Element &element,
const ElementVolumeVariables &elemDat)
{
#if 0
const VolumeVariables &vDat_i = elemDat[face().i];
const VolumeVariables &vDat_j = elemDat[face().j];
......@@ -261,6 +260,7 @@ private:
// calculate tortuosity at the nodes i and j needed
// for porous media diffusion coefficient
Scalar tau_i =
1.0/(vDat_i.porosity() * vDat_i.porosity()) *
pow(vDat_i.porosity() * vDat_i.saturation(phaseIdx), 7.0/3);
......@@ -272,8 +272,8 @@ private:
// -> harmonic mean
porousDiffCoeff_[phaseIdx] = harmonicMean(vDat_i.porosity() * vDat_i.saturation(phaseIdx) * tau_i * vDat_i.diffCoeff(phaseIdx),
vDat_j.porosity() * vDat_j.saturation(phaseIdx) * tau_j * vDat_j.diffCoeff(phaseIdx));
}
#endif
}
public:
......@@ -306,13 +306,11 @@ public:
int downstreamIdx(int phaseIdx) const
{ return downstreamIdx_[phaseIdx]; }
#if 0
/*!
* \brief The binary diffusion coefficient for each fluid phase.
*/
Scalar porousDiffCoeff(int phaseIdx) const
{ return porousDiffCoeff_[phaseIdx]; };
#endif
/*!
* \brief Return density \f$\mathrm{[kg/m^3]}\f$ of a phase at the integration
......@@ -364,10 +362,8 @@ protected:
// local index of the downwind vertex for each phase
int downstreamIdx_[numPhases];
/*
// the diffusion coefficient for the porous medium
Scalar porousDiffCoeff_[numPhases];
*/
};
} // end namepace
......
......@@ -74,6 +74,8 @@ protected:
typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
typedef TwoPTwoCFluidState<TypeTag> FluidState;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices)) Indices;
enum
......@@ -205,9 +207,7 @@ public:
flux = 0;
asImp_()->computeAdvectiveFlux(flux, vars);
Valgrind::CheckDefined(flux);
asImp_()->computeDiffusiveFlux(flux, vars);
Valgrind::CheckDefined(flux);
}
/*!
......@@ -248,13 +248,6 @@ public:
- mobilityUpwindAlpha) * (dn.density(phaseIdx)
* dn.mobility(phaseIdx) * dn.fluidState().massFrac(
phaseIdx, compIdx));
Valgrind::CheckDefined(vars.KmvpNormal(phaseIdx));
Valgrind::CheckDefined(up.density(phaseIdx));
Valgrind::CheckDefined(up.mobility(phaseIdx));
Valgrind::CheckDefined(up.fluidState().massFrac(phaseIdx, compIdx));
Valgrind::CheckDefined(dn.density(phaseIdx));
Valgrind::CheckDefined(dn.mobility(phaseIdx));
Valgrind::CheckDefined(dn.fluidState().massFrac(phaseIdx, compIdx));
}
}
}
......@@ -268,7 +261,6 @@ public:
*/
void computeDiffusiveFlux(PrimaryVariables &flux, const FluxVariables &vars) const
{
#if 0
// add diffusive flux of gas component in liquid phase
Scalar tmp =
- vars.porousDiffCoeff(lPhaseIdx) *
......@@ -284,7 +276,6 @@ public:
(vars.molarConcGrad(gPhaseIdx) * vars.face().normal);
flux[contiLEqIdx] += tmp * FluidSystem::molarMass(lCompIdx);
flux[contiGEqIdx] -= tmp * FluidSystem::molarMass(gCompIdx);
#endif
}
/*!
......
......@@ -142,6 +142,8 @@ class TwoPTwoCModel: public BoxModel<TypeTag>
formulation = GET_PROP_VALUE(TypeTag, PTAG(Formulation))
};
typedef TwoPTwoCFluidState<TypeTag> FluidState;
typedef typename GridView::template Codim<dim>::Entity Vertex;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
......
......@@ -33,15 +33,11 @@
#include <dumux/common/math.hh>
#include <dune/common/collectivecommunication.hh>
#include <dumux/material/fluidstates/equilibriumfluidstate.hh>
#include <dumux/material/constraintsolvers/compositionfromfugacities.hh>
#include <vector>
#include <iostream>
#include "2p2cproperties.hh"
#include "2p2cindices.hh"
#include "2p2cfluidstate.hh"
namespace Dumux
{
......@@ -73,34 +69,17 @@ class TwoPTwoCVolumeVariables : public BoxVolumeVariables<TypeTag>
numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)),
formulation = GET_PROP_VALUE(TypeTag, PTAG(Formulation)),
// component indices
lCompIdx = Indices::lCompIdx,
gCompIdx = Indices::gCompIdx,
// phase indices
lPhaseIdx = Indices::lPhaseIdx,
gPhaseIdx = Indices::gPhaseIdx,
// primary variable indices
pressureIdx = Indices::pressureIdx,
switchIdx = Indices::switchIdx,
// phase states
lPhaseOnly = Indices::lPhaseOnly,
gPhaseOnly = Indices::gPhaseOnly,
bothPhases = Indices::bothPhases,
// formulations
plSg = TwoPTwoCFormulation::plSg,
pgSl = TwoPTwoCFormulation::pgSl,
gPhaseIdx = Indices::gPhaseIdx
};
typedef typename GridView::template Codim<0>::Entity Element;
typedef TwoPTwoCFluidState<TypeTag> FluidState;
public:
//! The return type of the fluidState() method
typedef Dumux::EquilibriumFluidState<Scalar, FluidSystem> FluidState;
/*!
* \brief Update all quantities for a given control volume.
*
......@@ -125,205 +104,40 @@ public:
scvIdx,
isOldSol);
typename FluidSystem::MutableParameters mutParams;
typename FluidSystem::MutableParameters::FluidState &fs
= mutParams.fluidState();
int globalVertIdx = problem.model().vertexMapper().map(element, scvIdx, dim);
int phasePresence = problem.model().phasePresence(globalVertIdx, isOldSol);
/////////////
// set the phase saturations
/////////////
if (phasePresence == gPhaseOnly) {
fs.setSaturation(lPhaseIdx, 0.0);
fs.setSaturation(gPhaseIdx, 1.0);
}
else if (phasePresence == lPhaseOnly) {
fs.setSaturation(lPhaseIdx, 1.0);
fs.setSaturation(gPhaseIdx, 0.0);
}
else if (phasePresence == bothPhases) {
Scalar Sg;
if (formulation == plSg)
Sg = priVars[switchIdx];
else if (formulation == pgSl)
Sg = 1.0 - priVars[switchIdx];
fs.setSaturation(lPhaseIdx, 1 - Sg);
fs.setSaturation(gPhaseIdx, Sg);
}
/////////////
// set the phase temperatures
/////////////
// update the temperature part of the energy module
Scalar T = asImp_().getTemperature(priVars,
asImp().updateTemperature_(priVars,
element,
elemGeom,
scvIdx,
problem);
Valgrind::CheckDefined(T);
for (int i = 0; i < numPhases; ++i)
fs.setTemperature(i, T);
/////////////
// set the phase pressures
/////////////
// capillary pressure parameters
const MaterialLawParams &materialParams =
problem.spatialParameters().materialLawParams(element, elemGeom, scvIdx);
Scalar pC = MaterialLaw::pC(materialParams, fs.saturation(lPhaseIdx));
if (formulation == plSg) {
fs.setPressure(lPhaseIdx, priVars[pressureIdx]);
fs.setPressure(gPhaseIdx, priVars[pressureIdx] + pC);
}
else if (formulation == pgSl){
fs.setPressure(lPhaseIdx, priVars[pressureIdx] - pC);
fs.setPressure(gPhaseIdx, priVars[pressureIdx]);
}
Valgrind::CheckDefined(fs.pressure(lPhaseIdx));
Valgrind::CheckDefined(fs.pressure(gPhaseIdx));
// update the mutable parameters for the pure components
mutParams.updatePure(lPhaseIdx);
mutParams.updatePure(gPhaseIdx);
/////////////
// set the phase compositions.
/////////////
if (phasePresence == gPhaseOnly) {
// mass fractions
Scalar Xg1 = priVars[switchIdx];
Scalar Xg2 = 1 - Xg1;
// molar masses
Scalar M1 = FluidSystem::molarMass(lCompIdx);
Scalar M2 = FluidSystem::molarMass(gCompIdx);
// convert mass to mole fractions
Scalar meanM = M1*M2/(M2 + Xg2*(M1 - M2));
fs.setMoleFrac(gPhaseIdx, lCompIdx, Xg1 * M1/meanM);
fs.setMoleFrac(gPhaseIdx, gCompIdx, Xg2 * M2/meanM);
mutParams.updateMix(gPhaseIdx);
// calculate component fugacities in gas phase
Scalar fug1 = FluidSystem::computeFugacity(mutParams, gPhaseIdx, lCompIdx);
Scalar fug2 = FluidSystem::computeFugacity(mutParams, gPhaseIdx, gCompIdx);
fs.setFugacity(gPhaseIdx, lCompIdx, fug1);
fs.setFugacity(gPhaseIdx, gCompIdx, fug2);
// use same fugacities in liquid phase
fs.setFugacity(lPhaseIdx, lCompIdx, fug1);
fs.setFugacity(lPhaseIdx, gCompIdx, fug2);
// initial guess of liquid composition
fs.setMoleFrac(lPhaseIdx, lCompIdx, 0.98);
fs.setMoleFrac(lPhaseIdx, gCompIdx, 0.02);
// calculate liquid composition from fugacities
typedef Dumux::CompositionFromFugacities<Scalar, FluidSystem> CompFromFug;
CompFromFug::run(mutParams, lPhaseIdx);
Valgrind::CheckDefined(fs.moleFrac(gPhaseIdx, lCompIdx));
Valgrind::CheckDefined(fs.moleFrac(gPhaseIdx, gCompIdx));
Valgrind::CheckDefined(fs.moleFrac(lPhaseIdx, lCompIdx));
Valgrind::CheckDefined(fs.moleFrac(lPhaseIdx, gCompIdx));
// calculate molar volume of gas phase
Scalar Vmg = FluidSystem::computeMolarVolume(mutParams, gPhaseIdx);
fs.setMolarVolume(gPhaseIdx, Vmg);
}
else if (phasePresence == lPhaseOnly) {
// mass fractions
Scalar Xl2 = priVars[switchIdx];
Scalar Xl1 = 1 - Xl2;
// molar masses
Scalar M1 = FluidSystem::molarMass(lCompIdx);
Scalar M2 = FluidSystem::molarMass(gCompIdx);
// convert mass to mole fractions
Scalar meanM = M1*M2/(M2 + Xl2*(M1 - M2));
fs.setMoleFrac(lPhaseIdx, lCompIdx, Xl1 * M1/meanM);
fs.setMoleFrac(lPhaseIdx, gCompIdx, Xl2 * M2/meanM);
mutParams.updateMix(lPhaseIdx);
// calculate component fugacities in liquid phase
Scalar fug1 = FluidSystem::computeFugacity(mutParams, lPhaseIdx, lCompIdx);
Scalar fug2 = FluidSystem::computeFugacity(mutParams, lPhaseIdx, gCompIdx);
fs.setFugacity(lPhaseIdx, lCompIdx, fug1);
fs.setFugacity(lPhaseIdx, gCompIdx, fug2);
// use same fugacities in gas phase
fs.setFugacity(gPhaseIdx, lCompIdx, fug1);
fs.setFugacity(gPhaseIdx, gCompIdx, fug2);
// initial guess of gas composition
fs.setMoleFrac(gPhaseIdx, lCompIdx, 0.05);
fs.setMoleFrac(gPhaseIdx, gCompIdx, 0.95);
// calculate liquid composition from fugacities
typedef Dumux::CompositionFromFugacities<Scalar, FluidSystem> CompFromFug;
CompFromFug::run(mutParams, gPhaseIdx);
Valgrind::CheckDefined(fs.moleFrac(gPhaseIdx, lCompIdx));
Valgrind::CheckDefined(fs.moleFrac(gPhaseIdx, gCompIdx));
Valgrind::CheckDefined(fs.moleFrac(lPhaseIdx, lCompIdx));
Valgrind::CheckDefined(fs.moleFrac(lPhaseIdx, gCompIdx));
// calculate molar volume of liquid phase
Scalar Vml = FluidSystem::computeMolarVolume(mutParams, lPhaseIdx);
fs.setMolarVolume(lPhaseIdx, Vml);
}
else if (phasePresence == bothPhases) {
// HACK: assume both phases to be an ideal mixture,
// i.e. the fugacity coefficents do not depend on the
// composition
Scalar phi_l1 = FluidSystem::computeFugacityCoeff(mutParams, lPhaseIdx, lCompIdx);
Scalar phi_l2 = FluidSystem::computeFugacityCoeff(mutParams, lPhaseIdx, gCompIdx);
Scalar phi_g1 = FluidSystem::computeFugacityCoeff(mutParams, gPhaseIdx, lCompIdx);
Scalar phi_g2 = FluidSystem::computeFugacityCoeff(mutParams, gPhaseIdx, gCompIdx);
Scalar pl = fs.pressure(lPhaseIdx);
Scalar pg = fs.pressure(gPhaseIdx);
Valgrind::CheckDefined(phi_l1);
Valgrind::CheckDefined(phi_l2);
Valgrind::CheckDefined(phi_g1);
Valgrind::CheckDefined(phi_g2);
Scalar xg2 = (phi_g2/phi_l1 - pl/pg) / (phi_g1/phi_l1 - phi_g2/phi_l2);
Scalar xg1 = 1 - xg2;
Scalar xl2 = (xg2*pg*phi_g2)/(pl*phi_l2);
Scalar xl1 = 1 - xl2;
fs.setMoleFrac(lPhaseIdx, lCompIdx, xl1);
fs.setMoleFrac(lPhaseIdx, gCompIdx, xl2);
fs.setMoleFrac(gPhaseIdx, lCompIdx, xg1);
fs.setMoleFrac(gPhaseIdx, gCompIdx, xg2);
mutParams.updateMix(lPhaseIdx);
mutParams.updateMix(gPhaseIdx);
Scalar Vml = FluidSystem::computeMolarVolume(mutParams, lPhaseIdx);
Scalar Vmg = FluidSystem::computeMolarVolume(mutParams, gPhaseIdx);
fs.setMolarVolume(lPhaseIdx, Vml);
fs.setMolarVolume(gPhaseIdx, Vmg);
}
int globalVertIdx = problem.model().dofMapper().map(element, scvIdx, dim);
int phasePresence = problem.model().phasePresence(globalVertIdx, isOldSol);
// calculate phase state
fluidState_.update(priVars, materialParams, temperature(), phasePresence);
Valgrind::CheckDefined(fluidState_);
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
if (this->fluidState().saturation(phaseIdx) != 0.0) {
// Mobilities
Scalar muL = FluidSystem::computeViscosity(mutParams, lPhaseIdx);
Scalar krL = MaterialLaw::krw(materialParams, fs.saturation(lPhaseIdx));
mobility_[lPhaseIdx] = krL / muL;
Valgrind::CheckDefined(mobility_[lPhaseIdx]);
// ATTENTION: krn requires the liquid saturation as parameter!
Scalar muG = FluidSystem::computeViscosity(mutParams, gPhaseIdx);
Scalar krG = MaterialLaw::krn(materialParams, fs.saturation(lPhaseIdx));
mobility_[gPhaseIdx] = krG / muG;
Valgrind::CheckDefined(mobility_[gPhaseIdx]);
#if 0
const Scalar mu =
FluidSystem::phaseViscosity(phaseIdx,
fluidState().temperature(),
fluidState().phasePressure(lPhaseIdx),
fluidState());
Scalar kr;
if (phaseIdx == lPhaseIdx)
kr = MaterialLaw::krw(materialParams, saturation(lPhaseIdx));
else // ATTENTION: krn requires the liquid saturation
// as parameter!
kr = MaterialLaw::krn(materialParams, saturation(lPhaseIdx));
mobility_[phaseIdx] = kr / mu;
Valgrind::CheckDefined(mobility_[phaseIdx]);
// binary diffusion coefficents
diffCoeff_[phaseIdx] =
FluidSystem::diffCoeff(phaseIdx,
......@@ -333,18 +147,18 @@ public:
fluidState_.phasePressure(phaseIdx),
fluidState_);
Valgrind::CheckDefined(diffCoeff_[phaseIdx]);
#endif
}
else {
mobility_[phaseIdx] = 0;
diffCoeff_[phaseIdx] = 0;
}
}
// porosity
porosity_ = problem.spatialParameters().porosity(element,
elemGeom,
scvIdx);
Valgrind::CheckDefined(porosity_);
asImp_().updateEnergy(mutParams, priVars, element, elemGeom, scvIdx, problem);
// assign the equilibrium fluid state from the generic one
fluidState_.assign(fs);
}
/*!
......@@ -378,7 +192,7 @@ public:
* \param phaseIdx The phase index
*/
Scalar molarDensity(int phaseIdx) const
{ return fluidState_.molarDensity(phaseIdx); }
{ return fluidState_.density(phaseIdx) / fluidState_.averageMolarMass(phaseIdx); }
/*!
* \brief Returns the effective pressure of a given phase within
......@@ -387,7 +201,7 @@ public:
* \param phaseIdx The phase index
*/
Scalar pressure(int phaseIdx) const
{ return fluidState_.pressure(phaseIdx); }
{ return fluidState_.phasePressure(phaseIdx); }
/*!
* \brief Returns temperature inside the sub-control volume.
......@@ -397,7 +211,7 @@ public:
* identical.
*/
Scalar temperature() const
{ return fluidState_.temperature(); }
{ return temperature_; }
/*!
* \brief Returns the effective mobility of a given phase within
......@@ -414,7 +228,7 @@ public:
* \brief Returns the effective capillary pressure within the control volume.
*/
Scalar capillaryPressure() const
{ return fluidState_.pressure(lPhaseIdx) - fluidState_.pressure(gPhaseIdx); }
{ return fluidState_.capillaryPressure(); }
/*!
* \brief Returns the average porosity within the control volume.
......@@ -422,35 +236,35 @@ public:
Scalar porosity() const
{ return porosity_; }
#if 0
/*!
* \brief Returns the binary diffusion coefficients for a phase
*/
Scalar diffCoeff(int phaseIdx) const
{ return diffCoeff_[phaseIdx]; }
#endif
protected:
Scalar getTemperature_(const PrimaryVariables &priVars,
void updateTemperature_(const PrimaryVariables &priVars,
const Element &element,
const FVElementGeometry &elemGeom,
int scvIdx,
const Problem &problem)
{
return problem.temperature(element, elemGeom, scvIdx);
temperature_ = problem.temperature(element, elemGeom, scvIdx);
}
Scalar temperature_; //!< Temperature within the control volume
Scalar porosity_; //!< Effective porosity within the control volume
Scalar mobility_[numPhases]; //!< Effective mobility within the control volume
// Scalar diffCoeff_[numPhases]; //!< Binary diffusion coefficients for the phases
Scalar diffCoeff_[numPhases]; //!< Binary diffusion coefficients for the phases
FluidState fluidState_;
private:
Implementation &asImp_()
Implementation &asImp()
{ return *static_cast<Implementation*>(this); }
const Implementation &asImp_() const
const Implementation &asImp() const
{ return *static_cast<const Implementation*>(this); }
};
......
......@@ -112,11 +112,13 @@ public:
// compute the energy storage
result[energyEqIdx] =
vertDat.porosity()*(vertDat.density(lPhaseIdx) *
vertDat.fluidState().internalEnergy(lPhaseIdx) *
vertDat.internalEnergy(lPhaseIdx) *
//vertDat.enthalpy(lPhaseIdx) *
vertDat.saturation(lPhaseIdx)
+
vertDat.density(gPhaseIdx) *
vertDat.fluidState().internalEnergy(gPhaseIdx) *
vertDat.internalEnergy(gPhaseIdx) *
//vertDat.enthalpy(gPhaseIdx) *
vertDat.saturation(gPhaseIdx))
+
vertDat.temperature()*vertDat.heatCapacity();
......@@ -150,12 +152,12 @@ public:
mobilityUpwindAlpha * // upstream vertex
( up.density(phase) *
up.mobility(phase) *
up.fluidState().enthalpy(phase))
up.enthalpy(phase))
+
(1-mobilityUpwindAlpha) * // downstream vertex
( dn.density(phase) *
dn.mobility(phase) *
dn.fluidState().enthalpy(phase)) );
dn.enthalpy(phase)) );
}
}
......
......@@ -70,19 +70,6 @@ class TwoPTwoCNIVolumeVariables : public TwoPTwoCVolumeVariables<TypeTag>
//! \endcond
public:
/*!
* \brief Update the temperature of the sub-control volume.
*/
Scalar getTemperature(const PrimaryVariables &sol,
const Element &element,
const FVElementGeometry &elemGeom,
int scvIdx,
const Problem &problem) const
{
// retrieve temperature from solution vector
return sol[temperatureIdx];
}
/*!
* \brief Update all quantities for a given control volume.
*
......@@ -93,26 +80,61 @@ public:
* \param vertIdx The local index of the SCV (sub-control volume)
* \param isOldSol Evaluate function with solution of current or previous time step
*/
template <class MutableParams>
void updateEnergy(MutableParams &mutParams,
const PrimaryVariables &sol,
void update(const PrimaryVariables &sol,
const Problem &problem,
const Element &element,
const FVElementGeometry &elemGeom,
int scvIdx,
const Problem &problem)
int vertIdx,
bool isOldSol)
{
heatCapacity_ =
problem.spatialParameters().heatCapacity(element, elemGeom, scvIdx);
Valgrind::CheckDefined(heatCapacity_);
// vertex update data for the mass balance
ParentType::update(sol,
problem,
element,
elemGeom,
vertIdx,
isOldSol);
// the internal energies and the enthalpies
typename MutableParams::FluidState &fs = mutParams.fluidState();
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
Scalar h =
FluidSystem::computeEnthalpy(mutParams, phaseIdx);
fs.setEnthalpy(phaseIdx, h);
if (this->fluidState().saturation(phaseIdx) != 0.0) {
enthalpy_[phaseIdx] =
FluidSystem::phaseEnthalpy(phaseIdx,
this->fluidState().temperature(),
this->fluidState().phasePressure(phaseIdx),
this->fluidState());
internalEnergy_[phaseIdx] =
FluidSystem::phaseInternalEnergy(phaseIdx,
this->fluidState().temperature(),