| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- | ||
| 2 | // vi: set et ts=4 sw=4 sts=4: | ||
| 3 | // | ||
| 4 | // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder | ||
| 5 | // SPDX-License-Identifier: GPL-3.0-or-later | ||
| 6 | // | ||
| 7 | /*! | ||
| 8 | * \file | ||
| 9 | * \ingroup Flux | ||
| 10 | * \copydoc Dumux::ShallowWaterViscousFlux | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_FLUX_SHALLOW_WATER_VISCOUS_FLUX_HH | ||
| 13 | #define DUMUX_FLUX_SHALLOW_WATER_VISCOUS_FLUX_HH | ||
| 14 | |||
| 15 | #include <cmath> | ||
| 16 | #include <algorithm> | ||
| 17 | #include <utility> | ||
| 18 | #include <type_traits> | ||
| 19 | #include <array> | ||
| 20 | |||
| 21 | #include <dune/common/std/type_traits.hh> | ||
| 22 | #include <dune/common/exceptions.hh> | ||
| 23 | |||
| 24 | #include <dumux/common/parameters.hh> | ||
| 25 | #include <dumux/flux/fluxvariablescaching.hh> | ||
| 26 | #include <dumux/flux/shallowwater/fluxlimiterlet.hh> | ||
| 27 | |||
| 28 | namespace Dumux { | ||
| 29 | |||
| 30 | #ifndef DOXYGEN | ||
| 31 | namespace Detail { | ||
| 32 | // helper struct detecting if the user-defined spatial params class has a frictionLaw function | ||
| 33 | template <typename T, typename ...Ts> | ||
| 34 | using FrictionLawDetector = decltype(std::declval<T>().frictionLaw(std::declval<Ts>()...)); | ||
| 35 | |||
| 36 | template<class T, typename ...Args> | ||
| 37 | static constexpr bool implementsFrictionLaw() | ||
| 38 | { return Dune::Std::is_detected<FrictionLawDetector, T, Args...>::value; } | ||
| 39 | } // end namespace Detail | ||
| 40 | #endif | ||
| 41 | |||
| 42 | /*! | ||
| 43 | * \ingroup Flux | ||
| 44 | * \brief Compute the shallow water viscous momentum flux due to viscosity. | ||
| 45 | * | ||
| 46 | * The viscous momentum flux | ||
| 47 | * \f[ | ||
| 48 | * \int \int_{V} \mathbf{\nabla} \cdot \nu_t h \mathbf{\nabla} \mathbf{u} dV | ||
| 49 | * \f] | ||
| 50 | * is re-written using Gauss' divergence theorem to: | ||
| 51 | * \f[ | ||
| 52 | * \int_{S_f} \nu_t h \mathbf{\nabla} \mathbf{u} \cdot \mathbf{n_f} dS | ||
| 53 | * \f] | ||
| 54 | * | ||
| 55 | * The effective kinematic viscosity \f$ \nu_t \f$ | ||
| 56 | * can be calculated by adding a vertical (Elder-like) | ||
| 57 | * and a horizontal (Smagorinsky-like) part. This enabled by setting | ||
| 58 | * "ShallowWater.UseMixingLengthTurbulenceModel" to "true". | ||
| 59 | * | ||
| 60 | * For now the calculation of the shallow water viscous momentum flux is implemented | ||
| 61 | * strictly for 2D depth-averaged models (i.e. 3 equations). | ||
| 62 | */ | ||
| 63 | template<class NumEqVector, typename std::enable_if_t<NumEqVector::size() == 3, int> = 0> | ||
| 64 | class ShallowWaterViscousFlux | ||
| 65 | { | ||
| 66 | |||
| 67 | public: | ||
| 68 | |||
| 69 | using Cache = FluxVariablesCaching::EmptyDiffusionCache; | ||
| 70 | using CacheFiller = FluxVariablesCaching::EmptyCacheFiller; | ||
| 71 | /*! | ||
| 72 | * \ingroup Flux | ||
| 73 | * \brief Compute the viscous momentum flux contribution from the interface | ||
| 74 | * shear stress | ||
| 75 | */ | ||
| 76 | template<class Problem, class FVElementGeometry, class ElementVolumeVariables> | ||
| 77 | 4583245 | static NumEqVector flux(const Problem& problem, | |
| 78 | const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element, | ||
| 79 | const FVElementGeometry& fvGeometry, | ||
| 80 | const ElementVolumeVariables& elemVolVars, | ||
| 81 | const typename FVElementGeometry::SubControlVolumeFace& scvf) | ||
| 82 | { | ||
| 83 | using Scalar = typename NumEqVector::value_type; | ||
| 84 | using FluidSystem = typename ElementVolumeVariables::VolumeVariables::FluidSystem; | ||
| 85 | |||
| 86 | 4583245 | NumEqVector localFlux(0.0); | |
| 87 | |||
| 88 | // Get the inside and outside volume variables | ||
| 89 | 9166490 | const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()]; | |
| 90 | 9166490 | const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()]; | |
| 91 | |||
| 92 | 13749735 | const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); | |
| 93 | 9166490 | const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx()); | |
| 94 | |||
| 95 | 4583245 | const auto gradVelocity = [&]() | |
| 96 | { | ||
| 97 | // The left (inside) and right (outside) states | ||
| 98 | 4583245 | const auto velocityXLeft = insideVolVars.velocity(0); | |
| 99 | 4583245 | const auto velocityYLeft = insideVolVars.velocity(1); | |
| 100 | 4583245 | const auto velocityXRight = outsideVolVars.velocity(0); | |
| 101 | 4583245 | const auto velocityYRight = outsideVolVars.velocity(1); | |
| 102 | |||
| 103 | // Compute the velocity gradients normal to the face | ||
| 104 | // Factor that takes the direction of the unit vector into account | ||
| 105 | 13749735 | const auto& cellCenterToCellCenter = outsideScv.center() - insideScv.center(); | |
| 106 | 4583245 | const auto distance = cellCenterToCellCenter.two_norm(); | |
| 107 | 4583245 | const auto& unitNormal = scvf.unitOuterNormal(); | |
| 108 | 4583245 | const auto direction = (unitNormal*cellCenterToCellCenter)/distance; | |
| 109 | return std::array<Scalar, 2>{ | ||
| 110 | 4583245 | (velocityXRight-velocityXLeft)*direction/distance, | |
| 111 | 4583245 | (velocityYRight-velocityYLeft)*direction/distance | |
| 112 | 9166490 | }; | |
| 113 | 4583245 | }(); | |
| 114 | |||
| 115 | // Use a harmonic average of the depth at the interface. | ||
| 116 | 4583245 | const auto waterDepthLeft = insideVolVars.waterDepth(); | |
| 117 | 4583245 | const auto waterDepthRight = outsideVolVars.waterDepth(); | |
| 118 | 4583245 | const auto averageDepth = 2.0*(waterDepthLeft*waterDepthRight)/(waterDepthLeft + waterDepthRight); | |
| 119 | |||
| 120 | 4583245 | const Scalar effectiveKinematicViscosity = [&]() | |
| 121 | { | ||
| 122 | // The background viscosity, per default this is just the fluid viscosity (e.g. for the viscous flow regime) | ||
| 123 | // The default may be overwritten by the user, by setting the parameter "ShallowWater.TurbulentViscosity" to | ||
| 124 | // enable setting a background viscosity that already contains some effects of turbulence | ||
| 125 |
4/6✓ Branch 0 taken 7 times.
✓ Branch 1 taken 4583238 times.
✓ Branch 3 taken 7 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 7 times.
✗ Branch 7 not taken.
|
4583245 | static const auto backgroundKinematicViscosity = getParamFromGroup<Scalar>( |
| 126 |
3/6✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 7 times.
✗ Branch 8 not taken.
|
28 | problem.paramGroup(), "ShallowWater.TurbulentViscosity", insideVolVars.viscosity()/insideVolVars.density() |
| 127 | ); | ||
| 128 | |||
| 129 | // Check whether the mixing-length turbulence model is used | ||
| 130 |
4/6✓ Branch 0 taken 7 times.
✓ Branch 1 taken 4583238 times.
✓ Branch 3 taken 7 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 7 times.
✗ Branch 7 not taken.
|
4583245 | static const auto useMixingLengthTurbulenceModel = getParamFromGroup<bool>( |
| 131 |
1/2✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
|
14 | problem.paramGroup(), "ShallowWater.UseMixingLengthTurbulenceModel", false |
| 132 | ); | ||
| 133 | |||
| 134 | // constant eddy viscosity equal to the prescribed background eddy viscosity | ||
| 135 |
1/2✓ Branch 0 taken 4583245 times.
✗ Branch 1 not taken.
|
4583245 | if (!useMixingLengthTurbulenceModel) |
| 136 | 4583245 | return backgroundKinematicViscosity; | |
| 137 | |||
| 138 | using SpatialParams = typename Problem::SpatialParams; | ||
| 139 | using E = typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity; | ||
| 140 | using SCV = typename FVElementGeometry::SubControlVolume; | ||
| 141 | if constexpr (!Detail::implementsFrictionLaw<SpatialParams, E, SCV>()) | ||
| 142 | ✗ | DUNE_THROW(Dune::IOError, "Mixing length turbulence model enabled but spatial parameters do not implement the frictionLaw interface!"); | |
| 143 | else | ||
| 144 | { | ||
| 145 | // turbulence model based on mixing length | ||
| 146 | // Compute the turbulent viscosity using a combined horizontal/vertical mixing length approach | ||
| 147 | // Turbulence coefficients: vertical (Elder like) and horizontal (Smagorinsky like) | ||
| 148 | ✗ | static const auto turbConstV = getParamFromGroup<Scalar>(problem.paramGroup(), "ShallowWater.VerticalCoefficientOfMixingLengthModel", 1.0); | |
| 149 | ✗ | static const auto turbConstH = getParamFromGroup<Scalar>(problem.paramGroup(), "ShallowWater.HorizontalCoefficientOfMixingLengthModel", 0.1); | |
| 150 | |||
| 151 | /** The vertical (Elder-like) contribution to the turbulent viscosity scales with water depth \f[ h \f] and shear velocity \f[ u_{*} \f] : | ||
| 152 | * | ||
| 153 | * \f[ | ||
| 154 | * \nu_t^v = c^v \frac{\kappa}{6} u_{*} h | ||
| 155 | * \f] | ||
| 156 | */ | ||
| 157 | ✗ | constexpr Scalar kappa = 0.41; | |
| 158 | // Compute the average shear velocity on the face | ||
| 159 | ✗ | const Scalar ustar = [&]() | |
| 160 | { | ||
| 161 | // Get the bottom shear stress in the two adjacent cells | ||
| 162 | // Note that element is not needed in spatialParams().frictionLaw (should be removed). For now we simply pass the same element twice | ||
| 163 | ✗ | const auto& bottomShearStressInside = problem.spatialParams().frictionLaw(element, insideScv).bottomShearStress(insideVolVars); | |
| 164 | ✗ | const auto& bottomShearStressOutside = problem.spatialParams().frictionLaw(element, outsideScv).bottomShearStress(outsideVolVars); | |
| 165 | ✗ | const auto bottomShearStressInsideMag = bottomShearStressInside.two_norm(); | |
| 166 | ✗ | const auto bottomShearStressOutsideMag = bottomShearStressOutside.two_norm(); | |
| 167 | |||
| 168 | // Use a harmonic average of the viscosity and the depth at the interface. | ||
| 169 | using std::max; | ||
| 170 | ✗ | const auto averageBottomShearStress = 2.0*(bottomShearStressInsideMag*bottomShearStressOutsideMag) | |
| 171 | ✗ | / max(1.0e-8,bottomShearStressInsideMag + bottomShearStressOutsideMag); | |
| 172 | |||
| 173 | // Shear velocity possibly needed in mixing-length turbulence model in the computation of the diffusion flux | ||
| 174 | using std::sqrt; | ||
| 175 | static_assert(!FluidSystem::isCompressible(0), | ||
| 176 | "The shallow water model assumes incompressible fluids" | ||
| 177 | ); | ||
| 178 | // Since the shallow water equations are incompressible, the density is constant. | ||
| 179 | // Therefore, insideVolVars.density() is equal to outsideVolVars.density(). | ||
| 180 | ✗ | return sqrt(averageBottomShearStress/insideVolVars.density()); | |
| 181 | ✗ | }(); | |
| 182 | |||
| 183 | ✗ | const auto turbViscosityV = turbConstV * (kappa/6.0) * ustar * averageDepth; | |
| 184 | |||
| 185 | /** The horizontal (Smagorinsky-like) contribution to the turbulent viscosity scales with the water depth (squared) | ||
| 186 | * and the magnitude of the stress (rate-of-strain) tensor: | ||
| 187 | * | ||
| 188 | * \f[ | ||
| 189 | * \nu_t^h = (c^h h)^2 \sqrt{ 2\left(\frac{\partial u}{\partial x}\right)^2 + \left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}\right)^2 + 2\left(\frac{\partial v}{\partial y}\right)^2 } | ||
| 190 | * \f] | ||
| 191 | * | ||
| 192 | * However, based on the velocity vectors in the direct neighbours of the volume face, it is not possible to compute all components of the stress tensor. | ||
| 193 | * Only the gradients of u and v in the direction of the vector between the cell centres is available. | ||
| 194 | * To avoid the reconstruction of the full velocity gradient tensor based on a larger stencil, | ||
| 195 | * the horizontal contribution to the eddy viscosity (in the mixing-length model) is computed using only the velocity gradients normal to the face: | ||
| 196 | * | ||
| 197 | * \f[ | ||
| 198 | * \frac{\partial u}{\partial n} , \frac{\partial v}{\partial n} | ||
| 199 | * \f] | ||
| 200 | * | ||
| 201 | * In other words, the present approximation of the horizontal contribution to the turbulent viscosity reduces to: | ||
| 202 | * | ||
| 203 | * \f[ | ||
| 204 | * \nu_t^h = (c^h h)^2 \sqrt{ 2\left(\frac{\partial u}{\partial n}\right)^2 + 2\left(\frac{\partial v}{\partial n}\right)^2 } | ||
| 205 | * \f] | ||
| 206 | * | ||
| 207 | It should be noted that this simplified approach is formally inconsistent and will result in a turbulent viscosity that is dependent on the grid (orientation). | ||
| 208 | */ | ||
| 209 | using std::sqrt; | ||
| 210 | ✗ | const auto [gradU, gradV] = gradVelocity; | |
| 211 | ✗ | const auto mixingLengthSquared = turbConstH * turbConstH * averageDepth * averageDepth; | |
| 212 | ✗ | const auto turbViscosityH = mixingLengthSquared * sqrt(2.0*gradU*gradU + 2.0*gradV*gradV); | |
| 213 | |||
| 214 | // Total turbulent viscosity | ||
| 215 | ✗ | return backgroundKinematicViscosity + sqrt(turbViscosityV*turbViscosityV + turbViscosityH*turbViscosityH); | |
| 216 | } | ||
| 217 | 4583245 | }(); | |
| 218 | |||
| 219 | // Compute the viscous momentum fluxes with connecting water depth at the interface | ||
| 220 | using std::min; | ||
| 221 | using std::max; | ||
| 222 |
2/2✓ Branch 0 taken 2281176 times.
✓ Branch 1 taken 2302069 times.
|
4583245 | const auto freeSurfaceInside = insideVolVars.waterDepth() + insideVolVars.bedSurface(); |
| 223 |
2/2✓ Branch 0 taken 2281176 times.
✓ Branch 1 taken 2302069 times.
|
4583245 | const auto freeSurfaceOutside = outsideVolVars.waterDepth() + outsideVolVars.bedSurface(); |
| 224 |
8/10✓ Branch 0 taken 2281176 times.
✓ Branch 1 taken 2302069 times.
✓ Branch 2 taken 1743628 times.
✓ Branch 3 taken 2839617 times.
✓ Branch 4 taken 1743628 times.
✓ Branch 5 taken 2839617 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 4583245 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 4583245 times.
|
12632853 | const auto interfaceWaterDepth = max(min(freeSurfaceInside , freeSurfaceOutside) - max(insideVolVars.bedSurface(),outsideVolVars.bedSurface()),0.0); |
| 225 | 13749735 | const auto& [gradU, gradV] = gradVelocity; | |
| 226 | 4583245 | const auto uViscousFlux = effectiveKinematicViscosity * interfaceWaterDepth * gradU; | |
| 227 | 4583245 | const auto vViscousFlux = effectiveKinematicViscosity * interfaceWaterDepth * gradV; | |
| 228 | |||
| 229 | 9166490 | localFlux[0] = 0.0; | |
| 230 | 9166490 | localFlux[1] = -uViscousFlux * scvf.area(); | |
| 231 | 9166490 | localFlux[2] = -vViscousFlux * scvf.area(); | |
| 232 | |||
| 233 | 4583245 | return localFlux; | |
| 234 | } | ||
| 235 | }; | ||
| 236 | |||
| 237 | } // end namespace Dumux | ||
| 238 | |||
| 239 | #endif | ||
| 240 |