The Path Tracking Method as an alternative for tortuosity determination in granular beds

Granular Matter, Sep 2018

Tortuosity is one of the most elusive parameters of porous media. The fundamental question is whether it may be computed from the geometry only or the transport equations must be solved first. Elimination of the transport equations would significantly decrease the computation time and memory consumption and thus allow investigating larger samples. We compare the geometric to hydraulic tortuosity of a sphere-packed porous media. We applied the Discrete Element Method to generate a set of virtual beds based on experimental data taking into account the real porosity and particle distribution, the Lattice Boltzmann Method to compute the hydraulic tortuosity and geometrical approach, i.e. so-called Path Tracking Method, to calculate the geometrical tortuosity. Our study shows that the calculation time can be reduced from hours (if the LBM is used) to seconds (if the PTM is applied) without losing the accuracy of the final results. The relative error between average values of the tortuosity obtained for both used methods is less than 3%. We show that the applied geometrical method may serve as an attractive alternative to hydraulic tortuosity, particularly in large granular systems.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

The Path Tracking Method as an alternative for tortuosity determination in granular beds

Granular Matter November 2018, 20:72 | Cite as The Path Tracking Method as an alternative for tortuosity determination in granular beds AuthorsAuthors and affiliations Wojciech SobieskiMaciej MatykaJarosław GołembiewskiSeweryn Lipiński Open Access Original Paper First Online: 29 September 2018 Abstract Tortuosity is one of the most elusive parameters of porous media. The fundamental question is whether it may be computed from the geometry only or the transport equations must be solved first. Elimination of the transport equations would significantly decrease the computation time and memory consumption and thus allow investigating larger samples. We compare the geometric to hydraulic tortuosity of a sphere-packed porous media. We applied the Discrete Element Method to generate a set of virtual beds based on experimental data taking into account the real porosity and particle distribution, the Lattice Boltzmann Method to compute the hydraulic tortuosity and geometrical approach, i.e. so-called Path Tracking Method, to calculate the geometrical tortuosity. Our study shows that the calculation time can be reduced from hours (if the LBM is used) to seconds (if the PTM is applied) without losing the accuracy of the final results. The relative error between average values of the tortuosity obtained for both used methods is less than 3%. We show that the applied geometrical method may serve as an attractive alternative to hydraulic tortuosity, particularly in large granular systems. KeywordsTortuosity Granular beds Path Tracking Method Lattice Boltzmann Method Discrete Element Method  1 Introduction The fluid flow through porous media is a complex, multi-scale phenomena. Fluid particles penetrate pore space and draw paths of various shapes and lengths. If those paths are short, the drag of porous matrix is relatively low. If particles turn and travel on long and wiggly trajectories, the drag is high, what results in low permeability of the media. In order to take this effect into account, Kozeny [24] introduced a new physical quantity named tortuosity. He applied this parameter to correct the value of the hydraulic drop (which in turn depends on the path length) occurring during fluid flow through a porous body (Fig. 1a). The same parameter was then used by Carman [9] who used the tortuosity as a correction factor of the flow velocity in a porous medium (Fig. 1b). Additionally, Carman proposed a new formula for calculation of the so-called specific surface. The Kozeny formula modified by the Carman is now widely known as the Kozeny–Carman equation, which is destined to predict the pressure drop in porous media: $$\begin{aligned} \frac{dp}{dx} = C_{KC} T_f S^2_{0,Carman} \frac{(1-\phi )^2}{\phi ^3} (\mu v_f), \end{aligned}$$ (1) where p—pressure [Pa], x—coordinate along which the pressure drop occurs [m], \(C_{KC}\)—Kozeny–Carman pore shape factor (a model constant), which should be equal to 5.0 (no unit) [9], \(T_f\)—tortuosity factor (no unit) defined as the square of the tortuosity T (no unit), \(\phi \)—porosity (no unit), \(\mu \)—dynamic viscosity coefficient [kg/(m\(\cdot \)s)], \(v_f\)—filtration velocity [m/s]. Open image in new window Fig. 1 Visualisation of ideas proposed by: a Kozeny—the real path length inside pore channels (\(L_p\)) is higher than the thickness of the porous medium (\(L_0\)); b Carman—the real fluid velocity inside pore channels (\(v_p\)) is higher that the space averaged filtration velocity (\(v_f\)) Tortuosity T is defined as the ratio of the actual path length inside channels (\(L_p\)) to the thickness of the porous medium (\(L_0\)): $$\begin{aligned} T =\frac{L_p}{L_0}. \end{aligned}$$ (2) In a general case, tortuosity may be understood as a geometrical quantity (geometric tortuosity, \(T^g\)) or as a flow property (hydraulic tortuosity, \(T^h\)). Other kinds of tortuosity are also known, e.g. diffusional or electric tortuosity. Important is that obtaining the tortuosity value in a specific case is still very difficult and for this reason the application of the Kozeny–Carman equation (as well as many other formulas where tortuosity occurs) is limited. Due to such a state, all studies related to new ways of obtaining this physical quantity are needed. It should be mentioned that in the literature analytical formulas used to calculate the tortuosity in granular beds based on porosity may be found. However, the range of results obtained from these formulas is quite large [44], what leads to difficulties in chosing appropriate formula in a specific case. The hydraulic tortuosity may be calculated from a velocity field as follows: $$\begin{aligned} T^h = \frac{\sum v}{\sum v_x}, \end{aligned}$$ (3) where \(v_x\) [m/s] is the velocity component in direction of macroscopic flow in the porous material and v [m/s] is the velocity magnitude. The sum goes over pore space [13, 22]. The velocity field needed for the application of the formula (3) may be obtained experimentally (e.g. by the use of the Particle Image Velocimetry) or numerically using Computational Fluid Dynamics. In the last case, Finite Difference Method, Finite Volume Method, Immersed Boundary Method, Lattice Boltzmann Method or even Finite Element Method may be applied. In the context of porous media, Lattice Boltzmann Method is particularly convenient due it the feature of easy mapping of complex geometries [29, 30, 31, 38]. The methodology of calculating the hydraulic tortuosity on the basis of a velocity field was used before. It used directly velocity field [22] and path lines [32]. There, the hydraulic tortuosity was measured in 2D random generated pore structures consisting of squares and rectangles. Similar investigations were then performed by Nabovati and Sousa [35] and Duda et al. [13]. The results obtained by these Authors are quite consistent, but cannot be directly compared with results of this article due to differences in the porous media model. The same methodology was used by Aaltosalmi [3]. His work was focused on Lattice Boltzmann Method as such and on its application as for the simulation of fluid flows through materials such as paper, fabric, sandstone and others. Lane [26] performed a study on the hydraulic and the electric tortuosity in 2D sphere packs with a normal distribution. Papers directly related to the investigations of the hydraulic tortuosity in granular porous media are rare; e.g. Wang [49] investigated the relation between the porosity and the tortuosity in a channel filled with spherical particles with a constant diameter. The number of particles was not stated, but it can be estimated at a few thousands. In these simulations (performed by the use of the Discrete Element Method) the particles were randomly added over the bed in a calculation loop and next they fall down by gravity. In the second stage, the BGK D3Q19 Lattice Boltzmann Model was applied. On the basis of figures shown there it may by stated that for the porosity equal to 0.4 (which value corresponds approximately to the porosity of a bed sample used in this article), the hydraulic tortuosity is approx. 1.1975. Data on computational time needed for performing simulations is not given. Summarizing, papers related to the investigations of the hydraulic tortuosity in large porous beds created with the use of the Discrete Element Method, taking into account the particle size distribution, are not available. In the article, we propose such a methodology, however the relation between the particle size distribution and the tortuosity is not studied. The geometric tortuosity may be obtained in three ways: (1) through correlations with other parameters characterizing the geometry of a porous body; (2) experimentally and (3) by using algorithms that track the direction or shape of the porous space. In the first approach it is usually assumed that the geometric tortuosity is a direct function of the porosity. It is the most popular approach to calculate the tortuosity, and many functions (empirical or derived analytically), in which the relationships between these quantities are proposed, may be found in literature [2, 4, 21, 44, 48]. In the experimental investigations, the computed tomography and image analysis [15, 18, 50], acoustic methods [19, 20, 27], optical methods [36] as well as other methods [17] are used. In case of the last mentioned possibility, firstly the geometry of the pore part must be defined (or determined), and then its arrangement in the space is tracked. Yu and Li [52] analyzed a two-dimensional system of connected squares and proposed an analytical formula to calculate the tortuosity. Starley et al. [46] propose a three-dimensional tracer metric numerical model, based on mutually perpendicular cylinders, to estimate the tortuosity. Montes et al. [34] used the so-called Pathfinding A* algorithm to analyse electrical and thermal tortuosity in powder compacts (however, details concerning this algorithm were not given). In the literature, several other algoritms refering to the analysis the geometric tortuosity (or in fact: any path in a partially limited space) may be also found, however they are related to other research areas, such as medicine, geodesy, games theory, logistics, environmental science and other, and co they cannot be applied or even compared with the investigations described in the article. The numerical algorithm used here to calculate the geometric tortuosity in granular beds is called the Path Tracking Method (see Sect. 2.5) and was developed by Sobieski in 2009 [42, 45] (details in Sect. 2.5). In the paper, two different methods of determining the tortuosity with the use of the same virtual beds are compared—these are the methodology proposed by Koponen et al. [22, 23] and a methodology based on numerical analysis of the shapes of pore channels. Our investigations contain following steps: measurement of particle diameters and porosity in a sample of a real granular bed; specifying the distribution of diameters in this bed; division of the bed into fractions according to the specified distribution; generating the virtual bed representing the real bed sample (by simulations with the use of Discrete Element Method); obtaining the hydraulic tortuosity (by simulations with the use of Lattice Boltzmann Method); obtaining the geometric tortuosity (by the use of the above mentioned Path Tracking Method) and the comparison of results. All these steps are described below. We discuss not only the issues related directly to the tortuosity, but also with the porosity (as indicates the literature review it is the main factor impacting the tortuosity value) and with chosen features of the used methodologies. We show that the Path Tracking Method may be very useful due to the low demand for computing power (in opposition to Lattice Boltzmann Methods), relatively simple implementation (appropriate references are given) and the possibility of use in large granular beds. The last feature is particularly important, because the geometric method may be used even in cases, in which the application of the LBM approach is impossible. 2 Materials and methods 2.1 Materials In order to generate numerical samples of granular beds we need to know both particle sizes and the porosity of the sample. Both will be used as an input data for packing procedures. We used a real set of the SiLibeads Glass beads (Type S). First, we measured their diameters along two axes using the micrometer screw with \(\pm \,0.01\) [mm] accuracy. We found that diameters of 100 beads follow Gaussian distribution at a significance level 0.01 of the Kolmogorov–Smirnov test [28]. The mean diameter of 100 beads was found \({\bar{d}} = 6.072\) [mm] with standard deviation \(\sigma = 0.051\) [mm]. To test porosity of our sphere-packed media we loosely packed those glass beads into a graduated 200 [ml] cylinder and filled it slowly with water. We measured porosity by dividing the water volume (\(V_0\)) by the total volume of cylinder (\(V_c\)): \(\phi = V_0/V_c\). We repeated the procedure 15 times, each time with a new, dry set of marbles. We found that porosity of the sample is \(\phi = 0.41\pm \, 0.008\) what agrees with literature data [11, 40]. 2.2 Theoretical distribution of the particles forming the granular bed The YADE code used in further investigations allows introducing data from the above sieve analysis to generate particle clouds with non-uniform diameters. Due to inaccessibility of such data, to use this feature we need to convert the known distribution function to a discrete form with specific number of fractions. Such approach allows generating particle clouds with any distribution. The procedure of generating the theoretical discrete distribution from a Gaussian distribution is as follows. The range of \(\pm \,4 \sigma \) is divided into \(n_f\) (i.e. number of fractions) equal subranges. Each i-th (where \(i =1,2,\ldots , n_f\)) subrange has an attributed diameter: $$\begin{aligned} d_i = {\bar{d}} - 4 \sigma + \frac{4 \sigma }{n_f} (2 i+1). \end{aligned}$$ (4) Then, if N is the total number of spheres in the bed under creation, then the number of spheres \(n_i\) for the i-th bin is calculated as: $$\begin{aligned} n_i = N \int _{{\bar{d}} - 4 \sigma \frac{8 \sigma }{n_f} i}^{{\bar{d}} - 4 \sigma \frac{8 \sigma }{n_f} (i+1)} \frac{1}{\sigma \sqrt{2 \pi }} e^{- \left( \frac{x- {\bar{d}}}{2 \sigma } \right) ^2}. \end{aligned}$$ (5) Explanation for the above procedure is as follows. The sum of calculated integrals seen in the formula (5) obtained for each bin is equal to one. Multiplication of these values by the number of spheres that are supposed to form the granular bed as a whole, gives us the contribution of each fraction with an attributed diameter. Note, that the range \(\pm \,4 \sigma \) in case of normal distribution means that 99.994% of the population falls inside that range; in other words, 15,787 must be the minimum total population to expect the first case of value falling outside that range [16]. The result of the procedure described in this section is a cumulative curve of particle distribution given in a discrete form. 2.3 Discrete Elements Method (Radius Expansion Method) To generate a set of virtual beds for which the geometric and hydraulic tortuosity may be next calculated, Discrete Element Method (DEM) [12] implemented in the YADE code [25, 41] was used. This method allows modelling the dynamics of physical systems consisting of any number of separate bodies (e.g. particles). During the calculations, two types of equations are solved in a loop. Firstly, the displacement of every body in the system is calculated by the application of Newton’s laws, and then, the interaction forces between adjacent bodies are calculated. An important part of the algorithm is tracking contacts between objects. The main equations of the linear and the angular motion may be written as follows [12]: $$\begin{aligned} m_i \frac{d \mathbf {v}_i}{dt} = \sum _{k=1}^{n_c} \left( \mathbf {F}_{ij}^n + \mathbf {F}_{ij}^t \right) + \mathbf {F}_i^{ext} \end{aligned}$$ (6) and $$\begin{aligned} I_i \frac{d \varvec{\omega }_i}{dt} = \sum _{k=1}^{n_c} \left( r_i \times \mathbf {F}_{ij}^t \right) + \tau _{ij}^r, \end{aligned}$$ (7) where \(m_i\)—mass of the i-th body [kg], \(I_i\)—moment of inertia of the i-th body [kg\(\cdot \)m\(^2\)], \(\mathbf {v}_i\)—linear velocity of the i-th body [m/s], \(\varvec{\omega }_i\)—angular velocity of the i-th body [rad/s], \(\mathbf {F}_{ij}^n\)—normal forces between neighbouring bodies i and j [N], \(\mathbf {F}_{ij}^t\)—tangential forces between neighbouring bodies i and j [N], \(n_c\)—number of contacts between i-th body and neighbouring bodies \([-]\), \(\mathbf {F}_i^{ext}\)—external forces acting on the i-th body (e.g. gravity force) [N], \(r_i\)—distance between the contact point with the j-th body and the mass centre of the i-th body [m], \(t_{ij}^r\)—torque factor associated with the rolling friction [N\(\cdot \)m]. A virtual granular bed through DEM simulation may be obtained in three ways: (1) through generation of a particle cloud with proper diameters and size distribution, which is next compressed by moving walls until the target porosity is reached; (2) through generation of a particle cloud with proper size distribution and reduced diameters, that next are growing in a constant volume during the simulation until the target porosity is reached; (3) through adding particles one by one over the bed, which next fall down by gravity and form the bed. Independent on the approach, all interaction forces between neighbouring particles are calculated during a simulation and particles may move dependent on the local force equilibrium. In our study, the second approach (known in the literature as the Radius Expansion Method, REM) due to the much higher speed of operation was applied (Fig. 2). Some drawback of this approach is that the volume in which the particle cloud grows has to be very precisely defined. In other case, the final porosity or the final average particle diameter would be wrong. This method produces homogeneous beds with slightly varying diameters of particles (due to the differences in the local equilibrium states). We did not take the third method into account due to the conviction that the approach based on particle clouds gives more homogeneous beds. We refer here e.g. to the paper [55] in which this approach was applied in investigations of stratification effects in granular beds. Open image in new window Fig. 2 Visualisation of the Radius Expansion Method To calculate the bed volume, so-called characteristic length \(l_c\) is introduced. This variable is dependent on the number of particles, their size distribution and the whole sample porosity: $$\begin{aligned} l_c = \root 3 \of {\frac{V_s + V_p}{2}}, \end{aligned}$$ (8) where $$\begin{aligned} V_s = \sum _{i}^{n_f} {\bigg ( \sum _{j}^{n_j} {\frac{2}{3} \pi d_j^3} \bigg )} \end{aligned}$$ (9) and $$\begin{aligned} V_p = V - V_s = \frac{\phi }{1 - \phi } V_s. \end{aligned}$$ (10) \(V_s\) is the total volume of the solid part of the porous bed (particles) [m\(^3]\), \(V_p\) is the volume of the pore part of the porous bed [m\(^3]\), V is the total volume of the bed [m\(^3]\), \(n_j\) is the number of particles in the i-th fraction \([-]\) and \(d_j\) the diameter of the i-th particle [m]. 2.4 Lattice Boltzmann Method To compute the hydraulic tortuosity with Eq. (3) we need the complete velocity field in pore space [33]. We use Lattice Boltzmann Method (LBM) [7] and the Palabos numerical code [37] to reach this aim. The solver is based on kinetic theory of gases and uses the velocity distribution function \(f(\mathbf {x},t)\) in order to describe the state of the system. We solve the transport equation in the form of: $$\begin{aligned} f_i(\mathbf {x}+\mathbf {c}_i,t+1)-f_i(\mathbf {x},t)=\varOmega _i, \end{aligned}$$ (11) where \(f_i(\mathbf {x},t)\) is the discrete version of the velocity distribution function, \(\mathbf {c}_i\) is one of the lattice vectors, t is simulation time step and \(\varOmega _i\) is the collision term. We used the single relaxation time model for collision term (BGK) [47]. The resulting distribution function is time independent. Geometry obtained in a DEM simulation was converted to a grid of nodes needed for a LBM simulation. Some issues related to this step are further discussed. 2.5 Path Tracking Method To calculate the gometric toruosity in virtual beds created by the use of Discrete Element Method, so-called Path Tracking Method (PTM) [43] was applied. The algorithm is suitable for sphere-packed media. In PTM, the key role is played by tetrahedral structures formed by four adjacent particles. Such structures are used to approximate the pore space between spheres. The algorithm finds the shortest path in a given direction (usually chosen as the main direction of the transport in the media) starting from the bottom surface at various locations. The tetrahedrons are constructed one by one and the path length is calculated as a sum of distances between base triangles of all consecutive tetrahedrons. Its ratio to the sample height gives the tortuosity \(T^g\) [14, 43, 45]. The algorithm for calculating path length \(L_p\) may be summarized as follows (Fig. 3): arbitrarily choose the Initial Starting Point (ISP) at the bottom of the bed; find three particles (\(P_1\), \(P_2\) and \(P_3\)) which are the nearest to the ISP and which form a triangle in the space; calculate the coordinates of the gravity centre (GC) of the surface (in the triangle plane) through which fluid flows; move the Initial Starting Point to the Final Starting Point (FSP); in this way the first path section is perpendicular to the bottom surface of the bed; calculate the normal to the triangle; estimate coordinates of the so-called Ideal Location (IL), in which the next sphere (\(P_4\)) surrounding the path should be located; move the Ideal Location closer to the triangle plane (due to the fact, that particles forming the triangle basis may be separated from one another and in such case the fourth particle is located closer to the triangle plane); find the particle nearest to the Ideal Location; this is the Real Location (RL) of the 4th particle forming tetrahedron in the space; remove the lowest sphere from tetrahedron 1-2-3-4 to obtain the base triangle for the next tetrahedron; calculate the local length of the path inside the current tetrahedron (a local path is the line connecting centroids of base triangles of two neighbouring tetrahedrons); continue the calculations, until reaching the top surface of the bed. Open image in new window Fig. 3 Schematic representation of the PTM (particles are not in a scale) The algorithm stops when the distance between the gravity centre and the Ideal Location is smaller than the distance between the top surface of the bed and the Z coordinate of the Ideal Location. The last section is perpendicular to the top surface of the bed. The final path length \(L_p\) is calculated as the sum of all local path lengths inside tetrahedrons and two sections lying at the beginning, as well as at the end of the analyzed pore channel. The Path Tracking Method was implemented in the PathFinder code [43]. It is a freely available numerical code destined for the analysis of spatial structure of granular beds of a cylindrical or cubical shape, when data about locations and diameters of all particles in the bed are known. Currently two ways of obtaining the input data for the PathFinder code are possible: the indirect (with the use of simulation methods, e.g. DEM) and the direct (with the use of experimental methods, mainly CT/IA). The parameter set calculated by the PathFinder code contains: the bed volume (V), the volume of the porous part (\(V_p\)), the volume of the solid body (\(V_s\)), the porosity (\(\phi \)), the specific surface of the solid body (\(S_0\)) (with the use of two different definitions), the ratio of actual path length inside channels (\(L_p\)), the thickness of the porous medium (\(L_0\)) and the geometric tortuosity (\(T^g\)). Tortuosity may be calculated locally in a chosen ISP or for the whole bed in form of a scalar field [42]. Besides the geometric tortuosity, all other parameters are calculated based on analytical formulas. In the literature, alternative implementations of the PTM may be found [53, 54]. 3 Results and discussions 3.1 Virtual granular beds All virtual granular beds were created with the use of the Radius Expansion Method described in the Sect. 2.3. Every particle set has a cuboid shape with the size equal to \(2 l_c \times l_c \times l_c\). The bed size is doubled in the direction in which the geometric and the hydraulic tortuosity is calculated. The procedure of creating a virtual bed was integrated with scripts in Python programing language, needed for running the YADE code (YADE is in fact a set of libraries which has to be called from an other program - here written in Python), and contains following steps: reading the discrete cumulative curve of the particle distribution; calculating the characteristic lenght with the use of Eq. (8); performing DEM simulations; saving the results; converting the data to the other formats (for the Palabos and PathFinder codes). The generation process was performed for 12 cases, in which the number of fractions \(n_f\) was being changed. For every fraction, the bed was prepared three times with the use of different values of the seed in the random number generator. In such a way, three independent data sets for every \(n_f\) were obtained. The number of particles was set to 5000, as this is the number big enough for the reliable comparison of methods of tortuosity determination, allowing performing the analysis in acceptable computation time and not exceeding 15,787, which is the limit for the use of range in the process of generating the theoretical distribution of particles forming the granular bed (as described in the Sect. 2.2). For the purpose of performing all DEM simulations, the default so-called triaxial test was applied [41]. All settings of the triaxial stress controler were default, besides the final and the maximum multiplier factor of the REM method, which were set to 1.0004 and 1.004, respectively. To obtain a homogeneous bed, the gravity forces were turned off. In order to achieve the final porosity, the friction angle was decreased during a simulation. The starting value of this angle was set to 0.5 [rad]. The time of a typical simulation was about few minutes (Intel Inter Core 3.4 GHz processor, 32 GB RAM and Ubuntu 14.04 OS). An example of a starting and a final particle cloud recorded during a DEM simulation is visible in Fig. 4. Open image in new window Fig. 4 a Starting and b final configurations and normal forces in DEM packings in the growth algorithm. Dots visible on the right side show particle centers and connecting lines represent normal forces (\(\mathbf {F}_{ij}^n\) in Eq. 6) It can be seen (Table 1) that in every calculation, the obtained bed porosity is slightly different, even for the cases in which the number of particle fractions was constant. This is a characteristic feature of the used method. Note that the characteristic length is a little different in every case. Table 1 Results of calculations in the YADE code (all data sets) \(n_f\) \(l_c\) \(\phi _1\) \(\phi _2\) \(\phi _3\) \(\delta \phi _1 \) \(\delta \phi _2 \) \(\delta \phi _3 \) 3 0.07918935 0.4094 0.4071 0.4081 \(-\) 0.15 \(-\) 0.71 \(-\)0.46 5 0.07918937 0.4078 0.4089 0.4086 \(-\) 0.54 \(-\) 0.27 \(-\) 0.34 7 0.07920525 0.4079 0.4093 0.4081 \(-\) 0.51 \(-\) 0.17 \(-\) 0.46 9 0.07921055 0.4065 0.4081 0.4082 \(-\) 0.85 \(-\) 0.46 \(-\) 0.44 11 0.07920001 0.4069 0.4088 0.4073 \(-\) 0.76 \(-\) 0.29 \(-\) 0.66 13 0.07919474 0.4075 0.4082 0.4091 \(-\) 0.61 \(-\) 0.44 \(-\) 0.22 15 0.07920007 0.4085 0.4072 0.4098 \(-\) 0.37 \(-\) 0.68 \(-\) 0.05 17 0.07920543 0.4086 0.4075 0.4086 \(-\) 0.34 \(-\) 0.61 \(-\) 0.34 19 0.07920027 0.4074 0.4074 0.4087 \(-\) 0.63 \(-\) 0.63 \(-\) 0.32 21 0.07920579 0.4078 0.4095 0.4081 \(-\) 0.54 \(-\) 0.12 \(-\) 0.46 23 0.07919579 0.4087 0.4077 0.4098 \(-\) 0.32 \(-\) 0.56 \(-\) 0.05 25 0.07920155 0.4097 0.4080 0.4072 \(-\) 0.07 \(-\) 0.49 \(-\) 0.68 The porosity in every case is a little less than the target (experimental) porosity due to the condition of the calculations termination—the loop is finished when the current porosity is less than the target porosity. The relative errors (\(\delta \phi \)) between the obtained porosities and the experimental value are in all cases less than 1%. For this reason we can assume that the generated virtual beds represent the real bed sample very well. Open image in new window Fig. 5 Comparison between normal distribution and the particle distributions in the packings obtained with YADE at \(n_f=3,13\) and 25 (data set 1) Other important issue is to check the conformity of diameters in the real and the virtual bed. In Fig. 5 three examples of the diameter distribution (for \(n_f\) = 3, 13, 25 and data set 1) of particles included in the virtual bed, as well as the normal distribution, are shown. It can be seen that diameters in the virtual bed not coincide exactly with the theoretical distribution prescribed based on measurements of the real particles. It is the result of the way of generation of the virtual bed. The mean diameters of the spheres in generated beds are as follows: 6.0717 (for \(n_f\) = 3), 6.0774 (for \(n_f\) =13) and 6.0692 (for \(n_f\) = 25). As the average diameter of theoretical distribution is equal to 6.0720, it can be stated that the differences are negligible (0.0043%, 0.0884% and 0.0463%, respectively). The Pearson’s linear correlation coefficients [39] were also calculated between histograms for theoretical distribution and histograms for distributions of obtained virtual beds. The results (99.99% for \(n_f\) = 3, 97.81% for \(n_f\) =13 and 98.08% for \(n_f\) = 25) are satisfactory here as well. The same results were obtained for others data sets. 3.2 Hydraulic tortuosity (LBM simulations) To perform calculations in the Palabos code (see Sect. 2.4), the bed geometry obtained in a DEM simulation was transferred to the form required by the LBM. A regular, cubic lattice of nodes was used to approximate geometry and porosity of the discretized sample. The size of each grid was equal to \(2 L \times L \times L\). Due to the high demand for computing power only three cases of \(n_f\) were analysed. To find the minimal system size we discretize geometry of each sample and set each node of the lattice as solid if it lays entirely or at least by \(50\%\) inside any of the beads. Remaining nodes are left as fluid. The ratio between a number of fluid nodes to all lattice nodes gives porosity. We test various values of \(L\in {75,100,\ldots ,300}\) [lu] (lattice units) and find the minimal size \(L=200\) [lu] at which porosity does not change with an increasing grid size (Fig. 6). Open image in new window Fig. 6 Porosity of the discretized bed for LBM method in function of the system size L (numerical grid resolution in the transverse direction) Open image in new window Fig. 7 Exemplary streamlines in the fluid flow trough pore channels of the virtual granular bed (\(n_f=3\), data set 1) In the next step, we measured hydraulic tortuosity using the same samples. First, we simulated the fluid flow at the pore-level. We set periodic boundaries in the flow direction. The no-slip conditions were used at horizontal walls. Also, the no-slip conditions were used for solid nodes. The relaxation time was \(\tau =1\) resulting in the kinematic viscosity \(\nu = 1/6\) in lattice units [38]. The fluid was driven by the external gravity and the creeping flow was kept (the Reynolds number Re\( = d q/ \nu < 1\), where q is a Darcy flux, d is the mean grain diameter). An exemplary result of the simulation is given in Fig. 7. The data obtained for \(L = 200\) are collected in Table 2. Hydraulic tortuosity was calulated using Eq. 3. Table 2 Results of calculations in the Palabos code for \(L = 200\) (data averaged over three independent samples for chosen \(n_f\)) \(n_f\) \(\phi _1\) \(\phi _2\) \(\phi _3\) \(\delta \phi _1\) \(\delta \phi _2\) \(\delta \phi _3\) \(T^h_1\) \(T^h_2\) \(T^h_3\) 3 0.4141 0.4114 0.4123 1.00 0.33 0.55 1.2310 1.2305 1.2309 13 0.4120 0.4125 0.4134 0.50 0.60 0.83 1.2316 1.2338 1.2314 25 0.4114 0.4121 0.4116 0.33 0.51 0.39 1.2314 1.2331 1.2340 To check the finite-size effects we gradually increased the lattice size from \(L=75\) to \(L=300\). We found (Fig. 8) that tortuosity grows monotonically following an exponential function. Thus, in order to estimate the value of \(T^h\) for infinite lattice, we fit an exponential function of the form $$\begin{aligned} T^h(L) = T(\infty ) - a e^{(-b L)} \end{aligned}$$ (12) where \(T(\infty ) = 1.24\) and fit parameters a and b are equal to 0.13 and \(-\,0.014\), respectively. Details of the procedure are similar to those used i.e. in [32]. Please note that \(T^h \equiv {T}(\infty )\) from the above mentioned fitting procedure. 3.3 Geometric tortuosity (PTM calculations) The data obtained in DEM simulations was converted to format readable by PathFinder code (see Sect. 2.5), which was next used to calculate the porosity and the geometric tortuosity. In all calculations default settings were used [43]. The results of calculations are collected in Table 3. It can be seen (Fig. 9) that the porosity in every case is a little higher than the experimental porosity. The relative errors are usually slightly higher than these obtained in the direct calculations in the YADE code. However, we assume that these values are still very well representative. This is all the more justified that the differences in the porosity between particular ways of its calculation are less than the measurement error (\(\phi _{exp} \pm \, \varDelta \phi _{exp}\)). The geometric tortuosity was calculated 25 times for every virtual bed in the direction of the longest side of the calculation domain (in Table 2 only the average values are shown). Such a large number of values results from the fact that the path lengths inside channels were calculated in different zones of the virtual bed (Fig. 10) [43]. In this way, an average value was obtained, which better represents the tortuosity in this bed. In Fig. 11 the average values of geometric tortuosity for all three data sets and different particle fractions are shown. The average value of the tortuosity obtained with the use of PTM for all data sets and all values of \(n_f\) is equal to 1.2147. In this figure, the results reported by Wang (see Sect. 1) are also shown. The relative error between the above mentioned average value and the results shown by Wang is equal to 1.43 %. In such a way, the usefullnes of the PTM is validated not only by LBM simulations, but by an independent results as well. Open image in new window Fig. 8 Hydraulic tortuosity of the discretized bed for LBM method in function of the system size L (numerical grid resolution in the transverse direction) Table 3 Results of calculations in the PathFinder code (all data sets) \(n_f\) \(\phi _1\) \(\phi _2\) \(\phi _3\) \(\delta \phi _1\) \(\delta \phi _2\) \(\delta \phi _3\) \(T^g_1\) \(T^g_2\) \(T^g_3\) 3 0.4142 0.4117 0.4124 1.02 0.41 0.59 1.2123 1.2117 1.2059 5 0.4126 0.4130 0.4130 0.63 0.73 0.73 1.2087 1.2149 1.2066 7 0.4125 0.4136 0.4124 0.61 0.88 0.59 1.2090 1.2154 1.2118 9 0.4112 0.4125 0.4122 0.29 0.61 0.54 1.2143 1.2127 1.2125 11 0.4116 0.4133 0.4114 0.39 0.80 0.34 1.2338 1.2092 1.2179 13 0.4122 0.4127 0.4135 0.54 0.66 0.85 1.2121 1.2212 1.2192 15 0.4129 0.4118 0.4144 0.71 0.44 1.07 1.2106 1.2155 1.2097 17 0.4142 0.4112 0.4128 1.02 0.29 0.68 1.2167 1.2094 1.2231 19 0.4121 0.4122 0.4132 0.51 0.54 0.78 1.2132 1.2173 1.2338 21 0.4126 0.4138 0.4125 0.63 0.93 0.61 1.2099 1.2095 1.2164 23 0.4126 0.4119 0.4141 0.63 0.46 1.00 1.2166 1.2271 1.2095 25 0.4142 0.4122 0.4117 1.02 0.54 0.41 1.2128 1.2127 1.2148 No systematic changes with \(n_f\) were found, but, for \(n_f>10\) values of \(T^g\) are found to be more scattered, higher and encumbered with larger statistical errors of the mean. 4 Discussion 4.1 Direct data comparison In our simulations we show that the porosity in all used numerical models is very close to the experimental value. All the realative errors were close to 1% or (usually) even less than 1%. The particle distribution obtained by the use of REM follows the theoretical distributions at a satisfactory level. We noted that the higher values of the \(n_f\) give more representative virtual beds due to the better reconstruction of the Gaussian distribution. Finally we can assume that if the porosity value is correct, than the other, more complicated and porosity dependent, physical quantities may be also determined from the same virtual model (e.g. tortuosity). Open image in new window Fig. 9 Comparison of the porosity value obtained in the all used methods Open image in new window Fig. 10 Examples of path lines calculated by the use of the Path Tracking Method. Dots in the left figure represent the particle centres We observed that average values of hydraulic and geometric tortuosity are very similar and slightly higher than the results shown by Wang. Probably the reason is that Wang used particles with the same diameter. It suggests that tortuosity may increase with the increase of the standard deviation of the particle distribution. However, to confirm this, more investigations are needed. The main difference between the both used methodologies is the scattering of the obtained tortuosity values (Fig. 12). In the hydraulic approach, the tortuosity value is averaged over all nodes (available for the fluid) of the lattice grid. The local values of the tortuosity are not known. For this reason the average value is not very sensitive to the arrangement of the particles in the bed (i.e. the used data set). We suppose that the difference between average tortuosities obtained for different data sets would decrease with increase in the number of particles in the bed. This feature may be treated as an advantage (we obtain a very similar value every time) or as an drawback (we have no information about the possible local fluctuations of the tortuosity). In this context, the geometrical approach gives more information (e.g. in paper [42] two-dimensional tortuosity fields in the direction perpendicular to the main flow direction were calculated). We obtain not only one average value, but also the data on the range of the geometric tortuosity. Open image in new window Fig. 11 Relationships between the number of fractions and the tortuosity Open image in new window Fig. 12 The effect of overlapping the result sets It should be underlined that the range, as well as the average value of the geometric tortuosity, is in some way dependent on the number of the Initial Starting Points used in calculations (see Sect. 2.5). In the presented study, always 25 points were used. But, if we change the number of the Initial Starting Point to 9 (3 \(\times \) 3) or 10,000 (100 \(\times \) 100), we obtain very similar values of the geometric tortuosity: 1.2166 ± 0.0045 and 1.2047 ± 0.0002, respectively. These values vere averaged over all \(n_f\) and data sets, and are very similar to the previous presented value (1.2147). As it can be seen, the PTM is not very sensitive to the number of paths used to calculate the average tortuosity. However, in every case the geometrical tortuosity representative for the whole bed should be calculated from a larger number of individual values. 4.2 Indirect data comparison The observed relative differences between geometric and hydraulic tortuosity in granular beds may be expressed in terms of the anisotropy of an inclined capillary model of the flow [6]. In that model particles paths follow straight paths inclined by some angle \(\alpha \) (varying from path to path) with respect to the main axis of the flow. The tortuosity of each channel is equal to $$\begin{aligned} T=\cos ^{-1} \alpha =L_{eff}, \end{aligned}$$ (13) where \(\alpha \) is the inclination angle. To describe the whole medium, \(L_{eff}\) may be taken as the length of an average flow axis. Using the above relation and known values of \(T^g\) and \(T^h\) the inclination angles for geometric and hydraulic tortuosity are evaluated \(\alpha ^g=34^{\circ }\) and \(\alpha ^h=36^{\circ }\), respectively. From this, the absolute difference \(2^{\circ }\) and the relative difference \(\varDelta \alpha =3^{\circ }\) is calculated. Although this analysis is not perfect and suitable for simplified model of inclined channel flow, it gives some rough idea how big error would be introduced by taking \(T^g\) instead of \(T^h\) for the anisotropy angle estimations. In the real porous medium, due to a non-straight, wiggly paths of particles, one may expect that this difference will be decreased. Similarly, if we consider the empirical law for permeability (k) as a function of geometric parameters and tortuosity in which \(k\propto 1/T^2\) [23, 51], the maximum deviation $$\begin{aligned} \varDelta k=100|k^g - k^h|/k_h \approx 6\%. \end{aligned}$$ Here \(k^g\) is based on geometric and \(k^h\) is based on the hydraulic estimation of tortuosity. The same \(6\%\) deviation appears between effective diffusion coefficients in porous media related to diffusion in fluids with \(D_{eff}={D} T^{-2}\) [8]. Our results may be also used to analyze sphere-packed porous media in terms of the Archie law. This is an empirical relation for tortuosity as function of porosity: $$\begin{aligned} T=\phi ^{-b}, \end{aligned}$$ (14) where \(\phi \) is porosity and b is the free parameter related to a cementation factor in porous media [5, 10]. Using Eq. (14), for porosity \(\phi =0.41\) (Sect. 2.1) and tortuosity 1.24 (Sect. 3.2) and 1.21 (Sect. 3.3) we get \(b\approx 0.24\) and \(b\approx 0.21\). This is close to \(b=0.25\) calulated previously by numerical fitting of the Eq. 12 to tortuosity of the random, overlapping spheres three-dimensional porous medium [33]. This is an unexpected result which suggests that the cementation factor in the Archie law is independent of the way spheres are packed and arranged. We suggest that this is related to the shape of particles packed in the bed. This finding deserves further study and analysis but it is out of the scope of this work. 4.3 Possibility of investigations of the bed anisotropy In the paper, we discuss results of tortuosity calculations in only one space direction. However, both used methods may be used to determine the tortuosity in others directions. Having such data, any three-dimmensional model destined to analysis porous media may be used. The so-called Porous Media Model available in the popular ANSYS Fluent code may serve here as an example. The linear (Darcy) term in this model has a form [1]: $$\begin{aligned} s_i = - \sum _{j=1}^{3} D_{ij} \mu v_j, \end{aligned}$$ (15) where \(s_i\)—source of forces for the i-th space dimension [N/(kg \(\cdot \mathrm{m}^3)]\), \(\mu \)—dynamic viscosity coefficient [kg/(m \(\cdot \) s)], \(v_j\)—the j-th component of velocity [m/s], \(\mathbf {D}\)—a matrix with diagonal terms equal to 1 / k [1/m\(^2]\). Now, if we assume that the term \(\mathbf {D}\) is described e.g. by the Kozeny–Carman law (1), we can predict pressure distribution in a three-dimensional computational space. Note that this law does not define, which kind of the tortuosity should be used. We are aware that the anisotropy effects in granular beds consisting of spherical or quasi spherical particles are usually not very significant. However, in beds with high standard deviation, the particle diameters may be very different and the particles may tend to aggregation or form layers due to sizes. The way of generating such porous beds (real or virtual) may also affect in the level of the anisotropy. Except that, in a general case, the particle distribution may not follow Gaussian distributions and may by a mixture of many fractions with different distributions. The PTM may be applied for porous beds of each above kind. 4.4 Time and memory consumption We observed that computation time used to calculate the geometrical tortuosity was signicantly less than the time needed for LBM solver to simulate fluid flow required for \(T^h\) computation. For example, to determine single path in geometric methods using the earlier mentioned computer less than 6 [s] was necessary in average. In LBM, computation of fluid flow at reasonable accuracy required for \(T^h\) measurement varied with mesh size L and it was less than 20 [min] for L = 75, more than 1.5 [h] for L = 200 and more than 5 [h] for L = 300. Simulations in Palabos were performed in parallel mode with the use of 4 processors. Similarly, the memory consumption of LBM was large. It was 3GB for the D3Q27 model at L = 200 and 11GB for L = 300. For the same sample, PathFinder stores only grain positions and diameters, which required less than 0.2 MB for 5000 beds. Even using the sparse data structures in Palabos, that could drop the memory down to 45% [29], the memory required by LBM is orders of magnitude larger than required by geometrical algorithms. By taking into consideration lower time and memory resources required by geometric methods, we suggest that geometric methods will allow to estimate tortuosity in much larger, granular systems. This may also simplify measurement of other important parameters (i.e. permeability) in large scale samples. 5 Conclusions We found that geometric tortuosity in granular beds with slightly varying diameter is \(T^g \approx 1.21\), which is close to the hydraulic counterpart \(T^h \approx 1.24\). This value is also close to the hydraulic tortuosity given by Wang, who used a similar methodology of obtaining this quantity. We suggest that the applied geometrical method may serve as an attractive alternative to hydraulic tortuosity in large granular systems. The advantages of the Path Tracking Method are as follows: a very short time of calculations; the possibility of use in granular beds consisting of a large number of particles; the possibility of calculating tortuosities of individual geometric path lines in one virtual bed; no need of geometry conversion; no need of model calibration. In turn, in the LBM approach the demand for computing power is large, what limits its use. The main drawback of the PTM is that it is not related directly to the transport equations and that this method is appropriate only for beds consisting of spherical or quasi-spherical particles. However, the area in which such kinds of particles occurs is still relatively large. Notes Acknowledgements We would like to thank anonymous reviewers and dr hab. Zbigniew Koza for useful comments on the manuscript. Funding This study was funded by Polish Ministry of Science and Higher Education in the frames of the statutory research. Compliance with ethical standards Conflict of interest The authors declares that they have no conflict of interest. References 1. ANSYS Inc.: Porous media conditions. In: ANSYS Fluent User’s Guide, Release 13. (2010) 2. Allen, R., Sun, S.: Investigating the role of tortuosity in the Kozeny–Carman equation. In: International Conference on Numerical and Mathematical Modeling of Flow and Transport in Porous Media, Dubrovnik, Croatia, 29 September–3 October 2014 (2014)Google Scholar 3. Aaltosalmi, U.: Fluid flow in porous media with the lattice-Boltzman method. Ph.D. thesis, University of Jyväskylä, Finland (2005)Google Scholar 4. Ahmadi, M., Mohammadi, S., Hayati, A.: Analytical derivation of tortuosity and permeability of monosized spheres: a volume averaging approach. Phys. Rev. E 83(2), 026312 (2011)ADSCrossRefGoogle Scholar 5. Archie, G.: The electrical resistivity log as an aid in determining some reservoir characteristics. Trans. Am. Inst. Min. Metall. Pet. Eng. 146, 54–62 (1942)Google Scholar 6. Bear, J.: Dynamics of Fluids in Porous Media. Courier Dover Publications, New York (1972)zbMATHGoogle Scholar 7. Bhatnagar, P., Gross, E., Krook, M.: A model for collisional processes in gases. i. Small amplitude processes in charged and neutral one-component system. Phys. Rev. 94, 511–524 (1954)ADSCrossRefGoogle Scholar 8. Boudreau, B., Meysman, F.: Predicted tortuosity of muds. Geology 34, 693 (2006)ADSCrossRefGoogle Scholar 9. Carman, P.: Fluid flow through a granular bed. Trans. Inst. Chem. Eng. (Jubilee Supplement) 75, 32–48 (1997)CrossRefGoogle Scholar 10. Clennell, M.B.: Tortuosity: a guide through the maze. Geol. Soc. Spec. Publ. 122(1), 299–344 (1997). ADSCrossRefGoogle Scholar 11. Cooke, A., Rowe, R.: Extension of porosity and surface area models for uniform porous media. J. Environ. Eng. 125, 126–136 (1999)CrossRefGoogle Scholar 12. Cundall, P., Strack, O.: A discrete element model for granular assemblies. Géotechnique 29, 47–65 (1979)CrossRefGoogle Scholar 13. Duda, A., Koza, Z., Matyka, M.: Hydraulic tortuosity in arbitrary porous media flow. Phys. Rev. E 84(3), 036319 (2011)ADSCrossRefGoogle Scholar 14. Dudda, W., Sobieski, W.: Modification of the pathfinder algorithm for calculating granular beds with various particle size distributions. Tech. Sci. 17(2), 135–148 (2014)Google Scholar 15. Ebner, M., Chung, D.W., García, R., Wood, V.: Tortuosity anisotropy in lithium-ion battery electrodes. Adv. Energy Mater. 4, 1301278 (2013)CrossRefGoogle Scholar 16. El-Haik, B.: Axiomatic Quality: Integrating Axiomatic Design with Six-Sigma, Reliability, and Quality Engineering. Wiley, Hoboken (2005)CrossRefGoogle Scholar 17. Gao, H.Y., He, Y.H., Zou, J., Xu, N.P., Liu, C.: Tortuosity factor for porous feal intermetallics fabricated by reactive synthesis. Trans. Nonferrous Met. Soc. 22, 2179–2183 (2017)CrossRefGoogle Scholar 18. Gommes, C., Bons, A., Blacher, S., Dunsmuir, J., Tsou, A.: Practical methods for measuring the tortuosity of porous materials from binary or gray-tone tomographic reconstructions. AIChE J. 55, 2000–2012 (2009)CrossRefGoogle Scholar 19. Johnson, D., Plona, T., Scala, C., Pasierb, F., Kojima, H.: Tortuosity and acoustic slow waves. Phys. Rev. Lett. 49, 1840–1844 (1982)ADSCrossRefGoogle Scholar 20. Kochański, J., Kaczmarek, M., Kubik, J.: Determination of permeability and tortuosity of permeable media by ultrasonic method. Studies for sintered bronze. J. Theor. Appl. Mech. 39, 923–928 (2000)Google Scholar 21. Kong, W., Zhang, Q., Gao, X., Zhang, J., Chen, D., Su, S.: A method for predicting the tortuosity of pore phase in solid oxide fuel cells electrode. Int. J. Electrochem. Sci. 10, 5800–5811 (2015)Google Scholar 22. Koponen, A., Kataja, M., Timonen, J.: Tortuous flow in porous media. Phys. Rev. E 54, 406–410 (1996). ADSCrossRefGoogle Scholar 23. Koponen, A., Kataja, M., Timonen, J.: Permeability and effective porosity of porous media. Phys. Rev. E 56, 3319–3325 (1997). ADSCrossRefGoogle Scholar 24. Kozeny, J.: über kapillare leitung des wassers im boden (in German). Akademie der Wissenschaften in Wien (Sitzungsberichte 136/2a), 271–306 (1927)Google Scholar 25. Kozicki, J., Donze, F.: Yade-open dem: an open-source software using a discrete element method to simulate granular material. Eng. Comput. 26(7), 786–805 (2009). CrossRefzbMATHGoogle Scholar 26. Lane, N.: Numerical studies of flow in porous media using an unstructured approach. Ph.D. thesis, Louisiana State University and Agricultural and Mechanical College, Baton Rouge, USA (2011)Google Scholar 27. Le, L., Zhang, C., Ta, D., Lou, E.: Measurement of tortuosity in aluminum foams using airborne ultrasound. Ultrasonics 50, 1–5 (2010)CrossRefGoogle Scholar 28. Lilliefors, H.: On the Kolmogorov–Smirnov test for normality with mean and variance unknown. J. Am. Stat. Assoc. 62(318), 399–402. (1967). Accessed 21 Sept 2018CrossRefGoogle Scholar 29. Liu, H., Zhang, Y., Valocchi, A.: Lattice Boltzmann simulation of immiscible fluid displacement in porous media: homogeneous versus heterogeneous pore network. Phys. Fluids (2015). CrossRefGoogle Scholar 30. Maier, R., Kroll, D., Kutsovsky, Y., Davis, H., Bernard, R.: Simulation of flow through bead packs using the lattice Boltzmann method. Phys. Fluids 10(1), 60–74 (1998)ADSCrossRefGoogle Scholar 31. Manwart, C., Aaltosalmi, U., Koponen, A., Hilfer, R., Timonen, J.: Lattice-Boltzmann and finite-difference simulations for the permeability for three-dimensional porous media. Phys. Rev. E 66(1), 016702 (2002). ADSCrossRefGoogle Scholar 32. Matyka, M., Khalili, A., Koza, Z.: Tortuosity–porosity relation in porous media flow. Phys. Rev. E 78, 026306 (2008)ADSCrossRefGoogle Scholar 33. Matyka, M., Koza, Z.: How to calculate tortuosity easily? AIP Conf. Proc. 1453, 17–22 (2012)ADSCrossRefGoogle Scholar 34. Montes, J., Cuevas, F., Cintas, J.: Electrical and thermal tortuosity in powder compacts. Granul. Matter 9(6), 401–406 (2007). CrossRefGoogle Scholar 35. Nabovati, A., Sousa, A.: New Trends in Fluid Mechanics Research, Chap. Fluid Flow Simulation in Random Porous Media at Pore Level Using Lattice Boltzmann Method. Springer, Berlin (2007)Google Scholar 36. Nwaizu, C., Zhang, Q.: Characterizing tortuous airflow paths in a grain bulk using smoke visualization. Canad. Biosyst. Eng. 57(1), 3.13–3.22 (2015). CrossRefGoogle Scholar 37. Palabos: (2015). Accessed 21 Sept 2018 38. Pan, C., Luo, L.S., Miller, C.: An evaluation of lattice Boltzmann schemes for porous medium flow simulation. Comput. Fluids 35(8–9), 898–909 (2006)CrossRefGoogle Scholar 39. Pearson, K.: The problem of the random walk. Nature 72, 294 (1905)ADSCrossRefGoogle Scholar 40. Ribeiro, A., Neto, P., Pinho, C.: Mean porosity and pressure drop measurements in packed beds of monosized spheres. Int. Rev. Chem. Eng. 2, 40–46 (2010)Google Scholar 41. Smilauer, V., et al.: Yade Documentation, 2nd edn. The Yade Project. (2015). Accessed 21 Sept 2018 42. Sobieski, W.: The use of path tracking method for determining the tortuosity field in a porous bed. Granul. Matter 18(1–9), 72 (2016). CrossRefGoogle Scholar 43. Sobieski, W., Lipiński, S.: Pathfinder User’s Guide. (2015). Accessed 21 Sept 2018 44. Sobieski, W., Lipiński, S.: The analysis of the relation between porosity and tortuosity in granular beds. Tech. Sci. 20, 75–85 (2016)Google Scholar 45. Sobieski, W., Zhang, Q., Liu, C.: Predicting tortuosity for airflow through porous beds consisting of randomly packed spherical particles. Transp. Porous Med. 93, 431–445 (2012)CrossRefGoogle Scholar 46. Starly, B., Yildirim, E., Sun, W.: A tracer metric numerical model for predicting tortuosity factors in three-dimensional porous tissue scaffolds. Comput. Methods Prog. Bio. 87, 21–27 (2017)CrossRefGoogle Scholar 47. Succi, S.: The Lattice Boltzmann Equation for Fluid Dynamics and Beyond. Clarendon Press, New York (2001)zbMATHGoogle Scholar 48. Vallabh, R.: Tortuosity in fibrous porous media using computational fluid dynamics. Ph.D. thesis, North Carolina State University, USA (2009)Google Scholar 49. Wang, P.: Lattice Boltzmann simulation of permeability and tortuosity for flow through dense porous media. Math. Probl. Eng., 694350 (2014). ADSMathSciNetGoogle Scholar 50. Wu, Y., van Vliet, L., Frijlink, H., van der Voort Mssrschalk, K.: The determination of relative path length as a measure for tortuosity in compacts using image analysis. Eur. J. Pharm. Sci. 28, 433–440 (2006)CrossRefGoogle Scholar 51. Wyllie, M., Gregory, A.: Fluid flow through unconsolidated porous aggregates. Ind. Eng. Chem. 47(7), 1379–1388 (1955). CrossRefGoogle Scholar 52. Yu, B., Li, J.: A geometry model for tortuosity of flow path in porous media. Chin. Phys. Lett. 21, 1569–1571 (2004)ADSCrossRefGoogle Scholar 53. Yue, R.: Modeling pore structures and airflow in grain beds using discrete element method and pore-scale models. Ph.D. thesis, University of Manitoba (2017)Google Scholar 54. Yue, R., Zhang, Q.: A pore-scale model for predicting resistance to airflow in bulk grain. Biosyst. Eng. 155, 142–151 (2017)CrossRefGoogle Scholar 55. Zhang, D., Zhou, Z., Pinson, D.: Dem simulation of particle stratification and segregation in stockpile formation. EPJ Web Conf. 140, 1–4 (2017). CrossRefGoogle Scholar Copyright information © The Author(s) 2018 Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Authors and Affiliations Wojciech Sobieski1Email authorView author's OrcID profileMaciej Matyka2Jarosław Gołembiewski2Seweryn Lipiński1View author's OrcID profile1.Faculty of Technical SciencesUniversity of Warmia and MazuryOlsztynPoland2.Faculty of Physics and AstronomyUniversity of WrocławWrocławPoland

This is a preview of a remote PDF:

Wojciech Sobieski, Maciej Matyka, Jarosław Gołembiewski, Seweryn Lipiński. The Path Tracking Method as an alternative for tortuosity determination in granular beds, Granular Matter, 2018, 72, DOI: 10.1007/s10035-018-0842-x