| 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 FreeFlowPoreNetworkCoupling | ||
| 10 | * \copydoc Dumux::FreeFlowPoreNetworkCouplingConditions | ||
| 11 | */ | ||
| 12 | |||
| 13 | #ifndef DUMUX_MD_FREEFLOW_PORENETWORK_COUPLINGCONDITIONS_HH | ||
| 14 | #define DUMUX_MD_FREEFLOW_PORENETWORK_COUPLINGCONDITIONS_HH | ||
| 15 | |||
| 16 | #include <numeric> | ||
| 17 | |||
| 18 | #include <dune/common/exceptions.hh> | ||
| 19 | #include <dune/common/fvector.hh> | ||
| 20 | |||
| 21 | #include <dumux/common/properties.hh> | ||
| 22 | #include <dumux/common/math.hh> | ||
| 23 | #include <dumux/discretization/cellcentered/tpfa/computetransmissibility.hh> | ||
| 24 | #include <dumux/flux/fickslaw_fwd.hh> | ||
| 25 | #include <dumux/freeflow/navierstokes/momentum/velocitygradients.hh> | ||
| 26 | #include <dumux/multidomain/boundary/freeflowporousmedium/traits.hh> | ||
| 27 | #include <dumux/multidomain/boundary/freeflowporousmedium/indexhelper.hh> | ||
| 28 | |||
| 29 | namespace Dumux { | ||
| 30 | |||
| 31 | template<class MDTraits, class CouplingManager, bool enableEnergyBalance, bool isCompositional> | ||
| 32 | class FreeFlowPoreNetworkCouplingConditionsImplementation; | ||
| 33 | |||
| 34 | /*! | ||
| 35 | * \ingroup FreeFlowPoreNetworkCoupling | ||
| 36 | * \brief Coupling conditions for the coupling of a pore-network model | ||
| 37 | * with a (Navier-)Stokes model (staggerd grid). | ||
| 38 | */ | ||
| 39 | template<class MDTraits, class CouplingManager> | ||
| 40 | using FreeFlowPoreNetworkCouplingConditions | ||
| 41 | = FreeFlowPoreNetworkCouplingConditionsImplementation< | ||
| 42 | MDTraits, CouplingManager, | ||
| 43 | GetPropType<typename MDTraits::template SubDomain<1>::TypeTag, Properties::ModelTraits>::enableEnergyBalance(), | ||
| 44 | (GetPropType<typename MDTraits::template SubDomain<1>::TypeTag, Properties::ModelTraits>::numFluidComponents() > 1) | ||
| 45 | >; | ||
| 46 | |||
| 47 | /*! | ||
| 48 | * \ingroup FreeFlowPoreNetworkCoupling | ||
| 49 | * \brief A base class which provides some common methods used for free-flow/pore-network coupling. | ||
| 50 | */ | ||
| 51 | template<class MDTraits, class CouplingManager> | ||
| 52 | class FreeFlowPoreNetworkCouplingConditionsImplementationBase | ||
| 53 | { | ||
| 54 | using Scalar = typename MDTraits::Scalar; | ||
| 55 | |||
| 56 | template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag; | ||
| 57 | template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>; | ||
| 58 | template<std::size_t id> using Element = typename GridGeometry<id>::GridView::template Codim<0>::Entity; | ||
| 59 | template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView; | ||
| 60 | template<std::size_t id> using SubControlVolumeFace = typename GridGeometry<id>::LocalView::SubControlVolumeFace; | ||
| 61 | template<std::size_t id> using SubControlVolume = typename GridGeometry<id>::LocalView::SubControlVolume; | ||
| 62 | template<std::size_t id> using Indices = typename GetPropType<SubDomainTypeTag<id>, Properties::ModelTraits>::Indices; | ||
| 63 | template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView; | ||
| 64 | template<std::size_t id> using VolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::VolumeVariables; | ||
| 65 | template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>; | ||
| 66 | template<std::size_t id> using FluidSystem = GetPropType<SubDomainTypeTag<id>, Properties::FluidSystem>; | ||
| 67 | template<std::size_t id> using ModelTraits = GetPropType<SubDomainTypeTag<id>, Properties::ModelTraits>; | ||
| 68 | template<std::size_t id> using GlobalPosition = typename Element<id>::Geometry::GlobalCoordinate; | ||
| 69 | template<std::size_t id> using NumEqVector = typename Problem<id>::Traits::NumEqVector; | ||
| 70 | |||
| 71 | using VelocityGradients = StaggeredVelocityGradients; | ||
| 72 | |||
| 73 | public: | ||
| 74 | static constexpr auto freeFlowMomentumIndex = CouplingManager::freeFlowMomentumIndex; | ||
| 75 | static constexpr auto freeFlowMassIndex = CouplingManager::freeFlowMassIndex; | ||
| 76 | static constexpr auto poreNetworkIndex = CouplingManager::poreNetworkIndex; | ||
| 77 | private: | ||
| 78 | |||
| 79 | using AdvectionType = GetPropType<SubDomainTypeTag<poreNetworkIndex>, Properties::AdvectionType>; | ||
| 80 | |||
| 81 | static constexpr bool adapterUsed = ModelTraits<poreNetworkIndex>::numFluidPhases() > 1; | ||
| 82 | using IndexHelper = FreeFlowPorousMediumCoupling::IndexHelper<freeFlowMassIndex, poreNetworkIndex, FluidSystem<freeFlowMassIndex>, adapterUsed>; | ||
| 83 | |||
| 84 | static constexpr int enableEnergyBalance = GetPropType<SubDomainTypeTag<freeFlowMassIndex>, Properties::ModelTraits>::enableEnergyBalance(); | ||
| 85 | static_assert(GetPropType<SubDomainTypeTag<poreNetworkIndex>, Properties::ModelTraits>::enableEnergyBalance() == enableEnergyBalance, | ||
| 86 | "All submodels must both be either isothermal or non-isothermal"); | ||
| 87 | |||
| 88 | static_assert(FreeFlowPorousMediumCoupling::IsSameFluidSystem<FluidSystem<freeFlowMassIndex>, | ||
| 89 | FluidSystem<poreNetworkIndex>>::value, | ||
| 90 | "All submodels must use the same fluid system"); | ||
| 91 | |||
| 92 | using VelocityVector = GlobalPosition<freeFlowMomentumIndex>; | ||
| 93 | |||
| 94 | public: | ||
| 95 | /*! | ||
| 96 | * \brief Returns the corresponding phase index needed for coupling. | ||
| 97 | */ | ||
| 98 | template<std::size_t i> | ||
| 99 | static constexpr auto couplingPhaseIdx(Dune::index_constant<i> id, int coupledPhaseIdx = 0) | ||
| 100 | { return IndexHelper::couplingPhaseIdx(id, coupledPhaseIdx); } | ||
| 101 | |||
| 102 | /*! | ||
| 103 | * \brief Returns the corresponding component index needed for coupling. | ||
| 104 | */ | ||
| 105 | template<std::size_t i> | ||
| 106 | static constexpr auto couplingCompIdx(Dune::index_constant<i> id, int coupledCompIdx) | ||
| 107 | { return IndexHelper::couplingCompIdx(id, coupledCompIdx); } | ||
| 108 | |||
| 109 | /*! | ||
| 110 | * \brief Returns the momentum flux across the coupling boundary. | ||
| 111 | * | ||
| 112 | * For the normal momentum coupling, the porous medium side of the coupling condition | ||
| 113 | * is evaluated, i.e. -[p n]^pm. | ||
| 114 | * | ||
| 115 | */ | ||
| 116 | template<class Context> | ||
| 117 | 23798 | static NumEqVector<freeFlowMomentumIndex> momentumCouplingCondition(const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry, | |
| 118 | const SubControlVolumeFace<freeFlowMomentumIndex>& scvf, | ||
| 119 | const ElementVolumeVariables<freeFlowMomentumIndex>& elemVolVars, | ||
| 120 | const Context& context) | ||
| 121 | { | ||
| 122 | 23798 | NumEqVector<freeFlowMomentumIndex> momentumFlux(0.0); | |
| 123 | 23798 | const auto pnmPhaseIdx = couplingPhaseIdx(poreNetworkIndex); | |
| 124 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 23798 times.
|
23798 | const auto [pnmPressure, pnmViscosity] = [&] |
| 125 | { | ||
| 126 |
3/4✓ Branch 0 taken 23798 times.
✓ Branch 1 taken 23798 times.
✓ Branch 2 taken 47596 times.
✗ Branch 3 not taken.
|
47596 | for (const auto& scv : scvs(context.fvGeometry)) |
| 127 | { | ||
| 128 |
2/2✓ Branch 0 taken 23798 times.
✓ Branch 1 taken 23798 times.
|
47596 | if (scv.dofIndex() == context.poreNetworkDofIdx) |
| 129 | 23798 | return std::make_pair(context.elemVolVars[scv].pressure(pnmPhaseIdx), context.elemVolVars[scv].viscosity(pnmPhaseIdx)); | |
| 130 | } | ||
| 131 | ✗ | DUNE_THROW(Dune::InvalidStateException, "No coupled scv found"); | |
| 132 | 23798 | }(); | |
| 133 | |||
| 134 | // set p_freeflow = p_PNM | ||
| 135 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 23798 times.
|
23798 | momentumFlux[0] = pnmPressure; |
| 136 | |||
| 137 | // normalize pressure | ||
| 138 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 23798 times.
|
23798 | momentumFlux[0] -= elemVolVars.gridVolVars().problem().referencePressure(fvGeometry.element(), fvGeometry, scvf); |
| 139 | |||
| 140 | // Explicitly account for dv_i/dx_i, which is NOT part of the actual coupling condition. We do it here for convenience so | ||
| 141 | // we do not forget to set it in the problem. We assume that the velocity gradient at the boundary towards the interface is the same | ||
| 142 | // as the one in the center of the element. TODO check sign | ||
| 143 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 23798 times.
|
23798 | const auto& scv = fvGeometry.scv(scvf.insideScvIdx()); |
| 144 | 23798 | const auto& frontalInternalScvf = (*scvfs(fvGeometry, scv).begin()); | |
| 145 | 23798 | momentumFlux[0] -= 2*VelocityGradients::velocityGradII(fvGeometry, frontalInternalScvf, elemVolVars) * pnmViscosity; | |
| 146 | |||
| 147 | // We do NOT consider the inertia term here. If included, Newton convergence decreases drastically and the solution even does not converge to a reference solution. | ||
| 148 | // We furthermore assume creeping flow within the boundary layer thus neglecting this term is physically justified. | ||
| 149 | |||
| 150 | 23798 | momentumFlux[0] *= scvf.directionSign(); | |
| 151 | |||
| 152 | 23798 | return momentumFlux; | |
| 153 | } | ||
| 154 | |||
| 155 | template<class Context> | ||
| 156 | 93628 | static VelocityVector interfaceThroatVelocity(const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry, | |
| 157 | const SubControlVolumeFace<freeFlowMomentumIndex>& scvf, | ||
| 158 | const Context& context) | ||
| 159 | { | ||
| 160 |
1/2✓ Branch 0 taken 93628 times.
✗ Branch 1 not taken.
|
93628 | const auto& pnmElement = context.fvGeometry.element(); |
| 161 |
1/2✓ Branch 0 taken 93628 times.
✗ Branch 1 not taken.
|
93628 | const auto& pnmFVGeometry = context.fvGeometry; |
| 162 | 93628 | const auto& pnmScvf = pnmFVGeometry.scvf(0); | |
| 163 | 93628 | const auto& pnmElemVolVars = context.elemVolVars; | |
| 164 |
1/2✓ Branch 0 taken 93628 times.
✗ Branch 1 not taken.
|
93628 | const auto& pnmElemFluxVarsCache = context.elemFluxVarsCache; |
| 165 |
1/2✓ Branch 0 taken 93628 times.
✗ Branch 1 not taken.
|
93628 | const auto& pnmProblem = pnmElemVolVars.gridVolVars().problem(); |
| 166 | |||
| 167 |
1/2✓ Branch 0 taken 93628 times.
✗ Branch 1 not taken.
|
93628 | const auto pnmPhaseIdx = couplingPhaseIdx(poreNetworkIndex); |
| 168 |
1/2✓ Branch 0 taken 93628 times.
✗ Branch 1 not taken.
|
93628 | const Scalar area = pnmElemFluxVarsCache[pnmScvf].throatCrossSectionalArea(pnmPhaseIdx); |
| 169 | |||
| 170 | // only proceed if area > 0 in order to prevent division by zero (e.g., when the throat was not invaded yet) | ||
| 171 |
1/2✓ Branch 0 taken 93628 times.
✗ Branch 1 not taken.
|
93628 | if (area > 0.0) |
| 172 | { | ||
| 173 | using PNMFluxVariables = GetPropType<SubDomainTypeTag<poreNetworkIndex>, Properties::FluxVariables>; | ||
| 174 | 93628 | PNMFluxVariables fluxVars; | |
| 175 | 93628 | fluxVars.init(pnmProblem, pnmElement, pnmFVGeometry, pnmElemVolVars, pnmScvf, pnmElemFluxVarsCache); | |
| 176 | |||
| 177 | 187256 | const Scalar flux = fluxVars.advectiveFlux(pnmPhaseIdx, [pnmPhaseIdx](const auto& volVars){ return volVars.mobility(pnmPhaseIdx);}); | |
| 178 | |||
| 179 | // account for the orientation of the bulk face. | ||
| 180 | 187256 | VelocityVector velocity = (pnmElement.geometry().corner(1) - pnmElement.geometry().corner(0)); | |
| 181 | 93628 | velocity /= velocity.two_norm(); | |
| 182 | 93628 | velocity *= flux / area; | |
| 183 | |||
| 184 | // TODO: Multiple throats connected to the same pore? | ||
| 185 | 93628 | return velocity; | |
| 186 | } | ||
| 187 | else | ||
| 188 | ✗ | return VelocityVector(0.0); | |
| 189 | } | ||
| 190 | |||
| 191 | /*! | ||
| 192 | * \brief Evaluate an advective flux across the interface and consider upwinding. | ||
| 193 | */ | ||
| 194 | 46038 | static Scalar advectiveFlux(const Scalar insideQuantity, const Scalar outsideQuantity, const Scalar volumeFlow, bool insideIsUpstream) | |
| 195 | { | ||
| 196 | 46038 | const Scalar upwindWeight = 1.0; //TODO use Flux.UpwindWeight or something like Coupling.UpwindWeight? | |
| 197 | |||
| 198 | 46038 | if(insideIsUpstream) | |
| 199 | 25110 | return (upwindWeight * insideQuantity + (1.0 - upwindWeight) * outsideQuantity) * volumeFlow; | |
| 200 | else | ||
| 201 | 20928 | return (upwindWeight * outsideQuantity + (1.0 - upwindWeight) * insideQuantity) * volumeFlow; | |
| 202 | } | ||
| 203 | |||
| 204 | /*! | ||
| 205 | * \brief Evaluate the conductive energy flux across the interface. | ||
| 206 | */ | ||
| 207 | template<std::size_t i, std::size_t j, bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0> | ||
| 208 | 15461 | static Scalar conductiveEnergyFlux(Dune::index_constant<i> domainI, | |
| 209 | Dune::index_constant<j> domainJ, | ||
| 210 | const SubControlVolumeFace<freeFlowMassIndex>& scvf, | ||
| 211 | const SubControlVolume<i>& scvI, | ||
| 212 | const SubControlVolume<j>& scvJ, | ||
| 213 | const VolumeVariables<i>& volVarsI, | ||
| 214 | const VolumeVariables<j>& volVarsJ) | ||
| 215 | { | ||
| 216 | // use properties (distance and thermal conductivity) for transimissibillity coefficient | ||
| 217 | // only from FF side as vertex (pore body) of PNM grid lies on boundary | ||
| 218 | 15461 | const auto& freeFlowVolVars = std::get<const VolumeVariables<freeFlowMassIndex>&>(std::forward_as_tuple(volVarsI, volVarsJ)); | |
| 219 | 15461 | const auto& ffScv = std::get<const SubControlVolume<freeFlowMassIndex>&>(std::forward_as_tuple(scvI, scvJ)); | |
| 220 | // distance from FF cell center to interface | ||
| 221 | 15461 | const Scalar distance = getDistance_(ffScv, scvf); | |
| 222 | 15461 | const Scalar tij = freeFlowVolVars.fluidThermalConductivity() / distance; | |
| 223 | |||
| 224 | 15461 | const Scalar deltaT = volVarsJ.temperature() - volVarsI.temperature(); | |
| 225 | |||
| 226 | 15461 | return -deltaT * tij; | |
| 227 | } | ||
| 228 | |||
| 229 | protected: | ||
| 230 | /*! | ||
| 231 | * \brief Returns the distance between an scvf and the corresponding scv center. | ||
| 232 | */ | ||
| 233 | template<class Scv, class Scvf> | ||
| 234 | 25174 | static Scalar getDistance_(const Scv& scv, const Scvf& scvf) | |
| 235 | { | ||
| 236 | 50348 | return (scv.dofPosition() - scvf.ipGlobal()).two_norm(); | |
| 237 | } | ||
| 238 | }; | ||
| 239 | |||
| 240 | /*! | ||
| 241 | * \ingroup FreeFlowPoreNetworkCoupling | ||
| 242 | * \brief Coupling conditions specialization for non-compositional models. | ||
| 243 | */ | ||
| 244 | template<class MDTraits, class CouplingManager, bool enableEnergyBalance> | ||
| 245 | class FreeFlowPoreNetworkCouplingConditionsImplementation<MDTraits, CouplingManager, enableEnergyBalance, false> | ||
| 246 | : public FreeFlowPoreNetworkCouplingConditionsImplementationBase<MDTraits, CouplingManager> | ||
| 247 | { | ||
| 248 | using ParentType = FreeFlowPoreNetworkCouplingConditionsImplementationBase<MDTraits, CouplingManager>; | ||
| 249 | using Scalar = typename MDTraits::Scalar; | ||
| 250 | |||
| 251 | // the sub domain type tags | ||
| 252 | template<std::size_t id> | ||
| 253 | using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag; | ||
| 254 | |||
| 255 | template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>; | ||
| 256 | template<std::size_t id> using Element = typename GridGeometry<id>::GridView::template Codim<0>::Entity; | ||
| 257 | template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView; | ||
| 258 | template<std::size_t id> using SubControlVolumeFace = typename GridGeometry<id>::LocalView::SubControlVolumeFace; | ||
| 259 | template<std::size_t id> using SubControlVolume = typename GridGeometry<id>::LocalView::SubControlVolume; | ||
| 260 | template<std::size_t id> using Indices = typename GetPropType<SubDomainTypeTag<id>, Properties::ModelTraits>::Indices; | ||
| 261 | template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView; | ||
| 262 | template<std::size_t id> using VolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::VolumeVariables; | ||
| 263 | |||
| 264 | static_assert(GetPropType<SubDomainTypeTag<ParentType::poreNetworkIndex>, Properties::ModelTraits>::numFluidComponents() == GetPropType<SubDomainTypeTag<ParentType::poreNetworkIndex>, Properties::ModelTraits>::numFluidPhases(), | ||
| 265 | "Pore-network model must not be compositional"); | ||
| 266 | |||
| 267 | public: | ||
| 268 | using ParentType::ParentType; | ||
| 269 | using ParentType::couplingPhaseIdx; | ||
| 270 | |||
| 271 | /*! | ||
| 272 | * \brief Returns the mass flux across the coupling boundary as seen from the pore-network domain. | ||
| 273 | */ | ||
| 274 | template<class CouplingContext> | ||
| 275 | 896 | static Scalar massCouplingCondition(Dune::index_constant<ParentType::poreNetworkIndex> domainI, | |
| 276 | Dune::index_constant<ParentType::freeFlowMassIndex> domainJ, | ||
| 277 | const FVElementGeometry<ParentType::poreNetworkIndex>& fvGeometry, | ||
| 278 | const SubControlVolume<ParentType::poreNetworkIndex>& scv, | ||
| 279 | const ElementVolumeVariables<ParentType::poreNetworkIndex>& insideVolVars, | ||
| 280 | const CouplingContext& context) | ||
| 281 | { | ||
| 282 | 896 | Scalar massFlux(0.0); | |
| 283 | 896 | const auto& pnmVolVars = insideVolVars[scv.indexInElement()]; | |
| 284 | |||
| 285 |
2/2✓ Branch 0 taken 4198 times.
✓ Branch 1 taken 896 times.
|
5094 | for (const auto& c : context) |
| 286 | { | ||
| 287 | // positive values of normal free-flow velocity indicate flux leaving the free flow into the pore-network region | ||
| 288 | 4198 | const Scalar normalFFVelocity = c.velocity * c.scvf.unitOuterNormal(); | |
| 289 | |||
| 290 | // normal pnm velocity (correspondign to its normal vector) is in the opposite direction of normal free flow velocity | ||
| 291 | // positive values of normal pnm velocity indicate flux leaving the pnm into the free-flow region | ||
| 292 |
2/2✓ Branch 0 taken 3753 times.
✓ Branch 1 taken 445 times.
|
4198 | const Scalar normalPNMVelocity = -normalFFVelocity; |
| 293 | 4198 | const bool pnmIsUpstream = std::signbit(normalFFVelocity); | |
| 294 |
2/2✓ Branch 0 taken 3753 times.
✓ Branch 1 taken 445 times.
|
4198 | const Scalar area = c.scvf.area() * c.volVars.extrusionFactor(); |
| 295 | |||
| 296 | 8396 | auto flux = massFlux_(domainI, domainJ, c.scvf, scv, c.scv, pnmVolVars, c.volVars, normalPNMVelocity, pnmIsUpstream); | |
| 297 | // flux is used as source term (volumetric flux): positive values mean influx | ||
| 298 | // thus, it is multiplied with area and we flip the sign | ||
| 299 | 4198 | flux *= area; | |
| 300 | 4198 | flux *= -1.0; | |
| 301 | |||
| 302 | 4198 | massFlux += flux; | |
| 303 | } | ||
| 304 | |||
| 305 | 896 | return massFlux; | |
| 306 | } | ||
| 307 | |||
| 308 | /*! | ||
| 309 | * \brief Returns the mass flux across the coupling boundary as seen from the free-flow domain. | ||
| 310 | */ | ||
| 311 | template<class CouplingContext> | ||
| 312 | 2113 | static Scalar massCouplingCondition(Dune::index_constant<ParentType::freeFlowMassIndex> domainI, | |
| 313 | Dune::index_constant<ParentType::poreNetworkIndex> domainJ, | ||
| 314 | const FVElementGeometry<ParentType::freeFlowMassIndex>& fvGeometry, | ||
| 315 | const SubControlVolumeFace<ParentType::freeFlowMassIndex>& scvf, | ||
| 316 | const ElementVolumeVariables<ParentType::freeFlowMassIndex>& insideVolVars, | ||
| 317 | const CouplingContext& context) | ||
| 318 | { | ||
| 319 | // positive values indicate flux into pore-network region | ||
| 320 | 4226 | const Scalar normalFFVelocity = context.velocity * scvf.unitOuterNormal(); | |
| 321 | 2113 | const bool ffIsUpstream = !std::signbit(normalFFVelocity); | |
| 322 | 2113 | const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); | |
| 323 | 2113 | const auto& ffVolVars = insideVolVars[scvf.insideScvIdx()]; | |
| 324 | |||
| 325 | // flux is used in Neumann condition: positive values mean flux out of free-flow domain | ||
| 326 | 4226 | return massFlux_(domainI, domainJ, scvf, insideScv, context.scv, ffVolVars, context.volVars, normalFFVelocity, ffIsUpstream); | |
| 327 | } | ||
| 328 | |||
| 329 | /*! | ||
| 330 | * \brief Returns the energy flux across the coupling boundary as seen from the pore network. | ||
| 331 | */ | ||
| 332 | template<class CouplingContext, bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0> | ||
| 333 | 660 | static Scalar energyCouplingCondition(Dune::index_constant<ParentType::poreNetworkIndex> domainI, | |
| 334 | Dune::index_constant<ParentType::freeFlowMassIndex> domainJ, | ||
| 335 | const FVElementGeometry<ParentType::poreNetworkIndex>& fvGeometry, | ||
| 336 | const SubControlVolume<ParentType::poreNetworkIndex>& scv, | ||
| 337 | const ElementVolumeVariables<ParentType::poreNetworkIndex>& insideVolVars, | ||
| 338 | const CouplingContext& context) | ||
| 339 | { | ||
| 340 | 660 | Scalar energyFlux(0.0); | |
| 341 | |||
| 342 | //use VolumeVariables (same type as for context.volVars) instead of ElementVolumeVariables | ||
| 343 | 660 | const auto& pnmVolVars = insideVolVars[scv.indexInElement()]; | |
| 344 | |||
| 345 |
2/2✓ Branch 0 taken 3300 times.
✓ Branch 1 taken 660 times.
|
3960 | for(const auto& c : context) |
| 346 | { | ||
| 347 | // positive values indicate flux into pore-network region | ||
| 348 | 3300 | const Scalar normalFFVelocity = c.velocity * c.scvf.unitOuterNormal(); | |
| 349 | 3300 | const bool pnmIsUpstream = std::signbit(normalFFVelocity); | |
| 350 | 3300 | const Scalar normalPNMVelocity = -normalFFVelocity; | |
| 351 | 3300 | const Scalar area = c.scvf.area() * c.volVars.extrusionFactor(); | |
| 352 | |||
| 353 | 3300 | auto flux = energyFlux_(domainI, domainJ, c.scvf, scv, c.scv, pnmVolVars, c.volVars, normalPNMVelocity, pnmIsUpstream); | |
| 354 | |||
| 355 | // flux is used as source term (volumetric flux): positive values mean influx | ||
| 356 | // thus, it is multiplied with area and we flip the sign | ||
| 357 | 3300 | flux *= area; | |
| 358 | 3300 | flux *= -1.0; | |
| 359 | |||
| 360 | 3300 | energyFlux += flux; | |
| 361 | } | ||
| 362 | |||
| 363 | 660 | return energyFlux; | |
| 364 | } | ||
| 365 | |||
| 366 | /*! | ||
| 367 | * \brief Returns the energy flux across the coupling boundary as seen from the free-flow domain. | ||
| 368 | */ | ||
| 369 | template<class CouplingContext, bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0> | ||
| 370 | 1540 | static Scalar energyCouplingCondition(Dune::index_constant<ParentType::freeFlowMassIndex> domainI, | |
| 371 | Dune::index_constant<ParentType::poreNetworkIndex> domainJ, | ||
| 372 | const FVElementGeometry<ParentType::freeFlowMassIndex>& fvGeometry, | ||
| 373 | const SubControlVolumeFace<ParentType::freeFlowMassIndex>& scvf, | ||
| 374 | const ElementVolumeVariables<ParentType::freeFlowMassIndex>& insideVolVars, | ||
| 375 | const CouplingContext& context) | ||
| 376 | { | ||
| 377 | // positive values indicate flux into pore-network region | ||
| 378 | 3080 | const Scalar normalFFVelocity = context.velocity * scvf.unitOuterNormal(); | |
| 379 | 1540 | const bool ffIsUpstream = !std::signbit(normalFFVelocity); | |
| 380 | 1540 | const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); | |
| 381 | //use VolumeVariables (same type as for context.volVars) instead of ElementVolumeVariables | ||
| 382 | 1540 | const auto& ffVolVars = insideVolVars[scvf.insideScvIdx()]; | |
| 383 | |||
| 384 | 1540 | return energyFlux_(domainI, domainJ, scvf, insideScv, context.scv, ffVolVars, context.volVars, normalFFVelocity, ffIsUpstream); | |
| 385 | } | ||
| 386 | |||
| 387 | private: | ||
| 388 | /*! | ||
| 389 | * \brief Evaluate the mole/mass flux across the interface. | ||
| 390 | */ | ||
| 391 | template<std::size_t i, std::size_t j> | ||
| 392 | 6311 | static Scalar massFlux_(Dune::index_constant<i> domainI, | |
| 393 | Dune::index_constant<j> domainJ, | ||
| 394 | const SubControlVolumeFace<ParentType::freeFlowMassIndex>& scvf, | ||
| 395 | const SubControlVolume<i>& scvI, | ||
| 396 | const SubControlVolume<j>& scvJ, | ||
| 397 | const VolumeVariables<i>& insideVolVars, | ||
| 398 | const VolumeVariables<j>& outsideVolVars, | ||
| 399 | const Scalar velocity, | ||
| 400 | const bool insideIsUpstream) | ||
| 401 | { | ||
| 402 |
2/2✓ Branch 0 taken 269 times.
✓ Branch 1 taken 1844 times.
|
2113 | Scalar flux(0.0); |
| 403 | |||
| 404 |
2/2✓ Branch 0 taken 3753 times.
✓ Branch 1 taken 445 times.
|
4198 | const Scalar insideDensity = insideVolVars.density(couplingPhaseIdx(domainI)); |
| 405 |
2/2✓ Branch 0 taken 269 times.
✓ Branch 1 taken 1844 times.
|
2113 | const Scalar outsideDensity = outsideVolVars.density(couplingPhaseIdx(domainJ)); |
| 406 | |||
| 407 |
4/4✓ Branch 0 taken 3359 times.
✓ Branch 1 taken 514 times.
✓ Branch 2 taken 663 times.
✓ Branch 3 taken 1775 times.
|
10509 | flux = ParentType::advectiveFlux(insideDensity, outsideDensity, velocity, insideIsUpstream); |
| 408 | |||
| 409 | return flux; | ||
| 410 | } | ||
| 411 | |||
| 412 | /*! | ||
| 413 | * \brief Evaluate the energy flux across the interface. | ||
| 414 | */ | ||
| 415 | template<std::size_t i, std::size_t j, bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0> | ||
| 416 | 4840 | static Scalar energyFlux_(Dune::index_constant<i> domainI, | |
| 417 | Dune::index_constant<j> domainJ, | ||
| 418 | const SubControlVolumeFace<ParentType::freeFlowMassIndex>& scvf, | ||
| 419 | const SubControlVolume<i>& scvI, | ||
| 420 | const SubControlVolume<j>& scvJ, | ||
| 421 | const VolumeVariables<i>& insideVolVars, | ||
| 422 | const VolumeVariables<j>& outsideVolVars, | ||
| 423 | const Scalar velocity, | ||
| 424 | const bool insideIsUpstream) | ||
| 425 | { | ||
| 426 | 4840 | Scalar flux(0.0); | |
| 427 | |||
| 428 | // convective fluxes | ||
| 429 | 4840 | const Scalar insideTerm = insideVolVars.density(couplingPhaseIdx(domainI)) * insideVolVars.enthalpy(couplingPhaseIdx(domainI)); | |
| 430 | 4840 | const Scalar outsideTerm = outsideVolVars.density(couplingPhaseIdx(domainJ)) * outsideVolVars.enthalpy(couplingPhaseIdx(domainJ)); | |
| 431 | |||
| 432 | 4840 | flux += ParentType::advectiveFlux(insideTerm, outsideTerm, velocity, insideIsUpstream); | |
| 433 | |||
| 434 | // conductive energy fluxes | ||
| 435 | 4840 | flux += ParentType::conductiveEnergyFlux(domainI, domainJ, scvf, scvI, scvJ, insideVolVars, outsideVolVars); | |
| 436 | |||
| 437 | 4840 | return flux; | |
| 438 | } | ||
| 439 | }; | ||
| 440 | |||
| 441 | /*! | ||
| 442 | * \ingroup FreeFlowPoreNetworkCoupling | ||
| 443 | * \brief Coupling conditions specialization for compositional models. | ||
| 444 | */ | ||
| 445 | template<class MDTraits, class CouplingManager, bool enableEnergyBalance> | ||
| 446 | class FreeFlowPoreNetworkCouplingConditionsImplementation<MDTraits, CouplingManager, enableEnergyBalance, true> | ||
| 447 | : public FreeFlowPoreNetworkCouplingConditionsImplementationBase<MDTraits, CouplingManager> | ||
| 448 | { | ||
| 449 | using ParentType = FreeFlowPoreNetworkCouplingConditionsImplementationBase<MDTraits, CouplingManager>; | ||
| 450 | using Scalar = typename MDTraits::Scalar; | ||
| 451 | |||
| 452 | // the sub domain type tags | ||
| 453 | template<std::size_t id> | ||
| 454 | using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag; | ||
| 455 | |||
| 456 | template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>; | ||
| 457 | template<std::size_t id> using Element = typename GridGeometry<id>::GridView::template Codim<0>::Entity; | ||
| 458 | template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView; | ||
| 459 | template<std::size_t id> using SubControlVolumeFace = typename GridGeometry<id>::LocalView::SubControlVolumeFace; | ||
| 460 | template<std::size_t id> using SubControlVolume = typename GridGeometry<id>::LocalView::SubControlVolume; | ||
| 461 | template<std::size_t id> using Indices = typename GetPropType<SubDomainTypeTag<id>, Properties::ModelTraits>::Indices; | ||
| 462 | template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView; | ||
| 463 | template<std::size_t id> using VolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::VolumeVariables; | ||
| 464 | template<std::size_t id> using FluidSystem = typename VolumeVariables<id>::FluidSystem; | ||
| 465 | template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>; | ||
| 466 | template<std::size_t id> using NumEqVector = typename Problem<id>::Traits::NumEqVector; | ||
| 467 | |||
| 468 | static constexpr bool useMoles = GetPropType<SubDomainTypeTag<ParentType::freeFlowMassIndex>, Properties::ModelTraits>::useMoles(); | ||
| 469 | static constexpr auto numComponents = GetPropType<SubDomainTypeTag<ParentType::freeFlowMassIndex>, Properties::ModelTraits>::numFluidComponents(); | ||
| 470 | static constexpr auto referenceSystemFormulation = GetPropType<SubDomainTypeTag<ParentType::freeFlowMassIndex>, Properties::MolecularDiffusionType>::referenceSystemFormulation(); | ||
| 471 | static constexpr auto replaceCompEqIdx = GetPropType<SubDomainTypeTag<ParentType::freeFlowMassIndex>, Properties::ModelTraits>::replaceCompEqIdx(); | ||
| 472 | |||
| 473 | static_assert(GetPropType<SubDomainTypeTag<ParentType::poreNetworkIndex>, Properties::ModelTraits>::numFluidComponents() == numComponents, | ||
| 474 | "Models must have same number of components"); | ||
| 475 | |||
| 476 | public: | ||
| 477 | using ParentType::ParentType; | ||
| 478 | using ParentType::couplingPhaseIdx; | ||
| 479 | using ParentType::couplingCompIdx; | ||
| 480 | |||
| 481 | using NumCompVector = Dune::FieldVector<Scalar, numComponents>; | ||
| 482 | |||
| 483 | /*! | ||
| 484 | * \brief Returns the mass flux across the coupling boundary as seen from the pore-network domain. | ||
| 485 | */ | ||
| 486 | template<class CouplingContext> | ||
| 487 | 2895 | static NumCompVector massCouplingCondition(Dune::index_constant<ParentType::poreNetworkIndex> domainI, | |
| 488 | Dune::index_constant<ParentType::freeFlowMassIndex> domainJ, | ||
| 489 | const FVElementGeometry<ParentType::poreNetworkIndex>& fvGeometry, | ||
| 490 | const SubControlVolume<ParentType::poreNetworkIndex>& scv, | ||
| 491 | const ElementVolumeVariables<ParentType::poreNetworkIndex>& insideVolVars, | ||
| 492 | const CouplingContext& context) | ||
| 493 | { | ||
| 494 | 2895 | NumCompVector massFlux(0.0); | |
| 495 | 2895 | const auto& pnmVolVars = insideVolVars[scv.indexInElement()]; | |
| 496 | |||
| 497 |
2/2✓ Branch 0 taken 8685 times.
✓ Branch 1 taken 2895 times.
|
11580 | for (const auto& c : context) |
| 498 | { | ||
| 499 | // positive values indicate flux into pore-network region | ||
| 500 | 8685 | const Scalar normalFFVelocity = c.velocity * c.scvf.unitOuterNormal(); | |
| 501 | 8685 | const bool pnmIsUpstream = std::signbit(normalFFVelocity); | |
| 502 | 8685 | const Scalar normalPNMVelocity = -normalFFVelocity; | |
| 503 | 8685 | const Scalar area = c.scvf.area() * c.volVars.extrusionFactor(); | |
| 504 | |||
| 505 | 8685 | auto flux = massFlux_(domainI, domainJ, c.scvf, scv, c.scv, pnmVolVars, c.volVars, normalPNMVelocity, pnmIsUpstream); | |
| 506 | |||
| 507 | // flux is used as source term (volumetric flux): positive values mean influx | ||
| 508 | // thus, it is multiplied with area and we flip the sign | ||
| 509 |
2/2✓ Branch 0 taken 17370 times.
✓ Branch 1 taken 8685 times.
|
26055 | flux *= area; |
| 510 |
2/2✓ Branch 0 taken 17370 times.
✓ Branch 1 taken 8685 times.
|
26055 | flux *= -1.0; |
| 511 | |||
| 512 | 8685 | massFlux += flux; | |
| 513 | } | ||
| 514 | |||
| 515 | 2895 | return massFlux; | |
| 516 | } | ||
| 517 | |||
| 518 | /*! | ||
| 519 | * \brief Returns the mass flux across the coupling boundary as seen from the free-flow domain. | ||
| 520 | */ | ||
| 521 | template<class CouplingContext> | ||
| 522 | 5868 | static NumCompVector massCouplingCondition(Dune::index_constant<ParentType::freeFlowMassIndex> domainI, | |
| 523 | Dune::index_constant<ParentType::poreNetworkIndex> domainJ, | ||
| 524 | const FVElementGeometry<ParentType::freeFlowMassIndex>& fvGeometry, | ||
| 525 | const SubControlVolumeFace<ParentType::freeFlowMassIndex>& scvf, | ||
| 526 | const ElementVolumeVariables<ParentType::freeFlowMassIndex>& insideVolVars, | ||
| 527 | const CouplingContext& context) | ||
| 528 | { | ||
| 529 | // positive values indicate flux into pore-network region | ||
| 530 | 11736 | const Scalar normalFFVelocity = context.velocity * scvf.unitOuterNormal(); | |
| 531 | 5868 | const bool ffIsUpstream = !std::signbit(normalFFVelocity); | |
| 532 | 5868 | const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); | |
| 533 | 5868 | const auto& ffVolVars = insideVolVars[scvf.insideScvIdx()]; | |
| 534 | |||
| 535 | 5868 | return massFlux_(domainI, domainJ, scvf, insideScv, context.scv, ffVolVars, context.volVars, normalFFVelocity, ffIsUpstream); | |
| 536 | } | ||
| 537 | |||
| 538 | /*! | ||
| 539 | * \brief Returns the energy flux across the coupling boundary as seen from the pore network. | ||
| 540 | */ | ||
| 541 | template<class CouplingContext, bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0> | ||
| 542 | 1175 | static Scalar energyCouplingCondition(Dune::index_constant<ParentType::poreNetworkIndex> domainI, | |
| 543 | Dune::index_constant<ParentType::freeFlowMassIndex> domainJ, | ||
| 544 | const FVElementGeometry<ParentType::poreNetworkIndex>& fvGeometry, | ||
| 545 | const SubControlVolume<ParentType::poreNetworkIndex>& scv, | ||
| 546 | const ElementVolumeVariables<ParentType::poreNetworkIndex>& insideVolVars, | ||
| 547 | const CouplingContext& context) | ||
| 548 | { | ||
| 549 | 1175 | Scalar energyFlux(0.0); | |
| 550 | |||
| 551 | //use VolumeVariables (same type as for context.volVars) instead of ElementVolumeVariables | ||
| 552 | 1175 | const auto& pnmVolVars = insideVolVars[scv.indexInElement()]; | |
| 553 | |||
| 554 |
2/2✓ Branch 0 taken 3525 times.
✓ Branch 1 taken 1175 times.
|
4700 | for(const auto& c : context) |
| 555 | { | ||
| 556 | // positive values indicate flux into pore-network region | ||
| 557 | 3525 | const Scalar normalFFVelocity = c.velocity * c.scvf.unitOuterNormal(); | |
| 558 | 3525 | const bool pnmIsUpstream = std::signbit(normalFFVelocity); | |
| 559 | 3525 | const Scalar normalPNMVelocity = -normalFFVelocity; | |
| 560 | 3525 | const Scalar area = c.scvf.area() * c.volVars.extrusionFactor(); | |
| 561 | |||
| 562 | 3525 | auto flux = energyFlux_(domainI, domainJ, c.scvf, scv, c.scv, pnmVolVars, c.volVars, normalPNMVelocity, pnmIsUpstream); | |
| 563 | |||
| 564 | // flux is used as source term (volumetric flux): positive values mean influx | ||
| 565 | // thus, it is multiplied with area and we flip the sign | ||
| 566 | 3525 | flux *= area; | |
| 567 | 3525 | flux *= -1.0; | |
| 568 | 3525 | energyFlux += flux; | |
| 569 | } | ||
| 570 | |||
| 571 | 1175 | return energyFlux; | |
| 572 | } | ||
| 573 | |||
| 574 | /*! | ||
| 575 | * \brief Returns the energy flux across the coupling boundary as seen from the free-flow domain. | ||
| 576 | */ | ||
| 577 | template<class CouplingContext, bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0> | ||
| 578 | 2256 | static Scalar energyCouplingCondition(Dune::index_constant<ParentType::freeFlowMassIndex> domainI, | |
| 579 | Dune::index_constant<ParentType::poreNetworkIndex> domainJ, | ||
| 580 | const FVElementGeometry<ParentType::freeFlowMassIndex>& fvGeometry, | ||
| 581 | const SubControlVolumeFace<ParentType::freeFlowMassIndex>& scvf, | ||
| 582 | const ElementVolumeVariables<ParentType::freeFlowMassIndex>& insideVolVars, | ||
| 583 | const CouplingContext& context) | ||
| 584 | { | ||
| 585 | // positive values indicate flux into pore-network region | ||
| 586 | 4512 | const Scalar normalFFVelocity = context.velocity * scvf.unitOuterNormal(); | |
| 587 | 2256 | const bool ffIsUpstream = !std::signbit(normalFFVelocity); | |
| 588 | 2256 | const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); | |
| 589 | //use VolumeVariables (same type as for context.volVars) instead of ElementVolumeVariables | ||
| 590 | 2256 | const auto& ffVolVars = insideVolVars[scvf.insideScvIdx()]; | |
| 591 | |||
| 592 | 2256 | return energyFlux_(domainI, domainJ, scvf, insideScv, context.scv, ffVolVars, context.volVars, normalFFVelocity, ffIsUpstream); | |
| 593 | } | ||
| 594 | |||
| 595 | private: | ||
| 596 | /*! | ||
| 597 | * \brief Evaluate the compositional mole/mass flux across the interface. | ||
| 598 | */ | ||
| 599 | template<std::size_t i, std::size_t j> | ||
| 600 | 29106 | static NumCompVector massFlux_(Dune::index_constant<i> domainI, | |
| 601 | Dune::index_constant<j> domainJ, | ||
| 602 | const SubControlVolumeFace<ParentType::freeFlowMassIndex>& scvf, | ||
| 603 | const SubControlVolume<i>& scvI, | ||
| 604 | const SubControlVolume<j>& scvJ, | ||
| 605 | const VolumeVariables<i>& insideVolVars, | ||
| 606 | const VolumeVariables<j>& outsideVolVars, | ||
| 607 | const Scalar velocity, | ||
| 608 | const bool insideIsUpstream) | ||
| 609 | { | ||
| 610 | 29106 | NumCompVector flux(0.0); | |
| 611 | |||
| 612 | 87318 | auto moleOrMassFraction = [&](const auto& volVars, int phaseIdx, int compIdx) | |
| 613 | 116424 | { return useMoles ? volVars.moleFraction(phaseIdx, compIdx) : volVars.massFraction(phaseIdx, compIdx); }; | |
| 614 | |||
| 615 |
2/2✓ Branch 0 taken 5344 times.
✓ Branch 1 taken 6392 times.
|
87318 | auto moleOrMassDensity = [&](const auto& volVars, int phaseIdx) |
| 616 |
2/2✓ Branch 0 taken 5344 times.
✓ Branch 1 taken 6392 times.
|
87318 | { return useMoles ? volVars.molarDensity(phaseIdx) : volVars.density(phaseIdx); }; |
| 617 | |||
| 618 | // advective fluxes | ||
| 619 | 69948 | auto insideTerm = [&](int compIdx) | |
| 620 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 17370 times.
|
40842 | { return moleOrMassFraction(insideVolVars, couplingPhaseIdx(domainI), compIdx) * moleOrMassDensity(insideVolVars, couplingPhaseIdx(domainI)); }; |
| 621 | |||
| 622 | 75582 | auto outsideTerm = [&](int compIdx) | |
| 623 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 11736 times.
|
46476 | { return moleOrMassFraction(outsideVolVars, couplingPhaseIdx(domainJ), compIdx) * moleOrMassDensity(outsideVolVars, couplingPhaseIdx(domainJ)); }; |
| 624 | |||
| 625 | |||
| 626 |
2/2✓ Branch 0 taken 29106 times.
✓ Branch 1 taken 14553 times.
|
87318 | for (int compIdx = 0; compIdx < numComponents; ++compIdx) |
| 627 | { | ||
| 628 | 58212 | const int domainICompIdx = couplingCompIdx(domainI, compIdx); | |
| 629 |
2/2✓ Branch 0 taken 14553 times.
✓ Branch 1 taken 14553 times.
|
58212 | const int domainJCompIdx = couplingCompIdx(domainJ, compIdx); |
| 630 | |||
| 631 |
3/4✓ Branch 0 taken 14553 times.
✓ Branch 1 taken 14553 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 29106 times.
|
174636 | assert(FluidSystem<i>::componentName(domainICompIdx) == FluidSystem<j>::componentName(domainJCompIdx)); |
| 632 | |||
| 633 |
2/2✓ Branch 1 taken 14784 times.
✓ Branch 2 taken 14322 times.
|
116424 | flux[domainICompIdx] += ParentType::advectiveFlux(insideTerm(domainICompIdx), outsideTerm(domainJCompIdx), velocity, insideIsUpstream); |
| 634 | } | ||
| 635 | |||
| 636 | // diffusive fluxes | ||
| 637 | 29106 | NumCompVector diffusiveFlux = diffusiveMolecularFlux_(domainI, domainJ, scvf, scvI, scvJ, insideVolVars, outsideVolVars); | |
| 638 | |||
| 639 | //convert to correct units if necessary | ||
| 640 | if (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged && useMoles) | ||
| 641 | { | ||
| 642 |
2/2✓ Branch 0 taken 29106 times.
✓ Branch 1 taken 14553 times.
|
87318 | for (int compIdx = 0; compIdx < numComponents; ++compIdx) |
| 643 | { | ||
| 644 |
2/2✓ Branch 0 taken 14553 times.
✓ Branch 1 taken 14553 times.
|
58212 | const int domainICompIdx = couplingCompIdx(domainI, compIdx); |
| 645 | 58212 | diffusiveFlux[domainICompIdx] /= FluidSystem<i>::molarMass(domainICompIdx); | |
| 646 | } | ||
| 647 | } | ||
| 648 | if (referenceSystemFormulation == ReferenceSystemFormulation::molarAveraged && !useMoles) | ||
| 649 | { | ||
| 650 | for (int compIdx = 0; compIdx < numComponents; ++compIdx) | ||
| 651 | { | ||
| 652 | const int domainICompIdx = couplingCompIdx(domainI, compIdx); | ||
| 653 | diffusiveFlux[domainICompIdx] *= FluidSystem<i>::molarMass(domainICompIdx); | ||
| 654 | } | ||
| 655 | } | ||
| 656 | |||
| 657 | // total flux | ||
| 658 | 29106 | flux += diffusiveFlux; | |
| 659 | |||
| 660 | // convert to total mass/mole balance, if set be user | ||
| 661 | if (replaceCompEqIdx < numComponents) | ||
| 662 | flux[replaceCompEqIdx] = std::accumulate(flux.begin(), flux.end(), 0.0); | ||
| 663 | |||
| 664 | 29106 | return flux; | |
| 665 | } | ||
| 666 | |||
| 667 | /*! | ||
| 668 | * \brief Evaluate the energy flux (convective and conductive) across the interface. | ||
| 669 | */ | ||
| 670 | template<std::size_t i, std::size_t j, bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0> | ||
| 671 | 11562 | static Scalar energyFlux_(Dune::index_constant<i> domainI, | |
| 672 | Dune::index_constant<j> domainJ, | ||
| 673 | const SubControlVolumeFace<ParentType::freeFlowMassIndex>& scvf, | ||
| 674 | const SubControlVolume<i>& scvI, | ||
| 675 | const SubControlVolume<j>& scvJ, | ||
| 676 | const VolumeVariables<i>& insideVolVars, | ||
| 677 | const VolumeVariables<j>& outsideVolVars, | ||
| 678 | const Scalar velocity, | ||
| 679 | const bool insideIsUpstream) | ||
| 680 | { | ||
| 681 |
2/2✓ Branch 0 taken 781 times.
✓ Branch 1 taken 1475 times.
|
11562 | Scalar flux(0.0); |
| 682 | |||
| 683 | // convective fluxes | ||
| 684 |
2/2✓ Branch 0 taken 3084 times.
✓ Branch 1 taken 2697 times.
|
11562 | const Scalar insideTerm = insideVolVars.density(couplingPhaseIdx(domainI)) * insideVolVars.enthalpy(couplingPhaseIdx(domainI)); |
| 685 |
2/2✓ Branch 0 taken 3084 times.
✓ Branch 1 taken 2697 times.
|
11562 | const Scalar outsideTerm = outsideVolVars.density(couplingPhaseIdx(domainJ)) * outsideVolVars.enthalpy(couplingPhaseIdx(domainJ)); |
| 686 | |||
| 687 |
2/2✓ Branch 0 taken 3084 times.
✓ Branch 1 taken 2697 times.
|
23124 | flux += ParentType::advectiveFlux(insideTerm, outsideTerm, velocity, insideIsUpstream); |
| 688 | |||
| 689 | // conductive fluxes | ||
| 690 | 11562 | flux += ParentType::conductiveEnergyFlux(domainI, domainJ, scvf, scvI, scvJ, insideVolVars, outsideVolVars); | |
| 691 | |||
| 692 | // TODO: add contribution from diffusive fluxes | ||
| 693 | |||
| 694 | 11562 | return flux; | |
| 695 | } | ||
| 696 | |||
| 697 | /*! | ||
| 698 | * \brief Evaluate the diffusive mole/mass flux across the interface. | ||
| 699 | */ | ||
| 700 | template<std::size_t i, std::size_t j> | ||
| 701 | 14553 | static NumCompVector diffusiveMolecularFlux_(Dune::index_constant<i> domainI, | |
| 702 | Dune::index_constant<j> domainJ, | ||
| 703 | const SubControlVolumeFace<ParentType::freeFlowMassIndex>& scvf, | ||
| 704 | const SubControlVolume<i>& scvI, | ||
| 705 | const SubControlVolume<j>& scvJ, | ||
| 706 | const VolumeVariables<i>& volVarsI, | ||
| 707 | const VolumeVariables<j>& volVarsJ) | ||
| 708 | { | ||
| 709 | 14553 | NumCompVector diffusiveFlux(0.0); | |
| 710 | 14553 | const Scalar avgDensity = 0.5*(massOrMolarDensity(volVarsI, referenceSystemFormulation, couplingPhaseIdx(domainI)) | |
| 711 | 14553 | + massOrMolarDensity(volVarsJ, referenceSystemFormulation, couplingPhaseIdx(domainJ))); | |
| 712 | |||
| 713 | 14553 | const auto& freeFlowVolVars = std::get<const VolumeVariables<ParentType::freeFlowMassIndex>&>(std::forward_as_tuple(volVarsI, volVarsJ)); | |
| 714 | 14553 | const auto& ffScv = std::get<const SubControlVolume<ParentType::freeFlowMassIndex>&>(std::forward_as_tuple(scvI, scvJ)); | |
| 715 | 14553 | const Scalar distance = ParentType::getDistance_(ffScv, scvf); | |
| 716 | |||
| 717 | 29106 | for (int compIdx = 1; compIdx < numComponents; ++compIdx) | |
| 718 | { | ||
| 719 | 14553 | const int freeFlowMainCompIdx = couplingPhaseIdx(ParentType::freeFlowMassIndex); | |
| 720 | 14553 | const int domainICompIdx = couplingCompIdx(domainI, compIdx); | |
| 721 | 14553 | const int domainJCompIdx = couplingCompIdx(domainJ, compIdx); | |
| 722 | |||
| 723 | 29106 | assert(FluidSystem<i>::componentName(domainICompIdx) == FluidSystem<j>::componentName(domainJCompIdx)); | |
| 724 | |||
| 725 | 14553 | const Scalar massOrMoleFractionI = massOrMoleFraction(volVarsI, referenceSystemFormulation, couplingPhaseIdx(domainI), domainICompIdx); | |
| 726 | 14553 | const Scalar massOrMoleFractionJ = massOrMoleFraction(volVarsJ, referenceSystemFormulation, couplingPhaseIdx(domainJ), domainJCompIdx); | |
| 727 | 14553 | const Scalar deltaMassOrMoleFrac = massOrMoleFractionJ - massOrMoleFractionI; | |
| 728 | |||
| 729 | 14553 | const Scalar tij = freeFlowVolVars.effectiveDiffusionCoefficient(couplingPhaseIdx(ParentType::freeFlowMassIndex), | |
| 730 | freeFlowMainCompIdx, | ||
| 731 | couplingCompIdx(ParentType::freeFlowMassIndex, compIdx)) | ||
| 732 | / distance; | ||
| 733 | 14553 | diffusiveFlux[domainICompIdx] += -avgDensity * tij * deltaMassOrMoleFrac; | |
| 734 | } | ||
| 735 | |||
| 736 | 14553 | const Scalar cumulativeFlux = std::accumulate(diffusiveFlux.begin(), diffusiveFlux.end(), 0.0); | |
| 737 | 14553 | diffusiveFlux[couplingCompIdx(domainI, 0)] = -cumulativeFlux; | |
| 738 | |||
| 739 | 14553 | return diffusiveFlux; | |
| 740 | } | ||
| 741 | |||
| 742 | static Scalar getComponentEnthalpy_(const VolumeVariables<ParentType::freeFlowMassIndex>& volVars, int phaseIdx, int compIdx) | ||
| 743 | { | ||
| 744 | return FluidSystem<ParentType::freeFlowMassIndex>::componentEnthalpy(volVars.fluidState(), 0, compIdx); | ||
| 745 | } | ||
| 746 | |||
| 747 | static Scalar getComponentEnthalpy_(const VolumeVariables<ParentType::poreNetworkIndex>& volVars, int phaseIdx, int compIdx) | ||
| 748 | { | ||
| 749 | return FluidSystem<ParentType::poreNetworkIndex>::componentEnthalpy(volVars.fluidState(), phaseIdx, compIdx); | ||
| 750 | } | ||
| 751 | }; | ||
| 752 | |||
| 753 | } // end namespace Dumux | ||
| 754 | |||
| 755 | #endif | ||
| 756 |