The influence of particle size distribution on parameters characterizing the spatial structure of porous beds

The article proposes functions linking the standard deviation of a particle distribution in a porous bed consisting of spherical particles with various parameters characterising the spatial structure of the bed. The porosity, the inner surface, the specific surface and the geometrical tortuosity were analysed. In the first stage, a set of virtual beds was created with the use of the Discrete Element Method. The Radius Expansion Method was applied to generate virtual beds with different standard deviations. 150 virtual beds were created (25 standard deviations, 3 repetitions with different settings of the random number generator, 2 values of the radius expansion factor). In the second stage, the spatial structure of all virtual beds was analysed. The geometrical tortuosity was calculated with the use of the so-called Path Tracking Method; other parameters were calculated with the use of analytical formulas. The impact of the standard deviation on the parameters characterising the spatial structure of the granular bed was described by approximation functions, which can be used in order to obtain these parameters based on the particle size distribution for others porous beds.


Introduction
The prediction of the pressure drop, occurring during fluid flows through porous media, is one of the most important problems in the widely understood science and engineering. For example, this issue plays significant role in geology, civil engineering, agriculture, food industry, chemistry and many other areas. So-called porous beds, i.e. granular beds consisting of spherical or quasi-spherical particles, play an important role here-the article relates to such kind of the physical matter.
The pressure drop occurring during fluid flows through porous media in many cases can be described with the use of one general mathematical formula [43]: where p-pressure (Pa), x-a coordinate along which the pressure drop occurs (m), A( ) and B( )-two generalized parameters, dependent on the set of parameters characterizing the spatial structure of the porous medium,dynamic viscosity of the fluid [kg/(ms)], -density of the fluid (kg/m 3 ), v f -filtration velocity (m/s), T-temperature ( • C).
Equations proposed by Forchheimer [16], Ergun [14], Macdonald [26] and Kovács [21], as well as other similar equations, may serve as examples of mathematical formulas matched to the topology of the formula (1). If B( ) coefficient equals zero, formula (1) becomes Darcy-like topology. Besides Darcy [11], achievements of Hazen [19], Slichter [35], Terzhagi [44], Kozeny [22] or Carman [6] must be mentioned here. Detailed review of these formulas is not important in the context of this article (it can be found e.g. in [41]). The crucial thing is that all these formulas are functions of different parameters characterising the geometrical structure of the porous body. In the context of porous beds, such set of parameters (denoted here as ) may have the following form (Table 1) [37]: where d-the representative particle diameter (which may be defined in many ways [41]) (m), -the porosity (−), -the packing coefficient (−), e-the void ratio (−), V s -the volume of the solid part of the porous body (m 3 ), V p -the volume of the pore part of the porous body (m 3 ), V-the total volume of the porous body (m 3 ), -the tortuosity (−), L 0 -the depth of the porous bed (m), L p -the length of the pore channel (m), f -the tortuosity factor (−), S 0 -the specific surface of the solid body (in the Kozeny [22] or Carman meaning [6]) (1/m), S p -the inner surface of the porous body (m 2 ).
In the context of the described work it is important that the elements of the set are usually treated as constant values, representative for the whole porous bed. Such approach is used in all formulas cited in the description of the Eq. (1). However, these parameters may depend on other factors. In the investigations of granular beds consisting of spherical or quasi-spherical objects, the particle size distribution may be such factor. This is the issue that we study in this paper.
Most elements of the set may be calculated analytically, if only the number of particles in the bed and their sizes are known. The type of particle distribution is not directly considered in such calculations, but of course it is concealed in above-mentioned data. From the whole set, only the tortuosity has to be obtained in another way. In the literature, many formulas destined for calculation of this quantity may be found. In the [39] we checked 21 such relations and in most cases tortuosity is treated as a direct function of the porosity. An additional parameter (a shape factor) was introduced only in two cases. We did not found any paper investigating the relationships between the kind and parameters of the particle distribution and the tortuosity. It was the main motivation to perform research described in this article. Study like that has the potential to propose formulas allowing calculation of porous bed parameters based on variance of size of particles forming the investigated bed, with no need of experiments, that can be time-consuming, cost-intensive or simply impossible. It should be noted that in the literature investigations in which a particle distribution of a porous bed is taken into account may be found (e.g. [4,7,30,49]). However, the Authors focus usually on the general effect, like permeability or filtration coefficient (values of A( ) and B( ) terms), not on the elements of the set. Papers related to investigations of mono-and poly-disperse granular systems are known too, even related to their inner geometrical structure [28,29,33]. However, it is very difficult to compare the results of these investigations with our data on a quantitative level.
The need for the calculation of the tortuosity determines the methods used in the investigations. To reach the article aim, we use the so-called Path Tracking Method (PTM). This method needs the data on the locations and sizes of all particles forming a bed. Such data may be obtained experimentally (e.g. by the use of the computed tomograpfy and image analysis) or numerically (e.g. by the use of the Discrete Element Method, DEM). Due the fact that it is impossible to freely change the standard deviation in a real bed, in our study the second approach is used. In DEM simulations, the particle size distribution may be freely defined. To keep the resemblance to real granular beds, the porosity, the average particle diameter and the basic value of the standard deviation were taken from an earlier experiments [41].
The key question is whether the virtual beds generated in DEM simulations may be used to analyze the geometry of real granular beds. DEM is a very complicated technique and many different factors (the applied contact models, numerical schemes, time approach and many others) may influence the final spatial arrangement of the particles. To answer this question, all possible elements of the set were calculated in two ways: direct, based on the discrete form of the cumulative curve of particle distributions; and indirect, based on the data obtained from DEM simulations. This approach allows to show differences between these two ways.
The methodology of the presented investigations consists of four main steps: (a) preparation of a set of cumulative curves of particle distributions (with Gaussian distribution); (b) analytical calculations of all possible elements of the set in function of the standard deviation; (c) preparation of a set of virtual beds with the use of the Discrete Element Method (where the type and parameters of the particle distribution is specified); (d) use of our own algorithms destined for analysis Table 1 Parameters characterising the spatial structure of a porous bed [41] No.
Parameter Definition Inner surface

Materials
A set of virtual beds consisting of spherical particles and created with the use of the Discrete Element Method is the object of our investigations. We assume: (a) Gaussian distribution of the particle size; (b) one average particle diameter-equal to 6.072 (mm) ( d e ); (c) one basic value of the standard deviation equal to 0.05 (mm) ( e ); (d) one bed porosity (the so-called target porosity, e ) equal 0.413 (−). These values are taken from previous studies [41], in which real beds with such parameters were investigated.

Distributions of diameters
The main aim of this task is to create a set of cumulative curves of particle size distributions with a constant average diameter (equal to d e ) and different values of standard deviations ( ). Classic mathematical formulas for normal distribution are sufficient for generating theoretical, discrete distribution of the particles forming the bed. By the theoretical distribution, we mean the set of bins including numbers of spheres with a fixed diameter. Distribution is calculated in the range of ± 4 e from d e . In case of normal distribution, it means that ≈ 99.994 % of the population falls inside that range (15,787 must be the minimum total population to expect the first case falling outside that range [13]).
The procedure of generating the theoretical distribution is following. The range of ± 4 e is divided into n f (ie. number of fractions) equal subranges. Each i-th (where i = 1, 2, … , n f ) subrange has an attributed diameter: 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: Explanation for above procedure is as follows. The sum of calculated integrals seen in the formula (4) 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. On that basis, we can get the cumulative curve for each distribution and having the theoretical distribution of the particles forming the granular bed, we can create the virtual beds with attributed mean diameter and different variances.

Discrete Element Method
The main aim of this task is to obtain independent sets of virtual granular beds with one porosity, one average diameter and different values of the standard deviation. By the term "virtual bed" we mean here the data on position in the space (x, y and z coordinate) and the size (radius or diameter) of every particle forming the bed. The virtual beds were prepared with the use of the Discrete Element Method (DEM) proposed by Cundall and Strack [10]. This method is based on the analysis of dynamics of a set of solid bodies (e.g. spherical particles grouped in a cloud), carried out with the use of Newtonian laws. The main equations of the linear and the angular motion may be written as follows [1,12,15]: where m i -mass of the i-th body (kg), I i -moment of inertia of the i-th body (kg m 2 ), v i -linear velocity of the i-th body (m/s), -angular velocity of the i-th body (rad/s), n ijnormal forces between neighbouring bodies i and j (N), t ij -tangential forces between neighbouring bodies i and j (N), n c -number of contacts between i-th body and neighbouring bodies (−), ext i -external mass 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), r ij -torque factor associated with the rolling friction (N m).
To generate a virtual bed consisting of spherical particles, the so-called Radius Expansion Method (REM) implemented in the non-commercial YADE code was used [51]. This name denotes a specific way of generating a particle cloud in a DEM model. In this method, an initial cloud of small particles (without any contacts) with a given particle size distribution is randomly generated, and then all particles increase in size until the target porosity is reached [51]. The radiuses increases with the formula r t+dt = kr t , where k is the so-called radius expansion factor and t is time. The term "cloud" means here a set of objects (spherical in our model) located in the calculation domain. It should be emphasized that the random generation of objects in a cloud is applied very often [3,32,46,51]. The number of contacts between particles increases during calculations and it is fixed at a constant level after some time. The REM may be used to create clouds with particles having one specific diameter or different diameters. In the second case, a discrete cumulative curve of particle fractions must be defined (Sect. 2.2.1). A characteristic feature of the REM is that the volume of the space in which the initial cloud is generated, must be very precisely defined. To control this volume, an additional parameter ( l c ) has to be calculated in every single case. In other case, the final cloud of particles obtains the correct porosity as well as the features of the particle distribution, but the diameter values will be incorrect.
The virtual beds were prepared with an assumption that the porous body has a cuboid shape with the size of l c l c 2l c . The last direction is twice longer than the others to increase the L 0 dimension. This is the direction in which paths lengths ( L p ) are calculated to obtain the geometrical tortuosity. The volume of the porous body may be also calculated as follows: where l c is an characteristic dimension of the porous body [m]. This quantity may be obtained from the formula: where and Symbols from formula (9) have following denotations: V sum i -the total volume of all particles in the i-th bin (m 3 ), V i -the volume of a single particle belonging to the i-th bin (m 3 ), Y i -the fraction of particles in the i-th bin (−), d i -the diameter of particle belonging to the i-th bin (m). In our investigations we assume that the total number of spheres (N) is constant for all virtual beds.
The other possibility to create a virtual bed with a defined particle distribution is the use of a triaxial compression. In this approach, a cloud of particles is created (the particle sizes and distribution are correct in that moment) and then the cloud is compressed by the walls to obtain the target porosity. The particle distribution is introduced in the same way as in the previous method. The main drawback of this method is a long simulation time, counted in hours, not minutes, like in the REM method. Due to the fact that a big number of virtual beds had to be generated, we decided to use the REM approach. Since this method is often applied in many different studies [8,17,23,45,48,50], new conclusions related to the consequences of its use seems to be desirable.

Path Tracking Method
The main aim of this task is to obtain values of the geometrical tortuosity for every virtual bed prepared with the use of the Discrete Element Method. The geometrical tortuosity is calculated with the use of the Path Tracking Method (PTM) in a variant called Regular Grid Method (RGM) [36,42]. The other parameters are calculated analytically, but with the use of the same software, the so-called PathFinder code [31].
Path Tracking Method is an iterative method of determination of the length of a pore channel in a porous medium in the chosen space direction, which consists of analysis of the local structure of the pore space based on the information on the location and diameter of every particle forming the bed. The PTM method was developed by one of the authors in 2009 as a tool destined for investigating the spatial structure of granular beds. As a result, the numerical code called PathFinder was written, and it is currently freely available in the Internet [31]. The length of a pore channel ( L p ) in the PTM is determined between two parallel planes, based on the sum of the unitary lengths. In turn, these unitary lengths are calculated based on the so-called tetrahedral structures that establish the basis for the calculation algorithm ( Fig. 1). Tetrahedral structures are created based on the data on the location and diameter of each particle in the bed.
The Regular Grid Method defines a way of using the PTM algorithm [36].

Approximation functions
Six following models were fitted to data in each case ( Table 2). Symbols p 1 , p 2 , p 3 and p 4 are models parameters and y is the parameter of the porous bed being under consideration. Models were fitted to data using Levenberg-Marquardt algorithm [27]. Then, Akaike Information Criterion (AIC) was used to compare models [2,40]. Hence, formulas given below as linking porous bed parameters with variance of particle size distribution are those indicated by AIC as best describing analysed relationship. Moreover, residual analysis (whiteness test of residuals) was performed for each model indicated by AIC in order to check, whether chosen approximation function properly explains changes in analysed value [25].

Discrete form of particle size distribution
The set of standard deviations used in our investigations to generate cumulative curve of particle distribution was as follows: 0.0001, 0.001. 0.01, 0.10, 0.15, 0.20, ...,1.15. Each particle distribution was saved in a discrete form of 101 fractions ( n f = 101 ). The total number of particles (N) was set to 10,000 (−), what is less than the maximum safe value calculated in Sect. 2.2.1. Examples of such functions are visible in Fig. 2. The Y sum symbol denotes the cumulative fraction of all particles with a diameter equal or less than d i .

Direct analytical calculations
If the discrete function of the particle size distribution and the total number of particles are known, then the values of V s , V p , V, l c , S p , S 0,Carman , S 0,Kozeny as well as other quantities (see Eq. 2) for every standard deviation may be directly calculated. The mathematical formulas visible in Table 1 and Eqs. (7)-(10) were used. In this approach, firstly V s with the use of the Eq. (9) is calculated and next, taking into account the target porosity e , the volume of the pore part ( V p ) as well as the total volume of the bed (V) are determinated. After obtaining this data, the porosity is calculated. Knowing all diameters, the inner surface is also calculated. At the end, the fit functions were obtained for all parameters listed above. Results of the direct calculations are collected in Table 3. The "hat" symbol was added Fig. 1 Schemat of the tetrahedral structure used in the Path Tracking Method; the meaning of the abbreviations is as follows [38]: ISP-Initial Starting Point (a point in the bottom surface of the bed, from which the calculation starts), FSP-Final Stating Point (a point from which the path begins), GC-gravity center of the triangle formed by particles P 1 , P 2 and P 3 , IL-Ideal Location (a predicted point, in which the centre of the fourth particle forming the tetrahedral structure is located), RL-Real Location (a real point in which the centre of the particle P 4 is located)

Virtual beds
To prepare a set of virtual beds with the use of the Discrete Element Method, the following assumptions were made: (a) the average diameter is equal to d e ; (b) the standard deviation is changed in a range as wide as possible; (c) the total number of particles (N) is equal to 10000; (d) the target porosity is equal to e ; (e) every virtual bed is generated three times with the use of different parameters of the random function (the repetitions are identified by numbers 1, 2 and 3); (f) all test are performed two times with different values of the radius expansion factor. When the radius expansion factor is less, then the time of DEM simulation is longer but the target porosity may be obtained more precisely. These tests are identified as A and B. The radius expansion factor was equal to 0.0001 and 0.001 in tests A and B, respectively. The visco-elastic contact model was applied. The damping coefficient was set to 0.2. Mass forces in each direction were set to zero due to the desire to obtain more homogeneous beds. In every simulation the friction angle (denoted by ) was firstly set to 0.5 [rad], and next it was decreased during calculations. In Fig. 3, examples of virtual beds from DEM simulations are visible. It should be stressed that when the Radius Expansion Method is used, the final particle distribution differs from the normal distribution (Fig. 4). Based on additional experiment, we noticed that this effect is independent of initial value of the friction angle . The n p symbol in this figure denotes the number of particles with a specific diameter. The above mentioned effect occurs despite the fact that the initial cumulative curves are symmetric (see Fig. 2) and they all were created for the normal distributions. The reason is that in many cases particles are blocked by neighbouring particles and cannot grow evenly during simulation. This effect is significant mainly for particles with medium and higher diameters. In turn, to obtain the target porosity, the particles with smaller diameters have to grow more than they are supposed to. Confirmation of this phenomenon is visible in Fig. 5, where higher normal stress magnitude is noticeable for higher diameters (The mean normal stress was calculated using the formula given by [18] and [47]). This trend is not very strong but it is visible in every case. It should be noted that for some particles with small diameters  the normal stress magnitude is equal to zero. It indicates that no normal forces affected these particles (they may have no contact with other particles). Such effect is not visible for the biggest particles. It also happens that big single particles, which have space to growth, increases excessively in size (see the right end of the chart visible in Fig. 4). The other consequence of uneven growth of particles is the decrease of the average diameter when the standard deviation increases (Fig. 6). However, the relative error (where d e is taken as the precise value) is small and do not exceed 2%. On this stage we do not arbitrate whether the virtual beds can form the basis for further investigations of the geometry of such kind of porous media. This issue is discussed in the next section.
It should be emphasized that the investigations of virtual beds generated rendomly with the use of the REM and the cumulative curves are not new in the literature [5,9,20,24]. However, in the available literature we did not found comments concerning the changes in the particle distribution if this approach to modelling is applied.

Indirect analytical calculations
If the data on the location and size of every particle forming the virtual bed is known, then the PathFinder code destined to calculate all elements of the set may be used. In this section we focus on these elements of this set which may be calculated analitically. Additionally, we discuss some issues related to the data obtained with the use of the Discrete Element Method.
In Fig. 7, changes in the porosity in function of standard deviation are visible. Important is that these values were collected directly during DEM simulation in the main calculation loop (similar relationship is shown in the article further ( Fig. 9), but obtained in the PathFinder code). This time, at the beginning of every DEM simulation, a calculation domain with a specific volume (V) dependent on the parameter l c is created. The parameter l c is calculated analytically, in the same way as in the direct approach. The total volume of the particles ( V s ) and the volume of the pore space ( V p ) is calculated in every time step on the basis of the current data to obtain the current porosity. The main loop stops when the current porosity is equal or less than the target porosity ( e ). This explains why the values visible in Fig. 7 are always a bit lower than this target porosity. As it can be seen, the porosity value is much more constant in the A test. Besides, the porosity decreases in A test if the standard deviation increases, otherwise than in the B test, where the porosity changes are more non-uniform. The minimum and maximum relative errors (where the target porosity is used as the reference value) in A test were equal to 0.06% and 0.68%, in B test in turn 0.38% and 2.63%. None of considered models was able to describe changes in porosity for the B test. In case of the A test, changes in porosity were best explained with the use of the second degree polynomial model: Our results showed that the porosity slightly decreases when the standard deviation increases. It coincides with information from the literature (e.g. [34]), that the packing fraction in the polydisperse beds is higher than in the monodisperse system. The higher packing fraction denotes the lower porosity. Similar information is given in [47], where the influence of number of granulometric fractions on structure and micromechanics of compressed granular packings was investigated.
During investigations it turned out that the characteristic dimension increases slightly with increase of the standard deviation (Fig. 8). The porosity changes are caused by differences in the local arrangement of particles in the bed. Increase of a variance of particle sizes causes small changes in the structure of contacts between particles and in consequence in structure of forces and momentums. In case of characteristic dimension, the results were best explained with the use of the second degree polynomial model in a form of: In Fig. 9, the porosity in function of standard deviation for data obtained from the PathFinder code is presented. In opposition to porosity obtained in the YADE code, the values from the PathFinder code are slightly higher. The procedure of calculation the porosity in the YADE code was explained earlier. The bed size (volume V) in the PathFinder is calculated based on the data of locations and sizes of all particles (details in [38]). The volume of the solid part ( V s ) is calculated by the use of Eq. (9). The (12) l c = 0.0998 + 0.00008 + 0.0025 2 .
minimum and maximum relative errors (where the target porosity is used as the reference value) in A test were equal to 0.23% and 0.65% , in B test in turn 0.19% and 1.82% . Important is that all calculations related to the set are performed in the PathFinder code, thus the values presented earlier in Fig. 7 are not taken into account. Similarly like in a case of relationship between porosity and standard deviation obtained for data obtained through DEM simulations (shown in Fig. 7), none of considered models was able to properly describe changes in porosity for the B test, where obtained results are clearly nonlinear. In case of the A test, changes in porosity are best explained with the use of the second degree polynomial model:   The results visible in Figs. 7 and 9 show that virtual beds obtained in the A set have higher quality. By the term "quality" we mean here the level of the data scattering. It leads to the conclusion that the radius expansion factor should be as small as possible. Therefore, we assume, that results obtained in A test are good enough to reach the aim of the article. In the next parts, the approximation functions are shown only for data obtained in the A test. However, to continue the discussion results for the B test are also shown.
The porosity is dependent of the characteristic volumes of the bed and its parts. In Figs. 10, 11 and 12, the impact of standard deviation on the solid part of the bed, the pore part of the bed and the total volume of the bed are shown. Approximation functions for these relationships are as follows: The changes in the characteristic volumes of the bed and its parts shows that the structure of the pore space changes slightly in function of the standard deviation. The key role plays the volume of the solid part of the porous body (Eq. 9), which is dependent on the particle distribution in the bed. Such parameters like V p , V and l c depend directly on this volume.
In Fig. 13, the inner surface for all defined data sets is visible in function of the standard deviation. Important is that the inner surface was calculated with an assumption that all particles are non-deformable, what means that they have point contacts [see the formula no. 6 in Table 1  (16) V A = 0.00199 + 0.00016 2 , standard deviation due the fact that the relationship between the diameter and the inner surface is quadratic. Thus, the bigger particles have more significant impact on the total inner surface than the smallest particles. In Fig. 13 a clearly trend is visible with a maximum for ≈ 1.05 . This disortion results from the way used to generating the virtual bed and don't should appears. This is the only one case, for which the indirect analysis gives different shape of the fit function than in the direct calculations. The approximation function for specific surface is the third degree polynomial of following form: In Figs. 14 and 15, the specific surfaces in meanings of Kozeny Fig. 12 Changes in the total volume in function of the standard deviation functions for these figures are given in following exponential form: In Table 4 the mean and maximum relative errors between fit functions obtained from the indirect analysis in relation to the fit functions from the direct calculations are shown.
As it can be seen, these errors are very small in all cases. On this basis we can assume that virtual beds created with the Radius Expansion Method properly represent the real beds and so they may be used to further investigations of their geometrical features.

Geometrical tortuosity
The values of the geometrical tortuosity were obtained with the use of the Path Tracking Method [39] in the variant called the Regular Grid Method. The bottom surface of the bed was divided into a regular grid of Initial Starting Points ( 100 × 100 ), from which 10000 paths lengths for every virtual bed were calculated. In Fig. 16, the geometrical tortuosity maps for two chosen cases were shown. One can see that the geometrical tortuosity value differs depending on the location of the point, in which the path begins (details in [39]). It is characteristic that small areas in which the geometrical tortuosity has the same value are clearly visible. This is due to the fact that paths starting at different points lying close to each other lead often in whole or in part through the same pore channel. This effect was described in [36]. In Fig. 17, average values of geometrical tortuosity for three sets of virtual beds and both A and B tests are visible. The approximation function for A set is as follows:   It should be stressed here, that in case of geometrical tortuosity, obtained approximation function is proper for the B set as well, so it does not depend on bed quality. It leads to the conclusion that details related to the metod used to generate a virtual bed may have relatively little meaning in calculation this quantity. It is a very desirable conclusion due to the fact that from the set, the tortuosity is the most difficult parameter to obtain. Note, that if the standard deviation is equal to zero, then the tortuosity equals 1.208. We may compare this value with data shown by Wang [46], who used coupled DEM-LBM (Lattice Boltzmann Method) simulations to calculate the so-called hydraulic tortuosity in granular beds consisting (20) A = 1.208 + 0.006e 1.54 .
of spherical particles. The hydraulic tortuosity obtained by him was equal to 1.1975 for the porosity equal to 0.4. Such compliance is very satisfying, especially since in the both cases the used methodologies were different.
In Fig. 18 geometrical tortuosity maps of two chosen standard deviations and three sets of virtual beds are shown. These maps are different in every case what confirms that the trend visible in Fig. 17 is independent on the arrangement of local values of the tortuosity.
It was mentioned that in the literature the tortuosity is usually defined as a non-linear function of the porosity. Here we want to check, whether the trend visible in Fig. 17 is not caused by changes of the porosity shown in Fig. 9. For this purpose the available analytical formulas for calculating the tortuosity in granular beds were used. A review of such formulas is available in [39]. In Fig. 19 it can be seen that such analytical formulas do not respond to small changes of the porosity reported earlier (Fig. 9). It confirms that the trend visible in Fig. 17 is a general feature, dependent on the standard deviation and not on the porosity fluctuations. Moreover, since the porosity decreases slightly in function of the standard deviation, the changes in the tortuosity may by even bigger than it can result from the formula (20).
The performed analysis shows that the impact of the standard deviation on the geometrical tortuosity cannot be assessed on the basis of analytical formulas, in which = f ( ).

Sensitivity analysis
To obtain more information about the influence of the standard deviation on the particular elements of the set, the

Conclusions
The following conclusions can be formulated based on the above-discussed topics: (1) Geometrical parameters characterising granular beds (identified in the paper as the set ) may be treated as functions of the standard deviation of the particle distribution. These functions may be obtained on the basis of Discrete Element Method and statististical methods; (2) The obtained functions may be used as replacements of constant values in such formulas like e.g. Kozeny-Carman equation, in which the porosity, the tortuosity factor and the specific surface have to be stated; (3) The functions proposed in the paper are developed on the basis of virtual beds and may be treated only as an estimation. However, we hope that the obtained trends and intensivities of changes are close to what is in fact, at least in relation to granular beds consisting of spherical or quasi-spherical particle; (4) Due the fact, that the DEM is an increasingly popular method of modeling different granular systems, the investigations of features, possibilities and limitations of this approach is very important and in our opinion fully justified; (5) A drawback of the Radius Expansion Method is distortion of an originally given particle distribution during the simulations (not commented in the available literature). However, this method seems to be sufficient to perform investigations related to the spatial structure of a porous bed; (6) The quality of a virtual bed depends on the radius expansion factor. This parameter should be relatively small; (7) The sensitivity analysis performed in the article allows estimating the intensity and the direction of changes of all elements of the set when the standard deviation increases; (8) Geometrical tortuosity in granular beds increases slightly and nonlinearly with an increase in the standard deviation of the particle size distribution. Additionally, this relationship seem to be independent on the quality of the virtual bed.

Compliance with ethical standards
Funding This study was funded by Polish Ministry of Science and Higher Education in the frames of the statutory research.  Fig. 19 Relationship between the porosity and tortuosity calculated with the use of different analytical formulas [39] as well as the Eq. (20)