| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- | ||
| 2 | // vi: set et ts=4 sw=4 sts=4: | ||
| 3 | // | ||
| 4 | // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder | ||
| 5 | // SPDX-License-Identifier: GPL-3.0-or-later | ||
| 6 | // | ||
| 7 | /*! | ||
| 8 | * \file | ||
| 9 | * \ingroup PoreNetworkFlux | ||
| 10 | * \brief This file contains the data which is required to calculate | ||
| 11 | * the fluxes of the pore network model over a face of a finite volume. | ||
| 12 | */ | ||
| 13 | #ifndef DUMUX_FLUX_PNM_ADVECTION_HH | ||
| 14 | #define DUMUX_FLUX_PNM_ADVECTION_HH | ||
| 15 | |||
| 16 | #include <array> | ||
| 17 | #include <dumux/common/parameters.hh> | ||
| 18 | |||
| 19 | namespace Dumux::PoreNetwork::Detail { | ||
| 20 | |||
| 21 | template<class... TransmissibilityLawTypes> | ||
| 22 | struct Transmissibility : public TransmissibilityLawTypes... {}; | ||
| 23 | |||
| 24 | } // end namespace Dumux::PoreNetwork::Detail | ||
| 25 | |||
| 26 | namespace Dumux::PoreNetwork { | ||
| 27 | |||
| 28 | /*! | ||
| 29 | * \ingroup PoreNetworkFlux | ||
| 30 | * \brief Hagen–Poiseuille-type flux law to describe the advective flux for pore-network models. | ||
| 31 | */ | ||
| 32 | template<class ScalarT, class... TransmissibilityLawTypes> | ||
| 33 | class CreepingFlow | ||
| 34 | { | ||
| 35 | |||
| 36 | public: | ||
| 37 | //! Export the Scalar type | ||
| 38 | using Scalar = ScalarT; | ||
| 39 | |||
| 40 | //! Export the transmissibility law | ||
| 41 | using Transmissibility = Detail::Transmissibility<TransmissibilityLawTypes...>; | ||
| 42 | |||
| 43 | /*! | ||
| 44 | * \brief Returns the advective flux of a fluid phase | ||
| 45 | * across the given sub-control volume face (corresponding to a pore throat). | ||
| 46 | * \note The flux is given in N*m, and can be converted | ||
| 47 | * into a volume flux (m^3/s) or mass flux (kg/s) by applying an upwind scheme | ||
| 48 | * for the mobility (1/viscosity) or the product of density and mobility, respectively. | ||
| 49 | */ | ||
| 50 | template<class Problem, class Element, class FVElementGeometry, | ||
| 51 | class ElementVolumeVariables, class SubControlVolumeFace, class ElemFluxVarsCache> | ||
| 52 | 35151003 | static Scalar flux(const Problem& problem, | |
| 53 | const Element& element, | ||
| 54 | const FVElementGeometry& fvGeometry, | ||
| 55 | const ElementVolumeVariables& elemVolVars, | ||
| 56 | const SubControlVolumeFace& scvf, | ||
| 57 | const int phaseIdx, | ||
| 58 | const ElemFluxVarsCache& elemFluxVarsCache) | ||
| 59 | { | ||
| 60 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 35151003 times.
|
35151003 | const auto& fluxVarsCache = elemFluxVarsCache[scvf]; |
| 61 | |||
| 62 | // Get the inside and outside volume variables | ||
| 63 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 35151003 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 35151003 times.
|
70302006 | const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); |
| 64 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 35151003 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 35151003 times.
|
70302006 | const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx()); |
| 65 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 35151003 times.
|
35151003 | const auto& insideVolVars = elemVolVars[insideScv]; |
| 66 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 35151003 times.
|
35151003 | const auto& outsideVolVars = elemVolVars[outsideScv]; |
| 67 | |||
| 68 | // calculate the pressure difference | ||
| 69 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 35151003 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 35151003 times.
|
70302006 | const Scalar deltaP = insideVolVars.pressure(phaseIdx) - outsideVolVars.pressure(phaseIdx); |
| 70 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1641413 times.
|
35151003 | const Scalar transmissibility = fluxVarsCache.transmissibility(phaseIdx); |
| 71 | using std::isfinite; | ||
| 72 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 35151003 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 35151003 times.
|
70302006 | assert(isfinite(transmissibility)); |
| 73 | |||
| 74 | 35151003 | Scalar volumeFlow = transmissibility*deltaP; | |
| 75 | |||
| 76 | // add gravity term | ||
| 77 |
5/8✓ Branch 0 taken 14 times.
✓ Branch 1 taken 35150989 times.
✓ Branch 3 taken 14 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 14 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 14 times.
✗ Branch 10 not taken.
|
35151003 | static const bool enableGravity = getParamFromGroup<bool>(problem.paramGroup(), "Problem.EnableGravity"); |
| 78 |
2/2✓ Branch 0 taken 1872 times.
✓ Branch 1 taken 35149131 times.
|
35151003 | if (enableGravity) |
| 79 | { | ||
| 80 | 3744 | const Scalar rho = 0.5*insideVolVars.density(phaseIdx) + 0.5*outsideVolVars.density(phaseIdx); | |
| 81 | 9360 | const Scalar g = problem.spatialParams().gravity(scvf.center()) * scvf.unitOuterNormal(); | |
| 82 | |||
| 83 | // The transmissibility is with respect to the effective throat length (potentially dropping the pore body radii). | ||
| 84 | // For gravity, we need to consider the total throat length (i.e., the cell-center to cell-center distance). | ||
| 85 | // This might cause some inconsistencies TODO: is there a better way? | ||
| 86 | 1872 | volumeFlow += transmissibility * fluxVarsCache.poreToPoreDistance() * rho * g; | |
| 87 | } | ||
| 88 | |||
| 89 | 35151003 | return volumeFlow; | |
| 90 | } | ||
| 91 | |||
| 92 | /*! | ||
| 93 | * \brief Returns the throat conductivity | ||
| 94 | */ | ||
| 95 | template<class Problem, class Element, class FVElementGeometry, class ElementVolumeVariables, class FluxVariablesCache> | ||
| 96 | 1919640 | static Scalar calculateTransmissibility(const Problem& problem, | |
| 97 | const Element& element, | ||
| 98 | const FVElementGeometry& fvGeometry, | ||
| 99 | const typename FVElementGeometry::SubControlVolumeFace& scvf, | ||
| 100 | const ElementVolumeVariables& elemVolVars, | ||
| 101 | const FluxVariablesCache& fluxVarsCache, | ||
| 102 | const int phaseIdx) | ||
| 103 | { | ||
| 104 | if constexpr (ElementVolumeVariables::VolumeVariables::numFluidPhases() == 1) | ||
| 105 | 37204427 | return Transmissibility::singlePhaseTransmissibility(problem, element, fvGeometry, scvf, elemVolVars, fluxVarsCache, phaseIdx); | |
| 106 | else | ||
| 107 | { | ||
| 108 | static_assert(ElementVolumeVariables::VolumeVariables::numFluidPhases() == 2); | ||
| 109 | |||
| 110 | 1919640 | const auto& spatialParams = problem.spatialParams(); | |
| 111 | using FluidSystem = typename ElementVolumeVariables::VolumeVariables::FluidSystem; | ||
| 112 | 1919640 | const int wPhaseIdx = spatialParams.template wettingPhase<FluidSystem>(element, elemVolVars); | |
| 113 | 1919640 | const bool invaded = fluxVarsCache.invaded(); | |
| 114 | |||
| 115 |
2/2✓ Branch 0 taken 959820 times.
✓ Branch 1 taken 959820 times.
|
1919640 | if (phaseIdx == wPhaseIdx) |
| 116 | { | ||
| 117 |
2/2✓ Branch 0 taken 92457 times.
✓ Branch 1 taken 867363 times.
|
959820 | return invaded ? Transmissibility::wettingLayerTransmissibility(element, fvGeometry, scvf, fluxVarsCache) |
| 118 | 867363 | : Transmissibility::singlePhaseTransmissibility(problem, element, fvGeometry, scvf, elemVolVars, fluxVarsCache, phaseIdx); | |
| 119 | } | ||
| 120 | else // non-wetting phase | ||
| 121 | { | ||
| 122 |
2/2✓ Branch 0 taken 92457 times.
✓ Branch 1 taken 867363 times.
|
959820 | return invaded ? Transmissibility::nonWettingPhaseTransmissibility(element, fvGeometry, scvf, fluxVarsCache) |
| 123 | : 0.0; | ||
| 124 | } | ||
| 125 | } | ||
| 126 | } | ||
| 127 | |||
| 128 | template<class Problem, class Element, class FVElementGeometry, class ElementVolumeVariables, class FluxVariablesCache> | ||
| 129 | ✗ | static std::array<Scalar, 2> calculateTransmissibilities(const Problem& problem, | |
| 130 | const Element& element, | ||
| 131 | const FVElementGeometry& fvGeometry, | ||
| 132 | const ElementVolumeVariables& elemVolVars, | ||
| 133 | const typename FVElementGeometry::SubControlVolumeFace& scvf, | ||
| 134 | const FluxVariablesCache& fluxVarsCache) | ||
| 135 | { | ||
| 136 | static_assert(ElementVolumeVariables::VolumeVariables::numFluidPhases() == 1); | ||
| 137 | 605 | const Scalar t = calculateTransmissibility(problem, element, fvGeometry, scvf, elemVolVars, fluxVarsCache, 0); | |
| 138 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 603 times.
|
605 | return {t, -t}; |
| 139 | } | ||
| 140 | }; | ||
| 141 | |||
| 142 | /*! | ||
| 143 | * \ingroup PoreNetworkFlux | ||
| 144 | * \brief Non-creeping flow flux law to describe the advective flux for pore-network models based on | ||
| 145 | * El-Zehairy et al.(2019). | ||
| 146 | * https://doi.org/10.1016/j.advwatres.2019.103378 | ||
| 147 | */ | ||
| 148 | template <class ScalarT, class... TransmissibilityLawTypes> | ||
| 149 | class NonCreepingFlow | ||
| 150 | { | ||
| 151 | public: | ||
| 152 | //! Export the Scalar type | ||
| 153 | using Scalar = ScalarT; | ||
| 154 | |||
| 155 | //! Export the creeping flow transmissibility law | ||
| 156 | using TransmissibilityCreepingFlow = Detail::Transmissibility<TransmissibilityLawTypes...>; | ||
| 157 | |||
| 158 | //! Inherit transmissibility from creeping flow transmissibility law to cache non-creeping flow-related parameters | ||
| 159 | class Transmissibility: public TransmissibilityCreepingFlow | ||
| 160 | { | ||
| 161 | public: | ||
| 162 | class SinglePhaseCache : public TransmissibilityCreepingFlow::SinglePhaseCache | ||
| 163 | { | ||
| 164 | using ParentType = typename TransmissibilityCreepingFlow::SinglePhaseCache; | ||
| 165 | public: | ||
| 166 | using Scalar = ScalarT; | ||
| 167 | |||
| 168 | template<class Problem, class Element, class FVElementGeometry, class ElementVolumeVariables, class FluxVariablesCache> | ||
| 169 | 2554623 | void fill(const Problem& problem, | |
| 170 | const Element& element, | ||
| 171 | const FVElementGeometry& fvGeometry, | ||
| 172 | const typename FVElementGeometry::SubControlVolumeFace& scvf, | ||
| 173 | const ElementVolumeVariables& elemVolVars, | ||
| 174 | const FluxVariablesCache& fluxVarsCache, | ||
| 175 | const int phaseIdx) | ||
| 176 | { | ||
| 177 | 2554623 | ParentType::fill(problem, element, fvGeometry,scvf, elemVolVars, fluxVarsCache, phaseIdx); | |
| 178 | 2554623 | const auto elemSol = elementSolution(element, elemVolVars, fvGeometry); | |
| 179 | |||
| 180 |
2/2✓ Branch 0 taken 5109246 times.
✓ Branch 1 taken 2554623 times.
|
7663869 | for (const auto& scv : scvs(fvGeometry)) |
| 181 | { | ||
| 182 | 5109246 | const auto localIdx = scv.indexInElement(); | |
| 183 | 10218492 | const Scalar throatToPoreAreaRatio = fluxVarsCache.throatCrossSectionalArea() / problem.spatialParams().poreCrossSectionalArea(element, scv, elemSol); | |
| 184 | |||
| 185 | // dimensionless momentum coefficient | ||
| 186 | 5109246 | const Scalar momentumCoefficient = problem.spatialParams().momentumCoefficient(element, scv, elemSol); | |
| 187 | |||
| 188 | // dimensionless kinetik-energy coefficient | ||
| 189 | 5109246 | const Scalar kineticEnergyCoefficient = problem.spatialParams().kineticEnergyCoefficient(element, scv, elemSol); | |
| 190 | |||
| 191 | 15327738 | expansionCoefficient_[localIdx] = getExpansionCoefficient_(throatToPoreAreaRatio, momentumCoefficient, kineticEnergyCoefficient); | |
| 192 | 5109246 | contractionCoefficient_[localIdx] = getContractionCoefficient_(throatToPoreAreaRatio, momentumCoefficient, kineticEnergyCoefficient); | |
| 193 | } | ||
| 194 | 2554623 | } | |
| 195 | |||
| 196 | Scalar expansionCoefficient(int downstreamIdx) const | ||
| 197 | 5109246 | { return expansionCoefficient_[downstreamIdx]; } | |
| 198 | |||
| 199 | Scalar contractionCoefficient(int upstreamIdx) const | ||
| 200 | 5109246 | { return contractionCoefficient_[upstreamIdx]; } | |
| 201 | |||
| 202 | private: | ||
| 203 | |||
| 204 | ✗ | Scalar getExpansionCoefficient_(const Scalar throatToPoreAreaRatio, const Scalar momentumCoefficient, const Scalar kineticEnergyCoefficient) const | |
| 205 | { | ||
| 206 | 10218492 | Scalar expansionCoefficient = throatToPoreAreaRatio * throatToPoreAreaRatio * (2 * momentumCoefficient - kineticEnergyCoefficient) | |
| 207 | 5109246 | + kineticEnergyCoefficient - 2 * momentumCoefficient * throatToPoreAreaRatio | |
| 208 | 5109246 | - (1 - throatToPoreAreaRatio * throatToPoreAreaRatio); | |
| 209 | |||
| 210 | ✗ | return expansionCoefficient; | |
| 211 | } | ||
| 212 | |||
| 213 | ✗ | Scalar getContractionCoefficient_(const Scalar throatToPoreAreaRatio, const Scalar momentumCoefficient, const Scalar kineticEnergyCoefficient) const | |
| 214 | { | ||
| 215 | 10218492 | const Scalar contractionAreaRatio = getContractionAreaRatio_(throatToPoreAreaRatio); | |
| 216 | 15327738 | Scalar contractionCoefficient = (1 - (throatToPoreAreaRatio * throatToPoreAreaRatio * kineticEnergyCoefficient - 2 * momentumCoefficient /*+1-throatToPoreAreaRatio*throatToPoreAreaRatio*/) | |
| 217 | 10218492 | * contractionAreaRatio * contractionAreaRatio - 2 * contractionAreaRatio) / (contractionAreaRatio * contractionAreaRatio); | |
| 218 | |||
| 219 | ✗ | return contractionCoefficient; | |
| 220 | } | ||
| 221 | |||
| 222 | ✗ | Scalar getContractionAreaRatio_(const Scalar throatToPoreAreaRatio) const | |
| 223 | { | ||
| 224 | 5109246 | return 1.0-(1.0-throatToPoreAreaRatio)/(2.08*(1.0-throatToPoreAreaRatio)+0.5371); | |
| 225 | } | ||
| 226 | |||
| 227 | std::array<Scalar, 2> expansionCoefficient_; | ||
| 228 | std::array<Scalar, 2> contractionCoefficient_; | ||
| 229 | }; | ||
| 230 | }; | ||
| 231 | |||
| 232 | /*! | ||
| 233 | * \brief Returns the advective flux of a fluid phase | ||
| 234 | * across the given sub-control volume face (corresponding to a pore throat). | ||
| 235 | * \note The flux is given in N*m, and can be converted | ||
| 236 | * into a volume flux (m^3/s) or mass flux (kg/s) by applying an upwind scheme | ||
| 237 | * for the mobility (1/viscosity) or the product of density and mobility, respectively. | ||
| 238 | */ | ||
| 239 | template<class Problem, class Element, class FVElementGeometry, | ||
| 240 | class ElementVolumeVariables, class SubControlVolumeFace, class ElemFluxVarsCache> | ||
| 241 | 2554623 | static Scalar flux(const Problem& problem, | |
| 242 | const Element& element, | ||
| 243 | const FVElementGeometry& fvGeometry, | ||
| 244 | const ElementVolumeVariables& elemVolVars, | ||
| 245 | const SubControlVolumeFace& scvf, | ||
| 246 | const int phaseIdx, | ||
| 247 | const ElemFluxVarsCache& elemFluxVarsCache) | ||
| 248 | { | ||
| 249 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2554619 times.
|
2554623 | const auto& fluxVarsCache = elemFluxVarsCache[scvf]; |
| 250 | |||
| 251 | // Get the inside and outside volume variables | ||
| 252 |
4/4✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2554619 times.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 2554619 times.
|
5109246 | const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); |
| 253 |
4/4✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2554619 times.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 2554619 times.
|
5109246 | const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx()); |
| 254 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2554619 times.
|
2554623 | const auto& insideVolVars = elemVolVars[insideScv]; |
| 255 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2554619 times.
|
2554623 | const auto& outsideVolVars = elemVolVars[outsideScv]; |
| 256 | |||
| 257 | // calculate the pressure difference | ||
| 258 |
4/4✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2554619 times.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 2554619 times.
|
5109246 | Scalar deltaP = insideVolVars.pressure(phaseIdx) - outsideVolVars.pressure(phaseIdx); |
| 259 | |||
| 260 | // add gravity term | ||
| 261 |
6/8✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2554619 times.
✓ Branch 3 taken 3 times.
✓ Branch 4 taken 1 times.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 3 times.
✗ Branch 10 not taken.
|
2554623 | static const bool enableGravity = getParamFromGroup<bool>(problem.paramGroup(), "Problem.EnableGravity"); |
| 262 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2554623 times.
|
2554623 | if (enableGravity) |
| 263 | { | ||
| 264 | ✗ | const Scalar rho = 0.5*insideVolVars.density(phaseIdx) + 0.5*outsideVolVars.density(phaseIdx); | |
| 265 | |||
| 266 | // projected distance between pores in gravity field direction | ||
| 267 | ✗ | const auto poreToPoreVector = fluxVarsCache.poreToPoreDistance() * scvf.unitOuterNormal(); | |
| 268 | ✗ | const auto& gravityVector = problem.spatialParams().gravity(scvf.center()); | |
| 269 | |||
| 270 | ✗ | deltaP += rho * (gravityVector * poreToPoreVector); | |
| 271 | } | ||
| 272 | 2554623 | const Scalar creepingFlowTransmissibility = fluxVarsCache.transmissibility(phaseIdx); | |
| 273 | 2554623 | const Scalar throatCrossSectionalArea = fluxVarsCache.throatCrossSectionalArea(); | |
| 274 | using std::isfinite; | ||
| 275 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 2554623 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2554623 times.
|
5109246 | assert(isfinite(creepingFlowTransmissibility)); |
| 276 | |||
| 277 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 2554623 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2554623 times.
|
5109246 | assert(scvf.insideScvIdx() == 0); |
| 278 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 2554623 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2554623 times.
|
5109246 | assert(scvf.outsideScvIdx() == 1); |
| 279 | |||
| 280 | // determine the flow direction to predict contraction and expansion | ||
| 281 |
6/6✓ Branch 0 taken 955632 times.
✓ Branch 1 taken 1598991 times.
✓ Branch 2 taken 955632 times.
✓ Branch 3 taken 1598991 times.
✓ Branch 4 taken 955632 times.
✓ Branch 5 taken 1598991 times.
|
2554623 | const auto [upstreamIdx, downstreamIdx] = deltaP > 0 ? std::pair(0, 1) : std::pair(1, 0); |
| 282 | |||
| 283 |
4/4✓ Branch 0 taken 955632 times.
✓ Branch 1 taken 1598991 times.
✓ Branch 2 taken 955632 times.
✓ Branch 3 taken 1598991 times.
|
5109246 | const Scalar contractionCoefficient = fluxVarsCache.singlePhaseFlowVariables().contractionCoefficient(upstreamIdx); |
| 284 |
4/4✓ Branch 0 taken 955632 times.
✓ Branch 1 taken 1598991 times.
✓ Branch 2 taken 955632 times.
✓ Branch 3 taken 1598991 times.
|
5109246 | const Scalar expansionCoefficient = fluxVarsCache.singlePhaseFlowVariables().expansionCoefficient(downstreamIdx); |
| 285 |
4/4✓ Branch 0 taken 955632 times.
✓ Branch 1 taken 1598991 times.
✓ Branch 2 taken 955632 times.
✓ Branch 3 taken 1598991 times.
|
5109246 | const Scalar mu = elemVolVars[upstreamIdx].viscosity(); |
| 286 | |||
| 287 | //! El-Zehairy et al.(2019): Eq.(14): | ||
| 288 | //! A0Q^2 + B0Q + C0 = 0 where Q is volumetric flowrate, A0 = [Kc + Ke] * rho /[2aij^2], Kc and Ke are contraction and expansion coefficient respectively, | ||
| 289 | //! aij is cross sectional area of the throat, B0 = 1/K0 (K0 is trnasmissibility used in paper) and C0 = -deltaP | ||
| 290 | //! In our implementation viscosity of fluid will be included using upwinding later. Thus, creepingFlowTransmissibility = mu * K0 and | ||
| 291 | //! the volumetric flowrate calculated here q = mu * Q. Substitution of new variables into El-Zehairy et al.(2019), Eq.(14) gives: | ||
| 292 | //! Aq^2 + Bq + C = 0 where | ||
| 293 | //! A = A0 / mu^2, B = B0/mu = 1/creepingFlowTransmissibility and C = C0 | ||
| 294 | //! attention: the q, volumetric flowrate, calculated here is always positive and its sign needs to be determined based on flow direction | ||
| 295 | //! this approach is taken to prevent the term under the square root (discriminant) becoming negative | ||
| 296 |
8/8✓ Branch 0 taken 955632 times.
✓ Branch 1 taken 1598991 times.
✓ Branch 2 taken 955632 times.
✓ Branch 3 taken 1598991 times.
✓ Branch 4 taken 955632 times.
✓ Branch 5 taken 1598991 times.
✓ Branch 6 taken 955632 times.
✓ Branch 7 taken 1598991 times.
|
10218492 | const Scalar A = (contractionCoefficient * elemVolVars[upstreamIdx].density() + expansionCoefficient * elemVolVars[downstreamIdx].density()) |
| 297 | 2554623 | / (2.0 * mu * mu * throatCrossSectionalArea * throatCrossSectionalArea); | |
| 298 | 2554623 | const Scalar B = 1.0/creepingFlowTransmissibility; | |
| 299 | using std::abs; | ||
| 300 |
2/2✓ Branch 0 taken 955632 times.
✓ Branch 1 taken 1598991 times.
|
2554623 | const Scalar C = -abs(deltaP); |
| 301 | |||
| 302 | using std::sqrt; | ||
| 303 | 2554623 | const auto discriminant = B*B - 4*A*C; | |
| 304 | 2554623 | const auto q = (-B + sqrt(discriminant)) / (2*A); | |
| 305 | |||
| 306 | //! give the volume flowrate proper sign based on flow direction. | ||
| 307 |
2/2✓ Branch 0 taken 955632 times.
✓ Branch 1 taken 1598991 times.
|
2554623 | if (upstreamIdx == 0) |
| 308 | return q; | ||
| 309 | else | ||
| 310 | 955632 | return -q; | |
| 311 | |||
| 312 | } | ||
| 313 | |||
| 314 | template<class Problem, class Element, class FVElementGeometry, class ElementVolumeVariables, class FluxVariablesCache> | ||
| 315 | ✗ | static Scalar calculateTransmissibility(const Problem& problem, | |
| 316 | const Element& element, | ||
| 317 | const FVElementGeometry& fvGeometry, | ||
| 318 | const typename FVElementGeometry::SubControlVolumeFace& scvf, | ||
| 319 | const ElementVolumeVariables& elemVolVars, | ||
| 320 | const FluxVariablesCache& fluxVarsCache, | ||
| 321 | const int phaseIdx) | ||
| 322 | { | ||
| 323 | static_assert(ElementVolumeVariables::VolumeVariables::numFluidPhases() == 1); | ||
| 324 | 2554623 | return Transmissibility::singlePhaseTransmissibility(problem, element, fvGeometry, scvf, elemVolVars, fluxVarsCache, phaseIdx); | |
| 325 | } | ||
| 326 | |||
| 327 | template<class Problem, class Element, class FVElementGeometry, class ElementVolumeVariables, class FluxVariablesCache> | ||
| 328 | static std::array<Scalar, 2> calculateTransmissibilities(const Problem& problem, | ||
| 329 | const Element& element, | ||
| 330 | const FVElementGeometry& fvGeometry, | ||
| 331 | const ElementVolumeVariables& elemVolVars, | ||
| 332 | const typename FVElementGeometry::SubControlVolumeFace& scvf, | ||
| 333 | const FluxVariablesCache& fluxVarsCache) | ||
| 334 | { | ||
| 335 | static_assert(ElementVolumeVariables::VolumeVariables::numFluidPhases() == 1); | ||
| 336 | const Scalar t = calculateTransmissibility(problem, element, fvGeometry, scvf, elemVolVars, fluxVarsCache, 0); | ||
| 337 | return {t, -t}; | ||
| 338 | } | ||
| 339 | |||
| 340 | }; | ||
| 341 | |||
| 342 | } // end namespace Dumux | ||
| 343 | |||
| 344 | #endif | ||
| 345 |