// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* Copyright (C) 2008-2010 by Andreas Lauser *
* Copyright (C) 2008-2009 by Klaus Mosthaf *
* 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 Definition of the spatial parameters for the water-air problem.
*/
#ifndef DUMUX_WATER_AIR_SPATIAL_PARAMETERS_HH
#define DUMUX_WATER_AIR_SPATIAL_PARAMETERS_HH
#include <dumux/material/spatialparameters/boxspatialparameters.hh>
#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
namespace Dumux
{
//forward declaration
template<class TypeTag>
class WaterAirSpatialParameters;
namespace Properties
{
// The spatial parameters TypeTag
NEW_TYPE_TAG(WaterAirSpatialParameters);
// Set the spatial parameters
SET_TYPE_PROP(WaterAirSpatialParameters, SpatialParameters, Dumux::WaterAirSpatialParameters<TypeTag>);
// Set the material Law
SET_PROP(WaterAirSpatialParameters, MaterialLaw)
{
private:
// define the material law which is parameterized by effective
// saturations
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef RegularizedBrooksCorey<Scalar> EffMaterialLaw;
public:
// define the material law parameterized by absolute saturations
typedef EffToAbsLaw<EffMaterialLaw> type;
};
}
/*!
* \ingroup TwoPTwoCNIModel
* \ingroup BoxTestProblems
* \brief Definition of the spatial parameters for the water-air problem
*/
template<class TypeTag>
class WaterAirSpatialParameters : public BoxSpatialParameters<TypeTag>
{
typedef BoxSpatialParameters<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename Grid::ctype CoordScalar;
enum {
dim=GridView::dimension,
dimWorld=GridView::dimensionworld
};
typedef typename GET_PROP_TYPE(TypeTag, TwoPTwoCIndices) Indices;
enum {
lPhaseIdx = Indices::lPhaseIdx,
gPhaseIdx = Indices::gPhaseIdx
};
typedef Dune::FieldVector<CoordScalar,dim> LocalPosition;
typedef Dune::FieldVector<CoordScalar,dimWorld> GlobalPosition;
typedef Dune::FieldVector<CoordScalar,dimWorld> Vector;
typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
public:
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
typedef typename MaterialLaw::Params MaterialLawParams;
/*!
* \brief The constructor
*
* \param gv The grid view
*/
WaterAirSpatialParameters(const GridView &gv)
: ParentType(gv)
{
layerBottom_ = 22.0;
// intrinsic permeabilities
fineK_ = 1e-13;
coarseK_ = 1e-12;
// porosities
finePorosity_ = 0.3;
coarsePorosity_ = 0.3;
// residual saturations
fineMaterialParams_.setSwr(0.2);
fineMaterialParams_.setSnr(0.0);
coarseMaterialParams_.setSwr(0.2);
coarseMaterialParams_.setSnr(0.0);
// parameters for the Brooks-Corey law
fineMaterialParams_.setPe(1e4);
coarseMaterialParams_.setPe(1e4);
fineMaterialParams_.setLambda(2.0);
coarseMaterialParams_.setLambda(2.0);
}
~WaterAirSpatialParameters()
{}
/*!
* \brief Update the spatial parameters with the flow solution
* after a timestep.
*
* \param globalSolution The global solution vector
*/
void update(const SolutionVector &globalSolution)
{
};
/*!
* \brief Apply the intrinsic permeability tensor to a pressure
* potential gradient.
*
* \param element The current finite element
* \param fvElemGeom The current finite volume geometry of the element
* \param scvIdx The index of the sub-control volume
*/
const Scalar intrinsicPermeability(const Element &element,
const FVElementGeometry &fvElemGeom,
int scvIdx) const
{
const GlobalPosition &pos = fvElemGeom.subContVol[scvIdx].global;
if (isFineMaterial_(pos))
return fineK_;
return coarseK_;
}
/*!
* \brief Define the porosity \f$[-]\f$ of the spatial parameters
*
* \param element The finite element
* \param fvElemGeom The finite volume geometry
* \param scvIdx The local index of the sub-control volume where
* the porosity needs to be defined
*/
double porosity(const Element &element,
const FVElementGeometry &fvElemGeom,
int scvIdx) const
{
const GlobalPosition &pos = fvElemGeom.subContVol[scvIdx].global;
if (isFineMaterial_(pos))
return finePorosity_;
else
return coarsePorosity_;
}
/*!
* \brief return the parameter object for the Brooks-Corey material law which depends on the position
*
* \param element The current finite element
* \param fvElemGeom The current finite volume geometry of the element
* \param scvIdx The index of the sub-control volume
*/
const MaterialLawParams& materialLawParams(const Element &element,
const FVElementGeometry &fvElemGeom,
int scvIdx) const
{
const GlobalPosition &pos = fvElemGeom.subContVol[scvIdx].global;
if (isFineMaterial_(pos))
return fineMaterialParams_;
else
return coarseMaterialParams_;
}
/*!
* \brief Returns the heat capacity \f$[J/m^3 K]\f$ of the rock matrix.
*
* This is only required for non-isothermal models.
*
* \param element The finite element
* \param fvElemGeom The finite volume geometry
* \param scvIdx The local index of the sub-control volume where
* the heat capacity needs to be defined
*/
double heatCapacity(const Element &element,
const FVElementGeometry &fvElemGeom,
int scvIdx) const
{
return
790 // specific heat capacity of granite [J / (kg K)]
* 2700 // density of granite [kg/m^3]
* (1 - porosity(element, fvElemGeom, scvIdx));
}
/*!
* \brief Calculate the heat flux \f$[W/m^2]\f$ through the
* rock matrix based on the temperature gradient \f$[K / m]\f$
*
* This is only required for non-isothermal models.
*
* \param heatFlux The resulting heat flux vector
* \param fluxDat The flux variables
* \param vDat The volume variables
* \param tempGrad The temperature gradient
* \param element The current finite element
* \param fvElemGeom The finite volume geometry of the current element
* \param scvfIdx The local index of the sub-control volume face where
* the matrix heat flux should be calculated
*/
void matrixHeatFlux(Vector &heatFlux,
const FluxVariables &fluxDat,
const ElementVolumeVariables &vDat,
const Vector &tempGrad,
const Element &element,
const FVElementGeometry &fvElemGeom,
int scvfIdx) const
{
static const Scalar lWater = 0.6;
static const Scalar lGranite = 2.8;
// arithmetic mean of the liquid saturation and the porosity
const int i = fvElemGeom.subContVolFace[scvfIdx].i;
const int j = fvElemGeom.subContVolFace[scvfIdx].j;
Scalar Sl = std::max<Scalar>(0.0, (vDat[i].saturation(lPhaseIdx) +
vDat[j].saturation(lPhaseIdx)) / 2);
Scalar poro = (porosity(element, fvElemGeom, i) +
porosity(element, fvElemGeom, j)) / 2;
Scalar lsat = pow(lGranite, (1-poro)) * pow(lWater, poro);
Scalar ldry = pow(lGranite, (1-poro));
// the heat conductivity of the matrix. in general this is a
// tensorial value, but we assume isotropic heat conductivity.
Scalar heatCond = ldry + sqrt(Sl) * (ldry - lsat);
// the matrix heat flux is the negative temperature gradient
// times the heat conductivity.
heatFlux = tempGrad;
heatFlux *= -heatCond;
}
private:
bool isFineMaterial_(const GlobalPosition &pos) const
{ return pos[dim-1] > layerBottom_; };
Scalar fineK_;
Scalar coarseK_;
Scalar layerBottom_;
Scalar finePorosity_;
Scalar coarsePorosity_;
MaterialLawParams fineMaterialParams_;
MaterialLawParams coarseMaterialParams_;
};
}
#endif