| 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 FreeFlowPoreNetworkCoupling | ||
| 10 | * \copydoc Dumux::FreeFlowPoreNetworkCouplingConditions | ||
| 11 | */ | ||
| 12 | |||
| 13 | #ifndef DUMUX_MD_FREEFLOW_PORENETWORK_COUPLINGCONDITIONS_HH | ||
| 14 | #define DUMUX_MD_FREEFLOW_PORENETWORK_COUPLINGCONDITIONS_HH | ||
| 15 | |||
| 16 | #include <numeric> | ||
| 17 | |||
| 18 | #include <dune/common/exceptions.hh> | ||
| 19 | |||
| 20 | #include <dumux/common/properties.hh> | ||
| 21 | #include <dumux/common/math.hh> | ||
| 22 | #include <dumux/discretization/cellcentered/tpfa/computetransmissibility.hh> | ||
| 23 | #include <dumux/flux/fickslaw_fwd.hh> | ||
| 24 | #include <dumux/freeflow/navierstokes/momentum/velocitygradients.hh> | ||
| 25 | #include <dumux/multidomain/boundary/freeflowporousmedium/traits.hh> | ||
| 26 | #include <dumux/multidomain/boundary/freeflowporousmedium/indexhelper.hh> | ||
| 27 | |||
| 28 | namespace Dumux { | ||
| 29 | |||
| 30 | template<class MDTraits, class CouplingManager, bool enableEnergyBalance, bool isCompositional> | ||
| 31 | class FreeFlowPoreNetworkCouplingConditionsImplementation; | ||
| 32 | |||
| 33 | /*! | ||
| 34 | * \ingroup FreeFlowPoreNetworkCoupling | ||
| 35 | * \brief Coupling conditions for the coupling of a pore-network model | ||
| 36 | * with a (Navier-)Stokes model (staggerd grid). | ||
| 37 | */ | ||
| 38 | template<class MDTraits, class CouplingManager> | ||
| 39 | using FreeFlowPoreNetworkCouplingConditions | ||
| 40 | = FreeFlowPoreNetworkCouplingConditionsImplementation< | ||
| 41 | MDTraits, CouplingManager, | ||
| 42 | GetPropType<typename MDTraits::template SubDomain<0>::TypeTag, Properties::ModelTraits>::enableEnergyBalance(), | ||
| 43 | (GetPropType<typename MDTraits::template SubDomain<0>::TypeTag, Properties::ModelTraits>::numFluidComponents() > 1) | ||
| 44 | >; | ||
| 45 | |||
| 46 | /*! | ||
| 47 | * \ingroup FreeFlowPoreNetworkCoupling | ||
| 48 | * \brief A base class which provides some common methods used for free-flow/pore-network coupling. | ||
| 49 | */ | ||
| 50 | template<class MDTraits, class CouplingManager> | ||
| 51 | class FreeFlowPoreNetworkCouplingConditionsImplementationBase | ||
| 52 | { | ||
| 53 | using Scalar = typename MDTraits::Scalar; | ||
| 54 | |||
| 55 | template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag; | ||
| 56 | template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>; | ||
| 57 | template<std::size_t id> using Element = typename GridGeometry<id>::GridView::template Codim<0>::Entity; | ||
| 58 | template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView; | ||
| 59 | template<std::size_t id> using SubControlVolumeFace = typename GridGeometry<id>::LocalView::SubControlVolumeFace; | ||
| 60 | template<std::size_t id> using SubControlVolume = typename GridGeometry<id>::LocalView::SubControlVolume; | ||
| 61 | template<std::size_t id> using Indices = typename GetPropType<SubDomainTypeTag<id>, Properties::ModelTraits>::Indices; | ||
| 62 | template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView; | ||
| 63 | template<std::size_t id> using VolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::VolumeVariables; | ||
| 64 | template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>; | ||
| 65 | template<std::size_t id> using FluidSystem = GetPropType<SubDomainTypeTag<id>, Properties::FluidSystem>; | ||
| 66 | template<std::size_t id> using ModelTraits = GetPropType<SubDomainTypeTag<id>, Properties::ModelTraits>; | ||
| 67 | template<std::size_t id> using GlobalPosition = typename Element<id>::Geometry::GlobalCoordinate; | ||
| 68 | template<std::size_t id> using NumEqVector = typename Problem<id>::Traits::NumEqVector; | ||
| 69 | |||
| 70 | using VelocityGradients = StaggeredVelocityGradients; | ||
| 71 | |||
| 72 | public: | ||
| 73 | static constexpr auto freeFlowMomentumIndex = CouplingManager::freeFlowMomentumIndex; | ||
| 74 | static constexpr auto freeFlowMassIndex = CouplingManager::freeFlowMassIndex; | ||
| 75 | static constexpr auto poreNetworkIndex = CouplingManager::poreNetworkIndex; | ||
| 76 | private: | ||
| 77 | |||
| 78 | using AdvectionType = GetPropType<SubDomainTypeTag<poreNetworkIndex>, Properties::AdvectionType>; | ||
| 79 | |||
| 80 | static constexpr bool adapterUsed = ModelTraits<poreNetworkIndex>::numFluidPhases() > 1; | ||
| 81 | using IndexHelper = FreeFlowPorousMediumCoupling::IndexHelper<freeFlowMassIndex, poreNetworkIndex, FluidSystem<freeFlowMassIndex>, adapterUsed>; | ||
| 82 | |||
| 83 | static constexpr int enableEnergyBalance = GetPropType<SubDomainTypeTag<freeFlowMassIndex>, Properties::ModelTraits>::enableEnergyBalance(); | ||
| 84 | static_assert(GetPropType<SubDomainTypeTag<poreNetworkIndex>, Properties::ModelTraits>::enableEnergyBalance() == enableEnergyBalance, | ||
| 85 | "All submodels must both be either isothermal or non-isothermal"); | ||
| 86 | |||
| 87 | static_assert(FreeFlowPorousMediumCoupling::IsSameFluidSystem<FluidSystem<freeFlowMassIndex>, | ||
| 88 | FluidSystem<poreNetworkIndex>>::value, | ||
| 89 | "All submodels must use the same fluid system"); | ||
| 90 | |||
| 91 | using VelocityVector = GlobalPosition<freeFlowMomentumIndex>; | ||
| 92 | |||
| 93 | public: | ||
| 94 | /*! | ||
| 95 | * \brief Returns the corresponding phase index needed for coupling. | ||
| 96 | */ | ||
| 97 | template<std::size_t i> | ||
| 98 | static constexpr auto couplingPhaseIdx(Dune::index_constant<i> id, int coupledPhaseIdx = 0) | ||
| 99 | { return IndexHelper::couplingPhaseIdx(id, coupledPhaseIdx); } | ||
| 100 | |||
| 101 | /*! | ||
| 102 | * \brief Returns the corresponding component index needed for coupling. | ||
| 103 | */ | ||
| 104 | template<std::size_t i> | ||
| 105 | static constexpr auto couplingCompIdx(Dune::index_constant<i> id, int coupledCompIdx) | ||
| 106 | { return IndexHelper::couplingCompIdx(id, coupledCompIdx); } | ||
| 107 | |||
| 108 | |||
| 109 | /*! | ||
| 110 | * \brief Returns the momentum flux across the coupling boundary. | ||
| 111 | * | ||
| 112 | * For the normal momentum coupling, the porous medium side of the coupling condition | ||
| 113 | * is evaluated, i.e. -[p n]^pm. | ||
| 114 | * | ||
| 115 | */ | ||
| 116 | template<class Context> | ||
| 117 | 885 | static NumEqVector<freeFlowMomentumIndex> momentumCouplingCondition(const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry, | |
| 118 | const SubControlVolumeFace<freeFlowMomentumIndex>& scvf, | ||
| 119 | const ElementVolumeVariables<freeFlowMomentumIndex>& elemVolVars, | ||
| 120 | const Context& context) | ||
| 121 | { | ||
| 122 | 885 | NumEqVector<freeFlowMomentumIndex> momentumFlux(0.0); | |
| 123 | 885 | const auto pnmPhaseIdx = couplingPhaseIdx(poreNetworkIndex); | |
| 124 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 885 times.
|
885 | const auto [pnmPressure, pnmViscosity] = [&] |
| 125 | { | ||
| 126 |
1/2✓ Branch 0 taken 1770 times.
✗ Branch 1 not taken.
|
2655 | for (const auto& scv : scvs(context.fvGeometry)) |
| 127 | { | ||
| 128 |
2/2✓ Branch 0 taken 885 times.
✓ Branch 1 taken 885 times.
|
1770 | if (scv.dofIndex() == context.poreNetworkDofIdx) |
| 129 | 5310 | return std::make_pair(context.elemVolVars[scv].pressure(pnmPhaseIdx), context.elemVolVars[scv].viscosity(pnmPhaseIdx)); | |
| 130 | } | ||
| 131 | ✗ | DUNE_THROW(Dune::InvalidStateException, "No coupled scv found"); | |
| 132 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 885 times.
|
885 | }(); |
| 133 | |||
| 134 | // set p_freeflow = p_PNM | ||
| 135 | 885 | momentumFlux[scvf.normalAxis()] = pnmPressure; | |
| 136 | |||
| 137 | // normalize pressure | ||
| 138 | 885 | momentumFlux[scvf.normalAxis()] -= elemVolVars.gridVolVars().problem().referencePressure(fvGeometry.element(), fvGeometry, scvf); | |
| 139 | |||
| 140 | // Explicitly account for dv_i/dx_i, which is NOT part of the actual coupling condition. We do it here for convenience so | ||
| 141 | // we do not forget to set it in the problem. We assume that the velocity gradient at the boundary towards the interface is the same | ||
| 142 | // as the one in the center of the element. TODO check sign | ||
| 143 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 885 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 885 times.
|
1770 | const auto& scv = fvGeometry.scv(scvf.insideScvIdx()); |
| 144 | 885 | const auto& frontalInternalScvf = (*scvfs(fvGeometry, scv).begin()); | |
| 145 | 885 | momentumFlux[scvf.normalAxis()] -= 2*VelocityGradients::velocityGradII(fvGeometry, frontalInternalScvf, elemVolVars) * pnmViscosity; | |
| 146 | |||
| 147 | // We do NOT consider the inertia term here. If included, Newton convergence decreases drastically and the solution even does not converge to a reference solution. | ||
| 148 | // We furthermore assume creeping flow within the boundary layer thus neglecting this term is physically justified. | ||
| 149 | |||
| 150 | 885 | momentumFlux[scvf.normalAxis()] *= scvf.directionSign(); | |
| 151 | |||
| 152 | 885 | return momentumFlux; | |
| 153 | } | ||
| 154 | |||
| 155 | template<class Context> | ||
| 156 | 3606 | static VelocityVector interfaceThroatVelocity(const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry, | |
| 157 | const SubControlVolumeFace<freeFlowMomentumIndex>& scvf, | ||
| 158 | const Context& context) | ||
| 159 | { | ||
| 160 |
1/2✓ Branch 0 taken 3606 times.
✗ Branch 1 not taken.
|
3606 | const auto& pnmElement = context.fvGeometry.element(); |
| 161 | 3606 | const auto& pnmFVGeometry = context.fvGeometry; | |
| 162 |
1/2✓ Branch 0 taken 3606 times.
✗ Branch 1 not taken.
|
3606 | const auto& pnmScvf = pnmFVGeometry.scvf(0); |
| 163 | 3606 | const auto& pnmElemVolVars = context.elemVolVars; | |
| 164 | 3606 | const auto& pnmElemFluxVarsCache = context.elemFluxVarsCache; | |
| 165 | 3606 | const auto& pnmProblem = pnmElemVolVars.gridVolVars().problem(); | |
| 166 | |||
| 167 | 3606 | const auto pnmPhaseIdx = couplingPhaseIdx(poreNetworkIndex); | |
| 168 |
1/2✓ Branch 0 taken 3606 times.
✗ Branch 1 not taken.
|
3606 | const Scalar area = pnmElemFluxVarsCache[pnmScvf].throatCrossSectionalArea(pnmPhaseIdx); |
| 169 | |||
| 170 | // only proceed if area > 0 in order to prevent division by zero (e.g., when the throat was not invaded yet) | ||
| 171 |
1/2✓ Branch 0 taken 3606 times.
✗ Branch 1 not taken.
|
3606 | if (area > 0.0) |
| 172 | { | ||
| 173 | using PNMFluxVariables = GetPropType<SubDomainTypeTag<poreNetworkIndex>, Properties::FluxVariables>; | ||
| 174 | 3606 | PNMFluxVariables fluxVars; | |
| 175 | 3606 | fluxVars.init(pnmProblem, pnmElement, pnmFVGeometry, pnmElemVolVars, pnmScvf, pnmElemFluxVarsCache); | |
| 176 | |||
| 177 | 18030 | const Scalar flux = fluxVars.advectiveFlux(pnmPhaseIdx, [pnmPhaseIdx](const auto& volVars){ return volVars.mobility(pnmPhaseIdx);}); | |
| 178 | |||
| 179 | // account for the orientation of the bulk face. | ||
| 180 | 3606 | VelocityVector velocity = (pnmElement.geometry().corner(1) - pnmElement.geometry().corner(0)); | |
| 181 | 7212 | velocity /= velocity.two_norm(); | |
| 182 | 3606 | velocity *= flux / area; | |
| 183 | |||
| 184 | // TODO: Multiple throats connected to the same pore? | ||
| 185 | 3606 | return velocity; | |
| 186 | } | ||
| 187 | else | ||
| 188 | ✗ | return VelocityVector(0.0); | |
| 189 | } | ||
| 190 | |||
| 191 | /*! | ||
| 192 | * \brief Evaluate an advective flux across the interface and consider upwinding. | ||
| 193 | */ | ||
| 194 | static Scalar advectiveFlux(const Scalar insideQuantity, const Scalar outsideQuantity, const Scalar volumeFlow, bool insideIsUpstream) | ||
| 195 | { | ||
| 196 | 648 | const Scalar upwindWeight = 1.0; //TODO use Flux.UpwindWeight or something like Coupling.UpwindWeight? | |
| 197 | |||
| 198 | 648 | if(insideIsUpstream) | |
| 199 | 344 | return (upwindWeight * insideQuantity + (1.0 - upwindWeight) * outsideQuantity) * volumeFlow; | |
| 200 | else | ||
| 201 | 304 | return (upwindWeight * outsideQuantity + (1.0 - upwindWeight) * insideQuantity) * volumeFlow; | |
| 202 | } | ||
| 203 | |||
| 204 | protected: | ||
| 205 | |||
| 206 | /*! | ||
| 207 | * \brief Returns the distance between an scvf and the corresponding scv center. | ||
| 208 | */ | ||
| 209 | template<class Scv, class Scvf> | ||
| 210 | Scalar getDistance_(const Scv& scv, const Scvf& scvf) const | ||
| 211 | { | ||
| 212 | return (scv.dofPosition() - scvf.ipGlobal()).two_norm(); | ||
| 213 | } | ||
| 214 | }; | ||
| 215 | |||
| 216 | /*! | ||
| 217 | * \ingroup FreeFlowPoreNetworkCoupling | ||
| 218 | * \brief Coupling conditions specialization for non-compositional models. | ||
| 219 | */ | ||
| 220 | template<class MDTraits, class CouplingManager, bool enableEnergyBalance> | ||
| 221 | class FreeFlowPoreNetworkCouplingConditionsImplementation<MDTraits, CouplingManager, enableEnergyBalance, false> | ||
| 222 | : public FreeFlowPoreNetworkCouplingConditionsImplementationBase<MDTraits, CouplingManager> | ||
| 223 | { | ||
| 224 | using ParentType = FreeFlowPoreNetworkCouplingConditionsImplementationBase<MDTraits, CouplingManager>; | ||
| 225 | using Scalar = typename MDTraits::Scalar; | ||
| 226 | |||
| 227 | // the sub domain type tags | ||
| 228 | template<std::size_t id> | ||
| 229 | using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag; | ||
| 230 | |||
| 231 | template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>; | ||
| 232 | template<std::size_t id> using Element = typename GridGeometry<id>::GridView::template Codim<0>::Entity; | ||
| 233 | template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView; | ||
| 234 | template<std::size_t id> using SubControlVolumeFace = typename GridGeometry<id>::LocalView::SubControlVolumeFace; | ||
| 235 | template<std::size_t id> using SubControlVolume = typename GridGeometry<id>::LocalView::SubControlVolume; | ||
| 236 | template<std::size_t id> using Indices = typename GetPropType<SubDomainTypeTag<id>, Properties::ModelTraits>::Indices; | ||
| 237 | template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView; | ||
| 238 | template<std::size_t id> using VolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::VolumeVariables; | ||
| 239 | |||
| 240 | static_assert(GetPropType<SubDomainTypeTag<ParentType::poreNetworkIndex>, Properties::ModelTraits>::numFluidComponents() == GetPropType<SubDomainTypeTag<ParentType::poreNetworkIndex>, Properties::ModelTraits>::numFluidPhases(), | ||
| 241 | "Pore-network model must not be compositional"); | ||
| 242 | |||
| 243 | public: | ||
| 244 | using ParentType::ParentType; | ||
| 245 | using ParentType::couplingPhaseIdx; | ||
| 246 | |||
| 247 | /*! | ||
| 248 | * \brief Returns the mass flux across the coupling boundary as seen from the pore-network domain. | ||
| 249 | */ | ||
| 250 | template<class CouplingContext> | ||
| 251 | 168 | static Scalar massCouplingCondition(Dune::index_constant<ParentType::poreNetworkIndex> domainI, | |
| 252 | Dune::index_constant<ParentType::freeFlowMassIndex> domainJ, | ||
| 253 | const FVElementGeometry<ParentType::poreNetworkIndex>& fvGeometry, | ||
| 254 | const SubControlVolume<ParentType::poreNetworkIndex>& scv, | ||
| 255 | const ElementVolumeVariables<ParentType::poreNetworkIndex>& insideVolVars, | ||
| 256 | const CouplingContext& context) | ||
| 257 | { | ||
| 258 | 168 | Scalar massFlux(0.0); | |
| 259 | |||
| 260 |
4/4✓ Branch 0 taken 648 times.
✓ Branch 1 taken 168 times.
✓ Branch 2 taken 648 times.
✓ Branch 3 taken 168 times.
|
1800 | for (const auto& c : context) |
| 261 | { | ||
| 262 | // positive values indicate flux into pore-network region | ||
| 263 | 648 | const Scalar normalFFVelocity = c.velocity * c.scvf.unitOuterNormal(); | |
| 264 |
2/2✓ Branch 0 taken 344 times.
✓ Branch 1 taken 304 times.
|
648 | const bool pnmIsUpstream = std::signbit(normalFFVelocity); |
| 265 | |||
| 266 |
4/4✓ Branch 0 taken 344 times.
✓ Branch 1 taken 304 times.
✓ Branch 2 taken 344 times.
✓ Branch 3 taken 304 times.
|
1296 | const Scalar pnmDensity = insideVolVars[scv].density(couplingPhaseIdx(ParentType::poreNetworkIndex)); |
| 267 |
2/2✓ Branch 0 taken 344 times.
✓ Branch 1 taken 304 times.
|
648 | const Scalar ffDensity = c.volVars.density(couplingPhaseIdx(ParentType::freeFlowMassIndex)); |
| 268 | 648 | const Scalar area = c.scvf.area() * c.volVars.extrusionFactor(); | |
| 269 | |||
| 270 | // flux is used as source term: positive values mean influx | ||
| 271 |
2/2✓ Branch 0 taken 344 times.
✓ Branch 1 taken 304 times.
|
1296 | massFlux += ParentType::advectiveFlux(pnmDensity, ffDensity, normalFFVelocity*area, pnmIsUpstream); |
| 272 | } | ||
| 273 | |||
| 274 | 168 | return massFlux; | |
| 275 | } | ||
| 276 | |||
| 277 | /*! | ||
| 278 | * \brief Returns the mass flux across the coupling boundary as seen from the free-flow domain. | ||
| 279 | */ | ||
| 280 | template<class CouplingContext> | ||
| 281 | ✗ | static Scalar massCouplingCondition(Dune::index_constant<ParentType::freeFlowMassIndex> domainI, | |
| 282 | Dune::index_constant<ParentType::poreNetworkIndex> domainJ, | ||
| 283 | const FVElementGeometry<ParentType::freeFlowMassIndex>& fvGeometry, | ||
| 284 | const SubControlVolumeFace<ParentType::freeFlowMassIndex>& scvf, | ||
| 285 | const ElementVolumeVariables<ParentType::freeFlowMassIndex>& insideVolVars, | ||
| 286 | const CouplingContext& context) | ||
| 287 | { | ||
| 288 | // positive values indicate flux into pore-network region | ||
| 289 | ✗ | const Scalar normalFFVelocity = context.velocity * scvf.unitOuterNormal(); | |
| 290 | ✗ | const bool ffIsUpstream = !std::signbit(normalFFVelocity); | |
| 291 | |||
| 292 | ✗ | const Scalar ffDensity = insideVolVars[scvf.insideScvIdx()].density(couplingPhaseIdx(ParentType::freeFlowMassIndex)); | |
| 293 | ✗ | const Scalar pnmDensity = context.volVars.density(couplingPhaseIdx(ParentType::poreNetworkIndex)); | |
| 294 | |||
| 295 | // flux is used in Neumann condition: positive values mean flux out of free-flow domain | ||
| 296 | ✗ | return ParentType::advectiveFlux(ffDensity, pnmDensity, normalFFVelocity, ffIsUpstream); | |
| 297 | } | ||
| 298 | |||
| 299 | private: | ||
| 300 | |||
| 301 | }; | ||
| 302 | |||
| 303 | } // end namespace Dumux | ||
| 304 | |||
| 305 | #endif | ||
| 306 |