**1. Introduction**

Sand production is one of the common problems in the oil field during the production process, especially in weakly consolidated or unconsolidated sandstone reservoirs. The stress and pore pressure near the perforation cavity would be redistributed and the rock would change from an equilibrium state to another after the perforation process. The stress state around the perforation cavities region changed due to the considerable force exerted on the rock by the fluid as the fluids flow into the wellbore. Thus, rock failure will occur in the weakly consolidated or unconsolidated reservoir consequently. In this condition, sand particles will be produced along with the reservoir fluids, which lead to the evolution of perforation cavities and deteriorate the formation condition. A lot of money was spent to deal with problems including borehole instability, casing collapse, and well cleaning. Research shows that almost 30% of these sandstone reservoirs may produce sand easily [1]. Besides this, it is reported that the weakly consolidated or unconsolidated reservoir contains 70% of oil and gas in the world [2,3]. However, the sand production prediction is full of uncertainty. Therefore, accurate simulation and prediction of the sand production process are very important and helpful in sand control and enhancing oil recovery.

The commonly used methods to predict sand production consist of experiments and numerical simulations. Terzaghi firstly used a box with a movable lid at the bottom to conduct the sand production experiment and laid a foundation for the later laboratory experiments [4]. Up to now, experiments under more complex conditions such as uniaxial or tri-axial stress environments with axial and multiphase fluid injections can be conducted [5–8]. In these experiments, different categories of influencing factors including the characteristics of fluids, the state of the formation, the effect of confining pressure, and the effect of outlet size are discussed. However, the sand production experiment needs some assumptions since the scale is small and its result cannot be used in field practice due to the boundary effects. Besides this, it is hard to obtain the in-situ rock sample in some reservoirs that are weakly consolidated or unconsolidated.

Various numerical methods considering the plasticity and micromechanics of the rock, erosion effect, and the mechanical mechanisms due to waterflooding have been developed to accurately predict sand production. These methods consist of continuum method and discrete method. The continuum method considers the rock a continuous material and is widely used in sand production prediction because of its simplicity. It includes models based on plasticity theory, mixture theory, and the micromechanics approach. Methods based on plasticity theory considered the yielding behavior of the rock as the result of tensile and shear failures due to the depletion of well pressure [9]. The effect of hydrostatic stresses due to fluid flow is considered and many models are proposed to predict sand production [10]. Methods based on mixture theory considered the effect of hydrodynamic forces on the sand production process [11]. In these methods, the erosion phenomenon was studied and considered as the cause of the rock failure [12]. More influencing factors including flow rate, viscosity, density, grain size, and cavity height are employed to propose the model [13]. The micromechanics approach combined the poro-mechanical and erosion models extensively rather than employing mechanical yielding models and the erosion models independently. Two semi-analytical models, namely SPAM and SLAM model, were presented based on this concept and study the sand production due to the sheer pressure on the rock [14]. However, the continuum methods cannot reflect the real behavior of the rock such as nonlinear characteristics of mechanical response due to its simplicity. The discrete method treats the rock as a combination of many individual sand particles that have their characteristics. O'Connor et al. firstly used the DEM method to simulate sand production [15]. Based on their research, different shapes of particles such as the n-sided polygon were considered [16]. The growth of a micro-failure leading to rock mass failure is concluded by employing PFC2D for simulating the sanding process in hollow cylinder tests with fluid flow [17]. The coupled LBM-DEM method precisely presented the solid-fluid interface and enhanced the reliability of the prediction of sand production and the lattice Boltzmann method (LBM) has been proved to be the most appropriate method in investigating the fluid flow mechanism in the porous media [18–22]. However, these methods are mostly developed on two dimensions and the bond between the particles is not fully considered in the simulation. Besides this, the characteristic of force chain network in the sand production process is not studied in these research works.

Compared with particle methods, it is more appropriate to use the fluid-solid coupling method since sand production is a process in which fluid and solid interact with each other. We adopted the coupled lattice Boltzmann method and discrete element method (LBM-DEM) in this paper since it is effective to deal with the complicated boundary conditions and it is more feasible in pore-scale (mesoscale) fluid flow. Besides this, the CFD-DEM method is computationally expensive when dealing with fine mesh and solving the nonlinear Navier–Stokes equation in the pore scale. Then we validate the method by simulating the sedimentation of a single particle. Finally, we conduct the simulation of the sand production process and analyze the force network evolvement. Absolute and relative permeability before and after the sand production process are calculated and the effect of injection flow rate, cementation strength, and confining pressure are investigated. Compared with previous studies, we conducted our simulation on three dimensions and fully considered the bond between particles. Besides this, the complex network analysis of the force chain network evolution in the sand production process, which was not studied in previous research, is fully discussed in this paper.

#### **2. Methodology**

The coupled LBM-DEM method was firstly proposed by Cook, et al. [23]. This method was further improved by a great number of researchers subsequently [24–27]. The basic principles of each method and how the two methods are coupled with each other are given as follows.

#### *2.1. Bonded DEM*

DEM was developed to deal with rock mechanics problems in the 1970s [28]. It has been regarded as a useful tool for analyzing the behavior of materials with a discontinuous characteristic. In this method, the porous media are considered a combination of separate particulates. Due to the simplicity and numerical efficiency, we often use spherical particles.

In this method, Newton's second law is used to describe the motion of each particle. The equation that governs the translation motion of a particle can be described as [28]

$$m\frac{d\mathbf{v}}{dt} = \mathbf{F}\_b + \mathbf{F}\_c + \mathbf{F}\_h \tag{1}$$

where *m* and **v** denote the mass and the velocity of the particle, respectively. **F***<sup>b</sup>* indicates the body force, **F***c* refers to the total contact force which is calculated by summing up all the contact force on the particle. **F***<sup>h</sup>* represents the hydrodynamic force exerted by the flowing fluid. The equation that governs the rotation of a particle is [28]

$$I\frac{d\mathbf{w}}{dt} = \mathbf{T}\_c + \mathbf{T}\_h\tag{2}$$

where *I* is the moment of inertia of the particle, and ω is the angular of the particle. **T***<sup>c</sup>* refers to the total torque generated by the contact force and **T***<sup>h</sup>* indicates the hydrodynamic torque exerted by the flowing fluid.

In the DEM simulation, particles moved freely in the domain. We defined that contact exists when the distance between two particles centers is shorter than the summation of the particle radius. Figure 1 illustrates two particles in contact. The difference between the summation of the radius and the particle center distance was called overlapping length. By using the Hertz contact model, the contact force can be calculated by [28]

$$F\_n = k\_n \delta\_n^{3/2} \tag{3}$$

where *kn* denotes normal stiffness, which is influenced by the elastic modulus of the material. *δ<sup>n</sup>* refers to the overlapping length.

**Figure 1.** The sketch of two contact particles.

For each particle at each time step, the accelerations can be computed by using Equations (1) and (2). Particle velocity and position would be calculated by employing equations of motion consequently. In this way, the information of any single particle at any time step can be obtained and it is convenient for us to do further analyses.

For the cementation of the rock, we used the contact bond model to simulate the bond characteristics [29,30]. In this model, the grain particles are approximated as spherical particles and the cementation bonds are considered as sticks joining the two adjacent particles. The bond element can be regarded as two springs with normal and shear stiffness acting between the contacting particles. The contact bond model can be shown as [31]

$$F\_{\mathfrak{n}}^{b} = \begin{cases} \begin{array}{c} k\_{\mathfrak{n}}^{b} \delta \left( F\_{\mathfrak{n}}^{b} \leq F\_{\text{max}} \right) \\ 0 \text{ (} F\_{\mathfrak{n}}^{b} > F\_{\text{max}} \text{)} \end{array} \tag{4}$$

where *k<sup>b</sup> <sup>n</sup>* is the normal stiffness of the cement, *F*max is the critical tensile force. The magnitude of the critical force can be calculated as

$$F\_{\text{max}} = \sigma\_{\text{max}} S\_c \tag{5}$$

where *σ*max is the bond strength and the *S*<sup>c</sup> is the cross-sectional area of the bond element. The bond strength of the bond element can be determined by conducting a uniaxial tensile test [32,33]. The bond element can be regarded as a beam and the radius of the cross-section *rb* = (*r*<sup>1</sup> + *r*2)/2, where *r*<sup>1</sup> and *r*<sup>2</sup> are the radii of the bonded particles and the cross-sectional area can be calculated subsequently [34]. If the tensile force acting on the cement exceeds or equals the normal bond tensile strength, the bond breaks and the normal contact force should be zero.

#### *2.2. Lattice Boltzmann Method*

Fluid flow is often modeled by Navier–Stokes equations, and LBM is a potential tool for the solution of these equations. The BGK model is most widely used in LBM modeling and streaming and collision is given by [35]

$$f\_d(\mathbf{x} + \mathbf{e}\_d \Delta t, t + \Delta t) - f\_d(\mathbf{x}, t) = -\frac{\Delta t}{\pi} \left[ f\_d(\mathbf{x}, t) - f\_d^{eq}(\mathbf{x}, t) \right] \tag{6}$$

where **x** is the position in the lattice space, *τ* is the collision relaxation time, Δ*t* is the time interval, **e***<sup>a</sup>* refers to the lattice velocity, *a* indicates the discretized direction, *f* and *f eq* denote the local particle distribution function and equilibrium distribution function. In the equation, the left side is the streaming part and the right side is the collision part.

Collision of the fluid particles is regarded as a relaxation towards a local equilibrium and the equilibrium distribution function *f eq* is defined as [35]

$$f\_a^{eq}(\mathbf{x}) = w\_a \rho(\mathbf{x}) \left[ 1 + 3 \frac{\mathbf{e}\_d \cdot \mathbf{u}}{c^2} + \frac{9}{2} \frac{(\mathbf{e}\_d \cdot \mathbf{u})^2}{c^4} - \frac{3}{2} \frac{\mathbf{u}^2}{c^2} \right] \tag{7}$$

where the weights *wa* are 1/3 for the rest particles (*a* = 0), 1/18 for *a* = 1~6, and 1/36 for *a* = 7~18, and *c* is the basic speed on the lattice given by *c* = Δ*x*/Δ*t*, Δ*x* denotes the lattice spacing. The configuration of D3Q19 lattice cell is presented in Figure 2.

The distribution function represents the frequency of occurrence in a certain direction and the summation of frequencies of all directions can be considered as fluid densities. Therefore, the macroscopic fluid density is [35]

$$\rho \quad = \sum\_{a=-0}^{18} f\_a \tag{8}$$

**Figure 2.** The configuration of D3Q19 showing the 19 discrete velocities.

The macroscopic velocity **u** is an average of the microscopic velocities **e***a* weighted by the directional densities *fa* and can be expressed as [35]

$$\mathbf{u}\_{\parallel} = \frac{1}{\rho} \sum\_{a=-0}^{18} f\_a \mathbf{e}\_a \tag{9}$$

Equation (9) allows us to connect the discrete microscopic velocities in LBM with the continuous macroscopic velocities representing the fluid's motion.

Viscosity and pressure can be described by [35]

$$\nu = c\_s^2 (\tau - \frac{1}{2}) \Delta t \tag{10}$$

$$
\rho = \rho \mathbf{c}\_s^2 \tag{11}
$$

where *c*s is the sound velocity given by [36]

$$c\_s = \sqrt{\frac{c^2}{3}}\tag{12}$$

Shan-Chen model, also known as the pseudo-potential model, is one of the most widely used multiphase Lattice Boltzmann models. In this paper, we employed this method to simulate two-phase flow in the porous media and obtain the relative permeabilities in different stages of sand production.

When the fluid consists of multiple components, a pseudo-potential model can be established by incorporating interaction force between fluid components. For two-phase flow, two distribution functions represent the flow of wetting phase and non-wetting phase and the evolution function of the *σ*th component can be described as [37]

$$f\_a^{\sigma}(\mathbf{x} + \mathbf{e}\_d \Delta t, t + \Delta t) - f\_a^{\sigma}(\mathbf{x}, t) = -\frac{\Delta t}{\tau\_{\sigma}} [f\_a^{\sigma}(\mathbf{x}, t) - f\_a^{\sigma(c\eta)}(\mathbf{x}, t)] \tag{13}$$

where *f <sup>σ</sup> <sup>a</sup>* denote the local particle distribution function of the *σ*th component, *τσ* is the relaxation time of the *σ*th component and the corresponding kinematic viscosity can be expressed as [37]

$$\nu\_{\sigma} = c\_s^2 (\tau\_{\sigma} - \frac{1}{2}) \Delta t \tag{14}$$

In order to simulate two-phase flow in porous media, many researchers proposed different forms to incorporate the interaction force between two fluids in the fluid flow [38–42]. In this paper, we incorporate the interaction force in the collision process by adding a term into the velocity field and the macroscopic velocity of the *σ*th component can be calculated by [43]

$$\mathbf{u}\_{\sigma}^{\alpha \eta} = \tilde{\mathbf{u}} + \frac{\mathbf{r}\_{\sigma} \mathbf{F}\_{\sigma}}{\rho\_{\sigma}} \tag{15}$$

$$
\tilde{\mathbf{u}}^{\prime} = \frac{\sum\_{\sigma} \frac{\rho\_{\nu} \mathbf{u}\_{\sigma}}{\mathbf{\tilde{r}}\_{\sigma}}}{\sum\_{\sigma} \frac{\rho\_{\nu}}{\mathbf{\tilde{r}}\_{\sigma}}} \tag{16}
$$

where **<sup>~</sup> u** is the compound velocity of all components without considering additional forces. **F***<sup>σ</sup>* denotes the total force exerted on the *σ*th component, which includes the fluid–fluid interaction **F**f, the fluid–solid interaction **F***ads*, and the external force **F**<sup>e</sup> [44].

Therefore, the whole fluid velocity can be defined as [45,46]

$$\mathbf{u} \,\, = \, \frac{1}{\rho} \left( \sum\_{\sigma} \rho\_{\sigma} \mathbf{u}\_{\sigma} + \frac{1}{2} \sum\_{\sigma} \mathbf{F}\_{\sigma} \right) \tag{17}$$

where *ρ* is the total density of all the fluid components and can be expressed as [45]

$$
\rho = \sum\_{\sigma} \rho\_{\sigma} \tag{18}
$$

The fluid–fluid interaction only considers the interaction between the nearest neighbors and can be expressed by [44]

$$\mathbf{F}\_{\mathbf{f}}^{\tau}(\mathbf{x},t) = -\boldsymbol{\psi}^{\sigma}(\mathbf{x},t) \sum\_{\sigma'}^{S} G\_{\sigma\sigma'} \sum\_{a} \omega\_{i} \boldsymbol{\psi}^{\sigma'}(\mathbf{x} + \mathbf{e}\_{a} \Delta t, t) \mathbf{e}\_{a} \tag{19}$$

where *Gσσ'* denotes the Green function reflecting the strength of the intercomponent potential between the *σ*th component and the *σ* th component, *ψσ*(**x**, *t*) is the effective density of the *σ*th component and is commonly taken as the fluid density.

At the fluid-solid interface, the solid wall can be treated as a component with constant density. Therefore, the interaction force between the fluid and solid can be written as [44]

$$\mathbf{F}\_{ads}^{\sigma}(\mathbf{x},t) = \left. -\boldsymbol{\upvarphi}^{\sigma}(\mathbf{x},t) \sum\_{\sigma}^{S} \mathbf{G}\_{\sigma}^{\sigma} \sum\_{a} \boldsymbol{\omega}\_{i} \boldsymbol{\upvarphi}(\boldsymbol{\uprho}\_{\sigma}) \mathbf{s} \left( \mathbf{x} + \mathbf{e}\_{a} \boldsymbol{\Delta t}, t \right) \mathbf{e}\_{a} \tag{20}$$

where *G<sup>σ</sup> <sup>s</sup>* denotes the Green function reflecting the strength of potential between the fluid and the solid. *s*(**x** + **e**aΔ*t*, *t*) is a function that is equal to 1 when the node is solid and zero when the node is fluid. *ρ<sup>w</sup>* is the density of the solid wall.

#### *2.3. LBM-DEM Coupling*

The moving boundary condition is very important for the coupled LBM and DEM methods. Researchers have come up with numerous moving boundary conditions for coupling. In this study, we employ a method proposed by Noble and Torczynski [25,26,47]. In the coupling method, the collision operator in the LBM equation is expressed with the solid fraction *γ* in each node (Figure 3). Thus, Equation (6) can be described as [18]

$$f\_d(\mathbf{x} + \mathbf{e}\_d \Delta t, t + \Delta t) = f\_d(\mathbf{x}, t) - \frac{1}{\tau}(1 - \beta) \left[ f\_d(\mathbf{x}, t) - f\_d^{eq}(\mathbf{x}, t) \right] + \beta f\_d^m \tag{21}$$

where *β* is a weighting function which depends on the solid fraction *γ* in each node [18],

$$\beta = \frac{\gamma(\tau - 0.5)}{(1 - \gamma) + (\tau - 0.5)} \tag{22}$$

**Figure 3.** Noble and Torczynski's scheme.

*f m <sup>a</sup>* is the collision term that represents the bounce-back of the non-equilibrium part of the distribution function. It can be described by [18]

$$f\_a^{m} = f\_{-a}(\mathbf{x}, t) - f\_a(\mathbf{x}, t) + f\_a^{eq}(\rho, \mathbf{v}\_b) - f\_{-a}^{eq}(\rho, \mathbf{v}) \tag{23}$$

where −*a* is the direction opposite of *a.* The hydrodynamic forces and torque acted on a particle can be calculated as [18]

$$\mathbf{F}\_{fluid} = \frac{h^2}{\Delta t} \left[ \sum\_{n} \left( \beta\_n \sum\_{a} f\_a^m \mathbf{e}\_a \right) \right] \tag{24}$$

$$\mathbf{T}\_{fluid} = \frac{h^2}{\Delta t} \left[ \sum\_{n} (\mathbf{x}\_n - \mathbf{x}\_c) \times \sum\_{n} \left( \beta\_n \sum\_{a} f\_a^m \mathbf{e}\_a \right) \right] \tag{25}$$

where *n* is the number of nodes covered by a particle. The process of coupled LBM and DEM is shown in Figure 4.

#### *2.4. Method Validation*

In this section, we will verify the boundary condition for a solid that moves across the lattice by simulating a single sphere settling into a low particle Reynold number. For this case, it has a specific analytical expression for both the force and the torque [3]. The force for this case is given by

$$F = \frac{-6\pi\rho vrv}{1 - 0.6526(r/l) + 0.1475(r/l)^3 - 0.131(r/l)^4 - 0.0644(r/l)^5} \tag{26}$$

where *ρ* is the fluid density, *ν* is the fluid kinetic viscosity, *v* is the velocity of the particle, *r* denotes the particle radius, and *l* represents the distance between the particle center and the wall.

**Figure 4.** The cycle of computation for coupled LBM and DEM.

Our simulation was set up similarly to that by Strack et al. (2007). Two fixed walls 2*l* distance apart were employed in the transverse direction. The width and depth of the computational domain were set to a relatively larger value of 8*l* to decrease the dependence of the simulation results on the domain size. The walls were modeled with a bounce-back boundary condition, and periodic boundary conditions were adopted in the other two directions. A particle was initially located in the computational domain with a distance of *l* to both the left and the top. The particle diameter was covered by 10 lattices, so the computational domain was 40 × 160 × 160 lattices in the simulation. No body force was applied to the fluid and the particle settled under the influence of gravity. Due to the buoyancy force exerted on the particle, the gravitational force on the particle was considered based on the density difference between the particle and the fluid. The number of CPU cores is 4 in our simulation.

Figure 5a depicted the velocity of the particle and the force acting on the particle in the simulation. The velocity was normalized by the final settling velocity. The final setting velocity denotes the particle velocity in the steady state, which is defined by the relative error of settling velocity between n and n − 1 time steps smaller than the convergence limit <sup>1</sup> × <sup>10</sup>−6. The force was normalized by the analytical values based on the final settling velocity and the vertical displacements were normalized by the lattice spacing. From the figure, we can see that the velocity of the particle and the force exerted on the particle changed synchronously. When the particle reached the final settling velocity, the force exerted on the particle did not change anymore as well. The particle reached an equilibrium state and kept settling at a constant speed.

We studied the computational error of this method and the cumulative error of simulations under different Reynold numbers and lattice spacing can be seen from Figure 5b. The cumulative error was calculated by

$$\varepsilon = \sum \left| \frac{F\_{\text{sim}} - F\_{\text{a}}}{F\_{\text{a}}} \right| \tag{27}$$

where *Fsim* is the simulation results calculated by the LBM-DEM coupling method, Fa denotes the analytical method. As we can see from the figure, the simulation results become more and more accurate with the decrease of the lattice spacing. This is because small lattice spacing means more lattices covered on the particle, which will reduce the oscillation. It also can be spotted that the cumulative computational error is larger when the Reynold number is larger. A larger Reynold number indicates larger particle velocity and computational error is accumulated due to the update of the particle position and recalculation of the solid fraction. Therefore, it should be noticed the displacement of each step for particles should be shorter than the lattice spacing when employing this method to ensure the accuracy and stability of the computation.

**Figure 5.** (**a**) Normalized velocity and Force for a spherical particle settling between infinite walls. (**b**) Cumulative computational error under different Reynold numbers with different lattice spacing.

#### *2.5. Model Setup*

In this section, we used a 3D LBM-DEM approach to simulate the sand production process. The cross-section of our model for the simulation of the sand production process from a cemented reservoir is shown in Figure 6. Simulations were conducted to investigate the migration of sand particles in different circumstances. In this simulation, the fluid flows through a porous media from one end, which causes sand production at the other end through a small outlet. As it is shown in Figure 6, the sand particles were generated in the simulation region. The total length of the region was 60 mm and the length of the porous media was 20 mm. Both the width and height of the computational domain were 20 mm. The outlet of the simulation region is a circle in the middle of the right plane and its diameter equals 4 mm.

**Figure 6.** Cross-section of the simulation region.

Boundaries in both the y-direction and z-direction for the fluid were set as periodic boundaries. The fluid was assumed to flow into the computational domain with a specific velocity and a zero-pressure boundary condition was employed at the right boundary of the computational domain. Particles were allowed to leave the simulation domain from the small outlet. Body force was ignored in the simulation. The gravitational force on the particle was set according to the density difference between the particle and the fluid. The fluid density was 1140 kg/m3, and the dynamic viscosity was 1.14 mPa·s. The simulation parameters of sand particles can be seen in Table 1.

**Table 1.** Parameters for the sand particles in the simulation.


First, we conduct a simulation with fluid flow through the porous media to study the change of force during waterflooding. Then, absolute and relative permeability before and after the sand production process are calculated and the changes in the properties of porous media are studied. Finally, the effect of injection flow rate, cementation strength, and confining pressure are investigated.

#### **3. Results and Discussion**

#### *3.1. Force Chain Network Analysis*

Figure 7 shows the force chain of the porous media at different time steps after a period of waterflooding. From the figure, we can see the evolution of the force chain with the fluid injection. Comparing the force chain circled by the blue line in Figures 7a and 7b, we can spot strong force chain rupture during this process and the sand particles detach from the matrix and migrate through the porous media. Besides this, weak force chains circled by the red line in Figure 7a transform into strong force chains, which can be clearly recognized in Figure 7b. During sand production, the force chain network is reorganized continuously until it reaches a steady state.

**Figure 7.** Force chain of the porous media. ((**a**) t = 0.05 s; (**b**) t = 0.1 s).

The force chain network can be analyzed by using complex network theory. The complex network theory is one of the approaches to describe a granular system. Many researchers used it to study granular packing in the past few decades [48,49]. Since the unconsolidated rock can be regarded as sand packing, it is appropriate to use a contact network to clarify the intrinsic mechanisms in the sand production process.

In our view, the network structure is full of randomicity and we can hardly see any principles in it. However, the certainty of the network structure is often hidden in the randomicity and we can reveal it by studying the statistical properties. Therefore, the description and study of the statistical properties of complex networks are of great significance for the development and application of complex network theory. In order to effectively characterize the mathematical and statistical characteristics of complex network structures, various scholars have proposed many methods and tried to use the network parameters to describe the statistical characteristics of complex networks. The three most basic and important parameters are the degree distribution, the mean shortest-path distance, and the clustering coefficient.

The adjacency matrix is often used in describing the information of a complex network. For the weighted and undirected network, the weighted adjacency matrix **W** can be expressed as [50]

$$\mathbf{W}\_{ij} = \begin{cases} w\_{ij\prime} & \text{if there is an edge between node } i \text{ and node } j \\ 0, & \text{otherwise} \end{cases} \tag{28}$$

where *wij* is the corresponding weight between node *i* and node *j*. For the rock system, *wij* is the force between two particles and the adjacency matrix should be symmetric.

In the weighted network, the connectivity of nodes is also worth noticing regardless of how strongly two particles interact with each other. Therefore, the binary adjacency matrix **A** which only focuses on the link between nodes can be expressed as [50]

$$\mathbf{A}\_{ij} = \begin{cases} 1, & \text{if } w\_{\text{ij}} \neq 0 \\ 0, & \text{otherwise} \end{cases} \tag{29}$$

The degree of a node can be described as the number of edges that are attached to it in an undirected network. Therefore, the degree of the node *i* can be described as [50]

$$d\mathbf{c}\_i = \sum\_{j=-1}^{N} \mathbf{A}\_{ij} \tag{30}$$

where *N* is the number of nodes in the network. Therefore, the mean degree of the network can be computed as [50]

$$
\langle d\varepsilon \rangle \, \, \, \, \, \, \frac{1}{N} \sum\_{i} d\varepsilon\_{i} \, \, \, \tag{31}
$$

For the porous media system, it is easy to find that the node degree is exactly the coordination number of the particle from its definition. Since the degree of nodes in the network follows a certain probability distribution, the degree distribution is an important topological factor to describe and analyze complex networks. The degree distribution function *p*(*k*) represents the ratio of the number of nodes with degree equals to k to the number of all nodes in the network, that is, the proportion of nodes with a degree of k in the network.

The mean shortest-path distance *L* can be defined as the average of the distance between any two nodes and can be calculated by [50]

$$L = \frac{2}{N(N-1)} \sum\_{i,j, i \neq j} d\_{ij} \tag{32}$$

where *dij* denotes the shortest distance between node *i* and node *j*. It should be noticed that the distance should be calculated along edges in a path. The mean shortest-path distance, also known as the characteristic path distance of a network, can be used to measure the efficiency of information transfer between network nodes.

The local clustering coefficient *Ci* can be defined as the number of triangles containing *i* divided by the number of triples centered at node *i* [51,52]. It can be expressed as [50]

$$\mathbf{C}\_{i} = \begin{cases} \frac{\sum\_{\bar{h}} A\_{\bar{h}\bar{i}} A\_{\bar{i}\bar{h}} A\_{\bar{i}\bar{l}}}{k\_{\bar{i}}(k\_{\bar{i}} - 1)}, & k\_{\bar{i}} \ge 2\\ 0, & k\_{\bar{i}} = 0, 1 \end{cases} \tag{33}$$

Therefore, the global clustering coefficient *C* can be calculated as [50]

$$\mathbf{C} = \frac{1}{N} \sum\_{i}^{N} \mathbf{C}\_{i} \tag{34}$$

The clustering coefficient of a network describes the degree of the collectivity of a complex network. The larger the clustering coefficient of the network, the higher collectivity of the network, which means that the relationship between nodes is closer. The clustering coefficient can also reflect the density of triangles from the definition. The triangles, also known as three-cycles, are related to the stability of the network.

The mean degree change during the simulation and the comparison of degree distribution between the initial and final time steps can be seen in Figure 8. The mean degree of the porous media continues to increase in the simulation but the slope of the curve is constantly decreasing and finally tends to be 0 at the end of the simulation. The larger mean degree means more contact between particles and the porous media tends to reach a steady state. This conclusion can also be revealed from the degree distribution since the number of particles with a larger degree has increased after the simulation.

**Figure 8.** (**a**) Mean degree of porous media over time. (**b**) Comparison of degree distribution between the initial and final time steps.

The mean shortest-path distance and the clustering coefficient during the simulation process are presented in Figure 9. The mean shortest-path distance reduces gradually after a sharp initial decrease. At the end of the simulation, the curve of the simulation tends to be a horizontal line. The mean shortest-path distance indicates the efficiency of information transmission within the network. In the contact force network, it represents the efficiency of force transmission. Shorter mean shortest-path distance means higher force transmission efficiency between particles, which results in a more stable porous media. The tendency of the clustering coefficient is similar to that of the mean degree in the simulation. The clustering coefficient shows the cohesive tendency of the sand particles and the measure of the number of three-cycles in the network. A larger clustering coefficient means more three-cycles in the network, which reflects the network is more stable. All the network parameters indicate the porous media continues to evolve as simulations and gradually becomes stable.

#### *3.2. Absolute Permeability and Relative Permeability*

In this section, we present the analysis of absolute permeability and relative permeability of the porous media. For the absolute permeability, fluid flow was simulated by employing the single-phase Lattice Boltzmann method. First, we extract the structure of porous media at the different time steps of the sand production simulation for absolute permeability calculation. Then the fluid with a lattice density of 1 and zero velocity can be initialized throughout the simulation domain. A constant external body force F = 1 × <sup>10</sup>−<sup>5</sup> is applied along the x-direction to drive the flows. The periodic boundary condition is employed in all directions of the simulation domain, while the boundary condition for the solid particle is the bounce-back boundary condition. The simulation keeps calculating

until reaching a steady state, which is defined by the relative error of permeability between n and n-1 time steps and can be expressed as

**Figure 9.** The network parameters change over time for porous media. (**a**) mean shortest-path distance. (**b**) clustering coefficient.

For the relative permeability, fluid flow was simulated by employing the two-phase Shen–Chen Lattice Boltzmann model. Only the initial and the final structure of the porous media in the sand production simulation was employed in the relative permeability calculation. Different water saturation was assigned in the simulation domain by initializing the nodes with wetting fluid of different fractions. A constant external body force F = 1 × <sup>10</sup>−<sup>5</sup> is assigned to drive the fluid flow along the x-direction. The periodic boundary condition is employed in the x-direction, while the bounce-back boundary condition is adopted for the y and z directions. As with the absolute permeability calculation, the simulation should reach a steady state to obtain reliable results. In relative permeability calculation, the relative error of permeability for both fluids should be less than the convergence limit.

Figure 10a shows the porosity and the absolute permeability as the sand production simulation. We can see both the porosity and permeability increase during this process and finally becomes almost unchanged in the later stage of simulation. This can be explained by the gradually steady state after a period of sand production. The relative permeability of the porous media at the initial and final steps can be seen in Figure 10b. It can be spotted that the connate water saturation and residual oil saturation become smaller at the end of the simulation, which means wider saturation for two-phase flow than that at the initial time step. Besides this, the figure also shows higher relative permeability at the endpoints. This indicates that permeability and connectivity have improved after the sand production process. Moreover, the figure also reveals that the points at which *k*rw = 0 and *k*ro = *k*rw shift to the right side of the relative permeability curve. This can reflect that the degree of preferential wettability for water has increased after the sand production process.

#### *3.3. The Effect of Injection Flow Rate*

Sand production simulations with different injection rates at different confining pressures are conducted in this section. The injection rates employed in the simulations are <sup>2</sup> × <sup>10</sup>−<sup>6</sup> m3/s, 4 × <sup>10</sup>−<sup>6</sup> m3/s, 2 × <sup>10</sup>−<sup>5</sup> <sup>m</sup>3/s, and 4 × <sup>10</sup>−<sup>5</sup> m3/s, respectively. The confining pressures in the simulation are 0.8 MPa and 3 MPa. The cementation strength employed in the simulation is 1 × <sup>10</sup><sup>6</sup> Pa. The cumulative weight of produced sand with different injection rates at different confining pressures can be seen in Figure 11.

**Figure 10.** (**a**) porosity and normalized permeability change over time for porous media. (**b**) comparison of relative permeability for porous media at the initial and final steps.

**Figure 11.** Cumulative weight of produced sand at different injection rates. (**a**) confining pressure equals 0.8 MPa (**b**) confining pressure equals 3 MPa.

At lower confining pressure, we can clearly see that the curves of cumulative produced sand weight with larger flow rate are steeper than that with smaller flow rate from Figure 11a. This is a direct reflection of more produced sand for the rock in larger injection rate conditions. Due to the larger injection rate, the larger hydrodynamic force will exert on the particles and lead to more bond failure. Furthermore, it can be spotted from the figure that the rock under a larger injection rate produces sand earlier than the smaller ones. This reveals that the bond in the rock is more easily broken with the larger injection rate. Although the overall trend of the cumulative weight of produced sand keeps increasing for all simulations, the increase becomes smaller in the later stage for the simulation with a smaller injection flow rate. This indicates that the porous media will gradually tend to reach a steady state. From the figure, we can find that the increment of cumulative produced sand weight between different flow rates becomes larger as the flow rate increases. Therefore, we can conclude that productivity can be improved by allowing a limited amount of produced sand in the oil reservoirs, which corresponds with previous research [53,54].

However, the curves of cumulative produced sand weight reach a steady state after the initial slow increase for smaller flow rate conditions at a higher confining pressure, which can be clearly seen in Figure 11b. For injection rates equal to 2 × <sup>10</sup>−<sup>5</sup> m3/s and <sup>4</sup> × <sup>10</sup>−<sup>5</sup> m3/s, the cumulative weight of produced sand continuously increases as the process goes by for simulations due to their larger injection rate. Compared with curves at the lower confining pressure, their overall slopes are smaller at the higher confining pressure. For injection rates equal to 2 × <sup>10</sup>−<sup>6</sup> m3/s and 4 × <sup>10</sup>−<sup>6</sup> m3/s, the curves of cumulative produced sand weight tend to overlap each other as seen in the figure. This

phenomenon can be explained by the existence of a critical flow rate. For the conditions with a flow rate smaller than the critical value at a certain confining pressure, the total amount of produced sand is a definite value no matter what the flow rate is. It is important to determine the critical flow rate to avoid producing too much sand when the critical flow rate is exceeded. Therefore, many researchers proposed methods for predicting the critical flow rate [55–57]. In our simulation, the critical flow rate can be defined as the flow rate below which a steady state can be reached. From the figure, the weight of produced sand rapidly increases for flow rate larger than 4 × <sup>10</sup>−<sup>6</sup> m3/s, which reveals the flow rate exceeds the critical flow rate at this confining pressure.

#### *3.4. The Effect of Cementation Strength*

We conduct sand production simulations with different cementation strengths at different confining pressures in this section. The cementation strengths are 5 × 105 Pa, <sup>1</sup> × 106 Pa, and 5 × 106 Pa, respectively. The confining pressures in the simulation are 0.8 MPa and 3 MPa. The injection rate employed in the simulation is 2 × <sup>10</sup>−<sup>5</sup> m3/s. The cumulative weight of produced sand with different cementation strengths at different confining pressures can be seen in Figure 12.

**Figure 12.** Cumulative weight of produced sand at different cementation strength (**a**) confining pres-sure equals 0.8 MPa (**b**) confining pressure equals 3 MPa.

For both higher and lower confining pressure, the overall trend of the cumulative weight of produced sand keeps increasing for all simulations. The rock with smaller cementation strength will produce more sand because weaker bonds are more easily broken by the hydrodynamic force at the same injection rate. The curves of cumulative produced sand weight at lower confining pressure are much steeper than that at higher confining pressure. This is due to higher confining pressure which keeps the porous media more stable. Besides this, the increase of cumulative weight of the produced sand tends to be smaller in the later stage of simulation than that at the beginning at higher confining pressure. This indicates that the porous media gradually evolve to steady conditions. In this stage, the rock with larger cementation strength changes less because fewer bonds can be broken due to the hydrodynamic force in this flow rate.

#### *3.5. The Effect of Confining Pressure*

Simulations with different confining pressure are conducted in this section. The confining pressures employed in the simulations are 0.1 MPa, 1 Mpa, and 3 Mpa, respectively. The injection rate employed in the simulations is 4 × <sup>10</sup>−<sup>6</sup> <sup>m</sup>3/s. The cementation strength employed in the simulations is 1 × 106 Pa. The cumulative weight of the produced sand at different confining pressure can be seen in Figure 13.

**Figure 13.** Cumulative weight of produced sand at different confining pressure.

It can be spotted that the cumulative weight of produced sand continuously increases as the process goes by for simulations with confining pressure equal to 1 MPa and 0.1 MPa. For the simulation with injection rate equals to 3 MPa, the cumulative weight of produced sand first slowly increases and then reaches a steady state in the simulation. Moreover, the rock under higher confining pressure will produce less sand. This is due to the better stability caused by the high confining pressure since it will increase the shear strength and the friction resistance between sand particles. Therefore, there will be enough resistance in the simulation with a smaller injection rate to avoid sand production at higher confining pressure conditions.

#### **4. Conclusions**

This paper presents a coupled lattice Boltzmann method and discrete element method (LBM-DEM) to study the sand production process of the formation. The method has been validated by simulating the sedimentation of a single particle. Simulation of the sand production process is conducted and the force chain network evolvement is analyzed. Absolute and relative permeability before and after the sand production process are calculated and the changes in the properties of porous media are studied. The effect of injection flow rate, cementation strength, and confining pressure are investigated and conclusions are listed in the following.

Strong force chain rupture can be recognized and force chain reorganization can be identified. Sand particles will detach from the reservoir and the formation will reach a steady state after a period. The total amount of produced sand first increases and remains constant after the force chain reorganization. Meanwhile, the mean degree and the clustering coefficient of the porous media continue to increase in the simulation but the slope of the curve is constantly decreasing and finally tends to be 0 at the end of the simulation. However, the mean shortest-path distance reduces gradually after a sharply initial decrease. This reveals that the state of the porous media becomes more stable in the sand production process.

By analyzing the change of relative permeability, wider saturation for two-phase flow and higher relative permeability at endpoints can be spotted at the end of the process. This indicates that permeability and connectivity have improved after the sand production process. The degree of preferential wettability for water has increased after the sand production process since the points at which *k*rw = 0 and *k*ro = *k*rw shift to the right side of the relative permeability curve.

The flow rate, cementation strength, and confining pressure have an important effect on the sand production process. A higher injection rate leads to more sand production than a lower injection rate. Besides this, a critical flow rate below which porous media can reach a steady state exists. Moreover, productivity can be improved by allowing a limited amount of produced sand in the oil reservoirs. A reservoir with lower cementation strength produces more sand than the higher one due to weak bonds. Results show that porous media under higher confining pressure will be more stable due to the higher friction resistance between particles to prevent sand production.

**Author Contributions:** Conceptualization, Q.F. and T.X; methodology, S.W.; simulation, T.X. and J.Z.; validation, W.Z.; formal analysis, T.X. and X.Z.; investigation, S.W.; writing—original draft preparation, T.X. and W.Z.; writing—review and editing, J.Z. and S.W.; supervision, Q.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (Grant No. 51904319) and the National Science and Technology Major Project of China (Grant No. 2016ZX05011-001).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare that they have no known competing financial interest or personal relationship that could have appeared to influence the work reported in this paper.

### **Nomenclature**

Symbols


#### **References**

