Skip to content
Snippets Groups Projects
Commit 79a23fcd authored by Andreas Lauser's avatar Andreas Lauser
Browse files

richards box model: eliminate custom fluid state 

was replaced by ImmiscibleFluidState

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@7056 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 8bb2a761
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dumux/boxmodels/richards/richardsfluidstate.hh deleted 100644 → 0
+ 0
− 236
View file @ 8bb2a761
/*****************************************************************************
* Copyright (C) 2009 by Markus Wolff *
* 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 Calcultes the fluid state from the primary variables in the
* Richards model.
*/
#ifndef DUMUX_RICHARDS_FLUID_STATE_HH
#define DUMUX_RICHARDS_FLUID_STATE_HH
#include <dumux/material/fluidstate.hh>
#include "richardsproperties.hh"
namespace Dumux
{
/*!
* \brief Calcultes the fluid state from the primary variables in the
* Richards model.
*/
template <class TypeTag>
class RichardsFluidState : public FluidState<typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)),
RichardsFluidState<TypeTag> >
{
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLaw)) MaterialLaw;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLawParams)) MaterialLawParams;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(RichardsIndices)) Indices;
enum {
wPhaseIdx = Indices::wPhaseIdx,
nPhaseIdx = Indices::nPhaseIdx,
pwIdx= Indices::pwIdx,
};
public:
//! The number of fluid phases used by the fluid state
enum { numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)) };
/*!
* \brief Updates the fluid quantities from the primary variables
* of the Richards model.
*
* \param pnRef The reference pressure of the non-wetting fluid phase \f$\mathrm{[Pa]}\f$
* \param matParams The parameters for the capillary pressure/relative permeability material law
* \param priVars The primary variables for which the fluid state ought to be calculated
* \param temperature The temperature which should be used
*/
void update(Scalar pnRef,
const MaterialLawParams &matParams,
const PrimaryVariables &priVars,
Scalar temperature)
{
Scalar minPc = MaterialLaw::pC(matParams, 1.0);
pressure_[wPhaseIdx] = priVars[pwIdx];
pressure_[nPhaseIdx] = std::max(pnRef, priVars[pwIdx] + minPc);
Sn_ = 1.0 - MaterialLaw::Sw(matParams, capillaryPressure());
temperature_=temperature;
typename FluidSystem::ParameterCache paramCache;
paramCache.updateAll(*this);
densityWetting_ = FluidSystem::density(*this, paramCache, wPhaseIdx);
viscosityWetting_ = FluidSystem::viscosity(*this, paramCache, wPhaseIdx);
}
/*!
* \brief Returns the saturation [] of a phase.
*
* \param phaseIdx The index of the fluid phase
*/
Scalar saturation(int phaseIdx) const
{
if (phaseIdx == wPhaseIdx)
return Scalar(1.0) - Sn_;
else
return Sn_;
};
/*!
* \brief Returns the mass fraction [] of a component in a phase.
*
* \param phaseIdx The index of the fluid phase
* \param compIdx The index of the chemical species
*/
Scalar massFraction(int phaseIdx, int compIdx) const
{
if (compIdx == phaseIdx)
return 1.0;
return 0;
}
/*!
* \brief Returns the molar fraction [] of a component in a fluid phase.
*
* \param phaseIdx The index of the fluid phase
* \param compIdx The index of the chemical species
*/
Scalar moleFraction(int phaseIdx, int compIdx) const
{
return massFraction(phaseIdx, compIdx);
}
/*!
* \brief Returns the molar density of a phase \f$\mathrm{[mol/m^3]}\f$.
*
* \param phaseIdx The index of the fluid phase
* \param compIdx The index of the chemical species
*/
Scalar molarDensity(int phaseIdx, int compIdx) const
{
return density(phaseIdx)/FluidSystem::molarMass(compIdx);
};
/*!
* \brief Returns the molar density of a phase \f$\mathrm{[mol/m^3]}\f$.
*
* \param phaseIdx The index of the fluid phase
* \param compIdx The index of the chemical species
*/
Scalar molarity(int phaseIdx, int compIdx) const
{
if (phaseIdx != compIdx)
return 0;
return molarDensity(phaseIdx);
}
/*!
* \brief Returns the density of a phase \f$\mathrm{[kg/m^3]}\f$.
*
* \param phaseIdx The index of the fluid phase
*/
Scalar density(int phaseIdx) const
{
if (phaseIdx == wPhaseIdx)
return densityWetting_;
return 1e-10;
}
/*!
* \brief Returns the dynamic viscosity of a phase [Pa s].
*
* \param phaseIdx The index of the fluid phase
*/
Scalar viscosity(int phaseIdx) const
{
if (phaseIdx == wPhaseIdx)
return viscosityWetting_;
return 1e-20; // something quite small
}
/*!
* \brief Returns mean molar mass of a phase \f$\mathrm{[kg/mol]}\f$.
*
* This is equivalent to the sum of all component molar masses
* weighted by their respective mole fraction.
*
* \param phaseIdx The index of the fluid phase
*/
Scalar averageMolarMass(int phaseIdx) const
{ return FluidSystem::molarMass(phaseIdx); };
/*!
* \brief Returns the partial pressure of a component in the gas phase \f$\mathrm{[Pa]}\f$.
*
* \param compIdx The index of the chemical species
*/
Scalar fugacity(int compIdx) const
{
if (compIdx == wPhaseIdx)
return 0;
return pressure_[nPhaseIdx];
}
/*!
* \brief Returns the pressure of a fluid phase \f$\mathrm{[Pa]}\f$.
*
* \param phaseIdx The index of the fluid phase
*/
Scalar pressure(int phaseIdx) const
{ return pressure_[phaseIdx]; }
/*!
* \brief Returns the capillary pressure \f$\mathrm{[Pa]}\f$
*/
Scalar capillaryPressure() const
{ return pressure_[nPhaseIdx] - pressure_[wPhaseIdx]; }
/*!
* \brief Returns the temperature a single fluid \f$\mathrm{[K]}\f$.
*/
Scalar temperature(int phaseIdx) const
{ return temperature_; };
/*!
* \brief Returns the temperature of the fluids \f$\mathrm{[K]}\f$.
*
* Note that we assume thermodynamic equilibrium, so all fluids
* and the rock matrix exhibit the same temperature.
*/
Scalar temperature() const
{ return temperature_; };
public:
Scalar pressure_[numPhases];
Scalar densityWetting_;
Scalar viscosityWetting_;
Scalar temperature_;
Scalar Sn_;
};
} // end namepace
#endif
dumux/boxmodels/richards/richardsproperties.hh
+ 0
− 1
View file @ 79a23fcd
...@@ -58,7 +58,6 @@ NEW_PROP_TAG(MaterialLawParams); //!< The type of the parameter object for the m ...@@ -58,7 +58,6 @@ NEW_PROP_TAG(MaterialLawParams); //!< The type of the parameter object for the m
NEW_PROP_TAG(FluidSystem); //!< The fluid system to be used for the Richards model NEW_PROP_TAG(FluidSystem); //!< The fluid system to be used for the Richards model
NEW_PROP_TAG(WettingPhase); //!< Fluid which represents the wetting phase NEW_PROP_TAG(WettingPhase); //!< Fluid which represents the wetting phase
NEW_PROP_TAG(NonWettingPhase); //!< Fluid which represents the non-wetting phase NEW_PROP_TAG(NonWettingPhase); //!< Fluid which represents the non-wetting phase
NEW_PROP_TAG(FluidState); //!< The fluid state to be used for the Richards model
NEW_PROP_TAG(EnableGravity); //!< Returns whether gravity is considered in the problem NEW_PROP_TAG(EnableGravity); //!< Returns whether gravity is considered in the problem
NEW_PROP_TAG(MassUpwindWeight); //!< The value of the weight of the upwind direction in the mass conservation equations NEW_PROP_TAG(MassUpwindWeight); //!< The value of the weight of the upwind direction in the mass conservation equations
// \} // \}
......
dumux/boxmodels/richards/richardspropertydefaults.hh
+ 0
− 4
View file @ 79a23fcd
...@@ -34,7 +34,6 @@ ...@@ -34,7 +34,6 @@
#include "richardsindices.hh" #include "richardsindices.hh"
#include "richardsfluxvariables.hh" #include "richardsfluxvariables.hh"
#include "richardsvolumevariables.hh" #include "richardsvolumevariables.hh"
#include "richardsfluidstate.hh"
#include "richardsproperties.hh" #include "richardsproperties.hh"
#include "richardsnewtoncontroller.hh" #include "richardsnewtoncontroller.hh"
...@@ -146,9 +145,6 @@ public: ...@@ -146,9 +145,6 @@ public:
NonWettingPhase> type; NonWettingPhase> type;
}; };
//! The fluid state class
SET_TYPE_PROP(BoxRichards, FluidState, RichardsFluidState<TypeTag>);
// \} // \}
}; };
......
dumux/boxmodels/richards/richardsvolumevariables.hh
+ 69
− 15
View file @ 79a23fcd
...@@ -27,7 +27,8 @@ ...@@ -27,7 +27,8 @@
#define DUMUX_RICHARDS_VOLUME_VARIABLES_HH #define DUMUX_RICHARDS_VOLUME_VARIABLES_HH
#include "richardsproperties.hh" #include "richardsproperties.hh"
#include "richardsfluidstate.hh"
#include <dumux/material/MpNcfluidstates/immisciblefluidstate.hh>
namespace Dumux namespace Dumux
{ {
...@@ -50,7 +51,6 @@ class RichardsVolumeVariables : public BoxVolumeVariables<TypeTag> ...@@ -50,7 +51,6 @@ class RichardsVolumeVariables : public BoxVolumeVariables<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView; typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem; typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem; typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidState)) FluidState;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLaw)) MaterialLaw; typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLaw)) MaterialLaw;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLawParams)) MaterialLawParams; typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLawParams)) MaterialLawParams;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry; typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
...@@ -70,6 +70,9 @@ class RichardsVolumeVariables : public BoxVolumeVariables<TypeTag> ...@@ -70,6 +70,9 @@ class RichardsVolumeVariables : public BoxVolumeVariables<TypeTag>
typedef typename GridView::template Codim<0>::Entity Element; typedef typename GridView::template Codim<0>::Entity Element;
public: public:
//! The type returned by the fluidState() method
typedef Dumux::ImmiscibleFluidState<Scalar, FluidSystem> FluidState;
/*! /*!
* \brief Update all quantities for a given control volume. * \brief Update all quantities for a given control volume.
* *
...@@ -99,23 +102,55 @@ public: ...@@ -99,23 +102,55 @@ public:
scvIdx, scvIdx,
isOldSol); isOldSol);
this->asImp_().updateTemperature_(priVars, asImp_().updateTemperature_(priVars,
element, element,
elemGeom, elemGeom,
scvIdx, scvIdx,
problem); problem);
//////////
// specify the fluid state
//////////
// material law parameters // pressures
Scalar pnRef = problem.referencePressure(element, elemGeom, scvIdx); Scalar pnRef = problem.referencePressure(element, elemGeom, scvIdx);
const MaterialLawParams &materialParams = const MaterialLawParams &matParams =
problem.spatialParameters().materialLawParams(element, elemGeom, scvIdx); problem.spatialParameters().materialLawParams(element, elemGeom, scvIdx);
fluidState_.update(pnRef, materialParams, priVars, temperature_); Scalar minPc = MaterialLaw::pC(matParams, 1.0);
fluidState_.setPressure(wPhaseIdx, priVars[pwIdx]);
fluidState_.setPressure(nPhaseIdx, std::max(pnRef, priVars[pwIdx] + minPc));
// saturations
Scalar Sw = MaterialLaw::Sw(matParams, capillaryPressure());
fluidState_.setSaturation(wPhaseIdx, Sw);
fluidState_.setSaturation(nPhaseIdx, 1 - Sw);
relativePermeabilityWetting_ = // density and viscosity
MaterialLaw::krw(materialParams, fluidState_.saturation(wPhaseIdx)); typename FluidSystem::ParameterCache paramCache;
paramCache.updateAll(fluidState_);
fluidState_.setDensity(wPhaseIdx, FluidSystem::density(fluidState_, paramCache, wPhaseIdx));
fluidState_.setDensity(nPhaseIdx, 1e-10);
fluidState_.setViscosity(wPhaseIdx, FluidSystem::viscosity(fluidState_, paramCache, wPhaseIdx));
fluidState_.setViscosity(nPhaseIdx, 1e-10);
//////////
// specify the other parameters
//////////
relativePermeabilityWetting_ = MaterialLaw::krw(matParams, Sw);
porosity_ = problem.spatialParameters().porosity(element, elemGeom, scvIdx); porosity_ = problem.spatialParameters().porosity(element, elemGeom, scvIdx);
// energy related quantities
asImp_().updateEnergy_(paramCache, priVars, problem, element, elemGeom, scvIdx, isOldSol);
} }
/*!
* \brief Return the fluid configuration at the given primary
* variables
*/
const FluidState &fluidState() const
{ return fluidState_; }
/*! /*!
* \brief Returns the average porosity [] within the control volume. * \brief Returns the average porosity [] within the control volume.
...@@ -210,7 +245,7 @@ public: ...@@ -210,7 +245,7 @@ public:
* \f[ p_c = p_n - p_w \f] * \f[ p_c = p_n - p_w \f]
*/ */
Scalar capillaryPressure() const Scalar capillaryPressure() const
{ return fluidState_.capillaryPressure(); } { return fluidState_.pressure(nPhaseIdx) - fluidState_.pressure(wPhaseIdx); }
protected: protected:
void updateTemperature_(const PrimaryVariables &priVars, void updateTemperature_(const PrimaryVariables &priVars,
...@@ -219,13 +254,32 @@ protected: ...@@ -219,13 +254,32 @@ protected:
int scvIdx, int scvIdx,
const Problem &problem) const Problem &problem)
{ {
temperature_ = problem.boxTemperature(element, elemGeom, scvIdx); fluidState_.setTemperature(problem.boxTemperature(element, elemGeom, scvIdx));
} }
/*!
* \brief Called by update() to compute the energy related quantities
*/
template <class ParameterCache>
void updateEnergy_(ParameterCache &paramCache,
const PrimaryVariables &sol,
const Problem &problem,
const Element &element,
const FVElementGeometry &elemGeom,
int vertIdx,
bool isOldSol)
{ }
FluidState fluidState_; FluidState fluidState_;
Scalar relativePermeabilityWetting_; Scalar relativePermeabilityWetting_;
Scalar porosity_; Scalar porosity_;
Scalar temperature_;
private:
Implementation &asImp_()
{ return *static_cast<Implementation*>(this); }
const Implementation &asImp_() const
{ return *static_cast<const Implementation*>(this); }
}; };
} }
......
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