3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
couplingdata.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_STOKES_DARCY_COUPLINGDATA_HH
26#define DUMUX_STOKES_DARCY_COUPLINGDATA_HH
27
28#include <numeric>
29
31#include <dumux/common/math.hh>
36
37namespace Dumux {
38
45{
52 {
53 harmonic, arithmethic, ffOnly, pmOnly
54 };
55
60 static DiffusionCoefficientAveragingType stringToEnum(DiffusionCoefficientAveragingType, const std::string& diffusionCoefficientAveragingType)
61 {
62 if (diffusionCoefficientAveragingType == "Harmonic")
64 else if (diffusionCoefficientAveragingType == "Arithmethic")
66 else if (diffusionCoefficientAveragingType == "FreeFlowOnly")
68 else if (diffusionCoefficientAveragingType == "PorousMediumOnly")
70 else
71 DUNE_THROW(Dune::IOError, "Unknown DiffusionCoefficientAveragingType");
72 }
73
74};
75
83template<class FFFS, class PMFS>
85{
86 static_assert(FFFS::numPhases == 1, "Only single-phase fluidsystems may be used for free flow.");
87 static constexpr bool value = std::is_same<typename FFFS::MultiPhaseFluidSystem, PMFS>::value;
88};
89
95template<class FS>
96struct IsSameFluidSystem<FS, FS>
97{
98 static_assert(FS::numPhases == 1, "Only single-phase fluidsystems may be used for free flow.");
99 static constexpr bool value = std::is_same<FS, FS>::value; // always true
100};
101
102// forward declaration
103template <class TypeTag, DiscretizationMethod discMethod, ReferenceSystemFormulation referenceSystem>
105
111template<class DiffLaw>
112struct IsFicksLaw : public std::false_type {};
113
119template<class T, DiscretizationMethod discMethod, ReferenceSystemFormulation referenceSystem>
120struct IsFicksLaw<FicksLawImplementation<T, discMethod, referenceSystem>> : public std::true_type {};
121
131template<std::size_t stokesIdx, std::size_t darcyIdx, class FFFS, bool hasAdapter>
133
142template<std::size_t stokesIdx, std::size_t darcyIdx, class FFFS>
143struct IndexHelper<stokesIdx, darcyIdx, FFFS, false>
144{
148 template<std::size_t i>
149 static constexpr auto couplingPhaseIdx(Dune::index_constant<i>, int coupledPhaseIdx = 0)
150 { return coupledPhaseIdx; }
151
155 template<std::size_t i>
156 static constexpr auto couplingCompIdx(Dune::index_constant<i>, int coupledCompdIdx)
157 { return coupledCompdIdx; }
158};
159
168template<std::size_t stokesIdx, std::size_t darcyIdx, class FFFS>
169struct IndexHelper<stokesIdx, darcyIdx, FFFS, true>
170{
174 static constexpr auto couplingPhaseIdx(Dune::index_constant<stokesIdx>, int coupledPhaseIdx = 0)
175 { return 0; }
176
180 static constexpr auto couplingPhaseIdx(Dune::index_constant<darcyIdx>, int coupledPhaseIdx = 0)
181 { return FFFS::multiphaseFluidsystemPhaseIdx; }
182
186 static constexpr auto couplingCompIdx(Dune::index_constant<stokesIdx>, int coupledCompdIdx)
187 { return coupledCompdIdx; }
188
192 static constexpr auto couplingCompIdx(Dune::index_constant<darcyIdx>, int coupledCompdIdx)
193 { return FFFS::compIdx(coupledCompdIdx); }
194};
195
197template <class TypeTag, DiscretizationMethod discMethod>
198class DarcysLawImplementation;
199
201template <class TypeTag, DiscretizationMethod discMethod>
202class ForchheimersLawImplementation;
203
204
205template<class MDTraits, class CouplingManager, bool enableEnergyBalance, bool isCompositional>
207
213template<class MDTraits, class CouplingManager>
217
222template<class MDTraits, class CouplingManager>
224{
225 using Scalar = typename MDTraits::Scalar;
226
227 template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
228 template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
229 template<std::size_t id> using Element = typename GridGeometry<id>::GridView::template Codim<0>::Entity;
230 template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView;
231 template<std::size_t id> using SubControlVolumeFace = typename GridGeometry<id>::LocalView::SubControlVolumeFace;
232 template<std::size_t id> using SubControlVolume = typename GridGeometry<id>::LocalView::SubControlVolume;
233 template<std::size_t id> using Indices = typename GetPropType<SubDomainTypeTag<id>, Properties::ModelTraits>::Indices;
234 template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView;
235 template<std::size_t id> using VolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::VolumeVariables;
236 template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>;
237 template<std::size_t id> using FluidSystem = GetPropType<SubDomainTypeTag<id>, Properties::FluidSystem>;
238 template<std::size_t id> using ModelTraits = GetPropType<SubDomainTypeTag<id>, Properties::ModelTraits>;
239 template<std::size_t id> using GlobalPosition = typename Element<id>::Geometry::GlobalCoordinate;
240 static constexpr auto stokesIdx = CouplingManager::stokesIdx;
241 static constexpr auto darcyIdx = CouplingManager::darcyIdx;
242
244 using DarcysLaw = DarcysLawImplementation<SubDomainTypeTag<darcyIdx>, GridGeometry<darcyIdx>::discMethod>;
245 using ForchheimersLaw = ForchheimersLawImplementation<SubDomainTypeTag<darcyIdx>, GridGeometry<darcyIdx>::discMethod>;
246
247 static constexpr bool adapterUsed = ModelTraits<darcyIdx>::numFluidPhases() > 1;
249
250 static constexpr int enableEnergyBalance = GetPropType<SubDomainTypeTag<stokesIdx>, Properties::ModelTraits>::enableEnergyBalance();
251 static_assert(GetPropType<SubDomainTypeTag<darcyIdx>, Properties::ModelTraits>::enableEnergyBalance() == enableEnergyBalance,
252 "All submodels must both be either isothermal or non-isothermal");
253
255 FluidSystem<darcyIdx>>::value,
256 "All submodels must use the same fluid system");
257
258 using DiffusionCoefficientAveragingType = typename StokesDarcyCouplingOptions::DiffusionCoefficientAveragingType;
259
260public:
261 StokesDarcyCouplingDataImplementationBase(const CouplingManager& couplingmanager): couplingManager_(couplingmanager) {}
262
266 template<std::size_t i>
267 static constexpr auto couplingPhaseIdx(Dune::index_constant<i> id, int coupledPhaseIdx = 0)
268 { return IndexHelper::couplingPhaseIdx(id, coupledPhaseIdx); }
269
273 template<std::size_t i>
274 static constexpr auto couplingCompIdx(Dune::index_constant<i> id, int coupledCompdIdx)
275 { return IndexHelper::couplingCompIdx(id, coupledCompdIdx); }
276
281 { return couplingManager_; }
282
286 auto darcyPermeability(const Element<stokesIdx>& element, const SubControlVolumeFace<stokesIdx>& scvf) const
287 {
288 const auto& stokesContext = couplingManager().stokesCouplingContext(element, scvf);
289 return stokesContext.volVars.permeability();
290 }
291
299 template<class ElementFaceVariables>
300 Scalar momentumCouplingCondition(const Element<stokesIdx>& element,
301 const FVElementGeometry<stokesIdx>& fvGeometry,
302 const ElementVolumeVariables<stokesIdx>& stokesElemVolVars,
303 const ElementFaceVariables& stokesElemFaceVars,
304 const SubControlVolumeFace<stokesIdx>& scvf) const
305 {
306 static constexpr auto numPhasesDarcy = GetPropType<SubDomainTypeTag<darcyIdx>, Properties::ModelTraits>::numFluidPhases();
307
308 Scalar momentumFlux(0.0);
309 const auto& stokesContext = couplingManager_.stokesCouplingContext(element, scvf);
310 const auto darcyPhaseIdx = couplingPhaseIdx(darcyIdx);
311
312 // - p_pm * n_pm = p_pm * n_ff
313 const Scalar darcyPressure = stokesContext.volVars.pressure(darcyPhaseIdx);
314
315 if(numPhasesDarcy > 1)
316 momentumFlux = darcyPressure;
317 else // use pressure reconstruction for single phase models
318 momentumFlux = pressureAtInterface_(element, scvf, stokesElemFaceVars, stokesContext);
319 // TODO: generalize for permeability tensors
320
321 // normalize pressure
322 if(getPropValue<SubDomainTypeTag<stokesIdx>, Properties::NormalizePressure>())
323 momentumFlux -= couplingManager_.problem(stokesIdx).initial(scvf)[Indices<stokesIdx>::pressureIdx];
324
325 momentumFlux *= scvf.directionSign();
326
327 return momentumFlux;
328 }
329
333 Scalar advectiveFlux(const Scalar insideQuantity, const Scalar outsideQuantity, const Scalar volumeFlow, bool insideIsUpstream) const
334 {
335 const Scalar upwindWeight = 1.0; //TODO use Flux.UpwindWeight or something like Coupling.UpwindWeight?
336
337 if(insideIsUpstream)
338 return (upwindWeight * insideQuantity + (1.0 - upwindWeight) * outsideQuantity) * volumeFlow;
339 else
340 return (upwindWeight * outsideQuantity + (1.0 - upwindWeight) * insideQuantity) * volumeFlow;
341 }
342
343protected:
344
348 template<std::size_t i, std::size_t j>
349 Scalar transmissibility_(Dune::index_constant<i> domainI,
350 Dune::index_constant<j> domainJ,
351 const Scalar insideDistance,
352 const Scalar outsideDistance,
353 const Scalar avgQuantityI,
354 const Scalar avgQuantityJ,
355 const DiffusionCoefficientAveragingType diffCoeffAvgType) const
356 {
357 const Scalar totalDistance = insideDistance + outsideDistance;
358 if(diffCoeffAvgType == DiffusionCoefficientAveragingType::harmonic)
359 {
360 return harmonicMean(avgQuantityI, avgQuantityJ, insideDistance, outsideDistance)
361 / totalDistance;
362 }
363 else if(diffCoeffAvgType == DiffusionCoefficientAveragingType::arithmethic)
364 {
365 return arithmeticMean(avgQuantityI, avgQuantityJ, insideDistance, outsideDistance)
366 / totalDistance;
367 }
368 else if(diffCoeffAvgType == DiffusionCoefficientAveragingType::ffOnly)
369 return domainI == stokesIdx
370 ? avgQuantityI / totalDistance
371 : avgQuantityJ / totalDistance;
372
373 else // diffCoeffAvgType == DiffusionCoefficientAveragingType::pmOnly)
374 return domainI == darcyIdx
375 ? avgQuantityI / totalDistance
376 : avgQuantityJ / totalDistance;
377 }
378
382 template<class Scv, class Scvf>
383 Scalar getDistance_(const Scv& scv, const Scvf& scvf) const
384 {
385 return (scv.dofPosition() - scvf.ipGlobal()).two_norm();
386 }
387
391 template<std::size_t i, std::size_t j, bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0>
392 Scalar conductiveEnergyFlux_(Dune::index_constant<i> domainI,
393 Dune::index_constant<j> domainJ,
394 const FVElementGeometry<i>& fvGeometryI,
395 const FVElementGeometry<j>& fvGeometryJ,
396 const SubControlVolumeFace<i>& scvfI,
397 const SubControlVolume<i>& scvI,
398 const SubControlVolume<j>& scvJ,
399 const VolumeVariables<i>& volVarsI,
400 const VolumeVariables<j>& volVarsJ,
401 const DiffusionCoefficientAveragingType diffCoeffAvgType) const
402 {
403 const Scalar insideDistance = getDistance_(scvI, scvfI);
404 const Scalar outsideDistance = getDistance_(scvJ, scvfI);
405
406 const Scalar deltaT = volVarsJ.temperature() - volVarsI.temperature();
407 const Scalar tij = transmissibility_(domainI,
408 domainJ,
409 insideDistance,
410 outsideDistance,
411 volVarsI.effectiveThermalConductivity(),
412 volVarsJ.effectiveThermalConductivity(),
413 diffCoeffAvgType);
414
415 return -tij * deltaT;
416 }
417
421 template<class ElementFaceVariables, class CouplingContext>
422 Scalar pressureAtInterface_(const Element<stokesIdx>& element,
423 const SubControlVolumeFace<stokesIdx>& scvf,
424 const ElementFaceVariables& elemFaceVars,
425 const CouplingContext& context) const
426 {
427 GlobalPosition<stokesIdx> velocity(0.0);
428 velocity[scvf.directionIndex()] = elemFaceVars[scvf].velocitySelf();
429 const auto& darcyScvf = context.fvGeometry.scvf(context.darcyScvfIdx);
430 return computeCouplingPhasePressureAtInterface_(context.element, context.fvGeometry, darcyScvf, context.volVars, velocity, AdvectionType());
431 }
432
436 Scalar computeCouplingPhasePressureAtInterface_(const Element<darcyIdx>& element,
437 const FVElementGeometry<darcyIdx>& fvGeometry,
438 const SubControlVolumeFace<darcyIdx>& scvf,
439 const VolumeVariables<darcyIdx>& volVars,
440 const typename Element<stokesIdx>::Geometry::GlobalCoordinate& couplingPhaseVelocity,
441 ForchheimersLaw) const
442 {
443 const auto darcyPhaseIdx = couplingPhaseIdx(darcyIdx);
444 const Scalar cellCenterPressure = volVars.pressure(darcyPhaseIdx);
445 using std::sqrt;
446
447 // v + (cF*sqrt(K)*rho/mu*|v|) * v = - K/mu grad(p - rho g)
448 // multiplying with n and using a tpfa for the right-hand side yields
449 // v*n + (cF*sqrt(K)*rho/mu*|v|) * (v*n) = 1/mu * (ti*(p_center - p_interface) + rho*n^TKg)
450 // --> p_interface = (-mu*v*n + (cF*sqrt(K)*rho*|v|) * (-v*n) + rho*n^TKg)/ti + p_center
451 const auto velocity = couplingPhaseVelocity;
452 const Scalar mu = volVars.viscosity(darcyPhaseIdx);
453 const Scalar rho = volVars.density(darcyPhaseIdx);
454 const auto K = volVars.permeability();
455 const auto alpha = vtmv(scvf.unitOuterNormal(), K, couplingManager_.problem(darcyIdx).spatialParams().gravity(scvf.center()));
456
457 const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
458 const auto ti = computeTpfaTransmissibility(scvf, insideScv, K, 1.0);
459
460 // get the Forchheimer coefficient
461 Scalar cF = couplingManager_.problem(darcyIdx).spatialParams().forchCoeff(scvf);
462
463 const Scalar interfacePressure = ((-mu*(scvf.unitOuterNormal() * velocity))
464 + (-(scvf.unitOuterNormal() * velocity) * velocity.two_norm() * rho * sqrt(darcyPermeability(element, scvf)) * cF)
465 + rho * alpha)/ti + cellCenterPressure;
466 return interfacePressure;
467 }
468
472 Scalar computeCouplingPhasePressureAtInterface_(const Element<darcyIdx>& element,
473 const FVElementGeometry<darcyIdx>& fvGeometry,
474 const SubControlVolumeFace<darcyIdx>& scvf,
475 const VolumeVariables<darcyIdx>& volVars,
476 const typename Element<stokesIdx>::Geometry::GlobalCoordinate& couplingPhaseVelocity,
477 DarcysLaw) const
478 {
479 const auto darcyPhaseIdx = couplingPhaseIdx(darcyIdx);
480 const Scalar couplingPhaseCellCenterPressure = volVars.pressure(darcyPhaseIdx);
481 const Scalar couplingPhaseMobility = volVars.mobility(darcyPhaseIdx);
482 const Scalar couplingPhaseDensity = volVars.density(darcyPhaseIdx);
483 const auto K = volVars.permeability();
484
485 // A tpfa approximation yields (works if mobility != 0)
486 // v*n = -kr/mu*K * (gradP - rho*g)*n = mobility*(ti*(p_center - p_interface) + rho*n^TKg)
487 // -> p_interface = (1/mobility * (-v*n) + rho*n^TKg)/ti + p_center
488 // where v is the free-flow velocity (couplingPhaseVelocity)
489 const auto alpha = vtmv(scvf.unitOuterNormal(), K, couplingManager_.problem(darcyIdx).spatialParams().gravity(scvf.center()));
490
491 const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
492 const auto ti = computeTpfaTransmissibility(scvf, insideScv, K, 1.0);
493
494 return (-1/couplingPhaseMobility * (scvf.unitOuterNormal() * couplingPhaseVelocity) + couplingPhaseDensity * alpha)/ti
495 + couplingPhaseCellCenterPressure;
496 }
497
498
499private:
500 const CouplingManager& couplingManager_;
501
502};
503
508template<class MDTraits, class CouplingManager, bool enableEnergyBalance>
509class StokesDarcyCouplingDataImplementation<MDTraits, CouplingManager, enableEnergyBalance, false>
510: public StokesDarcyCouplingDataImplementationBase<MDTraits, CouplingManager>
511{
513 using Scalar = typename MDTraits::Scalar;
514 static constexpr auto stokesIdx = CouplingManager::stokesIdx;
515 static constexpr auto darcyIdx = CouplingManager::darcyIdx;
516 static constexpr auto stokesFaceIdx = CouplingManager::stokesFaceIdx;
517 static constexpr auto stokesCellCenterIdx = CouplingManager::stokesCellCenterIdx;
518
519 // the sub domain type tags
520 template<std::size_t id>
521 using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
522
523 template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
524 template<std::size_t id> using Element = typename GridGeometry<id>::GridView::template Codim<0>::Entity;
525 template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView;
526 template<std::size_t id> using SubControlVolumeFace = typename GridGeometry<id>::LocalView::SubControlVolumeFace;
527 template<std::size_t id> using SubControlVolume = typename GridGeometry<id>::LocalView::SubControlVolume;
528 template<std::size_t id> using Indices = typename GetPropType<SubDomainTypeTag<id>, Properties::ModelTraits>::Indices;
529 template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView;
530 template<std::size_t id> using ElementFaceVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridFaceVariables>::LocalView;
531 template<std::size_t id> using VolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::VolumeVariables;
532
534 "Darcy Model must not be compositional");
535
536 using DiffusionCoefficientAveragingType = typename StokesDarcyCouplingOptions::DiffusionCoefficientAveragingType;
537
538public:
539 using ParentType::ParentType;
540 using ParentType::couplingPhaseIdx;
541
545 Scalar massCouplingCondition(const Element<darcyIdx>& element,
546 const FVElementGeometry<darcyIdx>& fvGeometry,
547 const ElementVolumeVariables<darcyIdx>& darcyElemVolVars,
548 const SubControlVolumeFace<darcyIdx>& scvf) const
549 {
550 const auto& darcyContext = this->couplingManager().darcyCouplingContext(element, scvf);
551 const Scalar velocity = darcyContext.velocity * scvf.unitOuterNormal();
552 const Scalar darcyDensity = darcyElemVolVars[scvf.insideScvIdx()].density(couplingPhaseIdx(darcyIdx));
553 const Scalar stokesDensity = darcyContext.volVars.density();
554 const bool insideIsUpstream = velocity > 0.0;
555
556 return massFlux_(velocity, darcyDensity, stokesDensity, insideIsUpstream);
557 }
558
562 Scalar massCouplingCondition(const Element<stokesIdx>& element,
563 const FVElementGeometry<stokesIdx>& fvGeometry,
564 const ElementVolumeVariables<stokesIdx>& stokesElemVolVars,
565 const ElementFaceVariables<stokesIdx>& stokesElemFaceVars,
566 const SubControlVolumeFace<stokesIdx>& scvf) const
567 {
568 const auto& stokesContext = this->couplingManager().stokesCouplingContext(element, scvf);
569 const Scalar velocity = stokesElemFaceVars[scvf].velocitySelf();
570 const Scalar stokesDensity = stokesElemVolVars[scvf.insideScvIdx()].density();
571 const Scalar darcyDensity = stokesContext.volVars.density(couplingPhaseIdx(darcyIdx));
572 const bool insideIsUpstream = sign(velocity) == scvf.directionSign();
573
574 return massFlux_(velocity * scvf.directionSign(), stokesDensity, darcyDensity, insideIsUpstream);
575 }
576
580 template<bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0>
581 Scalar energyCouplingCondition(const Element<darcyIdx>& element,
582 const FVElementGeometry<darcyIdx>& fvGeometry,
583 const ElementVolumeVariables<darcyIdx>& darcyElemVolVars,
584 const SubControlVolumeFace<darcyIdx>& scvf,
585 const DiffusionCoefficientAveragingType diffCoeffAvgType = DiffusionCoefficientAveragingType::ffOnly) const
586 {
587 const auto& darcyContext = this->couplingManager().darcyCouplingContext(element, scvf);
588 const auto& darcyVolVars = darcyElemVolVars[scvf.insideScvIdx()];
589 const auto& stokesVolVars = darcyContext.volVars;
590
591 const Scalar velocity = darcyContext.velocity * scvf.unitOuterNormal();
592 const bool insideIsUpstream = velocity > 0.0;
593
594 return energyFlux_(darcyIdx, stokesIdx, fvGeometry, darcyContext.fvGeometry, scvf,
595 darcyVolVars, stokesVolVars, velocity, insideIsUpstream, diffCoeffAvgType);
596 }
597
601 template<bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0>
602 Scalar energyCouplingCondition(const Element<stokesIdx>& element,
603 const FVElementGeometry<stokesIdx>& fvGeometry,
604 const ElementVolumeVariables<stokesIdx>& stokesElemVolVars,
605 const ElementFaceVariables<stokesIdx>& stokesElemFaceVars,
606 const SubControlVolumeFace<stokesIdx>& scvf,
607 const DiffusionCoefficientAveragingType diffCoeffAvgType = DiffusionCoefficientAveragingType::ffOnly) const
608 {
609 const auto& stokesContext = this->couplingManager().stokesCouplingContext(element, scvf);
610 const auto& stokesVolVars = stokesElemVolVars[scvf.insideScvIdx()];
611 const auto& darcyVolVars = stokesContext.volVars;
612
613 const Scalar velocity = stokesElemFaceVars[scvf].velocitySelf();
614 const bool insideIsUpstream = sign(velocity) == scvf.directionSign();
615
616 return energyFlux_(stokesIdx, darcyIdx, fvGeometry, stokesContext.fvGeometry, scvf,
617 stokesVolVars, darcyVolVars, velocity * scvf.directionSign(), insideIsUpstream, diffCoeffAvgType);
618 }
619
620private:
621
625 Scalar massFlux_(const Scalar velocity,
626 const Scalar insideDensity,
627 const Scalar outSideDensity,
628 bool insideIsUpstream) const
629 {
630 return this->advectiveFlux(insideDensity, outSideDensity, velocity, insideIsUpstream);
631 }
632
636 template<std::size_t i, std::size_t j, bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0>
637 Scalar energyFlux_(Dune::index_constant<i> domainI,
638 Dune::index_constant<j> domainJ,
639 const FVElementGeometry<i>& insideFvGeometry,
640 const FVElementGeometry<j>& outsideFvGeometry,
641 const SubControlVolumeFace<i>& scvf,
642 const VolumeVariables<i>& insideVolVars,
643 const VolumeVariables<j>& outsideVolVars,
644 const Scalar velocity,
645 const bool insideIsUpstream,
646 const DiffusionCoefficientAveragingType diffCoeffAvgType) const
647 {
648 Scalar flux(0.0);
649
650 const auto& insideScv = (*scvs(insideFvGeometry).begin());
651 const auto& outsideScv = (*scvs(outsideFvGeometry).begin());
652
653 // convective fluxes
654 const Scalar insideTerm = insideVolVars.density(couplingPhaseIdx(domainI)) * insideVolVars.enthalpy(couplingPhaseIdx(domainI));
655 const Scalar outsideTerm = outsideVolVars.density(couplingPhaseIdx(domainJ)) * outsideVolVars.enthalpy(couplingPhaseIdx(domainJ));
656
657 flux += this->advectiveFlux(insideTerm, outsideTerm, velocity, insideIsUpstream);
658
659 flux += this->conductiveEnergyFlux_(domainI, domainJ, insideFvGeometry, outsideFvGeometry, scvf, insideScv, outsideScv, insideVolVars, outsideVolVars, diffCoeffAvgType);
660
661 return flux;
662 }
663
664};
665
670template<class MDTraits, class CouplingManager, bool enableEnergyBalance>
671class StokesDarcyCouplingDataImplementation<MDTraits, CouplingManager, enableEnergyBalance, true>
672: public StokesDarcyCouplingDataImplementationBase<MDTraits, CouplingManager>
673{
675 using Scalar = typename MDTraits::Scalar;
676 static constexpr auto stokesIdx = CouplingManager::stokesIdx;
677 static constexpr auto darcyIdx = CouplingManager::darcyIdx;
678 static constexpr auto stokesFaceIdx = CouplingManager::stokesFaceIdx;
679 static constexpr auto stokesCellCenterIdx = CouplingManager::stokesCellCenterIdx;
680
681 // the sub domain type tags
682 template<std::size_t id>
683 using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
684
685 template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
686 template<std::size_t id> using Element = typename GridGeometry<id>::GridView::template Codim<0>::Entity;
687 template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView;
688 template<std::size_t id> using SubControlVolumeFace = typename FVElementGeometry<id>::SubControlVolumeFace;
689 template<std::size_t id> using SubControlVolume = typename GridGeometry<id>::LocalView::SubControlVolume;
690 template<std::size_t id> using Indices = typename GetPropType<SubDomainTypeTag<id>, Properties::ModelTraits>::Indices;
691 template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView;
692 template<std::size_t id> using ElementFaceVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridFaceVariables>::LocalView;
693 template<std::size_t id> using VolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::VolumeVariables;
694 template<std::size_t id> using FluidSystem = GetPropType<SubDomainTypeTag<id>, Properties::FluidSystem>;
695
696 static constexpr auto numComponents = GetPropType<SubDomainTypeTag<stokesIdx>, Properties::ModelTraits>::numFluidComponents();
697 static constexpr auto replaceCompEqIdx = GetPropType<SubDomainTypeTag<stokesIdx>, Properties::ModelTraits>::replaceCompEqIdx();
698 static constexpr bool useMoles = GetPropType<SubDomainTypeTag<stokesIdx>, Properties::ModelTraits>::useMoles();
699 static constexpr auto referenceSystemFormulation = GetPropType<SubDomainTypeTag<stokesIdx>, Properties::MolecularDiffusionType>::referenceSystemFormulation();
700
701 static_assert(GetPropType<SubDomainTypeTag<darcyIdx>, Properties::ModelTraits>::numFluidComponents() == numComponents, "Both submodels must use the same number of components");
702 static_assert(getPropValue<SubDomainTypeTag<darcyIdx>, Properties::UseMoles>() == useMoles, "Both submodels must either use moles or not");
703 static_assert(getPropValue<SubDomainTypeTag<darcyIdx>, Properties::ReplaceCompEqIdx>() == replaceCompEqIdx, "Both submodels must use the same replaceCompEqIdx");
704 static_assert(GetPropType<SubDomainTypeTag<darcyIdx>, Properties::MolecularDiffusionType>::referenceSystemFormulation() == referenceSystemFormulation,
705 "Both submodels must use the same reference system formulation for diffusion");
706
707 using NumEqVector = Dune::FieldVector<Scalar, numComponents>;
708
709 using DiffusionCoefficientAveragingType = typename StokesDarcyCouplingOptions::DiffusionCoefficientAveragingType;
710
713 "Both submodels must use the same diffusion law.");
714
715 using ReducedComponentVector = Dune::FieldVector<Scalar, numComponents-1>;
716 using ReducedComponentMatrix = Dune::FieldMatrix<Scalar, numComponents-1, numComponents-1>;
717
719public:
720 using ParentType::ParentType;
721 using ParentType::couplingPhaseIdx;
722 using ParentType::couplingCompIdx;
723
727 NumEqVector massCouplingCondition(const Element<darcyIdx>& element,
728 const FVElementGeometry<darcyIdx>& fvGeometry,
729 const ElementVolumeVariables<darcyIdx>& darcyElemVolVars,
730 const SubControlVolumeFace<darcyIdx>& scvf,
731 const DiffusionCoefficientAveragingType diffCoeffAvgType = DiffusionCoefficientAveragingType::ffOnly) const
732 {
733 NumEqVector flux(0.0);
734 const auto& darcyContext = this->couplingManager().darcyCouplingContext(element, scvf);
735 const auto& darcyVolVars = darcyElemVolVars[scvf.insideScvIdx()];
736 const auto& stokesVolVars = darcyContext.volVars;
737 const auto& outsideScv = (*scvs(darcyContext.fvGeometry).begin());
738
739 const Scalar velocity = darcyContext.velocity * scvf.unitOuterNormal();
740 const bool insideIsUpstream = velocity > 0.0;
741
742 return massFlux_(darcyIdx, stokesIdx, fvGeometry,
743 scvf, darcyVolVars, stokesVolVars,
744 outsideScv, velocity, insideIsUpstream,
745 diffCoeffAvgType);
746 }
747
751 NumEqVector massCouplingCondition(const Element<stokesIdx>& element,
752 const FVElementGeometry<stokesIdx>& fvGeometry,
753 const ElementVolumeVariables<stokesIdx>& stokesElemVolVars,
754 const ElementFaceVariables<stokesIdx>& stokesElemFaceVars,
755 const SubControlVolumeFace<stokesIdx>& scvf,
756 const DiffusionCoefficientAveragingType diffCoeffAvgType = DiffusionCoefficientAveragingType::ffOnly) const
757 {
758 NumEqVector flux(0.0);
759 const auto& stokesContext = this->couplingManager().stokesCouplingContext(element, scvf);
760 const auto& stokesVolVars = stokesElemVolVars[scvf.insideScvIdx()];
761 const auto& darcyVolVars = stokesContext.volVars;
762 const auto& outsideScv = (*scvs(stokesContext.fvGeometry).begin());
763
764 const Scalar velocity = stokesElemFaceVars[scvf].velocitySelf();
765 const bool insideIsUpstream = sign(velocity) == scvf.directionSign();
766
767 return massFlux_(stokesIdx, darcyIdx, fvGeometry,
768 scvf, stokesVolVars, darcyVolVars,
769 outsideScv, velocity * scvf.directionSign(),
770 insideIsUpstream, diffCoeffAvgType);
771 }
772
776 template<bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0>
777 Scalar energyCouplingCondition(const Element<darcyIdx>& element,
778 const FVElementGeometry<darcyIdx>& fvGeometry,
779 const ElementVolumeVariables<darcyIdx>& darcyElemVolVars,
780 const SubControlVolumeFace<darcyIdx>& scvf,
781 const DiffusionCoefficientAveragingType diffCoeffAvgType = DiffusionCoefficientAveragingType::ffOnly) const
782 {
783 const auto& darcyContext = this->couplingManager().darcyCouplingContext(element, scvf);
784 const auto& darcyVolVars = darcyElemVolVars[scvf.insideScvIdx()];
785 const auto& stokesVolVars = darcyContext.volVars;
786
787 const Scalar velocity = darcyContext.velocity * scvf.unitOuterNormal();
788 const bool insideIsUpstream = velocity > 0.0;
789
790 return energyFlux_(darcyIdx, stokesIdx, fvGeometry, darcyContext.fvGeometry, scvf,
791 darcyVolVars, stokesVolVars, velocity, insideIsUpstream, diffCoeffAvgType);
792 }
793
797 template<bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0>
798 Scalar energyCouplingCondition(const Element<stokesIdx>& element,
799 const FVElementGeometry<stokesIdx>& fvGeometry,
800 const ElementVolumeVariables<stokesIdx>& stokesElemVolVars,
801 const ElementFaceVariables<stokesIdx>& stokesElemFaceVars,
802 const SubControlVolumeFace<stokesIdx>& scvf,
803 const DiffusionCoefficientAveragingType diffCoeffAvgType = DiffusionCoefficientAveragingType::ffOnly) const
804 {
805 const auto& stokesContext = this->couplingManager().stokesCouplingContext(element, scvf);
806 const auto& stokesVolVars = stokesElemVolVars[scvf.insideScvIdx()];
807 const auto& darcyVolVars = stokesContext.volVars;
808
809 const Scalar velocity = stokesElemFaceVars[scvf].velocitySelf();
810 const bool insideIsUpstream = sign(velocity) == scvf.directionSign();
811
812 return energyFlux_(stokesIdx, darcyIdx, fvGeometry, stokesContext.fvGeometry, scvf,
813 stokesVolVars, darcyVolVars, velocity * scvf.directionSign(), insideIsUpstream, diffCoeffAvgType);
814 }
815
816protected:
817
821 template<std::size_t i, std::size_t j>
822 NumEqVector massFlux_(Dune::index_constant<i> domainI,
823 Dune::index_constant<j> domainJ,
824 const FVElementGeometry<i>& insideFvGeometry,
825 const SubControlVolumeFace<i>& scvf,
826 const VolumeVariables<i>& insideVolVars,
827 const VolumeVariables<j>& outsideVolVars,
828 const SubControlVolume<j>& outsideScv,
829 const Scalar velocity,
830 const bool insideIsUpstream,
831 const DiffusionCoefficientAveragingType diffCoeffAvgType) const
832 {
833 NumEqVector flux(0.0);
834 NumEqVector diffusiveFlux(0.0);
835
836 auto moleOrMassFraction = [&](const auto& volVars, int phaseIdx, int compIdx)
837 { return useMoles ? volVars.moleFraction(phaseIdx, compIdx) : volVars.massFraction(phaseIdx, compIdx); };
838
839 auto moleOrMassDensity = [&](const auto& volVars, int phaseIdx)
840 { return useMoles ? volVars.molarDensity(phaseIdx) : volVars.density(phaseIdx); };
841
842 // treat the advective fluxes
843 auto insideTerm = [&](int compIdx)
844 { return moleOrMassFraction(insideVolVars, couplingPhaseIdx(domainI), compIdx) * moleOrMassDensity(insideVolVars, couplingPhaseIdx(domainI)); };
845
846 auto outsideTerm = [&](int compIdx)
847 { return moleOrMassFraction(outsideVolVars, couplingPhaseIdx(domainJ), compIdx) * moleOrMassDensity(outsideVolVars, couplingPhaseIdx(domainJ)); };
848
849
850 for (int compIdx = 0; compIdx < numComponents; ++compIdx)
851 {
852 const int domainICompIdx = couplingCompIdx(domainI, compIdx);
853 const int domainJCompIdx = couplingCompIdx(domainJ, compIdx);
854 flux[domainICompIdx] += this->advectiveFlux(insideTerm(domainICompIdx), outsideTerm(domainJCompIdx), velocity, insideIsUpstream);
855 }
856
857 // treat the diffusive fluxes
858 const auto& insideScv = insideFvGeometry.scv(scvf.insideScvIdx());
859
860 if (isFicksLaw)
861 diffusiveFlux += diffusiveMolecularFluxFicksLaw_(domainI, domainJ, scvf, insideScv, outsideScv, insideVolVars, outsideVolVars, diffCoeffAvgType);
862 else //maxwell stefan
863 diffusiveFlux += diffusiveMolecularFluxMaxwellStefan_(domainI, domainJ, scvf, insideScv, outsideScv, insideVolVars, outsideVolVars);
864
865 //convert to correct units if necessary
866 if (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged && useMoles)
867 {
868 for (int compIdx = 0; compIdx < numComponents; ++compIdx)
869 {
870 const int domainICompIdx = couplingCompIdx(domainI, compIdx);
871 diffusiveFlux[domainICompIdx] *= 1/FluidSystem<i>::molarMass(domainICompIdx);
872 }
873 }
874 if (referenceSystemFormulation == ReferenceSystemFormulation::molarAveraged && !useMoles)
875 {
876 for (int compIdx = 0; compIdx < numComponents; ++compIdx)
877 {
878 const int domainICompIdx = couplingCompIdx(domainI, compIdx);
879 diffusiveFlux[domainICompIdx] *= FluidSystem<i>::molarMass(domainICompIdx);
880 }
881 }
882
883 flux += diffusiveFlux;
884 // convert to total mass/mole balance, if set be user
885 if (replaceCompEqIdx < numComponents)
886 flux[replaceCompEqIdx] = std::accumulate(flux.begin(), flux.end(), 0.0);
887
888 return flux;
889 }
890
891 Scalar getComponentEnthalpy(const VolumeVariables<stokesIdx>& volVars, int phaseIdx, int compIdx) const
892 {
893 return FluidSystem<stokesIdx>::componentEnthalpy(volVars.fluidState(), 0, compIdx);
894 }
895
896 Scalar getComponentEnthalpy(const VolumeVariables<darcyIdx>& volVars, int phaseIdx, int compIdx) const
897 {
898 return FluidSystem<darcyIdx>::componentEnthalpy(volVars.fluidState(), phaseIdx, compIdx);
899 }
900
904 template<std::size_t i, std::size_t j>
905 NumEqVector diffusiveMolecularFluxMaxwellStefan_(Dune::index_constant<i> domainI,
906 Dune::index_constant<j> domainJ,
907 const SubControlVolumeFace<i>& scvfI,
908 const SubControlVolume<i>& scvI,
909 const SubControlVolume<j>& scvJ,
910 const VolumeVariables<i>& volVarsI,
911 const VolumeVariables<j>& volVarsJ) const
912 {
913 NumEqVector diffusiveFlux(0.0);
914
915 const Scalar insideDistance = this->getDistance_(scvI, scvfI);
916 const Scalar outsideDistance = this->getDistance_(scvJ, scvfI);
917
918 ReducedComponentVector moleFracInside(0.0);
919 ReducedComponentVector moleFracOutside(0.0);
920 ReducedComponentVector reducedFlux(0.0);
921 ReducedComponentMatrix reducedDiffusionMatrixInside(0.0);
922 ReducedComponentMatrix reducedDiffusionMatrixOutside(0.0);
923
924 //calculate the mole fraction vectors
925 for (int compIdx = 0; compIdx < numComponents-1; compIdx++)
926 {
927 const int domainICompIdx = couplingCompIdx(domainI, compIdx);
928 const int domainJCompIdx = couplingCompIdx(domainJ, compIdx);
929
930 assert(FluidSystem<i>::componentName(domainICompIdx) == FluidSystem<j>::componentName(domainJCompIdx));
931
932 //calculate x_inside
933 const Scalar xInside = volVarsI.moleFraction(couplingPhaseIdx(domainI), domainICompIdx);
934 //calculate outside molefraction with the respective transmissibility
935 const Scalar xOutside = volVarsJ.moleFraction(couplingPhaseIdx(domainJ), domainJCompIdx);
936 moleFracInside[domainICompIdx] = xInside;
937 moleFracOutside[domainICompIdx] = xOutside;
938 }
939
940 //now we have to do the tpfa: J_i = -J_j which leads to: J_i = -rho_i Bi^-1 omegai(x*-xi) with x* = (omegai rho_i Bi^-1 + omegaj rho_j Bj^-1)^-1 (xi omegai rho_i Bi^-1 + xj omegaj rho_j Bj^-1) with i inside and j outside.
941
942 //first set up the matrices containing the binary diffusion coefficients and mole fractions
943
944 //inside matrix. KIdx and LIdx are the indices for the k and l-th component, N for the n-th component
945 for (int compKIdx = 0; compKIdx < numComponents-1; compKIdx++)
946 {
947 const int domainICompKIdx = couplingCompIdx(domainI, compKIdx);
948 const Scalar xk = volVarsI.moleFraction(couplingPhaseIdx(domainI), domainICompKIdx);
949 const Scalar avgMolarMass = volVarsI.averageMolarMass(couplingPhaseIdx(domainI));
950 const Scalar Mn = FluidSystem<i>::molarMass(numComponents-1);
951 const Scalar tkn = volVarsI.effectiveDiffusionCoefficient(couplingPhaseIdx(domainI), domainICompKIdx, couplingCompIdx(domainI, numComponents-1));
952
953 // set the entries of the diffusion matrix of the diagonal
954 reducedDiffusionMatrixInside[domainICompKIdx][domainICompKIdx] += xk*avgMolarMass/(tkn*Mn);
955
956 for (int compLIdx = 0; compLIdx < numComponents; compLIdx++)
957 {
958 const int domainICompLIdx = couplingCompIdx(domainI, compLIdx);
959
960 // we don't want to calculate e.g. water in water diffusion
961 if (domainICompKIdx == domainICompLIdx)
962 continue;
963
964 const Scalar xl = volVarsI.moleFraction(couplingPhaseIdx(domainI), domainICompLIdx);
965 const Scalar Mk = FluidSystem<i>::molarMass(domainICompKIdx);
966 const Scalar Ml = FluidSystem<i>::molarMass(domainICompLIdx);
967 const Scalar tkl = volVarsI.effectiveDiffusionCoefficient(couplingPhaseIdx(domainI), domainICompKIdx, domainICompLIdx);
968 reducedDiffusionMatrixInside[domainICompKIdx][domainICompKIdx] += xl*avgMolarMass/(tkl*Mk);
969 reducedDiffusionMatrixInside[domainICompKIdx][domainICompLIdx] += xk*(avgMolarMass/(tkn*Mn) - avgMolarMass/(tkl*Ml));
970 }
971 }
972 //outside matrix
973 for (int compKIdx = 0; compKIdx < numComponents-1; compKIdx++)
974 {
975 const int domainJCompKIdx = couplingCompIdx(domainJ, compKIdx);
976 const int domainICompKIdx = couplingCompIdx(domainI, compKIdx);
977
978 const Scalar xk = volVarsJ.moleFraction(couplingPhaseIdx(domainJ), domainJCompKIdx);
979 const Scalar avgMolarMass = volVarsJ.averageMolarMass(couplingPhaseIdx(domainJ));
980 const Scalar Mn = FluidSystem<j>::molarMass(numComponents-1);
981 const Scalar tkn = volVarsJ.effectiveDiffusionCoefficient(couplingPhaseIdx(domainJ), domainJCompKIdx, couplingCompIdx(domainJ, numComponents-1));
982
983 // set the entries of the diffusion matrix of the diagonal
984 reducedDiffusionMatrixOutside[domainICompKIdx][domainICompKIdx] += xk*avgMolarMass/(tkn*Mn);
985
986 for (int compLIdx = 0; compLIdx < numComponents; compLIdx++)
987 {
988 const int domainJCompLIdx = couplingCompIdx(domainJ, compLIdx);
989 const int domainICompLIdx = couplingCompIdx(domainI, compLIdx);
990
991 // we don't want to calculate e.g. water in water diffusion
992 if (domainJCompLIdx == domainJCompKIdx)
993 continue;
994
995 const Scalar xl = volVarsJ.moleFraction(couplingPhaseIdx(domainJ), domainJCompLIdx);
996 const Scalar Mk = FluidSystem<j>::molarMass(domainJCompKIdx);
997 const Scalar Ml = FluidSystem<j>::molarMass(domainJCompLIdx);
998 const Scalar tkl = volVarsJ.effectiveDiffusionCoefficient(couplingPhaseIdx(domainJ), domainJCompKIdx, domainJCompLIdx);
999 reducedDiffusionMatrixOutside[domainICompKIdx][domainICompKIdx] += xl*avgMolarMass/(tkl*Mk);
1000 reducedDiffusionMatrixOutside[domainICompKIdx][domainICompLIdx] += xk*(avgMolarMass/(tkn*Mn) - avgMolarMass/(tkl*Ml));
1001 }
1002 }
1003
1004 const Scalar omegai = 1/insideDistance;
1005 const Scalar omegaj = 1/outsideDistance;
1006
1007 reducedDiffusionMatrixInside.invert();
1008 reducedDiffusionMatrixInside *= omegai*volVarsI.density(couplingPhaseIdx(domainI));
1009 reducedDiffusionMatrixOutside.invert();
1010 reducedDiffusionMatrixOutside *= omegaj*volVarsJ.density(couplingPhaseIdx(domainJ));
1011
1012 //in the helpervector we store the values for x*
1013 ReducedComponentVector helperVector(0.0);
1014 ReducedComponentVector gradientVectori(0.0);
1015 ReducedComponentVector gradientVectorj(0.0);
1016
1017 reducedDiffusionMatrixInside.mv(moleFracInside, gradientVectori);
1018 reducedDiffusionMatrixOutside.mv(moleFracOutside, gradientVectorj);
1019
1020 auto gradientVectorij = (gradientVectori + gradientVectorj);
1021
1022 //add the two matrixes to each other
1023 reducedDiffusionMatrixOutside += reducedDiffusionMatrixInside;
1024
1025 reducedDiffusionMatrixOutside.solve(helperVector, gradientVectorij);
1026
1027 //Bi^-1 omegai rho_i (x*-xi). As we previously multiplied rho_i and omega_i wit the insidematrix, this does not need to be done again
1028 helperVector -=moleFracInside;
1029 reducedDiffusionMatrixInside.mv(helperVector, reducedFlux);
1030
1031 reducedFlux *= -1;
1032
1033 for (int compIdx = 0; compIdx < numComponents-1; compIdx++)
1034 {
1035 const int domainICompIdx = couplingCompIdx(domainI, compIdx);
1036 diffusiveFlux[domainICompIdx] = reducedFlux[domainICompIdx];
1037 diffusiveFlux[couplingCompIdx(domainI, numComponents-1)] -= reducedFlux[domainICompIdx];
1038 }
1039 return diffusiveFlux;
1040 }
1041
1042 template<std::size_t i, std::size_t j>
1043 NumEqVector diffusiveMolecularFluxFicksLaw_(Dune::index_constant<i> domainI,
1044 Dune::index_constant<j> domainJ,
1045 const SubControlVolumeFace<i>& scvfI,
1046 const SubControlVolume<i>& scvI,
1047 const SubControlVolume<j>& scvJ,
1048 const VolumeVariables<i>& volVarsI,
1049 const VolumeVariables<j>& volVarsJ,
1050 const DiffusionCoefficientAveragingType diffCoeffAvgType) const
1051 {
1052 NumEqVector diffusiveFlux(0.0);
1053
1054 const Scalar rhoInside = massOrMolarDensity(volVarsI, referenceSystemFormulation, couplingPhaseIdx(domainI));
1055 const Scalar rhoOutside = massOrMolarDensity(volVarsJ, referenceSystemFormulation, couplingPhaseIdx(domainJ));
1056 const Scalar avgDensity = 0.5 * rhoInside + 0.5 * rhoOutside;
1057
1058 const Scalar insideDistance = this->getDistance_(scvI, scvfI);
1059 const Scalar outsideDistance = this->getDistance_(scvJ, scvfI);
1060
1061 for (int compIdx = 1; compIdx < numComponents; ++compIdx)
1062 {
1063 const int domainIMainCompIdx = couplingPhaseIdx(domainI);
1064 const int domainJMainCompIdx = couplingPhaseIdx(domainJ);
1065 const int domainICompIdx = couplingCompIdx(domainI, compIdx);
1066 const int domainJCompIdx = couplingCompIdx(domainJ, compIdx);
1067
1068 assert(FluidSystem<i>::componentName(domainICompIdx) == FluidSystem<j>::componentName(domainJCompIdx));
1069
1070 const Scalar massOrMoleFractionInside = massOrMoleFraction(volVarsI, referenceSystemFormulation, couplingPhaseIdx(domainI), domainICompIdx);
1071 const Scalar massOrMoleFractionOutside = massOrMoleFraction(volVarsJ, referenceSystemFormulation, couplingPhaseIdx(domainJ), domainJCompIdx);
1072
1073 const Scalar deltaMassOrMoleFrac = massOrMoleFractionOutside - massOrMoleFractionInside;
1074 const Scalar tij = this->transmissibility_(domainI,
1075 domainJ,
1076 insideDistance,
1077 outsideDistance,
1078 volVarsI.effectiveDiffusionCoefficient(couplingPhaseIdx(domainI), domainIMainCompIdx, domainICompIdx),
1079 volVarsJ.effectiveDiffusionCoefficient(couplingPhaseIdx(domainJ), domainJMainCompIdx, domainJCompIdx),
1080 diffCoeffAvgType);
1081 diffusiveFlux[domainICompIdx] += -avgDensity * tij * deltaMassOrMoleFrac;
1082 }
1083
1084 const Scalar cumulativeFlux = std::accumulate(diffusiveFlux.begin(), diffusiveFlux.end(), 0.0);
1085 diffusiveFlux[couplingCompIdx(domainI, 0)] = -cumulativeFlux;
1086
1087 return diffusiveFlux;
1088 }
1089
1093 template<std::size_t i, std::size_t j, bool isNI = enableEnergyBalance, typename std::enable_if_t<isNI, int> = 0>
1094 Scalar energyFlux_(Dune::index_constant<i> domainI,
1095 Dune::index_constant<j> domainJ,
1096 const FVElementGeometry<i>& insideFvGeometry,
1097 const FVElementGeometry<j>& outsideFvGeometry,
1098 const SubControlVolumeFace<i>& scvf,
1099 const VolumeVariables<i>& insideVolVars,
1100 const VolumeVariables<j>& outsideVolVars,
1101 const Scalar velocity,
1102 const bool insideIsUpstream,
1103 const DiffusionCoefficientAveragingType diffCoeffAvgType) const
1104 {
1105 Scalar flux(0.0);
1106
1107 const auto& insideScv = (*scvs(insideFvGeometry).begin());
1108 const auto& outsideScv = (*scvs(outsideFvGeometry).begin());
1109
1110 // convective fluxes
1111 const Scalar insideTerm = insideVolVars.density(couplingPhaseIdx(domainI)) * insideVolVars.enthalpy(couplingPhaseIdx(domainI));
1112 const Scalar outsideTerm = outsideVolVars.density(couplingPhaseIdx(domainJ)) * outsideVolVars.enthalpy(couplingPhaseIdx(domainJ));
1113
1114 flux += this->advectiveFlux(insideTerm, outsideTerm, velocity, insideIsUpstream);
1115
1116 flux += this->conductiveEnergyFlux_(domainI, domainJ, insideFvGeometry, outsideFvGeometry, scvf, insideScv, outsideScv, insideVolVars, outsideVolVars, diffCoeffAvgType);
1117
1118 auto diffusiveFlux = isFicksLaw ? diffusiveMolecularFluxFicksLaw_(domainI, domainJ, scvf, insideScv, outsideScv, insideVolVars, outsideVolVars, diffCoeffAvgType)
1119 : diffusiveMolecularFluxMaxwellStefan_(domainI, domainJ, scvf, insideScv, outsideScv, insideVolVars, outsideVolVars);
1120
1121
1122 for (int compIdx = 0; compIdx < diffusiveFlux.size(); ++compIdx)
1123 {
1124 const int domainICompIdx = couplingCompIdx(domainI, compIdx);
1125 const int domainJCompIdx = couplingCompIdx(domainJ, compIdx);
1126
1127 const bool insideDiffFluxIsUpstream = diffusiveFlux[domainICompIdx] > 0;
1128 const Scalar componentEnthalpy = insideDiffFluxIsUpstream ?
1129 getComponentEnthalpy(insideVolVars, couplingPhaseIdx(domainI), domainICompIdx)
1130 : getComponentEnthalpy(outsideVolVars, couplingPhaseIdx(domainJ), domainJCompIdx);
1131
1132 if (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged)
1133 flux += diffusiveFlux[domainICompIdx] * componentEnthalpy;
1134 else
1135 flux += diffusiveFlux[domainICompIdx] * FluidSystem<i>::molarMass(domainICompIdx) * componentEnthalpy;
1136 }
1137
1138 return flux;
1139 }
1140};
1141
1142} // end namespace Dumux
1143
1144#endif
Define some often used mathematical functions.
The available discretization methods in Dumux.
The reference frameworks and formulations available for splitting total fluxes into a advective and d...
Tensor::field_type computeTpfaTransmissibility(const SubControlVolumeFace &scvf, const SubControlVolume &scv, const Tensor &T, typename SubControlVolume::Traits::Scalar extrusionFactor)
Free function to evaluate the Tpfa transmissibility associated with the flux (in the form of flux = T...
Definition: tpfa/computetransmissibility.hh:47
VolumeVariables::PrimaryVariables::value_type massOrMoleFraction(const VolumeVariables &volVars, ReferenceSystemFormulation referenceSys, const int phaseIdx, const int compIdx)
returns the mass or mole fraction to be used in Fick's law based on the reference system
Definition: referencesystemformulation.hh:66
VolumeVariables::PrimaryVariables::value_type massOrMolarDensity(const VolumeVariables &volVars, ReferenceSystemFormulation referenceSys, const int phaseIdx)
evaluates the density to be used in Fick's law based on the reference system
Definition: referencesystemformulation.hh:55
constexpr Scalar harmonicMean(Scalar x, Scalar y, Scalar wx=1.0, Scalar wy=1.0) noexcept
Calculate the (weighted) harmonic mean of two scalar values.
Definition: math.hh:68
Dune::DenseMatrix< MAT >::value_type vtmv(const Dune::DenseVector< V1 > &v1, const Dune::DenseMatrix< MAT > &M, const Dune::DenseVector< V2 > &v2)
Evaluates the scalar product of a vector v2, projected by a matrix M, with a vector v1.
Definition: math.hh:849
constexpr int sign(const ValueType &value) noexcept
Sign or signum function.
Definition: math.hh:627
constexpr Scalar arithmeticMean(Scalar x, Scalar y, Scalar wx=1.0, Scalar wy=1.0) noexcept
Calculate the (weighted) arithmetic mean of two scalar values.
Definition: math.hh:49
Definition: adapt.hh:29
constexpr auto getPropValue()
get the value data member of a property
Definition: propertysystem.hh:153
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:149
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
Property whether to use moles or kg as amount unit for balance equations.
Definition: common/properties.hh:88
The component balance index that should be replaced by the total mass/mole balance.
Definition: common/properties.hh:90
Definition: common/properties.hh:101
The type for a global container for the volume variables.
Definition: common/properties.hh:108
The type for the calculation the advective fluxes.
Definition: common/properties.hh:140
The type for the calculation of the molecular diffusion fluxes.
Definition: common/properties.hh:144
The type of the fluid system to use.
Definition: common/properties.hh:155
Global vector containing face-related data.
Definition: common/properties.hh:217
Returns whether to normalize the pressure term in the momentum balance or not.
Definition: common/properties.hh:278
forward declaration of the method-specific implementation
Definition: flux/darcyslaw.hh:38
forward declaration of the method-specific implemetation
Definition: flux/box/fickslaw.hh:43
forward declare
Definition: forchheimerslaw.hh:38
This structs holds a set of options which allow to modify the Stokes-Darcy coupling mechanism during ...
Definition: couplingdata.hh:45
DiffusionCoefficientAveragingType
Defines which kind of averanging of diffusion coefficiencients (moleculat diffusion or thermal conduc...
Definition: couplingdata.hh:52
static DiffusionCoefficientAveragingType stringToEnum(DiffusionCoefficientAveragingType, const std::string &diffusionCoefficientAveragingType)
Convenience function to convert user input given as std::string to the corresponding enum class used ...
Definition: couplingdata.hh:60
This structs helps to check if the two sub models use the same fluidsystem. Specialization for the ca...
Definition: couplingdata.hh:85
static constexpr bool value
Definition: couplingdata.hh:87
This structs indicates that Fick's law is not used for diffusion.
Definition: couplingdata.hh:112
Helper struct to choose the correct index for phases and components. This is need if the porous-mediu...
Definition: couplingdata.hh:132
static constexpr auto couplingCompIdx(Dune::index_constant< i >, int coupledCompdIdx)
No adapter is used, just return the input index.
Definition: couplingdata.hh:156
static constexpr auto couplingPhaseIdx(Dune::index_constant< i >, int coupledPhaseIdx=0)
No adapter is used, just return the input index.
Definition: couplingdata.hh:149
static constexpr auto couplingPhaseIdx(Dune::index_constant< darcyIdx >, int coupledPhaseIdx=0)
The phase index of the porous-medium-flow model is given by the adapter fluidsytem (i....
Definition: couplingdata.hh:180
static constexpr auto couplingCompIdx(Dune::index_constant< stokesIdx >, int coupledCompdIdx)
The free-flow model does not need any change of the component index.
Definition: couplingdata.hh:186
static constexpr auto couplingCompIdx(Dune::index_constant< darcyIdx >, int coupledCompdIdx)
The component index of the porous-medium-flow model is mapped by the adapter fluidsytem.
Definition: couplingdata.hh:192
static constexpr auto couplingPhaseIdx(Dune::index_constant< stokesIdx >, int coupledPhaseIdx=0)
The free-flow model always uses phase index 0.
Definition: couplingdata.hh:174
Definition: couplingdata.hh:206
A base class which provides some common methods used for Stokes-Darcy coupling.
Definition: couplingdata.hh:224
auto darcyPermeability(const Element< stokesIdx > &element, const SubControlVolumeFace< stokesIdx > &scvf) const
Returns the intrinsic permeability of the coupled Darcy element.
Definition: couplingdata.hh:286
static constexpr auto couplingPhaseIdx(Dune::index_constant< i > id, int coupledPhaseIdx=0)
Returns the corresponding phase index needed for coupling.
Definition: couplingdata.hh:267
Scalar conductiveEnergyFlux_(Dune::index_constant< i > domainI, Dune::index_constant< j > domainJ, const FVElementGeometry< i > &fvGeometryI, const FVElementGeometry< j > &fvGeometryJ, const SubControlVolumeFace< i > &scvfI, const SubControlVolume< i > &scvI, const SubControlVolume< j > &scvJ, const VolumeVariables< i > &volVarsI, const VolumeVariables< j > &volVarsJ, const DiffusionCoefficientAveragingType diffCoeffAvgType) const
Returns the conductive energy flux acorss the interface.
Definition: couplingdata.hh:392
static constexpr auto couplingCompIdx(Dune::index_constant< i > id, int coupledCompdIdx)
Returns the corresponding component index needed for coupling.
Definition: couplingdata.hh:274
Scalar pressureAtInterface_(const Element< stokesIdx > &element, const SubControlVolumeFace< stokesIdx > &scvf, const ElementFaceVariables &elemFaceVars, const CouplingContext &context) const
Returns the pressure at the interface.
Definition: couplingdata.hh:422
Scalar advectiveFlux(const Scalar insideQuantity, const Scalar outsideQuantity, const Scalar volumeFlow, bool insideIsUpstream) const
Evaluate an advective flux across the interface and consider upwinding.
Definition: couplingdata.hh:333
Scalar computeCouplingPhasePressureAtInterface_(const Element< darcyIdx > &element, const FVElementGeometry< darcyIdx > &fvGeometry, const SubControlVolumeFace< darcyIdx > &scvf, const VolumeVariables< darcyIdx > &volVars, const typename Element< stokesIdx >::Geometry::GlobalCoordinate &couplingPhaseVelocity, DarcysLaw) const
Returns the pressure at the interface using Darcy's law for reconstruction.
Definition: couplingdata.hh:472
Scalar getDistance_(const Scv &scv, const Scvf &scvf) const
Returns the distance between an scvf and the corresponding scv center.
Definition: couplingdata.hh:383
Scalar momentumCouplingCondition(const Element< stokesIdx > &element, const FVElementGeometry< stokesIdx > &fvGeometry, const ElementVolumeVariables< stokesIdx > &stokesElemVolVars, const ElementFaceVariables &stokesElemFaceVars, const SubControlVolumeFace< stokesIdx > &scvf) const
Returns the momentum flux across the coupling boundary.
Definition: couplingdata.hh:300
Scalar computeCouplingPhasePressureAtInterface_(const Element< darcyIdx > &element, const FVElementGeometry< darcyIdx > &fvGeometry, const SubControlVolumeFace< darcyIdx > &scvf, const VolumeVariables< darcyIdx > &volVars, const typename Element< stokesIdx >::Geometry::GlobalCoordinate &couplingPhaseVelocity, ForchheimersLaw) const
Returns the pressure at the interface using Forchheimers's law for reconstruction.
Definition: couplingdata.hh:436
const CouplingManager & couplingManager() const
Returns a reference to the coupling manager.
Definition: couplingdata.hh:280
StokesDarcyCouplingDataImplementationBase(const CouplingManager &couplingmanager)
Definition: couplingdata.hh:261
Scalar transmissibility_(Dune::index_constant< i > domainI, Dune::index_constant< j > domainJ, const Scalar insideDistance, const Scalar outsideDistance, const Scalar avgQuantityI, const Scalar avgQuantityJ, const DiffusionCoefficientAveragingType diffCoeffAvgType) const
Returns the transmissibility used for either molecular diffusion or thermal conductivity.
Definition: couplingdata.hh:349
Scalar energyCouplingCondition(const Element< stokesIdx > &element, const FVElementGeometry< stokesIdx > &fvGeometry, const ElementVolumeVariables< stokesIdx > &stokesElemVolVars, const ElementFaceVariables< stokesIdx > &stokesElemFaceVars, const SubControlVolumeFace< stokesIdx > &scvf, const DiffusionCoefficientAveragingType diffCoeffAvgType=DiffusionCoefficientAveragingType::ffOnly) const
Returns the energy flux across the coupling boundary as seen from the free-flow domain.
Definition: couplingdata.hh:602
Scalar massCouplingCondition(const Element< stokesIdx > &element, const FVElementGeometry< stokesIdx > &fvGeometry, const ElementVolumeVariables< stokesIdx > &stokesElemVolVars, const ElementFaceVariables< stokesIdx > &stokesElemFaceVars, const SubControlVolumeFace< stokesIdx > &scvf) const
Returns the mass flux across the coupling boundary as seen from the free-flow domain.
Definition: couplingdata.hh:562
Scalar energyCouplingCondition(const Element< darcyIdx > &element, const FVElementGeometry< darcyIdx > &fvGeometry, const ElementVolumeVariables< darcyIdx > &darcyElemVolVars, const SubControlVolumeFace< darcyIdx > &scvf, const DiffusionCoefficientAveragingType diffCoeffAvgType=DiffusionCoefficientAveragingType::ffOnly) const
Returns the energy flux across the coupling boundary as seen from the Darcy domain.
Definition: couplingdata.hh:581
Scalar massCouplingCondition(const Element< darcyIdx > &element, const FVElementGeometry< darcyIdx > &fvGeometry, const ElementVolumeVariables< darcyIdx > &darcyElemVolVars, const SubControlVolumeFace< darcyIdx > &scvf) const
Returns the mass flux across the coupling boundary as seen from the Darcy domain.
Definition: couplingdata.hh:545
NumEqVector massCouplingCondition(const Element< stokesIdx > &element, const FVElementGeometry< stokesIdx > &fvGeometry, const ElementVolumeVariables< stokesIdx > &stokesElemVolVars, const ElementFaceVariables< stokesIdx > &stokesElemFaceVars, const SubControlVolumeFace< stokesIdx > &scvf, const DiffusionCoefficientAveragingType diffCoeffAvgType=DiffusionCoefficientAveragingType::ffOnly) const
Returns the mass flux across the coupling boundary as seen from the free-flow domain.
Definition: couplingdata.hh:751
Scalar energyCouplingCondition(const Element< darcyIdx > &element, const FVElementGeometry< darcyIdx > &fvGeometry, const ElementVolumeVariables< darcyIdx > &darcyElemVolVars, const SubControlVolumeFace< darcyIdx > &scvf, const DiffusionCoefficientAveragingType diffCoeffAvgType=DiffusionCoefficientAveragingType::ffOnly) const
Returns the energy flux across the coupling boundary as seen from the Darcy domain.
Definition: couplingdata.hh:777
Scalar getComponentEnthalpy(const VolumeVariables< darcyIdx > &volVars, int phaseIdx, int compIdx) const
Definition: couplingdata.hh:896
Scalar getComponentEnthalpy(const VolumeVariables< stokesIdx > &volVars, int phaseIdx, int compIdx) const
Definition: couplingdata.hh:891
NumEqVector massFlux_(Dune::index_constant< i > domainI, Dune::index_constant< j > domainJ, const FVElementGeometry< i > &insideFvGeometry, const SubControlVolumeFace< i > &scvf, const VolumeVariables< i > &insideVolVars, const VolumeVariables< j > &outsideVolVars, const SubControlVolume< j > &outsideScv, const Scalar velocity, const bool insideIsUpstream, const DiffusionCoefficientAveragingType diffCoeffAvgType) const
Evaluate the compositional mole/mass flux across the interface.
Definition: couplingdata.hh:822
Scalar energyFlux_(Dune::index_constant< i > domainI, Dune::index_constant< j > domainJ, const FVElementGeometry< i > &insideFvGeometry, const FVElementGeometry< j > &outsideFvGeometry, const SubControlVolumeFace< i > &scvf, const VolumeVariables< i > &insideVolVars, const VolumeVariables< j > &outsideVolVars, const Scalar velocity, const bool insideIsUpstream, const DiffusionCoefficientAveragingType diffCoeffAvgType) const
Evaluate the energy flux across the interface.
Definition: couplingdata.hh:1094
NumEqVector massCouplingCondition(const Element< darcyIdx > &element, const FVElementGeometry< darcyIdx > &fvGeometry, const ElementVolumeVariables< darcyIdx > &darcyElemVolVars, const SubControlVolumeFace< darcyIdx > &scvf, const DiffusionCoefficientAveragingType diffCoeffAvgType=DiffusionCoefficientAveragingType::ffOnly) const
Returns the mass flux across the coupling boundary as seen from the Darcy domain.
Definition: couplingdata.hh:727
Scalar energyCouplingCondition(const Element< stokesIdx > &element, const FVElementGeometry< stokesIdx > &fvGeometry, const ElementVolumeVariables< stokesIdx > &stokesElemVolVars, const ElementFaceVariables< stokesIdx > &stokesElemFaceVars, const SubControlVolumeFace< stokesIdx > &scvf, const DiffusionCoefficientAveragingType diffCoeffAvgType=DiffusionCoefficientAveragingType::ffOnly) const
Returns the energy flux across the coupling boundary as seen from the free-flow domain.
Definition: couplingdata.hh:798
NumEqVector diffusiveMolecularFluxMaxwellStefan_(Dune::index_constant< i > domainI, Dune::index_constant< j > domainJ, const SubControlVolumeFace< i > &scvfI, const SubControlVolume< i > &scvI, const SubControlVolume< j > &scvJ, const VolumeVariables< i > &volVarsI, const VolumeVariables< j > &volVarsJ) const
Evaluate the diffusive mole/mass flux across the interface.
Definition: couplingdata.hh:905
NumEqVector diffusiveMolecularFluxFicksLaw_(Dune::index_constant< i > domainI, Dune::index_constant< j > domainJ, const SubControlVolumeFace< i > &scvfI, const SubControlVolume< i > &scvI, const SubControlVolume< j > &scvJ, const VolumeVariables< i > &volVarsI, const VolumeVariables< j > &volVarsJ, const DiffusionCoefficientAveragingType diffCoeffAvgType) const
Definition: couplingdata.hh:1043
Definition: multidomain/couplingmanager.hh:46
const Problem< i > & problem(Dune::index_constant< i > domainIdx) const
Return a reference to the sub problem.
Definition: multidomain/couplingmanager.hh:264
Declares all properties used in Dumux.
Free functions to evaluate the transmissibilities associated with flux evaluations across sub-control...
The interface of the coupling manager for multi domain problems.