3.5-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
freeflowporenetwork/couplingconditions.hh
Go to the documentation of this file.
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 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
25#ifndef DUMUX_MD_FREEFLOW_PORENETWORK_COUPLINGCONDITIONS_HH
26#define DUMUX_MD_FREEFLOW_PORENETWORK_COUPLINGCONDITIONS_HH
27
28#include <numeric>
29
30#include <dune/common/exceptions.hh>
31
33#include <dumux/common/math.hh>
39
40namespace Dumux {
41
42template<class MDTraits, class CouplingManager, bool enableEnergyBalance, bool isCompositional>
44
50template<class MDTraits, class CouplingManager>
53 MDTraits, CouplingManager,
56 >;
57
62template<class MDTraits, class CouplingManager>
64{
65 using Scalar = typename MDTraits::Scalar;
66
67 template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
68 template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
69 template<std::size_t id> using Element = typename GridGeometry<id>::GridView::template Codim<0>::Entity;
70 template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView;
71 template<std::size_t id> using SubControlVolumeFace = typename GridGeometry<id>::LocalView::SubControlVolumeFace;
72 template<std::size_t id> using SubControlVolume = typename GridGeometry<id>::LocalView::SubControlVolume;
73 template<std::size_t id> using Indices = typename GetPropType<SubDomainTypeTag<id>, Properties::ModelTraits>::Indices;
74 template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView;
75 template<std::size_t id> using VolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::VolumeVariables;
76 template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>;
77 template<std::size_t id> using FluidSystem = GetPropType<SubDomainTypeTag<id>, Properties::FluidSystem>;
78 template<std::size_t id> using ModelTraits = GetPropType<SubDomainTypeTag<id>, Properties::ModelTraits>;
79 template<std::size_t id> using GlobalPosition = typename Element<id>::Geometry::GlobalCoordinate;
80 template<std::size_t id> using NumEqVector = typename Problem<id>::Traits::NumEqVector;
81
83
84public:
88private:
89
91
92 static constexpr bool adapterUsed = ModelTraits<poreNetworkIndex>::numFluidPhases() > 1;
94
95 static constexpr int enableEnergyBalance = GetPropType<SubDomainTypeTag<freeFlowMassIndex>, Properties::ModelTraits>::enableEnergyBalance();
96 static_assert(GetPropType<SubDomainTypeTag<poreNetworkIndex>, Properties::ModelTraits>::enableEnergyBalance() == enableEnergyBalance,
97 "All submodels must both be either isothermal or non-isothermal");
98
100 FluidSystem<poreNetworkIndex>>::value,
101 "All submodels must use the same fluid system");
102
103 using VelocityVector = GlobalPosition<freeFlowMomentumIndex>;
104
105public:
109 template<std::size_t i>
110 static constexpr auto couplingPhaseIdx(Dune::index_constant<i> id, int coupledPhaseIdx = 0)
111 { return IndexHelper::couplingPhaseIdx(id, coupledPhaseIdx); }
112
116 template<std::size_t i>
117 static constexpr auto couplingCompIdx(Dune::index_constant<i> id, int coupledCompIdx)
118 { return IndexHelper::couplingCompIdx(id, coupledCompIdx); }
119
120
128 template<class Context>
129 static NumEqVector<freeFlowMomentumIndex> momentumCouplingCondition(const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry,
130 const SubControlVolumeFace<freeFlowMomentumIndex>& scvf,
131 const ElementVolumeVariables<freeFlowMomentumIndex>& elemVolVars,
132 const Context& context)
133 {
134 NumEqVector<freeFlowMomentumIndex> momentumFlux(0.0);
135 const auto pnmPhaseIdx = couplingPhaseIdx(poreNetworkIndex);
136 const auto [pnmPressure, pnmViscosity] = [&]
137 {
138 for (const auto& scv : scvs(context.fvGeometry))
139 {
140 if (scv.dofIndex() == context.poreNetworkDofIdx)
141 return std::make_pair(context.elemVolVars[scv].pressure(pnmPhaseIdx), context.elemVolVars[scv].viscosity(pnmPhaseIdx));
142 }
143 DUNE_THROW(Dune::InvalidStateException, "No coupled scv found");
144 }();
145
146 // set p_freeflow = p_PNM
147 momentumFlux[scvf.normalAxis()] = pnmPressure;
148
149 // normalize pressure
150 momentumFlux[scvf.normalAxis()] -= elemVolVars.gridVolVars().problem().referencePressure(fvGeometry.element(), fvGeometry, scvf);
151
152 // Explicitly account for dv_i/dx_i, which is NOT part of the actual coupling condition. We do it here for convenience so
153 // 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
154 // as the one in the center of the element. TODO check sign
155 const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
156 const auto& frontalInternalScvf = (*scvfs(fvGeometry, scv).begin());
157 momentumFlux[scvf.normalAxis()] -= 2*VelocityGradients::velocityGradII(fvGeometry, frontalInternalScvf, elemVolVars) * pnmViscosity;
158
159 // We do NOT consider the intertia term here. If included, Newton convergence decreases drastically and the solution even does not converge to a reference solution.
160 // We furthermore assume creeping flow within the boundary layer thus neglecting this term is physically justified.
161
162 momentumFlux[scvf.normalAxis()] *= scvf.directionSign();
163
164 return momentumFlux;
165 }
166
167 template<class Context>
168 static VelocityVector interfaceThroatVelocity(const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry,
169 const SubControlVolumeFace<freeFlowMomentumIndex>& scvf,
170 const Context& context)
171 {
172 const auto& pnmElement = context.fvGeometry.element();
173 const auto& pnmFVGeometry = context.fvGeometry;
174 const auto& pnmScvf = pnmFVGeometry.scvf(0);
175 const auto& pnmElemVolVars = context.elemVolVars;
176 const auto& pnmElemFluxVarsCache = context.elemFluxVarsCache;
177 const auto& pnmProblem = pnmElemVolVars.gridVolVars().problem();
178
179 const auto pnmPhaseIdx = couplingPhaseIdx(poreNetworkIndex);
180 const Scalar area = pnmElemFluxVarsCache[pnmScvf].throatCrossSectionalArea(pnmPhaseIdx);
181
182 // only proceed if area > 0 in order to prevent division by zero (e.g., when the throat was not invaded yet)
183 if (area > 0.0)
184 {
186 PNMFluxVariables fluxVars;
187 fluxVars.init(pnmProblem, pnmElement, pnmFVGeometry, pnmElemVolVars, pnmScvf, pnmElemFluxVarsCache);
188
189 const Scalar flux = fluxVars.advectiveFlux(pnmPhaseIdx, [pnmPhaseIdx](const auto& volVars){ return volVars.mobility(pnmPhaseIdx);});
190
191 // account for the orientation of the bulk face.
192 VelocityVector velocity = (pnmElement.geometry().corner(1) - pnmElement.geometry().corner(0));
193 velocity /= velocity.two_norm();
194 velocity *= flux / area;
195
196 // TODO: Multiple throats connected to the same pore?
197 return velocity;
198 }
199 else
200 return VelocityVector(0.0);
201 }
202
206 static Scalar advectiveFlux(const Scalar insideQuantity, const Scalar outsideQuantity, const Scalar volumeFlow, bool insideIsUpstream)
207 {
208 const Scalar upwindWeight = 1.0; //TODO use Flux.UpwindWeight or something like Coupling.UpwindWeight?
209
210 if(insideIsUpstream)
211 return (upwindWeight * insideQuantity + (1.0 - upwindWeight) * outsideQuantity) * volumeFlow;
212 else
213 return (upwindWeight * outsideQuantity + (1.0 - upwindWeight) * insideQuantity) * volumeFlow;
214 }
215
216protected:
217
221 template<class Scv, class Scvf>
222 Scalar getDistance_(const Scv& scv, const Scvf& scvf) const
223 {
224 return (scv.dofPosition() - scvf.ipGlobal()).two_norm();
225 }
226};
227
232template<class MDTraits, class CouplingManager, bool enableEnergyBalance>
234: public FreeFlowPoreNetworkCouplingConditionsImplementationBase<MDTraits, CouplingManager>
235{
237 using Scalar = typename MDTraits::Scalar;
238
239 // the sub domain type tags
240 template<std::size_t id>
241 using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
242
243 template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
244 template<std::size_t id> using Element = typename GridGeometry<id>::GridView::template Codim<0>::Entity;
245 template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView;
246 template<std::size_t id> using SubControlVolumeFace = typename GridGeometry<id>::LocalView::SubControlVolumeFace;
247 template<std::size_t id> using SubControlVolume = typename GridGeometry<id>::LocalView::SubControlVolume;
248 template<std::size_t id> using Indices = typename GetPropType<SubDomainTypeTag<id>, Properties::ModelTraits>::Indices;
249 template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView;
250 template<std::size_t id> using VolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::VolumeVariables;
251
253 "Pore-network model must not be compositional");
254
255public:
256 using ParentType::ParentType;
257 using ParentType::couplingPhaseIdx;
258
262 template<class CouplingContext>
263 static Scalar massCouplingCondition(Dune::index_constant<ParentType::poreNetworkIndex> domainI,
264 Dune::index_constant<ParentType::freeFlowMassIndex> domainJ,
265 const FVElementGeometry<ParentType::poreNetworkIndex>& fvGeometry,
266 const SubControlVolume<ParentType::poreNetworkIndex>& scv,
267 const ElementVolumeVariables<ParentType::poreNetworkIndex>& insideVolVars,
268 const CouplingContext& context)
269 {
270 Scalar massFlux(0.0);
271
272 for (const auto& c : context)
273 {
274 // positive values indicate flux into pore-network region
275 const Scalar normalFFVelocity = c.velocity * c.scvf.unitOuterNormal();
276 const bool pnmIsUpstream = std::signbit(normalFFVelocity);
277
278 const Scalar pnmDensity = insideVolVars[scv].density(couplingPhaseIdx(ParentType::poreNetworkIndex));
279 const Scalar ffDensity = c.volVars.density(couplingPhaseIdx(ParentType::freeFlowMassIndex));
280 const Scalar area = c.scvf.area() * c.volVars.extrusionFactor();
281
282 // flux is used as source term: positive values mean influx
283 massFlux += ParentType::advectiveFlux(pnmDensity, ffDensity, normalFFVelocity*area, pnmIsUpstream);
284 }
285
286 return massFlux;
287 }
288
292 template<class CouplingContext>
293 static Scalar massCouplingCondition(Dune::index_constant<ParentType::freeFlowMassIndex> domainI,
294 Dune::index_constant<ParentType::poreNetworkIndex> domainJ,
295 const FVElementGeometry<ParentType::freeFlowMassIndex>& fvGeometry,
296 const SubControlVolumeFace<ParentType::freeFlowMassIndex>& scvf,
297 const ElementVolumeVariables<ParentType::freeFlowMassIndex>& insideVolVars,
298 const CouplingContext& context)
299 {
300 // positive values indicate flux into pore-network region
301 const Scalar normalFFVelocity = context.velocity * scvf.unitOuterNormal();
302 const bool ffIsUpstream = !std::signbit(normalFFVelocity);
303
304 const Scalar ffDensity = insideVolVars[scvf.insideScvIdx()].density(couplingPhaseIdx(ParentType::freeFlowMassIndex));
305 const Scalar pnmDensity = context.volVars.density(couplingPhaseIdx(ParentType::poreNetworkIndex));
306
307 // flux is used in Neumann condition: positive values mean flux out of free-flow domain
308 return ParentType::advectiveFlux(ffDensity, pnmDensity, normalFFVelocity, ffIsUpstream);
309 }
310
311private:
312
313};
314
315} // end namespace Dumux
316
317#endif
Define some often used mathematical functions.
Fick's law specilized for different discretization schemes. This file contains the data which is requ...
Index helpers for the free-flow/porous-medium-flow coupling.
Definition: adapt.hh:29
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property
Definition: propertysystem.hh:150
static constexpr auto freeFlowMassIndex
Definition: multidomain/boundary/freeflowporenetwork/couplingmanager.hh:47
static constexpr auto freeFlowMomentumIndex
Definition: multidomain/boundary/freeflowporenetwork/couplingmanager.hh:46
static constexpr auto poreNetworkIndex
Definition: multidomain/boundary/freeflowporenetwork/couplingmanager.hh:48
Traits class encapsulating model specifications.
Definition: common/properties.hh:53
Property to specify the type of a problem which has to be solved.
Definition: common/properties.hh:57
Definition: common/properties.hh:102
The type for a global container for the volume variables.
Definition: common/properties.hh:109
Container storing the different types of flux variables.
Definition: common/properties.hh:113
The type for the calculation the advective fluxes.
Definition: common/properties.hh:141
The type of the fluid system to use.
Definition: common/properties.hh:162
Helper class for calculating the velocity gradients for the Navier-Stokes model using the staggered g...
Definition: momentum/velocitygradients.hh:58
static auto velocityGradII(const FVElementGeometry &fvGeometry, const typename FVElementGeometry::SubControlVolumeFace &scvf, const ElemVolVars &elemVolVars)
Returns the in-axis velocity gradient.
Definition: momentum/velocitygradients.hh:158
Definition: freeflowporenetwork/couplingconditions.hh:43
A base class which provides some common methods used for free-flow/pore-network coupling.
Definition: freeflowporenetwork/couplingconditions.hh:64
static NumEqVector< freeFlowMomentumIndex > momentumCouplingCondition(const FVElementGeometry< freeFlowMomentumIndex > &fvGeometry, const SubControlVolumeFace< freeFlowMomentumIndex > &scvf, const ElementVolumeVariables< freeFlowMomentumIndex > &elemVolVars, const Context &context)
Returns the momentum flux across the coupling boundary.
Definition: freeflowporenetwork/couplingconditions.hh:129
static Scalar advectiveFlux(const Scalar insideQuantity, const Scalar outsideQuantity, const Scalar volumeFlow, bool insideIsUpstream)
Evaluate an advective flux across the interface and consider upwinding.
Definition: freeflowporenetwork/couplingconditions.hh:206
static VelocityVector interfaceThroatVelocity(const FVElementGeometry< freeFlowMomentumIndex > &fvGeometry, const SubControlVolumeFace< freeFlowMomentumIndex > &scvf, const Context &context)
Definition: freeflowporenetwork/couplingconditions.hh:168
static constexpr auto freeFlowMassIndex
Definition: freeflowporenetwork/couplingconditions.hh:86
static constexpr auto couplingPhaseIdx(Dune::index_constant< i > id, int coupledPhaseIdx=0)
Returns the corresponding phase index needed for coupling.
Definition: freeflowporenetwork/couplingconditions.hh:110
static constexpr auto poreNetworkIndex
Definition: freeflowporenetwork/couplingconditions.hh:87
Scalar getDistance_(const Scv &scv, const Scvf &scvf) const
Returns the distance between an scvf and the corresponding scv center.
Definition: freeflowporenetwork/couplingconditions.hh:222
static constexpr auto freeFlowMomentumIndex
Definition: freeflowporenetwork/couplingconditions.hh:85
static constexpr auto couplingCompIdx(Dune::index_constant< i > id, int coupledCompIdx)
Returns the corresponding component index needed for coupling.
Definition: freeflowporenetwork/couplingconditions.hh:117
static Scalar 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)
Returns the mass flux across the coupling boundary as seen from the pore-network domain.
Definition: freeflowporenetwork/couplingconditions.hh:263
static Scalar 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)
Returns the mass flux across the coupling boundary as seen from the free-flow domain.
Definition: freeflowporenetwork/couplingconditions.hh:293
Helper struct to choose the correct index for phases and components. This is need if the porous-mediu...
Definition: indexhelper.hh:42
This structs helps to check if the two sub models use the same fluidsystem. Specialization for the ca...
Definition: multidomain/boundary/freeflowporousmedium/traits.hh:42
The interface of the coupling manager for multi domain problems.
Definition: multidomain/couplingmanager.hh:60
Declares all properties used in Dumux.
Free functions to evaluate the transmissibilities associated with flux evaluations across sub-control...
Traits for the free-flow/porous-medium-flow coupling.