From 351ed2c1bdec653c700d15f2ecb40ebd47b7fa68 Mon Sep 17 00:00:00 2001
From: Anna Mareike Kostelecky <anmako96@web.de>
Date: Wed, 26 Feb 2025 14:21:09 +0100
Subject: [PATCH] [md][boundary][freeflowporenetwork] add energy and
compositional coupling conditions
---
.../freeflowporenetwork/couplingconditions.hh | 480 +++++++++++++++++-
.../freeflowporenetwork/couplingmanager.hh | 64 ++-
2 files changed, 512 insertions(+), 32 deletions(-)
diff --git a/dumux/multidomain/boundary/freeflowporenetwork/couplingconditions.hh b/dumux/multidomain/boundary/freeflowporenetwork/couplingconditions.hh
index 000e79d6f6..aade0712f0 100644
--- a/dumux/multidomain/boundary/freeflowporenetwork/couplingconditions.hh
+++ b/dumux/multidomain/boundary/freeflowporenetwork/couplingconditions.hh
@@ -16,6 +16,7 @@
#include <numeric>
#include <dune/common/exceptions.hh>
+#include <dune/common/fvector.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/math.hh>
@@ -39,8 +40,8 @@ template<class MDTraits, class CouplingManager>
using FreeFlowPoreNetworkCouplingConditions
= FreeFlowPoreNetworkCouplingConditionsImplementation<
MDTraits, CouplingManager,
- GetPropType<typename MDTraits::template SubDomain<0>::TypeTag, Properties::ModelTraits>::enableEnergyBalance(),
- (GetPropType<typename MDTraits::template SubDomain<0>::TypeTag, Properties::ModelTraits>::numFluidComponents() > 1)
+ GetPropType<typename MDTraits::template SubDomain<1>::TypeTag, Properties::ModelTraits>::enableEnergyBalance(),
+ (GetPropType<typename MDTraits::template SubDomain<1>::TypeTag, Properties::ModelTraits>::numFluidComponents() > 1)
>;
/*!
@@ -105,7 +106,6 @@ public:
static constexpr auto couplingCompIdx(Dune::index_constant<i> id, int coupledCompIdx)
{ return IndexHelper::couplingCompIdx(id, coupledCompIdx); }
-
/*!
* \brief Returns the momentum flux across the coupling boundary.
*
@@ -201,13 +201,37 @@ public:
return (upwindWeight * outsideQuantity + (1.0 - upwindWeight) * insideQuantity) * volumeFlow;
}
-protected:
+ /*!
+ * \brief Evaluate the conductive energy flux across the interface.
+ */
+ template<std::size_t i, std::size_t j, bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0>
+ static Scalar conductiveEnergyFlux(Dune::index_constant<i> domainI,
+ Dune::index_constant<j> domainJ,
+ const SubControlVolumeFace<freeFlowMassIndex>& scvf,
+ const SubControlVolume<i>& scvI,
+ const SubControlVolume<j>& scvJ,
+ const VolumeVariables<i>& volVarsI,
+ const VolumeVariables<j>& volVarsJ)
+ {
+ // use properties (distance and thermal conductivity) for transimissibillity coefficient
+ // only from FF side as vertex (pore body) of PNM grid lies on boundary
+ const auto& freeFlowVolVars = std::get<const VolumeVariables<freeFlowMassIndex>&>(std::forward_as_tuple(volVarsI, volVarsJ));
+ const auto& ffScv = std::get<const SubControlVolume<freeFlowMassIndex>&>(std::forward_as_tuple(scvI, scvJ));
+ // distance from FF cell center to interface
+ const Scalar distance = getDistance_(ffScv, scvf);
+ const Scalar tij = freeFlowVolVars.fluidThermalConductivity() / distance;
+
+ const Scalar deltaT = volVarsJ.temperature() - volVarsI.temperature();
+ return -deltaT * tij;
+ }
+
+protected:
/*!
* \brief Returns the distance between an scvf and the corresponding scv center.
*/
template<class Scv, class Scvf>
- Scalar getDistance_(const Scv& scv, const Scvf& scvf) const
+ static Scalar getDistance_(const Scv& scv, const Scvf& scvf)
{
return (scv.dofPosition() - scvf.ipGlobal()).two_norm();
}
@@ -256,19 +280,26 @@ public:
const CouplingContext& context)
{
Scalar massFlux(0.0);
+ const auto& pnmVolVars = insideVolVars[scv.indexInElement()];
for (const auto& c : context)
{
- // positive values indicate flux into pore-network region
+ // positive values of normal free-flow velocity indicate flux leaving the free flow into the pore-network region
const Scalar normalFFVelocity = c.velocity * c.scvf.unitOuterNormal();
- const bool pnmIsUpstream = std::signbit(normalFFVelocity);
- const Scalar pnmDensity = insideVolVars[scv].density(couplingPhaseIdx(ParentType::poreNetworkIndex));
- const Scalar ffDensity = c.volVars.density(couplingPhaseIdx(ParentType::freeFlowMassIndex));
+ // normal pnm velocity (correspondign to its normal vector) is in the opposite direction of normal free flow velocity
+ // positive values of normal pnm velocity indicate flux leaving the pnm into the free-flow region
+ const Scalar normalPNMVelocity = -normalFFVelocity;
+ const bool pnmIsUpstream = std::signbit(normalFFVelocity);
const Scalar area = c.scvf.area() * c.volVars.extrusionFactor();
- // flux is used as source term: positive values mean influx
- massFlux += ParentType::advectiveFlux(pnmDensity, ffDensity, normalFFVelocity*area, pnmIsUpstream);
+ auto flux = massFlux_(domainI, domainJ, c.scvf, scv, c.scv, pnmVolVars, c.volVars, normalPNMVelocity, pnmIsUpstream);
+ // flux is used as source term (volumetric flux): positive values mean influx
+ // thus, it is multiplied with area and we flip the sign
+ flux *= area;
+ flux *= -1.0;
+
+ massFlux += flux;
}
return massFlux;
@@ -288,16 +319,435 @@ public:
// positive values indicate flux into pore-network region
const Scalar normalFFVelocity = context.velocity * scvf.unitOuterNormal();
const bool ffIsUpstream = !std::signbit(normalFFVelocity);
-
- const Scalar ffDensity = insideVolVars[scvf.insideScvIdx()].density(couplingPhaseIdx(ParentType::freeFlowMassIndex));
- const Scalar pnmDensity = context.volVars.density(couplingPhaseIdx(ParentType::poreNetworkIndex));
+ const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+ const auto& ffVolVars = insideVolVars[scvf.insideScvIdx()];
// flux is used in Neumann condition: positive values mean flux out of free-flow domain
- return ParentType::advectiveFlux(ffDensity, pnmDensity, normalFFVelocity, ffIsUpstream);
+ return massFlux_(domainI, domainJ, scvf, insideScv, context.scv, ffVolVars, context.volVars, normalFFVelocity, ffIsUpstream);
+ }
+
+ /*!
+ * \brief Returns the energy flux across the coupling boundary as seen from the pore network.
+ */
+ template<class CouplingContext, bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0>
+ static Scalar energyCouplingCondition(Dune::index_constant<ParentType::poreNetworkIndex> domainI,
+ Dune::index_constant<ParentType::freeFlowMassIndex> domainJ,
+ const FVElementGeometry<ParentType::poreNetworkIndex>& fvGeometry,
+ const SubControlVolume<ParentType::poreNetworkIndex>& scv,
+ const ElementVolumeVariables<ParentType::poreNetworkIndex>& insideVolVars,
+ const CouplingContext& context)
+ {
+ Scalar energyFlux(0.0);
+
+ //use VolumeVariables (same type as for context.volVars) instead of ElementVolumeVariables
+ const auto& pnmVolVars = insideVolVars[scv.indexInElement()];
+
+ for(const auto& c : context)
+ {
+ // positive values indicate flux into pore-network region
+ const Scalar normalFFVelocity = c.velocity * c.scvf.unitOuterNormal();
+ const bool pnmIsUpstream = std::signbit(normalFFVelocity);
+ const Scalar normalPNMVelocity = -normalFFVelocity;
+ const Scalar area = c.scvf.area() * c.volVars.extrusionFactor();
+
+ auto flux = energyFlux_(domainI, domainJ, c.scvf, scv, c.scv, pnmVolVars, c.volVars, normalPNMVelocity, pnmIsUpstream);
+
+ // flux is used as source term (volumetric flux): positive values mean influx
+ // thus, it is multiplied with area and we flip the sign
+ flux *= area;
+ flux *= -1.0;
+
+ energyFlux += flux;
+ }
+
+ return energyFlux;
+ }
+
+ /*!
+ * \brief Returns the energy flux across the coupling boundary as seen from the free-flow domain.
+ */
+ template<class CouplingContext, bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0>
+ static Scalar energyCouplingCondition(Dune::index_constant<ParentType::freeFlowMassIndex> domainI,
+ Dune::index_constant<ParentType::poreNetworkIndex> domainJ,
+ const FVElementGeometry<ParentType::freeFlowMassIndex>& fvGeometry,
+ const SubControlVolumeFace<ParentType::freeFlowMassIndex>& scvf,
+ const ElementVolumeVariables<ParentType::freeFlowMassIndex>& insideVolVars,
+ const CouplingContext& context)
+ {
+ // positive values indicate flux into pore-network region
+ const Scalar normalFFVelocity = context.velocity * scvf.unitOuterNormal();
+ const bool ffIsUpstream = !std::signbit(normalFFVelocity);
+ const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+ //use VolumeVariables (same type as for context.volVars) instead of ElementVolumeVariables
+ const auto& ffVolVars = insideVolVars[scvf.insideScvIdx()];
+
+ return energyFlux_(domainI, domainJ, scvf, insideScv, context.scv, ffVolVars, context.volVars, normalFFVelocity, ffIsUpstream);
+ }
+
+private:
+ /*!
+ * \brief Evaluate the mole/mass flux across the interface.
+ */
+ template<std::size_t i, std::size_t j>
+ static Scalar massFlux_(Dune::index_constant<i> domainI,
+ Dune::index_constant<j> domainJ,
+ const SubControlVolumeFace<ParentType::freeFlowMassIndex>& scvf,
+ const SubControlVolume<i>& scvI,
+ const SubControlVolume<j>& scvJ,
+ const VolumeVariables<i>& insideVolVars,
+ const VolumeVariables<j>& outsideVolVars,
+ const Scalar velocity,
+ const bool insideIsUpstream)
+ {
+ Scalar flux(0.0);
+
+ const Scalar insideDensity = insideVolVars.density(couplingPhaseIdx(domainI));
+ const Scalar outsideDensity = outsideVolVars.density(couplingPhaseIdx(domainJ));
+
+ flux = ParentType::advectiveFlux(insideDensity, outsideDensity, velocity, insideIsUpstream);
+
+ return flux;
+ }
+
+ /*!
+ * \brief Evaluate the energy flux across the interface.
+ */
+ template<std::size_t i, std::size_t j, bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0>
+ static Scalar energyFlux_(Dune::index_constant<i> domainI,
+ Dune::index_constant<j> domainJ,
+ const SubControlVolumeFace<ParentType::freeFlowMassIndex>& scvf,
+ const SubControlVolume<i>& scvI,
+ const SubControlVolume<j>& scvJ,
+ const VolumeVariables<i>& insideVolVars,
+ const VolumeVariables<j>& outsideVolVars,
+ const Scalar velocity,
+ const bool insideIsUpstream)
+ {
+ Scalar flux(0.0);
+
+ // convective fluxes
+ const Scalar insideTerm = insideVolVars.density(couplingPhaseIdx(domainI)) * insideVolVars.enthalpy(couplingPhaseIdx(domainI));
+ const Scalar outsideTerm = outsideVolVars.density(couplingPhaseIdx(domainJ)) * outsideVolVars.enthalpy(couplingPhaseIdx(domainJ));
+
+ flux += ParentType::advectiveFlux(insideTerm, outsideTerm, velocity, insideIsUpstream);
+
+ // conductive energy fluxes
+ flux += ParentType::conductiveEnergyFlux(domainI, domainJ, scvf, scvI, scvJ, insideVolVars, outsideVolVars);
+
+ return flux;
+ }
+};
+
+/*!
+ * \ingroup FreeFlowPoreNetworkCoupling
+ * \brief Coupling conditions specialization for compositional models.
+ */
+template<class MDTraits, class CouplingManager, bool enableEnergyBalance>
+class FreeFlowPoreNetworkCouplingConditionsImplementation<MDTraits, CouplingManager, enableEnergyBalance, true>
+: public FreeFlowPoreNetworkCouplingConditionsImplementationBase<MDTraits, CouplingManager>
+{
+ using ParentType = FreeFlowPoreNetworkCouplingConditionsImplementationBase<MDTraits, CouplingManager>;
+ using Scalar = typename MDTraits::Scalar;
+
+ // the sub domain type tags
+ template<std::size_t id>
+ using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
+
+ template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
+ template<std::size_t id> using Element = typename GridGeometry<id>::GridView::template Codim<0>::Entity;
+ template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView;
+ template<std::size_t id> using SubControlVolumeFace = typename GridGeometry<id>::LocalView::SubControlVolumeFace;
+ template<std::size_t id> using SubControlVolume = typename GridGeometry<id>::LocalView::SubControlVolume;
+ template<std::size_t id> using Indices = typename GetPropType<SubDomainTypeTag<id>, Properties::ModelTraits>::Indices;
+ template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView;
+ template<std::size_t id> using VolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::VolumeVariables;
+ template<std::size_t id> using FluidSystem = typename VolumeVariables<id>::FluidSystem;
+ template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>;
+ template<std::size_t id> using NumEqVector = typename Problem<id>::Traits::NumEqVector;
+
+ static constexpr bool useMoles = GetPropType<SubDomainTypeTag<ParentType::freeFlowMassIndex>, Properties::ModelTraits>::useMoles();
+ static constexpr auto numComponents = GetPropType<SubDomainTypeTag<ParentType::freeFlowMassIndex>, Properties::ModelTraits>::numFluidComponents();
+ static constexpr auto referenceSystemFormulation = GetPropType<SubDomainTypeTag<ParentType::freeFlowMassIndex>, Properties::MolecularDiffusionType>::referenceSystemFormulation();
+ static constexpr auto replaceCompEqIdx = GetPropType<SubDomainTypeTag<ParentType::freeFlowMassIndex>, Properties::ModelTraits>::replaceCompEqIdx();
+
+ static_assert(GetPropType<SubDomainTypeTag<ParentType::poreNetworkIndex>, Properties::ModelTraits>::numFluidComponents() == numComponents,
+ "Models must have same number of components");
+
+public:
+ using ParentType::ParentType;
+ using ParentType::couplingPhaseIdx;
+ using ParentType::couplingCompIdx;
+
+ using NumCompVector = Dune::FieldVector<Scalar, numComponents>;
+
+ /*!
+ * \brief Returns the mass flux across the coupling boundary as seen from the pore-network domain.
+ */
+ template<class CouplingContext>
+ static NumCompVector massCouplingCondition(Dune::index_constant<ParentType::poreNetworkIndex> domainI,
+ Dune::index_constant<ParentType::freeFlowMassIndex> domainJ,
+ const FVElementGeometry<ParentType::poreNetworkIndex>& fvGeometry,
+ const SubControlVolume<ParentType::poreNetworkIndex>& scv,
+ const ElementVolumeVariables<ParentType::poreNetworkIndex>& insideVolVars,
+ const CouplingContext& context)
+ {
+ NumCompVector massFlux(0.0);
+ const auto& pnmVolVars = insideVolVars[scv.indexInElement()];
+
+ for (const auto& c : context)
+ {
+ // positive values indicate flux into pore-network region
+ const Scalar normalFFVelocity = c.velocity * c.scvf.unitOuterNormal();
+ const bool pnmIsUpstream = std::signbit(normalFFVelocity);
+ const Scalar normalPNMVelocity = -normalFFVelocity;
+ const Scalar area = c.scvf.area() * c.volVars.extrusionFactor();
+
+ auto flux = massFlux_(domainI, domainJ, c.scvf, scv, c.scv, pnmVolVars, c.volVars, normalPNMVelocity, pnmIsUpstream);
+
+ // flux is used as source term (volumetric flux): positive values mean influx
+ // thus, it is multiplied with area and we flip the sign
+ flux *= area;
+ flux *= -1.0;
+
+ massFlux += flux;
+ }
+
+ return massFlux;
+ }
+
+ /*!
+ * \brief Returns the mass flux across the coupling boundary as seen from the free-flow domain.
+ */
+ template<class CouplingContext>
+ static NumCompVector massCouplingCondition(Dune::index_constant<ParentType::freeFlowMassIndex> domainI,
+ Dune::index_constant<ParentType::poreNetworkIndex> domainJ,
+ const FVElementGeometry<ParentType::freeFlowMassIndex>& fvGeometry,
+ const SubControlVolumeFace<ParentType::freeFlowMassIndex>& scvf,
+ const ElementVolumeVariables<ParentType::freeFlowMassIndex>& insideVolVars,
+ const CouplingContext& context)
+ {
+ // positive values indicate flux into pore-network region
+ const Scalar normalFFVelocity = context.velocity * scvf.unitOuterNormal();
+ const bool ffIsUpstream = !std::signbit(normalFFVelocity);
+ const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+ const auto& ffVolVars = insideVolVars[scvf.insideScvIdx()];
+
+ return massFlux_(domainI, domainJ, scvf, insideScv, context.scv, ffVolVars, context.volVars, normalFFVelocity, ffIsUpstream);
+ }
+
+ /*!
+ * \brief Returns the energy flux across the coupling boundary as seen from the pore network.
+ */
+ template<class CouplingContext, bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0>
+ static Scalar energyCouplingCondition(Dune::index_constant<ParentType::poreNetworkIndex> domainI,
+ Dune::index_constant<ParentType::freeFlowMassIndex> domainJ,
+ const FVElementGeometry<ParentType::poreNetworkIndex>& fvGeometry,
+ const SubControlVolume<ParentType::poreNetworkIndex>& scv,
+ const ElementVolumeVariables<ParentType::poreNetworkIndex>& insideVolVars,
+ const CouplingContext& context)
+ {
+ Scalar energyFlux(0.0);
+
+ //use VolumeVariables (same type as for context.volVars) instead of ElementVolumeVariables
+ const auto& pnmVolVars = insideVolVars[scv.indexInElement()];
+
+ for(const auto& c : context)
+ {
+ // positive values indicate flux into pore-network region
+ const Scalar normalFFVelocity = c.velocity * c.scvf.unitOuterNormal();
+ const bool pnmIsUpstream = std::signbit(normalFFVelocity);
+ const Scalar normalPNMVelocity = -normalFFVelocity;
+ const Scalar area = c.scvf.area() * c.volVars.extrusionFactor();
+
+ auto flux = energyFlux_(domainI, domainJ, c.scvf, scv, c.scv, pnmVolVars, c.volVars, normalPNMVelocity, pnmIsUpstream);
+
+ // flux is used as source term (volumetric flux): positive values mean influx
+ // thus, it is multiplied with area and we flip the sign
+ flux *= area;
+ flux *= -1.0;
+ energyFlux += flux;
+ }
+
+ return energyFlux;
+ }
+
+ /*!
+ * \brief Returns the energy flux across the coupling boundary as seen from the free-flow domain.
+ */
+ template<class CouplingContext, bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0>
+ static Scalar energyCouplingCondition(Dune::index_constant<ParentType::freeFlowMassIndex> domainI,
+ Dune::index_constant<ParentType::poreNetworkIndex> domainJ,
+ const FVElementGeometry<ParentType::freeFlowMassIndex>& fvGeometry,
+ const SubControlVolumeFace<ParentType::freeFlowMassIndex>& scvf,
+ const ElementVolumeVariables<ParentType::freeFlowMassIndex>& insideVolVars,
+ const CouplingContext& context)
+ {
+ // positive values indicate flux into pore-network region
+ const Scalar normalFFVelocity = context.velocity * scvf.unitOuterNormal();
+ const bool ffIsUpstream = !std::signbit(normalFFVelocity);
+ const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+ //use VolumeVariables (same type as for context.volVars) instead of ElementVolumeVariables
+ const auto& ffVolVars = insideVolVars[scvf.insideScvIdx()];
+
+ return energyFlux_(domainI, domainJ, scvf, insideScv, context.scv, ffVolVars, context.volVars, normalFFVelocity, ffIsUpstream);
}
private:
+ /*!
+ * \brief Evaluate the compositional mole/mass flux across the interface.
+ */
+ template<std::size_t i, std::size_t j>
+ static NumCompVector massFlux_(Dune::index_constant<i> domainI,
+ Dune::index_constant<j> domainJ,
+ const SubControlVolumeFace<ParentType::freeFlowMassIndex>& scvf,
+ const SubControlVolume<i>& scvI,
+ const SubControlVolume<j>& scvJ,
+ const VolumeVariables<i>& insideVolVars,
+ const VolumeVariables<j>& outsideVolVars,
+ const Scalar velocity,
+ const bool insideIsUpstream)
+ {
+ NumCompVector flux(0.0);
+
+ auto moleOrMassFraction = [&](const auto& volVars, int phaseIdx, int compIdx)
+ { return useMoles ? volVars.moleFraction(phaseIdx, compIdx) : volVars.massFraction(phaseIdx, compIdx); };
+
+ auto moleOrMassDensity = [&](const auto& volVars, int phaseIdx)
+ { return useMoles ? volVars.molarDensity(phaseIdx) : volVars.density(phaseIdx); };
+
+ // advective fluxes
+ auto insideTerm = [&](int compIdx)
+ { return moleOrMassFraction(insideVolVars, couplingPhaseIdx(domainI), compIdx) * moleOrMassDensity(insideVolVars, couplingPhaseIdx(domainI)); };
+ auto outsideTerm = [&](int compIdx)
+ { return moleOrMassFraction(outsideVolVars, couplingPhaseIdx(domainJ), compIdx) * moleOrMassDensity(outsideVolVars, couplingPhaseIdx(domainJ)); };
+
+
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ {
+ const int domainICompIdx = couplingCompIdx(domainI, compIdx);
+ const int domainJCompIdx = couplingCompIdx(domainJ, compIdx);
+
+ assert(FluidSystem<i>::componentName(domainICompIdx) == FluidSystem<j>::componentName(domainJCompIdx));
+
+ flux[domainICompIdx] += ParentType::advectiveFlux(insideTerm(domainICompIdx), outsideTerm(domainJCompIdx), velocity, insideIsUpstream);
+ }
+
+ // diffusive fluxes
+ NumCompVector diffusiveFlux = diffusiveMolecularFlux_(domainI, domainJ, scvf, scvI, scvJ, insideVolVars, outsideVolVars);
+
+ //convert to correct units if necessary
+ if (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged && useMoles)
+ {
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ {
+ const int domainICompIdx = couplingCompIdx(domainI, compIdx);
+ diffusiveFlux[domainICompIdx] /= FluidSystem<i>::molarMass(domainICompIdx);
+ }
+ }
+ if (referenceSystemFormulation == ReferenceSystemFormulation::molarAveraged && !useMoles)
+ {
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ {
+ const int domainICompIdx = couplingCompIdx(domainI, compIdx);
+ diffusiveFlux[domainICompIdx] *= FluidSystem<i>::molarMass(domainICompIdx);
+ }
+ }
+
+ // total flux
+ flux += diffusiveFlux;
+
+ // convert to total mass/mole balance, if set be user
+ if (replaceCompEqIdx < numComponents)
+ flux[replaceCompEqIdx] = std::accumulate(flux.begin(), flux.end(), 0.0);
+
+ return flux;
+ }
+
+ /*!
+ * \brief Evaluate the energy flux (convective and conductive) across the interface.
+ */
+ template<std::size_t i, std::size_t j, bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0>
+ static Scalar energyFlux_(Dune::index_constant<i> domainI,
+ Dune::index_constant<j> domainJ,
+ const SubControlVolumeFace<ParentType::freeFlowMassIndex>& scvf,
+ const SubControlVolume<i>& scvI,
+ const SubControlVolume<j>& scvJ,
+ const VolumeVariables<i>& insideVolVars,
+ const VolumeVariables<j>& outsideVolVars,
+ const Scalar velocity,
+ const bool insideIsUpstream)
+ {
+ Scalar flux(0.0);
+
+ // convective fluxes
+ const Scalar insideTerm = insideVolVars.density(couplingPhaseIdx(domainI)) * insideVolVars.enthalpy(couplingPhaseIdx(domainI));
+ const Scalar outsideTerm = outsideVolVars.density(couplingPhaseIdx(domainJ)) * outsideVolVars.enthalpy(couplingPhaseIdx(domainJ));
+
+ flux += ParentType::advectiveFlux(insideTerm, outsideTerm, velocity, insideIsUpstream);
+
+ // conductive fluxes
+ flux += ParentType::conductiveEnergyFlux(domainI, domainJ, scvf, scvI, scvJ, insideVolVars, outsideVolVars);
+
+ // TODO: add contribution from diffusive fluxes
+
+ return flux;
+ }
+
+ /*!
+ * \brief Evaluate the diffusive mole/mass flux across the interface.
+ */
+ template<std::size_t i, std::size_t j>
+ static NumCompVector diffusiveMolecularFlux_(Dune::index_constant<i> domainI,
+ Dune::index_constant<j> domainJ,
+ const SubControlVolumeFace<ParentType::freeFlowMassIndex>& scvf,
+ const SubControlVolume<i>& scvI,
+ const SubControlVolume<j>& scvJ,
+ const VolumeVariables<i>& volVarsI,
+ const VolumeVariables<j>& volVarsJ)
+ {
+ NumCompVector diffusiveFlux(0.0);
+ const Scalar avgDensity = 0.5*(massOrMolarDensity(volVarsI, referenceSystemFormulation, couplingPhaseIdx(domainI))
+ + massOrMolarDensity(volVarsJ, referenceSystemFormulation, couplingPhaseIdx(domainJ)));
+
+ const auto& freeFlowVolVars = std::get<const VolumeVariables<ParentType::freeFlowMassIndex>&>(std::forward_as_tuple(volVarsI, volVarsJ));
+ const auto& ffScv = std::get<const SubControlVolume<ParentType::freeFlowMassIndex>&>(std::forward_as_tuple(scvI, scvJ));
+ const Scalar distance = ParentType::getDistance_(ffScv, scvf);
+
+ for (int compIdx = 1; compIdx < numComponents; ++compIdx)
+ {
+ const int freeFlowMainCompIdx = couplingPhaseIdx(ParentType::freeFlowMassIndex);
+ const int domainICompIdx = couplingCompIdx(domainI, compIdx);
+ const int domainJCompIdx = couplingCompIdx(domainJ, compIdx);
+
+ assert(FluidSystem<i>::componentName(domainICompIdx) == FluidSystem<j>::componentName(domainJCompIdx));
+
+ const Scalar massOrMoleFractionI = massOrMoleFraction(volVarsI, referenceSystemFormulation, couplingPhaseIdx(domainI), domainICompIdx);
+ const Scalar massOrMoleFractionJ = massOrMoleFraction(volVarsJ, referenceSystemFormulation, couplingPhaseIdx(domainJ), domainJCompIdx);
+ const Scalar deltaMassOrMoleFrac = massOrMoleFractionJ - massOrMoleFractionI;
+
+ const Scalar tij = freeFlowVolVars.effectiveDiffusionCoefficient(couplingPhaseIdx(ParentType::freeFlowMassIndex),
+ freeFlowMainCompIdx,
+ couplingCompIdx(ParentType::freeFlowMassIndex, compIdx))
+ / distance;
+ diffusiveFlux[domainICompIdx] += -avgDensity * tij * deltaMassOrMoleFrac;
+ }
+
+ const Scalar cumulativeFlux = std::accumulate(diffusiveFlux.begin(), diffusiveFlux.end(), 0.0);
+ diffusiveFlux[couplingCompIdx(domainI, 0)] = -cumulativeFlux;
+
+ return diffusiveFlux;
+ }
+
+ static Scalar getComponentEnthalpy_(const VolumeVariables<ParentType::freeFlowMassIndex>& volVars, int phaseIdx, int compIdx)
+ {
+ return FluidSystem<ParentType::freeFlowMassIndex>::componentEnthalpy(volVars.fluidState(), 0, compIdx);
+ }
+
+ static Scalar getComponentEnthalpy_(const VolumeVariables<ParentType::poreNetworkIndex>& volVars, int phaseIdx, int compIdx)
+ {
+ return FluidSystem<ParentType::poreNetworkIndex>::componentEnthalpy(volVars.fluidState(), phaseIdx, compIdx);
+ }
};
} // end namespace Dumux
diff --git a/dumux/multidomain/boundary/freeflowporenetwork/couplingmanager.hh b/dumux/multidomain/boundary/freeflowporenetwork/couplingmanager.hh
index 0c42137141..b032ea77cc 100644
--- a/dumux/multidomain/boundary/freeflowporenetwork/couplingmanager.hh
+++ b/dumux/multidomain/boundary/freeflowporenetwork/couplingmanager.hh
@@ -235,12 +235,12 @@ public:
/*!
* \brief Returns the mass flux across the coupling boundary.
*/
- auto massCouplingCondition(Dune::index_constant<poreNetworkIndex> domainI, Dune::index_constant<freeFlowMassIndex> domainJ,
+ auto massCouplingCondition(Dune::index_constant<poreNetworkIndex> domainI,
+ Dune::index_constant<freeFlowMassIndex> domainJ,
const FVElementGeometry<poreNetworkIndex>& fvGeometry,
const typename FVElementGeometry<poreNetworkIndex>::SubControlVolume& scv,
const ElementVolumeVariables<poreNetworkIndex>& elemVolVars) const
{
-
const auto& couplingContext = this->subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj) -> const auto& {
return cm.couplingContext(ii, fvGeometry, scv);
});
@@ -261,12 +261,12 @@ public:
/*!
* \brief Returns the mass flux across the coupling boundary.
*/
- auto massCouplingCondition(Dune::index_constant<freeFlowMassIndex> domainI, Dune::index_constant<poreNetworkIndex> domainJ,
+ auto massCouplingCondition(Dune::index_constant<freeFlowMassIndex> domainI,
+ Dune::index_constant<poreNetworkIndex> domainJ,
const FVElementGeometry<freeFlowMassIndex>& fvGeometry,
const typename FVElementGeometry<freeFlowMassIndex>::SubControlVolumeFace& scvf,
const ElementVolumeVariables<freeFlowMassIndex>& elemVolVars) const
{
-
const auto& couplingContext = this->subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj) -> const auto& {
return cm.couplingContext(ii, fvGeometry, scvf);
});
@@ -275,6 +275,49 @@ public:
return CouplingConditions::massCouplingCondition(domainI, domainJ, fvGeometry, scvf, elemVolVars, couplingContext);
}
+ /*!
+ * \brief Returns the energy flux across the coupling boundary.
+ */
+ auto energyCouplingCondition(Dune::index_constant<poreNetworkIndex> domainI,
+ Dune::index_constant<freeFlowMassIndex> domainJ,
+ const FVElementGeometry<poreNetworkIndex>& fvGeometry,
+ const typename FVElementGeometry<poreNetworkIndex>::SubControlVolume& scv,
+ const ElementVolumeVariables<poreNetworkIndex>& elemVolVars) const
+ {
+ const auto& couplingContext = this->subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj) -> const auto& {
+ return cm.couplingContext(ii, fvGeometry, scv);
+ });
+
+ const auto& freeFlowMassGridGeometry = this->subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj) -> const auto& {
+ return cm.problem(jj).gridGeometry();
+ });
+
+ for (auto& c : couplingContext)
+ {
+ const auto& freeFlowElement = freeFlowMassGridGeometry.element(c.scv.elementIndex());
+ c.velocity = this->subCouplingManager(freeFlowMomentumIndex, freeFlowMassIndex).faceVelocity(freeFlowElement, c.scvf);
+ }
+
+ return CouplingConditions::energyCouplingCondition(domainI, domainJ, fvGeometry, scv, elemVolVars, couplingContext);
+ }
+
+ /*!
+ * \brief Returns the energy flux across the coupling boundary.
+ */
+ auto energyCouplingCondition(Dune::index_constant<freeFlowMassIndex> domainI,
+ Dune::index_constant<poreNetworkIndex> domainJ,
+ const FVElementGeometry<freeFlowMassIndex>& fvGeometry,
+ const typename FVElementGeometry<freeFlowMassIndex>::SubControlVolumeFace& scvf,
+ const ElementVolumeVariables<freeFlowMassIndex>& elemVolVars) const
+ {
+ const auto& couplingContext = this->subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj) -> const auto& {
+ return cm.couplingContext(ii, fvGeometry, scvf);
+ });
+
+ couplingContext.velocity = this->subCouplingManager(freeFlowMomentumIndex, freeFlowMassIndex).faceVelocity(fvGeometry.element(), scvf);
+ return CouplingConditions::energyCouplingCondition(domainI, domainJ, fvGeometry, scvf, elemVolVars, couplingContext);
+ }
+
//////////////////////// Conditions for FreeFlowMomentum - PoreNetwork coupling //////////
///////////////////////////////////////////////////////////////////////////////////////////
@@ -444,19 +487,6 @@ public:
});
}
- // /*!
- // * \brief Returns whether a given scvf is coupled to the other domain
- // */
- // bool isCoupledLateralScvf(Dune::index_constant<freeFlowMomentumIndex> domainI,
- // Dune::index_constant<poreNetworkIndex> domainJ,
- // const SubControlVolumeFace<freeFlowMomentumIndex>& scvf) const
- // {
- // return this->subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj){
- // return cm.isCoupledLateralScvf(ii, scvf);
- // });
- // }
-
-
using ParentType::couplingStencil;
/*!
* \brief returns an iterable container of all indices of degrees of freedom of domain j
--
GitLab