Next Article in Journal
Investigating Machine Learning Techniques for Predicting the Process Characteristics of Stencil Printing
Previous Article in Journal
Fe/Mg/Fe Multilayer Composite Sheet Fabricated by Roll Cladding
Previous Article in Special Issue
Preparation and Characterization of Porous Materials from Pineapple Peel at Elevated Pyrolysis Temperatures
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Stochastic Filling and Modeling Algorithm of Non-Equal Diameter Particles with Specified Probability Density for Porous Permeable Materials

1
School of Mechanical and Electrical Engineering, Xi’an University of Architecture and Technology, Xi’an 710055, China
2
State Key Laboratory of Metal Extrusion and Forging Equipment Technology, China National Heavy Machinery Research Institute Co., Ltd., Xi’an 710018, China
*
Author to whom correspondence should be addressed.
Materials 2022, 15(14), 4733; https://doi.org/10.3390/ma15144733
Submission received: 13 May 2022 / Revised: 16 June 2022 / Accepted: 3 July 2022 / Published: 6 July 2022
(This article belongs to the Special Issue Feature Collection on Porous Materials)

Abstract

:
In this paper, a model generation algorithm for non-equal diameter particles with a specified probability density distribution is proposed. Based on considering the randomness of the size and distribution of the particles, the compact stacking of the particles is realized by the compactness algorithm, and then the spatial distribution of the tightly compacted particles is made to meet the random distribution of the specified probability density and the specified volume fraction by the filtering algorithm. The computational efficiency and effectiveness of the algorithm are verified, and the effects of the particle size and volume fraction on the distribution are analyzed. Finally, the proposed model has been used to study the permeability of a titanium porous filter cartridge. The results show that the size and location of the particle samples that are generated by the proposed algorithm follow specified probability distributions according to the requirements, and the volume fraction can be adjusted. Compared with the traditional algorithm, the computational effort and complexity are reduced. The resultant model can be used to study the permeability of porous materials and provide modeling support for structural optimization and further simulation of porous materials.

1. Introduction

Many new functional materials are bonded and compressed by granular materials with a pore structure that follows a certain probability distribution, such as titanium porous filter cartridge and metal-matrix composites. In order to obtain a mesoscopic mechanism for materials with a specific random distribution of particles, it is necessary to study the relationship between the properties of granular materials and the microstructure morphological characteristics of the formed materials, so as to carry out an effective structural performance design of the material [1,2]. However, because of the complex mesoscopic pore structure, it is difficult to realize the visualization of mesoscopic mechanism, for example, the seepage properties, through experiments [3]. Therefore, a meso-model with a specified probability density distribution is built to provide reference for improving the accuracy and efficiency of porous media meso-modeling.
As for the generation algorithm of particle distribution, Primera et al. [4] proposed the method of triangulation to describe the pore size, but this method was completed through two-dimensional slices and could not fully represent the three-dimensional morphology of the pores. Tory E M et al. [5] continuously filled the fixed space with the sequence addition method and completed the sequence filling of the powder layer based on considering the stability of the powder particles and the interaction force between the particles. Jodrey W S et al. [6] used the non-sequence rearrangement algorithm to start position allocation with large-diameter particles first, allocating the position with small-diameter particles, and then reducing the overlap degree by moving particles. Cundall P A et al. [7] made the stacking process of particles similar to the flow process of dynamic fluid based on the discrete element method and realized the stacking of powder particles. A. Bezrukov et al. [8] describes two algorithms for the generation of random packings of spheres with arbitrary diameter distribution. The first algorithm is the force-biased algorithm. It produces isotropic packings of a very high density. The second algorithm is the Jodrey–Tory sedimentation algorithm, which simulates the successive packing of a container with spheres following gravitation. It yields packings of a lower density and of weak anisotropy. Although the above algorithms could achieve the stacking of powder particles, there existed some shortcomings in the randomness of spatial distribution, and the filling efficiency was low. At the same time, the above algorithms inevitably need to constantly detect the interference between particles, which means tedious calculation.
In addition, Bailakanavar M et al. [9] used the random order adsorption algorithm to determine whether the current particles interfered with the existing particles, retained the data of the non-interfered particles, and then established the relevant model. Xin Zhenyang et al. [10] used the perturbation algorithm to assign regular distribution positions to all particles and applied the perturbation to the amount of random distance between each particle. Through interference judgment, the position information of particles was retained. Note that through the above, the generation algorithms could ensure the randomness of the space to a certain extent. Due to the need to compare the position relationship between the specified particles and all the particles in the distribution process, the computation burden is somewhat high, which affects the efficiency of the method.
The transmission characteristics of spherical packing include electrical conductivity and permeability, among which permeability is one of the most studied properties [11,12,13]. Experimental and theoretical studies on permeability have achieved certain results, such as Darcy’s Rule for fluid flow in porous media and the empirical formula Kozeny relation for permeability [14]. As the pore structure of sphere accumulation is complex like that of other porous media, it is difficult to study fluid flow in porous media by traditional methods [15].
To solve the above problems, this paper proposes a method to generate the meso-model of non-equal diameter particles with a specified probability density distribution. The filling method of particles avoids the tedious calculation that is caused by interference detection, and since the movement of particles is only related to a small number of particles in the neighborhood, the modeling efficiency is further improved to a significant degree. Through the comparison of computation with the traditional methods and a hypothesis test, the efficiency and effectiveness of the algorithm are verified, illustrating that the proposed generation algorithm could effectively build a meso-model for porous media and titanium foam with a specified probability density distribution. The lattice Boltzmann method is used to calculate the permeability of porous material with different probability density distributions. The results show that the porous material meso-model that is established by the proposed algorithm in this paper can provide a basis for further numerical simulation of fluid permeability.

2. Generation Algorithm for Meso-Distribution Particles

Considering the randomness of the particles in size and distribution, Python is employed to build the algorithm, and gives the generation steps of the particles a meso-model, according to the specified probability density. The proposed algorithm consists of two parts. The first part is the compactness algorithm, which causes the particles with the size of a specified probability density distribution to stack closely together. The second part is the filtering algorithm. By using this method, the distribution of the particles in the spatial position according to the specified probability density can be realized, and the specified volume fraction can be achieved based on the requirements.

2.1. Compactness Algorithm for Particles

Let   i   represent the direction of the   x axis, j represent the direction of the y axis, and k represent the direction of the z axis, respectively, corresponding to columns, rows and layers. The particles of random sizes are generated in the cube area according to the direction of the column, row and layer successively, and the tight compactness is realized by adjusting the position between the particles.

2.1.1. Formation of Non-Equal Diameter Particles

According to the actual demand, the particle size is randomly distributed according to the specified probability density by defining an appropriate sphere diameter range.
The expectation of X _ r is
E ( X _ r ) = A μ f d μ
In practice, the algorithm generates particles with a radius that is randomly distributed on the interval [a, b]; the average of the radius is
r a v g = E ( X _ r ) = a b x f ( x ) d x

2.1.2. Determine Initial Position

At this step, the algorithm places the random sized particles roughly in a 3D matrix pattern, with particle positions vibrated by random noises. To place the first particle with radius r 1   near coordinates (0, 0, 0), the algorithm considers the particle radius, the cube area boundary restriction and the random noise. The resulting center coordinates of the particle are noted as ( x 1 ,   y 1 ,   z 1 ) , where x 1 = r 1 + d n o i s e ,  y 1 = r 1 + d n o i s e   and   z 1 = r 1 + d n o i s e ; the second particle with radius r 2 is placed at coordinates ( x 2 ,   y 2 ,   z 2 ) , where x 2 = x 1 + r 1 + r 2 + d n o i s e , y 2 = r 2 + d n o i s e , z 2 = r 2 + d n o i s e . the third particle with radius r 3 is placed at ( x 2 + r 2 + r 3 + d n o i s e ,   r 3 + d n o i s e ,   r 3 + d n o i s e ) , etc. The notation d n o i s e stands for the random noise that the algorithm introduces; it takes a different value each time it appears. The particle placement process continues until the first row on the x y -plane is filled. Next, the algorithm places the second row close to the first row on the x y -plane and the row placement process repeats until the first x y -layer is filled. Afterwards, the second particle layer is placed closely above the first layer and the layer placement process repeats until there are enough layers. Prior to the placement of any particles, the algorithm estimates the number of particles to generate by considering the scenario where the cube area is filled with particles of equal radius r a v g . Assume the length of the edges of the cube along the x , y , z -axis are   L x , L y and L z , respectively; the number of particles in each row, the number of rows in each layer and the number of layers are estimated as N x = L x / ( 2 r a v g + ρ ) ,   N y = L y / ( 2 r a v g + ρ ) and N z = L z / ( 2 r a v g + ρ ) , respectively. The notation ρ stands for the upper bound of the gap distance between adjacent particles and is enforced by the compactness step (Section 2.1.3). Since the density of the particles is increased by the compactness step and the goal is to distribute the particles in the entire cube area following a given probability density, the estimated values N x ,   N y and N z are further scaled up by an overfill factor f o v > 1. After the compactness step, the excessive particles that remain outside of the cube boundaries are easily pruned.
To ensure that each newly generated particle and the last generated particle are adjacent in the same row and layer, the adjacent scales vibrate within a very small range. The advantage of this method is that the detection of particle interference can be omitted, and thus the randomness of the position can be retained, which significantly simplifies the computational effort of the algorithm. The logical steps are as follows:
(1) By specifying a minimal disturbance i n i t _ n o i s e to add a segment of white noise, its power spectral density will be constant in the whole frequency domain, as shown in Figure 1 (when i = 1 , the x coordinate of the center is evaluated on [ x _ m i n , x _ m i n + i n i t _ n o i s e ]   in a uniformly distributed scale, where x_min is simply the radius of the particle and init_noise is a preselected constant to control the amplitude of the noise). Otherwise, make the minimum coordinate x _ m i n of the newly generated particles in row j   of the k t h layer equal to the p r e v i o u s   m a x i m u m   x coordinate p r e v _ x _ m a x , plus the radius of the new particle, before taking the corresponding disturbance distance.
(2) Similarly, when   j = 1 , the y coordinate of the center is evaluated at a uniformly distributed scale on [ y _ m i n , y _ m i n + i n i t _ n o i s e ] , where y_min is simply the particle radius. Otherwise, since the ( j 1 ) th   row particles in layer k   have been determined, the minimum coordinate y _ m i n of the newly generated particles in layer k should be equal to the maximum y coordinate p r e v _ y _ m a x   of the ( j 1 ) th   row, plus the radius of the new particle, before taking the corresponding disturbance distance.
(3) Similarly, when k = 1 , the z coordinate of the center is evaluated on [ z _ m i n , z _ m i n + i n i t _ n o i s e ] , where z_min is simply the radius of the particle. Otherwise, since the particles at the   ( k 1 ) th layer have been determined, the minimum coordinate z _ m i n of the newly generated particles should be equal to the maximum z coordinate p r e v _ z _ m a x   ( k 1 ) t h layer, plus the radius of the new particle, before taking the corresponding disturbance distance.
If the amplitude of the noise is 0.1 and the mean value is 0.05, when the new particle has been generated 1000 times, the graph of the noise can be obtained, as shown in Figure 2.
The pseudocode of the placement process is as follows in Algorithm 1.
Algorithm 1 Algorithm of initial position determination
for k in [1, fov*Nz]
  z_max = −∞
  for j in [1, fov*Ny]
    y_max = −∞
    prev_x_max = −∞
    for i in [1, fov*Nx]
      generate particle with random size r
      if i = 1 then
        x_min = r
      else
        x_min = prev_x_max + r
      if j = 1 then
        y_min = r
      else
        y_min = prev_y_max + r
      if k = 1 then
        z_min = r
      else
        z_min = prev_z_max + r
      x_new = u(x_min, x_min+init_noise)
      y_new = u(y_min, y_min+init_noise)
      y_new = u(z_min, z_min+init_noise)
      place the newly generated particle at (x_new, y_new, z_new)
      prev_x_max = max(prev_x_max, x_new+r)
      y_max = max(y_max, y_new+r)
      z_max = max(z_max, z_new+r)
    prev_y_max = y_max
  prev_z_max = z_max
In the pseudo code, u (a, b) generates a uniformly distributed random value on the interval [a, b].

2.1.3. Compactness of Particles

Assume that the randomly generated particles are stored in a 3-dimensional array M , whose size is I m a x × J m a x × K m a x , n n e i g h b o r is a small positive integer, and ρ is a small positive value. The position of the particles will be repositioned by coordinate translation.
By calculating the distance d z between the particle and the n n e i g h b o r neighboring particles in layer ( k 1 ) , local neighborhood compression can be achieved to speed up the iterative process and more tightly compact the particles. While d y is the distance between the particles and the particles in the n n e i g h b o r range that are located in row ( j 1 )   of layer k, d x is the distance between the particle and the particles in the n n e i g h b o r range that are located in row j, column ( i 1 )   of layer k . Let d w = M a x   ( d x ,   d y , d z ) ( w   x ,   y ,   z ) . Move the particle by d w / 3 in the w direction. Repeat this step until d x ,   d y , d z < ρ .
The center coordinate of the particle P = ( x , y , z ) is reduced by 1 / 3 d w , d w = M a x   ( d x ,   d y , d z ) ( w   x , y , z ) , and a new position P   = ( x , y   , z   )   is obtained. The coordinates of the new position are derived as follows.
[ x y z 1 ] = [ 1 0 0 0 1 0 0 0 0 0 1 0 x y z 1 ] · 1 3 [ d w _ x d w _ y d w _ z 3 ]
According to the aforementioned methods, the particles at the initial position can be further compressed in distance to be tightly compacted.
The pseudocode of the compactness process is as follows in Algorithm 2.
Algorithm 2 Compactness algorithm
fork in [0, Kmax−1]
  for j in [0, Jmax−1]
    for i in [0, Imax−1]
      xboundary = max (x|x is the x coordinate of points on the particles at M [i−1,j,k])
      yboundary = max (y|y is the y coordinate of points on the particles at M [i,jn,k],
           where jn in [jnneighbor, j+nneighbor])
      zboundary = max (z|z is the z coordinate of points on the particles at M [i,jn,kn],
           where jn in [jnneighbor, j+nneighbor] and kn in [knneighbor, k+nneighbor]
   Assume:
     Plane Px parallel to the yz-plane intersects the x-axis at xboundary.
     Plane Py parallel to the xz-plane intersects the y-axis at yboundary.
     Plane Pz parallel to the xy-plane intersects the z-axis at zboundary.
       dx = distance between M [i,j,k] and Px
       dy = distance between M [i,j,k] and Py
       dz = distance between M [i,j,k] and Pz
        while max (dx, dy, dz) ⩾ ρ
          if dx = max (dx, dy, dz)
            subtract dx/3 from the x coordinate of M [i,j,k]
          else if dy = max (dx, dy, dz)
            subtract dy/3 from the y coordinate of M [i,j,k]
          else dz = max (dx, dy, dz)
            subtract dz/3 from the z coordinate of M [i,j,k]
          Recalculate dx,dy,dz
          if loop has executed more than 1000 times:
          Report exception and quit
Figure 3 is the compactness model diagram that is formed after the implementation of the compactness algorithm. The particles whose size complies with the specified distribution law are packed tightly in the cube region.

2.2. Filtering Algorithm

After the compactness algorithm, the particles with a random size distribution can be tightly compacted and a higher volume fraction of bulk density can be obtained. However, the spatial distribution cannot meet the random distribution of the assigned probability density, and the volume fraction cannot be adjusted according to the demand. Filtering algorithms help with this process. This paper makes the following assumptions:
(1) Assuming that the volume of the cube in which the particles are located is v , such that the particles are uniformly distributed, the probability density is the same throughout the cube;
  u x ,   y ,   z = 1 / v
(2) After the execution of the compactness algorithm and before the execution of the filtering algorithm, the total volume of the particles is v g ;
(3) If the target volume fraction is α , the total volume of the particles corresponding to this volume fraction is α v ;
(4) The target probability density function is f x , y , z .
Given that a small zone with volume d v is centered at ( x , y , z ) in the cube, the probability of retaining the particles near   ( x , y , z ) by the filtering algorithm is
s x , y , z = f x , y , z α   v d v u x , y , z v g d v
Substitute the hypothesis, simplify it, and we can achieve
s x , y , z = α   v 2 v g f x , y , z
Assuming that the tightly packed particles are stored in a three-dimensional array I m a x × J m a x × K m a x , after using the aforementioned methods, the particles that do not conform to the given probability density distribution in space are filtered out.
The pseudocode is as follows in Algorithm 3.
Algorithm 3 Filtering algorithm
fork in [0, Kmax1]
for j in [0, Jmax1]
  for I in [0, Imax1]
    Let x,y,z be the coordinates of the object M [i,j,k]
    if s > α v2/vg∙∙fx,y,z
    Delete M [i,j,k]
Figure 4 shows examples of common probability density functions (uniform, normal and exponential), as generated by the filtering algorithm.

3. Influence of Parameters on the Algorithm

The parameters of the algorithm directly affect its performance. In this section, the influence of the design parameters on particle filling and distribution effect in the application of this algorithm are discussed. The design parameters mainly include the particle size and the volume fraction. This section presents the distribution and filling effect and summarizes the related influence.

3.1. Distribution of Particles with Different Size Ranges with Same Probability Density

Taking the uniform distribution in the x ,   y and   z directions as an example, the distribution of the particles with different size ranges is simulated by using the proposed algorithm. The volume fraction is 5%, while the volume of the cube is 125 mm. The total volume of the particles is equal. The distribution results are shown in Table 1.
As can be seen from Table 1, when the volume fraction remains unchanged, the number of particles decreases significantly with the increase in particle size, and the distribution becomes sparser.

3.2. Particles Distribution with Same Probability Density and Different Volume Fraction

In this section, the uniform distribution in the x ,   y and z directions is taken as an example to simulate the distribution of particles in different volume fraction ranges by using the proposed algorithm. The radius sizes of the particles are all between 0.15–0.25 mm. The total volume of the cube is 125 mm3. The distribution results are shown in Table 2.
As can be seen from Table 2, when the radius size of the particles remains unchanged, the volume of the particles increases with the increase in volume fraction, but the maximum volume fraction of the particles cannot exceed the volume fraction of the tightest compactness that is formed by the compactness algorithm.

4. Algorithm Analysis

The effectiveness of the algorithm is analyzed in two regards, namely, the efficiency and the randomness. First, the efficiency of the algorithm is demonstrated from the perspective of computational complexity. Secondly, the sample that is generated by the algorithm is used to verify the randomness of the random algorithm in particle size and location by the statistical hypothesis testing method.

4.1. Computational Efficiency

Taking the uniform distribution as an example, it is assumed that n   particles obey the random distribution, and the efficiency analysis of the take-and-place algorithm and the algorithm proposed in this paper is as follows:
(1) The major steps of the take-and-place algorithm are as follows: randomly generate the position of the new particle and determine whether there is any interference between the current particle and the existing particles. If there is no interference, the relevant data of the current particle are saved, and the particle is accepted; otherwise, the algorithm needs to try new locations again and again. Therefore, the take-and-place algorithm not only reduces the randomness of sample generation, but also comes with the computational complexity of   c 1   ( m 1 + 2 m 2 + , + ( k 1 ) m k + , + ( n 1 ) m n ) that is far greater than O   ( n 2 ) . To clarify, m k is the number of repeated generation positions when filling the k th particle and a positive constant c 1 is independent of n ;
(2) As mentioned in Section 2, the algorithm that is proposed in this paper consists of two parts. In the first part of the algorithm, the position of each particle is only related to a small number of particles in its neighborhood when the initial position is generated or the tightness of the particles is adjusted. Therefore, the computational complexity is c 1   n 1 , where   n 1 is the number of particles and a positive constant c 1 is independent of n 1 . In the second part of the algorithm, each particle is independently processed, and the complexity is   c 2   n 1 , where   c 2 is a constant independent of n 1 . Accordingly, the overall computational complexity of our algorithm is O ( n ) .
The computational complexity of the particle random generation algorithm that is proposed in this paper is significantly lower than that of the traditional algorithm, and it is highly efficient.

4.2. The Randomness of Generated Samples

4.2.1. Random Test of Particles Size

At present, most of the simulations focus on the case of uniform and random distribution, so the random characteristics of the samples with uniform and random distribution are studied, and an χ 2 test [16] is used to evaluate the uniformity of the samples. The statistic test formula is
χ 2 = m ϵ i = 1 m ( ϵ i ϵ m ) 2
where ϵ is the number of all random numbers, m   is the number of intervals, and ϵ i is the number of the i th interval. In the uniformity test of the algorithm, we generated 3485 samples whose radii are on the interval (0.05, 0.15). The percentages of each radius are shown in Figure 5.
The particles sizes are evenly distributed in each interval. Assuming that the degree of freedom is 9, the calculated χ 2 value of radius R is 11.72, and the asymptotic significance is 0.230, greater than 0.05. According to statistic theory, the χ 2 distribution table shows that the distribution of particle radius R that is obtained by the algorithm can be considered as an ideal uniform random distribution. The test statistical scale is shown in Table 3.

4.2.2. Random Test of Particle Position

Taking normal distribution as an example, the Shapiro–Wilk test method [17] and the Kolmogorov–Smirnov test method [18,19] are used to test the normality of the samples that are generated by the random algorithms, using the compactness algorithm to make the spherical center coordinates of 301 particles normally distributed along the x-axis. Table 4 is the analysis result of the normality test.
As can be seen from Table 4, the sample data did not show statistical significance (p > 0.05), which means that the hypothesis (hypothesis: the data are normally distributed) is accepted and the sample has the characteristic of normality. The histogram of the sample distribution along the x-axis is shown in Figure 6.

5. Permeability of Porous Media

Titanium porous filter cartridge is made of porous titanium metal filter material by the powder metallurgy method. Its internal pores are curved and crisscross, and the filtering mechanism is typical deep filtration. Permeability is the ability of a porous medium to allow fluid matter to pass through. In this section, the lattice Boltzmann method is used to study the flow of fluid under the condition of mass particle filling and the influence of initial particle distribution on fluid permeability by using the meso-model of porous media that is constructed by the above algorithm.

5.1. Lattice Boltzmann Method

The motion of a fluid can be described by a set of partial differential equations, such as the Navier–Stokes equations, which are highly nonlinear in most cases and find it very difficult to obtain analytical solutions. The Lattice Boltzmann method is used to solve the numerical solution of the fluid motion equation by means of the discrete method. The lattice Boltzmann method can be regarded as a special discrete scheme of a continuous Boltzmann equation, as shown in the following formula.
g i ( x + e i δ t , t + δ t ) g i ( x , t ) = Ω i ( x , t )
where g   is the discrete distribution function, e   is the velocity space, i is the type of velocity, δ t is the discrete time step,   t is the current time step, x is a point on the grid, and Ω is the change caused by collision.
According to the operator that is proposed by Bhatnagar, Gross, and Krook [20], Ω i ( x , t )   can be replaced, as shown in the following formula.
g i ( x + e i δ t , t + δ t ) g i ( x , t ) = 1 τ ( g i e q g i )
where g i e q is an equilibrium distribution function to be determined; τ is the relaxation time.
The D d Q m model that is proposed by Qian et al. [21] is the basic model of the lattice Boltzmann model, where d represents d-dimensional space and m represents the number of discrete velocities of the lattice. In this paper, the velocity is discretized into 19 directions in three-dimensional space, as shown in Figure 7.
The discrete component of the velocity is
E = [ e 0 , e 1 , e 2 , e 3 , e 4 , e 5 , e 6 , e 7 , , e 18 ]
The equilibrium equation [21] of this model is
g i e q = ρ m ω i [ 1 + e i · μ c s 2 + ( e i · σ ^ ) 2 2 c s 4 σ ^ 2 2 c s 2 ] i = 0 , 1 , 2 , 18
where C s is the lattice speed, ω is the weight coefficient, σ ^ is the fluid velocity, and ρ m is the fluid density. The lattice velocity C s is as follows.
C s = c 3
where c is the ratio of grid step to time step, and its value is 1. The calculation of ω i is shown in the following formula.
ω i = { 1 / 3 1 / 18 1 / 36 e i 2 = 0 e i 2 = c 2 e i 2 = 2 c 2
Let g i ( x t , t )   be the distribution function of time t , at lattice point   x , velocity e , then the evolution equation of the distribution function is
g i ( x + e i δ t , t + δ t ) g i ( x , t ) = 1 τ ( g i e q ( x t , t ) g i ( x t , t ) )

5.2. Darcy’s Law

When single-phase flow flows in porous media at a low Reynolds number, it follows Darcy’s Law, also known as the Darcy model, which is one of the most basic and commonly used mathematical models for macroscopic seepage, as shown in the following equation
σ = h ν p l
where σ is the Darcy velocity, h   is the permeability of porous media, ν is the universality coefficient of the fluid, p l is the fluid pressure, and   is the Hamiltonian operator.
According to the distribution function, the macroscopic parameters of the fluid density ρ l and fluid velocity σ can be obtained from formula
ρ l = i g i σ = 1 ρ l i g i · e i
The permeability of porous media can be calculated according to Darcy’s formula [22].
h = ν σ ¯ p l
where h is the permeability of the porous media, ν is the coefficient of motion universality, p l is the pressure difference, and σ ¯ is the average speed.

5.3. Lattice Boltzmann Method Simulation Procedures

The lattice Boltzmann method simulation program structure is collision-migration. The specific process is as follows:
(1)
Setting of initial conditions;
(2)
Execute collision at time t ;
(3)
Boundary processing;
(4)
Calculate macroscopic quantities;
(5)
Check whether convergence exists. If not, return to Step 2. Otherwise, go to the next step;
(6)
Output the result.
Suppose the fluid flows in a 4.5 × 4.5 × 4.5 mm cube tube, as shown in Figure 8. The fluid is the single-phase flow of a low Reynolds number and is incompressible. The same number of lattice points are used in the x, y and z directions in the grid division of the stacking space. The flow direction of the fluid is in the z direction, the top and bottom surfaces are the fluid outlet and inlet surfaces, respectively, and the other surfaces are boundary surfaces. When the relaxation coefficient is fixed at 1, the universality coefficient is 10 6 . At the beginning, the velocity in the whole flow field is set at 0 and the density is 1. The method to judge the convergence is that in 50 time steps; if the change of the sum of the absolute values of the velocity along the direction of fluid flow is less than 0.0001, convergence is realized.

5.4. Simulation Results

5.4.1. The Influence of the Distribution Law

The model is divided into 1283 grid points, and each grid point has 19 velocity directions. The algorithm that is proposed in this paper is used to generate uniform distribution, normal distribution (mean is 3, and the standard deviation is 0.5) and exponential distribution, respectively (rate parameter λ is 2). The simulation results are shown in Table 5.
The above table shows that the average velocity and permeability of exponential distribution are the largest. There is little difference between a uniform distribution and a normal distribution. It should be noted that this example only calculates permeability under different probability density distributions to verify the feasibility of the modeling algorithm. The parameters that are related to the probability density function, such as the direction of the distribution law and the parameters of the distribution function, have significant effects on permeability, and their effects are often coupled with each other. Therefore, further research is needed to establish a general strong law before these questions can be resolved. This modeling algorithm provides a model foundation and technical support for further study of these laws.

5.4.2. The Influence of the Numbers of Particles and Grids

In the case of the lattice Boltzmann method, both the number of lattice points and the volume fraction of the particles become significant influencing factors. On the premise that the particle distribution law is determined as uniform distribution, the permeability of the model is simulated for the following two cases: (1) the volume fraction of the particles is 0.55 0.57, 0.61,0.63, respectively, and the number of grid points stays the same as 1283; (2) the volume fraction is 0.60 and the number of grid points is 643, 1283, 2563, 5123, respectively. The simulation results are shown in Figure 9 and Figure 10.
Through analytical methods [23], the relationship between the permeability and volume fraction of porous materials can be predicted by specific formulas. The permeability of porous materials is related to the minimum and maximum radius of the cross-section through the transport phase, the shortest path of transport and the volume fraction of the pores. The specific formula is as follows.
κ ^ = ( 0.94   r m i n + 0.06   r m a x ) 2   ε 2.14 τ 2.44   / 8  
where κ ^ is the predicted value of permeability, r m i n and   r m a x are the minimum and maximum radius of the cross-section through the transport phase, ε is the volume fraction of the pores, and τ is the shortest path of transportation.
It can be seen from the formula that the permeability of porous materials increases with the increase in the volume fraction of the pores. It is worth noting that the volume fraction of pore ε and the volume fraction of particle α are complementary, and the sum is 1. This is consistent with the change trend that is shown by the simulation.
As shown in Figure 9, the simulation results show that with the same distribution law, the permeability of the fluid through this distribution decreases with the increase in the volume fraction. Figure 10 shows that under the same distribution law, the permeability of the fluid through the same amount of particle decreases with the increase in the number of lattice points. When the lattice points increase to 1283 or higher, the permeability is almost unchanged.

6. Conclusions

In this paper, an efficient algorithm is proposed to generate the particle meso-model by specifying the probability density distribution. The following conclusions and innovations can be obtained:
  • The filling process avoids the huge calculation burden that is caused by the continuous iterative interference detection;
  • In the process of generating the initial position of the particles or realizing the compact stacking, the changes of each particle are only related to its small neighborhood, so it has a high compactness efficiency;
  • The computational complexity of the algorithm is first order, while that of the traditional algorithm is much larger than second order, which illustrates the computational efficiency of the proposed algorithm;
  • With this algorithm, the size and position of the particles can be distributed according to the arbitrary probability density based on the requirements, and the specified volume fraction can also be realized according to the requirements.
The results provide a new method for the space filling and modeling of porous materials and provide technical support for 3D printing when optimizing porous material structures. This method provides a modeling basis for further exploring the influence of porous materials with a specified probability density on the mesoscopic and multi-field simulation research of sound absorption, heat absorption, radiation resistance, seepage and extrusion resistance, which are difficult to observe by experimental methods.

Author Contributions

Methodology, W.Z. and F.W.; Software, W.Z. and L.H.; Data Analysis, L.H. and G.Z.; Supervision, L.H. and G.Z.; Writing—original draft, W.Z.; Writing—review & editing, F.W. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by Key Research and Development Projects [grant numbers 2022GY-399] of Shaanxi Province, China.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

Symbols and Abbreviations

f probability density
M 3-dimensional array
ρ   a small positive value
n n e i g h b o r a small positive integer
v volume of the cube
u probability density
v g total volume of the particles
α   volume fraction
n particle number
M k number of repeated generation position
c 1 a positive constant
ϵ number of all random numbers
g   discrete distribution function
e velocity space
δ t discrete time step
t   current time step
Ω change caused by collision
C s lattice velocity
ω weight coefficient
σ ^ fluid velocity
σ Darcy velocity
h permeability of porous media
p l fluid pressure
Hamiltonian operator
κ ^   predicted value of permeability
λ rate parameter of the exponential distribution

References

  1. Noubactep, C.; Caré, S.; Togue-Kamga, F.; Schöner, A.; Woafo, P. Extending Service Life of Household Water Filters by Mixing Metallic Iron with Sand. Clean Soil Air Water 2010, 38, 951–959. [Google Scholar] [CrossRef] [Green Version]
  2. Wang, W.; Parker, K.H. Effect of deformable porous surface layers on the motion of a sphere in a narrow cylindrical tube. J. Fluid Mech. 1995, 283, 287–305. [Google Scholar] [CrossRef] [Green Version]
  3. Sobieski, W. Waterfall Algorithm as a tool of investigation the geometrical features of granular porous media. Comput. Part. Mech. 2021, 9, 551–567. [Google Scholar] [CrossRef]
  4. Primera, J.; Hasmy, A.; Woignier, T. Numerical study of pore sizes distribution in gels. J. Sol-Gel Sci. Technol. 2003, 26, 671–675. [Google Scholar] [CrossRef]
  5. Tory, E.M.; Cochrane, N.A.; Waddell, S.R. Anisotropy in Simulated Random Packing of Equal Spheres. Nature 1986, 220, 1023–1024. [Google Scholar] [CrossRef]
  6. Jodrey, W.S.; Tory, E.M. Computer Simulation of Isotropic, Homogeneous, Dense Random Packing of Equal Spheres. Powder Technol. 1981, 30, 111–118. [Google Scholar] [CrossRef]
  7. Cundall, P.; Strack, O. Discussion: A discrete numerical model for granular assemblies. Geotechnique 1980, 30, 331–336. [Google Scholar] [CrossRef] [Green Version]
  8. Bezrukov, A.; Bargiel, M.; Stoyan, D. Statistical analysis of simulated random packings of spheres. Part. Part. Syst. Charact. Meas. Descr. Part. Prop. Behav. Powders Other Disperse Syst. 2002, 19, 111–118. [Google Scholar] [CrossRef]
  9. Bailakanavar, M.; Liu, Y.; Fish, J.; Zheng, Y. Automated modeling of random inclusion composites. Eng. Comput. 2012, 30, 609–625. [Google Scholar] [CrossRef]
  10. Xin, Z.; Miao, W.; Wang, Y.; Chen, H. Simulation of Tensile Fracture of ZrO2 Toughened Al2O3 Particles Reinforced Fe45 Composites by Finite-Discrete Element Coupling Method. Acta Mater. Sin. 2019, 36, 1471–1479. (In Chinese) [Google Scholar]
  11. Stevenl, B.; Peter, R.K.; David, M. Network Model Evaluation of Permeability and Spatial Correlation in a Real Random Sphere Packing. Transp. Porous Media 1993, 11, 53–70. [Google Scholar]
  12. Heijs, A.W.J. Numerical evaluation of the permeability and the Kozeny constant for two types of porous media. Phys. Rev. E 1995, 51, 4346–4353. [Google Scholar] [CrossRef] [PubMed]
  13. Bryant, S.; Blunt, M. Prediction of relative permeability in simple porous media. Phys. Rev. A 1992, 46, 2004–2011. [Google Scholar] [CrossRef] [PubMed]
  14. Raats, P.A.C. Dynamics of Fluids in Porous Media. Soil Sci. Soc. Am. J. 1973, 37, vi. [Google Scholar] [CrossRef]
  15. Zhu, J.; Xi, Z.; Tang, H.; Tan, P. Characterization of porous structures and fractal theory. Rare Met. Mater. Eng. 2006, 35, 452–456. [Google Scholar]
  16. Hosmane, B.S. Improved likelihood ratio tests and Pearson chi-squared tests for independence in two dimensional contingency tables. Commun. Stat. 1986, 16, 1875–1888. [Google Scholar] [CrossRef]
  17. Kline, R.; Kline, R.B.; Kline, R. Principles and Practice of Structural Equation Modelling. J. Am. Stat. Assoc. 2005, 4, 1941–1947. [Google Scholar]
  18. Drezner, Z.; Zerom, O.T.D. A Modified Kolmogorov-Smirnov Test for Normality. Commun. Stat.–Simul. Comput. 2010, 39, 693–704. [Google Scholar] [CrossRef] [Green Version]
  19. Lilliefors, H.W. On The Kolmogorov-Smirnov Test for Normality with Mean and Variance Unknown. J. Am. Stat. Assoc. 1967, 62, 399–402. [Google Scholar] [CrossRef]
  20. Bhatnagar, P.L.; Gross, E.P.; Krook, M. A model for collision processes in gases: I. Small amplitude processes in charged and neutral one-component systems. Phys. Rev. 1953, 94, 511–525. [Google Scholar] [CrossRef]
  21. Qian, Y.H.; D‘Humières, D.; Lallemand, P. Lattice BGK Models for Navier-Stokes Equation. EPL (Europhys. Lett.) 1991, 17, 479–484. [Google Scholar] [CrossRef]
  22. Maier, R.S.; Kroll, D.M.; Kutsovsky, Y.E.; Davis, H.T.; Bernard, R.S. Simulation of flow through bead packs using the lattice Boltzmann method. Phys. Fluids 1998, 10, 60–74. [Google Scholar] [CrossRef]
  23. Neumann, M.; Stenzel, O.; Willot, F.; Holzer, L.; Schmidt, V. Quantifying the influence of microstructure on effective conductivity and permeability: Virtual materials testing. Int. J. Solid Struct. 2020, 184, 211–220. [Google Scholar] [CrossRef]
Figure 1. Principle of initial position determination. (a) The position of new particle generated (i = 1). (b) The position of new particle generated (i ≠ 1).
Figure 1. Principle of initial position determination. (a) The position of new particle generated (i = 1). (b) The position of new particle generated (i ≠ 1).
Materials 15 04733 g001
Figure 2. White noise waveform.
Figure 2. White noise waveform.
Materials 15 04733 g002
Figure 3. Results of compactness algorithm execution.
Figure 3. Results of compactness algorithm execution.
Materials 15 04733 g003
Figure 4. The result of the filtering algorithm execution. (a) Uniform distribution along x, y, and z, (b) Normal distribution along the x direction, (c) Exponential distribution along the z direction.
Figure 4. The result of the filtering algorithm execution. (a) Uniform distribution along x, y, and z, (b) Normal distribution along the x direction, (c) Exponential distribution along the z direction.
Materials 15 04733 g004
Figure 5. Distribution of particle radius in each interval. (a) Percentage of particles in each radius, (b) The number and standard deviation of particles within each radius.
Figure 5. Distribution of particle radius in each interval. (a) Percentage of particles in each radius, (b) The number and standard deviation of particles within each radius.
Materials 15 04733 g005aMaterials 15 04733 g005b
Figure 6. Particle coordinate distribution along the x direction.
Figure 6. Particle coordinate distribution along the x direction.
Materials 15 04733 g006
Figure 7. Discrete velocity diagram of lattice points.
Figure 7. Discrete velocity diagram of lattice points.
Materials 15 04733 g007
Figure 8. Diagram of flow field.
Figure 8. Diagram of flow field.
Materials 15 04733 g008
Figure 9. Effect of volume fraction on fluid permeability.
Figure 9. Effect of volume fraction on fluid permeability.
Materials 15 04733 g009
Figure 10. Effect of grid number on fluid permeability.
Figure 10. Effect of grid number on fluid permeability.
Materials 15 04733 g010
Table 1. Distribution effects of particles with different size ranges with the same probability density.
Table 1. Distribution effects of particles with different size ranges with the same probability density.
Range of Particle SizesThe Distribution of Particles
R = 0.05–0.15 mm Materials 15 04733 i001
R = 0.15–0.25 mm Materials 15 04733 i002
R = 0.25–0.35 mm Materials 15 04733 i003
Table 2. Distribution effects of particles with the same probability density and different volume fractions.
Table 2. Distribution effects of particles with the same probability density and different volume fractions.
Volume Fraction Total   Volume   ( m m 3 ) The Distribution of Particles
10%12.550 Materials 15 04733 i004
20%25.504 Materials 15 04733 i005
30%37.743 Materials 15 04733 i006
Table 3. χ 2 test statistical table.
Table 3. χ 2 test statistical table.
Test Statistics
StatisticValue
χ 2 11.720
Degrees of freedom9
Asymptotic significance0.230
Table 4. Normal test analysis results.
Table 4. Normal test analysis results.
Sample SizeAverageStandard DeviationPartial DegreesKurtosisKolmogorov–Smirnov TestShapiro–Wilk Test
The Statistic D ValuepThe Statistic W Valuep
3013.0140.511−0.028−0.2190.0290.7420.9960.610
Table 5. Permeability and average velocity of different distribution.
Table 5. Permeability and average velocity of different distribution.
DistributionAverage Velocity
( m / s )
Permeability   k   ( m 2 )
Uniform0.586 × 10−35.204 × 10−10
Normal0.533 × 10−35.335 × 10−10
Exponential0.681 × 10−36.757 × 10−10
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Zhang, W.; He, L.; Wang, F.; Zhang, G. A Stochastic Filling and Modeling Algorithm of Non-Equal Diameter Particles with Specified Probability Density for Porous Permeable Materials. Materials 2022, 15, 4733. https://doi.org/10.3390/ma15144733

AMA Style

Zhang W, He L, Wang F, Zhang G. A Stochastic Filling and Modeling Algorithm of Non-Equal Diameter Particles with Specified Probability Density for Porous Permeable Materials. Materials. 2022; 15(14):4733. https://doi.org/10.3390/ma15144733

Chicago/Turabian Style

Zhang, Wei, Lile He, Fazhan Wang, and Guangyong Zhang. 2022. "A Stochastic Filling and Modeling Algorithm of Non-Equal Diameter Particles with Specified Probability Density for Porous Permeable Materials" Materials 15, no. 14: 4733. https://doi.org/10.3390/ma15144733

APA Style

Zhang, W., He, L., Wang, F., & Zhang, G. (2022). A Stochastic Filling and Modeling Algorithm of Non-Equal Diameter Particles with Specified Probability Density for Porous Permeable Materials. Materials, 15(14), 4733. https://doi.org/10.3390/ma15144733

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop