| 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 FreeFlowPorousMediumCoupling | ||
| 10 | * \copydoc Dumux::FFPMCouplingConditionsStaggeredCCTpfa | ||
| 11 | */ | ||
| 12 | |||
| 13 | #ifndef DUMUX_MD_FREEFLOW_POROUSMEDIUM_COUPLINGCONDITIONS_STAGGERED_TPFA_HH | ||
| 14 | #define DUMUX_MD_FREEFLOW_POROUSMEDIUM_COUPLINGCONDITIONS_STAGGERED_TPFA_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 | |||
| 23 | #include <dumux/discretization/cellcentered/tpfa/computetransmissibility.hh> | ||
| 24 | |||
| 25 | #include <dumux/flux/darcyslaw_fwd.hh> | ||
| 26 | #include <dumux/flux/fickslaw_fwd.hh> | ||
| 27 | #include <dumux/flux/forchheimerslaw_fwd.hh> | ||
| 28 | |||
| 29 | #include <dumux/multidomain/boundary/freeflowporousmedium/traits.hh> | ||
| 30 | #include <dumux/multidomain/boundary/freeflowporousmedium/indexhelper.hh> | ||
| 31 | |||
| 32 | namespace Dumux { | ||
| 33 | |||
| 34 | /*! | ||
| 35 | * \ingroup FreeFlowPorousMediumCoupling | ||
| 36 | * \brief This structs holds a set of options which allow to modify the Stokes-Darcy | ||
| 37 | * coupling mechanism during runtime. | ||
| 38 | */ | ||
| 39 | struct FreeFlowPorousMediumCouplingOptions | ||
| 40 | { | ||
| 41 | /*! | ||
| 42 | * \brief Defines which kind of averanging of diffusion coefficients | ||
| 43 | * (moleculat diffusion or thermal conductance) at the interface | ||
| 44 | * between free flow and porous medium shall be used. | ||
| 45 | */ | ||
| 46 | enum class DiffusionCoefficientAveragingType | ||
| 47 | { | ||
| 48 | harmonic, arithmetic, ffOnly, pmOnly | ||
| 49 | }; | ||
| 50 | |||
| 51 | /*! | ||
| 52 | * \brief Convenience function to convert user input given as std::string to the corresponding enum class used for choosing the type | ||
| 53 | * of averaging of the diffusion/conduction parameter at the interface between the two domains. | ||
| 54 | */ | ||
| 55 | static DiffusionCoefficientAveragingType stringToEnum(DiffusionCoefficientAveragingType, const std::string& diffusionCoefficientAveragingType) | ||
| 56 | { | ||
| 57 | if (diffusionCoefficientAveragingType == "Harmonic") | ||
| 58 | return DiffusionCoefficientAveragingType::harmonic; | ||
| 59 | else if (diffusionCoefficientAveragingType == "Arithmetic") | ||
| 60 | return DiffusionCoefficientAveragingType::arithmetic; | ||
| 61 | else if (diffusionCoefficientAveragingType == "FreeFlowOnly") | ||
| 62 | return DiffusionCoefficientAveragingType::ffOnly; | ||
| 63 | else if (diffusionCoefficientAveragingType == "PorousMediumOnly") | ||
| 64 | return DiffusionCoefficientAveragingType::pmOnly; | ||
| 65 | else | ||
| 66 | DUNE_THROW(Dune::IOError, "Unknown DiffusionCoefficientAveragingType"); | ||
| 67 | } | ||
| 68 | |||
| 69 | }; | ||
| 70 | |||
| 71 | template<class MDTraits, class CouplingManager, bool enableEnergyBalance, bool isCompositional> | ||
| 72 | class FFPMCouplingConditionsStaggeredCCTpfaImpl; | ||
| 73 | |||
| 74 | /*! | ||
| 75 | * \ingroup FreeFlowPorousMediumCoupling | ||
| 76 | * \brief Data for the coupling of a Darcy model (cell-centered finite volume) | ||
| 77 | * with a (Navier-)Stokes model (staggerd grid). | ||
| 78 | */ | ||
| 79 | template<class MDTraits, class CouplingManager> | ||
| 80 | using FFPMCouplingConditionsStaggeredCCTpfa | ||
| 81 | = FFPMCouplingConditionsStaggeredCCTpfaImpl< | ||
| 82 | MDTraits, CouplingManager, | ||
| 83 | GetPropType<typename MDTraits::template SubDomain<0>::TypeTag, Properties::ModelTraits>::enableEnergyBalance(), | ||
| 84 | (GetPropType<typename MDTraits::template SubDomain<0>::TypeTag, Properties::ModelTraits>::numFluidComponents() > 1) | ||
| 85 | >; | ||
| 86 | |||
| 87 | /*! | ||
| 88 | * \ingroup FreeFlowPorousMediumCoupling | ||
| 89 | * \brief A base class which provides some common methods used for Stokes-Darcy coupling. | ||
| 90 | */ | ||
| 91 | template<class MDTraits, class CouplingManager> | ||
| 92 | class FFPMCouplingConditionsStaggeredCCTpfaImplBase | ||
| 93 | { | ||
| 94 | using Scalar = typename MDTraits::Scalar; | ||
| 95 | |||
| 96 | template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag; | ||
| 97 | template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>; | ||
| 98 | template<std::size_t id> using Element = typename GridGeometry<id>::GridView::template Codim<0>::Entity; | ||
| 99 | template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView; | ||
| 100 | template<std::size_t id> using SubControlVolumeFace = typename GridGeometry<id>::LocalView::SubControlVolumeFace; | ||
| 101 | template<std::size_t id> using SubControlVolume = typename GridGeometry<id>::LocalView::SubControlVolume; | ||
| 102 | template<std::size_t id> using Indices = typename GetPropType<SubDomainTypeTag<id>, Properties::ModelTraits>::Indices; | ||
| 103 | template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView; | ||
| 104 | template<std::size_t id> using VolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::VolumeVariables; | ||
| 105 | template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>; | ||
| 106 | template<std::size_t id> using FluidSystem = GetPropType<SubDomainTypeTag<id>, Properties::FluidSystem>; | ||
| 107 | template<std::size_t id> using ModelTraits = GetPropType<SubDomainTypeTag<id>, Properties::ModelTraits>; | ||
| 108 | template<std::size_t id> using GlobalPosition = typename Element<id>::Geometry::GlobalCoordinate; | ||
| 109 | template<std::size_t id> using NumEqVector = typename Problem<id>::Traits::NumEqVector; | ||
| 110 | |||
| 111 | public: | ||
| 112 | static constexpr auto freeFlowMomentumIndex = CouplingManager::freeFlowMomentumIndex; | ||
| 113 | static constexpr auto freeFlowMassIndex = CouplingManager::freeFlowMassIndex; | ||
| 114 | static constexpr auto porousMediumIndex = CouplingManager::porousMediumIndex; | ||
| 115 | private: | ||
| 116 | |||
| 117 | using AdvectionType = GetPropType<SubDomainTypeTag<porousMediumIndex>, Properties::AdvectionType>; | ||
| 118 | using DarcysLaw = Dumux::DarcysLaw<SubDomainTypeTag<porousMediumIndex>>; | ||
| 119 | using ForchheimersLaw = Dumux::ForchheimersLaw<SubDomainTypeTag<porousMediumIndex>>; | ||
| 120 | |||
| 121 | static constexpr bool adapterUsed = ModelTraits<porousMediumIndex>::numFluidPhases() > 1; | ||
| 122 | using IndexHelper = FreeFlowPorousMediumCoupling::IndexHelper<freeFlowMassIndex, porousMediumIndex, FluidSystem<freeFlowMassIndex>, adapterUsed>; | ||
| 123 | |||
| 124 | static constexpr int enableEnergyBalance = GetPropType<SubDomainTypeTag<freeFlowMassIndex>, Properties::ModelTraits>::enableEnergyBalance(); | ||
| 125 | static_assert(GetPropType<SubDomainTypeTag<porousMediumIndex>, Properties::ModelTraits>::enableEnergyBalance() == enableEnergyBalance, | ||
| 126 | "All submodels must both be either isothermal or non-isothermal"); | ||
| 127 | |||
| 128 | static_assert(FreeFlowPorousMediumCoupling::IsSameFluidSystem<FluidSystem<freeFlowMassIndex>, | ||
| 129 | FluidSystem<porousMediumIndex>>::value, | ||
| 130 | "All submodels must use the same fluid system"); | ||
| 131 | |||
| 132 | using DiffusionCoefficientAveragingType = typename FreeFlowPorousMediumCouplingOptions::DiffusionCoefficientAveragingType; | ||
| 133 | |||
| 134 | public: | ||
| 135 | /*! | ||
| 136 | * \brief Returns the corresponding phase index needed for coupling. | ||
| 137 | */ | ||
| 138 | template<std::size_t i> | ||
| 139 | static constexpr auto couplingPhaseIdx(Dune::index_constant<i> id, int coupledPhaseIdx = 0) | ||
| 140 | { return IndexHelper::couplingPhaseIdx(id, coupledPhaseIdx); } | ||
| 141 | |||
| 142 | /*! | ||
| 143 | * \brief Returns the corresponding component index needed for coupling. | ||
| 144 | */ | ||
| 145 | template<std::size_t i> | ||
| 146 | static constexpr auto couplingCompIdx(Dune::index_constant<i> id, int coupledCompIdx) | ||
| 147 | { return IndexHelper::couplingCompIdx(id, coupledCompIdx); } | ||
| 148 | |||
| 149 | /*! | ||
| 150 | * \brief Returns the intrinsic permeability of the coupled Darcy element. | ||
| 151 | */ | ||
| 152 | template<class Context> | ||
| 153 | ✗ | static auto darcyPermeability(const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry, | |
| 154 | const SubControlVolumeFace<freeFlowMomentumIndex>& scvf, | ||
| 155 | const Context& context) | ||
| 156 | 143704 | { return context.volVars.permeability(); } | |
| 157 | |||
| 158 | /*! | ||
| 159 | * \brief Returns the momentum flux across the coupling boundary. | ||
| 160 | * | ||
| 161 | * For the normal momentum coupling, the porous medium side of the coupling condition | ||
| 162 | * is evaluated, i.e. -[p n]^pm. | ||
| 163 | * | ||
| 164 | */ | ||
| 165 | template<class Context> | ||
| 166 | static NumEqVector<freeFlowMomentumIndex> momentumCouplingCondition(const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry, | ||
| 167 | const SubControlVolumeFace<freeFlowMomentumIndex>& scvf, | ||
| 168 | const ElementVolumeVariables<freeFlowMomentumIndex>& elemVolVars, | ||
| 169 | const Context& context) | ||
| 170 | { | ||
| 171 | static constexpr auto numPhasesDarcy = GetPropType<SubDomainTypeTag<porousMediumIndex>, Properties::ModelTraits>::numFluidPhases(); | ||
| 172 |
1/2✓ Branch 0 taken 34686 times.
✗ Branch 1 not taken.
|
34686 | NumEqVector<freeFlowMomentumIndex> momentumFlux(0.0); |
| 173 | |||
| 174 |
1/2✓ Branch 0 taken 34686 times.
✗ Branch 1 not taken.
|
34686 | if (!scvf.isFrontal()) |
| 175 | return momentumFlux; | ||
| 176 | |||
| 177 | 34686 | const auto pmPhaseIdx = couplingPhaseIdx(porousMediumIndex); | |
| 178 | |||
| 179 | if(numPhasesDarcy > 1) | ||
| 180 | momentumFlux[scvf.normalAxis()] = context.volVars.pressure(pmPhaseIdx); | ||
| 181 | else // use pressure reconstruction for single phase models | ||
| 182 | 34686 | momentumFlux[scvf.normalAxis()] = pressureAtInterface_(fvGeometry, scvf, elemVolVars, context); | |
| 183 | |||
| 184 | // TODO: generalize for permeability tensors | ||
| 185 | |||
| 186 | // normalize pressure | ||
| 187 | 34686 | momentumFlux[scvf.normalAxis()] -= elemVolVars.gridVolVars().problem().referencePressure(fvGeometry.element(), fvGeometry, scvf); | |
| 188 | |||
| 189 | 34686 | momentumFlux[scvf.normalAxis()] *= scvf.directionSign(); | |
| 190 | |||
| 191 | return momentumFlux; | ||
| 192 | } | ||
| 193 | |||
| 194 | /*! | ||
| 195 | * \brief Evaluate an advective flux across the interface and consider upwinding. | ||
| 196 | */ | ||
| 197 | static Scalar advectiveFlux(const Scalar insideQuantity, const Scalar outsideQuantity, const Scalar volumeFlow, bool insideIsUpstream) | ||
| 198 | { | ||
| 199 | ✗ | const Scalar upwindWeight = 1.0; //TODO use Flux.UpwindWeight or something like Coupling.UpwindWeight? | |
| 200 | |||
| 201 | ✗ | if(insideIsUpstream) | |
| 202 | ✗ | return (upwindWeight * insideQuantity + (1.0 - upwindWeight) * outsideQuantity) * volumeFlow; | |
| 203 | else | ||
| 204 | ✗ | return (upwindWeight * outsideQuantity + (1.0 - upwindWeight) * insideQuantity) * volumeFlow; | |
| 205 | } | ||
| 206 | |||
| 207 | protected: | ||
| 208 | |||
| 209 | /*! | ||
| 210 | * \brief Returns the transmissibility used for either molecular diffusion or thermal conductivity. | ||
| 211 | */ | ||
| 212 | template<std::size_t i, std::size_t j> | ||
| 213 | Scalar transmissibility_(Dune::index_constant<i> domainI, | ||
| 214 | Dune::index_constant<j> domainJ, | ||
| 215 | const Scalar insideDistance, | ||
| 216 | const Scalar outsideDistance, | ||
| 217 | const Scalar avgQuantityI, | ||
| 218 | const Scalar avgQuantityJ, | ||
| 219 | const DiffusionCoefficientAveragingType diffCoeffAvgType) const | ||
| 220 | { | ||
| 221 | const Scalar totalDistance = insideDistance + outsideDistance; | ||
| 222 | if(diffCoeffAvgType == DiffusionCoefficientAveragingType::harmonic) | ||
| 223 | { | ||
| 224 | return harmonicMean(avgQuantityI, avgQuantityJ, insideDistance, outsideDistance) | ||
| 225 | / totalDistance; | ||
| 226 | } | ||
| 227 | else if(diffCoeffAvgType == DiffusionCoefficientAveragingType::arithmetic) | ||
| 228 | { | ||
| 229 | return arithmeticMean(avgQuantityI, avgQuantityJ, insideDistance, outsideDistance) | ||
| 230 | / totalDistance; | ||
| 231 | } | ||
| 232 | else if(diffCoeffAvgType == DiffusionCoefficientAveragingType::ffOnly) | ||
| 233 | return domainI == freeFlowMassIndex | ||
| 234 | ? avgQuantityI / totalDistance | ||
| 235 | : avgQuantityJ / totalDistance; | ||
| 236 | |||
| 237 | else // diffCoeffAvgType == DiffusionCoefficientAveragingType::pmOnly) | ||
| 238 | return domainI == porousMediumIndex | ||
| 239 | ? avgQuantityI / totalDistance | ||
| 240 | : avgQuantityJ / totalDistance; | ||
| 241 | } | ||
| 242 | |||
| 243 | /*! | ||
| 244 | * \brief Returns the distance between an scvf and the corresponding scv center. | ||
| 245 | */ | ||
| 246 | template<class Scv, class Scvf> | ||
| 247 | Scalar getDistance_(const Scv& scv, const Scvf& scvf) const | ||
| 248 | { | ||
| 249 | return (scv.dofPosition() - scvf.ipGlobal()).two_norm(); | ||
| 250 | } | ||
| 251 | |||
| 252 | /*! | ||
| 253 | * \brief Returns the conductive energy flux across the interface. | ||
| 254 | */ | ||
| 255 | template<std::size_t i, std::size_t j, bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0> | ||
| 256 | Scalar conductiveEnergyFlux_(Dune::index_constant<i> domainI, | ||
| 257 | Dune::index_constant<j> domainJ, | ||
| 258 | const FVElementGeometry<i>& fvGeometryI, | ||
| 259 | const FVElementGeometry<j>& fvGeometryJ, | ||
| 260 | const SubControlVolumeFace<i>& scvfI, | ||
| 261 | const SubControlVolume<i>& scvI, | ||
| 262 | const SubControlVolume<j>& scvJ, | ||
| 263 | const VolumeVariables<i>& volVarsI, | ||
| 264 | const VolumeVariables<j>& volVarsJ, | ||
| 265 | const DiffusionCoefficientAveragingType diffCoeffAvgType) const | ||
| 266 | { | ||
| 267 | const Scalar insideDistance = getDistance_(scvI, scvfI); | ||
| 268 | const Scalar outsideDistance = getDistance_(scvJ, scvfI); | ||
| 269 | |||
| 270 | const Scalar deltaT = volVarsJ.temperature() - volVarsI.temperature(); | ||
| 271 | const Scalar tij = transmissibility_( | ||
| 272 | domainI, domainJ, | ||
| 273 | insideDistance, outsideDistance, | ||
| 274 | volVarsI.effectiveThermalConductivity(), volVarsJ.effectiveThermalConductivity(), | ||
| 275 | diffCoeffAvgType | ||
| 276 | ); | ||
| 277 | |||
| 278 | return -tij * deltaT; | ||
| 279 | } | ||
| 280 | |||
| 281 | /*! | ||
| 282 | * \brief Returns the pressure at the interface | ||
| 283 | */ | ||
| 284 | template<class CouplingContext> | ||
| 285 | 34686 | static Scalar pressureAtInterface_(const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry, | |
| 286 | const SubControlVolumeFace<freeFlowMomentumIndex>& scvf, | ||
| 287 | const ElementVolumeVariables<freeFlowMomentumIndex>& elemVolVars, | ||
| 288 | const CouplingContext& context) | ||
| 289 | { | ||
| 290 | 34686 | GlobalPosition<freeFlowMomentumIndex> velocity(0.0); | |
| 291 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 34686 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 34686 times.
|
69372 | const auto& scv = fvGeometry.scv(scvf.insideScvIdx()); |
| 292 | 34686 | velocity[scv.dofAxis()] = elemVolVars[scv].velocity(); | |
| 293 | 34686 | const auto& pmScvf = context.fvGeometry.scvf(context.porousMediumScvfIdx); | |
| 294 | 34686 | return computeCouplingPhasePressureAtInterface_(context.fvGeometry, pmScvf, context.volVars, context.gravity, velocity, AdvectionType()); | |
| 295 | } | ||
| 296 | |||
| 297 | /*! | ||
| 298 | * \brief Returns the pressure at the interface using Forchheimers's law for reconstruction | ||
| 299 | */ | ||
| 300 | static Scalar computeCouplingPhasePressureAtInterface_(const FVElementGeometry<porousMediumIndex>& fvGeometry, | ||
| 301 | const SubControlVolumeFace<porousMediumIndex>& scvf, | ||
| 302 | const VolumeVariables<porousMediumIndex>& volVars, | ||
| 303 | const typename Element<freeFlowMomentumIndex>::Geometry::GlobalCoordinate& couplingPhaseVelocity, | ||
| 304 | ForchheimersLaw) | ||
| 305 | { | ||
| 306 | DUNE_THROW(Dune::NotImplemented, "Forchheimer's law pressure reconstruction"); | ||
| 307 | } | ||
| 308 | |||
| 309 | /*! | ||
| 310 | * \brief Returns the pressure at the interface using Darcy's law for reconstruction | ||
| 311 | */ | ||
| 312 | 34686 | static Scalar computeCouplingPhasePressureAtInterface_(const FVElementGeometry<porousMediumIndex>& fvGeometry, | |
| 313 | const SubControlVolumeFace<porousMediumIndex>& scvf, | ||
| 314 | const VolumeVariables<porousMediumIndex>& volVars, | ||
| 315 | const typename Element<freeFlowMomentumIndex>::Geometry::GlobalCoordinate& gravity, | ||
| 316 | const typename Element<freeFlowMomentumIndex>::Geometry::GlobalCoordinate& couplingPhaseVelocity, | ||
| 317 | DarcysLaw) | ||
| 318 | { | ||
| 319 | 34686 | const auto pmPhaseIdx = couplingPhaseIdx(porousMediumIndex); | |
| 320 | 34686 | const Scalar couplingPhaseCellCenterPressure = volVars.pressure(pmPhaseIdx); | |
| 321 | 34686 | const Scalar couplingPhaseMobility = volVars.mobility(pmPhaseIdx); | |
| 322 | 34686 | const Scalar couplingPhaseDensity = volVars.density(pmPhaseIdx); | |
| 323 | 34686 | const auto K = volVars.permeability(); | |
| 324 | |||
| 325 | // A tpfa approximation yields (works if mobility != 0) | ||
| 326 | // v*n = -kr/mu*K * (gradP - rho*g)*n = mobility*(ti*(p_center - p_interface) + rho*n^TKg) | ||
| 327 | // -> p_interface = (1/mobility * (-v*n) + rho*n^TKg)/ti + p_center | ||
| 328 | // where v is the free-flow velocity (couplingPhaseVelocity) | ||
| 329 | 69372 | const auto alpha = vtmv(scvf.unitOuterNormal(), K, gravity); | |
| 330 | |||
| 331 |
2/4✓ Branch 0 taken 34686 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 34686 times.
✗ Branch 3 not taken.
|
69372 | const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); |
| 332 | 34686 | const auto ti = computeTpfaTransmissibility(fvGeometry, scvf, insideScv, K, 1.0); | |
| 333 | |||
| 334 | 69372 | return (-1/couplingPhaseMobility * (scvf.unitOuterNormal() * couplingPhaseVelocity) + couplingPhaseDensity * alpha)/ti | |
| 335 | 34686 | + couplingPhaseCellCenterPressure; | |
| 336 | } | ||
| 337 | }; | ||
| 338 | |||
| 339 | /*! | ||
| 340 | * \ingroup FreeFlowPorousMediumCoupling | ||
| 341 | * \brief Coupling data specialization for non-compositional models. | ||
| 342 | */ | ||
| 343 | template<class MDTraits, class CouplingManager, bool enableEnergyBalance> | ||
| 344 | class FFPMCouplingConditionsStaggeredCCTpfaImpl<MDTraits, CouplingManager, enableEnergyBalance, false> | ||
| 345 | : public FFPMCouplingConditionsStaggeredCCTpfaImplBase<MDTraits, CouplingManager> | ||
| 346 | { | ||
| 347 | using ParentType = FFPMCouplingConditionsStaggeredCCTpfaImplBase<MDTraits, CouplingManager>; | ||
| 348 | using Scalar = typename MDTraits::Scalar; | ||
| 349 | |||
| 350 | // the sub domain type tags | ||
| 351 | template<std::size_t id> | ||
| 352 | using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag; | ||
| 353 | |||
| 354 | template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>; | ||
| 355 | template<std::size_t id> using Element = typename GridGeometry<id>::GridView::template Codim<0>::Entity; | ||
| 356 | template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView; | ||
| 357 | template<std::size_t id> using SubControlVolumeFace = typename GridGeometry<id>::LocalView::SubControlVolumeFace; | ||
| 358 | template<std::size_t id> using SubControlVolume = typename GridGeometry<id>::LocalView::SubControlVolume; | ||
| 359 | template<std::size_t id> using Indices = typename GetPropType<SubDomainTypeTag<id>, Properties::ModelTraits>::Indices; | ||
| 360 | template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView; | ||
| 361 | template<std::size_t id> using VolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::VolumeVariables; | ||
| 362 | |||
| 363 | static_assert(GetPropType<SubDomainTypeTag<ParentType::porousMediumIndex>, Properties::ModelTraits>::numFluidComponents() == GetPropType<SubDomainTypeTag<ParentType::porousMediumIndex>, Properties::ModelTraits>::numFluidPhases(), | ||
| 364 | "Darcy Model must not be compositional"); | ||
| 365 | |||
| 366 | using DiffusionCoefficientAveragingType = typename FreeFlowPorousMediumCouplingOptions::DiffusionCoefficientAveragingType; | ||
| 367 | |||
| 368 | public: | ||
| 369 | using ParentType::ParentType; | ||
| 370 | using ParentType::couplingPhaseIdx; | ||
| 371 | |||
| 372 | /*! | ||
| 373 | * \brief Returns the mass flux across the coupling boundary as seen from the Darcy domain. | ||
| 374 | */ | ||
| 375 | template<std::size_t i, std::size_t j, class CouplingContext> | ||
| 376 | ✗ | static Scalar massCouplingCondition(Dune::index_constant<i> domainI, | |
| 377 | Dune::index_constant<j> domainJ, | ||
| 378 | const FVElementGeometry<i>& fvGeometry, | ||
| 379 | const SubControlVolumeFace<i>& scvf, | ||
| 380 | const ElementVolumeVariables<i>& insideVolVars, | ||
| 381 | const CouplingContext& context) | ||
| 382 | { | ||
| 383 | ✗ | const Scalar normalVelocity = context.velocity * scvf.unitOuterNormal(); | |
| 384 | ✗ | const Scalar darcyDensity = insideVolVars[scvf.insideScvIdx()].density(couplingPhaseIdx(ParentType::porousMediumIndex)); | |
| 385 | ✗ | const Scalar stokesDensity = context.volVars.density(); | |
| 386 | ✗ | const bool insideIsUpstream = normalVelocity > 0.0; | |
| 387 | ✗ | return ParentType::advectiveFlux(darcyDensity, stokesDensity, normalVelocity, insideIsUpstream); | |
| 388 | } | ||
| 389 | |||
| 390 | /*! | ||
| 391 | * \brief Returns the energy flux across the coupling boundary as seen from the Darcy domain. | ||
| 392 | */ | ||
| 393 | template<bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0> | ||
| 394 | Scalar energyCouplingCondition(const Element<ParentType::porousMediumIndex>& element, | ||
| 395 | const FVElementGeometry<ParentType::porousMediumIndex>& fvGeometry, | ||
| 396 | const ElementVolumeVariables<ParentType::porousMediumIndex>& darcyElemVolVars, | ||
| 397 | const SubControlVolumeFace<ParentType::porousMediumIndex>& scvf, | ||
| 398 | const DiffusionCoefficientAveragingType diffCoeffAvgType = DiffusionCoefficientAveragingType::ffOnly) const | ||
| 399 | { | ||
| 400 | DUNE_THROW(Dune::NotImplemented, "Energy coupling condition"); | ||
| 401 | } | ||
| 402 | |||
| 403 | /*! | ||
| 404 | * \brief Returns the energy flux across the coupling boundary as seen from the free-flow domain. | ||
| 405 | */ | ||
| 406 | template<bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0> | ||
| 407 | Scalar energyCouplingCondition(const Element<ParentType::freeFlowMassIndex>& element, | ||
| 408 | const FVElementGeometry<ParentType::freeFlowMassIndex>& fvGeometry, | ||
| 409 | const ElementVolumeVariables<ParentType::freeFlowMassIndex>& stokesElemVolVars, | ||
| 410 | const SubControlVolumeFace<ParentType::freeFlowMassIndex>& scvf, | ||
| 411 | const DiffusionCoefficientAveragingType diffCoeffAvgType = DiffusionCoefficientAveragingType::ffOnly) const | ||
| 412 | { | ||
| 413 | DUNE_THROW(Dune::NotImplemented, "Energy coupling condition"); | ||
| 414 | } | ||
| 415 | |||
| 416 | private: | ||
| 417 | |||
| 418 | |||
| 419 | |||
| 420 | /*! | ||
| 421 | * \brief Evaluate the energy flux across the interface. | ||
| 422 | */ | ||
| 423 | template<std::size_t i, std::size_t j, bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0> | ||
| 424 | Scalar energyFlux_(Dune::index_constant<i> domainI, | ||
| 425 | Dune::index_constant<j> domainJ, | ||
| 426 | const FVElementGeometry<i>& insideFvGeometry, | ||
| 427 | const FVElementGeometry<j>& outsideFvGeometry, | ||
| 428 | const SubControlVolumeFace<i>& scvf, | ||
| 429 | const VolumeVariables<i>& insideVolVars, | ||
| 430 | const VolumeVariables<j>& outsideVolVars, | ||
| 431 | const Scalar velocity, | ||
| 432 | const bool insideIsUpstream, | ||
| 433 | const DiffusionCoefficientAveragingType diffCoeffAvgType) const | ||
| 434 | { | ||
| 435 | Scalar flux(0.0); | ||
| 436 | |||
| 437 | const auto& insideScv = (*scvs(insideFvGeometry).begin()); | ||
| 438 | const auto& outsideScv = (*scvs(outsideFvGeometry).begin()); | ||
| 439 | |||
| 440 | // convective fluxes | ||
| 441 | const Scalar insideTerm = insideVolVars.density(couplingPhaseIdx(domainI)) * insideVolVars.enthalpy(couplingPhaseIdx(domainI)); | ||
| 442 | const Scalar outsideTerm = outsideVolVars.density(couplingPhaseIdx(domainJ)) * outsideVolVars.enthalpy(couplingPhaseIdx(domainJ)); | ||
| 443 | |||
| 444 | flux += this->advectiveFlux(insideTerm, outsideTerm, velocity, insideIsUpstream); | ||
| 445 | |||
| 446 | flux += this->conductiveEnergyFlux_(domainI, domainJ, insideFvGeometry, outsideFvGeometry, scvf, insideScv, outsideScv, insideVolVars, outsideVolVars, diffCoeffAvgType); | ||
| 447 | |||
| 448 | return flux; | ||
| 449 | } | ||
| 450 | |||
| 451 | }; | ||
| 452 | |||
| 453 | } // end namespace Dumux | ||
| 454 | |||
| 455 | #endif | ||
| 456 |