| 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-FileCopyrightText: 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 StokesDarcyCoupling | ||
| 10 | * \copydoc Dumux::StokesDarcyCouplingData | ||
| 11 | */ | ||
| 12 | |||
| 13 | #ifndef DUMUX_STOKES_DARCY_COUPLINGDATA_HH | ||
| 14 | #define DUMUX_STOKES_DARCY_COUPLINGDATA_HH | ||
| 15 | |||
| 16 | #include <numeric> | ||
| 17 | |||
| 18 | #include <dumux/common/properties.hh> | ||
| 19 | #include <dumux/common/math.hh> | ||
| 20 | #include <dumux/discretization/method.hh> | ||
| 21 | #include <dumux/discretization/cellcentered/tpfa/computetransmissibility.hh> | ||
| 22 | #include <dumux/flux/referencesystemformulation.hh> | ||
| 23 | #include <dumux/multidomain/couplingmanager.hh> | ||
| 24 | |||
| 25 | namespace Dumux { | ||
| 26 | |||
| 27 | /*! | ||
| 28 | * \ingroup StokesDarcyCoupling | ||
| 29 | * \brief This structs holds a set of options which allow to modify the Stokes-Darcy | ||
| 30 | * coupling mechanism during runtime. | ||
| 31 | */ | ||
| 32 | struct StokesDarcyCouplingOptions | ||
| 33 | { | ||
| 34 | /*! | ||
| 35 | * \brief Defines which kind of averanging of diffusion coefficiencients | ||
| 36 | * (moleculat diffusion or thermal conductance) at the interface | ||
| 37 | * between free flow and porous medium shall be used. | ||
| 38 | */ | ||
| 39 | enum class DiffusionCoefficientAveragingType | ||
| 40 | { | ||
| 41 | harmonic, arithmetic, ffOnly, pmOnly | ||
| 42 | }; | ||
| 43 | |||
| 44 | /*! | ||
| 45 | * \brief Convenience function to convert user input given as std::string to the corresponding enum class used for choosing the type | ||
| 46 | * of averaging of the diffusion/conduction parameter at the interface between the two domains. | ||
| 47 | */ | ||
| 48 | 4 | static DiffusionCoefficientAveragingType stringToEnum(DiffusionCoefficientAveragingType, const std::string& diffusionCoefficientAveragingType) | |
| 49 | { | ||
| 50 |
1/2✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
|
4 | if (diffusionCoefficientAveragingType == "Harmonic") |
| 51 | return DiffusionCoefficientAveragingType::harmonic; | ||
| 52 |
1/2✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
|
4 | else if (diffusionCoefficientAveragingType == "Arithmetic") |
| 53 | return DiffusionCoefficientAveragingType::arithmetic; | ||
| 54 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
|
4 | else if (diffusionCoefficientAveragingType == "FreeFlowOnly") |
| 55 | return DiffusionCoefficientAveragingType::ffOnly; | ||
| 56 | ✗ | else if (diffusionCoefficientAveragingType == "PorousMediumOnly") | |
| 57 | return DiffusionCoefficientAveragingType::pmOnly; | ||
| 58 | else | ||
| 59 | ✗ | DUNE_THROW(Dune::IOError, "Unknown DiffusionCoefficientAveragingType"); | |
| 60 | } | ||
| 61 | |||
| 62 | }; | ||
| 63 | |||
| 64 | /*! | ||
| 65 | * \ingroup StokesDarcyCoupling | ||
| 66 | * \brief This structs helps to check if the two sub models use the same fluidsystem. | ||
| 67 | * Specialization for the case of using an adapter only for the free-flow model. | ||
| 68 | * \tparam FFFS The free-flow fluidsystem | ||
| 69 | * \tparam PMFS The porous-medium flow fluidsystem | ||
| 70 | */ | ||
| 71 | template<class FFFS, class PMFS> | ||
| 72 | struct IsSameFluidSystem | ||
| 73 | { | ||
| 74 | static_assert(FFFS::numPhases == 1, "Only single-phase fluidsystems may be used for free flow."); | ||
| 75 | static constexpr bool value = std::is_same<typename FFFS::MultiPhaseFluidSystem, PMFS>::value; | ||
| 76 | }; | ||
| 77 | |||
| 78 | /*! | ||
| 79 | * \ingroup StokesDarcyCoupling | ||
| 80 | * \brief This structs helps to check if the two sub models use the same fluidsystem. | ||
| 81 | * \tparam FS The fluidsystem | ||
| 82 | */ | ||
| 83 | template<class FS> | ||
| 84 | struct IsSameFluidSystem<FS, FS> | ||
| 85 | { | ||
| 86 | static_assert(FS::numPhases == 1, "Only single-phase fluidsystems may be used for free flow."); | ||
| 87 | static constexpr bool value = std::is_same<FS, FS>::value; // always true | ||
| 88 | }; | ||
| 89 | |||
| 90 | // forward declaration | ||
| 91 | template <class TypeTag, class DiscretizationMethod, ReferenceSystemFormulation referenceSystem> | ||
| 92 | class FicksLawImplementation; | ||
| 93 | |||
| 94 | /*! | ||
| 95 | * \ingroup StokesDarcyCoupling | ||
| 96 | * \brief This structs indicates that Fick's law is not used for diffusion. | ||
| 97 | * \tparam DiffLaw The diffusion law. | ||
| 98 | */ | ||
| 99 | template<class DiffLaw> | ||
| 100 | struct IsFicksLaw : public std::false_type {}; | ||
| 101 | |||
| 102 | /*! | ||
| 103 | * \ingroup StokesDarcyCoupling | ||
| 104 | * \brief This structs indicates that Fick's law is used for diffusion. | ||
| 105 | * \tparam DiffLaw The diffusion law. | ||
| 106 | */ | ||
| 107 | template<class T, class DiscretizationMethod, ReferenceSystemFormulation referenceSystem> | ||
| 108 | struct IsFicksLaw<FicksLawImplementation<T, DiscretizationMethod, referenceSystem>> : public std::true_type {}; | ||
| 109 | |||
| 110 | /*! | ||
| 111 | * \ingroup StokesDarcyCoupling | ||
| 112 | * \brief Helper struct to choose the correct index for phases and components. This is need if the porous-medium-flow model | ||
| 113 | features more fluid phases than the free-flow model. | ||
| 114 | * \tparam stokesIdx The domain index of the free-flow model. | ||
| 115 | * \tparam darcyIdx The domain index of the porous-medium-flow model. | ||
| 116 | * \tparam FFFS The free-flow fluidsystem. | ||
| 117 | * \tparam hasAdapter Specifies whether an adapter class for the fluidsystem is used. | ||
| 118 | */ | ||
| 119 | template<std::size_t stokesIdx, std::size_t darcyIdx, class FFFS, bool hasAdapter> | ||
| 120 | struct IndexHelper; | ||
| 121 | |||
| 122 | /*! | ||
| 123 | * \ingroup StokesDarcyCoupling | ||
| 124 | * \brief Helper struct to choose the correct index for phases and components. This is need if the porous-medium-flow model | ||
| 125 | features more fluid phases than the free-flow model. Specialization for the case that no adapter is used. | ||
| 126 | * \tparam stokesIdx The domain index of the free-flow model. | ||
| 127 | * \tparam darcyIdx The domain index of the porous-medium-flow model. | ||
| 128 | * \tparam FFFS The free-flow fluidsystem. | ||
| 129 | */ | ||
| 130 | template<std::size_t stokesIdx, std::size_t darcyIdx, class FFFS> | ||
| 131 | struct IndexHelper<stokesIdx, darcyIdx, FFFS, false> | ||
| 132 | { | ||
| 133 | /*! | ||
| 134 | * \brief No adapter is used, just return the input index. | ||
| 135 | */ | ||
| 136 | template<std::size_t i> | ||
| 137 | static constexpr auto couplingPhaseIdx(Dune::index_constant<i>, int coupledPhaseIdx = 0) | ||
| 138 | { return coupledPhaseIdx; } | ||
| 139 | |||
| 140 | /*! | ||
| 141 | * \brief No adapter is used, just return the input index. | ||
| 142 | */ | ||
| 143 | template<std::size_t i> | ||
| 144 | static constexpr auto couplingCompIdx(Dune::index_constant<i>, int coupledCompdIdx) | ||
| 145 | { return coupledCompdIdx; } | ||
| 146 | }; | ||
| 147 | |||
| 148 | /*! | ||
| 149 | * \ingroup StokesDarcyCoupling | ||
| 150 | * \brief Helper struct to choose the correct index for phases and components. This is need if the porous-medium-flow model | ||
| 151 | features more fluid phases than the free-flow model. Specialization for the case that a adapter is used. | ||
| 152 | * \tparam stokesIdx The domain index of the free-flow model. | ||
| 153 | * \tparam darcyIdx The domain index of the porous-medium-flow model. | ||
| 154 | * \tparam FFFS The free-flow fluidsystem. | ||
| 155 | */ | ||
| 156 | template<std::size_t stokesIdx, std::size_t darcyIdx, class FFFS> | ||
| 157 | struct IndexHelper<stokesIdx, darcyIdx, FFFS, true> | ||
| 158 | { | ||
| 159 | /*! | ||
| 160 | * \brief The free-flow model always uses phase index 0. | ||
| 161 | */ | ||
| 162 | static constexpr auto couplingPhaseIdx(Dune::index_constant<stokesIdx>, int coupledPhaseIdx = 0) | ||
| 163 | { return 0; } | ||
| 164 | |||
| 165 | /*! | ||
| 166 | * \brief The phase index of the porous-medium-flow model is given by the adapter fluidsytem (i.e., user input). | ||
| 167 | */ | ||
| 168 | static constexpr auto couplingPhaseIdx(Dune::index_constant<darcyIdx>, int coupledPhaseIdx = 0) | ||
| 169 | { return FFFS::multiphaseFluidsystemPhaseIdx; } | ||
| 170 | |||
| 171 | /*! | ||
| 172 | * \brief The free-flow model does not need any change of the component index. | ||
| 173 | */ | ||
| 174 | static constexpr auto couplingCompIdx(Dune::index_constant<stokesIdx>, int coupledCompdIdx) | ||
| 175 | { return coupledCompdIdx; } | ||
| 176 | |||
| 177 | /*! | ||
| 178 | * \brief The component index of the porous-medium-flow model is mapped by the adapter fluidsytem. | ||
| 179 | */ | ||
| 180 | 548454 | static constexpr auto couplingCompIdx(Dune::index_constant<darcyIdx>, int coupledCompdIdx) | |
| 181 | 548454 | { return FFFS::compIdx(coupledCompdIdx); } | |
| 182 | }; | ||
| 183 | |||
| 184 | //! forward declare | ||
| 185 | template <class TypeTag, class DiscretizationMethod> | ||
| 186 | class DarcysLawImplementation; | ||
| 187 | |||
| 188 | //! forward declare | ||
| 189 | template <class TypeTag, class DiscretizationMethod> | ||
| 190 | class ForchheimersLawImplementation; | ||
| 191 | |||
| 192 | |||
| 193 | template<class MDTraits, class CouplingManager, bool enableEnergyBalance, bool isCompositional> | ||
| 194 | class StokesDarcyCouplingDataImplementation; | ||
| 195 | |||
| 196 | /*! | ||
| 197 | * \ingroup StokesDarcyCoupling | ||
| 198 | * \brief Data for the coupling of a Darcy model (cell-centered finite volume) | ||
| 199 | * with a (Navier-)Stokes model (staggerd grid). | ||
| 200 | */ | ||
| 201 | template<class MDTraits, class CouplingManager> | ||
| 202 | using StokesDarcyCouplingData = StokesDarcyCouplingDataImplementation<MDTraits, CouplingManager, | ||
| 203 | GetPropType<typename MDTraits::template SubDomain<0>::TypeTag, Properties::ModelTraits>::enableEnergyBalance(), | ||
| 204 | (GetPropType<typename MDTraits::template SubDomain<0>::TypeTag, Properties::ModelTraits>::numFluidComponents() > 1)>; | ||
| 205 | |||
| 206 | /*! | ||
| 207 | * \ingroup StokesDarcyCoupling | ||
| 208 | * \brief A base class which provides some common methods used for Stokes-Darcy coupling. | ||
| 209 | */ | ||
| 210 | template<class MDTraits, class CouplingManager> | ||
| 211 | class StokesDarcyCouplingDataImplementationBase | ||
| 212 | { | ||
| 213 | using Scalar = typename MDTraits::Scalar; | ||
| 214 | |||
| 215 | template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag; | ||
| 216 | template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>; | ||
| 217 | template<std::size_t id> using Element = typename GridGeometry<id>::GridView::template Codim<0>::Entity; | ||
| 218 | template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView; | ||
| 219 | template<std::size_t id> using SubControlVolumeFace = typename GridGeometry<id>::LocalView::SubControlVolumeFace; | ||
| 220 | template<std::size_t id> using SubControlVolume = typename GridGeometry<id>::LocalView::SubControlVolume; | ||
| 221 | template<std::size_t id> using Indices = typename GetPropType<SubDomainTypeTag<id>, Properties::ModelTraits>::Indices; | ||
| 222 | template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView; | ||
| 223 | template<std::size_t id> using VolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::VolumeVariables; | ||
| 224 | template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>; | ||
| 225 | template<std::size_t id> using FluidSystem = GetPropType<SubDomainTypeTag<id>, Properties::FluidSystem>; | ||
| 226 | template<std::size_t id> using ModelTraits = GetPropType<SubDomainTypeTag<id>, Properties::ModelTraits>; | ||
| 227 | template<std::size_t id> using GlobalPosition = typename Element<id>::Geometry::GlobalCoordinate; | ||
| 228 | static constexpr auto stokesIdx = CouplingManager::stokesIdx; | ||
| 229 | static constexpr auto darcyIdx = CouplingManager::darcyIdx; | ||
| 230 | |||
| 231 | using AdvectionType = GetPropType<SubDomainTypeTag<darcyIdx>, Properties::AdvectionType>; | ||
| 232 | using DarcysLaw = DarcysLawImplementation<SubDomainTypeTag<darcyIdx>, typename GridGeometry<darcyIdx>::DiscretizationMethod>; | ||
| 233 | using ForchheimersLaw = ForchheimersLawImplementation<SubDomainTypeTag<darcyIdx>, typename GridGeometry<darcyIdx>::DiscretizationMethod>; | ||
| 234 | |||
| 235 | static constexpr bool adapterUsed = ModelTraits<darcyIdx>::numFluidPhases() > 1; | ||
| 236 | using IndexHelper = Dumux::IndexHelper<stokesIdx, darcyIdx, FluidSystem<stokesIdx>, adapterUsed>; | ||
| 237 | |||
| 238 | static constexpr int enableEnergyBalance = GetPropType<SubDomainTypeTag<stokesIdx>, Properties::ModelTraits>::enableEnergyBalance(); | ||
| 239 | static_assert(GetPropType<SubDomainTypeTag<darcyIdx>, Properties::ModelTraits>::enableEnergyBalance() == enableEnergyBalance, | ||
| 240 | "All submodels must both be either isothermal or non-isothermal"); | ||
| 241 | |||
| 242 | static_assert(IsSameFluidSystem<FluidSystem<stokesIdx>, | ||
| 243 | FluidSystem<darcyIdx>>::value, | ||
| 244 | "All submodels must use the same fluid system"); | ||
| 245 | |||
| 246 | using DiffusionCoefficientAveragingType = typename StokesDarcyCouplingOptions::DiffusionCoefficientAveragingType; | ||
| 247 | |||
| 248 | public: | ||
| 249 | 17 | StokesDarcyCouplingDataImplementationBase(const CouplingManager& couplingmanager): couplingManager_(couplingmanager) {} | |
| 250 | |||
| 251 | /*! | ||
| 252 | * \brief Returns the corresponding phase index needed for coupling. | ||
| 253 | */ | ||
| 254 | template<std::size_t i> | ||
| 255 | static constexpr auto couplingPhaseIdx(Dune::index_constant<i> id, int coupledPhaseIdx = 0) | ||
| 256 | { return IndexHelper::couplingPhaseIdx(id, coupledPhaseIdx); } | ||
| 257 | |||
| 258 | /*! | ||
| 259 | * \brief Returns the corresponding component index needed for coupling. | ||
| 260 | */ | ||
| 261 | template<std::size_t i> | ||
| 262 |
6/6✓ Branch 0 taken 21984 times.
✓ Branch 1 taken 21984 times.
✓ Branch 2 taken 21984 times.
✓ Branch 3 taken 21984 times.
✓ Branch 4 taken 12320 times.
✓ Branch 5 taken 12320 times.
|
112576 | static constexpr auto couplingCompIdx(Dune::index_constant<i> id, int coupledCompdIdx) |
| 263 | 548454 | { return IndexHelper::couplingCompIdx(id, coupledCompdIdx); } | |
| 264 | |||
| 265 | /*! | ||
| 266 | * \brief Returns a reference to the coupling manager. | ||
| 267 | */ | ||
| 268 | 1364080 | const CouplingManager& couplingManager() const | |
| 269 | 1364080 | { return couplingManager_; } | |
| 270 | |||
| 271 | /*! | ||
| 272 | * \brief Returns the intrinsic permeability of the coupled Darcy element. | ||
| 273 | */ | ||
| 274 | 863710 | auto darcyPermeability(const Element<stokesIdx>& element, const SubControlVolumeFace<stokesIdx>& scvf) const | |
| 275 | { | ||
| 276 | 863710 | const auto& stokesContext = couplingManager().stokesCouplingContext(element, scvf); | |
| 277 | 863710 | return stokesContext.volVars.permeability(); | |
| 278 | } | ||
| 279 | |||
| 280 | /*! | ||
| 281 | * \brief Returns the momentum flux across the coupling boundary. | ||
| 282 | * | ||
| 283 | * For the normal momentum coupling, the porous medium side of the coupling condition | ||
| 284 | * is evaluated, i.e. -[p n]^pm. | ||
| 285 | * | ||
| 286 | */ | ||
| 287 | template<class ElementFaceVariables> | ||
| 288 | 346770 | Scalar momentumCouplingCondition(const Element<stokesIdx>& element, | |
| 289 | const FVElementGeometry<stokesIdx>& fvGeometry, | ||
| 290 | const ElementVolumeVariables<stokesIdx>& stokesElemVolVars, | ||
| 291 | const ElementFaceVariables& stokesElemFaceVars, | ||
| 292 | const SubControlVolumeFace<stokesIdx>& scvf) const | ||
| 293 | { | ||
| 294 | static constexpr auto numPhasesDarcy = GetPropType<SubDomainTypeTag<darcyIdx>, Properties::ModelTraits>::numFluidPhases(); | ||
| 295 | |||
| 296 | 346770 | Scalar momentumFlux(0.0); | |
| 297 | 346770 | const auto& stokesContext = couplingManager_.stokesCouplingContext(element, scvf); | |
| 298 | 346770 | const auto darcyPhaseIdx = couplingPhaseIdx(darcyIdx); | |
| 299 | |||
| 300 | // - p_pm * n_pm = p_pm * n_ff | ||
| 301 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 124528 times.
|
346770 | const Scalar darcyPressure = stokesContext.volVars.pressure(darcyPhaseIdx); |
| 302 | |||
| 303 | if(numPhasesDarcy > 1) | ||
| 304 | 124528 | momentumFlux = darcyPressure; | |
| 305 | else // use pressure reconstruction for single phase models | ||
| 306 | 222242 | momentumFlux = pressureAtInterface_(element, scvf, stokesElemFaceVars, stokesContext); | |
| 307 | // TODO: generalize for permeability tensors | ||
| 308 | |||
| 309 | // normalize pressure | ||
| 310 | if(getPropValue<SubDomainTypeTag<stokesIdx>, Properties::NormalizePressure>()) | ||
| 311 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 346770 times.
|
346770 | momentumFlux -= couplingManager_.problem(stokesIdx).initial(scvf)[Indices<stokesIdx>::pressureIdx]; |
| 312 | |||
| 313 | 346770 | momentumFlux *= scvf.directionSign(); | |
| 314 | |||
| 315 | 346770 | return momentumFlux; | |
| 316 | } | ||
| 317 | |||
| 318 | /*! | ||
| 319 | * \brief Evaluate an advective flux across the interface and consider upwinding. | ||
| 320 | */ | ||
| 321 | 881080 | Scalar advectiveFlux(const Scalar insideQuantity, const Scalar outsideQuantity, const Scalar volumeFlow, bool insideIsUpstream) const | |
| 322 | { | ||
| 323 | 881080 | const Scalar upwindWeight = 1.0; //TODO use Flux.UpwindWeight or something like Coupling.UpwindWeight? | |
| 324 | |||
| 325 |
3/4✓ Branch 0 taken 43712 times.
✓ Branch 1 taken 256 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 171004 times.
|
820404 | if(insideIsUpstream) |
| 326 | 399208 | return (upwindWeight * insideQuantity + (1.0 - upwindWeight) * outsideQuantity) * volumeFlow; | |
| 327 | else | ||
| 328 | 481872 | return (upwindWeight * outsideQuantity + (1.0 - upwindWeight) * insideQuantity) * volumeFlow; | |
| 329 | } | ||
| 330 | |||
| 331 | protected: | ||
| 332 | |||
| 333 | /*! | ||
| 334 | * \brief Returns the transmissibility used for either molecular diffusion or thermal conductivity. | ||
| 335 | */ | ||
| 336 | template<std::size_t i, std::size_t j> | ||
| 337 | 731676 | Scalar transmissibility_(Dune::index_constant<i> domainI, | |
| 338 | Dune::index_constant<j> domainJ, | ||
| 339 | const Scalar insideDistance, | ||
| 340 | const Scalar outsideDistance, | ||
| 341 | const Scalar avgQuantityI, | ||
| 342 | const Scalar avgQuantityJ, | ||
| 343 | const DiffusionCoefficientAveragingType diffCoeffAvgType) const | ||
| 344 | { | ||
| 345 | 731676 | const Scalar totalDistance = insideDistance + outsideDistance; | |
| 346 |
2/2✓ Branch 0 taken 31600 times.
✓ Branch 1 taken 334238 times.
|
731676 | if(diffCoeffAvgType == DiffusionCoefficientAveragingType::harmonic) |
| 347 | { | ||
| 348 | 63200 | return harmonicMean(avgQuantityI, avgQuantityJ, insideDistance, outsideDistance) | |
| 349 | 63200 | / totalDistance; | |
| 350 | } | ||
| 351 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 334238 times.
|
668476 | else if(diffCoeffAvgType == DiffusionCoefficientAveragingType::arithmetic) |
| 352 | { | ||
| 353 | ✗ | return arithmeticMean(avgQuantityI, avgQuantityJ, insideDistance, outsideDistance) | |
| 354 | ✗ | / totalDistance; | |
| 355 | } | ||
| 356 |
1/2✓ Branch 0 taken 334238 times.
✗ Branch 1 not taken.
|
668476 | else if(diffCoeffAvgType == DiffusionCoefficientAveragingType::ffOnly) |
| 357 | return domainI == stokesIdx | ||
| 358 | ? avgQuantityI / totalDistance | ||
| 359 | 668476 | : avgQuantityJ / totalDistance; | |
| 360 | |||
| 361 | else // diffCoeffAvgType == DiffusionCoefficientAveragingType::pmOnly) | ||
| 362 | return domainI == darcyIdx | ||
| 363 | ? avgQuantityI / totalDistance | ||
| 364 | ✗ | : avgQuantityJ / totalDistance; | |
| 365 | } | ||
| 366 | |||
| 367 | /*! | ||
| 368 | * \brief Returns the distance between an scvf and the corresponding scv center. | ||
| 369 | */ | ||
| 370 | template<class Scv, class Scvf> | ||
| 371 | 1799000 | Scalar getDistance_(const Scv& scv, const Scvf& scvf) const | |
| 372 | { | ||
| 373 | 3598000 | return (scv.dofPosition() - scvf.ipGlobal()).two_norm(); | |
| 374 | } | ||
| 375 | |||
| 376 | /*! | ||
| 377 | * \brief Returns the conductive energy flux across the interface. | ||
| 378 | */ | ||
| 379 | template<std::size_t i, std::size_t j, bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0> | ||
| 380 | 121352 | Scalar conductiveEnergyFlux_(Dune::index_constant<i> domainI, | |
| 381 | Dune::index_constant<j> domainJ, | ||
| 382 | const FVElementGeometry<i>& fvGeometryI, | ||
| 383 | const FVElementGeometry<j>& fvGeometryJ, | ||
| 384 | const SubControlVolumeFace<i>& scvfI, | ||
| 385 | const SubControlVolume<i>& scvI, | ||
| 386 | const SubControlVolume<j>& scvJ, | ||
| 387 | const VolumeVariables<i>& volVarsI, | ||
| 388 | const VolumeVariables<j>& volVarsJ, | ||
| 389 | const DiffusionCoefficientAveragingType diffCoeffAvgType) const | ||
| 390 | { | ||
| 391 | 121352 | const Scalar insideDistance = getDistance_(scvI, scvfI); | |
| 392 | 121352 | const Scalar outsideDistance = getDistance_(scvJ, scvfI); | |
| 393 | |||
| 394 | 121352 | const Scalar deltaT = volVarsJ.temperature() - volVarsI.temperature(); | |
| 395 | 121352 | const Scalar tij = transmissibility_(domainI, | |
| 396 | domainJ, | ||
| 397 | insideDistance, | ||
| 398 | outsideDistance, | ||
| 399 | volVarsI.effectiveThermalConductivity(), | ||
| 400 | volVarsJ.effectiveThermalConductivity(), | ||
| 401 | diffCoeffAvgType); | ||
| 402 | |||
| 403 | 121352 | return -tij * deltaT; | |
| 404 | } | ||
| 405 | |||
| 406 | /*! | ||
| 407 | * \brief Returns the pressure at the interface | ||
| 408 | */ | ||
| 409 | template<class ElementFaceVariables, class CouplingContext> | ||
| 410 | 222242 | Scalar pressureAtInterface_(const Element<stokesIdx>& element, | |
| 411 | const SubControlVolumeFace<stokesIdx>& scvf, | ||
| 412 | const ElementFaceVariables& elemFaceVars, | ||
| 413 | const CouplingContext& context) const | ||
| 414 | { | ||
| 415 | 222242 | GlobalPosition<stokesIdx> velocity(0.0); | |
| 416 | 222242 | velocity[scvf.directionIndex()] = elemFaceVars[scvf].velocitySelf(); | |
| 417 | 222242 | const auto& darcyScvf = context.fvGeometry.scvf(context.darcyScvfIdx); | |
| 418 | 222242 | return computeCouplingPhasePressureAtInterface_(context.element, context.fvGeometry, darcyScvf, context.volVars, velocity, AdvectionType()); | |
| 419 | } | ||
| 420 | |||
| 421 | /*! | ||
| 422 | * \brief Returns the pressure at the interface using Forchheimers's law for reconstruction | ||
| 423 | */ | ||
| 424 | Scalar computeCouplingPhasePressureAtInterface_(const Element<darcyIdx>& element, | ||
| 425 | const FVElementGeometry<darcyIdx>& fvGeometry, | ||
| 426 | const SubControlVolumeFace<darcyIdx>& scvf, | ||
| 427 | const VolumeVariables<darcyIdx>& volVars, | ||
| 428 | const typename Element<stokesIdx>::Geometry::GlobalCoordinate& couplingPhaseVelocity, | ||
| 429 | ForchheimersLaw) const | ||
| 430 | { | ||
| 431 | const auto darcyPhaseIdx = couplingPhaseIdx(darcyIdx); | ||
| 432 | const Scalar cellCenterPressure = volVars.pressure(darcyPhaseIdx); | ||
| 433 | using std::sqrt; | ||
| 434 | |||
| 435 | // v + (cF*sqrt(K)*rho/mu*|v|) * v = - K/mu grad(p - rho g) | ||
| 436 | // multiplying with n and using a tpfa for the right-hand side yields | ||
| 437 | // v*n + (cF*sqrt(K)*rho/mu*|v|) * (v*n) = 1/mu * (ti*(p_center - p_interface) + rho*n^TKg) | ||
| 438 | // --> p_interface = (-mu*v*n + (cF*sqrt(K)*rho*|v|) * (-v*n) + rho*n^TKg)/ti + p_center | ||
| 439 | const auto velocity = couplingPhaseVelocity; | ||
| 440 | const Scalar mu = volVars.viscosity(darcyPhaseIdx); | ||
| 441 | const Scalar rho = volVars.density(darcyPhaseIdx); | ||
| 442 | const auto K = volVars.permeability(); | ||
| 443 | const auto alpha = vtmv(scvf.unitOuterNormal(), K, couplingManager_.problem(darcyIdx).spatialParams().gravity(scvf.center())); | ||
| 444 | |||
| 445 | const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); | ||
| 446 | const auto ti = computeTpfaTransmissibility(scvf, insideScv, K, 1.0); | ||
| 447 | |||
| 448 | // get the Forchheimer coefficient | ||
| 449 | Scalar cF = couplingManager_.problem(darcyIdx).spatialParams().forchCoeff(scvf); | ||
| 450 | |||
| 451 | const Scalar interfacePressure = ((-mu*(scvf.unitOuterNormal() * velocity)) | ||
| 452 | + (-(scvf.unitOuterNormal() * velocity) * velocity.two_norm() * rho * sqrt(darcyPermeability(element, scvf)) * cF) | ||
| 453 | + rho * alpha)/ti + cellCenterPressure; | ||
| 454 | return interfacePressure; | ||
| 455 | } | ||
| 456 | |||
| 457 | /*! | ||
| 458 | * \brief Returns the pressure at the interface using Darcy's law for reconstruction | ||
| 459 | */ | ||
| 460 | 222242 | Scalar computeCouplingPhasePressureAtInterface_(const Element<darcyIdx>& element, | |
| 461 | const FVElementGeometry<darcyIdx>& fvGeometry, | ||
| 462 | const SubControlVolumeFace<darcyIdx>& scvf, | ||
| 463 | const VolumeVariables<darcyIdx>& volVars, | ||
| 464 | const typename Element<stokesIdx>::Geometry::GlobalCoordinate& couplingPhaseVelocity, | ||
| 465 | DarcysLaw) const | ||
| 466 | { | ||
| 467 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 222242 times.
|
222242 | const auto darcyPhaseIdx = couplingPhaseIdx(darcyIdx); |
| 468 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 222242 times.
|
222242 | const Scalar couplingPhaseCellCenterPressure = volVars.pressure(darcyPhaseIdx); |
| 469 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 222242 times.
|
222242 | const Scalar couplingPhaseMobility = volVars.mobility(darcyPhaseIdx); |
| 470 | 222242 | const Scalar couplingPhaseDensity = volVars.density(darcyPhaseIdx); | |
| 471 | 222242 | const auto K = volVars.permeability(); | |
| 472 | |||
| 473 | // A tpfa approximation yields (works if mobility != 0) | ||
| 474 | // v*n = -kr/mu*K * (gradP - rho*g)*n = mobility*(ti*(p_center - p_interface) + rho*n^TKg) | ||
| 475 | // -> p_interface = (1/mobility * (-v*n) + rho*n^TKg)/ti + p_center | ||
| 476 | // where v is the free-flow velocity (couplingPhaseVelocity) | ||
| 477 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 222242 times.
✓ Branch 2 taken 222242 times.
✗ Branch 3 not taken.
|
399530 | const auto alpha = vtmv(scvf.unitOuterNormal(), K, couplingManager_.problem(darcyIdx).spatialParams().gravity(scvf.center())); |
| 478 | |||
| 479 |
1/2✓ Branch 0 taken 222242 times.
✗ Branch 1 not taken.
|
222242 | const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); |
| 480 | 222242 | const auto ti = computeTpfaTransmissibility(fvGeometry, scvf, insideScv, K, 1.0); | |
| 481 | |||
| 482 | 222242 | return (-1/couplingPhaseMobility * (scvf.unitOuterNormal() * couplingPhaseVelocity) + couplingPhaseDensity * alpha)/ti | |
| 483 | 222242 | + couplingPhaseCellCenterPressure; | |
| 484 | } | ||
| 485 | |||
| 486 | |||
| 487 | private: | ||
| 488 | const CouplingManager& couplingManager_; | ||
| 489 | |||
| 490 | }; | ||
| 491 | |||
| 492 | /*! | ||
| 493 | * \ingroup StokesDarcyCoupling | ||
| 494 | * \brief Coupling data specialization for non-compositional models. | ||
| 495 | */ | ||
| 496 | template<class MDTraits, class CouplingManager, bool enableEnergyBalance> | ||
| 497 | class StokesDarcyCouplingDataImplementation<MDTraits, CouplingManager, enableEnergyBalance, false> | ||
| 498 | : public StokesDarcyCouplingDataImplementationBase<MDTraits, CouplingManager> | ||
| 499 | { | ||
| 500 | using ParentType = StokesDarcyCouplingDataImplementationBase<MDTraits, CouplingManager>; | ||
| 501 | using Scalar = typename MDTraits::Scalar; | ||
| 502 | static constexpr auto stokesIdx = CouplingManager::stokesIdx; | ||
| 503 | static constexpr auto darcyIdx = CouplingManager::darcyIdx; | ||
| 504 | static constexpr auto stokesFaceIdx = CouplingManager::stokesFaceIdx; | ||
| 505 | static constexpr auto stokesCellCenterIdx = CouplingManager::stokesCellCenterIdx; | ||
| 506 | |||
| 507 | // the sub domain type tags | ||
| 508 | template<std::size_t id> | ||
| 509 | using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag; | ||
| 510 | |||
| 511 | template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>; | ||
| 512 | template<std::size_t id> using Element = typename GridGeometry<id>::GridView::template Codim<0>::Entity; | ||
| 513 | template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView; | ||
| 514 | template<std::size_t id> using SubControlVolumeFace = typename GridGeometry<id>::LocalView::SubControlVolumeFace; | ||
| 515 | template<std::size_t id> using SubControlVolume = typename GridGeometry<id>::LocalView::SubControlVolume; | ||
| 516 | template<std::size_t id> using Indices = typename GetPropType<SubDomainTypeTag<id>, Properties::ModelTraits>::Indices; | ||
| 517 | template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView; | ||
| 518 | template<std::size_t id> using ElementFaceVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridFaceVariables>::LocalView; | ||
| 519 | template<std::size_t id> using VolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::VolumeVariables; | ||
| 520 | |||
| 521 | static_assert(GetPropType<SubDomainTypeTag<darcyIdx>, Properties::ModelTraits>::numFluidComponents() == GetPropType<SubDomainTypeTag<darcyIdx>, Properties::ModelTraits>::numFluidPhases(), | ||
| 522 | "Darcy Model must not be compositional"); | ||
| 523 | |||
| 524 | using DiffusionCoefficientAveragingType = typename StokesDarcyCouplingOptions::DiffusionCoefficientAveragingType; | ||
| 525 | |||
| 526 | public: | ||
| 527 | 9 | using ParentType::ParentType; | |
| 528 | using ParentType::couplingPhaseIdx; | ||
| 529 | |||
| 530 | /*! | ||
| 531 | * \brief Returns the mass flux across the coupling boundary as seen from the Darcy domain. | ||
| 532 | */ | ||
| 533 | 25500 | Scalar massCouplingCondition(const Element<darcyIdx>& element, | |
| 534 | const FVElementGeometry<darcyIdx>& fvGeometry, | ||
| 535 | const ElementVolumeVariables<darcyIdx>& darcyElemVolVars, | ||
| 536 | const SubControlVolumeFace<darcyIdx>& scvf) const | ||
| 537 | { | ||
| 538 | 25500 | const auto& darcyContext = this->couplingManager().darcyCouplingContext(element, scvf); | |
| 539 | 51000 | const Scalar velocity = darcyContext.velocity * scvf.unitOuterNormal(); | |
| 540 |
2/2✓ Branch 1 taken 700 times.
✓ Branch 2 taken 24800 times.
|
25500 | const Scalar darcyDensity = darcyElemVolVars[scvf.insideScvIdx()].density(couplingPhaseIdx(darcyIdx)); |
| 541 |
2/2✓ Branch 0 taken 700 times.
✓ Branch 1 taken 24800 times.
|
25500 | const Scalar stokesDensity = darcyContext.volVars.density(); |
| 542 | 25500 | const bool insideIsUpstream = velocity > 0.0; | |
| 543 | |||
| 544 | 51000 | return massFlux_(velocity, darcyDensity, stokesDensity, insideIsUpstream); | |
| 545 | } | ||
| 546 | |||
| 547 | /*! | ||
| 548 | * \brief Returns the mass flux across the coupling boundary as seen from the free-flow domain. | ||
| 549 | */ | ||
| 550 | 85796 | Scalar massCouplingCondition(const Element<stokesIdx>& element, | |
| 551 | const FVElementGeometry<stokesIdx>& fvGeometry, | ||
| 552 | const ElementVolumeVariables<stokesIdx>& stokesElemVolVars, | ||
| 553 | const ElementFaceVariables<stokesIdx>& stokesElemFaceVars, | ||
| 554 | const SubControlVolumeFace<stokesIdx>& scvf) const | ||
| 555 | { | ||
| 556 | 85796 | const auto& stokesContext = this->couplingManager().stokesCouplingContext(element, scvf); | |
| 557 | 85796 | const Scalar velocity = stokesElemFaceVars[scvf].velocitySelf(); | |
| 558 |
2/2✓ Branch 1 taken 71757 times.
✓ Branch 2 taken 14039 times.
|
85796 | const Scalar stokesDensity = stokesElemVolVars[scvf.insideScvIdx()].density(); |
| 559 |
2/2✓ Branch 0 taken 71757 times.
✓ Branch 1 taken 14039 times.
|
85796 | const Scalar darcyDensity = stokesContext.volVars.density(couplingPhaseIdx(darcyIdx)); |
| 560 |
2/2✓ Branch 0 taken 71757 times.
✓ Branch 1 taken 14039 times.
|
85796 | const bool insideIsUpstream = sign(velocity) == scvf.directionSign(); |
| 561 | |||
| 562 | 171592 | return massFlux_(velocity * scvf.directionSign(), stokesDensity, darcyDensity, insideIsUpstream); | |
| 563 | } | ||
| 564 | |||
| 565 | /*! | ||
| 566 | * \brief Returns the energy flux across the coupling boundary as seen from the Darcy domain. | ||
| 567 | */ | ||
| 568 | template<bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0> | ||
| 569 | Scalar energyCouplingCondition(const Element<darcyIdx>& element, | ||
| 570 | const FVElementGeometry<darcyIdx>& fvGeometry, | ||
| 571 | const ElementVolumeVariables<darcyIdx>& darcyElemVolVars, | ||
| 572 | const SubControlVolumeFace<darcyIdx>& scvf, | ||
| 573 | const DiffusionCoefficientAveragingType diffCoeffAvgType = DiffusionCoefficientAveragingType::ffOnly) const | ||
| 574 | { | ||
| 575 | const auto& darcyContext = this->couplingManager().darcyCouplingContext(element, scvf); | ||
| 576 | const auto& darcyVolVars = darcyElemVolVars[scvf.insideScvIdx()]; | ||
| 577 | const auto& stokesVolVars = darcyContext.volVars; | ||
| 578 | |||
| 579 | const Scalar velocity = darcyContext.velocity * scvf.unitOuterNormal(); | ||
| 580 | const bool insideIsUpstream = velocity > 0.0; | ||
| 581 | |||
| 582 | return energyFlux_(darcyIdx, stokesIdx, fvGeometry, darcyContext.fvGeometry, scvf, | ||
| 583 | darcyVolVars, stokesVolVars, velocity, insideIsUpstream, diffCoeffAvgType); | ||
| 584 | } | ||
| 585 | |||
| 586 | /*! | ||
| 587 | * \brief Returns the energy flux across the coupling boundary as seen from the free-flow domain. | ||
| 588 | */ | ||
| 589 | template<bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0> | ||
| 590 | Scalar energyCouplingCondition(const Element<stokesIdx>& element, | ||
| 591 | const FVElementGeometry<stokesIdx>& fvGeometry, | ||
| 592 | const ElementVolumeVariables<stokesIdx>& stokesElemVolVars, | ||
| 593 | const ElementFaceVariables<stokesIdx>& stokesElemFaceVars, | ||
| 594 | const SubControlVolumeFace<stokesIdx>& scvf, | ||
| 595 | const DiffusionCoefficientAveragingType diffCoeffAvgType = DiffusionCoefficientAveragingType::ffOnly) const | ||
| 596 | { | ||
| 597 | const auto& stokesContext = this->couplingManager().stokesCouplingContext(element, scvf); | ||
| 598 | const auto& stokesVolVars = stokesElemVolVars[scvf.insideScvIdx()]; | ||
| 599 | const auto& darcyVolVars = stokesContext.volVars; | ||
| 600 | |||
| 601 | const Scalar velocity = stokesElemFaceVars[scvf].velocitySelf(); | ||
| 602 | const bool insideIsUpstream = sign(velocity) == scvf.directionSign(); | ||
| 603 | |||
| 604 | return energyFlux_(stokesIdx, darcyIdx, fvGeometry, stokesContext.fvGeometry, scvf, | ||
| 605 | stokesVolVars, darcyVolVars, velocity * scvf.directionSign(), insideIsUpstream, diffCoeffAvgType); | ||
| 606 | } | ||
| 607 | |||
| 608 | private: | ||
| 609 | |||
| 610 | /*! | ||
| 611 | * \brief Evaluate the mole/mass flux across the interface. | ||
| 612 | */ | ||
| 613 | 111296 | Scalar massFlux_(const Scalar velocity, | |
| 614 | const Scalar insideDensity, | ||
| 615 | const Scalar outSideDensity, | ||
| 616 | bool insideIsUpstream) const | ||
| 617 | { | ||
| 618 |
4/4✓ Branch 0 taken 71757 times.
✓ Branch 1 taken 14039 times.
✓ Branch 2 taken 700 times.
✓ Branch 3 taken 24800 times.
|
111296 | return this->advectiveFlux(insideDensity, outSideDensity, velocity, insideIsUpstream); |
| 619 | } | ||
| 620 | |||
| 621 | /*! | ||
| 622 | * \brief Evaluate the energy flux across the interface. | ||
| 623 | */ | ||
| 624 | template<std::size_t i, std::size_t j, bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0> | ||
| 625 | Scalar energyFlux_(Dune::index_constant<i> domainI, | ||
| 626 | Dune::index_constant<j> domainJ, | ||
| 627 | const FVElementGeometry<i>& insideFvGeometry, | ||
| 628 | const FVElementGeometry<j>& outsideFvGeometry, | ||
| 629 | const SubControlVolumeFace<i>& scvf, | ||
| 630 | const VolumeVariables<i>& insideVolVars, | ||
| 631 | const VolumeVariables<j>& outsideVolVars, | ||
| 632 | const Scalar velocity, | ||
| 633 | const bool insideIsUpstream, | ||
| 634 | const DiffusionCoefficientAveragingType diffCoeffAvgType) const | ||
| 635 | { | ||
| 636 | Scalar flux(0.0); | ||
| 637 | |||
| 638 | const auto& insideScv = (*scvs(insideFvGeometry).begin()); | ||
| 639 | const auto& outsideScv = (*scvs(outsideFvGeometry).begin()); | ||
| 640 | |||
| 641 | // convective fluxes | ||
| 642 | const Scalar insideTerm = insideVolVars.density(couplingPhaseIdx(domainI)) * insideVolVars.enthalpy(couplingPhaseIdx(domainI)); | ||
| 643 | const Scalar outsideTerm = outsideVolVars.density(couplingPhaseIdx(domainJ)) * outsideVolVars.enthalpy(couplingPhaseIdx(domainJ)); | ||
| 644 | |||
| 645 | flux += this->advectiveFlux(insideTerm, outsideTerm, velocity, insideIsUpstream); | ||
| 646 | |||
| 647 | flux += this->conductiveEnergyFlux_(domainI, domainJ, insideFvGeometry, outsideFvGeometry, scvf, insideScv, outsideScv, insideVolVars, outsideVolVars, diffCoeffAvgType); | ||
| 648 | |||
| 649 | return flux; | ||
| 650 | } | ||
| 651 | |||
| 652 | }; | ||
| 653 | |||
| 654 | /*! | ||
| 655 | * \ingroup StokesDarcyCoupling | ||
| 656 | * \brief Coupling data specialization for compositional models. | ||
| 657 | */ | ||
| 658 | template<class MDTraits, class CouplingManager, bool enableEnergyBalance> | ||
| 659 | class StokesDarcyCouplingDataImplementation<MDTraits, CouplingManager, enableEnergyBalance, true> | ||
| 660 | : public StokesDarcyCouplingDataImplementationBase<MDTraits, CouplingManager> | ||
| 661 | { | ||
| 662 | using ParentType = StokesDarcyCouplingDataImplementationBase<MDTraits, CouplingManager>; | ||
| 663 | using Scalar = typename MDTraits::Scalar; | ||
| 664 | static constexpr auto stokesIdx = CouplingManager::stokesIdx; | ||
| 665 | static constexpr auto darcyIdx = CouplingManager::darcyIdx; | ||
| 666 | static constexpr auto stokesFaceIdx = CouplingManager::stokesFaceIdx; | ||
| 667 | static constexpr auto stokesCellCenterIdx = CouplingManager::stokesCellCenterIdx; | ||
| 668 | |||
| 669 | // the sub domain type tags | ||
| 670 | template<std::size_t id> | ||
| 671 | using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag; | ||
| 672 | |||
| 673 | template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>; | ||
| 674 | template<std::size_t id> using Element = typename GridGeometry<id>::GridView::template Codim<0>::Entity; | ||
| 675 | template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView; | ||
| 676 | template<std::size_t id> using SubControlVolumeFace = typename FVElementGeometry<id>::SubControlVolumeFace; | ||
| 677 | template<std::size_t id> using SubControlVolume = typename GridGeometry<id>::LocalView::SubControlVolume; | ||
| 678 | template<std::size_t id> using Indices = typename GetPropType<SubDomainTypeTag<id>, Properties::ModelTraits>::Indices; | ||
| 679 | template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView; | ||
| 680 | template<std::size_t id> using ElementFaceVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridFaceVariables>::LocalView; | ||
| 681 | template<std::size_t id> using VolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::VolumeVariables; | ||
| 682 | template<std::size_t id> using FluidSystem = GetPropType<SubDomainTypeTag<id>, Properties::FluidSystem>; | ||
| 683 | |||
| 684 | static constexpr auto numComponents = GetPropType<SubDomainTypeTag<stokesIdx>, Properties::ModelTraits>::numFluidComponents(); | ||
| 685 | static constexpr auto replaceCompEqIdx = GetPropType<SubDomainTypeTag<stokesIdx>, Properties::ModelTraits>::replaceCompEqIdx(); | ||
| 686 | static constexpr bool useMoles = GetPropType<SubDomainTypeTag<stokesIdx>, Properties::ModelTraits>::useMoles(); | ||
| 687 | static constexpr auto referenceSystemFormulation = GetPropType<SubDomainTypeTag<stokesIdx>, Properties::MolecularDiffusionType>::referenceSystemFormulation(); | ||
| 688 | |||
| 689 | static_assert(GetPropType<SubDomainTypeTag<darcyIdx>, Properties::ModelTraits>::numFluidComponents() == numComponents, "Both submodels must use the same number of components"); | ||
| 690 | static_assert(getPropValue<SubDomainTypeTag<darcyIdx>, Properties::UseMoles>() == useMoles, "Both submodels must either use moles or not"); | ||
| 691 | static_assert(getPropValue<SubDomainTypeTag<darcyIdx>, Properties::ReplaceCompEqIdx>() == replaceCompEqIdx, "Both submodels must use the same replaceCompEqIdx"); | ||
| 692 | static_assert(GetPropType<SubDomainTypeTag<darcyIdx>, Properties::MolecularDiffusionType>::referenceSystemFormulation() == referenceSystemFormulation, | ||
| 693 | "Both submodels must use the same reference system formulation for diffusion"); | ||
| 694 | |||
| 695 | using NumEqVector = Dune::FieldVector<Scalar, numComponents>; | ||
| 696 | |||
| 697 | using DiffusionCoefficientAveragingType = typename StokesDarcyCouplingOptions::DiffusionCoefficientAveragingType; | ||
| 698 | |||
| 699 | static constexpr bool isFicksLaw = IsFicksLaw<GetPropType<SubDomainTypeTag<stokesIdx>, Properties::MolecularDiffusionType>>(); | ||
| 700 | static_assert(isFicksLaw == IsFicksLaw<GetPropType<SubDomainTypeTag<darcyIdx>, Properties::MolecularDiffusionType>>(), | ||
| 701 | "Both submodels must use the same diffusion law."); | ||
| 702 | |||
| 703 | using ReducedComponentVector = Dune::FieldVector<Scalar, numComponents-1>; | ||
| 704 | using ReducedComponentMatrix = Dune::FieldMatrix<Scalar, numComponents-1, numComponents-1>; | ||
| 705 | |||
| 706 | using MolecularDiffusionType = GetPropType<SubDomainTypeTag<stokesIdx>, Properties::MolecularDiffusionType>; | ||
| 707 | public: | ||
| 708 | 8 | using ParentType::ParentType; | |
| 709 | using ParentType::couplingPhaseIdx; | ||
| 710 | using ParentType::couplingCompIdx; | ||
| 711 | |||
| 712 | /*! | ||
| 713 | * \brief Returns the mass flux across the coupling boundary as seen from the Darcy domain. | ||
| 714 | */ | ||
| 715 | 67424 | NumEqVector massCouplingCondition(const Element<darcyIdx>& element, | |
| 716 | const FVElementGeometry<darcyIdx>& fvGeometry, | ||
| 717 | const ElementVolumeVariables<darcyIdx>& darcyElemVolVars, | ||
| 718 | const SubControlVolumeFace<darcyIdx>& scvf, | ||
| 719 | const DiffusionCoefficientAveragingType diffCoeffAvgType = DiffusionCoefficientAveragingType::ffOnly) const | ||
| 720 | { | ||
| 721 | 67424 | const auto& darcyContext = this->couplingManager().darcyCouplingContext(element, scvf); | |
| 722 | 67424 | const auto& darcyVolVars = darcyElemVolVars[scvf.insideScvIdx()]; | |
| 723 | 67424 | const auto& stokesVolVars = darcyContext.volVars; | |
| 724 | 67424 | const auto& outsideScv = (*scvs(darcyContext.fvGeometry).begin()); | |
| 725 | |||
| 726 | 67424 | const Scalar velocity = darcyContext.velocity * scvf.unitOuterNormal(); | |
| 727 | 67424 | const bool insideIsUpstream = velocity > 0.0; | |
| 728 | |||
| 729 | 67424 | return massFlux_(darcyIdx, stokesIdx, fvGeometry, | |
| 730 | scvf, darcyVolVars, stokesVolVars, | ||
| 731 | outsideScv, velocity, insideIsUpstream, | ||
| 732 | 67424 | diffCoeffAvgType); | |
| 733 | } | ||
| 734 | |||
| 735 | /*! | ||
| 736 | * \brief Returns the mass flux across the coupling boundary as seen from the free-flow domain. | ||
| 737 | */ | ||
| 738 | 260974 | NumEqVector massCouplingCondition(const Element<stokesIdx>& element, | |
| 739 | const FVElementGeometry<stokesIdx>& fvGeometry, | ||
| 740 | const ElementVolumeVariables<stokesIdx>& stokesElemVolVars, | ||
| 741 | const ElementFaceVariables<stokesIdx>& stokesElemFaceVars, | ||
| 742 | const SubControlVolumeFace<stokesIdx>& scvf, | ||
| 743 | const DiffusionCoefficientAveragingType diffCoeffAvgType = DiffusionCoefficientAveragingType::ffOnly) const | ||
| 744 | { | ||
| 745 | 260974 | const auto& stokesContext = this->couplingManager().stokesCouplingContext(element, scvf); | |
| 746 | 260974 | const auto& stokesVolVars = stokesElemVolVars[scvf.insideScvIdx()]; | |
| 747 | 260974 | const auto& darcyVolVars = stokesContext.volVars; | |
| 748 | 260974 | const auto& outsideScv = (*scvs(stokesContext.fvGeometry).begin()); | |
| 749 | |||
| 750 | 260974 | const Scalar velocity = stokesElemFaceVars[scvf].velocitySelf(); | |
| 751 | 260974 | const bool insideIsUpstream = sign(velocity) == scvf.directionSign(); | |
| 752 | |||
| 753 | 260974 | return massFlux_(stokesIdx, darcyIdx, fvGeometry, | |
| 754 | scvf, stokesVolVars, darcyVolVars, | ||
| 755 | 260974 | outsideScv, velocity * scvf.directionSign(), | |
| 756 | 260974 | insideIsUpstream, diffCoeffAvgType); | |
| 757 | } | ||
| 758 | |||
| 759 | /*! | ||
| 760 | * \brief Returns the energy flux across the coupling boundary as seen from the Darcy domain. | ||
| 761 | */ | ||
| 762 | template<bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0> | ||
| 763 | 12320 | Scalar energyCouplingCondition(const Element<darcyIdx>& element, | |
| 764 | const FVElementGeometry<darcyIdx>& fvGeometry, | ||
| 765 | const ElementVolumeVariables<darcyIdx>& darcyElemVolVars, | ||
| 766 | const SubControlVolumeFace<darcyIdx>& scvf, | ||
| 767 | const DiffusionCoefficientAveragingType diffCoeffAvgType = DiffusionCoefficientAveragingType::ffOnly) const | ||
| 768 | { | ||
| 769 | 12320 | const auto& darcyContext = this->couplingManager().darcyCouplingContext(element, scvf); | |
| 770 | 12320 | const auto& darcyVolVars = darcyElemVolVars[scvf.insideScvIdx()]; | |
| 771 | 12320 | const auto& stokesVolVars = darcyContext.volVars; | |
| 772 | |||
| 773 | 12320 | const Scalar velocity = darcyContext.velocity * scvf.unitOuterNormal(); | |
| 774 | 12320 | const bool insideIsUpstream = velocity > 0.0; | |
| 775 | |||
| 776 | 12320 | return energyFlux_(darcyIdx, stokesIdx, fvGeometry, darcyContext.fvGeometry, scvf, | |
| 777 | 12320 | darcyVolVars, stokesVolVars, velocity, insideIsUpstream, diffCoeffAvgType); | |
| 778 | } | ||
| 779 | |||
| 780 | /*! | ||
| 781 | * \brief Returns the energy flux across the coupling boundary as seen from the free-flow domain. | ||
| 782 | */ | ||
| 783 | template<bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0> | ||
| 784 | 48356 | Scalar energyCouplingCondition(const Element<stokesIdx>& element, | |
| 785 | const FVElementGeometry<stokesIdx>& fvGeometry, | ||
| 786 | const ElementVolumeVariables<stokesIdx>& stokesElemVolVars, | ||
| 787 | const ElementFaceVariables<stokesIdx>& stokesElemFaceVars, | ||
| 788 | const SubControlVolumeFace<stokesIdx>& scvf, | ||
| 789 | const DiffusionCoefficientAveragingType diffCoeffAvgType = DiffusionCoefficientAveragingType::ffOnly) const | ||
| 790 | { | ||
| 791 | 48356 | const auto& stokesContext = this->couplingManager().stokesCouplingContext(element, scvf); | |
| 792 | 48356 | const auto& stokesVolVars = stokesElemVolVars[scvf.insideScvIdx()]; | |
| 793 | 48356 | const auto& darcyVolVars = stokesContext.volVars; | |
| 794 | |||
| 795 | 48356 | const Scalar velocity = stokesElemFaceVars[scvf].velocitySelf(); | |
| 796 | 48356 | const bool insideIsUpstream = sign(velocity) == scvf.directionSign(); | |
| 797 | |||
| 798 | 48356 | return energyFlux_(stokesIdx, darcyIdx, fvGeometry, stokesContext.fvGeometry, scvf, | |
| 799 | 48356 | stokesVolVars, darcyVolVars, velocity * scvf.directionSign(), insideIsUpstream, diffCoeffAvgType); | |
| 800 | } | ||
| 801 | |||
| 802 | protected: | ||
| 803 | |||
| 804 | /*! | ||
| 805 | * \brief Evaluate the compositional mole/mass flux across the interface. | ||
| 806 | */ | ||
| 807 | template<std::size_t i, std::size_t j> | ||
| 808 | 656796 | NumEqVector massFlux_(Dune::index_constant<i> domainI, | |
| 809 | Dune::index_constant<j> domainJ, | ||
| 810 | const FVElementGeometry<i>& insideFvGeometry, | ||
| 811 | const SubControlVolumeFace<i>& scvf, | ||
| 812 | const VolumeVariables<i>& insideVolVars, | ||
| 813 | const VolumeVariables<j>& outsideVolVars, | ||
| 814 | const SubControlVolume<j>& outsideScv, | ||
| 815 | const Scalar velocity, | ||
| 816 | const bool insideIsUpstream, | ||
| 817 | const DiffusionCoefficientAveragingType diffCoeffAvgType) const | ||
| 818 | { | ||
| 819 | 656796 | NumEqVector flux(0.0); | |
| 820 | 656796 | NumEqVector diffusiveFlux(0.0); | |
| 821 | |||
| 822 | 1912352 | auto moleOrMassFraction = [&](const auto& volVars, int phaseIdx, int compIdx) | |
| 823 | 1912352 | { return useMoles ? volVars.moleFraction(phaseIdx, compIdx) : volVars.massFraction(phaseIdx, compIdx); }; | |
| 824 | |||
| 825 |
2/2✓ Branch 0 taken 270887 times.
✓ Branch 1 taken 165937 times.
|
1912352 | auto moleOrMassDensity = [&](const auto& volVars, int phaseIdx) |
| 826 |
2/2✓ Branch 0 taken 227175 times.
✓ Branch 1 taken 336685 times.
|
1912352 | { return useMoles ? volVars.molarDensity(phaseIdx) : volVars.density(phaseIdx); }; |
| 827 | |||
| 828 | // treat the advective fluxes | ||
| 829 | 1973732 | auto insideTerm = [&](int compIdx) | |
| 830 |
2/2✓ Branch 0 taken 43712 times.
✓ Branch 1 taken 272540 times.
|
1316936 | { return moleOrMassFraction(insideVolVars, couplingPhaseIdx(domainI), compIdx) * moleOrMassDensity(insideVolVars, couplingPhaseIdx(domainI)); }; |
| 831 | |||
| 832 | 1682156 | auto outsideTerm = [&](int compIdx) | |
| 833 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 392856 times.
|
1025360 | { return moleOrMassFraction(outsideVolVars, couplingPhaseIdx(domainJ), compIdx) * moleOrMassDensity(outsideVolVars, couplingPhaseIdx(domainJ)); }; |
| 834 | |||
| 835 | |||
| 836 |
2/2✓ Branch 0 taken 709108 times.
✓ Branch 1 taken 328398 times.
|
2075012 | for (int compIdx = 0; compIdx < numComponents; ++compIdx) |
| 837 | { | ||
| 838 |
2/2✓ Branch 0 taken 85502 times.
✓ Branch 1 taken 85502 times.
|
1418216 | const int domainICompIdx = couplingCompIdx(domainI, compIdx); |
| 839 | 1418216 | const int domainJCompIdx = couplingCompIdx(domainJ, compIdx); | |
| 840 |
3/3✓ Branch 1 taken 442051 times.
✓ Branch 2 taken 223345 times.
✓ Branch 0 taken 43712 times.
|
2836432 | flux[domainICompIdx] += this->advectiveFlux(insideTerm(domainICompIdx), outsideTerm(domainJCompIdx), velocity, insideIsUpstream); |
| 841 | } | ||
| 842 | |||
| 843 | // treat the diffusive fluxes | ||
| 844 |
1/2✓ Branch 0 taken 67424 times.
✗ Branch 1 not taken.
|
656796 | const auto& insideScv = insideFvGeometry.scv(scvf.insideScvIdx()); |
| 845 | |||
| 846 | if (isFicksLaw) | ||
| 847 | 488972 | diffusiveFlux += diffusiveMolecularFluxFicksLaw_(domainI, domainJ, scvf, insideScv, outsideScv, insideVolVars, outsideVolVars, diffCoeffAvgType); | |
| 848 | else //maxwell stefan | ||
| 849 | 167824 | diffusiveFlux += diffusiveMolecularFluxMaxwellStefan_(domainI, domainJ, scvf, insideScv, outsideScv, insideVolVars, outsideVolVars); | |
| 850 | |||
| 851 | //convert to correct units if necessary | ||
| 852 | if (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged && useMoles) | ||
| 853 | { | ||
| 854 |
2/2✓ Branch 0 taken 709108 times.
✓ Branch 1 taken 328398 times.
|
2075012 | for (int compIdx = 0; compIdx < numComponents; ++compIdx) |
| 855 | { | ||
| 856 |
2/2✓ Branch 0 taken 85502 times.
✓ Branch 1 taken 85502 times.
|
1418216 | const int domainICompIdx = couplingCompIdx(domainI, compIdx); |
| 857 | 1418216 | diffusiveFlux[domainICompIdx] *= 1/FluidSystem<i>::molarMass(domainICompIdx); | |
| 858 | } | ||
| 859 | } | ||
| 860 | if (referenceSystemFormulation == ReferenceSystemFormulation::molarAveraged && !useMoles) | ||
| 861 | { | ||
| 862 | for (int compIdx = 0; compIdx < numComponents; ++compIdx) | ||
| 863 | { | ||
| 864 | const int domainICompIdx = couplingCompIdx(domainI, compIdx); | ||
| 865 | diffusiveFlux[domainICompIdx] *= FluidSystem<i>::molarMass(domainICompIdx); | ||
| 866 | } | ||
| 867 | } | ||
| 868 | |||
| 869 | 656796 | flux += diffusiveFlux; | |
| 870 | // convert to total mass/mole balance, if set be user | ||
| 871 | if (replaceCompEqIdx < numComponents) | ||
| 872 | flux[replaceCompEqIdx] = std::accumulate(flux.begin(), flux.end(), 0.0); | ||
| 873 | |||
| 874 | 656796 | return flux; | |
| 875 | } | ||
| 876 | |||
| 877 | 60676 | Scalar getComponentEnthalpy(const VolumeVariables<stokesIdx>& volVars, int phaseIdx, int compIdx) const | |
| 878 | { | ||
| 879 | 60676 | return FluidSystem<stokesIdx>::componentEnthalpy(volVars.fluidState(), 0, compIdx); | |
| 880 | } | ||
| 881 | |||
| 882 | 60676 | Scalar getComponentEnthalpy(const VolumeVariables<darcyIdx>& volVars, int phaseIdx, int compIdx) const | |
| 883 | { | ||
| 884 | 60676 | return FluidSystem<darcyIdx>::componentEnthalpy(volVars.fluidState(), phaseIdx, compIdx); | |
| 885 | } | ||
| 886 | |||
| 887 | /*! | ||
| 888 | * \brief Evaluate the diffusive mole/mass flux across the interface. | ||
| 889 | */ | ||
| 890 | template<std::size_t i, std::size_t j> | ||
| 891 | 83912 | NumEqVector diffusiveMolecularFluxMaxwellStefan_(Dune::index_constant<i> domainI, | |
| 892 | Dune::index_constant<j> domainJ, | ||
| 893 | const SubControlVolumeFace<i>& scvfI, | ||
| 894 | const SubControlVolume<i>& scvI, | ||
| 895 | const SubControlVolume<j>& scvJ, | ||
| 896 | const VolumeVariables<i>& volVarsI, | ||
| 897 | const VolumeVariables<j>& volVarsJ) const | ||
| 898 | { | ||
| 899 | 83912 | NumEqVector diffusiveFlux(0.0); | |
| 900 | |||
| 901 | 83912 | const Scalar insideDistance = this->getDistance_(scvI, scvfI); | |
| 902 | 83912 | const Scalar outsideDistance = this->getDistance_(scvJ, scvfI); | |
| 903 | |||
| 904 | 83912 | ReducedComponentVector moleFracInside(0.0); | |
| 905 | 83912 | ReducedComponentVector moleFracOutside(0.0); | |
| 906 | 83912 | ReducedComponentVector reducedFlux(0.0); | |
| 907 | 83912 | ReducedComponentMatrix reducedDiffusionMatrixInside(0.0); | |
| 908 | 83912 | ReducedComponentMatrix reducedDiffusionMatrixOutside(0.0); | |
| 909 | |||
| 910 | //calculate the mole fraction vectors | ||
| 911 | 220136 | for (int compIdx = 0; compIdx < numComponents-1; compIdx++) | |
| 912 | { | ||
| 913 | 136224 | const int domainICompIdx = couplingCompIdx(domainI, compIdx); | |
| 914 | 136224 | const int domainJCompIdx = couplingCompIdx(domainJ, compIdx); | |
| 915 | |||
| 916 | 272448 | assert(FluidSystem<i>::componentName(domainICompIdx) == FluidSystem<j>::componentName(domainJCompIdx)); | |
| 917 | |||
| 918 | //calculate x_inside | ||
| 919 | 136224 | const Scalar xInside = volVarsI.moleFraction(couplingPhaseIdx(domainI), domainICompIdx); | |
| 920 | //calculate outside molefraction with the respective transmissibility | ||
| 921 | 136224 | const Scalar xOutside = volVarsJ.moleFraction(couplingPhaseIdx(domainJ), domainJCompIdx); | |
| 922 | 136224 | moleFracInside[domainICompIdx] = xInside; | |
| 923 | 136224 | moleFracOutside[domainICompIdx] = xOutside; | |
| 924 | } | ||
| 925 | |||
| 926 | //now we have to do the tpfa: J_i = -J_j which leads to: J_i = -rho_i Bi^-1 omegai(x*-xi) with x* = (omegai rho_i Bi^-1 + omegaj rho_j Bj^-1)^-1 (xi omegai rho_i Bi^-1 + xj omegaj rho_j Bj^-1) with i inside and j outside. | ||
| 927 | |||
| 928 | //first set up the matrices containing the binary diffusion coefficients and mole fractions | ||
| 929 | |||
| 930 | //inside matrix. KIdx and LIdx are the indices for the k and l-th component, N for the n-th component | ||
| 931 | 220136 | for (int compKIdx = 0; compKIdx < numComponents-1; compKIdx++) | |
| 932 | { | ||
| 933 | 136224 | const int domainICompKIdx = couplingCompIdx(domainI, compKIdx); | |
| 934 | 136224 | const Scalar xk = volVarsI.moleFraction(couplingPhaseIdx(domainI), domainICompKIdx); | |
| 935 | 136224 | const Scalar avgMolarMass = volVarsI.averageMolarMass(couplingPhaseIdx(domainI)); | |
| 936 | 136224 | const Scalar Mn = FluidSystem<i>::molarMass(numComponents-1); | |
| 937 | 136224 | const Scalar tkn = volVarsI.effectiveDiffusionCoefficient(couplingPhaseIdx(domainI), domainICompKIdx, couplingCompIdx(domainI, numComponents-1)); | |
| 938 | |||
| 939 | // set the entries of the diffusion matrix of the diagonal | ||
| 940 | 136224 | reducedDiffusionMatrixInside[domainICompKIdx][domainICompKIdx] += xk*avgMolarMass/(tkn*Mn); | |
| 941 | |||
| 942 | 513296 | for (int compLIdx = 0; compLIdx < numComponents; compLIdx++) | |
| 943 | { | ||
| 944 | 377072 | const int domainICompLIdx = couplingCompIdx(domainI, compLIdx); | |
| 945 | |||
| 946 | // we don't want to calculate e.g. water in water diffusion | ||
| 947 | 377072 | if (domainICompKIdx == domainICompLIdx) | |
| 948 | 136224 | continue; | |
| 949 | |||
| 950 | 240848 | const Scalar xl = volVarsI.moleFraction(couplingPhaseIdx(domainI), domainICompLIdx); | |
| 951 | 240848 | const Scalar Mk = FluidSystem<i>::molarMass(domainICompKIdx); | |
| 952 | 240848 | const Scalar Ml = FluidSystem<i>::molarMass(domainICompLIdx); | |
| 953 | 240848 | const Scalar tkl = volVarsI.effectiveDiffusionCoefficient(couplingPhaseIdx(domainI), domainICompKIdx, domainICompLIdx); | |
| 954 | 240848 | reducedDiffusionMatrixInside[domainICompKIdx][domainICompKIdx] += xl*avgMolarMass/(tkl*Mk); | |
| 955 | |||
| 956 | 209248 | if (domainICompLIdx < numComponents-1) | |
| 957 | 104624 | reducedDiffusionMatrixInside[domainICompKIdx][domainICompLIdx] += xk*(avgMolarMass/(tkn*Mn) - avgMolarMass/(tkl*Ml)); | |
| 958 | } | ||
| 959 | } | ||
| 960 | //outside matrix | ||
| 961 | 220136 | for (int compKIdx = 0; compKIdx < numComponents-1; compKIdx++) | |
| 962 | { | ||
| 963 | 136224 | const int domainJCompKIdx = couplingCompIdx(domainJ, compKIdx); | |
| 964 | 136224 | const int domainICompKIdx = couplingCompIdx(domainI, compKIdx); | |
| 965 | |||
| 966 | 136224 | const Scalar xk = volVarsJ.moleFraction(couplingPhaseIdx(domainJ), domainJCompKIdx); | |
| 967 | 136224 | const Scalar avgMolarMass = volVarsJ.averageMolarMass(couplingPhaseIdx(domainJ)); | |
| 968 | 136224 | const Scalar Mn = FluidSystem<j>::molarMass(numComponents-1); | |
| 969 | 136224 | const Scalar tkn = volVarsJ.effectiveDiffusionCoefficient(couplingPhaseIdx(domainJ), domainJCompKIdx, couplingCompIdx(domainJ, numComponents-1)); | |
| 970 | |||
| 971 | // set the entries of the diffusion matrix of the diagonal | ||
| 972 | 136224 | reducedDiffusionMatrixOutside[domainICompKIdx][domainICompKIdx] += xk*avgMolarMass/(tkn*Mn); | |
| 973 | |||
| 974 | 513296 | for (int compLIdx = 0; compLIdx < numComponents; compLIdx++) | |
| 975 | { | ||
| 976 | 377072 | const int domainJCompLIdx = couplingCompIdx(domainJ, compLIdx); | |
| 977 | 377072 | const int domainICompLIdx = couplingCompIdx(domainI, compLIdx); | |
| 978 | |||
| 979 | // we don't want to calculate e.g. water in water diffusion | ||
| 980 | 377072 | if (domainJCompLIdx == domainJCompKIdx) | |
| 981 | 136224 | continue; | |
| 982 | |||
| 983 | 240848 | const Scalar xl = volVarsJ.moleFraction(couplingPhaseIdx(domainJ), domainJCompLIdx); | |
| 984 | 240848 | const Scalar Mk = FluidSystem<j>::molarMass(domainJCompKIdx); | |
| 985 | 240848 | const Scalar Ml = FluidSystem<j>::molarMass(domainJCompLIdx); | |
| 986 | 240848 | const Scalar tkl = volVarsJ.effectiveDiffusionCoefficient(couplingPhaseIdx(domainJ), domainJCompKIdx, domainJCompLIdx); | |
| 987 | 240848 | reducedDiffusionMatrixOutside[domainICompKIdx][domainICompKIdx] += xl*avgMolarMass/(tkl*Mk); | |
| 988 | |||
| 989 | 209248 | if (domainICompLIdx < numComponents-1) | |
| 990 | 104624 | reducedDiffusionMatrixOutside[domainICompKIdx][domainICompLIdx] += xk*(avgMolarMass/(tkn*Mn) - avgMolarMass/(tkl*Ml)); | |
| 991 | } | ||
| 992 | } | ||
| 993 | |||
| 994 | 83912 | const Scalar omegai = 1/insideDistance; | |
| 995 | 83912 | const Scalar omegaj = 1/outsideDistance; | |
| 996 | |||
| 997 | 83912 | reducedDiffusionMatrixInside.invert(); | |
| 998 | 136224 | reducedDiffusionMatrixInside *= omegai*volVarsI.density(couplingPhaseIdx(domainI)); | |
| 999 | 83912 | reducedDiffusionMatrixOutside.invert(); | |
| 1000 | 83912 | reducedDiffusionMatrixOutside *= omegaj*volVarsJ.density(couplingPhaseIdx(domainJ)); | |
| 1001 | |||
| 1002 | //in the helpervector we store the values for x* | ||
| 1003 | 83912 | ReducedComponentVector helperVector(0.0); | |
| 1004 | 83912 | ReducedComponentVector gradientVectori(0.0); | |
| 1005 | 83912 | ReducedComponentVector gradientVectorj(0.0); | |
| 1006 | |||
| 1007 | 220136 | reducedDiffusionMatrixInside.mv(moleFracInside, gradientVectori); | |
| 1008 | 83912 | reducedDiffusionMatrixOutside.mv(moleFracOutside, gradientVectorj); | |
| 1009 | |||
| 1010 | 83912 | auto gradientVectorij = (gradientVectori + gradientVectorj); | |
| 1011 | |||
| 1012 | //add the two matrixes to each other | ||
| 1013 | 83912 | reducedDiffusionMatrixOutside += reducedDiffusionMatrixInside; | |
| 1014 | |||
| 1015 | 83912 | reducedDiffusionMatrixOutside.solve(helperVector, gradientVectorij); | |
| 1016 | |||
| 1017 | //Bi^-1 omegai rho_i (x*-xi). As we previously multiplied rho_i and omega_i with the insidematrix, this does not need to be done again | ||
| 1018 | 220136 | helperVector -=moleFracInside; | |
| 1019 | 220136 | reducedDiffusionMatrixInside.mv(helperVector, reducedFlux); | |
| 1020 | |||
| 1021 | 220136 | reducedFlux *= -1; | |
| 1022 | |||
| 1023 | 220136 | for (int compIdx = 0; compIdx < numComponents-1; compIdx++) | |
| 1024 | { | ||
| 1025 | 136224 | const int domainICompIdx = couplingCompIdx(domainI, compIdx); | |
| 1026 | 136224 | diffusiveFlux[domainICompIdx] = reducedFlux[domainICompIdx]; | |
| 1027 | 136224 | diffusiveFlux[couplingCompIdx(domainI, numComponents-1)] -= reducedFlux[domainICompIdx]; | |
| 1028 | } | ||
| 1029 | 83912 | return diffusiveFlux; | |
| 1030 | } | ||
| 1031 | |||
| 1032 | template<std::size_t i, std::size_t j> | ||
| 1033 | 426514 | NumEqVector diffusiveMolecularFluxFicksLaw_(Dune::index_constant<i> domainI, | |
| 1034 | Dune::index_constant<j> domainJ, | ||
| 1035 | const SubControlVolumeFace<i>& scvfI, | ||
| 1036 | const SubControlVolume<i>& scvI, | ||
| 1037 | const SubControlVolume<j>& scvJ, | ||
| 1038 | const VolumeVariables<i>& volVarsI, | ||
| 1039 | const VolumeVariables<j>& volVarsJ, | ||
| 1040 | const DiffusionCoefficientAveragingType diffCoeffAvgType) const | ||
| 1041 | { | ||
| 1042 | 426514 | NumEqVector diffusiveFlux(0.0); | |
| 1043 | |||
| 1044 | 426514 | const Scalar rhoInside = massOrMolarDensity(volVarsI, referenceSystemFormulation, couplingPhaseIdx(domainI)); | |
| 1045 | 426514 | const Scalar rhoOutside = massOrMolarDensity(volVarsJ, referenceSystemFormulation, couplingPhaseIdx(domainJ)); | |
| 1046 | 426514 | const Scalar avgDensity = 0.5 * rhoInside + 0.5 * rhoOutside; | |
| 1047 | |||
| 1048 | 426514 | const Scalar insideDistance = this->getDistance_(scvI, scvfI); | |
| 1049 | 426514 | const Scalar outsideDistance = this->getDistance_(scvJ, scvfI); | |
| 1050 | |||
| 1051 |
2/2✓ Branch 1 taken 121352 times.
✓ Branch 2 taken 121352 times.
|
853028 | for (int compIdx = 1; compIdx < numComponents; ++compIdx) |
| 1052 | { | ||
| 1053 | 426514 | const int domainIMainCompIdx = couplingPhaseIdx(domainI); | |
| 1054 | 426514 | const int domainJMainCompIdx = couplingPhaseIdx(domainJ); | |
| 1055 | 426514 | const int domainICompIdx = couplingCompIdx(domainI, compIdx); | |
| 1056 | 367570 | const int domainJCompIdx = couplingCompIdx(domainJ, compIdx); | |
| 1057 | |||
| 1058 |
3/4✓ Branch 1 taken 24640 times.
✓ Branch 2 taken 96712 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 24640 times.
|
1083598 | assert(FluidSystem<i>::componentName(domainICompIdx) == FluidSystem<j>::componentName(domainJCompIdx)); |
| 1059 | |||
| 1060 | 426514 | const Scalar massOrMoleFractionInside = massOrMoleFraction(volVarsI, referenceSystemFormulation, couplingPhaseIdx(domainI), domainICompIdx); | |
| 1061 | 426514 | const Scalar massOrMoleFractionOutside = massOrMoleFraction(volVarsJ, referenceSystemFormulation, couplingPhaseIdx(domainJ), domainJCompIdx); | |
| 1062 | |||
| 1063 | 426514 | const Scalar deltaMassOrMoleFrac = massOrMoleFractionOutside - massOrMoleFractionInside; | |
| 1064 | 426514 | const Scalar tij = this->transmissibility_(domainI, | |
| 1065 | domainJ, | ||
| 1066 | insideDistance, | ||
| 1067 | outsideDistance, | ||
| 1068 | volVarsI.effectiveDiffusionCoefficient(couplingPhaseIdx(domainI), domainIMainCompIdx, domainICompIdx), | ||
| 1069 | volVarsJ.effectiveDiffusionCoefficient(couplingPhaseIdx(domainJ), domainJMainCompIdx, domainJCompIdx), | ||
| 1070 | diffCoeffAvgType); | ||
| 1071 | 426514 | diffusiveFlux[domainICompIdx] += -avgDensity * tij * deltaMassOrMoleFrac; | |
| 1072 | } | ||
| 1073 | |||
| 1074 | 426514 | const Scalar cumulativeFlux = std::accumulate(diffusiveFlux.begin(), diffusiveFlux.end(), 0.0); | |
| 1075 | 426514 | diffusiveFlux[couplingCompIdx(domainI, 0)] = -cumulativeFlux; | |
| 1076 | |||
| 1077 | 426514 | return diffusiveFlux; | |
| 1078 | } | ||
| 1079 | |||
| 1080 | /*! | ||
| 1081 | * \brief Evaluate the energy flux across the interface. | ||
| 1082 | */ | ||
| 1083 | template<std::size_t i, std::size_t j, bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0> | ||
| 1084 | 60676 | Scalar energyFlux_(Dune::index_constant<i> domainI, | |
| 1085 | Dune::index_constant<j> domainJ, | ||
| 1086 | const FVElementGeometry<i>& insideFvGeometry, | ||
| 1087 | const FVElementGeometry<j>& outsideFvGeometry, | ||
| 1088 | const SubControlVolumeFace<i>& scvf, | ||
| 1089 | const VolumeVariables<i>& insideVolVars, | ||
| 1090 | const VolumeVariables<j>& outsideVolVars, | ||
| 1091 | const Scalar velocity, | ||
| 1092 | const bool insideIsUpstream, | ||
| 1093 | const DiffusionCoefficientAveragingType diffCoeffAvgType) const | ||
| 1094 | { | ||
| 1095 | 60676 | Scalar flux(0.0); | |
| 1096 | |||
| 1097 | 60676 | const auto& insideScv = (*scvs(insideFvGeometry).begin()); | |
| 1098 | 60676 | const auto& outsideScv = (*scvs(outsideFvGeometry).begin()); | |
| 1099 | |||
| 1100 | // convective fluxes | ||
| 1101 | 60676 | const Scalar insideTerm = insideVolVars.density(couplingPhaseIdx(domainI)) * insideVolVars.enthalpy(couplingPhaseIdx(domainI)); | |
| 1102 | 60676 | const Scalar outsideTerm = outsideVolVars.density(couplingPhaseIdx(domainJ)) * outsideVolVars.enthalpy(couplingPhaseIdx(domainJ)); | |
| 1103 | |||
| 1104 | 60676 | flux += this->advectiveFlux(insideTerm, outsideTerm, velocity, insideIsUpstream); | |
| 1105 | |||
| 1106 | 60676 | flux += this->conductiveEnergyFlux_(domainI, domainJ, insideFvGeometry, outsideFvGeometry, scvf, insideScv, outsideScv, insideVolVars, outsideVolVars, diffCoeffAvgType); | |
| 1107 | |||
| 1108 | 60676 | auto diffusiveFlux = isFicksLaw ? diffusiveMolecularFluxFicksLaw_(domainI, domainJ, scvf, insideScv, outsideScv, insideVolVars, outsideVolVars, diffCoeffAvgType) | |
| 1109 | : diffusiveMolecularFluxMaxwellStefan_(domainI, domainJ, scvf, insideScv, outsideScv, insideVolVars, outsideVolVars); | ||
| 1110 | |||
| 1111 | |||
| 1112 | 182028 | for (int compIdx = 0; compIdx < diffusiveFlux.size(); ++compIdx) | |
| 1113 | { | ||
| 1114 | 121352 | const int domainICompIdx = couplingCompIdx(domainI, compIdx); | |
| 1115 | 121352 | const int domainJCompIdx = couplingCompIdx(domainJ, compIdx); | |
| 1116 | |||
| 1117 | 121352 | const bool insideDiffFluxIsUpstream = diffusiveFlux[domainICompIdx] > 0; | |
| 1118 | 121352 | const Scalar componentEnthalpy = insideDiffFluxIsUpstream ? | |
| 1119 | 121352 | getComponentEnthalpy(insideVolVars, couplingPhaseIdx(domainI), domainICompIdx) | |
| 1120 | 121352 | : getComponentEnthalpy(outsideVolVars, couplingPhaseIdx(domainJ), domainJCompIdx); | |
| 1121 | |||
| 1122 | if (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged) | ||
| 1123 | 121352 | flux += diffusiveFlux[domainICompIdx] * componentEnthalpy; | |
| 1124 | else | ||
| 1125 | flux += diffusiveFlux[domainICompIdx] * FluidSystem<i>::molarMass(domainICompIdx) * componentEnthalpy; | ||
| 1126 | } | ||
| 1127 | |||
| 1128 | 60676 | return flux; | |
| 1129 | } | ||
| 1130 | }; | ||
| 1131 | |||
| 1132 | } // end namespace Dumux | ||
| 1133 | |||
| 1134 | #endif | ||
| 1135 |