| 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 Fluidmatrixinteractions | ||
| 10 | * \ingroup PoreNetworkModels | ||
| 11 | * \brief Implementation of the transmissibility laws for throats | ||
| 12 | */ | ||
| 13 | #ifndef DUMUX_PNM_THROAT_TRANSMISSIBILITY_2P_HH | ||
| 14 | #define DUMUX_PNM_THROAT_TRANSMISSIBILITY_2P_HH | ||
| 15 | |||
| 16 | #include <dumux/porenetwork/common/throatproperties.hh> | ||
| 17 | #include "emptycache.hh" | ||
| 18 | |||
| 19 | namespace Dumux::PoreNetwork { | ||
| 20 | |||
| 21 | namespace WettingLayerTransmissibility { | ||
| 22 | |||
| 23 | struct CreviceResistanceFactorZhou | ||
| 24 | { | ||
| 25 | /*! | ||
| 26 | * \brief Returns the crevice resistance factor used for calculating the w-phase conductance in an invaded pore throat | ||
| 27 | * | ||
| 28 | * \param alpha The corner half angle | ||
| 29 | * \param theta The contact angle | ||
| 30 | * \param f The boundary condition factor for the fluid interface (0 = free boundary) | ||
| 31 | * | ||
| 32 | * Zhou et al. (1997), eq. 19; Singh & Mohanty (2002), eq 8 | ||
| 33 | */ | ||
| 34 | template<class Scalar> | ||
| 35 | 7678560 | static Scalar beta(const Scalar alpha, const Scalar theta, const Scalar f = 0) | |
| 36 | { | ||
| 37 | using std::sin; | ||
| 38 | using std::cos; | ||
| 39 | using std::tan; | ||
| 40 | 7678560 | const Scalar sinAlpha = sin(alpha); | |
| 41 | 7678560 | const Scalar sinSum = sin(alpha + theta); | |
| 42 | 7678560 | const Scalar cosSum = cos(alpha + theta); | |
| 43 | 7678560 | const Scalar phi1 = cosSum*cosSum + cosSum*sinSum*tan(alpha); | |
| 44 | 7678560 | const Scalar phi2 = 1 - theta/(M_PI/2 - alpha); | |
| 45 | 7678560 | const Scalar phi3 = cosSum / cos(alpha); | |
| 46 | 7678560 | const Scalar B = (M_PI/2 - alpha)*tan(alpha); | |
| 47 | |||
| 48 | 7678560 | Scalar result = 12*sinAlpha*sinAlpha*(1-B)*(1-B)*(phi1 - B*phi2)*(phi3 + f*B*phi2)*(phi3 + f*B*phi2); | |
| 49 | 7678560 | result /= (1-sinAlpha)*(1-sinAlpha)*B*B*(phi1 - B*phi2)*(phi1 - B*phi2)*(phi1 - B*phi2); | |
| 50 | 7678560 | return result; | |
| 51 | } | ||
| 52 | }; | ||
| 53 | |||
| 54 | |||
| 55 | template<class Scalar, class CreviceResistanceFactor = CreviceResistanceFactorZhou> | ||
| 56 | struct RansohoffRadke | ||
| 57 | { | ||
| 58 |
1/4✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 6 taken 31352 times.
✗ Branch 7 not taken.
|
316534 | class WettingLayerCache |
| 59 | { | ||
| 60 | using NumCornerVector = Dune::ReservedVector<Scalar, 4>; | ||
| 61 | public: | ||
| 62 | template<class Problem, class Element, class FVElementGeometry, class ElementVolumeVariables, class FluxVariablesCache> | ||
| 63 | 1919640 | void fill(const Problem& problem, | |
| 64 | const Element& element, | ||
| 65 | const FVElementGeometry& fvGeometry, | ||
| 66 | const typename FVElementGeometry::SubControlVolumeFace& scvf, | ||
| 67 | const ElementVolumeVariables& elemVolVars, | ||
| 68 | const FluxVariablesCache& fluxVarsCache, | ||
| 69 | const int phaseIdx) | ||
| 70 | { | ||
| 71 | 1919640 | const auto& spatialParams = problem.spatialParams(); | |
| 72 | 1919640 | const auto& cornerHalfAngles = spatialParams.cornerHalfAngles(element); | |
| 73 | 1919640 | const Scalar contactAngle = spatialParams.contactAngle(element, elemVolVars); | |
| 74 | |||
| 75 | 3839280 | beta_.clear(); beta_.resize(cornerHalfAngles.size()); | |
| 76 |
2/2✓ Branch 0 taken 7678560 times.
✓ Branch 1 taken 1919640 times.
|
9598200 | for (int i = 0; i< cornerHalfAngles.size(); ++i) |
| 77 | 23035680 | beta_[i] = CreviceResistanceFactor::beta(cornerHalfAngles[i], contactAngle); | |
| 78 | 1919640 | } | |
| 79 | |||
| 80 | Scalar creviceResistanceFactor(const int cornerIdx) const | ||
| 81 | 739656 | { return beta_[cornerIdx]; } | |
| 82 | |||
| 83 | private: | ||
| 84 | NumCornerVector beta_; | ||
| 85 | }; | ||
| 86 | |||
| 87 | /*! | ||
| 88 | * \brief Returns the integral conductivity of all wetting layers occupying the corners of a throat | ||
| 89 | */ | ||
| 90 | template<class Element, class FVElementGeometry, class FluxVariablesCache> | ||
| 91 | 92457 | static Scalar wettingLayerTransmissibility(const Element& element, | |
| 92 | const FVElementGeometry& fvGeometry, | ||
| 93 | const typename FVElementGeometry::SubControlVolumeFace& scvf, | ||
| 94 | const FluxVariablesCache& fluxVarsCache) | ||
| 95 | { | ||
| 96 | 92457 | const Scalar throatLength = fluxVarsCache.throatLength(); | |
| 97 | 92457 | const Scalar rC = fluxVarsCache.curvatureRadius(); | |
| 98 | |||
| 99 | 184914 | const auto eIdx = fvGeometry.gridGeometry().elementMapper().index(element); | |
| 100 |
1/2✓ Branch 0 taken 92457 times.
✗ Branch 1 not taken.
|
92457 | const auto shape = fvGeometry.gridGeometry().throatCrossSectionShape(eIdx); |
| 101 | 92457 | const auto numCorners = Throat::numCorners(shape); | |
| 102 | |||
| 103 | // treat the wetting film layer in each corner of the throat individually (might have different corner half-angle and beta) | ||
| 104 | 92457 | Scalar result = 0.0; | |
| 105 |
2/2✓ Branch 1 taken 369828 times.
✓ Branch 2 taken 92457 times.
|
462285 | for (int i = 0; i < numCorners; ++i) |
| 106 | 1479312 | result += fluxVarsCache.wettingLayerCrossSectionalArea(i) * rC*rC / (throatLength*fluxVarsCache.wettingLayerFlowVariables().creviceResistanceFactor(i)); | |
| 107 | |||
| 108 | 92457 | return result; | |
| 109 | } | ||
| 110 | }; | ||
| 111 | } // end namespace WettingLayerTransmissibility | ||
| 112 | |||
| 113 | namespace NonWettingPhaseTransmissibility { | ||
| 114 | |||
| 115 | //! TODO: evaluate and maybe remove | ||
| 116 | template<class Scalar> | ||
| 117 | struct Mogensen | ||
| 118 | { | ||
| 119 | using NonWettingPhaseCache = EmptyCache; | ||
| 120 | |||
| 121 | template<class Element, class FVElementGeometry, class FluxVariablesCache> | ||
| 122 | static Scalar nonWettingPhaseTransmissibility(const Element& element, | ||
| 123 | const FVElementGeometry& fvGeometry, | ||
| 124 | const typename FVElementGeometry::SubControlVolumeFace& scvf, | ||
| 125 | const FluxVariablesCache& fluxVarsCache) | ||
| 126 | { | ||
| 127 | // Mogensen et al. (1999), does not really revover the single phase value | ||
| 128 | using std::sqrt; | ||
| 129 | const Scalar throatLength = fluxVarsCache.throatLength(); | ||
| 130 | const auto nPhaseIdx = fluxVarsCache.nPhaseIdx(); | ||
| 131 | const Scalar aNw = fluxVarsCache.throatCrossSectionalArea(nPhaseIdx); | ||
| 132 | const Scalar rEff = 0.5*(sqrt(aNw / M_PI) + fluxVarsCache.throatInscribedRadius()); | ||
| 133 | const Scalar result = M_PI/(8*throatLength) * rEff*rEff*rEff*rEff; | ||
| 134 | return result; | ||
| 135 | } | ||
| 136 | }; | ||
| 137 | |||
| 138 | //! TODO: evaluate and maybe remove | ||
| 139 | template<class Scalar, class SinglePhaseTransmissibilityLaw> | ||
| 140 | struct Valvatne | ||
| 141 | { | ||
| 142 | using NonWettingPhaseCache = EmptyCache; | ||
| 143 | |||
| 144 | template<class Element, class FVElementGeometry, class FluxVariablesCache> | ||
| 145 | static Scalar nonWettingPhaseTransmissibility(const Element& element, | ||
| 146 | const FVElementGeometry& fvGeometry, | ||
| 147 | const typename FVElementGeometry::SubControlVolumeFace& scvf, | ||
| 148 | const FluxVariablesCache& fluxVarsCache) | ||
| 149 | { | ||
| 150 | // Tora et al. (2012), also does not fully recover single-phase value, but is closer | ||
| 151 | using std::sqrt; | ||
| 152 | const auto nPhaseIdx = fluxVarsCache.nPhaseIdx(); | ||
| 153 | const Scalar aNw = fluxVarsCache.throatCrossSectionalArea(nPhaseIdx); | ||
| 154 | const Scalar aTot = fluxVarsCache.throatCrossSectionalArea(); | ||
| 155 | |||
| 156 | const Scalar result = SinglePhaseTransmissibilityLaw::singlePhaseTransmissibility(element, fvGeometry, scvf, fluxVarsCache) | ||
| 157 | * aNw / aTot; | ||
| 158 | |||
| 159 | return result; | ||
| 160 | } | ||
| 161 | }; | ||
| 162 | |||
| 163 | template<class Scalar> | ||
| 164 | struct BakkeOren | ||
| 165 | { | ||
| 166 | using NonWettingPhaseCache = EmptyCache; | ||
| 167 | |||
| 168 | /*! | ||
| 169 | * \brief Returns the conductivity of a throat. | ||
| 170 | * | ||
| 171 | * See Bakke & Oren (1997), eq. 9 | ||
| 172 | */ | ||
| 173 | template<class Element, class FVElementGeometry, class FluxVariablesCache> | ||
| 174 | ✗ | static Scalar nonWettingPhaseTransmissibility(const Element& element, | |
| 175 | const FVElementGeometry& fvGeometry, | ||
| 176 | const typename FVElementGeometry::SubControlVolumeFace& scvf, | ||
| 177 | const FluxVariablesCache& fluxVarsCache) | ||
| 178 | { | ||
| 179 | // Tora et al. (2012), quite close for single-phase value of square | ||
| 180 | using std::sqrt; | ||
| 181 | ✗ | const Scalar throatLength = fluxVarsCache.throatLength(); | |
| 182 | ✗ | const auto nPhaseIdx = fluxVarsCache.nPhaseIdx(); | |
| 183 | ✗ | const Scalar aNw = fluxVarsCache.throatCrossSectionalArea(nPhaseIdx); | |
| 184 | ✗ | const Scalar rEff = 0.5*(sqrt(aNw / M_PI) + fluxVarsCache.throatInscribedRadius()); | |
| 185 | ✗ | const Scalar result = rEff*rEff*aNw / (8.0*throatLength); | |
| 186 | ✗ | return result; | |
| 187 | } | ||
| 188 | }; | ||
| 189 | } // end namespace NonWettingPhaseTransmissibility | ||
| 190 | } // end namespace Dumux::PoreNetwork | ||
| 191 | |||
| 192 | #endif | ||
| 193 |