Simulation of Particle-Fluid Interaction in Fractal Fractures Based on the Immersed Boundary-Lattice Boltzmann Method

In order to reveal the influence of particles on fluid flow characteristics in rough fractures under fluid-solid coupling, a range of fracture systems with varying roughness were generated using the Weierstrass-Mandelbrot function. Fluid-particle interactions in rough fractal fractures were simulated using the immersed boundary-lattice Boltzmann method. In this paper, the effects of fluid viscosity, particle size, particle quantity, fracture fractal dimension, and particle grading composition are studied. Results illustrate that increasing fluid viscosity hinders the movement of particles, resulting in the decreasing of particle velocity. As particle size and particle quantity increase, the particle occupation of the channel area grows larger, which lead to lower permeability of the channel. Increasing fracture fractal dimension surges the curvature of the fluid channel, but permeability has a negative exponential correlation to fractal dimension. With increasing particle grading composition, the blocking effect of particles on fracture flow also increases with increasing particle proportion.


Introduction
Fractures are widespread in rock masses, and fluid flow through fractures affects their mechanical properties by inducing deformation, fracturing, softening, argillization, and corrosion [1]. For oil shale and coal reservoirs, grand and undetected fractures may lead to significant engineering problems because they prominently alter groundwater movement [2,3]. As fractures increase the transport distance of contaminants, a better understanding of fracture geometry and fracture occurrence will contribute to a more accurate determination of water protection areas [4,5]. Particle migration and the consequent permeability reduction can cause reservoir damage in numerous petroleum, environmental, and water resource technologies [6,7]. It is of great scientific significance and engineering application value to accurately and quantitatively describe the structure characteristics of fractured rock and evaluate the influence of internal particle migration.
The permeability performance and seepage direction of rock masses are not only related to the expansion and cutting characteristics of the fracture network [8] but are also closely related to the geometric characteristics, aperture, roughness, connectivity, and contact degree of a single fracture [9][10][11]. Fracture roughness is the essential element that causes deviation between the actual fracture permeability and theoretical permeability.
Sedimentary transport and adhesion plugging in the fracture space considerably change the permeability characteristics of fractured rock masses, resulting in a change in the groundwater level and the redistribution of pore water pressure within the formation, eventually inducing geological disasters [12][13][14]. Presently, the porous media capillary bundle model has been widely studied to explore the sediment transport mechanism [15,16]. The model simplifies the coal medium into a uniform combination of capillary units, and each unit contains certain amount of capillaries. Happel [17] proposed an ideal model that considers hydrodynamic and gravity effects. The method is based on a granular porous medium model, which assumes that the porous medium is composed of numerous spherical particle elements. Yao [18] proposed an analysis method for determining particle motion paths and determined that the control factors of suspended particle motion are diffusion, sedimentation, and interception. For the sake of simplicity, the coarse-grained model ignores the dynamics of fluids and the Van Der Waals forces between molecules. Tufenkji and Elimelech [19] revised the governing equations of diffusion, sedimentation, and interception in the particle motion path analysis method.
James and Chrysiko [20] believe that suspended particles migrate at high speed with central fluid in the pores. Raizah and Ahmed [21] examined the influence of two crosssection heated ellipses on the fluid flow and heat transfer inside porous enclosures. Tien, Sharma, and Mattson et al. [22,23] established a 3D fracture simulation system. In this method, the Lagrangian method was used to track solid particles, and the influence of the density, size, and flow rate of solid particles on fluid migration velocity was considered. Oseen [24] modified the resistance expression of Stokes formula by comprehensively examining the influence of the inertia term. In porous media, Civan and Shamar [25] divided particle sedimentation into surface precipitation, pore throat blockage, pore filling, and filter cake formation.
In the past decades, there have been numerous direct numerical simulation methods for particle-fluid coupling, such as the finite element-arbitrary Lagrangian Euler method, fictitious domain method, lattice Boltzmann method, immersed boundary method, and smoothed particle hydrodynamics method [26][27][28][29][30][31]. Feng et al. [29] proposed the immersed boundary-lattice Boltzmann method (IB-LBM) that couples the computational fluid dynamics and immersion boundary method. In this method, fluid flow is simulated using the Navier-Stokes equation, and the hydrodynamic force applied to the solid particles can be calculated using the structure and velocity of the solid particles. Unlike the lattice Boltzmann method, the IB-LBM method uses a Cartesian grid and a set of Lagrangian points to discretize the boundaries of solid obstacles. Discrete Lagrangian lattice points are closer to the curve boundary than those in a staircase fashion. Simultaneously, the IB-LBM solves the Navier-Stokes equations in a rectangular domain, where flexible boundaries can move or change shapes in a complicated fashion. The boundary force is applied to the fluid using a chosen form of the Dirac function. The IB-LBM method has a higher efficiency than bodyfitted methods, e.g., finite-element method [32].
In order to understand the flow mechanism through microchannels with varying morphological structures and the influence of sedimentary particles during drainage, the IB-LBM is utilized to analyze and calculate the fluid seepage behavior in a rock mass with a two-dimensional fracture structure in this contribution. Using the Weierstrass-Mandelbrot function, a series of fractal curves with varying roughness and microscale rock fracture structure model are established.
Through the numerical simulations of the steady-state seepage within the fractured structure, the velocity field distribution and the equivalent permeability coefficient of the fluid are calculated. Through the analysis of simulation results, the relationship between the fracture fractal dimension and equivalent permeability coefficient is established. The seepage behavior of fluid in complex fractures with par-ticles is explored, and the influence of particle size, particles quantities, fracture fractal dimension, and particle grading composition on fluid passage is considered.

Numerical Methods
2.1. Lattice Boltzmann Method. The Navier-Stokes equations are the fundamental partial differential equations that describe incompressible fluid flow. Using the rate of stress and rate of strain tensors, the velocity-pressure formulas can be written as where u is the fluid velocity, t is the time, ν is the kinematic viscosity, p is the pressure, ρ is the fluid density, and F is the forcing term. In the IB-LB method, the effect of a solid boundary pushing or blocking a fluid is convert into a force term, which is the most essential component of F. The other part of F is the conservative force, such as gravity or a magnetic force.
The fluid velocity must satisfy the Mach number Ma = μ f /c s < <0:1 in order to meet the condition of an incompressible fluid, where μ f is the fluid velocity and c s is the lattice speed of sound.
In the LB method, the D2Q9 model is the commonly used discrete format for solving Eqs. (1) and (2). The fluid flow can be determined by resolving the particle collision and streaming processes, and the lattice Boltzmann equation is used to determine the streaming and collision processes of fluid particles. In this model, velocity is microscopically discretized into 9 directions and the discrete velocity c i = ½0, 0 ; 1, 0 ; 0, 1;−1, 0 ; 0,−1 ; 1, 1;−1, 1;−1,−1 ; 1,−1 ( Figure 1). The density distribution function is the only quantity of concern, and the discrete lattice Boltzmann equation can be written as where f i is the density distribution function of moving populations. Here, τ is the dimensionless mean relaxation coefficient, and δt is the time step. The equilibrium distribution function f eq i ðx, tÞ is obtained using a Taylor series expansion using the Maxwell-Boltzmann distribution function with velocity u up to the second order, which can be written as where the ω i is the weight coefficient associated with velocity c i , and the nondimensional lattice speed of sound is modeldependent. In the D2Q9 model, c s commonly equals 1/ ffiffi ffi 3 p and ω i = ð4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36Þ.

Geofluids
Applying the Chapman-Enskog multiscale expansion method, the Navier-Stokes equations can be recovered using the LB equation with first-order time precision and secondorder space precision. Density and local momentum are determined using the summation and the first order moment of the density distribution function: Performing the Chapman-Enskog expansion on the LB dynamics with a low Mach number limit gives the following expression for kinematic viscosity: where δr is the lattice spacing.
In order to implement IB in the LB method, Eq. (3) can be separated into collision and streaming components. Thus, in the IB-LB method, the collision and streaming components are Streaming : As mentioned above, F i ðx, tÞ is the discrete force term. In this paper, the discrete method for the force term proposed by Guo et al. [33] is used, and the discrete force term is However, the existence of the external force term changes the density and momentum; hence, a velocity correction is needed: There are numerous methods for the dispersion of force terms, but the accuracy of this method is relatively high, and a second-order precision in space can be achieved. In addition, the implementation of the method is relatively simple.

Immersed Boundary Method.
Coupling LB and IB requires a procedure to implement the effect of the body force between the Eulerian grids and Lagrangian grids. Currently, there are two methods to couple force terms, the direct force method (DF) [34] and interpolation force method (IF) [35]. The DF method aims at directly obtaining the force term on the Euler grid using the Dirac δ function under the condition that the solid boundary stress is already known. However, the IF method aims to determine the velocity and force terms on the Euler grid through high-precision interpolation. Dupuis et al. compare the differences between the two methods in detail and found that the difference between the two methods was not notable [36]. The DF is adopted in this paper to exert the force term. The computational scheme for the implementation of the IB-LBM in the particle-fluid interaction problem can be found in the paper by Liu et al. [37].
In order to model the physical surface of a particle, we consider a system with N physical particles with M discrete points. We assume that the reference point of the particle is stationary, and all boundary points are on a circle. At time t, the center of the particle is at x i ðtÞ, and the instantaneous particle rotational matrix is R i ðtÞ. The position of a reference point j in a Lagrangian coordinate system with respect to the particle i may be determined using the following formula: where x r ij ðtÞ is the boundary reference point of the particles at timestep t. The linear restoring force will generate when the reference point and the boundary point are not in the same position. We suppose that the particle deviation is ξ ij ðtÞ = x ij ðtÞ − x r ij , where the subscripts indicate the ith particle and the jth discrete point. The restoration force can be calculated as [37] where k is the spring constant. The internal link force in the simulation of the rigid particle motion is simply the balance force that was mentioned above. Hence, the total force F i and the total torque T i exerted on the ith particle are where m i is the mass of the ith particle.

Geofluids
In order to satisfy the no-slip boundary condition on the particle-fluid interface, the velocity on the particle point must strictly be the same as on its neighboring fluid. We can obtain the velocity of the Lagrangian points from the velocity of the adjacent fluid.
where Xðs, tÞ is the displacement of solid particles. We can also use the Dirac δ function to composite the fluid and particle: When the IB-LBM is used to simulate particle motions in fluid, the include external force (such as buoyancy), internal force (such as the particle collision force), and restoration force can be calculated. Hence, the total force exert on the particle surface node can be written as [37] A collision technique is applied if the gap between the particles is less than a given thresholdh. The functional form of the repulsive force is given by Glowinski et al. [38]: The force density fðx, tÞ is obtained from the immersed The Dirac δ function is a continuous interpolation function. We cannot interpolate all the points one by one in the interpolation process; therefore, the Dirac δ function must be approximated at the grid points. Numerous methods exist for approximating the Dirac δ function [39]. We use a more continuous method: 2.3. Characterization of Fracture Roughness. The rock damage process is actually the appearance, growth, bifurcation, and intersection of microcracks in the rock to form a macroscopic fracture. The seepage structure of fractures commonly displays fractal characteristics [40][41][42]. The three-dimensional fracture structure network of rock masses has complex structural characteristics composed of multiple intersecting and combined two-dimensional fracture surfaces; however, the two-dimensional fracture surface is composed of multiple groups of single fractures. The structure curve of a single fracture obtained by cutting the crosssection along a specific section presents an irregular fractal structure.
A quantitative description method of single fracture can determine the influence of the rough structure fracture shape on fluid seepage and clarify the fluid-particle coupling mechanism in the rough structure. Firstly, we use the Weierstrass-

Geofluids
Mandelbrot function to obtain fractal curves for different fractal dimensions D, which is where γ is a parameter greater than 1, and phase ϕ n is any angle that makes W exhibit deterministic or stochastic behavior. D ∈ ð1, 2Þ is the Hausdorff fractal dimension of the graph of WðtÞ. The fractal governing function CðtÞ is the real part of WðtÞ.
The fractal dimension D is equal to 1.0, 1.1, 1.2, 1.3, 1.4, and 1.5, respectively (Figure 2). During the fracture curve generation process, γ is fixed at 1.4, and the time control parameter is t ∈ ð0 : 0001 : 1Þ and n ∈ ð0 : 100Þ. The numerical simulation model is 10 mm long by 4 mm wide, and the average gap width of the fracture is 1 mm.
Judging from the shape of the fractal structure curves, with increasing fractal dimension, the crack surface gradually becomes steeper, and wrinkle degree surges. Conversely, the smaller the fractal roughness index, the smoother the fractal surface.
The influence of roughness on fluid flow with particlefluid interaction is investigated using a two-dimensional incompressible model. The 2D profile of the fracture surface is imported into the IB-LBM in order to generate a fractured model (Figure 3). The length of the fracture model is set to 10 mm to be consistent with a realistic profile. The aperture between the two parallel rough profiles is 1 mm. The pressure is applied at the inlet and outlet of the fracture, and the bounce-back scheme is used to simulate a nonslip boundary at the fluid-solid interface.

Benchmark-Single Particle Sedimentation in Fluid.
In order to validate the accuracy of the immersed boundarylattice Boltzmann method, we adopt a model that considers the sedimentation of a circular particle under gravity and buoyancy in a channel filled with a viscous fluid. In 2006, Wan and Turek [43] simulated this single-channel model using FEM, which considered the initial state of water flow in the pipeline as static. In order to verify the correctness of the model in this paper, we compared the results of the immersion boundary-lattice Boltzmann simulation with the FEM results of Wan and Turek. The dimension of the box is 2 cm (in the X-direction) by 6 cm (in the Y-direction). The fluid domain is divided into 200 × 600 square lattices with spacing of h = 100 μm. The viscosity and density of the fluid are 0.01 Pa·s and 1 g/cm 3 , respectively. The radius of the particle is 0.125 cm, the particle density is 1250 kg/m 3 , and the relaxation time in the simulation is 0.53. The four boundaries of the model are stationary walls, which are noslip boundaries for the fluid. Initially, one stationary solid particle is generated at the position (1 cm, 4 cm). Owing to gravity, a solid particle will move downward gradually.
During the simulation, the particle position coordinates and flow field velocity distribution were recorded. Figure 4 shows the curve of particle velocity in three cases with varying viscosities as time increases. Particle velocity increases continuously in the initial state due to gravity, and at about 0.3 s, the speed reaches a steady state (Figure 4). The particles produce a vertical upward water resistance during acceleration, and the resistance is positively correlated to the velocity. When the velocity reaches a certain value, the gravitational force and the resistance force on the particle remain the D A H L P out P in Figure 6: Schematic of flow around a circular particle close to an irregular wall.  : Velocity contours for a circular particle in a corrugated wall channel. 6 Geofluids same, and the particle reaches the ultimate velocity. The ultimate velocity of the particle will gradually decrease as viscosity increases because the viscosity of the fluid increases the viscous resistance of the particles (Figure 4). Our simulation results are consistent with the FEM results Wan and Turek in 2006, which indicate the correctness of the IB-LBM. Particle movement and fluid velocity contours at times of 0, 0.2, 0.4, 0.6, 0.8, 0.9, and 1.0 seconds are given in Figure 5. Particle sedimentation impacts the four channel walls. When the particles pass through a certain area, a large impact force is generated on the vertical pipe wall of the area. The particles squeeze the fluid, causing it to spread out on both sides, and the fluid moves in the opposite direction relative to the particle movement. A smaller impact area is formed in front of the particle movement, and a larger vortex area is formed behind the particle. The difference between these two areas will gradually increase as particle speed increases.

Particle Disturbance in Rough Fractures
. This section will focus on the effects of the particulate matter on fluid flow in irregular fractures, including the influence of the particle size on channel permeability, the effect of particle number on fluid flow, the filling behavior of particles in rough fractures with varying fractal dimensions, and the influence of particle gradation on seepage effects. Figure 6 shows a schematic diagram of the flow configuration with our coordinate system and the geometrical parameters adopted in this study. In order to examine the effect of the particle size on the seepage in the fractured channel, we set up a sinusoidal uplift wall with a height of A, and three circular particles with varying diameters were set up in the trough of the wall.

Effect of the Particle Size on Fluid Flow.
Refer to the parameters of the benchmark program above, the density and viscosity of the fluid are set to 1000 kg/m 3 and 0.01 Pa·s, respectively. The particle density (1250 kg/m 3 ) is slightly heavier than that of the fluid. The inlet and outlet are pressure boundary conditions, the pressure difference is set to 0.01 kPa, and the upper and lower walls are set to complete bound back boundary conditions to ensure that the speed at the side wall is zero.
The spacing of each grid point is 10 μm, the length of the entire field is 5 mm × 1 mm, the number of grids is 500 × 100, and the relaxation time in the simulation is set to 0.53. In order to examine the effect of the particle diameter on fluid passage, we set the height A of the wall protrusion to 0.2 mm. The velocity distribution in the fluid channel under four different cases is analyzed. In the first case, fluid flow through a sinusoidal channel containing no particles was analyzed. The other three cases all examined flow through a sinusoidal channel containing a circular particle with diameters of 0.12, 0.24, and 0.48 mm.
The IB-LBM is used to simulate the fluid flow field and particle motion in this model to obtain the velocity contours of the flow field with particles of different sizes. Figure 6 shows the fluid flow path of different particles in the same channel. When the diameter of the particles is smaller than the height A of the corrugated curved channel, there is almost  7 Geofluids no effect on the seepage velocity in the channel, but as the diameter of the particles increases, the passage of the fluid in the channel is blocked by the particles (Figure 7). The fluid passage deteriorates, thereby increasing the velocity of this part of the fluid and simultaneously reducing the seepage area of the particles as well as the seepage flux. Extending this concept to a complex fracture network indicates that the particle size will hinder the channel only when the particle height reaches the height of the rough protrusions on the fractures, and this hindering effect will increase as particle size increases.
In order to quantitatively determine the influence of particle radius on the channel, the particle drag coefficient and the channel permeability are analyzed. The ratio of the drag force F D exerted by the fluid on the particle to the product of the projected area of the particle in the direction of movement and the fluid dynamic pressure is where F D represents the resistance of the particle surface along the incoming flow direction, U in represents the incident velocity of the fluid, and d represents the particle diameter. Figure 8 shows the drag coefficients of particles with varying diameters. From the stable value of the drag coefficient, we can find that particles with a large radius have greater resistance to the particles because the particles hinder the smooth flow of the fluid and increase the path of movement of the fluid, thereby increasing the velocity of the fluid and further increasing the drag coefficient. The fluid movement in the microscopic seepage channel can also be characterized by permeability, which represents the difficulty for the fluid to pass through complex fractures. Permeability at the microscopic scale is given by where Q is the flow quantity, which is estimated using a numerical integration algorithm, L is the partial path length, and A is the cross-sectional area of channel. Table 1 shows the flow rate and permeability change with particle radius. The existence of large particles changes the dominant path of the flow field, thereby changing the flux and permeability of the flow field (Figure 7). We can make the same observation using the change law of permeability. As particle size increases, the permeability of the flow field decreases. Owing to the obstruction effect of particles, the fluid cannot deliver effectively; so, the pressure drops, and the flux reduces, which presents the decline of permeability.

Effect of the Particle Number on the Flow Channel.
In order to study the influence of the number of particles on the fluid movement in complex fractal fractures, we set up a curved and complex fractal fracture keeping the structure of the fractal fracture unchanged. The simulation discussed the influence of the presence of particles on fluid flow when the number of particles is 0, 5, 10, and 15. The fractal dimension of the fracture structure is 1.1, the length of the channel is set to 10 mm × 4 mm, and the aperture of the fracture is set to 1 mm. The spacing of each grid point is 50 μm; so, the number of grids is 200 × 80. The density and viscosity of the fluid are 1000 kg/m 3 and 0.01 Pa·s, respectively. The density of the particles is slightly heavier than that of the fluid (1250 kg/m 3 ), and the radius of the particles is 0.1 mm. The inlet and outlet are pressure boundary conditions, the  Figure 9: Flow contours with particles in a fractal fracture with varying particle amounts. 8 Geofluids pressure difference is set to 0.01 kPa, and the upper and lower fracture surfaces are set to bounce back boundary conditions. Figure 9 shows the seepage process of fractures with a fractal dimension of 1.1 containing different numbers of particles. From numerical analysis of the contour map, it can be inferred that as the number of particles increases, the average seepage speed decreases (Figure 9). The presence of particles decreases the cross-sectional area of the fluid flow and hinders fluid passage. The fluid will form a velocity increase zone around the particles, but this has a lesser impact than the permeability reduction caused by the decrease in the fluid passage cross-sectional area. The more particles, the larger the fluid channel area occupied by the particles, and the stronger the hindrance. As the number of particles increases, the area of the fluid channel occupied by the particles becomes increases, and the blocking effect becomes stronger. Table 2 shows the seepage velocity, flux, permeability, and permeability reduction rate of complex fractures with varying particle numbers. The presence of more particles will significantly reduce the passage of fluid; that is, reduce the permeability of the fracture space (Table 2). Simultaneously, the permeability reduction coefficient positively correlates with the space occupied by the particles. The larger the space occupied by particles, the greater the permeability reduction  9 Geofluids coefficient. The fluid-particle interaction flow not only needs to overcome the hindrance of irregular boundaries but also needs to entrain the solid particles to move together, so it causes more energy expenditures. The greater the number of particles, the greater the energy expenditure. This is the main factor of permeability reduction caused by particles.

Particle Motion in Fractures with Varying Fractal
Dimensions. Considering the seepage process of complex fractures with varying fractal dimensions, this section constructs six fractal curves with varying fractal dimensions, which correspond to six fracture structures with varying roughness. The fractal dimensions of the six different fractal curves are 1.0, 1.1, 1.2, 1.3, 1.4, and 1.5, and the parameter settings of fluid and particles are exactly the same as those in Section 2.3. Figure 10 shows the contour map of fluidparticle coupled seepage velocity under different fractal dimensions. Fractal dimension has a significant influence on the flow. Increasing the fractal dimension means that the crack roughness increases, the tortuosity of the fluid passing through the crack increases, and the channel permeability decreases. Figure 11 shows the permeability curves under different fractal dimensions. The larger the fractal dimension of the fracture, the smaller the permeability of the channel. The increase of the fracture fractal dimension means that the structure is more complex, which indicates that the uplift and depression amplitudes are greater. When fluid flows through a fracture with a larger fractal dimension, the boundary impact that needs to be overcome is stronger, resulting in a decrease in permeability. From a numerical point of view, this relationship approximately satisfies the negative exponential correlation between permeability and fractal dimension, which has also been experimentally verified by Jafari and Babadagli [44]. The permeability of complex fractures has a positive correlation with the fractal dimension, and the negative exponential law is also approximately satisfied for fluids containing particles.

Blocking
Effect of the Particle Size Distribution. In order to simulate the effect of the particle size distribution on fluid seepage, the following four sets of seepage channel diagrams of size particles with varying proportions are set up ( Figure 12). The cavity size of the four models is 81 × 81 mm, the radius of the large particles is 9 mm, and the radius of the small particles is 4.5 mm. The large to small particle radius ratio is 2 : 1, and the ratio of the occupied area is 4 : 1. In order to eliminate the error caused by the position of the particles to the fluid channel, we attempt to make the position of the particles as symmetrical and regular as possible, and the area of the total cavity occupied by large and small particles is roughly equal. When the fluid flows through a complex seepage channel, it will autonomously flow through the superior channel; that is, flow in the place where the flow path is the shortest and the seepage width is wide ( Figure 12). Figure 12 shows the velocity contour map with varying the particle size distribution. With increasing amounts of small particles, that is, the ratio of large particles to small particles increases, the fluid seepage channel becomes complex, the advantageous seepage path incorporates more twists and turns, and fluid seepage velocity increases significantly, but the overall permeability decreases. Increasing the particle gradation ratio makes the fluid path more dispersed, which is  10 Geofluids not enough to form a larger main seepage path and reduces permeability.

Conclusion
A series of particle-fluid coupling models are constructed using the IB-LBM and Weierstrass -Mandelbrot function.
Numerical simulation results show that the IB-LBM has good numerical stability and convergence in simulating complex fluid seepage flow in fractures with particles, and this method can accurately simulate fluid-particle coupling. In order to explore the seepage behavior of fluid in complex fractures with particles, the effects of the particle size, particle number, fractal dimension, and particle gradation composition on fluid passage are examined. The coupling of fluid particles is analyzed at the microstructure level.
With increasing particle diameter, fluid passage in the channel is blocked by particles, which decreases the seepage area, resulting in a decrease in seepage flux and permeability. With increasing particle number, the average seepage velocity decreases. Fluid velocity will increase around the particles, but this effect is smaller than the overall permeability reduction caused by the decrease in fluid cross-sectional area. With increasing fractal dimension, the roughness of fracture increases, and the complexity of corresponding fluid flow path through the fracture increases. Increasing particle gradation composition makes the fluid channel more dispersed but is not enough to form a large main seepage channel, resulting in a permeability decrease.

u:
Fluid velocity t: Time ν: Kinematic viscosity p: Pressure ρ: Fluid density F: Force term Ma: Mach number μ f : Fluid velocity c s : Lattice speed of sound c i : Discrete velocity f i : Density distribution function τ: Dimensionless mean relaxation coefficient δt: Time step f eq i ðx, tÞ: Equilibrium distribution function F i ðx, tÞ: Discrete force N: Number of physical particles M: Number of discrete points in each particle x i ðtÞ: Particle center R i ðtÞ: Instantaneous particle rotational matrix x r ij ðtÞ: Boundary reference point ξ ij ðtÞ: Particle deviation F ij : Restoration force k: Spring constant F i : Total force T i : Total torque m i : Mass of the ith particle Xðs, tÞ: Displacement of solid particles δ: Dirac function F coll : Collision force Fðs, tÞ: Immersed boundary force density D: Fractal dimension γ: Fractal parameter ϕ n : Fractal angle W: Fractal stochastic index CðtÞ: Fractal governing function F D : Resistance coefficient U in : Fluid incident velocity d: Particle diameter K: Permeability Q: Flow quantity L: Partial path length A: Channel cross-sectional area.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
No potential conflicts of interest were reported by the authors.