From c1e3fe548943fd1d48ef052b99325f0ce70ac385 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= <dennis.glaeser@iws.uni-stuttgart.de> Date: Mon, 7 May 2018 21:44:29 +0200 Subject: [PATCH] reintroduce geomechanics framework for box --- dumux/CMakeLists.txt | 1 + dumux/common/properties.hh | 5 + dumux/discretization/box/hookeslaw.hh | 126 +++++++++++++++ dumux/discretization/hookeslaw.hh | 45 ++++++ dumux/geomechanics/CMakeLists.txt | 9 ++ dumux/geomechanics/elastic/CMakeLists.txt | 6 + dumux/geomechanics/elastic/indices.hh | 58 +++++++ dumux/geomechanics/elastic/localresidual.hh | 151 ++++++++++++++++++ dumux/geomechanics/elastic/model.hh | 113 +++++++++++++ dumux/geomechanics/elastic/volumevariables.hh | 121 ++++++++++++++ dumux/geomechanics/fvproblem.hh | 49 ++++++ dumux/geomechanics/lameparams.hh | 57 +++++++ dumux/geomechanics/properties.hh | 78 +++++++++ dumux/geomechanics/stressvariablescache.hh | 66 ++++++++ dumux/geomechanics/velocityoutput.hh | 60 +++++++ dumux/material/spatialparams/fvelastic.hh | 143 +++++++++++++++++ 16 files changed, 1088 insertions(+) create mode 100644 dumux/discretization/box/hookeslaw.hh create mode 100644 dumux/discretization/hookeslaw.hh create mode 100644 dumux/geomechanics/CMakeLists.txt create mode 100644 dumux/geomechanics/elastic/CMakeLists.txt create mode 100644 dumux/geomechanics/elastic/indices.hh create mode 100644 dumux/geomechanics/elastic/localresidual.hh create mode 100644 dumux/geomechanics/elastic/model.hh create mode 100644 dumux/geomechanics/elastic/volumevariables.hh create mode 100644 dumux/geomechanics/fvproblem.hh create mode 100644 dumux/geomechanics/lameparams.hh create mode 100644 dumux/geomechanics/properties.hh create mode 100644 dumux/geomechanics/stressvariablescache.hh create mode 100644 dumux/geomechanics/velocityoutput.hh create mode 100644 dumux/material/spatialparams/fvelastic.hh diff --git a/dumux/CMakeLists.txt b/dumux/CMakeLists.txt index 696feee855..61918c6fcf 100644 --- a/dumux/CMakeLists.txt +++ b/dumux/CMakeLists.txt @@ -2,6 +2,7 @@ add_subdirectory("adaptive") add_subdirectory("assembly") add_subdirectory("common") add_subdirectory("discretization") +add_subdirectory("geomechanics") add_subdirectory("freeflow") add_subdirectory("io") add_subdirectory("linear") diff --git a/dumux/common/properties.hh b/dumux/common/properties.hh index 78a0ac7c6a..08a752fc93 100644 --- a/dumux/common/properties.hh +++ b/dumux/common/properties.hh @@ -138,6 +138,11 @@ NEW_PROP_TAG(EnableWaterDiffusionInAir); //!< Property for turning Richards into ////////////////////////////////////////////////////////////// NEW_PROP_TAG(OnlyGasPhaseCanDisappear); //!< reduces the phasestates to threePhases and wnPhaseOnly +///////////////////////////////////////////////////////////// +// Properties used by geomechanical models: +///////////////////////////////////////////////////////////// +NEW_PROP_TAG(StressType); //!< The type used for the evaluation of stress tensors and forces + ///////////////////////////////////////////////////////////// // Properties used by the staggered-grid discretization method ///////////////////////////////////////////////////////////// diff --git a/dumux/discretization/box/hookeslaw.hh b/dumux/discretization/box/hookeslaw.hh new file mode 100644 index 0000000000..0749d4b05e --- /dev/null +++ b/dumux/discretization/box/hookeslaw.hh @@ -0,0 +1,126 @@ +// -*- 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 Specialization of Hooke's law for the box scheme. This computes + * the stress tensor and surface forces resulting from mechanical deformation. + */ +#ifndef DUMUX_DISCRETIZATION_BOX_HOOKES_LAW_HH +#define DUMUX_DISCRETIZATION_BOX_HOOKES_LAW_HH + +#include <dune/common/fmatrix.hh> + +#include <dumux/common/math.hh> +#include <dumux/common/parameters.hh> +#include <dumux/discretization/methods.hh> + +namespace Dumux { + +/*! + * \ingroup BoxDiscretization + * \brief Hooke's law for box scheme + * \tparam Scalar the scalar type for scalar physical quantities + * \tparam FVGridGeometry the grid geometry + */ +template<class Scalar, class FVGridGeometry> +class HookesLaw<Scalar, FVGridGeometry, DiscretizationMethod::box> +{ + using FVElementGeometry = typename FVGridGeometry::LocalView; + using SubControlVolume = typename FVElementGeometry::SubControlVolume; + using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; + + using GridView = typename FVGridGeometry::GridView; + using Element = typename GridView::template Codim<0>::Entity; + using CoordScalar = typename GridView::ctype; + + static constexpr int dim = GridView::dimension; + static constexpr int dimWorld = GridView::dimensionworld; + static_assert(dim == dimWorld, "Hookes Law not implemented for network/surface grids"); + +public: + //! export the type used for the stress tensor + using StressTensor = Dune::FieldMatrix<Scalar, dim, dimWorld>; + //! export the type used for force vectors + using ForceVector = typename StressTensor::row_type; + + //! computes the force acting on a sub-control volume face + template<class Problem, class ElementVolumeVariables, class ElementFluxVarsCache> + static ForceVector force(const Problem& problem, + const Element& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const SubControlVolumeFace& scvf, + const ElementFluxVarsCache& elemFluxVarCache) + { + const auto sigma = stressTensor(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarCache); + + ForceVector scvfForce(0.0); + sigma.mv(scvf.unitOuterNormal(), scvfForce); + scvfForce *= scvf.area(); + + return scvfForce; + } + + //! assembles the stress tensor at the integration point of an scvf + template<class Problem, class ElementVolumeVariables, class ElementFluxVarsCache> + static StressTensor stressTensor(const Problem& problem, + const Element& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const SubControlVolumeFace& scvf, + const ElementFluxVarsCache& elemFluxVarCache) + { + const auto& fluxVarCache = elemFluxVarCache[scvf]; + const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); + const auto& insideVolVars = elemVolVars[insideScv]; + + // TODO What to do with extrusion factor?? + // TODO this only works for element-wise constant lame params + // How could this be done for heterogeneous params within the element? + const auto& lameParams = insideVolVars.lameParams(); + + // evaluate displacement gradient + StressTensor gradU(0.0); + for (int dir = 0; dir < dim; ++dir) + for (const auto& scv : scvs(fvGeometry)) + gradU[dir].axpy(elemVolVars[scv.indexInElement()].displacement(dir), fluxVarCache.gradN(scv.indexInElement())); + + // evaluate strain tensor + StressTensor epsilon; + for (int i = 0; i < dim; ++i) + for (int j = 0; j < dimWorld; ++j) + epsilon[i][j] = 0.5*(gradU[i][j] + gradU[j][i]); + + // calculate sigma + StressTensor sigma(0.0); + const auto traceEpsilon = trace(epsilon); + for (int i = 0; i < dim; ++i) + { + sigma[i][i] = lameParams.lambda()*traceEpsilon; + for (int j = 0; j < dimWorld; ++j) + sigma[i][j] += 2.0*lameParams.mu()*epsilon[i][j]; + } + + return sigma; + } +}; + +} // end namespace Dumux + +#endif diff --git a/dumux/discretization/hookeslaw.hh b/dumux/discretization/hookeslaw.hh new file mode 100644 index 0000000000..f3abc81ec7 --- /dev/null +++ b/dumux/discretization/hookeslaw.hh @@ -0,0 +1,45 @@ +// -*- 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 Discretization + * \brief Hooke's law specialized for different discretization schemes. + * This computes the stress tensor and surface forces resulting from mechanical deformation. + */ +#ifndef DUMUX_DISCRETIZATION_HOOKES_LAW_HH +#define DUMUX_DISCRETIZATION_HOOKES_LAW_HH + +#include <dumux/discretization/methods.hh> + +namespace Dumux { + +/*! + * \ingroup Discretization + * \brief Evaluates the normal component of the Darcy velocity on a (sub)control volume face. + * \note Specializations are provided for the different discretization methods. + * These specializations are found in the headers included below. + */ +template <class Scalar, class FVGridGeometry, DiscretizationMethod dm = FVGridGeometry::discMethod> +class HookesLaw; + +} // end namespace Dumux + +#include <dumux/discretization/box/hookeslaw.hh> + +#endif diff --git a/dumux/geomechanics/CMakeLists.txt b/dumux/geomechanics/CMakeLists.txt new file mode 100644 index 0000000000..6cf6436721 --- /dev/null +++ b/dumux/geomechanics/CMakeLists.txt @@ -0,0 +1,9 @@ +add_subdirectory("elastic") + +install(FILES +fvproblem.hh +lameparams.hh +properties.hh +stressvariablescache.hh +velocityoutput.hh +DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/geomechanics) diff --git a/dumux/geomechanics/elastic/CMakeLists.txt b/dumux/geomechanics/elastic/CMakeLists.txt new file mode 100644 index 0000000000..d37ff68eec --- /dev/null +++ b/dumux/geomechanics/elastic/CMakeLists.txt @@ -0,0 +1,6 @@ +install(FILES +indices.hh +localresidual.hh +model.hh +volumevariables.hh +DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/geomechanics/elastic) diff --git a/dumux/geomechanics/elastic/indices.hh b/dumux/geomechanics/elastic/indices.hh new file mode 100644 index 0000000000..782d06c20e --- /dev/null +++ b/dumux/geomechanics/elastic/indices.hh @@ -0,0 +1,58 @@ +// -*- 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 Geomechanics + * \ingroup Elastic + * \brief Defines the indices for the elastic model + */ +#ifndef DUMUX_ELASTIC_INDICES_HH +#define DUMUX_ELASTIC_INDICES_HH + +namespace Dumux { +// \{ + +/*! + * \ingroup Geomechanics + * \ingroup Elastic + * \brief The indices for the linear elasticity model. + */ +struct ElasticIndices +{ + // returns the equation index for a given space direction + static constexpr int momentum(int dirIdx) { return dirIdx; }; + + // returns the primary variable index for a given space direction + static constexpr int u(int dirIdx) { return dirIdx; }; + + // Equation indices + static constexpr int momentumXEqIdx = 0; + static constexpr int momentumYEqIdx = 1; + static constexpr int momentumZEqIdx = 2; + + // primary variable indices + static constexpr int uxIdx = 0; + static constexpr int uyIdx = 1; + static constexpr int uzIdx = 2; +}; + +// \} +} // end namespace + +#endif diff --git a/dumux/geomechanics/elastic/localresidual.hh b/dumux/geomechanics/elastic/localresidual.hh new file mode 100644 index 0000000000..99f6cf0fc7 --- /dev/null +++ b/dumux/geomechanics/elastic/localresidual.hh @@ -0,0 +1,151 @@ +// -*- 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 Geomechanics + * \ingroup Elastic + * \brief Element-wise calculation of the local residual for problems + * using the elastic model considering linear elasticity. + */ +#ifndef DUMUX_GEOMECHANICS_ELASTIC_LOCAL_RESIDUAL_HH +#define DUMUX_GEOMECHANICS_ELASTIC_LOCAL_RESIDUAL_HH + +#include <vector> +#include <dumux/common/properties.hh> + +namespace Dumux { + +/*! + * \ingroup Geomechanics + * \ingroup Elastic + * \brief Element-wise calculation of the local residual for problems + * using the elastic model considering linear elasticity. + */ +template<class TypeTag> +class ElasticLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual) +{ + using ParentType = typename GET_PROP_TYPE(TypeTag, BaseLocalResidual); + using Implementation = typename GET_PROP_TYPE(TypeTag, LocalResidual); + + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using Element = typename GridView::template Codim<0>::Entity; + + using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices; + using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView; + using SubControlVolume = typename FVElementGeometry::SubControlVolume; + using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; + using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector); + using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, GridFluxVariablesCache)::LocalView; + using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView; + using VolumeVariables = typename ElementVolumeVariables::VolumeVariables; + + // class assembling the stress tensor + using StressType = typename GET_PROP_TYPE(TypeTag, StressType); + +public: + using ParentType::ParentType; + + /*! + * \brief For the elastic model the storage term is zero since + * we neglect inertia forces. + * + * \param problem The problem + * \param scv The sub control volume + * \param volVars The current or previous volVars + */ + NumEqVector computeStorage(const Problem& problem, + const SubControlVolume& scv, + const VolumeVariables& volVars) const + { + return NumEqVector(0.0); + } + + /*! + * \brief Evaluates the force in all grid directions acting + * on a face of a sub-control volume. + * + * \param problem The problem + * \param element The current element. + * \param fvGeometry The finite-volume geometry + * \param elemVolVars The volume variables of the current element + * \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 + { + // obtain force on the face from stress type + return StressType::force(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache); + } + + /*! + * \brief Calculate the source term of the equation + * + * \param problem The problem to solve + * \param element The DUNE Codim<0> entity for which the residual + * ought to be calculated + * \param fvGeometry The finite-volume geometry of the element + * \param elemVolVars The volume variables associated with the element stencil + * \param scv The sub-control volume over which we integrate the source term + * \note This is the default implementation for geomechanical models adding to + * the user defined sources the source stemming from the gravitational acceleration. + * + */ + NumEqVector computeSource(const Problem& problem, + const Element& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const SubControlVolume &scv) const + { + NumEqVector source(0.0); + + // add contributions from volume flux sources + source += problem.source(element, fvGeometry, elemVolVars, scv); + + // add contribution from possible point sources + source += problem.scvPointSources(element, fvGeometry, elemVolVars, scv); + + // maybe add gravitational acceleration + static const bool gravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity"); + if (gravity) + { + const auto& g = problem.gravityAtPos(scv.center()); + for (int dir = 0; dir < GridView::dimensionworld; ++dir) + source[Indices::momentum(dir)] += elemVolVars[scv].solidDensity()*g[dir]; + } + + return source; + } + +protected: + Implementation *asImp_() + { return static_cast<Implementation *> (this); } + + const Implementation *asImp_() const + { return static_cast<const Implementation *> (this); } +}; + +} // end namespace Dumux + +#endif diff --git a/dumux/geomechanics/elastic/model.hh b/dumux/geomechanics/elastic/model.hh new file mode 100644 index 0000000000..5a6dad5b41 --- /dev/null +++ b/dumux/geomechanics/elastic/model.hh @@ -0,0 +1,113 @@ +// -*- 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 Properties + * \ingroup Geomechanics + * \ingroup Elastic + * \brief Defines a type tag and some properties for the elastic geomechanical model + */ +#ifndef DUMUX_GEOMECHANICS_ELASTIC_MODEL_HH +#define DUMUX_GEOMECHANICS_ELASTIC_MODEL_HH + +#include <dumux/common/properties.hh> +#include <dumux/common/properties/model.hh> + +#include <dumux/geomechanics/properties.hh> +#include <dumux/discretization/hookeslaw.hh> + +#include "indices.hh" +#include "localresidual.hh" +#include "volumevariables.hh" + +namespace Dumux { + +/*! + * \ingroup Geomechanics + * \ingroup Elastic + * \brief Specifies a number properties of the elastic model + */ +template< int dim, int numSolidComp > +struct ElasticModelTraits +{ + //! export the type encapsulating indices + using Indices = ElasticIndices; + //! the number of equations is equal to grid dimension + static constexpr int numEq() { return dim; } + //! This model does not consider fluid phases + static constexpr int numFluidPhases() { return 0; } + //! This model does not consider fluid phases + static constexpr int numFluidComponents() { return 0; } + //! We have one solid phase here + static constexpr int numSolidComponents() { return numSolidComp; } + + //! Energy balance not yet implemented + static constexpr bool enableEnergyBalance() { return false; } +}; + +/*! + * \ingroup Elastic + * \brief Traits class for the volume variables of the elastic model. + * + * \tparam PV The type used for primary variables + * \tparam MT The model traits + * \tparam LP The class used for storing the lame parameters + * \tparam SST The solid state + * \tparam SSY The solid system + */ +template<class PV, class MT, class LP, class SST, class SSY> +struct ElasticVolumeVariablesTraits +{ + using PrimaryVariables = PV; + using ModelTraits = MT; + using LameParams = LP; + using SolidState = SST; + using SolidSystem = SSY; +}; + +namespace Properties { + +//! Type tag for the elastic geomechanical model +NEW_TYPE_TAG(Elastic, INHERITS_FROM(Geomechanics)); + +//! Use the local residual of the elastic model +SET_TYPE_PROP(Elastic, LocalResidual, ElasticLocalResidual<TypeTag>); + +//! By default, we use hooke's law for stress evaluations +SET_TYPE_PROP(Elastic, ModelTraits, ElasticModelTraits< GET_PROP_TYPE(TypeTag, GridView)::dimension, + GET_PROP_TYPE(TypeTag, SolidSystem)::numComponents >); + +//! Set the volume variables property +SET_PROP(Elastic, VolumeVariables) +{ +private: + using PV = typename GET_PROP_TYPE(TypeTag, PrimaryVariables); + using MT = typename GET_PROP_TYPE(TypeTag, ModelTraits); + using LP = typename GET_PROP_TYPE(TypeTag, SpatialParams)::LameParams; + using SST = typename GET_PROP_TYPE(TypeTag, SolidState); + using SSY = typename GET_PROP_TYPE(TypeTag, SolidSystem); + using Traits = ElasticVolumeVariablesTraits<PV, MT, LP, SST, SSY>; +public: + using type = ElasticVolumeVariables<Traits>; +}; + +} // namespace Properties +} // namespace Dumux + + #endif diff --git a/dumux/geomechanics/elastic/volumevariables.hh b/dumux/geomechanics/elastic/volumevariables.hh new file mode 100644 index 0000000000..80ef4b0241 --- /dev/null +++ b/dumux/geomechanics/elastic/volumevariables.hh @@ -0,0 +1,121 @@ +// -*- 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 Geomechanics + * \ingroup Elastic + * \brief Quantities required by the elastic model defined on a sub-control volume. + */ +#ifndef DUMUX_ELASTIC_VOLUME_VARIABLES_HH +#define DUMUX_ELASTIC_VOLUME_VARIABLES_HH + +#include <dumux/material/solidstates/updatesolidvolumefractions.hh> + +namespace Dumux { + +/*! +* \ingroup Geomechanics +* \ingroup Elastic + * \brief Contains the quantities which are constant within a + * finite volume in the elastic model. + * + * \tparam Traits Class encapsulating types to be used by the vol vars + */ +template<class Traits> +class ElasticVolumeVariables +{ + using Scalar = typename Traits::PrimaryVariables::value_type; + using LameParams = typename Traits::LameParams; + using ModelTraits = typename Traits::ModelTraits; + + //! The elastic model only makes sense with inert solid systems + static_assert(Traits::SolidSystem::isInert(), "Elastic model can only be used with inert solid systems"); +public: + //! export the type used for the primary variables + using PrimaryVariables = typename Traits::PrimaryVariables; + //! export the type encapsulating primary variable indices + using Indices = typename ModelTraits::Indices; + //! export type of solid state + using SolidState = typename Traits::SolidState; + //! export the solid system used + 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) + { + priVars_ = elemSol[scv.localDofIndex()]; + extrusionFactor_ = problem.extrusionFactor(element, scv, elemSol); + lameParams_ = problem.spatialParams().lameParams(element, elemSol); + + //! set the volume fractions of the solid components + updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, /*numFluidComps=*/0); + // set the temperature of the solid phase + setSolidTemperature_(problem, elemSol); + // update the density of the solid phase + solidState_.setDensity(SolidSystem::density(solidState_)); + }; + + //! Return the average porosity \f$\mathrm{[-]}\f$ within the control volume. + Scalar solidDensity() const { return solidState_.density(); } + //! Returns the permeability within the control volume in \f$[m^2]\f$. + Scalar displacement(unsigned int dir) const { return priVars_[ Indices::momentum(dir) ]; } + //! Return a component of primary variable vector for a given index + Scalar priVar(const int pvIdx) const { return priVars_[pvIdx]; } + //! Return the vector of primary variables + const PrimaryVariables& priVars() const { return priVars_; } + //! Return the object containing the lame params of this scv + const LameParams& lameParams() const { return lameParams_; } + //! TODO We don't know yet how to interpret extrusion for mechanics + static constexpr Scalar extrusionFactor() { return 1.0; } + +private: + //! sets the temperature in the solid state for non-isothermal models + template< class Problem, class ElemSol, + bool EB = ModelTraits::enableEnergyBalance(), std::enable_if_t<EB, int> = 0 > + void setSolidTemperature_(const Problem& problem, const ElemSol& elemSol) + { DUNE_THROW(Dune::InvalidStateException, "Non-isothermal elastic model."); } + + //! sets the temperature in the solid state for isothermal models + template< class Problem, class ElemSol, + bool EB = ModelTraits::enableEnergyBalance(), std::enable_if_t<!EB, int> = 0 > + void setSolidTemperature_(const Problem& problem, const ElemSol& elemSol) + { solidState_.setTemperature(problem.temperature()); } + + // data members + Scalar extrusionFactor_; + PrimaryVariables priVars_; + LameParams lameParams_; + SolidState solidState_; +}; + +} + +#endif diff --git a/dumux/geomechanics/fvproblem.hh b/dumux/geomechanics/fvproblem.hh new file mode 100644 index 0000000000..dd180b1e8e --- /dev/null +++ b/dumux/geomechanics/fvproblem.hh @@ -0,0 +1,49 @@ +// -*- 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 Base class for all geomechanical problems + */ +#ifndef DUMUX_GEOMECHANICS_FV_PROBLEM_HH +#define DUMUX_GEOMECHANICS_FV_PROBLEM_HH + +#include <dumux/porousmediumflow/problem.hh> + +namespace Dumux { + +/*! + * \ingroup Geomechanics + * \brief Base class for all geomechanical problems + * \note We actually require the same functionality as the + * porous medium flow problem, which is why we simply + * inherit from that here. + */ +template<class TypeTag> +class GeomechanicsFVProblem : public PorousMediumFlowProblem<TypeTag> +{ + using ParentType = PorousMediumFlowProblem<TypeTag>; + +public: + //! pull up the constructor of the parent class + using ParentType::ParentType; +}; + +} // end namespace Dumux + +#endif diff --git a/dumux/geomechanics/lameparams.hh b/dumux/geomechanics/lameparams.hh new file mode 100644 index 0000000000..04b9fea971 --- /dev/null +++ b/dumux/geomechanics/lameparams.hh @@ -0,0 +1,57 @@ +// -*- 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 \copydoc Dumux::LameParams + */ +#ifndef DUMUX_GEOMECHANICS_LAME_PARAMS_HH +#define DUMUX_GEOMECHANICS_LAME_PARAMS_HH + +namespace Dumux +{ + +/*! + * \ingroup Geomechanics + * \ingroup Elastic + * \brief Structure encapsulating the lame parameters + */ +template<class Scalar> +struct LameParams +{ + //! Default constructor + LameParams() = default; + //! Constructor taking lambda and mu directly + LameParams(Scalar lambda, Scalar mu) : lambda_(lambda), mu_(mu) {} + + //! Return the first lame parameter + Scalar lambda() const { return lambda_; } + //! Return the second lame parameter + Scalar mu() const { return mu_; } + + //! set the first lame parameter + void setLambda(Scalar lambda) { lambda_ = lambda; } + //! set the second lame parameter + void setMu(Scalar mu) { mu_ = mu; } + +private: + Scalar lambda_; + Scalar mu_; +}; +} +#endif diff --git a/dumux/geomechanics/properties.hh b/dumux/geomechanics/properties.hh new file mode 100644 index 0000000000..34b73e92e9 --- /dev/null +++ b/dumux/geomechanics/properties.hh @@ -0,0 +1,78 @@ +// -*- 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 Properties + * \ingroup Geomechanics + * \brief Defines a type tag and some properties for geomechanical DuMuX models. + */ + +#ifndef DUMUX_GEOMECHANICS_PROPERTIES_HH +#define DUMUX_GEOMECHANICS_PROPERTIES_HH + +#include <dumux/common/properties.hh> +#include <dumux/common/properties/model.hh> +#include <dumux/material/components/constant.hh> +#include <dumux/material/solidstates/inertsolidstate.hh> +#include <dumux/material/solidsystems/inertsolidphase.hh> +#include <dumux/discretization/hookeslaw.hh> + +#include "stressvariablescache.hh" +#include "velocityoutput.hh" + +namespace Dumux { +namespace Properties { + +//! Type tag for geomechanical models +NEW_TYPE_TAG(Geomechanics, INHERITS_FROM(ModelProperties)); + +//! The flux variables cache class for models involving flow in porous media +SET_TYPE_PROP(Geomechanics, FluxVariablesCache, StressVariablesCache< typename GET_PROP_TYPE(TypeTag, Scalar), + typename GET_PROP_TYPE(TypeTag, FVGridGeometry) >); + +//! By default, we use hooke's law for stress evaluations +SET_TYPE_PROP(Geomechanics, StressType, HookesLaw< typename GET_PROP_TYPE(TypeTag, Scalar), + typename GET_PROP_TYPE(TypeTag, FVGridGeometry) >); + +//! The (currently empty) velocity output +SET_TYPE_PROP(Geomechanics, VelocityOutput, GeomechanicsVelocityOutput); + +//! The solid state must be inert +SET_PROP(Geomechanics, SolidState) +{ +private: + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using SolidSystem = typename GET_PROP_TYPE(TypeTag, SolidSystem); +public: + using type = InertSolidState<Scalar, SolidSystem>; +}; + +//! Per default we use one constant component in the inert solid system +SET_PROP(Geomechanics, SolidSystem) +{ +private: + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using InertComponent = Components::Constant<1, Scalar>; +public: + using type = SolidSystems::InertSolidPhase<Scalar, InertComponent>; +}; +} // namespace Properties +} // namespace Dumux + + #endif diff --git a/dumux/geomechanics/stressvariablescache.hh b/dumux/geomechanics/stressvariablescache.hh new file mode 100644 index 0000000000..cd192aa4e7 --- /dev/null +++ b/dumux/geomechanics/stressvariablescache.hh @@ -0,0 +1,66 @@ +// -*- 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 Geomechanics + * \brief Base class for the stress variables cache + */ +#ifndef DUMUX_GEOMECHANICS_STRESSVARIABLESCACHE_HH +#define DUMUX_GEOMECHANICS_STRESSVARIABLESCACHE_HH + +#include <dune/common/exceptions.hh> + +#include <dumux/discretization/methods.hh> +#include <dumux/discretization/fluxvariablescaching.hh> +#include <dumux/discretization/box/fluxvariablescache.hh> + +namespace Dumux { + +/*! + * \ingroup Geomechanics + * \brief The stress variables cache classes for models involving geomechanics. + * Store data required for stress calculation. + */ +template< class Scalar, class FVGridGeometry, DiscretizationMethod dm = FVGridGeometry::discMethod > +class StressVariablesCache; + +//! We only store discretization-related quantities for the box method. +template< class Scalar, class FVGridGeometry > +class StressVariablesCache<Scalar, FVGridGeometry, DiscretizationMethod::box> +: public BoxFluxVariablesCache< Scalar, FVGridGeometry > +{}; + +// specialization for the cell centered tpfa method +template< class Scalar, class FVGridGeometry > +class StressVariablesCache<Scalar, FVGridGeometry, DiscretizationMethod::cctpfa> : public FluxVariablesCaching::_EmptyCache +{ +public: + template<typename... Args> + void update(Args&&... args) { DUNE_THROW(Dune::NotImplemented, "Geomechanics with cell-centered schemes"); } +}; + +// specialization for the cell centered mpfa method +template< class Scalar, class FVGridGeometry > +class StressVariablesCache<Scalar, FVGridGeometry, DiscretizationMethod::ccmpfa> +: public StressVariablesCache<Scalar, FVGridGeometry, DiscretizationMethod::cctpfa> +{}; + +} // end namespace Dumux + +#endif diff --git a/dumux/geomechanics/velocityoutput.hh b/dumux/geomechanics/velocityoutput.hh new file mode 100644 index 0000000000..dd4571eb49 --- /dev/null +++ b/dumux/geomechanics/velocityoutput.hh @@ -0,0 +1,60 @@ +// -*- 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 Velocity output for geomechanical models + */ +#ifndef DUMUX_GEOMECHANICS_VELOCITYOUTPUT_HH +#define DUMUX_GEOMECHANICS_VELOCITYOUTPUT_HH + +namespace Dumux { + +/*! + * \brief Velocity output for geomechanical models. + * This class could be used to compute the temporal derivative + * of the displacement. Currently this is not implemented and + * we simply define this here in order to be able to reuse the + * VtkOutputModule which expects a VelocityOutput class. + */ +class GeomechanicsVelocityOutput +{ +public: + //! The constructor + template< typename... Args > + GeomechanicsVelocityOutput(Args&&... args) {} + + //! Output is currently disabled (not implemented) + static constexpr bool enableOutput() { return false; } + + //! There is always only one solid phase + static constexpr int numPhaseVelocities() { return 1; } + + //! Returns the name of phase for which velocity is computed + static std::string phaseName(int phaseIdx) + { DUNE_THROW(Dune::NotImplemented, "Velocity output for geomechanical models."); } + + //! Calculate the velocities for the scvs in the element + template< typename... Args > + void calculateVelocity(Args... args) + { DUNE_THROW(Dune::NotImplemented, "Velocity output for geomechanical models."); } +}; + +} // end namespace Dumux + +#endif diff --git a/dumux/material/spatialparams/fvelastic.hh b/dumux/material/spatialparams/fvelastic.hh new file mode 100644 index 0000000000..2a2430cdb3 --- /dev/null +++ b/dumux/material/spatialparams/fvelastic.hh @@ -0,0 +1,143 @@ +// -*- 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 The base class for spatial parameters of linear elastic geomechanical problems + */ +#ifndef DUMUX_GEOMECHANICS_ELASTIC_FV_SPATIAL_PARAMS_HH +#define DUMUX_GEOMECHANICS_ELASTIC_FV_SPATIAL_PARAMS_HH + +#include <dumux/common/typetraits/isvalid.hh> +#include <dumux/geomechanics/lameparams.hh> + +namespace Dumux { + +#ifndef DOXYGEN +namespace Detail { +// helper struct detecting if the user-defined spatial params class has a lameParamsAtPos function +// for g++ > 5.3, this can be replaced by a lambda +template<class GlobalPosition> +struct hasLameParamsAtPos +{ + template<class SpatialParams> + auto operator()(const SpatialParams& a) + -> decltype(a.lameParamsAtPos(std::declval<GlobalPosition>())) + {}; +}; + +} // end namespace Detail +#endif + +/*! + * \ingroup Geomechanics + * \brief The base class for spatial parameters of linear elastic geomechanical problems + */ +template<class Scalar, class FVGridGeometry, class Implementation> +class FVSpatialParamsElastic +{ + using SubControlVolume = typename FVGridGeometry::SubControlVolume; + using GridView = typename FVGridGeometry::GridView; + using Element = typename GridView::template Codim<0>::Entity; + using GlobalPosition = typename Element::Geometry::GlobalCoordinate; + +public: + //! The constructor + FVSpatialParamsElastic(std::shared_ptr<const FVGridGeometry> fvGridGeometry) + : fvGridGeometry_(fvGridGeometry) + {} + + /*! + * \brief Function for defining the solid volume fraction. + * That is possibly solution dependent. + * + * \param element The current element + * \param scv The sub-control volume inside the element. + * \param elemSol The solution at the dofs connected to the element. + * \param compIdx The solid component index + * \return the volume fraction of the inert solid component with index compIdx + */ + template<class SolidSystem, class ElementSolution> + Scalar inertVolumeFraction(const Element& element, + const SubControlVolume& scv, + const ElementSolution& elemSol, + int compIdx) const + { + static_assert(SolidSystem::isInert(), "Elastic model can only be used with inert solid systems"); + + // when there is only one component, the volume fraction must be one + if (SolidSystem::numInertComponents == 1) + return 1.0; + + // otherwise we require the user to define the solid composition + return asImp_().template inertVolumeFractionAtPos<SolidSystem>(scv.center(), compIdx); + } + + /*! + * \brief Function for defining the solid volume fraction. + * That is possibly solution dependent. + * + * \param globalPos The global position + * \param compIdx The solid component index + * \return the volume fraction of the inert solid component with index compIdx + */ + template<class SolidSystem> + Scalar inertVolumeFractionAtPos(const GlobalPosition& globalPos, int compIdx) const + { DUNE_THROW(Dune::InvalidStateException, "The spatial parameters do not provide inertVolumeFractionAtPos() method."); } + + /*! + * \brief Define the Lame parameters + * \note It is possibly solution dependent. + * + * \param element The current element + * \param elemSol The solution at the dofs connected to the element. + * \return lame parameters + * \todo TODO Could the lame parameters also be heterogeneous inside element + * (i.e. pass scv as additional argument to this function)? + * This would need appropriate adjustments in Hooke's law. + */ + template<class ElementSolution> + decltype(auto) lameParams(const Element& element, + const ElementSolution& elemSol) const + { + static_assert(decltype(isValid(Detail::hasLameParamsAtPos<GlobalPosition>())(this->asImp_()))::value," \n\n" + " Your spatial params class has to either implement\n\n" + " const LameParams& lameParamsAtPos(const GlobalPosition& globalPos) const\n\n" + " or overload this function\n\n" + " template<class ElementSolution>\n" + " const LameParams& lameParams(const Element& element,\n" + " const ElementSolution& elemSol) const\n\n"); + + return asImp_().lameParamsAtPos(element.geometry().center()); + } + + //! The finite volume grid geometry + const FVGridGeometry& fvGridGeometry() const { return *fvGridGeometry_; } + +protected: + Implementation &asImp_() + { return *static_cast<Implementation*>(this); } + + const Implementation &asImp_() const + { return *static_cast<const Implementation*>(this); } + +private: + std::shared_ptr<const FVGridGeometry> fvGridGeometry_; +}; +} +#endif -- GitLab