| 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 | /*! | ||
| 9 | * \file | ||
| 10 | * \ingroup NavierStokesModel | ||
| 11 | * \copydoc Dumux::NavierStokesMomentumBoundaryFlux | ||
| 12 | */ | ||
| 13 | #ifndef DUMUX_FREEFLOW_NAVIERSTOKES_MOMENTUM_BOUNDARY_FLUXHELPER_HH | ||
| 14 | #define DUMUX_FREEFLOW_NAVIERSTOKES_MOMENTUM_BOUNDARY_FLUXHELPER_HH | ||
| 15 | |||
| 16 | #include <dumux/common/math.hh> | ||
| 17 | #include <dumux/common/parameters.hh> | ||
| 18 | #include <dumux/discretization/method.hh> | ||
| 19 | #include <dumux/freeflow/navierstokes/momentum/velocitygradients.hh> | ||
| 20 | #include <dumux/freeflow/navierstokes/slipcondition.hh> | ||
| 21 | |||
| 22 | namespace Dumux { | ||
| 23 | |||
| 24 | /*! | ||
| 25 | * \ingroup NavierStokesModel | ||
| 26 | * \brief Class to compute the boundary flux for the momentum balance of the Navier-Stokes model | ||
| 27 | * This helper class is typically used in the Neumann function of the momentum problem. | ||
| 28 | * | ||
| 29 | * \tparam DiscretizationMethod The discretization method used for the momentum problem | ||
| 30 | * \tparam SlipVelocityPolicy The policy to compute the slip velocity at the boundary | ||
| 31 | * | ||
| 32 | * The slip velocity policy is a class with a static member function `velocity` | ||
| 33 | * that computes the slip velocity at the boundary. | ||
| 34 | */ | ||
| 35 | template<class DiscretizationMethod, class SlipVelocityPolicy = void> | ||
| 36 | struct NavierStokesMomentumBoundaryFlux; | ||
| 37 | |||
| 38 | /*! | ||
| 39 | * \ingroup NavierStokesModel | ||
| 40 | * \brief Class to compute the boundary flux for the momentum balance of the Navier-Stokes model | ||
| 41 | * This helper class is typically used in the Neumann function of the momentum problem. | ||
| 42 | * | ||
| 43 | * \tparam SlipVelocityPolicy The policy to compute the slip velocity at the boundary | ||
| 44 | * | ||
| 45 | * The slip velocity policy is a class with a static member function `velocity` | ||
| 46 | * that computes the slip velocity at the boundary. | ||
| 47 | */ | ||
| 48 | template<class SlipVelocityPolicy> | ||
| 49 | struct NavierStokesMomentumBoundaryFlux<DiscretizationMethods::FCStaggered, SlipVelocityPolicy> | ||
| 50 | { | ||
| 51 | /*! | ||
| 52 | * \brief Returns the momentum flux a fixed-pressure boundary. | ||
| 53 | * \param problem The problem | ||
| 54 | * \param fvGeometry The finite-volume geometry | ||
| 55 | * \param scvf The sub control volume face | ||
| 56 | * \param elemVolVars The volume variables for the element | ||
| 57 | * \param elemFluxVarsCache The flux variables cache for the element | ||
| 58 | * \param pressure The pressure given at the boundary | ||
| 59 | * \param zeroNormalVelocityGradient If set to false, this yields a "natural" outflow condition where the | ||
| 60 | * flow tends to diverge at the boundary. If set to true (default), this | ||
| 61 | * yields a "do-nothing" outflow condition where \f$\mathbf{v} \cdot \mathbf{n} = 0\f$ | ||
| 62 | * is satisfied and the general flow field is preserved in many cases. | ||
| 63 | * See, e.g., https://www.math.uni-magdeburg.de/~richter/WS17/numns/fem.pdf for details. | ||
| 64 | */ | ||
| 65 | template<class Problem, class FVElementGeometry, class ElementVolumeVariables, class ElementFluxVariablesCache, class Scalar> | ||
| 66 | 721439 | static auto fixedPressureMomentumFlux(const Problem& problem, | |
| 67 | const FVElementGeometry& fvGeometry, | ||
| 68 | const typename FVElementGeometry::SubControlVolumeFace& scvf, | ||
| 69 | const ElementVolumeVariables& elemVolVars, | ||
| 70 | const ElementFluxVariablesCache& elemFluxVarsCache, | ||
| 71 | const Scalar pressure, | ||
| 72 | const bool zeroNormalVelocityGradient = true) | ||
| 73 | { | ||
| 74 | // TODO density upwinding? | ||
| 75 | using MomentumFlux = typename Problem::MomentumFluxType; | ||
| 76 | 721439 | constexpr std::size_t dim = static_cast<std::size_t>(FVElementGeometry::GridGeometry::GridView::dimension); | |
| 77 | static_assert( | ||
| 78 | MomentumFlux::size() == dim, | ||
| 79 | "Expects momentum flux to have as many entries as dimensions." | ||
| 80 | ); | ||
| 81 | |||
| 82 | 721439 | const auto& element = fvGeometry.element(); | |
| 83 | |||
| 84 | static_assert( | ||
| 85 | std::decay_t<decltype(problem.dirichlet(element, scvf))>::size() == dim, | ||
| 86 | "Expects problem.dirichlet to return an array with as many entries as dimensions." | ||
| 87 | ); | ||
| 88 | |||
| 89 |
2/2✓ Branch 0 taken 116055 times.
✓ Branch 1 taken 605384 times.
|
721439 | if (scvf.isFrontal()) |
| 90 | { | ||
| 91 | 116055 | MomentumFlux flux(0.0); | |
| 92 | |||
| 93 | // pressure contribution | ||
| 94 |
2/2✓ Branch 0 taken 85740 times.
✓ Branch 1 taken 30315 times.
|
116055 | flux[scvf.normalAxis()] = (pressure - problem.referencePressure(element, fvGeometry, scvf)) * scvf.directionSign(); |
| 95 | |||
| 96 | // advective terms | ||
| 97 |
2/2✓ Branch 0 taken 85740 times.
✓ Branch 1 taken 30315 times.
|
116055 | if (problem.enableInertiaTerms()) |
| 98 | { | ||
| 99 | 171480 | const auto v = elemVolVars[scvf.insideScvIdx()].velocity(); | |
| 100 | 85740 | flux[scvf.normalAxis()] += v*v * problem.density(element, fvGeometry, scvf) * scvf.directionSign(); | |
| 101 | } | ||
| 102 | |||
| 103 | 116055 | return flux; | |
| 104 | } | ||
| 105 |
1/2✓ Branch 0 taken 605384 times.
✗ Branch 1 not taken.
|
605384 | else if (scvf.isLateral()) |
| 106 | { | ||
| 107 | // if no zero velocity gradient is desired the lateral flux is zero | ||
| 108 |
2/2✓ Branch 0 taken 5556 times.
✓ Branch 1 taken 599828 times.
|
605384 | if (!zeroNormalVelocityGradient) |
| 109 | 11112 | return MomentumFlux(0.0); | |
| 110 | |||
| 111 | // otherwise, make sure the flow does not diverge by accounting for the off-diagonal entries of the stress tensor | ||
| 112 | |||
| 113 | // if the lateral scvf is adjacent to a slip boundary, call the helper function for this special case | ||
| 114 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 567796 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 567796 times.
|
1199656 | const auto& scv = fvGeometry.scv(scvf.insideScvIdx()); |
| 115 |
4/4✓ Branch 0 taken 298490 times.
✓ Branch 1 taken 301338 times.
✓ Branch 4 taken 80 times.
✓ Branch 5 taken 2614 times.
|
599828 | if (scv.boundary() && problem.onSlipBoundary(fvGeometry, fvGeometry.frontalScvfOnBoundary(scv))) |
| 116 | 80 | return slipVelocityMomentumFlux(problem, fvGeometry, scvf, elemVolVars, elemFluxVarsCache); | |
| 117 | |||
| 118 | // viscous terms | ||
| 119 | 599748 | MomentumFlux flux(0.0); | |
| 120 | 599748 | const Scalar mu = problem.effectiveViscosity(element, fvGeometry, scvf); | |
| 121 | |||
| 122 | // lateral face normal to boundary (integration point touches boundary) | ||
| 123 |
2/2✓ Branch 0 taken 298410 times.
✓ Branch 1 taken 301338 times.
|
599748 | if (scv.boundary()) |
| 124 | 298410 | flux[scv.dofAxis()] -= mu * StaggeredVelocityGradients::velocityGradIJ(fvGeometry, scvf, elemVolVars) | |
| 125 | 298410 | * scvf.directionSign(); | |
| 126 | // lateral face coinciding with boundary | ||
| 127 |
1/2✓ Branch 0 taken 301338 times.
✗ Branch 1 not taken.
|
301338 | else if (scvf.boundary()) |
| 128 | 301338 | flux[scv.dofAxis()] -= mu * StaggeredVelocityGradients::velocityGradJI(fvGeometry, scvf, elemVolVars) | |
| 129 | 301338 | * scvf.directionSign(); | |
| 130 | |||
| 131 | // advective terms | ||
| 132 |
2/2✓ Branch 0 taken 493960 times.
✓ Branch 1 taken 105788 times.
|
599748 | if (problem.enableInertiaTerms()) // TODO revise advection terms! Take care with upwinding! |
| 133 | { | ||
| 134 | 987920 | const auto transportingVelocity = [&]() | |
| 135 | { | ||
| 136 | 493960 | const auto& orthogonalScvf = fvGeometry.lateralOrthogonalScvf(scvf); | |
| 137 |
2/2✓ Branch 2 taken 10 times.
✓ Branch 3 taken 493950 times.
|
987920 | const auto innerTransportingVelocity = elemVolVars[orthogonalScvf.insideScvIdx()].velocity(); |
| 138 | |||
| 139 |
4/6✓ Branch 0 taken 10 times.
✓ Branch 1 taken 493950 times.
✓ Branch 3 taken 10 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 10 times.
✗ Branch 7 not taken.
|
493960 | static const bool useOldScheme = getParam<bool>("FreeFlow.UseOldTransportingVelocity", true); // TODO how to deprecate? |
| 140 | |||
| 141 |
1/2✓ Branch 0 taken 493960 times.
✗ Branch 1 not taken.
|
493960 | if (useOldScheme) |
| 142 | { | ||
| 143 |
2/2✓ Branch 0 taken 267456 times.
✓ Branch 1 taken 226504 times.
|
493960 | if (scvf.boundary()) |
| 144 | { | ||
| 145 |
2/2✓ Branch 1 taken 17648 times.
✓ Branch 2 taken 249808 times.
|
517264 | if (problem.boundaryTypes(element, scvf).isDirichlet(scvf.normalAxis())) |
| 146 | 45808 | return problem.dirichlet(element, scvf)[scvf.normalAxis()]; | |
| 147 | } | ||
| 148 | else | ||
| 149 | return innerTransportingVelocity; | ||
| 150 | } | ||
| 151 | |||
| 152 | // use the Dirichlet velocity as transporting velocity if the lateral face is on a Dirichlet boundary | ||
| 153 |
1/2✓ Branch 0 taken 249808 times.
✗ Branch 1 not taken.
|
249808 | if (scvf.boundary()) |
| 154 | { | ||
| 155 |
1/2✓ Branch 1 taken 249808 times.
✗ Branch 2 not taken.
|
499616 | if (problem.boundaryTypes(element, scvf).isDirichlet(scvf.normalAxis())) |
| 156 | ✗ | return 0.5*(problem.dirichlet(element, scvf)[scvf.normalAxis()] + innerTransportingVelocity); | |
| 157 | } | ||
| 158 | |||
| 159 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 249808 times.
|
249808 | if (orthogonalScvf.boundary()) |
| 160 | { | ||
| 161 | ✗ | if (problem.boundaryTypes(element, orthogonalScvf).isDirichlet(scvf.normalAxis())) | |
| 162 | ✗ | return 0.5*(problem.dirichlet(element, scvf)[scvf.normalAxis()] + innerTransportingVelocity); | |
| 163 | else | ||
| 164 | ✗ | return innerTransportingVelocity; // fallback value, should actually never be called | |
| 165 | } | ||
| 166 | |||
| 167 | // average the transporting velocity by weighting with the scv volumes | ||
| 168 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 235024 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 235024 times.
|
499616 | const auto insideVolume = fvGeometry.scv(orthogonalScvf.insideScvIdx()).volume(); |
| 169 | 499616 | const auto outsideVolume = fvGeometry.scv(orthogonalScvf.outsideScvIdx()).volume(); | |
| 170 | 499616 | const auto outerTransportingVelocity = elemVolVars[orthogonalScvf.outsideScvIdx()].velocity(); | |
| 171 | 249808 | return (insideVolume*innerTransportingVelocity + outsideVolume*outerTransportingVelocity) / (insideVolume + outsideVolume); | |
| 172 | 493960 | }(); | |
| 173 | |||
| 174 | // lateral face normal to boundary (integration point touches boundary) | ||
| 175 |
2/2✓ Branch 0 taken 244152 times.
✓ Branch 1 taken 249808 times.
|
493960 | if (scv.boundary()) |
| 176 | { | ||
| 177 | 488304 | const auto innerVelocity = elemVolVars[scvf.insideScvIdx()].velocity(); | |
| 178 | 488304 | const auto outerVelocity = elemVolVars[scvf.outsideScvIdx()].velocity(); | |
| 179 | |||
| 180 | 244152 | const auto rho = problem.insideAndOutsideDensity(element, fvGeometry, scvf); | |
| 181 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 244142 times.
|
244152 | const bool selfIsUpstream = scvf.directionSign() == sign(transportingVelocity); |
| 182 | |||
| 183 | 244152 | const auto insideMomentum = innerVelocity * rho.first; | |
| 184 | 244152 | const auto outsideMomentum = outerVelocity * rho.second; | |
| 185 | |||
| 186 |
5/8✓ Branch 0 taken 10 times.
✓ Branch 1 taken 244142 times.
✓ Branch 3 taken 10 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 10 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 10 times.
✗ Branch 10 not taken.
|
244152 | static const auto upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(), "Flux.UpwindWeight"); |
| 187 | |||
| 188 |
2/2✓ Branch 0 taken 97155 times.
✓ Branch 1 taken 146997 times.
|
244152 | const auto transportedMomentum = selfIsUpstream ? (upwindWeight * insideMomentum + (1.0 - upwindWeight) * outsideMomentum) |
| 189 | 146997 | : (upwindWeight * outsideMomentum + (1.0 - upwindWeight) * insideMomentum); | |
| 190 | |||
| 191 | 488304 | flux[scv.dofAxis()] += transportingVelocity * transportedMomentum * scvf.directionSign(); | |
| 192 | } | ||
| 193 | |||
| 194 | // lateral face coinciding with boundary | ||
| 195 |
1/2✓ Branch 0 taken 249808 times.
✗ Branch 1 not taken.
|
249808 | else if (scvf.boundary()) |
| 196 | { | ||
| 197 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 235024 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 235024 times.
|
499616 | const auto insideDensity = problem.density(element, fvGeometry.scv(scvf.insideScvIdx())); |
| 198 | 499616 | const auto innerVelocity = elemVolVars[scvf.insideScvIdx()].velocity(); | |
| 199 | 499616 | flux[scv.dofAxis()] += innerVelocity * transportingVelocity * insideDensity * scvf.directionSign(); | |
| 200 | } | ||
| 201 | } | ||
| 202 | |||
| 203 | 599748 | return flux; | |
| 204 | } | ||
| 205 | |||
| 206 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Face is neither lateral nor frontal."); | |
| 207 | } | ||
| 208 | |||
| 209 | /*! | ||
| 210 | * \brief Returns the momentum flux for a boundary with a slip condition. | ||
| 211 | * \param problem The problem | ||
| 212 | * \param fvGeometry The finite-volume geometry | ||
| 213 | * \param scvf The sub control volume face | ||
| 214 | * \param elemVolVars The volume variables for the element | ||
| 215 | * \param elemFluxVarsCache The flux variables cache for the element | ||
| 216 | * \param tangentialVelocityGradient A user-specified velocity gradient in tangential direction of the boundary. Only used for certain situations | ||
| 217 | * where this gradient cannot be determined automatically (e.g., at lateral Neumann boundaries). | ||
| 218 | */ | ||
| 219 | template<class Problem, class FVElementGeometry, class ElementVolumeVariables, class ElementFluxVariablesCache, class TangentialVelocityGradient = double> | ||
| 220 | 182881 | static auto slipVelocityMomentumFlux(const Problem& problem, | |
| 221 | const FVElementGeometry& fvGeometry, | ||
| 222 | const typename FVElementGeometry::SubControlVolumeFace& scvf, | ||
| 223 | const ElementVolumeVariables& elemVolVars, | ||
| 224 | const ElementFluxVariablesCache& elemFluxVarsCache, | ||
| 225 | const TangentialVelocityGradient& tangentialVelocityGradient = TangentialVelocityGradient(0.0)) | ||
| 226 | { | ||
| 227 | using MomentumFlux = typename Problem::MomentumFluxType; | ||
| 228 | 182881 | constexpr std::size_t dim = static_cast<std::size_t>(FVElementGeometry::GridGeometry::GridView::dimension); | |
| 229 | static_assert( | ||
| 230 | MomentumFlux::size() == dim, | ||
| 231 | "Expects momentum flux to have as many entries as dimensions." | ||
| 232 | ); | ||
| 233 | |||
| 234 | 182881 | MomentumFlux flux(0.0); | |
| 235 | |||
| 236 | // slip velocity momentum contribution only makes sense for lateral scvfs | ||
| 237 |
2/2✓ Branch 0 taken 35571 times.
✓ Branch 1 taken 147310 times.
|
182881 | if (scvf.isFrontal()) |
| 238 | 35571 | return flux; | |
| 239 | |||
| 240 | using Scalar = std::decay_t<decltype(flux[0])>; | ||
| 241 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 147310 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 147310 times.
|
294620 | const auto& scv = fvGeometry.scv(scvf.insideScvIdx()); |
| 242 | 147310 | const auto& orthogonalScvf = fvGeometry.lateralOrthogonalScvf(scvf); | |
| 243 |
4/4✓ Branch 0 taken 76908 times.
✓ Branch 1 taken 70402 times.
✓ Branch 2 taken 76908 times.
✓ Branch 3 taken 70402 times.
|
294620 | const auto& orthogonalScv = fvGeometry.scv(orthogonalScvf.insideScvIdx()); |
| 244 | |||
| 245 |
5/5✓ Branch 0 taken 76908 times.
✓ Branch 1 taken 70402 times.
✓ Branch 3 taken 73128 times.
✓ Branch 4 taken 3700 times.
✓ Branch 5 taken 80 times.
|
147310 | if (scvf.boundary() && problem.onSlipBoundary(fvGeometry, fvGeometry.frontalScvfOnBoundary(orthogonalScv))) |
| 246 | { | ||
| 247 | /* | ||
| 248 | * ---------#######::::::::: x dof position ***** porous boundary at bottom | ||
| 249 | * slip | || | :: | ||
| 250 | * gradient | || | velI :: -- element | ||
| 251 | * -------> | || scv x~~~~> :: | ||
| 252 | * ------> | || | :: O position at which gradient is evaluated (integration point) | ||
| 253 | * -----> | velJ ^ | :^ | ||
| 254 | * ----> -------|-######O::::::::| <----This velocity dof (outer velJ) does not exist if the scv itself lies on a | ||
| 255 | * ***************~~~>****** non-Dirichlet boundary. In that case, use the given tangentialVelocityGradient. | ||
| 256 | * frontal scvf v_slip | ||
| 257 | * on porous || and # staggered half-control-volume (own element) | ||
| 258 | * boundary | ||
| 259 | * :: staggered half-control-volume (neighbor element) | ||
| 260 | * | ||
| 261 | */ | ||
| 262 | |||
| 263 | if constexpr (!std::is_void_v<SlipVelocityPolicy>) | ||
| 264 | { | ||
| 265 | 153656 | const Scalar velI = elemVolVars[scvf.insideScvIdx()].velocity(); | |
| 266 | |||
| 267 | // viscous terms | ||
| 268 | 76828 | const Scalar mu = problem.effectiveViscosity(fvGeometry.element(), fvGeometry, scvf); | |
| 269 | 307312 | const Scalar distance = (scv.dofPosition()- scvf.ipGlobal()).two_norm(); | |
| 270 | 76828 | const Scalar velGradJI = [&] | |
| 271 | { | ||
| 272 |
2/2✓ Branch 2 taken 76788 times.
✓ Branch 3 taken 40 times.
|
153656 | if (elemVolVars.hasVolVars(orthogonalScvf.outsideScvIdx())) |
| 273 | 76788 | return StaggeredVelocityGradients::velocityGradJI(fvGeometry, scvf, elemVolVars); | |
| 274 | else | ||
| 275 | 40 | return tangentialVelocityGradient; | |
| 276 | 76828 | }(); | |
| 277 | |||
| 278 |
2/2✓ Branch 1 taken 33240 times.
✓ Branch 2 taken 43588 times.
|
76828 | const Scalar slipVelocity = SlipVelocityPolicy::velocity(problem, fvGeometry, scvf, elemVolVars, velGradJI)[scv.dofAxis()]; |
| 279 | 76828 | const Scalar velGradIJ = (slipVelocity - velI) / distance * scvf.directionSign(); | |
| 280 | |||
| 281 |
2/2✓ Branch 0 taken 33240 times.
✓ Branch 1 taken 43588 times.
|
76828 | flux[scv.dofAxis()] -= (mu * (velGradIJ + velGradJI))*scvf.directionSign(); |
| 282 | |||
| 283 | // advective terms | ||
| 284 |
2/2✓ Branch 0 taken 33240 times.
✓ Branch 1 taken 43588 times.
|
76828 | if (problem.enableInertiaTerms()) |
| 285 | { | ||
| 286 | // transporting velocity corresponds to velJ | ||
| 287 | 33240 | const auto transportingVelocity = [&]() | |
| 288 | { | ||
| 289 | 66480 | const auto innerTransportingVelocity = elemVolVars[orthogonalScvf.insideScvIdx()].velocity(); | |
| 290 | |||
| 291 |
1/2✓ Branch 2 taken 33240 times.
✗ Branch 3 not taken.
|
66480 | if (!elemVolVars.hasVolVars(orthogonalScvf.outsideScvIdx())) |
| 292 | return innerTransportingVelocity; // fallback | ||
| 293 | |||
| 294 |
1/2✓ Branch 2 taken 33240 times.
✗ Branch 3 not taken.
|
66480 | const auto outerTransportingVelocity = elemVolVars[orthogonalScvf.outsideScvIdx()].velocity(); |
| 295 | |||
| 296 | // if the orthogonal scvf lies on a boundary and if there are outside volvars, we assume that these come from a Dirichlet condition | ||
| 297 |
1/2✓ Branch 0 taken 33240 times.
✗ Branch 1 not taken.
|
33240 | if (orthogonalScvf.boundary()) |
| 298 | return outerTransportingVelocity; | ||
| 299 | |||
| 300 |
4/6✓ Branch 0 taken 3 times.
✓ Branch 1 taken 33237 times.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
|
33240 | static const bool useOldScheme = getParam<bool>("FreeFlow.UseOldTransportingVelocity", true); // TODO how to deprecate? |
| 301 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 33240 times.
|
33240 | if (useOldScheme) |
| 302 | return innerTransportingVelocity; | ||
| 303 | else | ||
| 304 | { | ||
| 305 | // average the transporting velocity by weighting with the scv volumes | ||
| 306 | ✗ | const auto insideVolume = fvGeometry.scv(orthogonalScvf.insideScvIdx()).volume(); | |
| 307 | ✗ | const auto outsideVolume = fvGeometry.scv(orthogonalScvf.outsideScvIdx()).volume(); | |
| 308 | ✗ | const auto outerTransportingVelocity = elemVolVars[orthogonalScvf.outsideScvIdx()].velocity(); | |
| 309 | ✗ | return (insideVolume*innerTransportingVelocity + outsideVolume*outerTransportingVelocity) / (insideVolume + outsideVolume); | |
| 310 | } | ||
| 311 | 33240 | }(); | |
| 312 | |||
| 313 | // Do not use upwinding here but directly take the slip velocity located on the boundary. Upwinding with a weight of 0.5 | ||
| 314 | // would actually prevent second order grid convergence. | ||
| 315 | 33240 | const auto rho = problem.insideAndOutsideDensity(fvGeometry.element(), fvGeometry, scvf); | |
| 316 | 33240 | const auto transportedMomentum = slipVelocity * rho.second; | |
| 317 | |||
| 318 | 66480 | flux[scv.dofAxis()] += transportingVelocity * transportedMomentum * scvf.directionSign(); | |
| 319 | } | ||
| 320 | } | ||
| 321 | else | ||
| 322 | DUNE_THROW(Dune::InvalidStateException, "SlipVelocityPolicy needs to be specified as template argument"); | ||
| 323 | } | ||
| 324 |
3/5✓ Branch 0 taken 70482 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 66968 times.
✓ Branch 4 taken 3514 times.
✗ Branch 5 not taken.
|
70482 | else if (scv.boundary() && problem.onSlipBoundary(fvGeometry, fvGeometry.frontalScvfOnBoundary(scv))) |
| 325 | { | ||
| 326 | /* * | ||
| 327 | * ---------#######* x dof position ***** porous boundary at right side | ||
| 328 | * | || |* | ||
| 329 | * | || |* velI -- element | ||
| 330 | * | || scv x~~~~> | ||
| 331 | * | || |* O position at which gradient is evaluated (integration point) | ||
| 332 | * | velJ ^ |*^ | ||
| 333 | * -------|-######O*| v_slip || and # staggered half-control-volume (own element) | ||
| 334 | * | |* | ||
| 335 | * | neighbor |* | ||
| 336 | * | element ~~~~> velI outside (Does not exist if lower later lateral scvf lies on non-Dirichlet | ||
| 337 | * | |* boundary. In that case, use the given tangentialVelocityGradient.) | ||
| 338 | * | |* | ||
| 339 | * ----------------* :: staggered half-control-volume (neighbor element) | ||
| 340 | * ^ | ||
| 341 | * | ^ | ||
| 342 | * slip gradient | | ^ | ||
| 343 | * | | | | ||
| 344 | * | ||
| 345 | */ | ||
| 346 | |||
| 347 | if constexpr (!std::is_void_v<SlipVelocityPolicy>) | ||
| 348 | { | ||
| 349 | 140964 | const Scalar velJ = elemVolVars[orthogonalScvf.insideScvIdx()].velocity(); | |
| 350 | |||
| 351 | // viscous terms | ||
| 352 | 70482 | const Scalar mu = problem.effectiveViscosity(fvGeometry.element(), fvGeometry, scvf); | |
| 353 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 70482 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 70482 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 70482 times.
|
281928 | const Scalar distance = (fvGeometry.scv(orthogonalScvf.insideScvIdx()).dofPosition()- scvf.ipGlobal()).two_norm(); |
| 354 | |||
| 355 | 70482 | const Scalar velGradIJ = [&] | |
| 356 | { | ||
| 357 |
2/2✓ Branch 2 taken 70442 times.
✓ Branch 3 taken 40 times.
|
140964 | if (elemVolVars.hasVolVars(scvf.outsideScvIdx())) |
| 358 | 70442 | return StaggeredVelocityGradients::velocityGradIJ(fvGeometry, scvf, elemVolVars); | |
| 359 | else | ||
| 360 | 40 | return tangentialVelocityGradient; | |
| 361 | 70482 | }(); | |
| 362 | |||
| 363 |
2/2✓ Branch 1 taken 30440 times.
✓ Branch 2 taken 40042 times.
|
70482 | const Scalar slipVelocity = SlipVelocityPolicy::velocity(problem, fvGeometry, orthogonalScvf, elemVolVars, velGradIJ)[scvf.normalAxis()]; |
| 364 | 70482 | const Scalar velGradJI = (slipVelocity - velJ) / distance * orthogonalScvf.directionSign(); | |
| 365 | |||
| 366 |
2/2✓ Branch 0 taken 30440 times.
✓ Branch 1 taken 40042 times.
|
70482 | flux[scv.dofAxis()] -= (mu * (velGradIJ + velGradJI))*scvf.directionSign(); |
| 367 | |||
| 368 | // advective terms | ||
| 369 |
2/2✓ Branch 0 taken 30440 times.
✓ Branch 1 taken 40042 times.
|
70482 | if (problem.enableInertiaTerms()) |
| 370 | { | ||
| 371 | // transporting velocity corresponds to velJ | ||
| 372 | 30440 | const auto transportingVelocity = slipVelocity; | |
| 373 | |||
| 374 | // if the scvf lies on a boundary and if there are outside volvars, we assume that these come from a Dirichlet condition | ||
| 375 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 30440 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
30440 | if (scvf.boundary() && elemVolVars.hasVolVars(scvf.outsideScvIdx())) |
| 376 | { | ||
| 377 | ✗ | flux[scv.dofAxis()] += problem.density(fvGeometry.element(), scv) | |
| 378 | ✗ | * elemVolVars[scvf.outsideScvIdx()].velocity() * transportingVelocity * scvf.directionSign(); // TODO revise density | |
| 379 | ✗ | return flux; | |
| 380 | } | ||
| 381 | |||
| 382 | 60880 | const auto innerVelocity = elemVolVars[scvf.insideScvIdx()].velocity(); | |
| 383 | 30440 | const auto outerVelocity = [&] | |
| 384 | { | ||
| 385 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 30440 times.
|
60880 | if (!elemVolVars.hasVolVars(scvf.outsideScvIdx())) |
| 386 | ✗ | return innerVelocity; // fallback | |
| 387 | else | ||
| 388 | 60880 | return elemVolVars[scvf.outsideScvIdx()].velocity(); | |
| 389 | 30440 | }(); | |
| 390 | |||
| 391 | 30440 | const auto rho = problem.insideAndOutsideDensity(fvGeometry.element(), fvGeometry, scvf); | |
| 392 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 30437 times.
|
30440 | const bool selfIsUpstream = scvf.directionSign() == sign(transportingVelocity); |
| 393 | |||
| 394 | 30440 | const auto insideMomentum = innerVelocity * rho.first; | |
| 395 | 30440 | const auto outsideMomentum = outerVelocity * rho.second; | |
| 396 | |||
| 397 |
5/8✓ Branch 0 taken 3 times.
✓ Branch 1 taken 30437 times.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 3 times.
✗ Branch 10 not taken.
|
30440 | static const auto upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(), "Flux.UpwindWeight"); |
| 398 | |||
| 399 |
2/2✓ Branch 0 taken 12453 times.
✓ Branch 1 taken 17987 times.
|
30440 | const auto transportedMomentum = selfIsUpstream ? (upwindWeight * insideMomentum + (1.0 - upwindWeight) * outsideMomentum) |
| 400 | 17987 | : (upwindWeight * outsideMomentum + (1.0 - upwindWeight) * insideMomentum); | |
| 401 | |||
| 402 | 60880 | flux[scv.dofAxis()] += transportingVelocity * transportedMomentum * scvf.directionSign(); | |
| 403 | } | ||
| 404 | } | ||
| 405 | else | ||
| 406 | DUNE_THROW(Dune::InvalidStateException, "SlipVelocityPolicy needs to be specified as template argument"); | ||
| 407 | } | ||
| 408 | |||
| 409 | 147310 | return flux; | |
| 410 | } | ||
| 411 | }; | ||
| 412 | |||
| 413 | using NavierStokesMomentumBoundaryFluxHelper | ||
| 414 | [[deprecated("Replace with implementation class `NavierStokesMomentumBoundaryFlux`with template arguments `DiscretizationMethod` and `SlipVelocityPolicy`. This will be removed after 3.9.")]] | ||
| 415 | = NavierStokesMomentumBoundaryFlux<DiscretizationMethods::FCStaggered, | ||
| 416 | NavierStokesSlipVelocity<DiscretizationMethods::FCStaggered, NavierStokes::SlipConditions::BJ>>; | ||
| 417 | |||
| 418 | } // end namespace Dumux | ||
| 419 | |||
| 420 | #endif | ||
| 421 |