- Email: [email protected]
A stochastic approach for predicting tortuosity in porous media via pore network modeling
Computers and Geotechnics 120 (2020) 103406
Contents lists available at ScienceDirect
Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo
Research Paper
A stochastic approach for predicting tortuosity in porous media via pore network modeling
T
⁎
Ramin Nemati, Javad Rahbar Shahrouzi , Reza Alizadeh Faculty of Chemical Engineering, Sahand University of Technology, Sahand New Town, 53318-17634 Tabriz, Iran
A R T I C LE I N FO
A B S T R A C T
Keywords: Tortuosity Pore network modeling Stochastic Probability Porous media
In this study, a novel procedure is introduced and employed to calculate the tortuosity of porous media using pore network modeling. For this purpose, a 3D pore network model for porous media is constructed using network elements (pores and throats) along with its geometrical and topological characteristics. A unique method is introduced for connecting each pore to neighbors in an arbitrary direction with a variable pore coordination number. Next, fluid flow is simulated and the pressure distribution is achieved. Then, a stochastic procedure is applied to determine the tortuosity of the porous media. Since the probability of each path is calculated based on fluid flow rates, the stochastic approach makes it possible to calculate the tortuosity for effective paths, and there is no need for total path determination. The tortuosity of Berea sandstone, as a case study, is computed to be 2.07 using the described procedure, which is in good agreement with the experimental data. Finally, the impact of the repetition number and network parameters, including network size, pore coordination number, throat radius distribution, pore radius distribution, and pore/throat shape on the tortuosity calculation is investigated in detail.
1. Introduction The geometric and topological complexity of porous media, particularly natural ones, makes the flow deviate from straight conductive pathways to curved and chaotic ones. The important intrinsic property of porous media that offers a means of quantifying this deflection is called tortuosity. Its significance can be observed in various applications of engineering and science ranging from chemical and electrochemical science (catalyst [1], membrane [2,3], biofilm [4], fuel cell [5], and Li-ion battery [6,7]), soil science [8,9] and petroleum engineering [10] to acoustics applications [11,12]. The concept of tortuosity was first introduced in the theory of porous media by Carman in 1937 to compensate for more fundamental errors in the bundle of capillary tubes model [13–15]. The objectionable fundamental assumption of the simple capillary model was that the streamlines were considered completely straight and parallel to each other while a real porous medium supports flow in any direction [13,16]. Later, Carman used the tortuosity as an additional parameter to approximate the actual length of the paths that fluid travels [14]. However, it was still an unknown question of whether the tortuosity is an adjustable parameter or a function of some intrinsic properties of pore space [17]. To find the answer to this question, several attempts were made by some researchers, which in most cases several empirical ⁎
relations were suggested to relate tortuosity to the porosity of the porous media. Some of these relations were reviewed by Shen and Chen [18]. Despite providing a better prediction of tortuosity, one of the unknown aspects regarding these empirical relations is the fact that they contain adjustable parameters that their nature is not well addressed. Discovering the relationship between tortuosity and connectivity was a fundamental step towards understanding the nature of tortuosity. The connectivity, number of channels connected to a pore, is one of the well-known structural parameters that has a significant effect on the transport properties, particularly on tortuosity. This parameter, with various interpretations and definitions, has been investigated in some studies [19–22]. For example, Jarvis et al. [19] expressed that in the image-based methods, considering a different number of nearest neighbors for a voxel, in a cubic grid, results in the distinct definitions of connection and consequently can strongly affect the degree of connectivity and percolation thresholds. Cheon et al. [20], by comparing the obtained data for Berea and Otway sandstones, showed that a larger number of connecting paths results in higher permeability and affects tortuosity. In another work, Promentilla et al. [21], using micro-tomographic images, showed that better connection of pore space in segmented pores results in a decrease in tortuosity and vice versa. In their work, a two-pore voxel is considered to be connected when both of
Corresponding author. E-mail address: [email protected] (J. Rahbar Shahrouzi).
https://doi.org/10.1016/j.compgeo.2019.103406 Received 20 July 2019; Received in revised form 25 September 2019; Accepted 17 December 2019 0266-352X/ © 2019 Elsevier Ltd. All rights reserved.
Computers and Geotechnics 120 (2020) 103406
R. Nemati, et al.
Nomenclature A fk g gij g Pi gTij G l Lij L Pi LTij m N Ncoi NTout Nu pi Q qji qr qtot qv
r R R Pi t xi, yi, zi Z
Cross-sectional area of pore/throat perpendicular to the flow direction Frequency of each path Conductance Conductance between pore i and j Conductance of pore i Conductance of throat ij Shape factor of pore/throat Length of system Distance between pore i and j Half-length of pore i Length of throat ij Number of inlet pore Repetition number of algorithm Pore coordination number of pore i Number of throats that the fluid can flow out of Pc Number of unique paths Pressure of pore i Total volumetric flow of network Volumetric flow between pore i and j Volumetric flow of rth inlet pore Sum of the volumetric flow passing through the throats out of Pc Cumulative volumetric flow passing through the throats 1
Greek symbols λ μ τ τk τr
λ l
Average length of real flow path Viscosity of fluid Tortuosity of network Tortuosity of kth unique path Tortuosity of network for rth inlet pore
Abbreviations CFD DIC LBM PNM P Pc Pn T
Computational Fluid Dynamic Digital Image Correlation Lattice Boltzmann Model Pore Network Model Pore Current pore Next pore Throat
inexpensive and can be readily used for a wide range of porous media [32,33]. However, most of the numerical methods need experimental results as input data and also for verification of simulation results. This means that combining the experimental and numerical methods can be a powerful tool to estimate tortuosity and improve our understanding of fluid flow in porous media. Several numerical methods, such as Lattice Boltzmann Model (LBM) [16,17,33], Computational Fluid Dynamic (CFD) [34,35], and Pore Network Models (PNMs) are used to analyze the fluid flow in the porous media. Among these methods, LBM is more popular because of solving discrete Boltzmann transport equation instead of hard nonlinear equations (e.g., Navier-Stokes), and its ability to simulate the flow for a vast range of conditions and flow regimes [16]. For instance, Matyka et al. [16] used LBM to find a relationship between hydraulic tortuosity and porosity. They constructed 2D porous media with equal size square obstacles and then solved flow streamlines for laminar fluid flow, which resulted in a logarithmic relationship between tortuosity and porosity. In another study [33], the impact of aspect ratio on the tortuosity value of irregular porous media was investigated using a 2D stochastic model. Despite the increasing use of LBM, several limitations of this method cannot be ignored: discontinuity of streamlines, applicability only for two-dimensional models, and porosity dependency of presented correlations [14,24,32]. Among numerical methods, PNMs have many desirable features, such as providing a realistic model, having the ability to adapt to real pore space, using simpler equations, and having relatively low-computational costs that attract researchers. The mentioned advantages of PNMs have led to these methods being among the most popular methods to investigate fluid flow in porous media in recent decades. For instance, these methods are used in various applications of science and engineering, including investigation of multiphase flow, adsorption, biomass growth, hydrate dissociation process, and catalyst deactivation [36–40]. In spite of the aforementioned interesting features and broad applicability of PNMs, there is no evidence of using this model for fluxweighted average tortuosity calculation [33]. However, it appears that PNMs can be proper models for estimating tortuosity in combination with stochastic methods, which is proposed in this research. The aim of
them are orthogonally in face-to-face contact with each other. They also confirmed that the degree of connectivity and measured tortuosity depend on the digital resolution of images. Besides, Bellini et al. [22], by defining connectivity as a characteristic parameter related to the diffusional accessibility of the structure voids, demonstrated that the connectivity and tortuosity are inversely proportional to each other. Due to the unclear concept of tortuosity, because of its various methods of measurement, there is not a consistent scientific definition in the literature that uniquely describes it [14,23,24]; thus, in the literature, four different types of definition including geometric, hydraulic, electrical and diffusive tortuosity can be observed, which are often used interchangeably [14,16]. Although providing the consistent definition of tortuosity may be challenging, perhaps the most perceptible and simplistic description is that of the average ratio of the lengths of real flow paths to the length of the system as Eq. (1) [25].
τ=
to v Random number between 0 and 1 Inscribed radius of pore/throat Inscribed radius of pore i Index of selected throat that the fluid flow out of Pc Coordinates of the center of pore i Number of nearest neighboring pores
(1)
It should be noted that the calculation of the average length of flow paths can be carried out in two different approaches: in the first approach, different volumetric flow rates of fluid in each path are neglected, and only the actual length of flow paths is used to calculate the average length. However, in the second approach, by considering the volumetric flow rate, the weighted length is employed [14]. The later approach that is used in this study is more realistic when fluid flow in porous media is considered [25]. In the literature, four types of methods, including experimental [26], analytical [27,28], numerical [16,17,29], and imaging methods [30,31] have been mainly conducted to estimate the value of different types of tortuosity. Despite providing an acceptable estimation of tortuosity, experimental, analytical, and imaging methods have some limitations. For example, experimental methods are expensive and time-consuming, analytical methods are used for ideal porous media and they are not applicable to natural ones, and in the case of imagebased techniques, the 2D image analysis methods (2D DIC) are limited to in-plane processes, and the 3D X-ray tomography methods are expensive. On the other hand, numerical methods are relatively 2
Computers and Geotechnics 120 (2020) 103406
R. Nemati, et al.
coordination number of 4 in two dimensional or 6 in three dimensional); however, most of the realistic porous media are originally random in structure and such simple models could not generate proper macroscopic prediction models. To construct an irregular pore network, the center of pores is distributed randomly in three-dimensional space and a radius is randomly assigned to each pore, but the completely random distribution of pores causes the network to become uneven and also the probability of overlapping of pores to become greater. To overcome these problems, a gridded network is considered so that the center of every pore can only be selected randomly inside one of these grids. In such a situation, after distributing the pores in the specified spaces, overlapping may be happened again for some pores with large radii. To deal with this problem, each overlapped pore is reselected randomly inside the same gird until overlapping disappears. If there is still overlapping, the pore is replaced by one of the small pores. After generating the network and assigning geometric values for each pore and assuming pore length to be twice as the pore radius [41], the next and main step is to present a procedure to establish connectivity between pores. This step that is of key importance in topologically adapting the constructed model with original pore space has been studied widely. However, in this study, a new procedure is introduced so that more compatibility could be gained between the model and the original pore space. After assigning a random coordination number for each pore, the following procedure is applied to establish connectivity between pores:
this study is to introduce a novel approach that uses the stochastic method to calculate the pathway of incompressible fluid flow inside the pore network. In this approach, the stochastic method makes it possible to calculate tortuosity for effective paths, and there is no need for total path determination. The stochastic algorithm helps PNM to be qualified as an interesting and time-saving method for calculating tortuosity. Another novelty of this article is devising a method to assign the topological properties of the PNM, i.e., pore placement and pore connectivity, in a flexible way so that each pore can connect to other pores without any restrictions in terms of number and direction of the connection. Also, as the assumptions, it is considered to have a single-phase and steady-state fluid flow within the network, in which all the connecting paths are completely saturated with the fluid. In the present work, using the aforementioned stochastic method and flexibly constructed PNM, the hydraulic tortuosity of Berea sandstone is estimated and compared with experimental data. Furthermore, the effect of network size on the estimation of the tortuosity is discussed. Finally, the effect of some microscopic parameters of the network, including pore coordination number, throat radius distribution, pore radius distribution, and pore/throat shape are also investigated. 2. Modeling procedure In this section, the construction method of the network and calculation of the network tortuosity is presented. Firstly, the proposed method of constructing the three-dimensional, topologically equivalent, pore network is explained. In the second step, the pressure distribution of each pore is calculated using continuity and hydrodynamic equations, and finally, the tortuosity of the constructed network is calculated by the roulette wheel selection algorithm. Detailed representation of the steps above is declared below, respectively.
a) The distance between the pore (i) and other existing pores (throat length) is calculated using Eq. (2).
LTij =
2.1. Network construction
(x i − x j )2 + (yi − yj )2 + (z i − zj )2 − (R Pi + R Pj )
(2)
b) The set of distances is arranged in ascending order, and a determined number (Z) of the nearest neighboring pores is selected. c) Ncoi numbers of these pores are randomly selected to connect to the ith pore. If none of the selected pores have the capacity to connect, another pore is selected and substituted. d) Steps (a) to (c) are repeated until all the pores are connected to neighboring pores.
The constructed pore network in this study, similar to conventional ones, consists of two elements: pores and throats with the triangular, square, and circular cross-sectional shapes. Pores are the main elements of the network that are connected together with throats, and the number of throats connected to each pore is called the pore coordination number of that pore. In order to easily describe the whole network and for being the proposed algorithm convenient to handle, an n × n × n cell array is used in which each cell corresponds to a pore. In this way, for a typical pore, necessary input data (pore radius, length, shape, shape factor, coordination number, etc.) along with the values obtained from calculations (cross-sectional area, volume, and pore conductance, pore pressure, etc.) will be assigned. As well as the information of connected pores (pore number, position, radius, length, shape, conductance, pore pressure, etc.) and connected throats (throat length, radius, shape, cross-sectional area, volume, conductance, volumetric flow, etc.) will be stored in the information cell of same pore during the network construction and the flow simulation steps. By considering this approach, each pore, in addition to its own properties, contains information about its neighboring connected pores and connected throats that makes pores interactive and causes path-finding calculations to be easier. It should be mentioned that in most cases, the resistance to the flow within the pore is neglected as a simplifying assumption, but in this study, both of these elements are assumed to resist against the flow and to affect pressure distribution. Fig. 1 schematically shows two pores connected to each other with a throat. The topology of porous media that shows the spatial positioning and connectivity of pores, especially in tortuosity studies, is of great importance in terms of the compatibility of the constructed model. Two main topological properties are the structure of porous media (regular or irregular structure) as well as the mean coordination number and its distribution. Most of the models constructed to this point are considered regular (i.e., the pores positioned in the nodes of regular lattice with a
Fig. 2 shows the aforementioned algorithm for connecting the pores to each other in detail. The advantage of this procedure is that it allows each pore to connect to the arbitrary number of neighboring pores in an arbitrary direction, i.e., each pore can connect to other pores without any restrictions in terms of number and direction of the connection. The 2D schematic representation of the regular and irregular network constructed according to the mentioned steps is shown in Fig. 3, in which red, blue, and white pores represent the inlet, outlet, and inner pores of the network, respectively, and gray lines are the throats. In the next step, the cross-sectional shape and inscribed radius of throats are assigned using existing data [42] while their length is resulted as a function of the distance between centers of two pores minus
Fig. 1. The schematic representation of the connected pores via a throat. 3
Computers and Geotechnics 120 (2020) 103406
R. Nemati, et al.
Fig. 3. The schematic representation of (a) regular network and (b) irregular network.
2.2. Simulation of flow The most important aspect of the flow simulation inside the network is to calculate the distribution of pressure in the network. This can be done by applying appropriate boundary conditions and using the mass conservation and Hagen-Poiseuille equation. This operation is described below in more detail. The first step in obtaining the distribution of pressure in the network is to write a mass balance for each pore. Because of the fluid incompressibility, volume conservation can be used instead of mass conservation, which is shown in Eq. (3). NCoi
∑
Fig. 2. The algorithm for establishing the connection between the pores.
qji = 0 j = 1, 2, ...Ncoi
j=1
their inscribed radiuses according to Eq. (2). After assigning geometrical (inscribed radius and shape of pores and throats) and topological properties (mean coordination number and its distribution), reported by Piri and Blunt [26], the cross-section area and the volume of circular, square and triangular elements (pore or throat) can be computed. So the porosity of the constructed network can be achieved from the sum of the volume of pores and throats divided into the total volume of the network. To make the network as consistent as possible with the input data in terms of throat porosity and total throat number (established connections), the network should be constructed over and over again until the actual values and the simulated results converge.
(3)
Pictorial representation of volume conservation in the 2D constructed network, for a pore with four coordination number, is shown in Fig. 4. For laminar flow, the volumetric flow rate between two connected pores can be described by Eq. (4):
qij =
gij Lij
(pj − pi )
(4)
In which conductance of each circular, square, and triangular element is calculated as Eqs. (5)–(7), respectively [42]:
circular: gc =
4
0.5GA2 μ
(5)
Computers and Geotechnics 120 (2020) 103406
R. Nemati, et al.
Fig. 5. Representation of the possible paths for the fluid flowing out of a new pore.
Fig. 4. Volumetric flow conservation for each pore.
square: gs =
0.5623GA2 μ
triangular: gt =
2.3. Calculation of hydraulic tortuosity (6) The essential parameter that must be determined to compute tortuosity is λ. For this purpose, all possible paths that connect the inlet to the outlet of the network should be detected and contribute to tortuosity calculation. However, not all these paths are equal in terms of fluid flow, i.e., some of them that allow passing more volumetric flow have more contribution to the tortuosity value and vice versa. Thus, the contribution of each path must be calculated on the flow rate basis to obtain a weighted mean tortuous path. However, it should be considered that, according to the previous researches, finding all possible paths between two nodes on a network graph is challenging, and in most cases, it is impossible [45]. To overcome this issue, a stochastic simulation algorithm is used to find effective paths and calculate tortuosity. To achieve this goal, as shown in Fig. 5, suppose that the fluid enters the network through the inlet pore P1 and exits via one of the outlet pores. In the first movement of the fluid, it fills the pore P2 through the throat T12 , because the only possible way that the fluid can continue to flow is through this throat. After filling the pore P2 , the fluid can continue flowing through one of the throats T23 , T24 , or T25 to fill one of the pores P3, P4 , or P5 , respectively. Selecting the next throat is implemented using a new method inspired from the Gillespie algorithm, which is used to simulate the time behavior of a chemical reaction system. The algorithm has been extended to new applications by other researchers (e.g., stochastic reaction network generation and kinetic modeling of conversion processes) [46–48].
0.15R2A μ
(7)
In these equations, shape factor (G) is the dimensionless geometric parameter of the cross-section area perpendicular to the fluid flow direction, which is defined as a ratio of area to perimeter squared and R is the radius of the inscribed circle in the pore/throat [42]. The harmonic mean of the conductance of pores i and j and the throat between them weighed by their lengths parallel to the fluid flow direction is calculated according to Eq. (8) [42,43]:
Lij gij
=
LTij gTij
+
L Pj L Pi + g Pi g Pj
(8)
By assuming two different pressure levels in the inlet and outlet (often zero pressure in the outlet and constant pressure in the inlet) as a boundary condition and also solving a set of linear algebraic equations obtained from a combination of Eqs. (3)–(8) and applying them to each pore, the fluid pressure in all the pores and sequentially pressure distribution through the network can be achieved [44]. After determining the pressure distribution, volumetric flow through each pore and total volumetric flow through the network can be calculated. 5
Computers and Geotechnics 120 (2020) 103406
R. Nemati, et al.
Berea sandstone, is presented in Fig. 7 in which, the inlet and outlet pores are shown with black, inner pores with dark gray and throats with light gray. The constructed network is a 23 × 23 × 23 gridded network with a physical size of 3 × 3 × 3 mm3, in which the distance between two grid points along all three x, y, and z orientations is about 130 μm. The network contains pores with cross-section shapes of 95.5% triangular, 4.3% square, and 0.2% circular. As well as, the minimum, maximum, and average inscribed radius of the pores are 3.623 μm, 73.539 μm, and 19.16 μm, respectively. In this study, the length of each pore is assumed twice as the pore inscribed radius, so the 3D shape of pore are triangular prism, square prism, or cylinder. The throats are also prisms, in which the base shapes are 90.7% triangular, 7.6% square, and 1.7% circular. Besides, the minimum, maximum, and average inscribed radius of throats are 0.903 μm, 56.850 μm, and 10.970 μm, each in turn. Despite the size of the network considered in this work is 23 × 23 × 23; however, to make the figure less crowded, a 12 × 12 × 12 network is presented instead. For calculating the tortuosity of the network, the proper number of repetitions must be determined for the mentioned stochastic algorithm (N). As stated above, each time the algorithm is repeated, one of the paths inside the network is randomly determined, and its tortuosity is obtained. By continuing the repetition of this algorithm, newer paths are specified and added to the list of previous paths and finally, the mean value of tortuosity is calculated. With increasing the repetition of the algorithm (N), this mean value converges to a constant value. The smallest amount of N, in which the mean value starts to converge, is the proper number. For a better understanding of this procedure, see Fig. 8(a), which is considered for the determination of the appropriate repetition number of the algorithm for an inlet pore. According to this figure, where the tortuosity of each path and mean tortuosity value of the network are shown, in the low number of N, the value of the mean tortuosity is not constant; however, by increasing the repetition, tortuosity converges to a constant value. As can be observed, after almost 600 repetitions of the procedure, the mean tortuosity does not alter considerably, and this number can be selected as a proper repetition number of the algorithm for all the other inlet pores. For more certainty, a detailed description of the statistical properties (maximum, minimum, and mean) for certain repetitions of the algorithm is shown in Fig. 8(b). As demonstrated in the figure, with increasing N, the mean tortuosity converges to a constant value, but convergence of mean value occurs in a lower number of N. This means that a limited number of repetition of the algorithm is sufficient to calculate the tortuosity of the network and determination of all possible paths is not necessary. In addition, Fig. 9 presents the histogram of the
The stochastic algorithm used in this study is as follows: 1. The number of throats (NTout) that the fluid can flow out of the current pore (Pc) is counted. 2. The sum of the volumetric flow (qtot ) passing through the throats out of the current pore is calculated. 3. A random number (r) between 0 and 1 is generated. 4. The tth throat is selected if
t−1 ν=1
5. 6. 7.
8.
t
∑ qν < rqtot ⩽ ∑ qν in which ν=1
1 ≤ t ≤ NTout. Through the throat t, the fluid moves to the next connected pore (Pn ). If Pn is one of the outlet pores, the next step appears, otherwise Pc = Pn and it proceeds to step 1. The length of the path, through which the fluid passes from the inlet pore to one of the outlet pores is calculated and tortuosity of this path is computed by Eq. (1). The steps (1) to (7) are repeated for N times to find other paths or repetitive ones. The determination of N is explained in the Result and Discussion section.
After completing repetitions, the number of unique paths and the frequency of each path are counted, and the average tortuosity of the network by assuming that the fluid enters from the rth inlet pore is calculated as follows: Nu
∑ fk τk τr =
k=1 Nu
∑ fk k=1
(9)
In Fig. 6, two paths amongst the thousands of possible paths, starting from the same inlet pore, which are generated using the mentioned algorithm, are shown. It should be stated that this algorithm does not handle simultaneous flow paths outgoing from the specified inlet pore, i.e., in each repetition only one unique path is detected, not two or more paths at the same repetition. To compensate for this limitation, we repeat the process of path detecting N times. If all the above steps are repeated for the other remaining inlet pores, the total tortuosity of the network can be calculated easily by weighted averaging of all the inlet pores, as Eq. (10): m
∑ qr τr τ=
r=1
Q
(10)
3. Result and discussion In this section, firstly, the simulation results for Berea sandstone are presented and compared with the experimental data reported in Refs. [20,49]. Since the comprehensive microscopic and macroscopic characteristics of Berea sandstone are fairly available (see [42]), they are used as input data to construct the PNM in the present work. After the representation of the constructed model, the analysis of the tortuosity calculation considering the fluid entering from one inlet pore is described. Afterward, the same analysis is performed to calculate the tortuosity of the entire network considering all inlet pores. In the second step, the effect of network size on tortuosity is discussed, and finally, the effect of some microscopic parameters of pore network (pore coordination number, pore radius distribution, and throat radius distribution along with pore/throat shape) on tortuosity is investigated. 3.1. Tortuosity calculation result Fig. 6. Representation of two different paths starting from the same inlet pore percolating through the medium and reaching one of the outlet pores.
The 3D schematic representation of the constructed regular and irregular pore network model, according to the existing properties of 6
Computers and Geotechnics 120 (2020) 103406
R. Nemati, et al.
the calculated tortuosity range is between 1.53 and 3.2. However, the tortuosity amounts ranging from 1.7 to 2.6 are more probable. The most important point about the stochastic model is the fact that this method must be repeated several times for the entire steps of a process. In this way, by averaging the results of all records, the stochastic model gives approximately the same result as the equivalent deterministic model. So in order to satisfy this, as well as checking the repeatability of the results, all the steps from network construction to total tortuosity calculation of network are repeated 25 times, which resulted in the overall mean tortuosity of 2.07 (see Fig. 12). The obtained overall tortuosity is almost in good agreement with the reported experimental data by Cheon et al. [20], who provided the average tortuosity in three different directions to be between 2 and 2.2 using experimental data from microfocus X-ray computed tomography system. Takahashi et al. [49] also reported the mean value of 2.14 for the tortuosity of Berea sandstone in a specified direction. 3.2. The effect of network size on tortuosity Due to the fact that the network size may affect tortuosity, in this section, the tortuosity of cubic networks of n × n × n with varying n from 15 to 25 is calculated and the results are shown in Fig. 13. According to this figure, for networks with n < 21, some oscillation is observed in the result of tortuosity, but for larger n, the results are almost the same. This result assures that a network size of 23 × 23 × 23, used in this work, is sufficient for calculating tortuosity. 3.3. The effect of the microscopic parameters of the network on tortuosity In order to study the effect of the microscopic parameters on tortuosity, several networks with the same porosity but different microscopic properties are created and the effect of pore coordination number, pore radius, and throat radius are investigated. Since changing these parameters changes the porosity, to compensate for the decrease or increase of porosity, the radius of the pores or throats are adjusted with the proper ratio. As the first case, to determine the effect of the pore coordination number on tortuosity, connections between pores in the constructed network are randomly eliminated to obtain networks with fewer pore coordination numbers. By applying this method, the tortuosity corresponding to the mean coordination numbers ranging from 2.95 to 4.2 is calculated. As can be seen in Fig. 14, with decreasing the pore coordination number, the tortuosity increases. This can be interpreted that by reducing the connections between the pores, a number of paths are also removed and the fluid is forced to flow in more tortuous paths. The second parameter examined here is the distribution of the throat radius. An attempt is carried out to find whether tortuosity depends on throat size distribution or not. To find the answer, tortuosity is calculated for networks with both non-uniformly distributed (min = 0.903 μm , max = 56.85 μm , average = 10.970 μm ) and uniformly distributed (10.970 μm ) throat radii and their comparison is shown in Fig. 15. According to this figure, the tortuosity range of the two cases is significantly different, and in the system with uniform throat radii, the average tortuosity of the network is 24% less than that in the network with non-uniformly distributed throat radii. To interpret this difference, it can be said that in the non-uniformly distributed state, the throats with radii less than the mean value have a higher resistance against the fluid flow, i.e., according to Eqs. (5)–(7), the conductance of a pore or a throat approximately scales as the fourth power of the radius. The presence of such a narrow throat in the fluid path causes the fluid to prefer paths with lower resistance and results in longer fluid paths. However, in the uniformly distributed radii, the resistance of throats is almost the same, so shorter paths are prioritized for fluid flow and, consequently, tortuosity decreases. Similar to the throat, the effect of non-uniformly distributed pore radius (min = 3.623μm , max = 73.539 μm , average = 19.167 μm ) and
Fig. 7. Representation of the constructed 3D pore network for Berea sandstone: (a) regular and (b) irregular. Inlet and outlet pores are shown in black color. Inner pores and throats are represented in dark and light gray, respectively.
tortuosity distribution of a given inlet pore for the repetition number of 10,000 and 1,000,000. It can be deduced from this figure that the tortuosity distribution in both repetition numbers is approximately the same and that 10,000 repetitions are a good representation of the network tortuosity. This lower number of repetitions results in lower computational cost. For the other starting pores, N is considered to be 10,000, and thus the total repetition number is 529 × 10000, in which 529 (23 × 23) is the number of all the inlet pores. It should be noted that when the stochastic algorithm is applied, some flow paths are selected several times. Multiple selections of some paths indicate the fact that these paths are more conductive and have a greater contribution to the fluid flow (Eq. (9)); therefore, they have a greater impact on the flux-weighted average tortuosity of the network. Fig. 10 illustrates the frequency of unique paths for 1,000,000 repetitions and highlights the importance of some paths in fluid flow. It should be noted that the unique paths can be obtained by assigning a matrix for each path according to the pore number sequences. If two matrices are the same, the corresponding two paths are not unique and both of them are the same. After selecting the proper repetition number, tortuosity corresponding to the remaining inlet pores (528 pores) is calculated by the same procedure and the total mean tortuosity of the network is computed to be equal to 2.06 using Eq. (10). Moreover, the tortuosity distribution for Berea sandstone is shown in Fig. 11. As seen in this figure, 7
Computers and Geotechnics 120 (2020) 103406
R. Nemati, et al.
Fig. 8. The impact of repetition number on the tortuosity of the flow path starting from a certain inlet pore: (a) tortuosity of each unique path and mean tortuosity, (b) mean, maximum, and minimum values of tortuosity for a selected number of repetitions for path selection.
uniformly distributed pore radius (19.167 μm ) is investigated. As can be seen in Fig. 16 both the non-uniformly and uniformly distributed pore have almost the same tortuosity. It is because of the higher mean value of pore radius that is nearly twice the throat mean radius. Therefore, they are still the throats that have the highest resistance to the flow and control it, so either non-uniformly or uniformly distributed pore radii has no significant effect on tortuosity. As the last parameter, to understand the effect of throat shapes on tortuosity, four networks with the same uniformly distributed throat radii but different throat shapes are constructed; in the network no. 1–3 throat shapes are triangular, square, and circular, respectively, but in the network no. 4 shapes are 90.7% triangular, 7.6% square, and 1.7% circular, according to the Berea sandstone data. It is observed that the networks with the uniform throat shapes have roughly equal tortuosity. It is because of the fact that in the same value of the radius, the ratio of conductances is a constant value, so by changing the shape of the throats, the conductances of the throats are increased or decreased with a certain ratio, and consequently, the fluid maintains the same pattern of motion. However, in the case of non-uniformly distributed throat shapes (network no. 4), tortuosity is about 4.5% higher than the networks with uniformly distributed throat shapes (see Fig. 17). This can be interpreted that in the uniform networks, all the throats have the same shapes and consequently, the same conductance, so the fluid prefers shorter paths, but in the case of non-uniformly distributed throat shapes, not all the conductance of throats are the same, and the fluid prefers paths with higher conductance (It should be noted that in the same inscribed radius, the conductance of triangular shapes is greater than square or circular throats). In addition, the effect of pore shapes on the tortuosity is negligible because, as discussed later, the main contribution to the value of tortuosity is related to the throats. Fig. 9. The tortuosity distribution of the network for a certain inlet pore: (a) N = 10,000 and (b) N = 1,000,000.
8
Computers and Geotechnics 120 (2020) 103406
R. Nemati, et al.
Fig. 10. Frequency of each unique path for 1,000,000 repetitions.
Fig. 14. The effect of the pore coordination number on tortuosity. Fig. 11. The tortuosity distribution for Berea sandstone.
Fig. 12. Repeatability of total tortuosity calculation.
Fig. 15. The effect of throat radius distribution on tortuosity.
Fig. 13. The effect of network size (n × n × n) on tortuosity.
4. Summary and conclusion In this work, using pore network modeling, a new stochastic approach is presented to calculate the tortuosity of porous media. For this goal, after assigning the geometric properties, an attempt is made to connect the pores in a new way to have better compatibility with the
Fig. 16. The effect of pore radius distribution on tortuosity.
9
Computers and Geotechnics 120 (2020) 103406
R. Nemati, et al.
[3] Karanikola V, Corral AF, Jiang H, Sáez AE, Ela WP, Arnold RG. Effects of membrane structure and operational variables on membrane distillation performance. J Membr Sci 2017;524:87–96. [4] Zhang TC, Bishop PL. Evaluation of tortuosity factors and effective diffusivities in biofilms. Water Res 1994;28:2279–87. [5] Asinari P, Quaglia MC, von Spakovsky MR, Kasula BV. Direct numerical calculation of the kinematic tortuosity of reactive mixture flow in the anode layer of solid oxide fuel cells by the lattice Boltzmann method. J Power Sources 2007;170:359–75. [6] Thorat IV, Stephenson DE, Zacharias NA, Zaghib K, Harb JN, Wheeler DR. Quantifying tortuosity in porous Li-ion battery materials. J Power Sources 2009;188:592–600. [7] Ebner M, Wood V. Tool for tortuosity estimation in lithium ion battery porous electrodes. J Electrochem Soc 2015;162:A3064–70. [8] Moldrup P, Olesen T, Komatsu T, Schjønning P, Rolston D. Tortuosity, diffusivity, and permeability in the soil liquid and gaseous phases. Soil Sci Soc Am J 2001;65:613–23. [9] Hernández, OleMF T, Gironás J, Braud I, Suárez F, Muñoz JF. Assessment of evaporation and water fluxes in a column of dry saline soil subject to different water table levels. Hydrol Process 2014;28:3655–69. [10] Azar JH, Javaherian A, Pishvaie MR, Nabi-Bidhendi M. An approach to defining tortuosity and cementation factor in carbonate reservoir rocks. J Petrol Sci Eng 2008;60:125–31. [11] Allard JF, Castagnede B, Henry M, Lauriks W. Evaluation of tortuosity in acoustic porous materials saturated by air. Rev Sci Instrum 1994;65:754–5. [12] Johnson DL, Plona T, Scala C, Pasierb F, Kojima H. Tortuosity and acoustic slow waves. Phys Rev Lett 1982;49:1840. [13] Hunt A, Ewing R, Horton R. What’s wrong with soil physics? Soil Sci Soc Am J 2013;77:1877–87. [14] Ghanbarian B, Hunt AG, Ewing RP, Sahimi M. Tortuosity in porous media: a critical review. Soil Sci Soc Am J 2013;77:1461–77. [15] Carman PC. Fluid flow through granular beds. Trans Inst Chem Eng 1937;15:150–66. [16] Matyka M, Khalili A, Koza Z. Tortuosity-porosity relation in porous media flow. Phys Rev E 2008;78. 026306. [17] Duda A, Koza Z, Matyka M. Hydraulic tortuosity in arbitrary porous media flow. Phys Rev E 2011;84. 036319. [18] Shen L, Chen Z. Critical review of the impact of tortuosity on diffusion. Chem Eng Sci 2007;62:3748–55. [19] Jarvis N, Larsbo M, Koestel J. Connectivity and percolation of structural pore networks in a cultivated silt loam soil quantified by X-ray tomography. Geoderma 2017;287:71–9. [20] Cheon D-S, Takahashi M, Park E-S. Permeability change and geometrical information of void spaces in Berea sandstone and Otway sandstone. Geosyst Eng 2017;20:1–8. [21] Promentilla MAB, Sugiyama T, Hitomi T, Takeda N. Quantification of tortuosity in hardened cement pastes using synchrotron-based X-ray computed microtomography. Cem Concr Res 2009;39:548–57. [22] Bellini S, Azzato G, Grandinetti M, Stellato V, De Marco G, Sun Y, et al. A novel connectivity factor for morphological characterization of membranes and porous media: A simulation study on structures of mono-sized spherical particles. Appl Sci 2018;8:573. [23] Thauvin F, Mohanty K. Network modeling of non-Darcy flow through porous media. Transp Porous Media 1998;31:19–37. [24] Clennell MB. Tortuosity: a guide through the maze. Geol Soc, London, Spe Publ 1997;122:299–344. [25] Koponen A, Kataja M, Timonen J. Tortuous flow in porous media. Phys Rev E 1996;54:406. [26] Soukup K, Hejtmánek V, Čapek P, Stanczyk K, Šolcová O. Modeling of contaminant migration through porous media after underground coal gasification in shallow coal seam. Fuel Process Technol 2015;140:188–97. [27] Du Plessis JP, Masliyah JH. Flow through isotropic granular porous media. Transp Porous Media 1991;6:207–21. [28] Ahmadi MM, Mohammadi S, Hayati AN. Analytical derivation of tortuosity and permeability of monosized spheres: A volume averaging approach. Phys Rev E 2011;83. 026312. [29] Sun Z, Tang X, Cheng G. Numerical simulation for tortuosity of porous media. Microporous Mesoporous Mater 2013;173:37–42. [30] Naveed M, Hamamoto S, Kawamoto K, Sakaki T, Takahashi M, Komatsu T, et al. Correlating gas transport parameters and X-ray computed tomography measurements in porous media. Soil Sci 2013;178:60–8. [31] Shanti NO, Chan VW, Stock SR, De Carlo F, Thornton K, Faber KT. X-ray microcomputed tomography and tortuosity calculations of percolating pore networks. Acta Mater 2014;71:126–35. [32] Al-Raoush RI, Madhoun IT. TORT3D: A MATLAB code to compute geometric tortuosity from 3D images of unconsolidated porous media. Powder Technol 2017;320:99–107. [33] Khabbazi AE, Hinebaugh J, Bazylak A. Determining the impact of rectangular grain aspect ratio on tortuosity–porosity correlations of two-dimensional stochastically generated porous media. Sci Bull 2016;61:601–11. [34] Azzato G, De Marco G, Stellato V, Sun Y, Caravella A. Tortuosity and connectivity evaluation by CFD simulation for morphological characterization of membranes and catalytic structures. Case study: CaF2-like structure. Chem Eng Sci 2018. [35] Pawlowski S, Nayak N, Meireles M, Portugal C, Velizarov S, Crespo J. CFD modelling of flow patterns, tortuosity and residence time distribution in monolithic porous columns reconstructed from X-ray tomography data. Chem Eng J 2018. [36] Blunt MJ, Jackson MD, Piri M, Valvatne PH. Detailed physics, predictive
Fig. 17. The effect of throat shape on tortuosity.
original pore space; then, the pressure distribution of the network is obtained posing proper equations to the pores and throats. Finally, to compute tortuosity, the stochastic approach is used to achieve more probable paths for fluid flow. The results show that a limited number of paths are a good representative of all paths to calculate flux-weighted tortuosity, i.e., there is no need to find all possible paths for fluid flow inside the network. As a result, it is found that both 3D regular and irregular pore networks with well-determined microscopic and macroscopic properties are appropriate models for computing the tortuosity of the porous media. A Comparison of the result of the proposed model and Berea sandstone data shows the validity of this claim. Furthermore, in order to show the effect of geometric and topological microscopic characteristics of the network on the tortuosity, the influence of network size, pore coordination number, throat radius distribution, pore radius distribution, and pore/throat shape are investigated. It is shown that networks with small sizes lead to uncertain results, but as network size increases, the confidence in tortuosity estimation via pore network increases. To investigate the effect of pore coordination number on tortuosity, a number of connections between the pores are randomly eliminated to obtain lower coordination numbers. It is shown that by reducing the number of connections, due to the loss of some conductive paths, tortuosity increases. Regarding the throat radius, the results show that, because of the presence of narrow throats, the tortuosity of network with non-uniformly distributed throat radii is 24% higher than that of the uniformly distributed throat radii. However, it is concluded that the size distribution of pores has no significant effect on tortuosity and the average pore radius can be used instead. This is due to the fact that the impact of throat size, which is almost half of the pore size, is much greater than the pore size in fluid flow resistance. This means, neither uniformly or non-uniformly distributed pore radii alter tortuosity. Moreover, it is revealed that the tortuosity of network with various shapes of throats, because of non-uniformly distribution of throat conductances, is 4.5% higher than the network with uniform throat shapes. Nevertheless, the shape of pores has a negligible effect on tortuosity. As a final result, it is concluded that the stochastic approach, since there is no need to determine all possible paths, is preferable to the deterministic method to calculate the tortuosity of the network. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. References [1] Prachayawarakorn S, Mann R. Effects of pore assembly architecture on catalyst particle tortuosity and reaction effectiveness. Catal Today 2007;128:88–99. [2] Manickam SS, Gelb J, McCutcheon JR. Pore structure characterization of asymmetric membranes: Non-destructive characterization of porosity and tortuosity. J Membr Sci 2014;454:549–54.
10
Computers and Geotechnics 120 (2020) 103406
R. Nemati, et al.
[37]
[38]
[39]
[40]
[41]
[42]
2005;71. 026301. [43] Hughes RG, Blunt MJ. Network modeling of multiphase flow in fractures. Adv Water Resour 2001;24:409–21. [44] Bakke S, Øren P-E. 3-D pore-scale modelling of sandstones and flow simulations in the pore networks. SPE J 1997;2:136–49. [45] Migliore M, Martorana V, Sciortino F. An algorithm to find all paths between two nodes in a graph. J Comput Phys 1990;87:231–6. [46] Gillespie DT. Exact stochastic simulation of coupled chemical reactions. J Phys Chem 1977;81:2340–61. [47] Shahrouzi JR, Guillaume D, Rouchon P, Da Costa P. Stochastic simulation and single events kinetic modeling: application to olefin oligomerization. Ind Eng Chem Res 2008;47:4308–16. [48] de Oliveira LP, Verstraete JJ, Kolb M. Molecule-based kinetic modeling by Monte Carlo methods for heavy petroleum conversion. Sci China Chem 2013;56:1608–22. [49] Takahashi M, Kato M, Lin W, Sato M. Three-dimensional pore geometry and permeability anisotropy of Berea sandstone under hydrostatic pressure: connecting path and tortuosity data obtained by microfocus X-ray CT. Geol Soc, London, Eng Geol Spec Publ 2016;27:207–15.
capabilities and macroscopic consequences for pore-network models of multiphase flow. Adv Water Resour 2002;25:1069–89. Raoof A, Nick H, Hassanizadeh SM, Spiers C. PoreFlow: A complex pore-network model for simulation of reactive transport in variably saturated porous media. Comput Geosci 2013;61:160–74. Rosenzweig R, Furman A, Shavit U. A channel network model as a framework for characterizing variably saturated flow in biofilm-affected soils. Vadose Zone J 2013;12. Abdoli S, Shafiei S, Raoof A, Ebadi A, Jafarzadeh Y. Insight into heterogeneity effects in methane hydrate dissociation via pore-scale modeling. Transp Porous Media 2018:1–19. Ye G, Wang H, Duan X, Sui Z, Zhou X, Coppens MO, et al. Pore network modeling of catalyst deactivation by coking, from single site to particle, during propane dehydrogenation. AIChE J 2018. Jamshidi S, Boozarjomehry RB, Pishvaie MR. Application of GA in optimization of pore network models generated by multi-cellular growth algorithms. Adv Water Resour 2009;32:1543–53. Piri M, Blunt MJ. Three-dimensional mixed-wet random pore-scale network modeling of two-and three-phase flow in porous media. I. Model description. Phys Rev E
11
Contents lists available at ScienceDirect
Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo
Research Paper
A stochastic approach for predicting tortuosity in porous media via pore network modeling
T
⁎
Ramin Nemati, Javad Rahbar Shahrouzi , Reza Alizadeh Faculty of Chemical Engineering, Sahand University of Technology, Sahand New Town, 53318-17634 Tabriz, Iran
A R T I C LE I N FO
A B S T R A C T
Keywords: Tortuosity Pore network modeling Stochastic Probability Porous media
In this study, a novel procedure is introduced and employed to calculate the tortuosity of porous media using pore network modeling. For this purpose, a 3D pore network model for porous media is constructed using network elements (pores and throats) along with its geometrical and topological characteristics. A unique method is introduced for connecting each pore to neighbors in an arbitrary direction with a variable pore coordination number. Next, fluid flow is simulated and the pressure distribution is achieved. Then, a stochastic procedure is applied to determine the tortuosity of the porous media. Since the probability of each path is calculated based on fluid flow rates, the stochastic approach makes it possible to calculate the tortuosity for effective paths, and there is no need for total path determination. The tortuosity of Berea sandstone, as a case study, is computed to be 2.07 using the described procedure, which is in good agreement with the experimental data. Finally, the impact of the repetition number and network parameters, including network size, pore coordination number, throat radius distribution, pore radius distribution, and pore/throat shape on the tortuosity calculation is investigated in detail.
1. Introduction The geometric and topological complexity of porous media, particularly natural ones, makes the flow deviate from straight conductive pathways to curved and chaotic ones. The important intrinsic property of porous media that offers a means of quantifying this deflection is called tortuosity. Its significance can be observed in various applications of engineering and science ranging from chemical and electrochemical science (catalyst [1], membrane [2,3], biofilm [4], fuel cell [5], and Li-ion battery [6,7]), soil science [8,9] and petroleum engineering [10] to acoustics applications [11,12]. The concept of tortuosity was first introduced in the theory of porous media by Carman in 1937 to compensate for more fundamental errors in the bundle of capillary tubes model [13–15]. The objectionable fundamental assumption of the simple capillary model was that the streamlines were considered completely straight and parallel to each other while a real porous medium supports flow in any direction [13,16]. Later, Carman used the tortuosity as an additional parameter to approximate the actual length of the paths that fluid travels [14]. However, it was still an unknown question of whether the tortuosity is an adjustable parameter or a function of some intrinsic properties of pore space [17]. To find the answer to this question, several attempts were made by some researchers, which in most cases several empirical ⁎
relations were suggested to relate tortuosity to the porosity of the porous media. Some of these relations were reviewed by Shen and Chen [18]. Despite providing a better prediction of tortuosity, one of the unknown aspects regarding these empirical relations is the fact that they contain adjustable parameters that their nature is not well addressed. Discovering the relationship between tortuosity and connectivity was a fundamental step towards understanding the nature of tortuosity. The connectivity, number of channels connected to a pore, is one of the well-known structural parameters that has a significant effect on the transport properties, particularly on tortuosity. This parameter, with various interpretations and definitions, has been investigated in some studies [19–22]. For example, Jarvis et al. [19] expressed that in the image-based methods, considering a different number of nearest neighbors for a voxel, in a cubic grid, results in the distinct definitions of connection and consequently can strongly affect the degree of connectivity and percolation thresholds. Cheon et al. [20], by comparing the obtained data for Berea and Otway sandstones, showed that a larger number of connecting paths results in higher permeability and affects tortuosity. In another work, Promentilla et al. [21], using micro-tomographic images, showed that better connection of pore space in segmented pores results in a decrease in tortuosity and vice versa. In their work, a two-pore voxel is considered to be connected when both of
Corresponding author. E-mail address: [email protected] (J. Rahbar Shahrouzi).
https://doi.org/10.1016/j.compgeo.2019.103406 Received 20 July 2019; Received in revised form 25 September 2019; Accepted 17 December 2019 0266-352X/ © 2019 Elsevier Ltd. All rights reserved.
Computers and Geotechnics 120 (2020) 103406
R. Nemati, et al.
Nomenclature A fk g gij g Pi gTij G l Lij L Pi LTij m N Ncoi NTout Nu pi Q qji qr qtot qv
r R R Pi t xi, yi, zi Z
Cross-sectional area of pore/throat perpendicular to the flow direction Frequency of each path Conductance Conductance between pore i and j Conductance of pore i Conductance of throat ij Shape factor of pore/throat Length of system Distance between pore i and j Half-length of pore i Length of throat ij Number of inlet pore Repetition number of algorithm Pore coordination number of pore i Number of throats that the fluid can flow out of Pc Number of unique paths Pressure of pore i Total volumetric flow of network Volumetric flow between pore i and j Volumetric flow of rth inlet pore Sum of the volumetric flow passing through the throats out of Pc Cumulative volumetric flow passing through the throats 1
Greek symbols λ μ τ τk τr
λ l
Average length of real flow path Viscosity of fluid Tortuosity of network Tortuosity of kth unique path Tortuosity of network for rth inlet pore
Abbreviations CFD DIC LBM PNM P Pc Pn T
Computational Fluid Dynamic Digital Image Correlation Lattice Boltzmann Model Pore Network Model Pore Current pore Next pore Throat
inexpensive and can be readily used for a wide range of porous media [32,33]. However, most of the numerical methods need experimental results as input data and also for verification of simulation results. This means that combining the experimental and numerical methods can be a powerful tool to estimate tortuosity and improve our understanding of fluid flow in porous media. Several numerical methods, such as Lattice Boltzmann Model (LBM) [16,17,33], Computational Fluid Dynamic (CFD) [34,35], and Pore Network Models (PNMs) are used to analyze the fluid flow in the porous media. Among these methods, LBM is more popular because of solving discrete Boltzmann transport equation instead of hard nonlinear equations (e.g., Navier-Stokes), and its ability to simulate the flow for a vast range of conditions and flow regimes [16]. For instance, Matyka et al. [16] used LBM to find a relationship between hydraulic tortuosity and porosity. They constructed 2D porous media with equal size square obstacles and then solved flow streamlines for laminar fluid flow, which resulted in a logarithmic relationship between tortuosity and porosity. In another study [33], the impact of aspect ratio on the tortuosity value of irregular porous media was investigated using a 2D stochastic model. Despite the increasing use of LBM, several limitations of this method cannot be ignored: discontinuity of streamlines, applicability only for two-dimensional models, and porosity dependency of presented correlations [14,24,32]. Among numerical methods, PNMs have many desirable features, such as providing a realistic model, having the ability to adapt to real pore space, using simpler equations, and having relatively low-computational costs that attract researchers. The mentioned advantages of PNMs have led to these methods being among the most popular methods to investigate fluid flow in porous media in recent decades. For instance, these methods are used in various applications of science and engineering, including investigation of multiphase flow, adsorption, biomass growth, hydrate dissociation process, and catalyst deactivation [36–40]. In spite of the aforementioned interesting features and broad applicability of PNMs, there is no evidence of using this model for fluxweighted average tortuosity calculation [33]. However, it appears that PNMs can be proper models for estimating tortuosity in combination with stochastic methods, which is proposed in this research. The aim of
them are orthogonally in face-to-face contact with each other. They also confirmed that the degree of connectivity and measured tortuosity depend on the digital resolution of images. Besides, Bellini et al. [22], by defining connectivity as a characteristic parameter related to the diffusional accessibility of the structure voids, demonstrated that the connectivity and tortuosity are inversely proportional to each other. Due to the unclear concept of tortuosity, because of its various methods of measurement, there is not a consistent scientific definition in the literature that uniquely describes it [14,23,24]; thus, in the literature, four different types of definition including geometric, hydraulic, electrical and diffusive tortuosity can be observed, which are often used interchangeably [14,16]. Although providing the consistent definition of tortuosity may be challenging, perhaps the most perceptible and simplistic description is that of the average ratio of the lengths of real flow paths to the length of the system as Eq. (1) [25].
τ=
to v Random number between 0 and 1 Inscribed radius of pore/throat Inscribed radius of pore i Index of selected throat that the fluid flow out of Pc Coordinates of the center of pore i Number of nearest neighboring pores
(1)
It should be noted that the calculation of the average length of flow paths can be carried out in two different approaches: in the first approach, different volumetric flow rates of fluid in each path are neglected, and only the actual length of flow paths is used to calculate the average length. However, in the second approach, by considering the volumetric flow rate, the weighted length is employed [14]. The later approach that is used in this study is more realistic when fluid flow in porous media is considered [25]. In the literature, four types of methods, including experimental [26], analytical [27,28], numerical [16,17,29], and imaging methods [30,31] have been mainly conducted to estimate the value of different types of tortuosity. Despite providing an acceptable estimation of tortuosity, experimental, analytical, and imaging methods have some limitations. For example, experimental methods are expensive and time-consuming, analytical methods are used for ideal porous media and they are not applicable to natural ones, and in the case of imagebased techniques, the 2D image analysis methods (2D DIC) are limited to in-plane processes, and the 3D X-ray tomography methods are expensive. On the other hand, numerical methods are relatively 2
Computers and Geotechnics 120 (2020) 103406
R. Nemati, et al.
coordination number of 4 in two dimensional or 6 in three dimensional); however, most of the realistic porous media are originally random in structure and such simple models could not generate proper macroscopic prediction models. To construct an irregular pore network, the center of pores is distributed randomly in three-dimensional space and a radius is randomly assigned to each pore, but the completely random distribution of pores causes the network to become uneven and also the probability of overlapping of pores to become greater. To overcome these problems, a gridded network is considered so that the center of every pore can only be selected randomly inside one of these grids. In such a situation, after distributing the pores in the specified spaces, overlapping may be happened again for some pores with large radii. To deal with this problem, each overlapped pore is reselected randomly inside the same gird until overlapping disappears. If there is still overlapping, the pore is replaced by one of the small pores. After generating the network and assigning geometric values for each pore and assuming pore length to be twice as the pore radius [41], the next and main step is to present a procedure to establish connectivity between pores. This step that is of key importance in topologically adapting the constructed model with original pore space has been studied widely. However, in this study, a new procedure is introduced so that more compatibility could be gained between the model and the original pore space. After assigning a random coordination number for each pore, the following procedure is applied to establish connectivity between pores:
this study is to introduce a novel approach that uses the stochastic method to calculate the pathway of incompressible fluid flow inside the pore network. In this approach, the stochastic method makes it possible to calculate tortuosity for effective paths, and there is no need for total path determination. The stochastic algorithm helps PNM to be qualified as an interesting and time-saving method for calculating tortuosity. Another novelty of this article is devising a method to assign the topological properties of the PNM, i.e., pore placement and pore connectivity, in a flexible way so that each pore can connect to other pores without any restrictions in terms of number and direction of the connection. Also, as the assumptions, it is considered to have a single-phase and steady-state fluid flow within the network, in which all the connecting paths are completely saturated with the fluid. In the present work, using the aforementioned stochastic method and flexibly constructed PNM, the hydraulic tortuosity of Berea sandstone is estimated and compared with experimental data. Furthermore, the effect of network size on the estimation of the tortuosity is discussed. Finally, the effect of some microscopic parameters of the network, including pore coordination number, throat radius distribution, pore radius distribution, and pore/throat shape are also investigated. 2. Modeling procedure In this section, the construction method of the network and calculation of the network tortuosity is presented. Firstly, the proposed method of constructing the three-dimensional, topologically equivalent, pore network is explained. In the second step, the pressure distribution of each pore is calculated using continuity and hydrodynamic equations, and finally, the tortuosity of the constructed network is calculated by the roulette wheel selection algorithm. Detailed representation of the steps above is declared below, respectively.
a) The distance between the pore (i) and other existing pores (throat length) is calculated using Eq. (2).
LTij =
2.1. Network construction
(x i − x j )2 + (yi − yj )2 + (z i − zj )2 − (R Pi + R Pj )
(2)
b) The set of distances is arranged in ascending order, and a determined number (Z) of the nearest neighboring pores is selected. c) Ncoi numbers of these pores are randomly selected to connect to the ith pore. If none of the selected pores have the capacity to connect, another pore is selected and substituted. d) Steps (a) to (c) are repeated until all the pores are connected to neighboring pores.
The constructed pore network in this study, similar to conventional ones, consists of two elements: pores and throats with the triangular, square, and circular cross-sectional shapes. Pores are the main elements of the network that are connected together with throats, and the number of throats connected to each pore is called the pore coordination number of that pore. In order to easily describe the whole network and for being the proposed algorithm convenient to handle, an n × n × n cell array is used in which each cell corresponds to a pore. In this way, for a typical pore, necessary input data (pore radius, length, shape, shape factor, coordination number, etc.) along with the values obtained from calculations (cross-sectional area, volume, and pore conductance, pore pressure, etc.) will be assigned. As well as the information of connected pores (pore number, position, radius, length, shape, conductance, pore pressure, etc.) and connected throats (throat length, radius, shape, cross-sectional area, volume, conductance, volumetric flow, etc.) will be stored in the information cell of same pore during the network construction and the flow simulation steps. By considering this approach, each pore, in addition to its own properties, contains information about its neighboring connected pores and connected throats that makes pores interactive and causes path-finding calculations to be easier. It should be mentioned that in most cases, the resistance to the flow within the pore is neglected as a simplifying assumption, but in this study, both of these elements are assumed to resist against the flow and to affect pressure distribution. Fig. 1 schematically shows two pores connected to each other with a throat. The topology of porous media that shows the spatial positioning and connectivity of pores, especially in tortuosity studies, is of great importance in terms of the compatibility of the constructed model. Two main topological properties are the structure of porous media (regular or irregular structure) as well as the mean coordination number and its distribution. Most of the models constructed to this point are considered regular (i.e., the pores positioned in the nodes of regular lattice with a
Fig. 2 shows the aforementioned algorithm for connecting the pores to each other in detail. The advantage of this procedure is that it allows each pore to connect to the arbitrary number of neighboring pores in an arbitrary direction, i.e., each pore can connect to other pores without any restrictions in terms of number and direction of the connection. The 2D schematic representation of the regular and irregular network constructed according to the mentioned steps is shown in Fig. 3, in which red, blue, and white pores represent the inlet, outlet, and inner pores of the network, respectively, and gray lines are the throats. In the next step, the cross-sectional shape and inscribed radius of throats are assigned using existing data [42] while their length is resulted as a function of the distance between centers of two pores minus
Fig. 1. The schematic representation of the connected pores via a throat. 3
Computers and Geotechnics 120 (2020) 103406
R. Nemati, et al.
Fig. 3. The schematic representation of (a) regular network and (b) irregular network.
2.2. Simulation of flow The most important aspect of the flow simulation inside the network is to calculate the distribution of pressure in the network. This can be done by applying appropriate boundary conditions and using the mass conservation and Hagen-Poiseuille equation. This operation is described below in more detail. The first step in obtaining the distribution of pressure in the network is to write a mass balance for each pore. Because of the fluid incompressibility, volume conservation can be used instead of mass conservation, which is shown in Eq. (3). NCoi
∑
Fig. 2. The algorithm for establishing the connection between the pores.
qji = 0 j = 1, 2, ...Ncoi
j=1
their inscribed radiuses according to Eq. (2). After assigning geometrical (inscribed radius and shape of pores and throats) and topological properties (mean coordination number and its distribution), reported by Piri and Blunt [26], the cross-section area and the volume of circular, square and triangular elements (pore or throat) can be computed. So the porosity of the constructed network can be achieved from the sum of the volume of pores and throats divided into the total volume of the network. To make the network as consistent as possible with the input data in terms of throat porosity and total throat number (established connections), the network should be constructed over and over again until the actual values and the simulated results converge.
(3)
Pictorial representation of volume conservation in the 2D constructed network, for a pore with four coordination number, is shown in Fig. 4. For laminar flow, the volumetric flow rate between two connected pores can be described by Eq. (4):
qij =
gij Lij
(pj − pi )
(4)
In which conductance of each circular, square, and triangular element is calculated as Eqs. (5)–(7), respectively [42]:
circular: gc =
4
0.5GA2 μ
(5)
Computers and Geotechnics 120 (2020) 103406
R. Nemati, et al.
Fig. 5. Representation of the possible paths for the fluid flowing out of a new pore.
Fig. 4. Volumetric flow conservation for each pore.
square: gs =
0.5623GA2 μ
triangular: gt =
2.3. Calculation of hydraulic tortuosity (6) The essential parameter that must be determined to compute tortuosity is λ. For this purpose, all possible paths that connect the inlet to the outlet of the network should be detected and contribute to tortuosity calculation. However, not all these paths are equal in terms of fluid flow, i.e., some of them that allow passing more volumetric flow have more contribution to the tortuosity value and vice versa. Thus, the contribution of each path must be calculated on the flow rate basis to obtain a weighted mean tortuous path. However, it should be considered that, according to the previous researches, finding all possible paths between two nodes on a network graph is challenging, and in most cases, it is impossible [45]. To overcome this issue, a stochastic simulation algorithm is used to find effective paths and calculate tortuosity. To achieve this goal, as shown in Fig. 5, suppose that the fluid enters the network through the inlet pore P1 and exits via one of the outlet pores. In the first movement of the fluid, it fills the pore P2 through the throat T12 , because the only possible way that the fluid can continue to flow is through this throat. After filling the pore P2 , the fluid can continue flowing through one of the throats T23 , T24 , or T25 to fill one of the pores P3, P4 , or P5 , respectively. Selecting the next throat is implemented using a new method inspired from the Gillespie algorithm, which is used to simulate the time behavior of a chemical reaction system. The algorithm has been extended to new applications by other researchers (e.g., stochastic reaction network generation and kinetic modeling of conversion processes) [46–48].
0.15R2A μ
(7)
In these equations, shape factor (G) is the dimensionless geometric parameter of the cross-section area perpendicular to the fluid flow direction, which is defined as a ratio of area to perimeter squared and R is the radius of the inscribed circle in the pore/throat [42]. The harmonic mean of the conductance of pores i and j and the throat between them weighed by their lengths parallel to the fluid flow direction is calculated according to Eq. (8) [42,43]:
Lij gij
=
LTij gTij
+
L Pj L Pi + g Pi g Pj
(8)
By assuming two different pressure levels in the inlet and outlet (often zero pressure in the outlet and constant pressure in the inlet) as a boundary condition and also solving a set of linear algebraic equations obtained from a combination of Eqs. (3)–(8) and applying them to each pore, the fluid pressure in all the pores and sequentially pressure distribution through the network can be achieved [44]. After determining the pressure distribution, volumetric flow through each pore and total volumetric flow through the network can be calculated. 5
Computers and Geotechnics 120 (2020) 103406
R. Nemati, et al.
Berea sandstone, is presented in Fig. 7 in which, the inlet and outlet pores are shown with black, inner pores with dark gray and throats with light gray. The constructed network is a 23 × 23 × 23 gridded network with a physical size of 3 × 3 × 3 mm3, in which the distance between two grid points along all three x, y, and z orientations is about 130 μm. The network contains pores with cross-section shapes of 95.5% triangular, 4.3% square, and 0.2% circular. As well as, the minimum, maximum, and average inscribed radius of the pores are 3.623 μm, 73.539 μm, and 19.16 μm, respectively. In this study, the length of each pore is assumed twice as the pore inscribed radius, so the 3D shape of pore are triangular prism, square prism, or cylinder. The throats are also prisms, in which the base shapes are 90.7% triangular, 7.6% square, and 1.7% circular. Besides, the minimum, maximum, and average inscribed radius of throats are 0.903 μm, 56.850 μm, and 10.970 μm, each in turn. Despite the size of the network considered in this work is 23 × 23 × 23; however, to make the figure less crowded, a 12 × 12 × 12 network is presented instead. For calculating the tortuosity of the network, the proper number of repetitions must be determined for the mentioned stochastic algorithm (N). As stated above, each time the algorithm is repeated, one of the paths inside the network is randomly determined, and its tortuosity is obtained. By continuing the repetition of this algorithm, newer paths are specified and added to the list of previous paths and finally, the mean value of tortuosity is calculated. With increasing the repetition of the algorithm (N), this mean value converges to a constant value. The smallest amount of N, in which the mean value starts to converge, is the proper number. For a better understanding of this procedure, see Fig. 8(a), which is considered for the determination of the appropriate repetition number of the algorithm for an inlet pore. According to this figure, where the tortuosity of each path and mean tortuosity value of the network are shown, in the low number of N, the value of the mean tortuosity is not constant; however, by increasing the repetition, tortuosity converges to a constant value. As can be observed, after almost 600 repetitions of the procedure, the mean tortuosity does not alter considerably, and this number can be selected as a proper repetition number of the algorithm for all the other inlet pores. For more certainty, a detailed description of the statistical properties (maximum, minimum, and mean) for certain repetitions of the algorithm is shown in Fig. 8(b). As demonstrated in the figure, with increasing N, the mean tortuosity converges to a constant value, but convergence of mean value occurs in a lower number of N. This means that a limited number of repetition of the algorithm is sufficient to calculate the tortuosity of the network and determination of all possible paths is not necessary. In addition, Fig. 9 presents the histogram of the
The stochastic algorithm used in this study is as follows: 1. The number of throats (NTout) that the fluid can flow out of the current pore (Pc) is counted. 2. The sum of the volumetric flow (qtot ) passing through the throats out of the current pore is calculated. 3. A random number (r) between 0 and 1 is generated. 4. The tth throat is selected if
t−1 ν=1
5. 6. 7.
8.
t
∑ qν < rqtot ⩽ ∑ qν in which ν=1
1 ≤ t ≤ NTout. Through the throat t, the fluid moves to the next connected pore (Pn ). If Pn is one of the outlet pores, the next step appears, otherwise Pc = Pn and it proceeds to step 1. The length of the path, through which the fluid passes from the inlet pore to one of the outlet pores is calculated and tortuosity of this path is computed by Eq. (1). The steps (1) to (7) are repeated for N times to find other paths or repetitive ones. The determination of N is explained in the Result and Discussion section.
After completing repetitions, the number of unique paths and the frequency of each path are counted, and the average tortuosity of the network by assuming that the fluid enters from the rth inlet pore is calculated as follows: Nu
∑ fk τk τr =
k=1 Nu
∑ fk k=1
(9)
In Fig. 6, two paths amongst the thousands of possible paths, starting from the same inlet pore, which are generated using the mentioned algorithm, are shown. It should be stated that this algorithm does not handle simultaneous flow paths outgoing from the specified inlet pore, i.e., in each repetition only one unique path is detected, not two or more paths at the same repetition. To compensate for this limitation, we repeat the process of path detecting N times. If all the above steps are repeated for the other remaining inlet pores, the total tortuosity of the network can be calculated easily by weighted averaging of all the inlet pores, as Eq. (10): m
∑ qr τr τ=
r=1
Q
(10)
3. Result and discussion In this section, firstly, the simulation results for Berea sandstone are presented and compared with the experimental data reported in Refs. [20,49]. Since the comprehensive microscopic and macroscopic characteristics of Berea sandstone are fairly available (see [42]), they are used as input data to construct the PNM in the present work. After the representation of the constructed model, the analysis of the tortuosity calculation considering the fluid entering from one inlet pore is described. Afterward, the same analysis is performed to calculate the tortuosity of the entire network considering all inlet pores. In the second step, the effect of network size on tortuosity is discussed, and finally, the effect of some microscopic parameters of pore network (pore coordination number, pore radius distribution, and throat radius distribution along with pore/throat shape) on tortuosity is investigated. 3.1. Tortuosity calculation result Fig. 6. Representation of two different paths starting from the same inlet pore percolating through the medium and reaching one of the outlet pores.
The 3D schematic representation of the constructed regular and irregular pore network model, according to the existing properties of 6
Computers and Geotechnics 120 (2020) 103406
R. Nemati, et al.
the calculated tortuosity range is between 1.53 and 3.2. However, the tortuosity amounts ranging from 1.7 to 2.6 are more probable. The most important point about the stochastic model is the fact that this method must be repeated several times for the entire steps of a process. In this way, by averaging the results of all records, the stochastic model gives approximately the same result as the equivalent deterministic model. So in order to satisfy this, as well as checking the repeatability of the results, all the steps from network construction to total tortuosity calculation of network are repeated 25 times, which resulted in the overall mean tortuosity of 2.07 (see Fig. 12). The obtained overall tortuosity is almost in good agreement with the reported experimental data by Cheon et al. [20], who provided the average tortuosity in three different directions to be between 2 and 2.2 using experimental data from microfocus X-ray computed tomography system. Takahashi et al. [49] also reported the mean value of 2.14 for the tortuosity of Berea sandstone in a specified direction. 3.2. The effect of network size on tortuosity Due to the fact that the network size may affect tortuosity, in this section, the tortuosity of cubic networks of n × n × n with varying n from 15 to 25 is calculated and the results are shown in Fig. 13. According to this figure, for networks with n < 21, some oscillation is observed in the result of tortuosity, but for larger n, the results are almost the same. This result assures that a network size of 23 × 23 × 23, used in this work, is sufficient for calculating tortuosity. 3.3. The effect of the microscopic parameters of the network on tortuosity In order to study the effect of the microscopic parameters on tortuosity, several networks with the same porosity but different microscopic properties are created and the effect of pore coordination number, pore radius, and throat radius are investigated. Since changing these parameters changes the porosity, to compensate for the decrease or increase of porosity, the radius of the pores or throats are adjusted with the proper ratio. As the first case, to determine the effect of the pore coordination number on tortuosity, connections between pores in the constructed network are randomly eliminated to obtain networks with fewer pore coordination numbers. By applying this method, the tortuosity corresponding to the mean coordination numbers ranging from 2.95 to 4.2 is calculated. As can be seen in Fig. 14, with decreasing the pore coordination number, the tortuosity increases. This can be interpreted that by reducing the connections between the pores, a number of paths are also removed and the fluid is forced to flow in more tortuous paths. The second parameter examined here is the distribution of the throat radius. An attempt is carried out to find whether tortuosity depends on throat size distribution or not. To find the answer, tortuosity is calculated for networks with both non-uniformly distributed (min = 0.903 μm , max = 56.85 μm , average = 10.970 μm ) and uniformly distributed (10.970 μm ) throat radii and their comparison is shown in Fig. 15. According to this figure, the tortuosity range of the two cases is significantly different, and in the system with uniform throat radii, the average tortuosity of the network is 24% less than that in the network with non-uniformly distributed throat radii. To interpret this difference, it can be said that in the non-uniformly distributed state, the throats with radii less than the mean value have a higher resistance against the fluid flow, i.e., according to Eqs. (5)–(7), the conductance of a pore or a throat approximately scales as the fourth power of the radius. The presence of such a narrow throat in the fluid path causes the fluid to prefer paths with lower resistance and results in longer fluid paths. However, in the uniformly distributed radii, the resistance of throats is almost the same, so shorter paths are prioritized for fluid flow and, consequently, tortuosity decreases. Similar to the throat, the effect of non-uniformly distributed pore radius (min = 3.623μm , max = 73.539 μm , average = 19.167 μm ) and
Fig. 7. Representation of the constructed 3D pore network for Berea sandstone: (a) regular and (b) irregular. Inlet and outlet pores are shown in black color. Inner pores and throats are represented in dark and light gray, respectively.
tortuosity distribution of a given inlet pore for the repetition number of 10,000 and 1,000,000. It can be deduced from this figure that the tortuosity distribution in both repetition numbers is approximately the same and that 10,000 repetitions are a good representation of the network tortuosity. This lower number of repetitions results in lower computational cost. For the other starting pores, N is considered to be 10,000, and thus the total repetition number is 529 × 10000, in which 529 (23 × 23) is the number of all the inlet pores. It should be noted that when the stochastic algorithm is applied, some flow paths are selected several times. Multiple selections of some paths indicate the fact that these paths are more conductive and have a greater contribution to the fluid flow (Eq. (9)); therefore, they have a greater impact on the flux-weighted average tortuosity of the network. Fig. 10 illustrates the frequency of unique paths for 1,000,000 repetitions and highlights the importance of some paths in fluid flow. It should be noted that the unique paths can be obtained by assigning a matrix for each path according to the pore number sequences. If two matrices are the same, the corresponding two paths are not unique and both of them are the same. After selecting the proper repetition number, tortuosity corresponding to the remaining inlet pores (528 pores) is calculated by the same procedure and the total mean tortuosity of the network is computed to be equal to 2.06 using Eq. (10). Moreover, the tortuosity distribution for Berea sandstone is shown in Fig. 11. As seen in this figure, 7
Computers and Geotechnics 120 (2020) 103406
R. Nemati, et al.
Fig. 8. The impact of repetition number on the tortuosity of the flow path starting from a certain inlet pore: (a) tortuosity of each unique path and mean tortuosity, (b) mean, maximum, and minimum values of tortuosity for a selected number of repetitions for path selection.
uniformly distributed pore radius (19.167 μm ) is investigated. As can be seen in Fig. 16 both the non-uniformly and uniformly distributed pore have almost the same tortuosity. It is because of the higher mean value of pore radius that is nearly twice the throat mean radius. Therefore, they are still the throats that have the highest resistance to the flow and control it, so either non-uniformly or uniformly distributed pore radii has no significant effect on tortuosity. As the last parameter, to understand the effect of throat shapes on tortuosity, four networks with the same uniformly distributed throat radii but different throat shapes are constructed; in the network no. 1–3 throat shapes are triangular, square, and circular, respectively, but in the network no. 4 shapes are 90.7% triangular, 7.6% square, and 1.7% circular, according to the Berea sandstone data. It is observed that the networks with the uniform throat shapes have roughly equal tortuosity. It is because of the fact that in the same value of the radius, the ratio of conductances is a constant value, so by changing the shape of the throats, the conductances of the throats are increased or decreased with a certain ratio, and consequently, the fluid maintains the same pattern of motion. However, in the case of non-uniformly distributed throat shapes (network no. 4), tortuosity is about 4.5% higher than the networks with uniformly distributed throat shapes (see Fig. 17). This can be interpreted that in the uniform networks, all the throats have the same shapes and consequently, the same conductance, so the fluid prefers shorter paths, but in the case of non-uniformly distributed throat shapes, not all the conductance of throats are the same, and the fluid prefers paths with higher conductance (It should be noted that in the same inscribed radius, the conductance of triangular shapes is greater than square or circular throats). In addition, the effect of pore shapes on the tortuosity is negligible because, as discussed later, the main contribution to the value of tortuosity is related to the throats. Fig. 9. The tortuosity distribution of the network for a certain inlet pore: (a) N = 10,000 and (b) N = 1,000,000.
8
Computers and Geotechnics 120 (2020) 103406
R. Nemati, et al.
Fig. 10. Frequency of each unique path for 1,000,000 repetitions.
Fig. 14. The effect of the pore coordination number on tortuosity. Fig. 11. The tortuosity distribution for Berea sandstone.
Fig. 12. Repeatability of total tortuosity calculation.
Fig. 15. The effect of throat radius distribution on tortuosity.
Fig. 13. The effect of network size (n × n × n) on tortuosity.
4. Summary and conclusion In this work, using pore network modeling, a new stochastic approach is presented to calculate the tortuosity of porous media. For this goal, after assigning the geometric properties, an attempt is made to connect the pores in a new way to have better compatibility with the
Fig. 16. The effect of pore radius distribution on tortuosity.
9
Computers and Geotechnics 120 (2020) 103406
R. Nemati, et al.
[3] Karanikola V, Corral AF, Jiang H, Sáez AE, Ela WP, Arnold RG. Effects of membrane structure and operational variables on membrane distillation performance. J Membr Sci 2017;524:87–96. [4] Zhang TC, Bishop PL. Evaluation of tortuosity factors and effective diffusivities in biofilms. Water Res 1994;28:2279–87. [5] Asinari P, Quaglia MC, von Spakovsky MR, Kasula BV. Direct numerical calculation of the kinematic tortuosity of reactive mixture flow in the anode layer of solid oxide fuel cells by the lattice Boltzmann method. J Power Sources 2007;170:359–75. [6] Thorat IV, Stephenson DE, Zacharias NA, Zaghib K, Harb JN, Wheeler DR. Quantifying tortuosity in porous Li-ion battery materials. J Power Sources 2009;188:592–600. [7] Ebner M, Wood V. Tool for tortuosity estimation in lithium ion battery porous electrodes. J Electrochem Soc 2015;162:A3064–70. [8] Moldrup P, Olesen T, Komatsu T, Schjønning P, Rolston D. Tortuosity, diffusivity, and permeability in the soil liquid and gaseous phases. Soil Sci Soc Am J 2001;65:613–23. [9] Hernández, OleMF T, Gironás J, Braud I, Suárez F, Muñoz JF. Assessment of evaporation and water fluxes in a column of dry saline soil subject to different water table levels. Hydrol Process 2014;28:3655–69. [10] Azar JH, Javaherian A, Pishvaie MR, Nabi-Bidhendi M. An approach to defining tortuosity and cementation factor in carbonate reservoir rocks. J Petrol Sci Eng 2008;60:125–31. [11] Allard JF, Castagnede B, Henry M, Lauriks W. Evaluation of tortuosity in acoustic porous materials saturated by air. Rev Sci Instrum 1994;65:754–5. [12] Johnson DL, Plona T, Scala C, Pasierb F, Kojima H. Tortuosity and acoustic slow waves. Phys Rev Lett 1982;49:1840. [13] Hunt A, Ewing R, Horton R. What’s wrong with soil physics? Soil Sci Soc Am J 2013;77:1877–87. [14] Ghanbarian B, Hunt AG, Ewing RP, Sahimi M. Tortuosity in porous media: a critical review. Soil Sci Soc Am J 2013;77:1461–77. [15] Carman PC. Fluid flow through granular beds. Trans Inst Chem Eng 1937;15:150–66. [16] Matyka M, Khalili A, Koza Z. Tortuosity-porosity relation in porous media flow. Phys Rev E 2008;78. 026306. [17] Duda A, Koza Z, Matyka M. Hydraulic tortuosity in arbitrary porous media flow. Phys Rev E 2011;84. 036319. [18] Shen L, Chen Z. Critical review of the impact of tortuosity on diffusion. Chem Eng Sci 2007;62:3748–55. [19] Jarvis N, Larsbo M, Koestel J. Connectivity and percolation of structural pore networks in a cultivated silt loam soil quantified by X-ray tomography. Geoderma 2017;287:71–9. [20] Cheon D-S, Takahashi M, Park E-S. Permeability change and geometrical information of void spaces in Berea sandstone and Otway sandstone. Geosyst Eng 2017;20:1–8. [21] Promentilla MAB, Sugiyama T, Hitomi T, Takeda N. Quantification of tortuosity in hardened cement pastes using synchrotron-based X-ray computed microtomography. Cem Concr Res 2009;39:548–57. [22] Bellini S, Azzato G, Grandinetti M, Stellato V, De Marco G, Sun Y, et al. A novel connectivity factor for morphological characterization of membranes and porous media: A simulation study on structures of mono-sized spherical particles. Appl Sci 2018;8:573. [23] Thauvin F, Mohanty K. Network modeling of non-Darcy flow through porous media. Transp Porous Media 1998;31:19–37. [24] Clennell MB. Tortuosity: a guide through the maze. Geol Soc, London, Spe Publ 1997;122:299–344. [25] Koponen A, Kataja M, Timonen J. Tortuous flow in porous media. Phys Rev E 1996;54:406. [26] Soukup K, Hejtmánek V, Čapek P, Stanczyk K, Šolcová O. Modeling of contaminant migration through porous media after underground coal gasification in shallow coal seam. Fuel Process Technol 2015;140:188–97. [27] Du Plessis JP, Masliyah JH. Flow through isotropic granular porous media. Transp Porous Media 1991;6:207–21. [28] Ahmadi MM, Mohammadi S, Hayati AN. Analytical derivation of tortuosity and permeability of monosized spheres: A volume averaging approach. Phys Rev E 2011;83. 026312. [29] Sun Z, Tang X, Cheng G. Numerical simulation for tortuosity of porous media. Microporous Mesoporous Mater 2013;173:37–42. [30] Naveed M, Hamamoto S, Kawamoto K, Sakaki T, Takahashi M, Komatsu T, et al. Correlating gas transport parameters and X-ray computed tomography measurements in porous media. Soil Sci 2013;178:60–8. [31] Shanti NO, Chan VW, Stock SR, De Carlo F, Thornton K, Faber KT. X-ray microcomputed tomography and tortuosity calculations of percolating pore networks. Acta Mater 2014;71:126–35. [32] Al-Raoush RI, Madhoun IT. TORT3D: A MATLAB code to compute geometric tortuosity from 3D images of unconsolidated porous media. Powder Technol 2017;320:99–107. [33] Khabbazi AE, Hinebaugh J, Bazylak A. Determining the impact of rectangular grain aspect ratio on tortuosity–porosity correlations of two-dimensional stochastically generated porous media. Sci Bull 2016;61:601–11. [34] Azzato G, De Marco G, Stellato V, Sun Y, Caravella A. Tortuosity and connectivity evaluation by CFD simulation for morphological characterization of membranes and catalytic structures. Case study: CaF2-like structure. Chem Eng Sci 2018. [35] Pawlowski S, Nayak N, Meireles M, Portugal C, Velizarov S, Crespo J. CFD modelling of flow patterns, tortuosity and residence time distribution in monolithic porous columns reconstructed from X-ray tomography data. Chem Eng J 2018. [36] Blunt MJ, Jackson MD, Piri M, Valvatne PH. Detailed physics, predictive
Fig. 17. The effect of throat shape on tortuosity.
original pore space; then, the pressure distribution of the network is obtained posing proper equations to the pores and throats. Finally, to compute tortuosity, the stochastic approach is used to achieve more probable paths for fluid flow. The results show that a limited number of paths are a good representative of all paths to calculate flux-weighted tortuosity, i.e., there is no need to find all possible paths for fluid flow inside the network. As a result, it is found that both 3D regular and irregular pore networks with well-determined microscopic and macroscopic properties are appropriate models for computing the tortuosity of the porous media. A Comparison of the result of the proposed model and Berea sandstone data shows the validity of this claim. Furthermore, in order to show the effect of geometric and topological microscopic characteristics of the network on the tortuosity, the influence of network size, pore coordination number, throat radius distribution, pore radius distribution, and pore/throat shape are investigated. It is shown that networks with small sizes lead to uncertain results, but as network size increases, the confidence in tortuosity estimation via pore network increases. To investigate the effect of pore coordination number on tortuosity, a number of connections between the pores are randomly eliminated to obtain lower coordination numbers. It is shown that by reducing the number of connections, due to the loss of some conductive paths, tortuosity increases. Regarding the throat radius, the results show that, because of the presence of narrow throats, the tortuosity of network with non-uniformly distributed throat radii is 24% higher than that of the uniformly distributed throat radii. However, it is concluded that the size distribution of pores has no significant effect on tortuosity and the average pore radius can be used instead. This is due to the fact that the impact of throat size, which is almost half of the pore size, is much greater than the pore size in fluid flow resistance. This means, neither uniformly or non-uniformly distributed pore radii alter tortuosity. Moreover, it is revealed that the tortuosity of network with various shapes of throats, because of non-uniformly distribution of throat conductances, is 4.5% higher than the network with uniform throat shapes. Nevertheless, the shape of pores has a negligible effect on tortuosity. As a final result, it is concluded that the stochastic approach, since there is no need to determine all possible paths, is preferable to the deterministic method to calculate the tortuosity of the network. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. References [1] Prachayawarakorn S, Mann R. Effects of pore assembly architecture on catalyst particle tortuosity and reaction effectiveness. Catal Today 2007;128:88–99. [2] Manickam SS, Gelb J, McCutcheon JR. Pore structure characterization of asymmetric membranes: Non-destructive characterization of porosity and tortuosity. J Membr Sci 2014;454:549–54.
10
Computers and Geotechnics 120 (2020) 103406
R. Nemati, et al.
[37]
[38]
[39]
[40]
[41]
[42]
2005;71. 026301. [43] Hughes RG, Blunt MJ. Network modeling of multiphase flow in fractures. Adv Water Resour 2001;24:409–21. [44] Bakke S, Øren P-E. 3-D pore-scale modelling of sandstones and flow simulations in the pore networks. SPE J 1997;2:136–49. [45] Migliore M, Martorana V, Sciortino F. An algorithm to find all paths between two nodes in a graph. J Comput Phys 1990;87:231–6. [46] Gillespie DT. Exact stochastic simulation of coupled chemical reactions. J Phys Chem 1977;81:2340–61. [47] Shahrouzi JR, Guillaume D, Rouchon P, Da Costa P. Stochastic simulation and single events kinetic modeling: application to olefin oligomerization. Ind Eng Chem Res 2008;47:4308–16. [48] de Oliveira LP, Verstraete JJ, Kolb M. Molecule-based kinetic modeling by Monte Carlo methods for heavy petroleum conversion. Sci China Chem 2013;56:1608–22. [49] Takahashi M, Kato M, Lin W, Sato M. Three-dimensional pore geometry and permeability anisotropy of Berea sandstone under hydrostatic pressure: connecting path and tortuosity data obtained by microfocus X-ray CT. Geol Soc, London, Eng Geol Spec Publ 2016;27:207–15.
capabilities and macroscopic consequences for pore-network models of multiphase flow. Adv Water Resour 2002;25:1069–89. Raoof A, Nick H, Hassanizadeh SM, Spiers C. PoreFlow: A complex pore-network model for simulation of reactive transport in variably saturated porous media. Comput Geosci 2013;61:160–74. Rosenzweig R, Furman A, Shavit U. A channel network model as a framework for characterizing variably saturated flow in biofilm-affected soils. Vadose Zone J 2013;12. Abdoli S, Shafiei S, Raoof A, Ebadi A, Jafarzadeh Y. Insight into heterogeneity effects in methane hydrate dissociation via pore-scale modeling. Transp Porous Media 2018:1–19. Ye G, Wang H, Duan X, Sui Z, Zhou X, Coppens MO, et al. Pore network modeling of catalyst deactivation by coking, from single site to particle, during propane dehydrogenation. AIChE J 2018. Jamshidi S, Boozarjomehry RB, Pishvaie MR. Application of GA in optimization of pore network models generated by multi-cellular growth algorithms. Adv Water Resour 2009;32:1543–53. Piri M, Blunt MJ. Three-dimensional mixed-wet random pore-scale network modeling of two-and three-phase flow in porous media. I. Model description. Phys Rev E
11