25#ifndef DUMUX_FLUX_PNM_ADVECTION_HH
26#define DUMUX_FLUX_PNM_ADVECTION_HH
33template<
class... TransmissibilityLawTypes>
45template<
class ScalarT,
class... TransmissibilityLawTypes>
63 template<
class Problem,
class Element,
class FVElementGeometry,
64 class ElementVolumeVariables,
class SubControlVolumeFace,
class ElemFluxVarsCache>
66 const Element& element,
67 const FVElementGeometry& fvGeometry,
68 const ElementVolumeVariables& elemVolVars,
69 const SubControlVolumeFace& scvf,
71 const ElemFluxVarsCache& elemFluxVarsCache)
73 const auto& fluxVarsCache = elemFluxVarsCache[scvf];
76 const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
77 const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
78 const auto& insideVolVars = elemVolVars[insideScv];
79 const auto& outsideVolVars = elemVolVars[outsideScv];
82 const Scalar deltaP = insideVolVars.pressure(phaseIdx) - outsideVolVars.pressure(phaseIdx);
83 const Scalar transmissibility = fluxVarsCache.transmissibility(phaseIdx);
85 Scalar volumeFlow = transmissibility*deltaP;
88 static const bool enableGravity = getParamFromGroup<bool>(problem.paramGroup(),
"Problem.EnableGravity");
91 const Scalar rho = 0.5*insideVolVars.density(phaseIdx) + 0.5*outsideVolVars.density(phaseIdx);
92 const Scalar g = problem.spatialParams().gravity(scvf.center()) * scvf.unitOuterNormal();
97 volumeFlow += transmissibility * fluxVarsCache.poreToPoreDistance() * rho * g;
106 template<
class Problem,
class Element,
class FVElementGeometry,
class ElementVolumeVariables,
class FluxVariablesCache>
108 const Element& element,
109 const FVElementGeometry& fvGeometry,
110 const typename FVElementGeometry::SubControlVolumeFace& scvf,
111 const ElementVolumeVariables& elemVolVars,
112 const FluxVariablesCache& fluxVarsCache,
115 if constexpr (ElementVolumeVariables::VolumeVariables::numFluidPhases() == 1)
116 return Transmissibility::singlePhaseTransmissibility(problem, element, fvGeometry, scvf, elemVolVars, fluxVarsCache, phaseIdx);
119 static_assert(ElementVolumeVariables::VolumeVariables::numFluidPhases() == 2);
121 const auto& spatialParams = problem.spatialParams();
122 using FluidSystem =
typename ElementVolumeVariables::VolumeVariables::FluidSystem;
123 const int wPhaseIdx = spatialParams.template wettingPhase<FluidSystem>(element, elemVolVars);
124 const bool invaded = fluxVarsCache.invaded();
126 if (phaseIdx == wPhaseIdx)
128 return invaded ? Transmissibility::wettingLayerTransmissibility(element, fvGeometry, scvf, fluxVarsCache)
129 : Transmissibility::singlePhaseTransmissibility(problem, element, fvGeometry, scvf, elemVolVars, fluxVarsCache, phaseIdx);
133 return invaded ? Transmissibility::nonWettingPhaseTransmissibility(element, fvGeometry, scvf, fluxVarsCache)
139 template<
class Problem,
class Element,
class FVElementGeometry,
class ElementVolumeVariables,
class FluxVariablesCache>
141 const Element& element,
142 const FVElementGeometry& fvGeometry,
143 const ElementVolumeVariables& elemVolVars,
144 const typename FVElementGeometry::SubControlVolumeFace& scvf,
145 const FluxVariablesCache& fluxVarsCache)
147 static_assert(ElementVolumeVariables::VolumeVariables::numFluidPhases() == 1);
160template <
class ScalarT,
class... TransmissibilityLawTypes>
176 using ParentType =
typename TransmissibilityCreepingFlow::SinglePhaseCache;
180 template<
class Problem,
class Element,
class FVElementGeometry,
class ElementVolumeVariables,
class FluxVariablesCache>
181 void fill(
const Problem& problem,
182 const Element& element,
183 const FVElementGeometry& fvGeometry,
184 const typename FVElementGeometry::SubControlVolumeFace& scvf,
185 const ElementVolumeVariables& elemVolVars,
186 const FluxVariablesCache& fluxVarsCache,
189 ParentType::fill(problem, element, fvGeometry,scvf, elemVolVars, fluxVarsCache, phaseIdx);
190 const auto elemSol =
elementSolution(element, elemVolVars, fvGeometry);
192 for (
const auto& scv : scvs(fvGeometry))
194 const auto localIdx = scv.indexInElement();
195 const Scalar throatToPoreAreaRatio = fluxVarsCache.throatCrossSectionalArea() / problem.spatialParams().poreCrossSectionalArea(element, scv, elemSol);
198 const Scalar momentumCoefficient = problem.spatialParams().momentumCoefficient(element, scv, elemSol);
201 const Scalar kineticEnergyCoefficient = problem.spatialParams().kineticEnergyCoefficient(element, scv, elemSol);
203 expansionCoefficient_[localIdx] = getExpansionCoefficient_(throatToPoreAreaRatio, momentumCoefficient, kineticEnergyCoefficient);
204 contractionCoefficient_[localIdx] = getContractionCoefficient_(throatToPoreAreaRatio, momentumCoefficient, kineticEnergyCoefficient);
209 {
return expansionCoefficient_[downstreamIdx]; }
212 {
return contractionCoefficient_[upstreamIdx]; }
216 Scalar getExpansionCoefficient_(
const Scalar throatToPoreAreaRatio,
const Scalar momentumCoefficient,
const Scalar kineticEnergyCoefficient)
const
218 Scalar expansionCoefficient = throatToPoreAreaRatio * throatToPoreAreaRatio * (2 * momentumCoefficient - kineticEnergyCoefficient)
219 + kineticEnergyCoefficient - 2 * momentumCoefficient * throatToPoreAreaRatio
220 - (1 - throatToPoreAreaRatio * throatToPoreAreaRatio);
225 Scalar getContractionCoefficient_(
const Scalar throatToPoreAreaRatio,
const Scalar momentumCoefficient,
const Scalar kineticEnergyCoefficient)
const
227 const Scalar contractionAreaRatio = getContractionAreaRatio_(throatToPoreAreaRatio);
228 Scalar contractionCoefficient = (1 - (throatToPoreAreaRatio * throatToPoreAreaRatio * kineticEnergyCoefficient - 2 * momentumCoefficient )
229 * contractionAreaRatio * contractionAreaRatio - 2 * contractionAreaRatio) / (contractionAreaRatio * contractionAreaRatio);
234 Scalar getContractionAreaRatio_(
const Scalar throatToPoreAreaRatio)
const
236 return 1.0-(1.0-throatToPoreAreaRatio)/(2.08*(1.0-throatToPoreAreaRatio)+0.5371);
239 std::array<Scalar, 2> expansionCoefficient_;
240 std::array<Scalar, 2> contractionCoefficient_;
251 template<
class Problem,
class Element,
class FVElementGeometry,
252 class ElementVolumeVariables,
class SubControlVolumeFace,
class ElemFluxVarsCache>
254 const Element& element,
255 const FVElementGeometry& fvGeometry,
256 const ElementVolumeVariables& elemVolVars,
257 const SubControlVolumeFace& scvf,
259 const ElemFluxVarsCache& elemFluxVarsCache)
261 const auto& fluxVarsCache = elemFluxVarsCache[scvf];
264 const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
265 const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
266 const auto& insideVolVars = elemVolVars[insideScv];
267 const auto& outsideVolVars = elemVolVars[outsideScv];
270 Scalar deltaP = insideVolVars.pressure(phaseIdx) - outsideVolVars.pressure(phaseIdx);
273 static const bool enableGravity = getParamFromGroup<bool>(problem.paramGroup(),
"Problem.EnableGravity");
276 const Scalar rho = 0.5*insideVolVars.density(phaseIdx) + 0.5*outsideVolVars.density(phaseIdx);
279 const auto poreToPoreVector = fluxVarsCache.poreToPoreDistance() * scvf.unitOuterNormal();
280 const auto& gravityVector = problem.spatialParams().gravity(scvf.center());
282 deltaP += rho * (gravityVector * poreToPoreVector);
284 const Scalar creepingFlowTransmissibility = fluxVarsCache.transmissibility(phaseIdx);
285 const Scalar throatCrossSectionalArea = fluxVarsCache.throatCrossSectionalArea();
287 assert(scvf.insideScvIdx() == 0);
288 assert(scvf.outsideScvIdx() == 1);
291 const auto [upstreamIdx, downstreamIdx] = deltaP > 0 ? std::pair(0, 1) : std::pair(1, 0);
293 const Scalar contractionCoefficient = fluxVarsCache.singlePhaseFlowVariables().contractionCoefficient(upstreamIdx);
294 const Scalar expansionCoefficient = fluxVarsCache.singlePhaseFlowVariables().expansionCoefficient(downstreamIdx);
295 const Scalar mu = elemVolVars[upstreamIdx].viscosity();
306 const Scalar A = (contractionCoefficient * elemVolVars[upstreamIdx].density() + expansionCoefficient * elemVolVars[downstreamIdx].density())
307 / (2.0 * mu * mu * throatCrossSectionalArea * throatCrossSectionalArea);
308 const Scalar B = 1.0/creepingFlowTransmissibility;
310 const Scalar C = -abs(deltaP);
313 const auto discriminant = B*B - 4*A*C;
314 const auto q = (-B + sqrt(discriminant)) / (2*A);
317 if (upstreamIdx == 0)
324 template<
class Problem,
class Element,
class FVElementGeometry,
class ElementVolumeVariables,
class FluxVariablesCache>
326 const Element& element,
327 const FVElementGeometry& fvGeometry,
328 const typename FVElementGeometry::SubControlVolumeFace& scvf,
329 const ElementVolumeVariables& elemVolVars,
330 const FluxVariablesCache& fluxVarsCache,
333 static_assert(ElementVolumeVariables::VolumeVariables::numFluidPhases() == 1);
334 return Transmissibility::singlePhaseTransmissibility(problem, element, fvGeometry, scvf, elemVolVars, fluxVarsCache, phaseIdx);
337 template<
class Problem,
class Element,
class FVElementGeometry,
class ElementVolumeVariables,
class FluxVariablesCache>
339 const Element& element,
340 const FVElementGeometry& fvGeometry,
341 const ElementVolumeVariables& elemVolVars,
342 const typename FVElementGeometry::SubControlVolumeFace& scvf,
343 const FluxVariablesCache& fluxVarsCache)
345 static_assert(ElementVolumeVariables::VolumeVariables::numFluidPhases() == 1);
The infrastructure to retrieve run-time parameters from Dune::ParameterTrees.
auto elementSolution(const Element &element, const SolutionVector &sol, const GridGeometry &gg) -> std::enable_if_t< GridGeometry::discMethod==DiscretizationMethod::box, BoxElementSolution< typename GridGeometry::LocalView, std::decay_t< decltype(std::declval< SolutionVector >()[0])> > >
Make an element solution for box schemes.
Definition: box/elementsolution.hh:118
Definition: discretization/porenetwork/fvelementgeometry.hh:33
Definition: advection.hh:31
Definition: advection.hh:34
Definition: advection.hh:47
static std::array< Scalar, 2 > calculateTransmissibilities(const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const typename FVElementGeometry::SubControlVolumeFace &scvf, const FluxVariablesCache &fluxVarsCache)
Definition: advection.hh:140
ScalarT Scalar
Export the Scalar type.
Definition: advection.hh:51
static Scalar calculateTransmissibility(const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const typename FVElementGeometry::SubControlVolumeFace &scvf, const ElementVolumeVariables &elemVolVars, const FluxVariablesCache &fluxVarsCache, const int phaseIdx)
Returns the throat conductivity.
Definition: advection.hh:107
static Scalar flux(const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const SubControlVolumeFace &scvf, const int phaseIdx, const ElemFluxVarsCache &elemFluxVarsCache)
Returns the advective flux of a fluid phase across the given sub-control volume face (corresponding t...
Definition: advection.hh:65
Definition: advection.hh:162
static Scalar flux(const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const SubControlVolumeFace &scvf, const int phaseIdx, const ElemFluxVarsCache &elemFluxVarsCache)
Returns the advective flux of a fluid phase across the given sub-control volume face (corresponding t...
Definition: advection.hh:253
ScalarT Scalar
Export the Scalar type.
Definition: advection.hh:165
static Scalar calculateTransmissibility(const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const typename FVElementGeometry::SubControlVolumeFace &scvf, const ElementVolumeVariables &elemVolVars, const FluxVariablesCache &fluxVarsCache, const int phaseIdx)
Definition: advection.hh:325
static std::array< Scalar, 2 > calculateTransmissibilities(const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const typename FVElementGeometry::SubControlVolumeFace &scvf, const FluxVariablesCache &fluxVarsCache)
Definition: advection.hh:338
Inherit transmissibility from creeping flow transmissibility law to cache non-creeping flow-related p...
Definition: advection.hh:172
Definition: advection.hh:175
void fill(const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const typename FVElementGeometry::SubControlVolumeFace &scvf, const ElementVolumeVariables &elemVolVars, const FluxVariablesCache &fluxVarsCache, const int phaseIdx)
Definition: advection.hh:181
ScalarT Scalar
Definition: advection.hh:178
Scalar contractionCoefficient(int upstreamIdx) const
Definition: advection.hh:211
Scalar expansionCoefficient(int downstreamIdx) const
Definition: advection.hh:208