On the Application of the Lattice Boltzmann Method to Predict Soil Meso Seepage Characteristics

: In this study, a two-dimensional approach is elaborated to study with the lattice Boltzmann method (LBM) the seepage of water in the pores of a soil. Firstly, the D2Q9 model is selected to account for the discrete velocity distribution of water ﬂ ow. In particular, impermeability is considered as macroscopic boundary condition for the left and right domain sides, while the upper and lower boundaries are assumed to behave as pressure boundaries controlled by different densities. The micro-boundary conditions are implemented through the standard rebound strategy and a non-equilibrium extrapolation scheme. Matlab is used for the development of the related algorithm. Finally, the in ﬂ uence of porosity, permeability, osmotic pressure and other factors is assessed with regard to seepage characteristics and the ensuing results are compared with Darcy ’ s law. The com-putations show that, for ﬁ xed initial conditions, the pore structure has a certain in ﬂ uence on the local velocity of seepage, but the overall state is stable, and the average velocity of each layer is the same. The larger the pore passage is, the faster the ﬂ ow velocity is, and vice versa. For low permeability, the numerical results are consistent with the Darcy's law. The greater the pressure difference between the inlet and outlet of seepage, the greater the seepage rate. The relationship between them is linear (yet in good agreement with Darcy ’ s law).


Introduction
Soil is a very important discontinuous porous medium material in geotechnical engineering, which is composed of solid, liquid, gas and pore. The pores can connect with each other in space, and the water in saturated soil flows under the influence of unbalanced potential energy, which is called seepage [1]. Seepage problem is an important subject of geotechnical engineering, such as the calculation of water inflow (seepage flow) of high-rise building foundation and deep foundation pit drainage [2]; Water conservancy engineering problems such as permeability of embankment materials are considered in the construction of embankment [3]. In daily life and production engineering, the exploitation of oil, groundwater, underground minerals and other energy sources, the prediction and prevention of natural field by common calculation methods. The Navier-Stokes equation can be derived from the discrete equation of the particle location r, the corresponding seepage time t and the particle distribution function z (r, t) in the flow field, so as to describe the basic law of fluid flow in the seepage field from the microscopic point of view [13]. The macroscopic equation of two-dimensional seepage can be written as the mass conservation equation.
Momentum conservation equation: In Eqs. (1) and (2), r is the Lap operator, and u is the volume average flow velocity of the fluid. ρ is the density of the fluid, and n is the porosity of the porous medium. P is the pressure of the fluid. m e is the effective viscosity coefficient. F is the resultant force of the fluid in the porous medium. Without external force, the expression is as follows: In Eq. (3), ν is the viscosity coefficient, which is related to m e but not necessarily equal. F n is the geometric factor of the porous medium. K is the permeability of the porous medium and indicates the difficulty of the fluid passing through the porous medium. Then, F n can be expressed as follows: In Eq. (5), d p is the effective diameter of the solid particles of the porous medium.
The lattice Boltzmann method solves the discrete Boltzmann equations of the particle distribution function at time t and position r to satisfy the law of conservation of mass and the law of conservation of momentum. Thereby, the law of fluid flow in porous medium on the meso level is calculated. BLE equation of discrete lattice can be used to express the time rule of particle distribution function f a (r, t), as shown in Eq. (6): In Eq. (1): f a (r, t) represents the distribution function of particles along the r direction at the time of t; e a represents the discrete speed of particles when they move in the flow field; δ t represents the time used to discrete the whole particle flow field; τ represents the particle relaxation time of Dimensionalization; faeq (r, t) represents the distribution function corresponding to the flow field in equilibrium.

LBM Mathematical Model in the Calculation of Pore Seepage
Generally speaking, the variable of lattice BLE algorithm is more complex. However, lattice LBM algorithm can be simplified by single relaxation discrete model. At present, the single relaxation discrete model (LBGK model) is the most widely used LBM model, especially the DnQb model proposed by Qian et al. [14]. For two-dimensional plane problems, the most commonly used are D2Q7 model and D2Q9 model. For three-dimensional problems, the most widely used models are D3Q19 and D3Q15. The seepage problem in this study is a two-dimensional plane problem, so D2Q9 model is adopted. D2Q9 model divides the seepage field of soil into a square grid with uniform size. A single particle has nine moving squares, each of which is adjacent to eight surrounding nodes, as shown in Fig. 1.
The set of 9 velocity vectors in different directions is as follows: In Eq. (7), c is the lattice velocity in the composition of D2Q9 model. c = δ x /δ t , δ x is the lattice step, and δ t is the time step. In general, the lattice steps in x and y directions are the same, that is, δ x = δ y . According to the law of conservation of mass and the law of conservation of momentum, the relationship between the distribution function of equilibrium state and macro variables of D2Q9 model can be obtained, as shown from Eq. (8) to Eq. (11): In Eq. (3): c s is the lattice velocity, c s 2 = c 2/3 . ω a is the weight coefficient. ρ is the model density. u is the macro variable. f a e a (11) By using Chapman-Enskog expansion method, evolution (Eq. (6)) can be reduced to Navier-Stokes equation of hydrodynamics, and the relationship between viscosity coefficient v of fluid motion and dimensionally unified relaxation time τ can be obtained.

Boundary Treatment of LBM Method
Boundary conditions play an important role in seepage problems. For example, for the steady-state problem, when the seepage in the model is stable, the flow field has nothing to do with the initial conditions, but only with the boundary conditions. For different methods and models, it is very important to choose the right numerical treatment method and the right boundary conditions in seepage simulation. Boundary treatment is a key content in the application of lattice BLE method to seepage simulation, which will have a great impact on the accuracy, calculation efficiency and numerical stability of numerical simulation results [15]. The boundary conditions of seepage simulation model are composed of two schemes: standard rebound scheme and non-equilibrium extrapolation scheme. Among them, the upper boundary (inflow boundary) and the lower boundary (outflow boundary) of the soil are extrapolated in non-equilibrium form, and the left and right impermeable boundary and the soil particle boundary of the model are rebounded in standard form.

Standard Bounce Format
The standard rebound scheme is a kind of heuristic scheme which is often used to deal with the static sliding free wall. The scheme assumes that the particles collide with the non-slip wall, and then the velocity reverses, that is, the velocity remains the same and rebounds in the opposite direction. Fig. 2 is a schematic diagram of the standard rebound format when the lower boundary is fixed.
In Fig. 2, O is the boundary node. f 1 , f 4 , f 5 , f 6 , f 7 , f 8 are the distribution functions of the velocity reversal after the particle collides with the boundary wall. From the fluid node to the boundary node O, it will return along the original path without collision. In Fig. 2, the collision distribution function is shown in Eq. (13):

Nonequilibrium Extrapolation Scheme
According to a new boundary treatment scheme proposed by Guo et al., the distribution function on the boundary node is divided into two parts: equilibrium state and non-equilibrium state, that is, non-equilibrium extrapolation scheme [16]. The equilibrium part is approximately obtained by the definition of boundary conditions, while the non-equilibrium part is determined by non-equilibrium extrapolation. Fig. 3 is a schematic diagram of the non-equilibrium extrapolation scheme. In the figure, A, O and C are boundary grid points, and grid points E, B and D are located in the flow field. f 0 , f 1 , f 2 , f 3 , f 5 , f 6 are the distribution functions of the velocity reversal after the particle collides with the boundary wall.
The non-equilibrium extrapolation method divides the distribution function of boundary nodes into two parts: equilibrium (f eq ) and non-equilibrium (f neq ), which is Eq. (14): As shown in Fig. 3, the equilibrium function can be obtained from the speed and density under macro conditions, while the non-equilibrium part usually has unknown macro physical quantities. In order to simplify the calculation, the corresponding value of point B is used instead. The non-equilibrium extrapolation scheme of boundary O points can be expressed as Eq. (15): The lattice B of evolution direction is as Eq. (16):

The Transformation between Physical Quantities and Lattice Units
In dealing with the problem of physical simulation, there are two ways to connect the physical quantity with the numerical simulation. One is to replace all the physical quantities in the program code in the simulation with the actual physical quantities. The other is to carry out dimensionless treatment of all the physical quantities and ensure that the flow criterion before and after dimensionless treatment is unchanged. The second method will be mainly introduced, that is, how to convert physical units into lattice units in numerical simulation. In order to realize the seepage simulation of water in soil, the actual unit is converted. Between the actual unit and the lattice unit, the relationship between them can be In the D2Q9 model, the basic parameters are length L, evolution time T and kinematic viscosity coefficient v. These are set as lattice units, and the macroscopic physical quantity is L′, v′, T′, and the intermediate conversion unit is Lr, as shown in Eqs. (17) and (18): In Eqs. (17) and (18), the lattice units L and v are used as set parameters, and the macroscopic physical quantities L′ can be obtained from relevant data, so Lr can be derived. According to the Reynolds number Re, the relationship is obtained as shown in Eqs. (19) and (20): In Eqs. (19) and (20), u is the lattice velocity and u′ is the macroscopic velocity. Through the above equations and known quantities, the required unknown quantities can be obtained.

Establishment of Soil Pore Model
The pore structure of soil is complex in space and distributed in disorder. In order to describe the pore structure characteristics of soil, it is assumed that the soil is only composed of solid soil particles and pores, and its spatial distribution can be expressed as the Eq. (16): In Eq. (21), Z (x) is a random variable, its statistical characteristics can reflect the distribution of soil pores, and the statistical average value of ZðxÞ = n. n is the porosity of soil.
In this study, the range of soil model is set in Matlab, and the rand command is introduced to generate a uniform random array N (x). Among them, rand () is between (0, 1). The size of N (x) and porosity n is Fig. 4 shows the structure of soil model with porosity n = 0.5. The black area is the soil particle skeleton, and the white area is the pore.
It can be seen from Fig. 4 soil pore structure that when the porosity is small, the soil pore connectivity is poor, and some of the pores cannot be penetrated. At this time, it is difficult to achieve the purpose of matching the actual application in seepage simulation, and the overall seepage work cannot be achieved. Therefore, in most cases, the method of this study is only applicable to the research of large-scale granular soil, such as sandy soil, while it is not applicable to the research of small porosity, such as finegrained cohesive soil.

Validation Study on Soil Meso Seepage Based on Lattice LBM Method
For the meso seepage model of saturated soil based on the LBM method, the D2Q9 model is used, and the corresponding calculation program is compiled using Matlab. Then, referring to the calculation method in reference [17], the soil seepage is numerically simulated. The weight factor is configured as: x 0 = 4/9, x 1 =1/ 9, x 5 = 1/36, and C s = l/ ffiffi ffi 3 p . The calculation time step is δ 1 = 1 and the grid step is δ x = δ y = 1. r oh_in = 1.001 is taken from the entrance and r oh_out = 0.999 is taken from the exit. The viscosity coefficient is set to 1.28 cst and the calculation area is 200 mm × 200 mm. The unit is the fragment. This numerical simulation uses mixed soil samples. In addition, the particles are classified as shown in the table. Using a self-editing program, the soil particle gradation is input to obtain randomly generated soil samples. The calculation results of different grid sizes fluctuate greatly. It can be seen that the calculation result is sensitive to the grid size. As shown in Tab. 1, when the grid size is 1~0.5 mm, the fluctuation of the calculation result is relatively small. Therefore, the model uses 1 mm as the grid distance.
To realize the LBM method to simulate soil see page, a 20 cm × 20 cm area is selected for soil reconstruction. Water seeps into the model from top to bottom. The initial velocity is 0.062 cm/s. The porosity n = 0.5 and 0.6. The time step δ 1 = 1. The lattice step δ x = δ y = 1. The upper and lower boundaries adopt non-equilibrium extrapolation scheme. The left and right boundaries are impermeable surfaces. Boundary rebound is the standard rebound format. The model is shown in Fig. 5.
According to Section 2.4, the unit transformation is carried out, and the macroscopic physical quantities used in this simulation and the transformed lattice units are shown in Tab. 2:   To facilitate calculations, the following basic assumptions are made on the soil seepage model during validation study: 1) The water seepage velocity in the soil is small, all of which are laminar; 2) The soil is fully saturated, and no gas phase affects water seepage; 3) All the water in the pores of the soil is free water, and there is no crystal water. During the seepage process, the free water content in the unit volume of soil remains unchanged; 4) For the seepage, the external factors such as gravity and the deformation of the soil body are not considered, that is, the shape of the pore channel in the soil body is unchanged; 5) The left and right boundaries of the model are impervious, that is, when rebound occurs at the boundary, there is no loss of speed and energy.
The calculation flow is shown in Fig. 6: According to Fig. 6, in Matlab, first, the seepage model of soil micro-structure is introduced, and the initialization parameters of the model is set, such as flow velocity v and density ρ. Then, the particles collide. For the seepage model, boundary conditions are applied. Then, the next seepage velocity v, soil density ρ and lattice particle distribution function are calculated to determine whether the model has converged. If there is no convergence, the process needs to be repeated until the convergence and calculation are completed. In this study, the relative error of the macro velocity of two adjacent time steps is used to judge the convergence of the calculation result. The basis for judgment is as follows: In Eq. (22): Given a small quantity m, when err < m, the calculation converges and the operation stops.

Darcy's Law of Seepage
According to the LBM method, calculation process and convergence basis proposed in this study, the relationship between porosity n, permeability K, osmotic pressure P and seepage velocity under pressure is discussed respectively. Moreover, the relevant calculation results are compared with the distribution law of Darcy seepage law. In the case of laminar flow, the effect of gravity is not considered. Darcy's law equation of soil seepage [18] is as Eq. (25): In the Eq. (25), u is the flow velocity of seepage and v is the viscosity coefficient of fluid.

Example Calculation Results of Meso Seepage Model of Soil
According to the calculation flow and model in Figs. 5 and 6, the velocity distribution image after seepage stabilization of porous n = 0.5 (a) and porous n = 0.6 (b) are obtained in Fig. 7. Fig. 8 is the vector diagram of local seepage velocity in the seepage field when n = 0.6. From Fig. 7, it can be seen that comparing two seepage images with different porosity, after the particle seepage in the fluid is stable, the connectivity of the pore channel determines the seepage velocity of the internal pore particles in the flow field. In the soil model, the flow velocity of seepage is uniform as a whole. In the channel with good pore connectivity, the flow velocity of seepage is faster, and in the area with smaller channel, the flow velocity of seepage is smaller. It can be seen from Fig. 8 that in the channel with good pore connectivity, the velocity vector direction in the seepage flow field is relatively constant, and all flows in the constant direction. In the channel with poor pore connectivity, the velocity distribution is disordered. It can be seen here that it is similar to the actual seepage situation, which proves that the LBM method can be used for seepage research.

The Influence of Porosity on Seepage Characteristics
In order to study the influence of soil porosity on the seepage characteristics, the seepage inlet density is set to ρ in = 1.1, the seepage outlet density ρout = 0.9, and the permeability K = 0.1. Under the same boundary conditions, the lattice LBM algorithm proposed in this study is used to calculate the distribution of seepage velocity in the internal region of the model when the porosity n of the model is 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 and 1.0 respectively. The results are compared with those of Darcy's law, as shown in Fig. 9.   From Fig. 9, it can be seen that the less dense the soil is, the larger the porosity of the soil is, the worse the connectivity of the pore channel is, and the smaller the corresponding seepage rate is. The smaller the pore of soil, the better the connectivity of pore channel, and the larger the corresponding seepage velocity. The porosity and the flow velocity are in inverse proportion.

Influence of Permeability on Seepage Characteristics
In order to study the influence of soil permeability on seepage characteristics, the seepage inlet density is set as ρ in = 1.1, the seepage outlet density as ρ out = 0.9, and the permeability as n = 0.5. Under the same boundary conditions, through the lattice LBM algorithm proposed in this study, the seepage velocity distribution in the internal area of the model is obtained when the permeability K of the soil microstructure seepage model is 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1.0, 2.0, 5.0 and 10.0, respectively. The results are compared with those of Darcy's law, as shown in Fig. 10.
It can be seen from the Fig. 10 that the larger the permeability K is, the greater the seepage rate is. When the permeability K < 0.5, in this study, the calculation results of seepage velocity are in good agreement with Darcy's seepage law. When permeability K ≥ 0.5, there is a big deviation between them. Moreover, the larger the permeability is, the larger the deviation is.

Influence of Seepage Pressure on Seepage Characteristics
In order to study the influence of seepage pressure on seepage characteristics, porosity is set as n = 0.5, permeability K = 0.1. Under the same boundary conditions, through the lattice LBM algorithm proposed in this study, the seepage velocity in the inner region of the model when the difference ρ in -ρ out between the density at the entrance and the density at the exit of the seepage model of soil meso-structure is 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 and 1.0 is calculated, respectively, and the results are compared with the calculation results of Darcy's law. The results are shown in Fig. 11.
It can be seen from Fig. 11 that the calculation results of this study are consistent with Darcy's seepage law, and the seepage velocity increases with the increase of the pressure difference between the inlet and outlet of seepage, and the two change linearly. Therefore, the variation law of seepage velocity is consistent with the actual situation, and the calculation results of this study are in good agreement with Darcy's law.

Discussion
As a porous medium material, soil contains many pore channels with different sizes and specifications. It is an important carrier and material in production and life. The internal structure in the soil is extremely complex and even exhibits greatly obvious multi-scale structural characteristics. The seepage of water in the soil will cause many environmental and engineering problems, which have a great impact on human life and the environment. It is closely related to engineering fields such as construction, transportation, mining, water conservancy, and petroleum. Therefore, based on the LBM method, the distribution and variation of the seepage field velocity in the saturated soil with a certain porosity are explored. From a meso perspective, the seepage mechanism of the pores in the soil under the action of external water flow is investigated. It can be seen from the velocity distribution and velocity vector distribution in the experimental soil seepage field obtained by numerical simulation that the established model simulates the seepage situation similar to the actual seepage situation. It is proved that the LBM method can be used for seepage research. Wang et al. [19] proposed a three-dimensional coupled bonded particle and lattice Boltzmann method (BPLBM) with an immersive moving boundary scheme. It is used to investigate the erosion process of soil particles in a particle filter in an earth dam. Three particle filters with different representative size ratios are simulated. The numerical results are consistent with empirical standards, which is consistent with the results of this research.
The influence of different pore structure, porosity, permeability, osmotic pressure and other factors on seepage characteristics is explored. It proves that under the same porosity and other initial conditions of seepage, only the pore structure is changed, and the local velocity of seepage changes, but the overall state is stable. The average flow velocity of each layer is equal, and is similar to the actual seepage situation. Under the premise of controlling other conditions unchanged, only the porosity is changed (the pore structure also changes with it). The velocity of seepage is faster in the area where the pore channel is larger, and the velocity is smaller in the area where the pore channel is smaller. It is consistent with the calculation results obtained by Darcy's law. It is ensured that other initial conditions remain unchanged, and only the permeability is changed. The seepage results indicate that the greater the permeability, the greater the seepage rate. When the permeability is small, the numerical calculation results in the research are in good agreement with Darcy's law. When the permeability is large, the deviation between the two is large. When considering the osmotic pressure, since the porosity is the same, the greater the pressure difference between the inlet and outlet of the seepage, the greater the seepage rate. There is a linear relationship between the two. Also, the numerical calculation results in the research have good consistency with Darcy's law. Wang et al. applied the lattice Boltzmann method to study the seepage behavior of large pore soils with good vegetation. It is proved that along the depth direction, the larger the porosity of the large pores, the faster the average seepage velocity of the profile. In closed or poorly connected large pores, the conductivity of water is greatly small, and the flow velocity is almost equal to zero [20]. It is almost consistent with the results of this research. The pore structure is generated by the four-parameter random growth method. When m ≤ 0.5, the connectivity of the pores is relatively poor. Therefore, based on the lattice LBM method, the soil method and the control variables are reconstructed, indicating the diversity of factors affecting the permeability of the random microscopic soil structure. The lattice BLE method has a good application to the theoretical study of the seepage direction, and it can make certain reference value and guiding significance to the seepage phenomenon in actual production and life.

Conclusion
In this study, based on the lattice BLE method, the D2Q9 model, the standard rebound scheme and the non-equilibrium extrapolation scheme are used, and the two-dimensional micro soil seepage characteristics are discussed through the random configuration method. The experiment of this study shows that the single relaxation method of lattice BLE method can simulate the seepage problem of micro soil very well, and has a good effect. Based on the lattice BLE method and the soil reconstruction method, the seepage simulation calculation is carried out and compared with the calculation results of Darcy's law. Among the factors such as pore structure, porosity, permeability and osmotic pressure, if one of them is changed and the other factors remain unchanged, the numerical results are in good agreement with Darcy's law. Therefore, the calculation in this study is suitable for cohesive soil with large porosity, but not for sand with small porosity. In this study, based on the lattice LBM method, the reconstruction method and control variables show that the random micro soil structure affects the diversity of permeability factors. The lattice BLE method has a good application in the theoretical research of seepage direction, and it can provide a certain reference value and guidance for the seepage phenomenon in actual production and life. The premise of this study is to assume that the soil is completely saturated and the water flow is always in the laminar state during the seepage process. In the actual production and life, the completely saturated soil is relatively rare. The water in the soil is also divided into combined water and free water. In the process of seepage, there is also the problem of energy loss, and the effect and influence of external forces should be considered. Therefore, in the future research, it is necessary to consider the influence of these factors, in order to be more closely linked with the reality, so as to achieve the purpose of theory guiding practice.