Commit 1dc11abb authored by Timo Koch's avatar Timo Koch
Browse files

[solidenergy] New model for energy equation of solid

parent 6abe298f
......@@ -5,6 +5,8 @@ Differences Between DuMuX 3.1 and DuMuX 3.0
### Improvements and Enhancements
- Added new porous medium model for the energy balance of a porous solid (general heat equation)
### Immediate interface changes not allowing/requiring a deprecation period
### Deprecated classes/files, to be removed after 3.1:
......
......@@ -166,6 +166,12 @@
* \brief Richards multi-component flow
* \copydetails ./porousmediumflow/richardsnc/model.hh
*/
/*!
* \ingroup PorousmediumflowModels
* \defgroup SolidEnergyModel Solid energy
* \brief Energy equation for the solid (general heat equation)
* \copydetails ./porousmediumflow/solidenergy/model.hh
*/
/*!
* \ingroup PorousmediumflowModels
* \defgroup TracerModel Tracer
......
......@@ -20,6 +20,7 @@ add_subdirectory(nonisothermal)
add_subdirectory(richards)
add_subdirectory(richardsnc)
add_subdirectory(sequential)
add_subdirectory(solidenergy)
add_subdirectory(tracer)
install(FILES
......
install(FILES
localresidual.hh
model.hh
volumevariables.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/solidenergy)
// -*- 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
* \ingroup SolidEnergyModel
* \brief Element-wise calculation of the residual
*/
#ifndef DUMUX_SOLID_ENERGY_LOCAL_RESIDUAL_HH
#define DUMUX_SOLID_ENERGY_LOCAL_RESIDUAL_HH
#include <dumux/common/properties.hh>
namespace Dumux {
/*!
* \ingroup SolidEnergyModel
* \brief Element-wise calculation of the residual
*/
template<class TypeTag>
class SolidEnergyLocalResidual : public GetPropType<TypeTag, Properties::BaseLocalResidual>
{
using ParentType = GetPropType<TypeTag, Properties::BaseLocalResidual>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Problem = GetPropType<TypeTag, Properties::Problem>;
using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
using FluxVariables = GetPropType<TypeTag, Properties::FluxVariables>;
using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using Element = typename GridView::template Codim<0>::Entity;
public:
using ParentType::ParentType;
/*!
* \brief Evaluate the rate of change of all conservation quantites
* \param problem The problem
* \param scv The sub control volume
* \param volVars The current or previous volVars
* \note This function should not include the source and sink terms.
* \note The volVars can be different to allow computing
* the implicit/explicit euler time derivative here
*/
NumEqVector computeStorage(const Problem& problem,
const SubControlVolume& scv,
const VolumeVariables& volVars) const
{
NumEqVector storage;
storage[0] = volVars.temperatureSolid()
* volVars.solidHeatCapacity()
* volVars.solidDensity()
* (1.0 - volVars.porosity());
return storage;
}
/*!
* \brief Evaluate the energy flux over a face of a sub control volume
*
* \param problem The problem
* \param element The element
* \param fvGeometry The finite volume geometry context
* \param elemVolVars The volume variables for all flux stencil elements
* \param scvf The sub control volume face to compute the flux on
* \param elemFluxVarsCache The cache related to flux compuation
*/
NumEqVector computeFlux(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf,
const ElementFluxVariablesCache& elemFluxVarsCache) const
{
FluxVariables fluxVars;
fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
return fluxVars.heatConductionFlux();
}
};
} // end namespace Dumux
#endif
// -*- 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
* \ingroup SolidEnergyModel
* \brief The energy balance equation for a porous solid
*
* The energy balance is described by the following equation:
\f[
\frac{ \partial n c_p \varrho T}{\partial t}
- \text{div} \left\lbrace \lambda_\text{pm} \textbf{grad} T \right\rbrace = q,
\f]
* where \f$n\f$ is the volume fraction of the conducting material, \f$c_p\f$ its specific heat capacity,
* \f$\varrho\f$ its density, \f$T\f$ the temperature, and \f$\lambda\f$ the heat conductivity of the porous solid.
*/
#ifndef DUMUX_SOLID_ENERGY_MODEL_HH
#define DUMUX_SOLID_ENERGY_MODEL_HH
#include <dumux/common/properties.hh>
#include <dumux/porousmediumflow/properties.hh>
#include <dumux/porousmediumflow/nonisothermal/iofields.hh>
#include <dumux/porousmediumflow/solidenergy/volumevariables.hh>
#include <dumux/porousmediumflow/solidenergy/localresidual.hh>
namespace Dumux {
/*!
* \ingroup SolidEnergyModel
* \brief The indices
*/
struct SolidEnergyIndices
{
static constexpr int temperatureIdx = 0;
static constexpr int energyEqIdx = 0;
};
/*!
* \ingroup SolidEnergyModel
* \brief The energy balance equation for a porous solid
*/
struct SolidEnergyModelTraits
{
using Indices = SolidEnergyIndices;
static constexpr int numEq() { return 1; }
static constexpr int numFluidPhases() { return 0; }
static constexpr int numFluidComponents() { return 0; }
static constexpr int numSolidComponents() { return 1; }
static constexpr bool enableAdvection() { return false; }
static constexpr bool enableMolecularDiffusion() { return false; }
static constexpr bool enableEnergyBalance() { return true; }
};
/*!
* \ingroup SolidEnergyModel
* \brief The volume variable traits
*/
template<class PV, class SSY, class SST, class MT>
struct SolidEnergyVolumeVariablesTraits
{
using PrimaryVariables = PV;
using SolidSystem = SSY;
using SolidState = SST;
using ModelTraits = MT;
};
struct ThermalConductivitySolidEnergy
{
/*!
* \brief Relation for a simple effective thermal conductivity \f$\mathrm{[W/(m K)]}\f$
*
* \param volVars volume variables
* \param spatialParams spatial parameters
* \param element element (to be passed to spatialParams)
* \param fvGeometry fvGeometry (to be passed to spatialParams)
* \param scv scv (to be passed to spatialParams)
*
* \return effective thermal conductivity \f$\mathrm{[W/(m K)]}\f$
*/
template<class VolumeVariables, class SpatialParams, class Element, class FVGeometry>
static typename VolumeVariables::PrimaryVariables::value_type
effectiveThermalConductivity(const VolumeVariables& volVars,
const SpatialParams& spatialParams,
const Element& element,
const FVGeometry& fvGeometry,
const typename FVGeometry::SubControlVolume& scv)
{
return volVars.solidThermalConductivity()*(1.0-volVars.porosity());
}
};
// \{
namespace Properties {
//////////////////////////////////////////////////////////////////
// Type tags
//////////////////////////////////////////////////////////////////
//! The type tags for the fully implicit tracer model.
// Create new type tags
namespace TTag {
struct SolidEnergy { using InheritsFrom = std::tuple<PorousMediumFlow>; };
} // end namespace TTag
///////////////////////////////////////////////////////////////////////////
// Properties
///////////////////////////////////////////////////////////////////////////
//! set the model traits
template<class TypeTag>
struct ModelTraits<TypeTag, TTag::SolidEnergy> { using type = SolidEnergyModelTraits; };
//! set the local residual
template<class TypeTag>
struct LocalResidual<TypeTag, TTag::SolidEnergy> { using type = SolidEnergyLocalResidual<TypeTag>; };
//! Set the vtk output fields specific to this model
template<class TypeTag>
struct IOFields<TypeTag, TTag::SolidEnergy> { using type = EnergyIOFields<>; };
//! Set the volume variables property
template<class TypeTag>
struct VolumeVariables<TypeTag, TTag::SolidEnergy>
{
private:
using PV = GetPropType<TypeTag, Properties::PrimaryVariables>;
using SSY = GetPropType<TypeTag, Properties::SolidSystem>;
using SST = GetPropType<TypeTag, Properties::SolidState>;
using MT = GetPropType<TypeTag, Properties::ModelTraits>;
using Traits = SolidEnergyVolumeVariablesTraits<PV, SSY, SST, MT>;
public:
using type = SolidEnergyVolumeVariables<Traits>;
};
//! Use the average for effective conductivities
template<class TypeTag>
struct ThermalConductivityModel<TypeTag, TTag::SolidEnergy>
{ using type = ThermalConductivitySolidEnergy; };
} // end namespace Properties
// \}
} // end namespace Dumux
#endif
// -*- 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
* \ingroup SolidEnergyModel
* \brief Class for computation of all volume averaged quantities
*/
#ifndef DUMUX_SOLID_ENERGY_VOLUME_VARIABLES_HH
#define DUMUX_SOLID_ENERGY_VOLUME_VARIABLES_HH
#include <type_traits>
#include <dumux/material/solidstates/updatesolidvolumefractions.hh>
#include <dumux/porousmediumflow/volumevariables.hh>
#include <dumux/porousmediumflow/nonisothermal/volumevariables.hh>
namespace Dumux {
/*!
* \ingroup SolidEnergyModel
* \brief Class for computation of all volume averaged quantities
*/
template<class Traits>
class SolidEnergyVolumeVariables : public PorousMediumFlowVolumeVariables<Traits>
{
using ParentType = PorousMediumFlowVolumeVariables<Traits>;
using Scalar = typename Traits::PrimaryVariables::value_type;
static constexpr int temperatureIdx = Traits::ModelTraits::Indices::temperatureIdx;
public:
//! export the type used for the solid state
using SolidState = typename Traits::SolidState;
//! export the type used for the solid system
using SolidSystem = typename Traits::SolidSystem;
/*!
* \brief Update all quantities for a given control volume
*
* \param elemSol A vector containing all primary variables connected to the element
* \param problem The object specifying the problem which ought to
* be simulated
* \param element An element which contains part of the control volume
* \param scv The sub-control volume
*/
template<class ElemSol, class Problem, class Element, class Scv>
void update(const ElemSol& elemSol,
const Problem& problem,
const Element& element,
const Scv& scv)
{
ParentType::update(elemSol, problem, element, scv);
updateTemperature(elemSol, problem, element, scv, solidState_);
updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, 0);
updateSolidEnergyParams(elemSol, problem, element, scv, solidState_);
}
//! Fill temperature in the solid state
template<class ElemSol, class Problem, class Element, class Scv>
void updateTemperature(const ElemSol& elemSol,
const Problem& problem,
const Element& element,
const Scv& scv,
SolidState& solidState)
{
const Scalar T = elemSol[scv.localDofIndex()][temperatureIdx];
solidState.setTemperature(T);
}
//! Fill solid matrix parameters in the solid state
template<class ElemSol, class Problem, class Element, class Scv>
void updateSolidEnergyParams(const ElemSol &elemSol,
const Problem& problem,
const Element &element,
const Scv &scv,
SolidState & solidState)
{
Scalar cs = solidHeatCapacity_(elemSol, problem, element, scv, solidState);
solidState.setHeatCapacity(cs);
Scalar rhos = solidDensity_(elemSol, problem, element, scv, solidState);
solidState.setDensity(rhos);
Scalar lambdas = solidThermalConductivity_(elemSol, problem, element, scv, solidState);
solidState.setThermalConductivity(lambdas);
}
/*!
* \brief Returns the temperature in the sub-control volume.
*/
Scalar temperatureSolid() const
{ return solidState_.temperature(); }
/*!
* \brief Returns the temperature in the sub-control volume.
*/
Scalar temperature() const
{ return solidState_.temperature(); }
/*!
* \brief Returns the total heat capacity \f$\mathrm{[J/(kg K)]}\f$ of the rock matrix in
* the sub-control volume.
*/
Scalar solidHeatCapacity() const
{ return solidState_.heatCapacity(); }
/*!
* \brief Returns the mass density \f$\mathrm{[kg/m^3]}\f$ of the rock matrix in
* the sub-control volume.
*/
Scalar solidDensity() const
{ return solidState_.density(); }
/*!
* \brief Returns the thermal conductivity \f$\mathrm{[W/(m*K)]}\f$ of the solid phase in the sub-control volume.
*/
Scalar solidThermalConductivity() const
{ return solidState_.thermalConductivity(); }
/*!
* \brief Return the average porosity \f$\mathrm{[-]}\f$ within the control volume.
*/
Scalar porosity() const
{ return solidState_.porosity(); }
private:
//! the solid state
SolidState solidState_;
/*!
* It has to be decided if the full solid system / solid state interface is used (general option, but more complicated),
* or the simple nonisothermal spatial params interface (simpler but less general).
* In the simple nonisothermal spatial params interface the functions solidHeatCapacity, solidDensity, and solidThermalConductivity
* in the spatial params overwrite the parameters given in the solid system. This only makes sense in combination
* with the simplest solid system InertSolidPhase, and can be used to quickly change parameters in certain domain regions.
* For setups with more general solids with several components these functions should not exist. Instead, the solid system
* determines the values for solidHeatCapacity, solidDensity, and solidThermalConductivity depending on the given composition.
*/
/*!
* \name Access functions for the solidsystem / solidstate interface
*/
// \{
/*!
* \brief get the solid heat capacity in an scv
* \param elemSol the element solution vector
* \param problem the problem to solve
* \param element the element (codim-0-entity) the scv belongs to
* \param scv the sub control volume
* \param solidState the solid state
* \note this gets selected if the user uses the solidsystem / solidstate interface
*/
template<class ElemSol, class Problem, class Element, class Scv,
std::enable_if_t<!Detail::hasSolidHeatCapacity<typename Problem::SpatialParams, Element, Scv, ElemSol, SolidState>(), int> = 0>
Scalar solidHeatCapacity_(const ElemSol& elemSol,
const Problem& problem,
const Element& element,
const Scv& scv,
const SolidState& solidState)
{
return SolidSystem::heatCapacity(solidState);
}
/*!
* \brief get the solid density in an scv
* \param elemSol the element solution vector
* \param problem the problem to solve
* \param element the element (codim-0-entity) the scv belongs to
* \param scv the sub control volume
* \param solidState the solid state
* \note this gets selected if the user uses the solidsystem / solidstate interface
*/
template<class ElemSol, class Problem, class Element, class Scv,
std::enable_if_t<!Detail::hasSolidDensity<typename Problem::SpatialParams, Element, Scv, ElemSol, SolidState>(), int> = 0>
Scalar solidDensity_(const ElemSol& elemSol,
const Problem& problem,
const Element& element,
const Scv& scv,
const SolidState& solidState)
{
return SolidSystem::density(solidState);
}
/*!
* \brief get the solid's thermal conductivity in an scv
* \param elemSol the element solution vector
* \param problem the problem to solve
* \param element the element (codim-0-entity) the scv belongs to
* \param scv the sub control volume
* \param solidState the solid state
* \note this gets selected if the user uses the solidsystem / solidstate interface
*/
template<class ElemSol, class Problem, class Element, class Scv,
std::enable_if_t<!Detail::hasSolidThermalConductivity<typename Problem::SpatialParams, Element, Scv, ElemSol, SolidState>(), int> = 0>
Scalar solidThermalConductivity_(const ElemSol& elemSol,
const Problem& problem,
const Element& element,
const Scv& scv,
const SolidState& solidState)
{
return SolidSystem::thermalConductivity(solidState);
}
// \}
/*!
* \name Access functions for the simple nonisothermal spatial params interface in
* combination with an InertSolidPhase as solid system
*/
// \{
/*!
* \brief get the solid heat capacity in an scv
* \param elemSol the element solution vector
* \param problem the problem to solve
* \param element the element (codim-0-entity) the scv belongs to
* \param scv the sub control volume
* \param solidState the solid state
* \note this gets selected if the user uses the simple spatial params interface in
* combination with an InertSolidPhase as solid system
*/
template<class ElemSol, class Problem, class Element, class Scv,
std::enable_if_t<Detail::hasSolidHeatCapacity<typename Problem::SpatialParams, Element, Scv, ElemSol, SolidState>(), int> = 0>
Scalar solidHeatCapacity_(const ElemSol& elemSol,
const Problem& problem,
const Element& element,
const Scv& scv,
const SolidState& solidState)
{
static_assert(Detail::isInertSolidPhase<SolidSystem>::value,
"solidHeatCapacity can only be overwritten in the spatial params when the solid system is a simple InertSolidPhase\n"
"If you select a proper solid system, the solid heat capacity will be computed as stated in the solid system!");
return problem.spatialParams().solidHeatCapacity(element, scv, elemSol, solidState);
}
/*!
* \brief get the solid density in an scv
* \param elemSol the element solution vector
* \param problem the problem to solve
* \param element the element (codim-0-entity) the scv belongs to
* \param scv the sub control volume
* \param solidState the solid state
* \note this gets selected if the user uses the simple spatial params interface in
* combination with an InertSolidPhase as solid system
*/
template<class ElemSol, class Problem, class Element, class Scv,
std::enable_if_t<Detail::hasSolidDensity<typename Problem::SpatialParams, Element, Scv, ElemSol, SolidState>(), int> = 0>
Scalar solidDensity_(const ElemSol& elemSol,
const Problem& problem,
const Element& element,
const Scv& scv,
const SolidState& solidState)
{
static_assert(Detail::isInertSolidPhase<SolidSystem>::value,
"solidDensity can only be overwritten in the spatial params when the solid system is a simple InertSolidPhase\n"
"If you select a proper solid system, the solid density will be computed as stated in the solid system!");
return problem.spatialParams().solidDensity(element, scv, elemSol, solidState);
}
/*!
* \brief get the solid's heat capacity in an scv
* \param elemSol the element solution vector
* \param problem the problem to solve
* \param element the element (codim-0-entity) the scv belongs to
* \param scv the sub control volume
* \param solidState the solid state
* \note this gets selected if the user uses the simple spatial params interface in
* combination with an InertSolidPhase as solid system
*/
template<class ElemSol, class Problem, class Element, class Scv,
std::enable_if_t<Detail::hasSolidThermalConductivity<typename Problem::SpatialParams, Element, Scv, ElemSol, SolidState>(), int> = 0>
Scalar solidThermalConductivity_(const ElemSol& elemSol,
const Problem& problem,