diff --git a/CHANGELOG.md b/CHANGELOG.md index 8b727bebfb4160afa21515317ad19706fb4a29f0..2d780782266d70ad754924cf54da32d60ebe1951 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 072a6c5d42aa627979de68b35bf38f68d9f15e32..c1f8d8ec21183b32b9eb4a890872b98629b01eba 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 12de3eacf7646951ef895c4d075e83a1bc800404..db984d34c95c68506bcda92d6eb3c4f2c8006b60 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 52450b5d743302787e428e4afef617837532118d..7d19df39cb525f113f2e8bf5addc32741ddb54a6 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 31499ee087108db916d50b3eed4db93de1037acc..b063d080ac8349e1e76670d01aa808245fd46e4f 100644 --- a/dumux/material/fluidmatrixinteractions/frictionlaws/frictionlaw.hh +++ b/dumux/material/fluidmatrixinteractions/frictionlaws/frictionlaw.hh @@ -43,25 +43,28 @@ 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. * * 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 ----------- @@ -82,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; }; diff --git a/dumux/material/fluidmatrixinteractions/frictionlaws/manning.hh b/dumux/material/fluidmatrixinteractions/frictionlaws/manning.hh index a0b8a369ecbe6d83612103047c1bc0dc19e83e7b..51e70977ee04a5b840c7d4d623c32425ded6fb57 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 10a8d15cc46e6f8ab91767ebebad1c37bc4891f4..258341e71642e650c136b8cb7cbfe68ca6d865de 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 b2072fc4172354261d3b2e16e56272c3c06fe126..851a0aa9252a92a8a9fddfa77797269a63c1d3de 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 aff44d7233fae9deeacc2026ce6d673ba3d00f55..eee68f25d2ace4dfaef371f250b0b9f4b4753116 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 a43f11fe0b9400e7d9cdd4166b4cb18b858abb64..837b7d83d28dc1b0ed5b53ea58404d5bcc4c995b 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; }