From 7eb6507a153c1ba4dc39776da6635876672f3a23 Mon Sep 17 00:00:00 2001 From: Martin Utz Date: Wed, 16 Mar 2022 11:33:30 +0100 Subject: [PATCH 1/2] [swe][friction] Fix several things about the 'shearStress' interface * shearstress is deprecated (it returned bottom shear stress normalized by negative density) * new interface 'bottomShearStress' that returns the (positive) bottom shear stress in N/m^2 * introduce FluidSystem for shallow water models * Make the assumption of incompressible fluids explicit by static_asserts --- CHANGELOG.md | 1 + dumux/flux/shallowwaterviscousflux.hh | 12 +++++-- dumux/freeflow/shallowwater/model.hh | 17 +++++++++- .../freeflow/shallowwater/volumevariables.hh | 17 ++++++++++ .../frictionlaws/frictionlaw.hh | 15 +++++---- .../frictionlaws/manning.hh | 31 ++++++++++--------- .../frictionlaws/nikuradse.hh | 24 +++++++------- .../frictionlaws/nofriction.hh | 13 ++++---- examples/shallowwaterfriction/problem.hh | 12 +++---- .../shallowwater/poiseuilleflow/params.input | 0 .../shallowwater/roughchannel/problem.hh | 6 ++-- 11 files changed, 98 insertions(+), 50 deletions(-) mode change 100755 => 100644 test/freeflow/shallowwater/poiseuilleflow/params.input diff --git a/CHANGELOG.md b/CHANGELOG.md index 8b727bebfb..2d78078226 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -101,6 +101,7 @@ The parameters describing your dispersion tensor can then be included in your `s - Porous medium flow models should now inherit from the new base spatial parameters that can be found in the folder `dumux/porousmediumflow/`, which allow users to overload the new `temperature` and `extrusionFactor` interfaces. - Free flow and pore network models now also expect the user problems to expose spatial parameters, in which `gravity`, `temperature` and `extrusionFactor` are defined. The respective problem interfaces have been deprecated. - The problem interfaces for fluid properties in the poroelastic model, namely `effectiveFluidDensity` and `effectivePorePressure`, have been deprecated and were moved to the spatial parameters. +- The function `shearStress` in the class `FrictionLaw` and its child classes has been renamed to `bottomShearStress` and the return value has been changed. The function returns now the actual bottom shear stress and not the bottom shear stress term as it used in the shallow water model. The bottom shear stress term of the shallow water model is the bottom shear stress multiplied with minus one and normalised by the water density. ### New experimental features (possibly subject to backwards-incompatible changes in the future) diff --git a/dumux/flux/shallowwaterviscousflux.hh b/dumux/flux/shallowwaterviscousflux.hh index 072a6c5d42..c1f8d8ec21 100644 --- a/dumux/flux/shallowwaterviscousflux.hh +++ b/dumux/flux/shallowwaterviscousflux.hh @@ -87,6 +87,7 @@ public: const typename FVElementGeometry::SubControlVolumeFace& scvf) { using Scalar = typename NumEqVector::value_type; + using FluidSystem = typename ElementVolumeVariables::VolumeVariables::FluidSystem; NumEqVector localFlux(0.0); @@ -160,8 +161,8 @@ public: { // Get the bottom shear stress in the two adjacent cells // Note that element is not needed in spatialParams().frictionLaw (should be removed). For now we simply pass the same element twice - const auto& bottomShearStressInside = problem.spatialParams().frictionLaw(element, insideScv).shearStress(insideVolVars); - const auto& bottomShearStressOutside = problem.spatialParams().frictionLaw(element, outsideScv).shearStress(outsideVolVars); + const auto& bottomShearStressInside = problem.spatialParams().frictionLaw(element, insideScv).bottomShearStress(insideVolVars); + const auto& bottomShearStressOutside = problem.spatialParams().frictionLaw(element, outsideScv).bottomShearStress(outsideVolVars); const auto bottomShearStressInsideMag = bottomShearStressInside.two_norm(); const auto bottomShearStressOutsideMag = bottomShearStressOutside.two_norm(); @@ -172,7 +173,12 @@ public: // Shear velocity possibly needed in mixing-length turbulence model in the computation of the diffusion flux using std::sqrt; - return sqrt(averageBottomShearStress); + static_assert(!FluidSystem::isCompressible(0), + "The shallow water model assumes incompressible fluids" + ); + // Since the shallow water equations are incompressible, the density is constant. + // Therefore, insideVolVars.density() is equal to outsideVolVars.density(). + return sqrt(averageBottomShearStress/insideVolVars.density()); }(); const auto turbViscosityV = turbConstV * (kappa/6.0) * ustar * averageDepth; diff --git a/dumux/freeflow/shallowwater/model.hh b/dumux/freeflow/shallowwater/model.hh index 12de3eacf7..db984d34c9 100644 --- a/dumux/freeflow/shallowwater/model.hh +++ b/dumux/freeflow/shallowwater/model.hh @@ -69,6 +69,8 @@ #include #include #include +#include +#include #include "localresidual.hh" #include "volumevariables.hh" @@ -98,13 +100,16 @@ struct ShallowWaterModelTraits * \brief Traits class for the volume variables of the shallow water model. * * \tparam PV The type used for primary variables + * \tparam FSY The fluid system type * \tparam MT The model traits */ template struct ShallowWaterVolumeVariablesTraits { using PrimaryVariables = PV; + using FluidSystem = FSY; using ModelTraits = MT; }; @@ -137,8 +142,9 @@ struct VolumeVariables { private: using PV = GetPropType; + using FSY = GetPropType; using MT = GetPropType; - using Traits = ShallowWaterVolumeVariablesTraits; + using Traits = ShallowWaterVolumeVariablesTraits; public: using type = ShallowWaterVolumeVariables; }; @@ -163,6 +169,15 @@ template struct ViscousFluxType { using type = ShallowWaterViscousFlux< Dumux::NumEqVector> >; }; +template +struct FluidSystem +{ +private: + using Scalar = GetPropType; +public: + using type = FluidSystems::OnePLiquid>; +}; + } // end properties } // end namespace Dumux diff --git a/dumux/freeflow/shallowwater/volumevariables.hh b/dumux/freeflow/shallowwater/volumevariables.hh index 52450b5d74..7d19df39cb 100644 --- a/dumux/freeflow/shallowwater/volumevariables.hh +++ b/dumux/freeflow/shallowwater/volumevariables.hh @@ -38,6 +38,8 @@ class ShallowWaterVolumeVariables public: using PrimaryVariables = typename Traits::PrimaryVariables; + //! export the underlying fluid system + using FluidSystem = typename Traits::FluidSystem; template void update(const ElemSol &elemSol, @@ -90,6 +92,21 @@ public: return bedSurface_; } + /*! + * \brief Return the fluid density + * + */ + Scalar density(int phaseIdx = 0) const + { + static_assert(!FluidSystem::isCompressible(0), + "The shallow water model assumes incompressible fluids" + ); + + // call with hard-coded sensible default values for water/river applications for now + return FluidSystem::density(283.15, 1e5); + + } + private: PrimaryVariables priVars_; Scalar bedSurface_; diff --git a/dumux/material/fluidmatrixinteractions/frictionlaws/frictionlaw.hh b/dumux/material/fluidmatrixinteractions/frictionlaws/frictionlaw.hh index 31499ee087..465eb4948a 100644 --- a/dumux/material/fluidmatrixinteractions/frictionlaws/frictionlaw.hh +++ b/dumux/material/fluidmatrixinteractions/frictionlaws/frictionlaw.hh @@ -43,18 +43,21 @@ class FrictionLaw using Scalar = typename VolumeVariables::PrimaryVariables::value_type; public: /*! - * \brief Compute the shear stress. + * \brief Compute the bottom shear stress. * - * \param volVar Volume Variables. + * \param volVars Volume variables * - * Compute the shear stress due to friction. The shear stress is not a tensor as know - * from contiuums mechanics, but a force projected on an area. Therefore it is a - * vector with two entries. + * Compute the bottom shear stress due to bottom friction. + * The bottom shear stress is a projection of the shear stress tensor onto the river bed. + * It can therefore be represented by a (tangent) vector with two entries. * * \return shear stress [N/m^2]. First entry is the x-component, the second the y-component. */ + virtual Dune::FieldVector bottomShearStress(const VolumeVariables& volVars) const = 0; - virtual Dune::FieldVector shearStress(const VolumeVariables& volVar) const = 0; + [[deprecated("Use bottomShearStress. Note that the unit and sign of the return value is different. Will be removed after release 3.5")]] + virtual Dune::FieldVector shearStress(const VolumeVariables& volVars) const + { auto bss = bottomShearStress(volVars); bss /= -volVars.density(); return bss; } /*! * \brief Limit the friction for small water depth. diff --git a/dumux/material/fluidmatrixinteractions/frictionlaws/manning.hh b/dumux/material/fluidmatrixinteractions/frictionlaws/manning.hh index a0b8a369ec..51e70977ee 100644 --- a/dumux/material/fluidmatrixinteractions/frictionlaws/manning.hh +++ b/dumux/material/fluidmatrixinteractions/frictionlaws/manning.hh @@ -45,24 +45,24 @@ public: /*! * \brief Constructor * - * \param gravity Gravity constant [m/s^2] - * \param manningN Manning friction coefficient [-] + * \param gravity Gravity constant (in m/s^2) + * \param manningN Manning friction coefficient (in s/m^(1/3)) */ FrictionLawManning(const Scalar gravity, const Scalar manningN) - : gravity_(gravity), manningN_(manningN) {} + : gravity_(gravity), manningN_(manningN) {} /*! - * \brief Compute the shear stress. + * \brief Compute the bottom shear stress. * * \param volVars Volume variables * - * Compute the shear stress due to friction. The shear stress is not a tensor as know - * from contiuums mechanics, but a force projected on an area. Therefore it is a - * vector with two entries. + * Compute the bottom shear stress due to bottom friction. + * The bottom shear stress is a projection of the shear stress tensor onto the river bed. + * It can therefore be represented by a (tangent) vector with two entries. * - * \return shear stress [N/m^2]. First entry is the x-component, the second the y-component. + * \return shear stress in N/m^2. First entry is the x-component, the second the y-component. */ - Dune::FieldVector shearStress(const VolumeVariables& volVars) const final + Dune::FieldVector bottomShearStress(const VolumeVariables& volVars) const final { using std::pow; using Dune::power; @@ -70,16 +70,19 @@ public: Dune::FieldVector shearStress(0.0); - Scalar roughnessHeight = power(25.68/(1.0/manningN_),6); + Scalar roughnessHeight = power(25.68/(1.0/manningN_), 6); roughnessHeight = this->limitRoughH(roughnessHeight, volVars.waterDepth()); - const Scalar c = pow((volVars.waterDepth() + roughnessHeight),1.0/6.0) * 1.0/(manningN_); - const Scalar uv = hypot(volVars.velocity(0),volVars.velocity(1)); + // c has units of m^(1/2)/s so c^2 has units of m/s^2 + const Scalar c = pow(volVars.waterDepth() + roughnessHeight, 1.0/6.0) * 1.0/(manningN_); + const Scalar uv = hypot(volVars.velocity(0), volVars.velocity(1)); - shearStress[0] = -gravity_/(c*c) * volVars.velocity(0) * uv; - shearStress[1] = -gravity_/(c*c) * volVars.velocity(1) * uv; + const Scalar dimensionlessFactor = gravity_/(c*c); + shearStress[0] = dimensionlessFactor * volVars.velocity(0) * uv * volVars.density(); + shearStress[1] = dimensionlessFactor * volVars.velocity(1) * uv * volVars.density(); return shearStress; } + private: Scalar gravity_; Scalar manningN_; diff --git a/dumux/material/fluidmatrixinteractions/frictionlaws/nikuradse.hh b/dumux/material/fluidmatrixinteractions/frictionlaws/nikuradse.hh index 10a8d15cc4..258341e716 100644 --- a/dumux/material/fluidmatrixinteractions/frictionlaws/nikuradse.hh +++ b/dumux/material/fluidmatrixinteractions/frictionlaws/nikuradse.hh @@ -45,23 +45,23 @@ public: /*! * \brief Constructor * - * \param ks Equivalent sand roughness [m] + * \param ks Equivalent sand roughness (in m) */ FrictionLawNikuradse(const Scalar ks) - : ks_(ks) {} + : ks_(ks) {} /*! - * \brief Compute the shear stress. + * \brief Compute the bottom shear stress. * * \param volVars Volume variables * - * Compute the shear stress due to friction. The shear stress is not a tensor as know - * from contiuums mechanics, but a force projected on an area. Therefore it is a - * vector with two entries. + * Compute the bottom shear stress due to bottom friction. + * The bottom shear stress is a projection of the shear stress tensor onto the river bed. + * It can therefore be represented by a (tangent) vector with two entries. * - * \return shear stress [N/m^2]. First entry is the x-component, the second the y-component. + * \return shear stress in N/m^2. First entry is the x-component, the second the y-component. */ - Dune::FieldVector shearStress(const VolumeVariables& volVars) const final + Dune::FieldVector bottomShearStress(const VolumeVariables& volVars) const final { using Dune::power; using std::log; @@ -71,14 +71,16 @@ public: Scalar roughnessHeight = ks_; roughnessHeight = this->limitRoughH(roughnessHeight, volVars.waterDepth()); - const Scalar ustarH = power(0.41,2)/power(log((12*(volVars.waterDepth() + roughnessHeight))/ks_),2); + const Scalar karmanConstant = 0.41; // Karman's constant is dimensionless + const Scalar dimensionlessFactor = power(karmanConstant, 2)/power(log((12*(volVars.waterDepth() + roughnessHeight))/ks_), 2); const Scalar uv = hypot(volVars.velocity(0),volVars.velocity(1)); - shearStress[0] = -ustarH * volVars.velocity(0) * uv; - shearStress[1] = -ustarH * volVars.velocity(1) * uv; + shearStress[0] = dimensionlessFactor * volVars.velocity(0) * uv * volVars.density(); + shearStress[1] = dimensionlessFactor * volVars.velocity(1) * uv * volVars.density(); return shearStress; } + private: Scalar ks_; }; diff --git a/dumux/material/fluidmatrixinteractions/frictionlaws/nofriction.hh b/dumux/material/fluidmatrixinteractions/frictionlaws/nofriction.hh index b2072fc417..851a0aa925 100644 --- a/dumux/material/fluidmatrixinteractions/frictionlaws/nofriction.hh +++ b/dumux/material/fluidmatrixinteractions/frictionlaws/nofriction.hh @@ -43,17 +43,18 @@ public: FrictionLawNoFriction() = default; /*! - * \brief Compute the shear stress. + * \brief Compute the bottom shear stress. * * \param volVars Volume variables * - * Compute the shear stress due to friction. The shear stress is not a tensor as know - * from contiuums mechanics, but a force projected on an area. Therefore it is a - * vector with two entries. For this law without friction, the shearStress is zero. + * Compute the bottom shear stress due to bottom friction. + * The bottom shear stress is a projection of the shear stress tensor onto the river bed. + * It can therefore be represented by a (tangent) vector with two entries. + * For this law without bottom friction, the bottom shear stress is zero. * - * \return shear stress [N/m^2]. First entry is the x-component, the second the y-component. + * \return shear stress in N/m^2. First entry is the x-component, the second the y-component. */ - Dune::FieldVector shearStress(const VolumeVariables& volVars) const final + Dune::FieldVector bottomShearStress(const VolumeVariables& volVars) const final { return {0.0, 0.0}; } diff --git a/examples/shallowwaterfriction/problem.hh b/examples/shallowwaterfriction/problem.hh index aff44d7233..eee68f25d2 100644 --- a/examples/shallowwaterfriction/problem.hh +++ b/examples/shallowwaterfriction/problem.hh @@ -161,20 +161,20 @@ public: // Accordingly, the third entry of the `bottomFrictionSource` is equal to the second component of the `bottomShearStress`. // [[codeblock]] NumEqVector bottomFrictionSource(const Element& element, - const FVElementGeometry& fvGeometry, - const ElementVolumeVariables& elemVolVars, - const SubControlVolume &scv) const + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const SubControlVolume &scv) const { NumEqVector bottomFrictionSource(0.0); const auto& volVars = elemVolVars[scv]; // bottom shear stress vector - Dune::FieldVector bottomShearStress = this->spatialParams().frictionLaw(element, scv).shearStress(volVars); + Dune::FieldVector bottomShearStress = this->spatialParams().frictionLaw(element, scv).bottomShearStress(volVars); // source term due to bottom friction bottomFrictionSource[0] = 0.0; - bottomFrictionSource[1] = bottomShearStress[0]; - bottomFrictionSource[2] = bottomShearStress[1]; + bottomFrictionSource[1] = -bottomShearStress[0] / volVars.density(); + bottomFrictionSource[2] = -bottomShearStress[1] / volVars.density(); return bottomFrictionSource; } diff --git a/test/freeflow/shallowwater/poiseuilleflow/params.input b/test/freeflow/shallowwater/poiseuilleflow/params.input old mode 100755 new mode 100644 diff --git a/test/freeflow/shallowwater/roughchannel/problem.hh b/test/freeflow/shallowwater/roughchannel/problem.hh index a43f11fe0b..837b7d83d2 100644 --- a/test/freeflow/shallowwater/roughchannel/problem.hh +++ b/test/freeflow/shallowwater/roughchannel/problem.hh @@ -202,11 +202,11 @@ public: NumEqVector bottomFrictionSource(0.0); const auto& volVars = elemVolVars[scv]; - Dune::FieldVector bottomShearStress = this->spatialParams().frictionLaw(element, scv).shearStress(volVars); + Dune::FieldVector bottomShearStress = this->spatialParams().frictionLaw(element, scv).bottomShearStress(volVars); bottomFrictionSource[0] = 0.0; - bottomFrictionSource[1] =bottomShearStress[0]; - bottomFrictionSource[2] =bottomShearStress[1]; + bottomFrictionSource[1] = -bottomShearStress[0] / volVars.density(); + bottomFrictionSource[2] = -bottomShearStress[1] / volVars.density(); return bottomFrictionSource; } -- GitLab From 73d0b6bd37165b4683babb69c3cce498acc49e70 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Thu, 7 Apr 2022 12:34:57 +0200 Subject: [PATCH 2/2] [swe] Fix documentation in friction law base class --- .../frictionlaws/frictionlaw.hh | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/dumux/material/fluidmatrixinteractions/frictionlaws/frictionlaw.hh b/dumux/material/fluidmatrixinteractions/frictionlaws/frictionlaw.hh index 465eb4948a..b063d080ac 100644 --- a/dumux/material/fluidmatrixinteractions/frictionlaws/frictionlaw.hh +++ b/dumux/material/fluidmatrixinteractions/frictionlaws/frictionlaw.hh @@ -63,8 +63,8 @@ public: * \brief Limit the friction for small water depth. * * We define a water depth minUpperH. If the water depth is - * smaller, we start to limit the friciton. - * So the friciton term get's not extreme large for small water + * smaller, we start to limit the friction. + * So the friction term get's not extreme large for small water * depths. * * ------------------------- minUpperH ----------- @@ -85,23 +85,22 @@ public: * For the limitation of the roughness height L = 0.0, T = 2.0 and E = 1.0 are choosen. * Therefore the calculation of the mobility is simplified significantly. * - * \param roughnessHeight roughness height of the representive structure (e.g. largest grain size). + * \param roughnessHeight roughness height of the representative structure (e.g. largest grain size). * \param waterDepth water depth. */ Scalar limitRoughH(const Scalar roughnessHeight, const Scalar waterDepth) const { - using std::min; - using std::max; + using std::clamp; - Scalar mobilityMax = 1.0; //!< maximal mobility + const Scalar mobilityMax = 1.0; //!< maximal mobility - Scalar minUpperH = roughnessHeight * 2.0; - Scalar sw = min(waterDepth * (1.0/minUpperH),1.0); - sw = max(0.0,sw); - auto mobility = mobilityMax /(1 + (1.0-sw)*(1.0-sw)); + const Scalar minUpperH = roughnessHeight * 2.0; + const Scalar sw = clamp(waterDepth * (1.0/minUpperH), 0.0, 1.0); + const Scalar mobility = mobilityMax /(1.0 + (1.0-sw)*(1.0-sw)); return roughnessHeight * (1.0 - mobility); } + // virtual base class needs a virtual destructor virtual ~FrictionLaw() = default; }; -- GitLab