Calculating the Binary Tortuosity in DEM-Generated Granular Beds

Abstract

In this paper, a methodology of calculating the tortuosity in three-dimensional granular beds saved in a form of binary geometry with the application of the A-Star Algorithm and the Path Searching Algorithm is presented. The virtual beds serving as examples are prepared with the use of the Discrete Element Method based on data of real, existing samples. The obtained results are compared with the results described in other papers (obtained by the use of the Lattice Boltzmann Method and the Path Tracking Method) as well as with the selected empirical formulas found in the literature. It was stated in the paper that the A-Star Algorithm gives values similar (but always slightly underestimated) to the values obtained via approaches based on the Lattice Boltzmann Method or the Path Tracking Method. In turn, the Path Searching Algorithm gives results in the same value range as popular empirical formulas and additionally it is approximately two times faster than the A-Star Algorithm.

1. Introduction

The prediction of the pressure drop occurring during fluid flows is one of the most important tasks in mass or heat transport in porous media. If the dynamics of the fluid flow is correctly mapped, then the more complicated phenomena or processes may be taken into account. Despite many years of investigations into the relationships between features of porous media and pressure losses, the problem is still not fully resolved. The main difficulty results from the huge variety of porous materials occurring in nature as well in industry.

In the literature, many laws for calculating the pressure drop in fluid flow systems with porous media can be found; however, most of them may by saved in one unified form [1]

$$-\frac{dp}{dx}=A\left(\mathsf{\Phi }\right)\left[\mu {\overrightarrow{v}}_f\right]+B\left(\mathsf{\Phi }\right)\left[\rho {\overrightarrow{v}}_f^2\right]$$

(1)

where p is the pressure [Pa], x is a coordinate along which the pressure drop occurs [m], $A\left(\mathsf{\Phi }\right)$ and $B\left(\mathsf{\Phi }\right)$ are two generalized terms, dependent on the set $\mathsf{\Phi }$ characterizing the spatial structure of the porous medium, $\mu$ is the dynamic viscosity of the fluid [kg/ms], $\rho$ is the density of the fluid [kg/m${}^3]$, and ${\overrightarrow{v}}_f$ is the filtration velocity [m/s].

The set $\mathsf{\Phi }$ may be defined as follows,

$$\mathsf{\Phi }=\{d,\varphi (V,V_f),ϵ(V,V_s),\tau (L_p,L_0),S_0(S_S,V,V_s),\mathsf{\Psi },...\}$$

(2)

where d is a representative (e.g., average) particle diameter [m], $\varphi$ is the porosity $[-]$, V is the volume of the porous body [m${}^3]$, $V_f$ is the volume of the pore space (filled by a fluid) [m${}^3]$, $\epsilon$ is the volume fraction $[-]$, $V_s$ is the volume of the solid part of the porous body [m${}^3]$, $\tau$ is the tortuosity $[-]$, $L_p$ is the actual path length inside pore space [m], $L_0$ is the thickness of the porous medium [m], $S_0$ is the specific surface of the porous body [1/m], $S_s$ is the inner surface of the solid body [m${}^2]$, and $\psi$ is a sphericity coefficient $[-]$ (less than 1 in general and equal to 1 when the particles are spherical in shape).

In the literature, one can find many formulas for $A\left(\mathsf{\Phi }\right)$ and $B\left(\mathsf{\Phi }\right)$ terms matching the topology of Equation (1), e.g., Darcy (1856) [2], Forchheimer (1901) [3], Kozeny-Carman (1937) [4,5], Ergun (1952) [6], Mian (1992) [7], Skjetne (1999) [8], Samsuri et al. (2000) [9], Belyadi (2006) [10], Belyadi (2006) [11], Lord et al. (2006) [12], and Wu et al. (2006) [13], but almost all of them are related to the geometrical structure of the porous media (elements of the set $\mathsf{\Phi }$).

In a porous body consisting of spherical particles (a granular bed), most elements of the set $\mathsf{\Phi }$ may be easy calculated in an analytical way on the basis of the particle locations and sizes. In fact, the tortuosity, defined as the ratio of the actual path length inside pore channels to the thickness of the porous medium, is the only parameter which is difficult to determine.

There are three main approaches to the calculation of tortuosity in porous media:

• Direct analysis of the pore channels geometry. The computational space may be expressed directly in a form of the vector geometry, or indirectly, by a grid of nodes, cells, elementary volumes, pixels, or voxels. Algorithms may have analytical origins or may be based on different discrete techniques. In this approach the tortuosity is expressed as follows [4,5],

  $$\tau =\frac{L_p}{L_0}$$

  (3)
  If the tortuosity is calculated directly, on the basis of shapes of pore channels, then it is called the geometric or geometrical tortuosity. Examples of models based on pore channel geometry may be found in [14,15,16].
• The analysis of transportation properties of fluids flowing through pore channels. In this approach [17,18]

  $$\tau =\frac{〈\left|v\right|〉}{〈v_x〉}$$

  (4)

  where $\left|v\right|$ is the the absolute value of local flow velocity obtained for a creeping flow, $v_x$ is the X-component of velocity (where X is the direction of main flow), and is a spatial average over the pore space. The tortuosity calculated with the use of Equation (4) is in the literature called the hydraulic tortuosity. If a velocity field is available, than the so-called streamline tortuosity may be also calculated [19]. Other examples of calculating the hydraulic tortuosity may be found in [20,21,22].
• The analysis of diffusional properties of porous media and the application of the following relationship [23,24],

  $$D_{eff}=D\frac{\epsilon }{{\tau }_f}$$

  (5)

  where $D_{eff}$ is the effective diffusivity [m${}^2$/s], D is the intrinsic diffusivity of the conductive phase [m${}^2$/s], $\epsilon$ is the volume fraction of the conductive phase $[-]$, and ${\tau }_f$ is the tortuosity factor [-] defined as the square of the tortuosity. The diffusional approach is applied, e.g., in [25,26,27].

Based on the above-mentioned methodologies, as well as on the experimental measurements that are not described here, different empirical formulas were proposed. In

Table 1. Examples of empirical formulas serving to calculate the tortuosity for packed beds.

No.	Source	Formula	$\mathit{\tau }(\mathit{\varphi }=0.413)$
1	Bartell & Osterhof (1928) [28,29]	$\tau =0.5\pi$	$1.5708$
2	Carman (1937) [4,5]	$\tau =\sqrt{2}$	$1.4142$
3	Weissberg (1963) [28,30]	$\tau =1-0.49ln\left(\varphi \right)$	1.4333
4	Du Plessis & Masliyah (1988) [28,30]	$\tau =\frac{\varphi }{1-{(1-\varphi )}^{2/3}}$	1.3816
5	Comiti & Renaud (1989) [28,30,31]	$\tau =1-Cln\left(\varphi \right)$	1.3626 for $C=0.41$ [29]
6	Boudreau (1996)	$\tau =\sqrt{1-ln\left({\varphi }^2\right)}$	1.6639
7	Lanfrey et al. (2010) [29]	$\tau =1.23\frac{{(1-\varphi )}^{4/3}}{\varphi k^2}$, where k is a shape factor	1.4638 for $k=1$

, examples of such formulas destined for packed beds are presented. It may be seen that the tortuosity is usually treated as a direct function of the porosity. In the last column, some examples of results illustrating the range of the tortuosity for a chosen porosity are shown (see Section 2.1).

In the context of the hereby investigation different methodologies based on the so-called Random Walk technique [32] are particularly interesting. In this approach a virtual object is created, which randomly moves in the available calculational space.

Nakashima and Watanabe (2002) [33] applied the Random Walk method to calculate the tortuosity in a granular bed consisting of spherical particles with the average diameter equal to 2.11 [mm] and the standard deviation 0.06 [mm]. The geometry of the porous medium is described by a discrete lattice of voxels. The walker executes a random jump to one of the six nearest unoccupied sites. If the randomly selected site or a voxel is occupied by an obstacle or solid, the jump is not performed. This model shares many features with the Path Searching Algorithm developed by the author and described in Section 2.4.

Boudreau and Meysman (2015) [34] proposed a geometrical tortuosity model which serves for predicting the tortuosity in marine muds represented by the nonoverlapping disks arranged in layers and located randomly. In this model a virtual walker is defined, attempts to move through the pore space in a direction parallel to the cylindrical axis of the disks. In their model, if the walker is in the pore space between the disks, it will simply move along the main direction. If the walker encounters a disk, it will choose a random direction and take a straight line along the disk surface, until it reaches an edge of the disk and the movement in the main direction will be again possible. The idea is quite similar to the Path Searching Algorithm described in Section 2.4, but the model developed by Boudreau and Meysman is analytical and destined only for obstacles having one shape. In this context the Path Searching Algorithm is more universal.

Huang et al. (2019) [35] propose the so-called Random Walking Particle Tracking (RWPT) method. The random walk of a particle in an advection velocity field is modeled by a stochastic differential equation containing certain terms responsible for the particle displacement, the advection displacement, and a random walk displacement. The geometry is saved in the form of a structural grid, in which each cell represents the pore space or the solid part of the porous medium. Walkers follow a straight trajectory and change this trajectory only as a result of collisions with walls.

Amien et al. (2019) [36] used the Simple Neurite Tracer (SNT) to analyze the tortuosity in four kind of digital samples of porous rock model with the size of 256×256 pixels. NST is a plugin of ImageJ software designed to allow easy semi-automatic tracing of neurons or other filament-like structures (e.g., microtubules and blood vessels) through either 2D images or 3D image stacks [37]. Details related to the method applied are not described, but on the home page of the software it is mentioned that A-Star algorithm is also implemented in this code. In the context of the paper, it is quite important to underline that this approach may be applied only in specific geometries.

Cao et al. (2019) [38] presented an analytical model of calculating the tortuosity in cement-based materials containing aggregates with spherical or quasi-spherical shape. Paths are expressed as straight lines. However, when the path meets an obstacle, it bypasses it and returns to the original trajectory. Their results seem to be debatable, because the obtained tortuosity values are very high (between 25 and 50) and this does not correspond with the schemas presented in this work, in which paths are relatively straight.

In the hereby paper a methodology, based on porous media in which the geometry is saved in a form of binary matrix of zeros and ones, is proposed. The idea was inspired by the Lattice Boltzmann Method [39]. Formula (3), in which the path length is calculated by the use of the A-Star Algorithm or the Path Searching Algorithm, is applied to calculate the tortuosity. In Figure 1, the main idea is presented. Black circles represent the solid part of the porous body. In turn, the white circles represent the grid nodes belonging to the pore space. An example of path is marked by red color. It is important to note that the trajectory depends on the assumptions of the algorithm used (for this reason in Figure 1 two variants of paths are shown). The tortuosity obtained in the described way is called here the binary tortuosity (it is a kind of the geometrical tortuosity). The main difference between binary algorithms described in the paper and Random Walk methods is that the binary algorithms are determined by rigid rules (except for one specific case in the Path Searching Algorithm) and not by random choices.

It should be underlined that the A-Star algorithm is relatively often applied to calculate paths in a space containing different obstacles or limitations. However, it is usually used in other research areas, such as mobile robots, vehicle traffic, transportation, or logistics [40,41]. Papers in which the A-Star algorithm is used directly in the context of porous media are very difficult to find. However, the mentioning of the implementation of A-Star algorithms in the ImageJ code suggests that such papers most probably exist. The implementation of both algorithms (for in the case of the A-Star Algorithm there are two variants) was performed by the author in 2019. The Patch Searching Algorithm was developed by the author the same year.

The proposed methodology covers the following steps.

• Preparing the geometry of a porous body with the use of any method. In the hereby paper, virtual beds generated by the use of the Discrete Element Method (Section 2.2) are applied. They represent the existing bed samples consisting of glass marbles. Such a choice was dictated by the fact that the alternative results for these bed samples are available and thus can be used as comparative data. At this stage, it is assumed that the porosity of virtual bed samples must be the same as the porosity of a real bed sample. The formulas collected in Table 1 clearly indicate that porosity is the most important factor influencing the tortuosity value.
• Converting the available vector geometry to a binary geometry.
• Calculating the path lengths (Section 3.3) and then calculating the binary tortuosity (Section 3.4) with the use of the A-Star Algorithm (Section 2.3) and the Path Searching Algorithm (Section 2.4). At this stage, it is assumed that all calculations will be repeated for three geometries (Section 3.1) and for a few different resolutions of the binary geometry (Section 3.2).
• Comparing the obtained results with other independent data as well as with the data obtained from empirical formulas collected in Table 1 (Section 3.4). At this stage, the hydraulic tortuosity, calculated with the use of the Lattice Boltzmann Method (which is also based on the geometry expressed by binary matrices), as well as the geometric tortuosity, determined by the Path Tracking Method, serve as the reference data [42]. These tortuosities were obtained for the same material as in the present study (Section 2.1). Additionally, the results obtained by Wang, who investigated tortuosity in very similar granular systems [22], are taken into account during the comparison.

It should be highlighted, however, that the proposed methodology may be used only in the pore system with all channels open. If some pore channels are closed (blind cavities appear), then the path length will be much too long and the tortuosity will be strongly overestimated [21]. There are two kinds of systems in which the binary tortuosity may be calculated: granular beds (in 3D only, where a particle has contact points and the space around these contact points is always available) and any system (in 2D or 3D) with relatively high porosity consisting of separate solid objects. The limitation mentioned here is related to both algorithms used.

The main aim of the study is to find a method for calculating the tortuosity in high scale granular beds (which means the beds consisting of a high number of particles) within a reasonable calculation time. In other words, it is a try to find a good alternative for the approach proposed by Koponen et al. [17], which may be applied only in relatively small systems and which demands a huge computational power [21,42,43]. The hereby paper may be also treated as a continuation of investigations described in the three cited papers.

2. Materials and Methods

2.1. Materials

It was assumed that the numerical investigations are related to a real sample of granular bed consisting of SiLibeads Glass beads (Type S) [44] (Figure 2). To obtain the data related to particle sizes, 100 beads were randomly selected and then their diameters were measured, along two axes, with the use of the micrometre screw of $\pm 0.01$ [mm] accuracy. It was stated that the diameters follow Gaussian distribution at significance level 0.01 of the Kolmogorov–Smirnov test [45]. The mean diameter was equal to 6.072 [mm] with the standard deviation of 0.051 [mm]. During the next measurements, based on a graduated 200 [mm] cylinder, it was found out that the average porosity (obtained from 15 repetitions) was equal to 0.413 $[-]$. Such a value agrees with the literature data [46,47].

2.2. Discrete Element Method

The Discrete Element Method (DEM) is a method of analyzing various solid systems, which is based on the classic Newton’s dynamic principles [48]. Typically, two kinds of solids are considered: particles, most often representing a granular matter, and walls, representing in turn various parts of technical infrastructure (tanks, bins, silos, pipes, chutes, conveyors components, etc.). Originally, the particles could only have a spherical shape, but nowadays other simple shapes (such as cylinders, cuboids, polyhedrons, etc.) or complex shapes consisting of rigidly connected bodies with basic shapes (the so-called clumps) are increasingly used. Due to the article context, the discussion is limited to analyzing spherical particles only. However, the methodology described may be also used in beds consisting of particles with a more complicated shape.

In general, a progression movement equation and a rotation movement equation in DEM is considered. In the case of i-th spherical body, which is in contact with other bodies (Figure 3), the equations may be saved as follows,

$$m_i\frac{d{\overrightarrow{v}}_i}{dt}=\sum _{j=1}^{n_c}\left({\overrightarrow{F}}_{n,ji}+{\overrightarrow{F}}_{t,ji}\right)+{\overrightarrow{F}}_e$$

(6)

and

$$I_i\frac{d{\overrightarrow{\omega }}_i}{dt}=\sum _{j=1}^{n_c}\left(r_i\times {\overrightarrow{F}}_{t,ji}+f_{ij}\times {\overrightarrow{F}}_{n,ji}\right)$$

(7)

where m is the mass [kg], I is the moment of inertia [kg m${}^2]$, $\overrightarrow{v}$ is the linear velocity [m/s], $\overrightarrow{\omega }$ is the angular velocity $[1/s]$, t is time $\left[s\right]$, $n_c$ is the number of contacts $[-]$, ${\overrightarrow{F}}_n$ is the normal force $\left[N\right]$, ${\overrightarrow{F}}_t$ is the tangential force $\left[N\right]$, ${\overrightarrow{F}}_e$ is the external force $\left[N\right]$, r is the radius [m], and f is the distance between the direction of acting the normal force and the mass centre [m].

If the touching circular or spherical bodies are non-deformable, the contact zone is limited to a point, and the distance $f_{ij}=f_{ji}$ is equal to zero (Figure 3a). The distance between the centers of the particles i and j is equal to $r_i+r_j$. This type of contact is called hard contact in the literature. If the bodies are deformable (Figure 3b), instead of a point contact, a linear (in 2D) or surface (in 3D) contact will be created. In this case, the point of acting of the resulting normal force depends on the distribution of normal stress on this line or surface. Consequently, normal force will generate an additional moment of force. Because the particles are deformable, the distance between their centers is smaller than the sum of the radiuses $r_i$ and $r_j$. In this way, the degree of the overlap of particles becomes a measure of their deformation. The contact of deformable particles is called soft contact in the literature. Looking at Figure 3, it can be seen that there are two key problems in the Discrete Element Method. The problems are related to (a) searching for all contact points and (b) calculating forces and moments of forces at all contact points based on appropriate contact models. The original Cundall and Strack model (1979) [48,49] or the Hertz theory (1881) [50] is very often used to calculate the normal forces in contact points [51,52,53,54]. In the literature, other contact models may be found, e.g., Kelvin-Voigt (1890) [55], Hunt-Crossley (1975) [56], Kuwabara-Kono (1987) [57], Lankarani-Nikravesh (1990) [58], Tsuji et al. (1992) [59], Zabulionis et al. (2012) [60], and Singh et al. (2014) [61]. Two ranges are usually considered when determining the tangential force. If the tangential force is relatively small, then the tangential deformation is calculated. In this range, similar formulas are used as for the calculation of normal force. If the tangential force exceeds the static friction force, it slips with a constant dynamic friction force.

2.3. A-Star Algorithm

The A-Star (A*) Algorithm is an iterative algorithm serving for finding paths in a graph or a grid of nodes based on a specific weighted function. The aim is to find a path from a chosen Start Point to a given Stop Point at the smallest cost (i.e., the least distance traveled). Points may be freely located in the space or arranged in a form of a structured grid. A binary geometry may serve as a calculation space for the A-Star Algorithm. In such a case, the path may go only through the nodes with a value representing the pore space. Another limitation is that a path cannot go twice through the same grid node. The A-Star Algorithm was described in 1968 by Peter Hart, Nils Nilsson, and Bertram Raphael [62]. In general, the algorithm may be used not only to search distances, but also, for example, to calculate the minimum time of traveling. In the traditional approach, the algorithm is applied for 2D systems, but in the hereby study a 3D implementation has been developed.

The main idea of the A-Star Algorithm acting may be saved as follows,

$$f(x,y,z)=g(x,y,z)+h(x,y,z)$$

(8)

where $f(x,y,z)$ is a function serving to choose the next point of the graph or grid of nodes, $g(x,y,z)$ is the movement cost (the distance between the Start Point and the current node), and $h(x,y,z)$ is a heuristic function (the distance between the current node and the Stop Point).

In the regularly structured grids (i.e., the grids with a uniform distance between the nodes in each direction), the function $g(x,y,z)$ may have only three values: $D·1$ for nodes having two common indexes with the current Start Point (point A in Figure 4); $D·\sqrt{2}$ for nodes having one common index with the current Start Point (point B in Figure 4); $D·\sqrt{3}$ for nodes without common indexes with the current Start Point (point C in Figure 4). The current Start Point is marked in Figure 4 by black color. In regular structured grids the space scale coefficient D is usually equal to 10. In such a way, the values of $g(x,y,z)$ function are approximately equal to 10, 14, or 17. In fact, the exact values are not important because what really matters is the relation between particular values.

In the regularly structured grids, the function $h(x,y,z)$ may be calculated with the use of a few different heuristic functions. The most popular of them are as follows.

• Manhattan function

  $$h(x,y,z)=D·\left(\left|x_i-x_{stop}\right|+\left|y_i-y_{stop}\right|+\left|z_i-z_{stop}\right|\right)$$

  (9)
• Diagonal function

  $$h(x,y,z)=D·max\left(\left|x_i-x_{stop}\right|+\left|y_i-y_{stop}\right|+\left|z_i-z_{stop}\right|\right)$$

  (10)
• Euklides function

  $$h(x,y,z)=D·\sqrt{{\left(x_i-x_{stop}\right)}^2+{\left(y_i-y_{stop}\right)}^2+{\left(z_i-z_{stop}\right)}^2}$$

  (11)
• Quasi-Euklides function

  $$h(x,y,z)=D·\left({\left(x_i-x_{stop}\right)}^2+{\left(y_i-y_{stop}\right)}^2+{\left(z_i-z_{stop}\right)}^2\right)$$

  (12)

  where D is a space scale coefficient; $x_i$, $y_i$, and $z_i$ are coordinates of i-th point located adjacent to the Start Point; and $x_{stop}$, $y_{stop}$, and $z_{stop}$ are coordinates of the assumed Stop Point of the path.

After calculating the value of $f(x,y,z)$ function for every node neighboring to the current Start Point (there are 26 such points in 3D space, see Figure 4), the node with the minimum value of this function is finally selected as the next point of the path. In the next iteration, this new node is treated as the next Start Point.

2.4. Path Searching Algorithm

Path Searching Algorithm (PSA) is destined to search free passages through pore structures saved in a form of binary geometry. The Path Searching Algorithm has the same limitation as the A-Star Algorithm: (1) the path may go only through the nodes with the value representing the pore space and (2) the path cannot go twice through the same grid node. The algorithm may be used in 2D as well as in 3D domains. The algorithm was developed by author in 2019 [21]. It was used in order to filter data in a set of pore structures consisting of 15,000 cases (with random geometries without a free passage omitted). In the present study, the algorithm is employed to calculate the tortuosity.

In the opposition to the A-Star Algorithm, PSA has a solely local character and the locations of Stop Points are not needed for its acting. Other differences are as follows. (a) The path may change the direction only at right angles; (b) the path length increases in each iteration by 1, which simplifies the way of calculating the total paths lengths; and (c) in some specific cases the next node is drawn in a random process (which means that after running a search twice, both paths may have different shapes and lengths).

The acting of the algorithm may be summarized as follows (see Figure 5).

• If the movement forward is possible (it is a so-called free node), than the path goes straight in this direction (Figure 5a). A free node means here a node located in the pore space and not belonging to the current path. In Figure 5, the main direction in which the path is determined is from the left to the right.
• If the movement forward is impossible (the next node in the main direction is not free), then the path turns at a right angle (Figure 5b). The movement is possible only to a free node. If there is more than one free node (Figure 5c), then the movement direction is drawn by the use of the random number generator.
• If the movement forward is impossible, the path cannot go perpendicular to the main direction, but the node located behind the current node is free, then the path turns back (Figure 5d).
• If no further movement is possible (such a case is called as a cavity, Figure 5e), then the path points are deleted in sequence until another free node (not used before) is found.

3. Results and Discussion

3.1. DEM Simulations

To generate a set of virtual beds, the Discrete Element Method implemented in the YADE code [49,63] was used. The so-called three-axial compression technique was applied, in which in the first phase a random particle cloud inside a cuboid is generated. In the second phase, the cuboid walls are getting closer and the particle cloud is compressed until a target porosity equals to the porosity of real bed samples. The YADE code was used two times in the presented investigations. In the first stage, the relation between the number of particles and the porosity was investigated. In this study, uniform particles were used with the diameter equal to $6.072$ [mm]. The number of particles was in the range between 100 and 1000 and during each simulation, one million time steps was performed. The friction angle was set at 0 to facilitate the particle relative movement and to accelerate compression. The traditional Cundall’s linear elastic-plastic contact model in the calculations was used (class Law2_ScGeom_FrictPhys_CundallStrack [49]). It was assumed that the particles are made of glass and the walls of steel.

The results of the calculations are presented in Figure 6. It may be seen that ~1000 particles have to be assumed to obtain the porosity at the similar level as in the experiment. Due to the random character of the particle cloud creation, the resulting porosity varies in some range. It means that not every simulation will fulfill the condition of obtaining an assumed experimental porosity.

In the second stage of the investigations, the number of particles was set to 1000 and the simulation was continued until the target porosity (equal to 0.413) was reached. The other difference was that the particle size was defined by a cumulative discrete distribution function consisting of 25 bins. In such a way the average diameter and the standard deviation of samples described in Section 2.1 have been mapped. The details related to this issue may be found in the paper [64]. The simulation had to be repeated several times in order to obtain three examples of virtual beds with appropriate porosity. In Figure 7, the resulting virtual beds are visible.

3.2. Geometry Conversion

The virtual beds generated by the Discrete Element Method are expressed in a form of a vector geometry. It is done by saving coordinates of all particle centers and their diameters. It is important to remember that these particles cannot be treated as elements of a graph or as nodes of a non-structural grid for the A-Star Algorithm or the Path Searching Algorithm (tortuosity is calculated between solid objects and not inside them). It means that for the following purposes, the geometry of the pore space must be mapped in some way. To reach this aim, first, the space in which the virtual bed is located, is covered by a structural grid of nodes. Next, a binary value is assigned to each node: TRUE (or 1) if the node is located inside a sphere, and FALSE (or 0) if the node is belonging to the pore space.

As the resolution of grids of cells or nodes in numerical techniques is usually very important, it is probably this factor (for the appropriate data are not available in the literature), which may have an impact on the results in the methodology proposed. For this reason, four numerical node grids with the following resolutions of 50 × 50 × 50 (Figure 8), 100 × 100 × 100 (Figure 9), 150 × 150 × 150 (Figure 10) and 200 × 200 × 200 (not possible to visualize with the available computer) were prepared for every virtual bed. The total number of nodes was limited by the available computational power. It may be seen (Figure 8, right) that the smallest resolution does not provide circular cross sections of the spherical particles. The geometry of the virtual beds used is represented better by higher resolutions. The bed geometry visible in Figure 8, Figure 9 and Figure 10 is rotated in comparison to the original geometry obtained from the YADE code, but this does not matter due to the isotropic character of the sample.

3.3. Determination of Path Lengths

Virtual beds in the binary representation were used to calculate lengths of paths on the basis of the A-Star Algorithm and the Path Tracking Method. Both algorithms were self-implemented in the Fortran programming language. All results were saved in VTK file format, which is very convenient to visualize, e.g., in the ParaView software. The A-Star Algorithm was applied to both versions, i.e., a standard one (here called static) and a modified one (here called dynamic). In the static approach, the coordinates of the Stop Point in the plane perpendicular to the main direction (here X direction) are the same as the coordinates of the Start Point:

$$i_{stop}=n_x,j_{stop}=j_{start},k_{stop}=k_{start}$$

(13)

In the modified version, the Stop Point changes their locations during the calculation as follows,

$$i_{stop}=n_x,j_{stop}=j_i,k_{stop}=k_i$$

(14)

where index i means the current path node and $n_x$ the number of grid nodes in the X direction. The modification is proposed due the fact that, in the transport processes in porous media, the fluid trajectory depends mainly on the local conditions (geometry). In other words, it is not strongly determined by a point located somewhere on the outflow surface far down the flow.

In each case, 25 different Start Points regularly arranged in the inflow plane were used. Their locations were calculated as follows,

$$i_{start}=1,j_{start}=s{p}_j·int\left(\frac{n_y}{ns{p}_j+1}\right),k_{start}=s{p}_k·int\left(\frac{n_z}{ns{p}_k+1}\right)$$

(15)

where ${sp}_j$ is the ordinal number of Start Points in Y direction $[-]$, ${sp}_k$ is the ordinal number of Start Points in Z direction $[-]$, $n_y$ is the number of grid nodes in the Y direction $[-]$, ${nsp}_j$ is the total number of Start Points in Y direction $[-]$, $n_z$ is the number of grid nodes in the Z direction $[-]$, and ${nsp}_k$ is the total number of Start Points in Z direction $[-]$. The function $int\left(\right)$ means that the resulting value is rounded to the nearest integer number.

In Figure 11, an example of using the A-Star Algorithm in the static variant for grid resolution equal to 200 × 200 × 200 is presented. On the left, all of the obtained paths for the first geometry are visible. Path points are colored according to the order of adding. On the right, four chosen paths for the same Start Points and three geometries are compared. It should be noted that in the static variant all paths are relatively straight due to assumed “attraction” to the Stop Point. In Figure 12, analogical results, but this time for the A-Star Algorithm in dynamic variant, are visible. It may be seen that the effect of the Stop Point “attraction” does not occur here. In Figure 13, the results of acting of the Path Searching Algorithm are visible. Paths are not as smooth as in the A-Star Algorithm, which is caused by the limitation when it comes to choosing the next path points (a diagonal direction is prohibited).

In Figure 14, four chosen paths for the first virtual sample and different resolutions are visible. Paths obtained for resolution 50 × 50 × 50, 100 × 100 × 100, 150 × 150 × 150 and 200 × 200 × 200 are marked by black, green, blue, and red color, respectively. It may be seen that the resulting number of paths depends on the grid resolution. All Start Points located in the solid space are omitted. If the resolution is low, the probability that a Start Point will be located in a solid node is higher. This probability decreases when the grid resolution increases. This relationship is illustrated in Figure 15.

In Figure 16, the total number of path points divided by the number of grid nodes in the main direction is visible. This number depends on the grid resolution and the method applied. It may be observed that the results obtained for the 50 × 50 × 50 resolution differs in all cases from the other results. Therefore, it may be concluded that the quality of the binary representation of the real geometry is in this case too low. This conclusion is obvious however: the minimum value may depend on the approach used in the investigations. For example, in his study based on the Lattice Boltzmann Method (LBM), Wang stated [22], that the radius of sphere in a binary representation should have at least ten grid nodes and the media length along flow direction should be more than twenty radii. Looking at Figure 7, it may be seen that within the domain, 9 particles are usually located along an edge. It means that the number of nodes per sphere diameter is approximately equal to 6, 11, 17, and 22 for the resolutions of 50 × 50 × 50, 100 × 100 × 100, 150 × 150 × 150 and 200 × 200 × 200, respectively. It seems that if the A-Star Algorithm and the Path Searching Algorithm (which are qualitatively similar) are used, the resolution may be two times smaller then the one mentioned by Wang. It is an important information, because the demand for computing power increases linearly with the number of grid nodes. Lower resolution means faster acting and shorter solution times. The conclusion related to the minimum size of a sample is similar to this, which was reported by Wang. In all virtual bed samples used in the investigations, the domain length was higher than 18 radiuses, which is very close to the value mentioned by Wang (20 radiuses).

One problem with the LBM method is that it works incorrect for the channels narrower than approximately four lattice units [19]. It means that the binary geometry obtained in the current study should be refined at least four times to fulfill this condition. In such a way the resolution of the highest grid should be equal to 800 × 800 × 800 (512 ×${10}^6$ nodes) instead of 200 × 200 × 200 (8 ×${10}^6$ nodes). It is impossible to analyze such big data with a PC computer. In this context, it is very comfortable that the channel width in the proposed methodology may be equal to only one node and a grid refinement is not needed at all.

In Figure 17, the real calculation time via the function of grid nodes is presented. It turns out that the Path Searching Algorithm works approximately 2 times faster than the A-Star Algorithm. Moreover, no significant differences between A-Star Algorithm in static and dynamic versions are visible. It should be also noted that in general the calculation time is relatively short. In previous investigations, in which the Lattice Boltzmann Method to calculate the hydraulic tortuosity in a virtual bed consisting of 50 particles was used [43], the calculation time was much longer. For example, the calculation time for a grid with the resolution of 32 × 32 × 64 and 8000 iterations was equal to 8990 s. It is 2.36 times more than the longest calculation time in this study (equal to 3815 second for the A-Star Algorithm in static version, the virtual bed sample number 1 and the grid resolution of 200 × 200 × 200). The difference seems to be not very significant, but previously the binary grid had only 65,536 nodes, while here the number of nodes is equal to 8 ×${10}^6$ (and the number of particles is 20 times higher). The conclusion is that the application of algorithm, such as the A-Star or Path Searching Algorithm, to calculate the path lengths (and tortuosity) is very attractive in comparison with the Lattice Boltzmann Method.

3.4. Tortuosity Calculation

In the last stage of investigations, the tortuosity was calculated according to the Formula (3), wherein $L_0$ is equal to $n_x-1$. In Figure 18, all the obtained results are collected. On the left, the data for all resolutions are visible, whereas on the right, only these for the highest grid resolution are presented. In

Table 2. Relative errors between the results obtained in the investigation and data from the works in [22,42].

Method	$\mathit{\tau }$	${\mathit{\delta }}_{\mathbf{Wang}}$	${\mathit{\delta }}_{\mathbf{LBM}}$	${\mathit{\delta }}_{\mathbf{PTM}}$
	$[-]$	$[\%]$	$[\%]$	$[\%]$
A-Star static	1.1464	−4.27	−7.01	−5.53
A-Star dynamic	1.1259	−5.98	−8.67	−7.21
PSA	1.4220	18.75	15.35	17.19

, a comparison between the average values of tortuosity calculated for all the method applied in this investigation and results from the literature [22,42]. It turned out that in comparison with hydraulic tortuosity calculated on the basis of the Lattice Boltzmann Method as well as with geometric tortuosity calculated with the use of the Path Tracking Method (details in [42]), the A-Star Algorithm underestimates the result and the Path Searching Algorithm overestimates it. Additionally, the Path Searching Algorithm gives in all cases higher relative error. The differences between static and dynamic variants of the A-Star Algorithm are very small, however, the static variant leads to obtaining results with the smallest relative errors.

The situation differs if the obtained data are compared with the results of calculations based on empirical formulas collected in Table 1. In this comparison (visible in Figure 19), the Path Searching Algorithm provides results in the same range. It is visible once more that the smallest resolution does not have enough good quality.

It may be stated that the binary algorithms such as the A-Star Algorithm or The Path Searching Algorithm may be a good alternative for other techniques of calculating the tortuosity, particularly these based on the numerical methods. They are not complex in theory, are relatively easy to implement and use, fast, and do not generate big data files.

4. Summary

The following conclusions can be formulated based on the results of the present study.

• To obtain the same porosity as in a real bed sample, the virtual bed has to contain approximately 1000 particles.
• The resolution of the binary representation of a real (or virtual) bed cannot be too small. However, this resolution may be still much lower than the resolution of grids required by other numerical techniques, in particular by the Lattice Boltzmann Method (used by Wang [22] or in the previous investigations described in the [42]).
• It is possible to calculate tortuosity with algorithms working on a geometry saved in a form of binary matrices.
• The Path Searching Algorithm gives always higher values of the tortuosity than the A-Star Algorithm.
• The A-Star Algorithm (in both implemented variants) underestimated the tortuosity in comparison with the tortuosity (hydraulic or geometric) calculated by alternative numerical techniques.
• The Path Searching Algorithm overestimates the tortuosity in comparison with the tortuosity (hydraulic or geometric) calculated by alternative numerical techniques. However, the results are in good agreement with the values calculated with the use of empirical formulas collected in Table 1.
• The Path Searching Algorithm is about two times faster than the A-Star Algorithm.