Commit ff1cd5bb authored by Alexander Kissinger's avatar Alexander Kissinger
Browse files

- Richards model can now be used with the nonisothermal model

- Included tests richardsniconduction and convection test for box and cc
- Fixed unit of thermalconductivitySolid() in all spatial parameters in the test folder from W/m^2 to W/(m K)

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@14054 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 27cfad4d
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* 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 Overwrites the volumeFlux() function from the ImplicitDarcyFluxVariables.
*
*/
#ifndef DUMUX_RICHARDS_FLUX_VARIABLES_HH
#define DUMUX_RICHARDS_FLUX_VARIABLES_HH
#include <dumux/common/math.hh>
#include <dumux/common/parameters.hh>
#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
#include "richardsproperties.hh"
namespace Dumux
{
/*!
* \ingroup RichardsModel
* \ingroup ImplicitFluxVariables
* \brief Overwrites the volumeFlux() function from the ImplicitDarcyFluxVariables.
*/
template <class TypeTag>
class RichardsFluxVariables : public ImplicitDarcyFluxVariables<TypeTag>
{
typedef ImplicitDarcyFluxVariables<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GridView::template Codim<0>::Entity Element;
enum { dim = GridView::dimension} ;
enum { dimWorld = GridView::dimensionworld} ;
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases)} ;
enum { nPhaseIdx = Indices::nPhaseIdx} ;
public:
/*
* \brief The constructor
*
* \param problem The problem
* \param element The finite element
* \param fvGeometry The finite-volume geometry
* \param faceIdx The local index of the SCV (sub-control-volume) face
* \param elemVolVars The volume variables of the current element
* \param onBoundary A boolean variable to specify whether the flux variables
* are calculated for interior SCV faces or boundary faces, default=false
*/
RichardsFluxVariables(const Problem &problem,
const Element &element,
const FVElementGeometry &fvGeometry,
const int fIdx,
const ElementVolumeVariables &elemVolVars,
const bool onBoundary = false)
: ParentType(problem, element, fvGeometry, fIdx, elemVolVars, onBoundary)
{ }
public:
/*!
* \brief Return the volumetric flux over a face of a given phase.
*
* This is the calculated velocity multiplied by the unit normal
* and the area of the face.
* face().normal
* has already the magnitude of the area.
* For the Richards model the velocity of the non-wetting phase
* is set to zero.
*
* \param phaseIdx index of the phase
*/
Scalar volumeFlux(const unsigned int phaseIdx) const
{
if(phaseIdx == nPhaseIdx)
return 0.;
else
return ParentType::volumeFlux(phaseIdx);
}
};
} // end namespace
#endif
......@@ -36,6 +36,7 @@ namespace Dumux
template<class TypeTag>
class RichardsLocalResidual : public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
{
typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
......@@ -113,12 +114,27 @@ public:
this->curVolVars_(),
onBoundary);
flux = 0;
asImp_()->computeAdvectiveFlux(flux, fluxVars);
asImp_()->computeDiffusiveFlux(flux, fluxVars);
}
/*!
* \brief Evaluates the advective mass flux of all components over
* a face of a sub-control volume.
*
* \param flux The advective flux over the sub-control-volume face for each component
* \param fluxVars The flux variables at the current SCV
*/
void computeAdvectiveFlux(PrimaryVariables &flux, const FluxVariables &fluxVars) const
{
// data attached to upstream and the downstream vertices
// of the current phase
const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx(wPhaseIdx));
const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx(wPhaseIdx));
flux[contiEqIdx] =
flux[contiEqIdx] +=
fluxVars.volumeFlux(wPhaseIdx)
*
(( massUpwindWeight_)*up.density(wPhaseIdx)
......@@ -126,6 +142,24 @@ public:
(1 - massUpwindWeight_)*dn.density(wPhaseIdx));
}
/*!
* \brief Adds the diffusive flux to the flux vector over
* the face of a sub-control volume.
*
* \param flux The diffusive flux over the sub-control-volume face for each phase
* \param fluxVars The flux variables at the current SCV
*
* This function doesn't do anything but may be used by the
* non-isothermal three-phase models to calculate diffusive heat
* fluxes
*/
void computeDiffusiveFlux(PrimaryVariables &flux, const FluxVariables &fluxVars) const
{
// diffusive fluxes
flux += 0.0;
}
/*!
* \brief Calculate the source term of the equation
*
......@@ -143,8 +177,21 @@ public:
this->curVolVars_());
}
protected:
Implementation *asImp_()
{
return static_cast<Implementation *> (this);
}
const Implementation *asImp_() const
{
return static_cast<const Implementation *> (this);
}
private:
Scalar massUpwindWeight_;
};
}
......
......@@ -30,6 +30,7 @@
#include <dumux/implicit/box/boxproperties.hh>
#include <dumux/implicit/cellcentered/ccproperties.hh>
#include <dumux/implicit/nonisothermal/niproperties.hh>
namespace Dumux
{
......@@ -43,11 +44,16 @@ namespace Properties {
// Type tags
//////////////////////////////////////////////////////////////////
//! The type tags for the implicit Richards problems
//! The type tags for the implicit isothermal one-phase two-component problems
NEW_TYPE_TAG(Richards);
NEW_TYPE_TAG(BoxRichards, INHERITS_FROM(BoxModel, Richards));
NEW_TYPE_TAG(CCRichards, INHERITS_FROM(CCModel, Richards));
//! The type tags for the corresponding non-isothermal problems
NEW_TYPE_TAG(RichardsNI, INHERITS_FROM(Richards, NonIsothermal));
NEW_TYPE_TAG(BoxRichardsNI, INHERITS_FROM(BoxModel, RichardsNI));
NEW_TYPE_TAG(CCRichardsNI, INHERITS_FROM(CCModel, RichardsNI));
//////////////////////////////////////////////////////////////////
// Property tags
//////////////////////////////////////////////////////////////////
......
......@@ -28,6 +28,7 @@
#ifndef DUMUX_RICHARDS_PROPERTY_DEFAULTS_HH
#define DUMUX_RICHARDS_PROPERTY_DEFAULTS_HH
#include "richardsfluxvariables.hh"
#include "richardsmodel.hh"
#include "richardsproblem.hh"
#include "richardsindices.hh"
......@@ -36,10 +37,12 @@
#include "richardsnewtoncontroller.hh"
#include "richardslocalresidual.hh"
#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
#include <dumux/implicit/nonisothermal/nipropertydefaults.hh>
#include <dumux/material/components/nullcomponent.hh>
#include <dumux/material/fluidmatrixinteractions/2p/thermalconductivitysomerton.hh>
#include <dumux/material/fluidsystems/2pimmisciblefluidsystem.hh>
#include <dumux/material/spatialparams/implicitspatialparams.hh>
#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
namespace Dumux
{
......@@ -50,23 +53,23 @@ namespace Properties {
// Properties values
//////////////////////////////////////////////////////////////////
//! Number of equations required by the model
SET_INT_PROP(Richards, NumEq, 1);
SET_INT_PROP(Richards, NumEq, GET_PROP_VALUE(TypeTag, IsothermalNumEq));
//! Number of fluid phases considered
SET_INT_PROP(Richards, NumPhases, 2);
//! The local residual operator
SET_TYPE_PROP(Richards,
LocalResidual,
RichardsLocalResidual<TypeTag>);
typename GET_PROP_TYPE(TypeTag, IsothermalLocalResidual));
//! The global model used
SET_TYPE_PROP(Richards, Model, RichardsModel<TypeTag>);
SET_TYPE_PROP(Richards, Model, typename GET_PROP_TYPE(TypeTag, IsothermalModel));
//! The class for the volume averaged quantities
SET_TYPE_PROP(Richards, VolumeVariables, RichardsVolumeVariables<TypeTag>);
SET_TYPE_PROP(Richards, VolumeVariables, typename GET_PROP_TYPE(TypeTag, IsothermalVolumeVariables));
//! The class for the quantities required for the flux calculation
SET_TYPE_PROP(Richards, FluxVariables, ImplicitDarcyFluxVariables<TypeTag>);
SET_TYPE_PROP(Richards, FluxVariables, typename GET_PROP_TYPE(TypeTag, IsothermalFluxVariables));
//! The class of the newton controller
SET_TYPE_PROP(Richards, NewtonController, RichardsNewtonController<TypeTag>);
......@@ -78,7 +81,7 @@ SET_SCALAR_PROP(Richards, ImplicitMassUpwindWeight, 1.0);
SET_SCALAR_PROP(Richards, ImplicitMobilityUpwindWeight, 1.0);
//! The class with all index definitions for the model
SET_TYPE_PROP(Richards, Indices, RichardsIndices<TypeTag>);
SET_TYPE_PROP(Richards, Indices, typename GET_PROP_TYPE(TypeTag, IsothermalIndices));
//! The spatial parameters to be employed.
//! Use ImplicitSpatialParams by default.
......@@ -164,6 +167,41 @@ SET_BOOL_PROP(Richards, ProblemEnableGravity, true);
// (Nield, Bejan, Convection in porous media, 2006, p. 10)
SET_SCALAR_PROP(BoxModel, SpatialParamsForchCoeff, 0.55);
//! Somerton is used as default model to compute the effective thermal heat conductivity
SET_PROP(NonIsothermal, ThermalConductivityModel)
{
private:
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
public:
typedef ThermalConductivitySomerton<Scalar, Indices> type;
};
//! temperature is already written by the isothermal model
SET_BOOL_PROP(RichardsNI, NiOutputLevel, 0);
//////////////////////////////////////////////////////////////////
// Property values for isothermal model required for the general non-isothermal model
//////////////////////////////////////////////////////////////////
// set isothermal Model
SET_TYPE_PROP(Richards, IsothermalModel, RichardsModel<TypeTag>);
// set isothermal FluxVariables
SET_TYPE_PROP(Richards, IsothermalFluxVariables, RichardsFluxVariables<TypeTag>);
//set isothermal VolumeVariables
SET_TYPE_PROP(Richards, IsothermalVolumeVariables, RichardsVolumeVariables<TypeTag>);
//set isothermal LocalResidual
SET_TYPE_PROP(Richards, IsothermalLocalResidual, RichardsLocalResidual<TypeTag>);
//set isothermal Indices
SET_TYPE_PROP(Richards, IsothermalIndices, RichardsIndices<TypeTag>);
//set isothermal NumEq
SET_INT_PROP(Richards, IsothermalNumEq, 1);
// \}
}
......
......@@ -140,6 +140,11 @@ public:
fluidState.setViscosity(wPhaseIdx, FluidSystem::viscosity(fluidState, paramCache, wPhaseIdx));
fluidState.setViscosity(nPhaseIdx, 1e-10);
// compute and set the enthalpy
fluidState.setEnthalpy(wPhaseIdx, Implementation::enthalpy_(fluidState, paramCache, wPhaseIdx));
fluidState.setEnthalpy(nPhaseIdx, Implementation::enthalpy_(fluidState, paramCache, nPhaseIdx));
}
/*!
......@@ -251,6 +256,14 @@ protected:
{
return problem.temperatureAtPos(fvGeometry.subContVol[scvIdx].global);
}
template<class ParameterCache>
static Scalar enthalpy_(const FluidState& fluidState,
const ParameterCache& paramCache,
int phaseIdx)
{
return 0;
}
/*!
* \brief Called by update() to compute the energy related quantities
......
......@@ -144,7 +144,7 @@ public:
}
/*!
* \brief Returns the thermal conductivity \f$[W/m^2]\f$ of the porous material.
* \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the porous material.
*
* \param element The finite element
* \param fvGeometry The finite volume geometry
......
......@@ -144,7 +144,7 @@ public:
}
/*!
* \brief Returns the thermal conductivity \f$[W/m^2]\f$ of the porous material.
* \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the porous material.
*
* \param element The finite element
* \param fvGeometry The finite volume geometry
......
......@@ -146,7 +146,7 @@ public:
}
/*!
* \brief Returns the thermal conductivity \f$[W/m^2]\f$ of the porous material.
* \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the porous material.
*
* \param element The finite element
* \param fvGeometry The finite volume geometry
......
......@@ -200,7 +200,7 @@ public:
}
/*!
* \brief Returns the thermal conductivity \f$[W/m^2]\f$ of the solid
* \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the solid
*
* This is only required for non-isothermal models.
*
......
......@@ -195,7 +195,7 @@ public:
}
/*!
* \brief Returns the thermal conductivity \f$[W/m^2]\f$ of the porous material.
* \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the porous material.
*
* \param element The finite element
* \param fvGeometry The finite volume geometry
......
......@@ -187,7 +187,7 @@ public:
}
/*!
* \brief Returns the thermal conductivity \f$[W/m^2]\f$ of the porous material.
* \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the porous material.
*
* \param element The finite element
* \param fvGeometry The finite volume geometry
......
......@@ -235,7 +235,7 @@ public:
/*!
* \brief Returns the thermal conductivity \f$[W/m^2]\f$ of the porous material.
* \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the porous material.
*
* \param element The finite element
* \param fvGeometry The finite volume geometry
......
......@@ -206,55 +206,11 @@ public:
* (1 - porosity(element, fvGeometry, 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 fluxVars The flux variables
* \param elemVolVars The volume variables
* \param tempGrad The temperature gradient
* \param element The current finite element
* \param fvGeometry 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(DimVector &heatFlux,
const FluxVariables &fluxVars,
const ElementVolumeVariables &elemVolVars,
const DimVector &tempGrad,
const Element &element,
const FVElementGeometry &fvGeometry,
int scvfIdx) const
{
static const Scalar lDry = 0.35;
static const Scalar lsw1 = 1.8;
static const Scalar lsn1 = 0.65;
// arithmetic mean of the liquid saturation and the porosity
const int i = fluxVars.face().i;
const int j = fluxVars.face().j;
Scalar sw = std::max(0.0, (elemVolVars[i].saturation(wPhaseIdx) +
elemVolVars[j].saturation(wPhaseIdx)) / 2);
Scalar sn = std::max(0.0, (elemVolVars[i].saturation(nPhaseIdx) +
elemVolVars[j].saturation(nPhaseIdx)) / 2);
// the heat conductivity of the matrix. in general this is a
// tensorial value, but we assume isotropic heat conductivity.
Scalar heatCond = lDry + sqrt(sw) * (lsw1-lDry) + sqrt(sn) * (lsn1-lDry);
// the matrix heat flux is the negative temperature gradient
// times the heat conductivity.
heatFlux = tempGrad;
heatFlux *= -heatCond;
}
const MaterialLawParams& materialLawParams() const
{
return MaterialParams_;
}
private:
bool isFineMaterial_(const GlobalPosition &globalPos) const
{ return
......
......@@ -226,52 +226,7 @@ public:
}
/*!
* \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 fluxVars The flux variables
* \param elemVolVars The volume variables
* \param tempGrad The temperature gradient
* \param element The current finite element
* \param fvGeometry The finite volume geometry of the current element
* \param fIdx The local index of the sub-control volume face where
* the matrix heat flux should be calculated
*/
void matrixHeatFlux(DimVector &heatFlux,
const FluxVariables &fluxVars,
const ElementVolumeVariables &elemVolVars,
const DimVector &tempGrad,
const Element &element,
const FVElementGeometry &fvGeometry,
const int fIdx) const
{
static const Scalar lDry = 0.35;
static const Scalar lsw1 = 1.8;
static const Scalar lsn1 = 0.65;
// arithmetic mean of the liquid saturation and the porosity
const int i = fluxVars.face().i;
const int j = fluxVars.face().j;
Scalar sw = std::max(0.0, (elemVolVars[i].saturation(wPhaseIdx) +
elemVolVars[j].saturation(wPhaseIdx)) / 2);
Scalar sn = std::max(0.0, (elemVolVars[i].saturation(nPhaseIdx) +
elemVolVars[j].saturation(nPhaseIdx)) / 2);
// the heat conductivity of the matrix. in general this is a
// tensorial value, but we assume isotropic heat conductivity.
Scalar heatCond = lDry + sqrt(sw) * (lsw1-lDry) + sqrt(sn) * (lsn1-lDry);
// the matrix heat flux is the negative temperature gradient
// times the heat conductivity.
heatFlux = tempGrad;
heatFlux *= -heatCond;
}
/*!
* \brief Returns the thermal conductivity \f$[W/m^2]\f$ of the porous material.
* \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the porous material.
*
* \param element The finite element
* \param fvGeometry The finite volume geometry
......
......@@ -231,7 +231,7 @@ public:
}
/*!
* \brief Returns the thermal conductivity \f$[W/m^2]\f$ of the solid
* \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the solid
*
* This is only required for non-isothermal models.
*
......
......@@ -6,6 +6,7 @@ if(MPI_CXX_FOUND)
${CMAKE_CURRENT_BINARY_DIR}/s0002-p0000-richardslensbox-00008.vtu
"${MPIEXEC} -np 2 ${CMAKE_CURRENT_BINARY_DIR}/test_boxrichards")
else()
# isothermal tests
add_dumux_test(test_boxrichards test_boxrichards test_boxrichards.cc
${CMAKE_SOURCE_DIR}/bin/runTest.sh
${CMAKE_SOURCE_DIR}/bin/fuzzycomparevtu.py
......@@ -20,3 +21,32 @@ add_dumux_test(test_ccrichards test_ccrichards test_ccrichards.cc
${CMAKE_SOURCE_DIR}/test/references/richardslenscc-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/richardslenscc-00008.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_ccrichards)
# non-isothermal tests
add_dumux_test(test_boxrichardsniconvection test_boxrichardsniconvection test_boxrichardsniconvection.cc
${CMAKE_SOURCE_DIR}/bin/runTest.sh
${CMAKE_SOURCE_DIR}/bin/fuzzycomparevtu.py
${CMAKE_SOURCE_DIR}/test/references/richardsniconvectionbox-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_boxrichardsniconvection-00011.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_boxrichardsniconvection)
add_dumux_test(test_ccrichardsniconvection test_ccrichardsniconvection test_ccrichardsniconvection.cc
${CMAKE_SOURCE_DIR}/bin/runTest.sh
${CMAKE_SOURCE_DIR}/bin/fuzzycomparevtu.py
${CMAKE_SOURCE_DIR}/test/references/richardsniconvectioncc-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_ccrichardsniconvection-00011.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_ccrichardsniconvection)
add_dumux_test(test_boxrichardsniconduction test_boxrichardsniconduction test_boxrichardsniconduction.cc
${CMAKE_SOURCE_DIR}/bin/runTest.sh
${CMAKE_SOURCE_DIR}/bin/fuzzycomparevtu.py
${CMAKE_SOURCE_DIR}/test/references/richardsniconductionbox-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_boxrichardsniconduction-00007.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_boxrichardsniconduction)
add_dumux_test(test_ccrichardsniconduction test_ccrichardsniconduction test_ccrichardsniconduction.cc
${CMAKE_SOURCE_DIR}/bin/runTest.sh
${CMAKE_SOURCE_DIR}/bin/fuzzycomparevtu.py
${CMAKE_SOURCE_DIR}/test/references/richardsniconductioncc-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_ccrichardsniconduction-00007.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_ccrichardsniconduction)
\ No newline at end of file
check_PROGRAMS = test_boxrichards test_ccrichards
check_PROGRAMS = \
test_boxrichards \
test_ccrichards \
test_boxrichardsniconvection \
test_ccrichardsniconvection \
test_boxrichardsniconduction \
test_ccrichardsniconduction
noinst_HEADERS := $(wildcard *.hh)
EXTRA_DIST:=$(wildcard *.input) $(wildcard grids/*.dgf) CMakeLists.txt
test_boxrichards_SOURCES = test_boxrichards.cc
test_ccrichards_SOURCES = test_ccrichards.cc
test_boxrichardsniconvection_SOURCES = test_boxrichardsniconvection.cc
test_ccrichardsniconvection_SOURCES = test_ccrichardsniconvection.cc
test_boxrichardsniconduction_SOURCES = test_boxrichardsniconduction.cc
test_ccrichardsniconduction_SOURCES = test_ccrichardsniconduction.cc
include $(top_srcdir)/am/global-rules