GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/porenetwork/common/throatproperties.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 52 83 62.7%
Functions: 9 9 100.0%
Branches: 47 260 18.1%
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 PoreNetworkModels
10 * \brief This file contains functions related to calculate pore-throat properties.
11 */
12 #ifndef DUMUX_PNM_BASE_THROAT_PROPERTIES_HH
13 #define DUMUX_PNM_BASE_THROAT_PROPERTIES_HH
14
15 #include <string>
16 #include <cmath>
17 #include <numeric>
18
19 #include <dune/common/exceptions.hh>
20 #include <dune/common/reservedvector.hh>
21
22 namespace Dumux::PoreNetwork::Throat {
23
24 //! Collection of different pore-throat shapes
25 enum class Shape
26 { scaleneTriangle, equilateralTriangle, square, rectangle, circle, twoPlates, polygon };
27
28 //! Get the shape from a string description of the shape
29 74 inline std::string shapeToString(Shape s)
30 {
31
6/8
✓ Branch 0 taken 20 times.
✓ Branch 1 taken 21 times.
✓ Branch 2 taken 21 times.
✓ Branch 3 taken 5 times.
✓ Branch 4 taken 5 times.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 74 switch (s)
32 {
33 40 case Shape::scaleneTriangle: return "ScaleneTriangle";
34 42 case Shape::equilateralTriangle: return "EquilateralTriangle";
35 42 case Shape::square: return "Square";
36 10 case Shape::rectangle: return "Rectangle";
37 10 case Shape::circle: return "Circle";
38 4 case Shape::twoPlates: return "TwoPlates";
39 ✗ case Shape::polygon: return "Polygon";
40 ✗ default: DUNE_THROW(Dune::InvalidStateException, "Unknown shape!");
41 }
42 }
43
44 //! Get the shape from a string description of the shape
45 20 inline Shape shapeFromString(const std::string& s)
46 {
47
3/6
✗ Branch 1 not taken.
✓ Branch 2 taken 20 times.
✓ Branch 3 taken 20 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 20 times.
✗ Branch 6 not taken.
 40 if (s == shapeToString(Shape::equilateralTriangle)) return Shape::equilateralTriangle;
48
3/6
✗ Branch 1 not taken.
✓ Branch 2 taken 20 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 20 times.
✓ Branch 5 taken 20 times.
✗ Branch 6 not taken.
 20 else if (s == shapeToString(Shape::scaleneTriangle)) return Shape::scaleneTriangle;
49
5/6
✓ Branch 1 taken 19 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 20 times.
✓ Branch 5 taken 4 times.
✓ Branch 6 taken 16 times.
 39 else if (s == shapeToString(Shape::square)) return Shape::square;
50
4/6
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 4 times.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
 5 else if (s == shapeToString(Shape::rectangle)) return Shape::rectangle;
51
5/6
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 4 times.
✓ Branch 5 taken 1 times.
✓ Branch 6 taken 3 times.
 7 else if (s == shapeToString(Shape::circle)) return Shape::circle;
52
3/6
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
 2 else if (s == shapeToString(Shape::twoPlates)) return Shape::twoPlates;
53 ✗ else if (s == shapeToString(Shape::polygon)) return Shape::polygon;
54 ✗ else DUNE_THROW(Dune::InvalidStateException, s << " is not a valid shape");
55 }
56
57 /*!
58 * \brief Returns the radius of a pore throat
59 *
60 * \param poreRadiusOne The radius of the first pore
61 * \param poreRadiusTwo The radius of the second pore
62 * \param centerTocenterDist The center-to-center distance between the pores
63 * \param n Fitting parameter
64 *
65 * Joekar-Niasar et al. (2008) https://doi.org/10.1007/s11242-007-9191-7
66 */
67 template<class Scalar>
68 43923 inline Scalar averagedRadius(const Scalar poreRadiusOne, const Scalar poreRadiusTwo, const Scalar centerTocenterDist, const Scalar n = 0.1)
69 {
70
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 43923 times.
 43923 assert(n > 0.0);
71 using std::sin; using std::cos; using std::pow;
72 43923 const Scalar rOneTilde = poreRadiusOne/centerTocenterDist;
73 43923 const Scalar rTwoTilde = poreRadiusTwo/centerTocenterDist;
74 43923 const Scalar a = sin(M_PI/4.0);
75 43923 const Scalar b = cos(M_PI/4.0);
76 43923 const Scalar rhoOne = rOneTilde*a / pow((1.0 - rOneTilde*b), n);
77 43923 const Scalar rhoTwo = rTwoTilde*a / pow((1.0 - rTwoTilde*b), n);
78 43923 const Scalar rTilde = rhoOne*rhoTwo * pow((pow(rhoOne, 1.0/n) + pow(rhoTwo, 1.0/n)), -n);
79 43923 return rTilde * centerTocenterDist;
80 }
81
82
83 //! Returns the corner half angle
84 template<class Scalar>
85 959950 inline Dune::ReservedVector<Scalar, 4> cornerHalfAngles(Shape shape)
86 {
87
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 959950 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 959950 switch(shape)
88 {
89 ✗ case Shape::equilateralTriangle:
90 {
91 ✗ const Scalar value = M_PI/6.0;
92 ✗ return Dune::ReservedVector<Scalar, 4>{value, value, value};
93 }
94 959950 case Shape::square:
95 {
96 959950 const Scalar value = M_PI/4.0;
97 1919900 return Dune::ReservedVector<Scalar, 4>{value, value, value, value};
98 }
99 ✗ case Shape::rectangle:
100 {
101 ✗ const Scalar value = M_PI/4.0;
102 ✗ return Dune::ReservedVector<Scalar, 4>{value, value, value, value};
103 }
104 ✗ case Shape::circle:
105 {
106 ✗ const Scalar value = 0.5*M_PI; // we define the (single) corner angle of a circle as 180°
107 ✗ return { value };
108 }
109 ✗ case Shape::polygon: DUNE_THROW(Dune::InvalidStateException, "Corner half angles for polygons must be calculated explicitly");
110 ✗ default: DUNE_THROW(Dune::InvalidStateException, "Unknown shape");
111 // TODO implement angles for scaleneTriangle as given by Valvatne & Blunt (2004)
112 }
113 }
114
115 //! Returns the value of the shape factor for an equilateral triangle
116 template<class Scalar>
117 inline constexpr Scalar shapeFactorEquilateralTriangle() noexcept
118 {
119 using std::sqrt;
120 return sqrt(3.0)/36.0;
121 }
122
123 //! Returns the value of the shape factor for a square
124 template<class Scalar>
125 inline constexpr Scalar shapeFactorSquare() noexcept
126 {
127 return 1.0/16.0;
128 }
129
130 //! Returns the value of the shape factor for a rectangle
131 template<class Scalar>
132 inline constexpr Scalar shapeFactorRectangle(const Scalar inscribedRadius, const Scalar height) noexcept
133 {
134 ✗ const Scalar a = 2.0*inscribedRadius; // shorter side length
135 ✗ const Scalar b = height; // longer side length
136 ✗ const Scalar area = a*b;
137 ✗ const Scalar perimeter = 2.0*a + 2.0*b;
138 ✗ return area / (perimeter*perimeter);
139 }
140
141 //! Returns the value of the shape factor for a circle
142 template<class Scalar>
143 inline constexpr Scalar shapeFactorCircle() noexcept
144 {
145 return 1.0/(4.0*M_PI);
146 }
147
148 //! Returns the value of the shape factor for a given shape
149 template<class Scalar>
150 35774 inline Scalar shapeFactor(Shape shape, const Scalar inscribedRadius)
151 {
152
4/6
✓ Branch 0 taken 33210 times.
✓ Branch 1 taken 2562 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 35774 switch(shape)
153 {
154 case Shape::equilateralTriangle: return shapeFactorEquilateralTriangle<Scalar>();
155 33210 case Shape::square: return shapeFactorSquare<Scalar>();
156 2562 case Shape::circle: return shapeFactorCircle<Scalar>();
157 1 case Shape::twoPlates: return 0.0; // TODO is this a good idea?
158 ✗ case Shape::polygon: DUNE_THROW(Dune::InvalidStateException, "The shape factor for a polygon has to be defined by the input data");
159 ✗ default: DUNE_THROW(Dune::InvalidStateException, "Unknown shape");
160 }
161 }
162
163 //! Returns the shape for a given shape factor
164 template<class Scalar>
165 inline constexpr Shape shape(const Scalar shapeFactor) noexcept
166 {
167
2/2
✓ Branch 0 taken 4122 times.
✓ Branch 1 taken 2260 times.
 6382 if (shapeFactor < shapeFactorEquilateralTriangle<Scalar>())
168 return Shape::scaleneTriangle;
169
1/2
✓ Branch 0 taken 4122 times.
✗ Branch 1 not taken.
 4122 else if (shapeFactor <= shapeFactorEquilateralTriangle<Scalar>())
170 return Shape::equilateralTriangle;
171
2/2
✓ Branch 0 taken 3262 times.
✓ Branch 1 taken 860 times.
 4122 else if (shapeFactor < shapeFactorSquare<Scalar>())
172 return Shape::rectangle;
173
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3262 times.
 3262 else if (shapeFactor == shapeFactorSquare<Scalar>())
174 return Shape::square;
175 ✗ else if (shapeFactor == shapeFactorCircle<Scalar>())
176 return Shape::circle;
177 else
178 return Shape::polygon;
179 }
180
181 //! Returns if a shape is regular
182 959946 inline bool isRegularShape(Shape shape)
183 {
184
1/5
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 959946 times.
 959946 switch (shape)
185 {
186 case Shape::equilateralTriangle: return true;
187 case Shape::square: return true;
188 case Shape::circle: return true;
189 ✗ case Shape::rectangle: return false;
190 ✗ case Shape::scaleneTriangle: return false;
191 ✗ case Shape::polygon: DUNE_THROW(Dune::InvalidStateException, "Equality of Corner half angles for polygons must be determined explicitly");
192 ✗ default:
193 ✗ DUNE_THROW(Dune::InvalidStateException, "Unknown shape");
194 }
195 }
196
197 //! Returns the cross-sectional area of a given geometry
198 template<class Scalar>
199 35778 inline Scalar totalCrossSectionalArea(const Shape shape, const Scalar inscribedRadius)
200 {
201 using std::sqrt;
202
4/5
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 33211 times.
✓ Branch 2 taken 2563 times.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
 35778 switch(shape)
203 {
204 2 case Shape::equilateralTriangle: return 3.0*sqrt(3.0)*inscribedRadius*inscribedRadius;
205 33211 case Shape::square: return 4.0*inscribedRadius*inscribedRadius;
206 2563 case Shape::circle: return M_PI*inscribedRadius*inscribedRadius;
207 2 case Shape::twoPlates: return 2.0*inscribedRadius;
208 ✗ default : DUNE_THROW(Dune::InvalidStateException, "Unsupported geometry: " << shapeToString(shape));
209 }
210 }
211
212 //! Returns the cross-sectional area of a rectangle
213 template<class Scalar>
214 inline constexpr Scalar totalCrossSectionalAreaForRectangle(const Scalar inscribedRadius, const Scalar height) noexcept
215 {
216 ✗ return 2.0*inscribedRadius * height;
217 }
218
219 //! Returns the number of corners of a given geometry
220 92457 inline std::size_t numCorners(Shape shape)
221 {
222
1/2
✓ Branch 0 taken 92457 times.
✗ Branch 1 not taken.
 92457 switch(shape)
223 {
224 case Shape::scaleneTriangle: return 3;
225 case Shape::equilateralTriangle: return 3;
226 case Shape::square: return 4;
227 case Shape::rectangle: return 4;
228 case Shape::circle: return 0;
229 ✗ default : DUNE_THROW(Dune::InvalidStateException, "Unsupported geometry: " << shapeToString(shape));
230 }
231 }
232
233 /*!
234 * \brief Return the cross-sectional area of a wetting layer residing in a corner of a throat
235 *
236 * \param curvatureRadius The radius of curvature within the throat (surfaceTension/pc)
237 * \param contactAngle The contact angle within the throat
238 * \param cornerHalfAngle The corner half-angle of the the throat's corner of interest
239 *
240 * See, e.g, Valvatne & Blunt (2004), eq. A7 https://doi.org/10.1029/2003WR002627
241 */
242 template<class Scalar>
243 369828 inline constexpr Scalar wettingLayerCrossSectionalArea(const Scalar curvatureRadius,
244 const Scalar contactAngle,
245 const Scalar cornerHalfAngle) noexcept
246 {
247 using std::sin;
248 using std::cos;
249 369828 return curvatureRadius*curvatureRadius *(cos(contactAngle) * cos(contactAngle + cornerHalfAngle) / sin(cornerHalfAngle)
250 369828 + cornerHalfAngle + contactAngle - M_PI/2.0);
251 }
252
253 } // end Dumux::PoreNetwork::Throat
254
255 #endif
256