Comparison of Concentration Transport Approach and MP-PIC Method for Simulating Proppant Transport Process

*Junsheng Zeng and Heng Li*

## **Abstract**

In this work, proppant transport process is studied based on two popular numerical methods: multiphase particle-in-cell method (MP-PIC) and concentration transport method. Derivations of governing equations in these two frameworks are reviewed, and then similarities and differences between these two methods are fully discussed. Several cases are designed to study the particle settling and conveying processes at different fluid Reynolds number. Simulation results indicate that two physical mechanisms become significant in the high Reynolds number cases, which leads to big differences between the simulation results of the two methods. One is the gravity convection effect in the early stage and the other is the particle packing, which determines the shape of sandbank. Above all, the MP-PIC method performs better than the concentration transport approach because more physical mechanisms are considered in the former framework. Besides, assumptions of ignoring unsteady terms and transient terms for the fluid governing equations in the concentration transport approach are only reasonable when Reynolds number is smaller than 100.

**Keywords:** proppant transport, two-phase flow, gravity convection, multiphase particle-in-cell method, concentration transport approach

## **1. Introduction**

In unconventional oil and gas industry, there exists a significant granular flow process, which is known as the proppant transport [1]. It is necessary to pump highstrength granular materials such as ceramic particles and sand into the stimulated fracture networks with carrying fluid. Eventually after the flow-back of fluid, the granular materials remain in the fractures and fracture networks are efficiently propped, which contributes to a high conductivity for gas/oil exploitation. Therefore, it is important to reveal the physical mechanisms in the proppant transport process.

Essentially, proppant transport process is a two-phase flow problem constrained in a channel with various widths. In previous works, concentration transport approach was very popular for simulating the proppant transport. In the approach, proppant is considered as a continuum, and is quantitatively described using

concentration. The motion of the particle phase is solved based on a concentration transport equation. This method is firstly established by Mobbs and Hammond [2]. They derived the governing equations for the fluid-particle mixture (i.e., slurry) by combining mass conservation laws of two phases and convection models. Then a Poisson equation for the fluid pressure is obtained. They proposed an important dimensionless number, that is, the Buoyancy number, to quantitatively describe the relative intensity of horizontal convection effect and the vertical settling effect. Their pioneering work is then adopted and extended in many later works. For example, Gadde and Sharma [3] and Gu and Mohanty [4] extended this framework by considering the effects of fracture propagation. Wang et al. [5] introduced a blocking function in order to consider the proppant bridging effect. Dontsov and Peirce [6] utilized a more accurate velocity retardation model based on their theoretical analysis. Roostaei et al. [7] applied the WENO (weighted essentially nonoscillation) scheme to solve the concentration transport equation, which greatly reduces the numerical diffusion.

The MP-PIC method [8, 9] is another numerical method for simulating largescale fluid-particle coupling system, which is popular in chemical engineering. In the MP-PIC method, fluid motion is governed by the volume-averaged Navier-Stokes equation, and particle motion is solved using the Newton's second law under the Lagrangian framework, which is different from those in the concentration approach. Due to its Lagrangian feature and high fidelity, the MP-PIC method is also shown to be a powerful tool for simulating proppant transport process.

where *w* is the fracture width. Then it is trivial to obtain the differential form [2]:

þ ∇ � ð Þ 1 � *C wu*!

*p*

*f*

*f*

*up*,*<sup>x</sup>* ¼ *u <sup>f</sup>*,*<sup>x</sup>*, *up*,*<sup>y</sup>* ¼ *u <sup>f</sup>*,*<sup>y</sup>* þ *usettling* , (5)

<sup>¼</sup> <sup>0</sup>*:* (2)

<sup>¼</sup> <sup>0</sup>*:* (3)

<sup>¼</sup> *wuleak*, (4)

<sup>¼</sup> *wuleak*, (6)

is the slurry velocity.

þ ∇ � *Cwu*!

*Comparison of Concentration Transport Approach and MP-PIC Method for Simulating…*

Considering fluid leak-off, it is trivial to add a source term in the RHS:

where *uleak* is the leak velocity, which has the dimension of 1/[TIME].

þ ∇ � ð Þ 1 � *C wu*!

In large-scale cases, we assume that in the horizontal direction, the velocities of the two phases are the same (homogeneous slurry flow), and in the vertical direction, the particle phase velocity differs from the fluid phase due to particle settling.

where *usettling* is the particle settling velocity. There are also other modifications indicating that in the horizontal direction particle velocity does not equal the fluid

*s*

*<sup>f</sup>* þ *Cusettling j*

Eq. (6) is necessary for solving fluid pressure and it is illustrated below. As we know, for pure Newtonian fluid, we have the constitutive laws in which viscosity is a significant parameter. In viscous case (low Re number), based on the Poiseuille's Law, we can derive the relationship between the pressure gradient and average

!

Using Eqs. (2) and (4), we can obtain the fluid/particle mixture (slurry)

<sup>þ</sup> <sup>∇</sup> � *wu*!

*<sup>p</sup>* ¼ *u* !

*<sup>∂</sup>*ð Þ *Cw ∂t*

*<sup>∂</sup>*ð Þ ð Þ <sup>1</sup> � *<sup>C</sup> <sup>w</sup> ∂t*

*<sup>∂</sup>*ð Þ ð Þ <sup>1</sup> � *<sup>C</sup> <sup>w</sup> ∂t*

velocity, which is called proppant retardation [3, 4].

*<sup>∂</sup>*ð Þ *<sup>w</sup> ∂t*

*<sup>f</sup>* þ *Cu*!

Then we have the following formula:

*Control volume in concentration transport approach.*

*DOI: http://dx.doi.org/10.5772/intechopen.92227*

*<sup>s</sup>* ¼ ð Þ 1 � *C u*

velocity in a channel/fracture:

!

continuity equation.

**Figure 1.**

where *u* !

**181**

Similarly, we have the mass balance equation for fluid phase:

In this work, the two above numerical methods are both applied in simulating the proppant transport process. Though the two numerical methods are built under different frameworks, there exist both similarities and differences between them. The hidden facts are revealed based on the analysis of the governing equation sets, as well as the numerical results. The remaining contents of this work are organized as follows. In Section 2, basic governing equation sets of the two methods are demonstrated, and relationship between the two methods are discussed. In Section 3, several numerical cases are designed to illustrate the performance of the two methods, and the numerical results are then discussed. Finally, conclusions are drawn in the Section 4.

## **2. Methodology**

#### **2.1 Concentration transport approach**

Assume that there exists only one kind of proppant (same density and size, or mono-disperse) in the fracture and the particle phase is well distributed so that in the large scale we can take the derivative of the particle concentration in most regions (except discontinuity), and the particle and fluid phase are both incompressible. Then we have following unknown variables: fluid velocity *u* ! *<sup>f</sup>* , particle velocity *u* ! *<sup>p</sup>*, particle concentration *C*, and fluid pressure *P*.

Let us start from the continuity equation. **Figure 1** shows the fluxes and accumulation in a control volume. It is clear that for the particle phase we have the mass balance equation:

$$\begin{split} \frac{\mathbf{C}^{n+1} w^{n+1} - \mathbf{C}^{n} w^{n}}{\Delta t} &= \frac{\left(\mathbf{C} w u\_{p, \mathbf{x}}\right)\_{W} \cdot dy - \left(\mathbf{C} w u\_{p, \mathbf{x}}\right)\_{E} \cdot dy}{d\mathbf{x}} \\ &+ \frac{\left(\mathbf{C} w u\_{p, \mathbf{y}}\right)\_{N} \cdot d\mathbf{x} - \left(\mathbf{C} w u\_{p, \mathbf{y}}\right)\_{S} \cdot d\mathbf{x}}{dy} \end{split} \tag{1}$$

*Comparison of Concentration Transport Approach and MP-PIC Method for Simulating… DOI: http://dx.doi.org/10.5772/intechopen.92227*

**Figure 1.** *Control volume in concentration transport approach.*

concentration. The motion of the particle phase is solved based on a concentration transport equation. This method is firstly established by Mobbs and Hammond [2]. They derived the governing equations for the fluid-particle mixture (i.e., slurry) by combining mass conservation laws of two phases and convection models. Then a Poisson equation for the fluid pressure is obtained. They proposed an important dimensionless number, that is, the Buoyancy number, to quantitatively describe the relative intensity of horizontal convection effect and the vertical settling effect. Their pioneering work is then adopted and extended in many later works. For example, Gadde and Sharma [3] and Gu and Mohanty [4] extended this framework by considering the effects of fracture propagation. Wang et al. [5] introduced a blocking function in order to consider the proppant bridging effect. Dontsov and Peirce [6] utilized a more accurate velocity retardation model based on their theoretical analysis. Roostaei et al. [7] applied the WENO (weighted essentially nonoscillation) scheme to solve the concentration transport equation, which greatly

The MP-PIC method [8, 9] is another numerical method for simulating largescale fluid-particle coupling system, which is popular in chemical engineering. In the MP-PIC method, fluid motion is governed by the volume-averaged Navier-Stokes equation, and particle motion is solved using the Newton's second law under the Lagrangian framework, which is different from those in the concentration approach. Due to its Lagrangian feature and high fidelity, the MP-PIC method is also shown to be a powerful tool for simulating proppant transport process.

In this work, the two above numerical methods are both applied in simulating the

Assume that there exists only one kind of proppant (same density and size, or mono-disperse) in the fracture and the particle phase is well distributed so that in the large scale we can take the derivative of the particle concentration in most regions (except discontinuity), and the particle and fluid phase are both incom-

Let us start from the continuity equation. **Figure 1** shows the fluxes and accumulation in a control volume. It is clear that for the particle phase we have the mass

*<sup>W</sup>* � *dy* � *Cwup*,*<sup>x</sup>*

*<sup>N</sup>* � *dx* � *Cwup*,*<sup>y</sup>*

*dy*

*dx*

*<sup>E</sup>* � *dy*

*<sup>S</sup>* � *dx*

*Cwup*,*<sup>y</sup>*  !

*<sup>f</sup>* , particle

, (1)

pressible. Then we have following unknown variables: fluid velocity *u*

*<sup>p</sup>*, particle concentration *C*, and fluid pressure *P*.

*<sup>Δ</sup><sup>t</sup>* <sup>¼</sup> *Cwup*,*<sup>x</sup>*

þ

proppant transport process. Though the two numerical methods are built under different frameworks, there exist both similarities and differences between them. The hidden facts are revealed based on the analysis of the governing equation sets, as well as the numerical results. The remaining contents of this work are organized as follows. In Section 2, basic governing equation sets of the two methods are demonstrated, and relationship between the two methods are discussed. In Section 3, several numerical cases are designed to illustrate the performance of the two methods, and the numerical results are then discussed. Finally, conclusions are drawn in the Section 4.

reduces the numerical diffusion.

*Progress in Fine Particle Plasmas*

**2. Methodology**

velocity *u* !

**180**

balance equation:

*C<sup>n</sup>*þ<sup>1</sup>

*<sup>w</sup><sup>n</sup>*þ<sup>1</sup> � *<sup>C</sup>nwn*

**2.1 Concentration transport approach**

where *w* is the fracture width. Then it is trivial to obtain the differential form [2]:

$$\frac{\partial (\mathcal{C}w)}{\partial t} + \nabla \cdot \left(\mathcal{C}w \vec{u}\_p\right) = \mathbf{0}.\tag{2}$$

Similarly, we have the mass balance equation for fluid phase:

$$\frac{\partial((\mathbf{1} - \mathbf{C})w)}{\partial t} + \nabla \cdot \left((\mathbf{1} - \mathbf{C})w\overrightarrow{u}\_f\right) = \mathbf{0}.\tag{3}$$

Considering fluid leak-off, it is trivial to add a source term in the RHS:

$$\frac{\partial ((\mathbf{1} - \mathbf{C})\boldsymbol{\omega})}{\partial \mathbf{t}} + \nabla \cdot \left( (\mathbf{1} - \mathbf{C})\boldsymbol{\omega}\vec{\boldsymbol{u}}\_f \right) = \boldsymbol{\omega}\boldsymbol{u}\_{label},\tag{4}$$

where *uleak* is the leak velocity, which has the dimension of 1/[TIME].

In large-scale cases, we assume that in the horizontal direction, the velocities of the two phases are the same (homogeneous slurry flow), and in the vertical direction, the particle phase velocity differs from the fluid phase due to particle settling. Then we have the following formula:

$$
\mathfrak{u}\_{p,\mathbf{x}} = \mathfrak{u}\_{\,f,\mathbf{x}}, \; \mathfrak{u}\_{p,\mathbf{y}} = \mathfrak{u}\_{\,f,\mathbf{y}} + \mathfrak{u}\_{\text{setting}}, \tag{5}
$$

where *usettling* is the particle settling velocity. There are also other modifications indicating that in the horizontal direction particle velocity does not equal the fluid velocity, which is called proppant retardation [3, 4].

Using Eqs. (2) and (4), we can obtain the fluid/particle mixture (slurry) continuity equation.

$$\frac{\partial(w)}{\partial t} + \nabla \cdot \left(w\overrightarrow{u}\_s\right) = wu\_{leak},\tag{6}$$

where *u* ! *<sup>s</sup>* ¼ ð Þ 1 � *C u* ! *<sup>f</sup>* þ *Cu*! *<sup>p</sup>* ¼ *u* ! *<sup>f</sup>* þ *Cusettling j* ! is the slurry velocity.

Eq. (6) is necessary for solving fluid pressure and it is illustrated below. As we know, for pure Newtonian fluid, we have the constitutive laws in which viscosity is a significant parameter. In viscous case (low Re number), based on the Poiseuille's Law, we can derive the relationship between the pressure gradient and average velocity in a channel/fracture:

$$
\overrightarrow{\boldsymbol{\mu}}\_f = -\frac{\boldsymbol{w}^2}{12\mu\_f} \nabla \left( \boldsymbol{P} + \rho\_f \mathbf{g} \mathbf{z} \right),
\tag{7}
$$

where *u* !

force in the term 1 � *<sup>ρ</sup> <sup>f</sup>*

hydrostatic pressure.

*∂w <sup>∂</sup><sup>t</sup>* � *<sup>∂</sup> ∂x*

8

>>>>>>>>><

>>>>>>>>>:

*u* !

*∂ α fw* � � *∂t*

1 *w*

8

>>>>>>>>>><

>>>>>>>>>>:

**183**

*du*! *p dt* <sup>¼</sup> *Dp <sup>u</sup>*

*<sup>∂</sup>*ð Þ *Cw ∂t*

*ρp* � �*<sup>g</sup>*

*DOI: http://dx.doi.org/10.5772/intechopen.92227*

**2.3 Comparison between the two methods**

Concentration transport approach:

*w*3 12*μ*ð Þ *C*

þ ∇ � *Cwu*!

!

þ ∇ � *α fwu*!

þ 1 *w*

and fluid-particle interaction force *f particle*.

sets, which are denoted as set I and set II, respectively:

then the second equation in set II is simplified as:

If the term *<sup>f</sup> wall* converges into the formula of <sup>12</sup>*μ α*ð Þ*<sup>p</sup>*

without hydrostatic pressure. Note that here *u*

*f* � �

> ! *<sup>f</sup>* � *u* ! *p* � � <sup>þ</sup> <sup>1</sup> � *<sup>ρ</sup> <sup>f</sup>*

*<sup>s</sup>* ¼ ð Þ 1 � *C u*

MP-PIC approach:

*∂ α fwu*!

*∂t*

� �

*∂P ∂x*

*p* � � <sup>¼</sup> <sup>0</sup>

*<sup>f</sup>* þ *Cu*!

*f* � � <sup>¼</sup> *wuleak*

∇ � *α fwu*!

*<sup>f</sup>* ⊗ *u* ! *f* � � ¼ � <sup>1</sup>

*ρp*

*g* ! � <sup>1</sup> *ρp*

In Eq. (13), we need to introduce models for settling velocity *usettling* and effective viscosity *μ*ð Þ *C* . In Eq. (14), we need to model the two terms: wall friction *f wall*

Next we will discuss the similarity and difference between these two equation

a. Fluid part: it is clear that in set II unsteady and convection/inertial terms are considered. Set I is suitable for homogeneous slurry flow, which indicates that particles settle pretty slowly and Reynolds number is low. In slick water cases, set II is preferred. Actually set II shall converge into set I in low Re number cases. In steady cases, the unsteady terms vanish and the fluid-particle interaction term converges into additional gravity force of particle phase,

∇h i *P* ¼ �*f wall* � *α<sup>p</sup> ρ<sup>p</sup>* � *ρ <sup>f</sup>*

!

*up*,*<sup>x</sup>* ¼ *u <sup>f</sup>*,*<sup>x</sup>*, *up*,*<sup>y</sup>* ¼ *u <sup>f</sup>*,*<sup>y</sup>* þ *usettling*

� *∂ ∂z*

*<sup>p</sup>* ¼ *u* !

*<sup>p</sup>* is the parcel velocity, *Dp* is the drag coefficient, *α<sup>p</sup>* is the volume fraction

*<sup>∂</sup> <sup>P</sup>* <sup>þ</sup> *<sup>ρ</sup><sup>s</sup>* ð Þ *gz ∂z*

!

*ρ f*

∇ � *PI* � *T <sup>f</sup>*

� �*<sup>g</sup>*

! *s* D E <sup>¼</sup> *<sup>α</sup> <sup>f</sup> <sup>u</sup>*

*<sup>ρ</sup> <sup>f</sup> <sup>w</sup>*<sup>2</sup> *u* !

> ! *f* D E <sup>þ</sup> <sup>1</sup> *V <sup>f</sup>*

¼ � *<sup>w</sup>*<sup>2</sup>

� �

! and the fluid pressure is the net pressure excluding the

¼ *wuleak*

<sup>12</sup>*μ*ð Þ *<sup>C</sup>* <sup>∇</sup> *<sup>P</sup>* <sup>þ</sup> *<sup>ρ</sup><sup>s</sup>* ð Þ *gz*

∇*P* � *f wall* þ *f particle*

*θpρ<sup>p</sup>*

! (15)

*<sup>s</sup>* we can recover the Eq. (8)

P*Vp u* ! *p*

<sup>∇</sup>*τ<sup>p</sup>* � <sup>1</sup> *τw u* ! *p*

� � � <sup>1</sup>

(13)

*:*

(14)

of the particle phase, *τ<sup>p</sup>* is the particle normal stress, and *τ<sup>w</sup>* is the damping relaxation time due to particle-wall interaction, *P* is the local average fluid pressure, and *τ <sup>f</sup>* is the local average fluid viscous stress tensor. Here we consider the Archimedes buoyancy

*Comparison of Concentration Transport Approach and MP-PIC Method for Simulating…*

First let us summarize the governing equations of the two approaches.

*w*3 12*μ*ð Þ *C*

*<sup>f</sup>* þ *Cusettling j*

where *μ <sup>f</sup>* is the fluid viscosity and *w* is the fracture width. The term *w*2*=*12 is also considered as the effective permeability of a fracture. This relationship connects the fluid pressure with the velocity, and if we substitute it into continuity equation, we will obtain a Poisson equation for pressure and it can be easily solved.

In the case of fluid/particle mixture, we also expect there exists a similar relationship between fluid pressure and slurry velocity. There are many previous literatures including experimental and numerical works revealing this relationship. It is well known that the apparent viscosity of the mixture is higher than that of the pure fluid and a formula similar to Eq. (7) can be obtained introducing the effective viscosity of the mixture [7]:

$$\overrightarrow{\boldsymbol{\mu}}\_{s} = -\frac{\boldsymbol{w}^{2}}{12\mu(\mathbf{C})}\nabla(\boldsymbol{P} + \rho\_{s}\mathbf{g}\mathbf{z}),\tag{8}$$

where *μ*ð Þ *C* is the effective viscosity, which depends on particle concentration *C* and obviously *μ*ð Þ¼ 0 *μ <sup>f</sup>* , *ρ<sup>s</sup>* ¼ *ρ <sup>f</sup>*ð Þþ 1 � *C ρpC* is the slurry density.

By substituting Eq. (8) into Eq. (6), we get the following pressure Eq. (7):

$$\frac{\partial w}{\partial t} - \frac{\partial}{\partial \mathbf{x}} \left( \frac{w^3}{12\mu(\mathbf{C})} \frac{\partial P}{\partial \mathbf{x}} \right) - \frac{\partial}{\partial \mathbf{z}} \left( \frac{w^3}{12\mu(\mathbf{C})} \frac{\partial (P + \rho\_\mathbf{g} \mathbf{z})}{\partial \mathbf{z}} \right) = w u\_{leak} \tag{9}$$

If the fracture width does not vary with time, then the time derivative vanishes, with the pressure Poisson equation remaining, and the sand concentration *C* can be solved explicitly or implicitly using Eq. (2) with high-order accuracy schemes.

#### **2.2 MP-PIC method**

The MP-PIC method is an Eulerian-Lagrangian method, in which fluid motion is solved in the Eulerian grids and particle motion is solved under the Lagrangian framework. The governing equation of fluid motion reads [10]:

$$\frac{\partial (a\_f w)}{\partial t} + \nabla \cdot \left(a\_f w \overrightarrow{u}\_f\right) = \omega u\_{leak},\tag{10}$$

$$\frac{1}{w}\frac{\partial \left(a\_f w \overrightarrow{u}\_f\right)}{\partial t} + \frac{1}{w} \nabla \cdot \left(a\_f w \overrightarrow{u}\_f \otimes \overrightarrow{u}\_f\right) = -\frac{1}{\rho\_f} \nabla P + \overrightarrow{f}\_w + \overrightarrow{f}\_p,\tag{11}$$

where *w* is the fracture width, *α <sup>f</sup>* is the volume fraction of fluid, *f* ! *<sup>w</sup>* is the wall friction term and *f* ! *<sup>p</sup>* is the fluid-particle coupling source term, *u* ! *<sup>f</sup>* is the fluid velocity, symbol " ⊗ " indicates the union product.

In the MP-PIC method, the particle phase is discretized into parcels, and every parcel represents several physical particles, which possess the same properties such as density and size, and also similar kinetic behavior such as the velocity and acceleration. Parcel motion is governed by the Newton's second law listed as follows:

$$\frac{d\overrightarrow{u}\_p}{dt} = D\_p \left(\overrightarrow{u}\_f - \overrightarrow{u}\_p\right) + \left(\mathbf{1} - \frac{\rho\_f}{\rho\_p}\right)\overrightarrow{\mathbf{g}} - \frac{\mathbf{1}}{\rho\_p}\nabla \cdot \left(\vec{P}\overleftarrow{\overline{I}} - \overrightarrow{\overline{\tau}}\_f\right) - \frac{\nabla \tau\_p}{a\_p \rho\_p} - \frac{\mathbf{1}}{\tau\_w}\overrightarrow{u}\_p \tag{12}$$

*Comparison of Concentration Transport Approach and MP-PIC Method for Simulating… DOI: http://dx.doi.org/10.5772/intechopen.92227*

where *u* ! *<sup>p</sup>* is the parcel velocity, *Dp* is the drag coefficient, *α<sup>p</sup>* is the volume fraction of the particle phase, *τ<sup>p</sup>* is the particle normal stress, and *τ<sup>w</sup>* is the damping relaxation time due to particle-wall interaction, *P* is the local average fluid pressure, and *τ <sup>f</sup>* is the local average fluid viscous stress tensor. Here we consider the Archimedes buoyancy force in the term 1 � *<sup>ρ</sup> <sup>f</sup> ρp* � �*<sup>g</sup>* ! and the fluid pressure is the net pressure excluding the hydrostatic pressure.

### **2.3 Comparison between the two methods**

First let us summarize the governing equations of the two approaches. Concentration transport approach:

$$\begin{cases} \frac{\partial w}{\partial t} - \frac{\partial}{\partial \mathbf{x}} \left( \frac{w^3}{12\mu(\mathbf{C})} \frac{\partial P}{\partial \mathbf{x}} \right) - \frac{\partial}{\partial \mathbf{z}} \left( \frac{w^3}{12\mu(\mathbf{C})} \frac{\partial (P + \rho\_t \mathbf{g} \mathbf{z})}{\partial \mathbf{z}} \right) = w u\_{lank} \\\\ \frac{\partial (\mathbf{C}w)}{\partial t} + \nabla \cdot \left( \mathbf{C}w \overrightarrow{u}\_p \right) = \mathbf{0} \\ u\_{p, \mathbf{x}} = u\_{f, \mathbf{x}}, \ u\_{p, \mathbf{y}} = u\_{f, \mathbf{y}} + u\_{\text{settling}} \\\\ \overrightarrow{u}\_t = (\mathbf{1} - \mathbf{C})\overrightarrow{u}\_f + \mathbf{C}\overrightarrow{u}\_p = \overrightarrow{u}\_f + \mathbf{C}u\_{\text{settling}} \overrightarrow{j} = -\frac{w^2}{12\mu(\mathbf{C})} \nabla (P + \rho\_\mathbf{g} \mathbf{z}) \end{cases} (\text{13})$$

MP-PIC approach:

*u* !

*u* !

*w*3 12*μ*ð Þ *C*

� �

viscosity of the mixture [7]:

*Progress in Fine Particle Plasmas*

*∂w <sup>∂</sup><sup>t</sup>* � *<sup>∂</sup> ∂x*

**2.2 MP-PIC method**

1 *w*

friction term and *f*

*du* ! *p dt* <sup>¼</sup> *Dp <sup>u</sup>*

**182**

*∂ α fwu*!

*∂t*

! *<sup>f</sup>* � *u* ! *p* � �

*f* � �

!

þ 1 *w*

velocity, symbol " ⊗ " indicates the union product.

*<sup>f</sup>* ¼ � *<sup>w</sup>*<sup>2</sup> 12*μ <sup>f</sup>*

will obtain a Poisson equation for pressure and it can be easily solved.

*<sup>s</sup>* ¼ � *<sup>w</sup>*<sup>2</sup>

and obviously *μ*ð Þ¼ 0 *μ <sup>f</sup>* , *ρ<sup>s</sup>* ¼ *ρ <sup>f</sup>*ð Þþ 1 � *C ρpC* is the slurry density.

*∂P ∂x*

framework. The governing equation of fluid motion reads [10]:

<sup>∇</sup> � *<sup>α</sup> fwu*!

<sup>þ</sup> <sup>1</sup> � *<sup>ρ</sup> <sup>f</sup>*

!

*ρp*

*g* ! � <sup>1</sup> *ρp*

*∂ α fw* � � *∂t*

∇ *P* þ *ρ <sup>f</sup> gz* � �

where *μ <sup>f</sup>* is the fluid viscosity and *w* is the fracture width. The term *w*2*=*12 is also considered as the effective permeability of a fracture. This relationship connects the fluid pressure with the velocity, and if we substitute it into continuity equation, we

In the case of fluid/particle mixture, we also expect there exists a similar relationship between fluid pressure and slurry velocity. There are many previous literatures including experimental and numerical works revealing this relationship. It is well known that the apparent viscosity of the mixture is higher than that of the pure fluid and a formula similar to Eq. (7) can be obtained introducing the effective

where *μ*ð Þ *C* is the effective viscosity, which depends on particle concentration *C*

*w*3 12*μ*ð Þ *C*

If the fracture width does not vary with time, then the time derivative vanishes, with the pressure Poisson equation remaining, and the sand concentration *C* can be solved explicitly or implicitly using Eq. (2) with high-order accuracy schemes.

The MP-PIC method is an Eulerian-Lagrangian method, in which fluid motion is

*f* � �

> ¼ � <sup>1</sup> *ρ f*

∇ � *PI* � *τ <sup>f</sup>* � �

∇*P* þ *f* ! *<sup>w</sup>* þ *f* !

solved in the Eulerian grids and particle motion is solved under the Lagrangian

<sup>þ</sup> <sup>∇</sup> � *<sup>α</sup> fwu*!

*<sup>f</sup>* ⊗ *u* ! *f*

*<sup>p</sup>* is the fluid-particle coupling source term, *u*

In the MP-PIC method, the particle phase is discretized into parcels, and every parcel represents several physical particles, which possess the same properties such as density and size, and also similar kinetic behavior such as the velocity and acceleration. Parcel motion is governed by the Newton's second law listed as follows:

� �

where *w* is the fracture width, *α <sup>f</sup>* is the volume fraction of fluid, *f*

By substituting Eq. (8) into Eq. (6), we get the following pressure Eq. (7):

¼ *wuleak* (9) .

� *∂ ∂z* , (7)

<sup>12</sup>*μ*ð Þ *<sup>C</sup>* <sup>∇</sup> *<sup>P</sup>* <sup>þ</sup> *<sup>ρ</sup><sup>s</sup>* ð Þ *gz* , (8)

¼ *wuleak*, (10)

!

� <sup>∇</sup>*τ<sup>p</sup> αpρ<sup>p</sup>* � 1 *τw u* !

*<sup>p</sup>*, (11)

*<sup>w</sup>* is the wall

*<sup>p</sup>* (12)

!

*<sup>f</sup>* is the fluid

*<sup>∂</sup> <sup>P</sup>* <sup>þ</sup> *<sup>ρ</sup><sup>s</sup>* ð Þ *gz ∂z*

� �

$$\begin{cases} \frac{\partial \left(a\_f w\right)}{\partial t} + \nabla \cdot \left(a\_f w \overrightarrow{u}\_f\right) = w u\_{leak} \\ \frac{1}{w} \frac{\partial \left(a\_f w \overrightarrow{u}\_f\right)}{\partial t} + \frac{1}{w} \nabla \cdot \left(a\_f w \overrightarrow{u}\_f \otimes \overrightarrow{u}\_f\right) = -\frac{1}{\rho\_f} \nabla P - f\_{wall} + f\_{particle} \\ \frac{d \overrightarrow{u}\_p}{dt} = D\_p \left(\overrightarrow{u}\_f - \overrightarrow{u}\_p\right) + \left(1 - \frac{\rho\_f}{\rho\_p}\right) \overrightarrow{g} - \frac{1}{\rho\_p} \nabla \cdot \left(P \overleftarrow{\overline{I}} - \overrightarrow{\overline{T}\_f}\right) - \frac{1}{\theta\_p \rho\_p} \nabla \tau\_p - \frac{1}{\tau\_w} \overrightarrow{u}\_p \end{cases} . \tag{14}$$

In Eq. (13), we need to introduce models for settling velocity *usettling* and effective viscosity *μ*ð Þ *C* . In Eq. (14), we need to model the two terms: wall friction *f wall* and fluid-particle interaction force *f particle*.

Next we will discuss the similarity and difference between these two equation sets, which are denoted as set I and set II, respectively:

a. Fluid part: it is clear that in set II unsteady and convection/inertial terms are considered. Set I is suitable for homogeneous slurry flow, which indicates that particles settle pretty slowly and Reynolds number is low. In slick water cases, set II is preferred. Actually set II shall converge into set I in low Re number cases. In steady cases, the unsteady terms vanish and the fluid-particle interaction term converges into additional gravity force of particle phase, then the second equation in set II is simplified as:

$$\nabla \langle P \rangle = -f\_{wall} - a\_p \left(\rho\_p - \rho\_f\right) \overrightarrow{\mathbf{g}} \tag{15}$$

If the term *<sup>f</sup> wall* converges into the formula of <sup>12</sup>*μ α*ð Þ*<sup>p</sup> <sup>ρ</sup> <sup>f</sup> <sup>w</sup>*<sup>2</sup> *u* ! *<sup>s</sup>* we can recover the Eq. (8) without hydrostatic pressure. Note that here *u* ! *s* D E <sup>¼</sup> *<sup>α</sup> <sup>f</sup> <sup>u</sup>* ! *f* D E <sup>þ</sup> <sup>1</sup> *V <sup>f</sup>* P*Vp u* ! *p*

according to definition or we can just consider it as the average fluid velocity for simplicity. However, how the wall friction term changes in high Re case is still unrevealed.

where Re *<sup>p</sup>* <sup>¼</sup> <sup>2</sup>*<sup>ρ</sup> <sup>f</sup> <sup>u</sup>*

**3. Case study**

respectively.

**Figure 2.**

**185**

*A sketch for the numerical simulation domain.*

**3.1 Parameter setting**

radius of proppant particle.

! *<sup>f</sup>* �*u* ! j j*<sup>p</sup> rp μ f*

*DOI: http://dx.doi.org/10.5772/intechopen.92227*

these parcels are modified following a bounce-back way.

the low Reynolds number case as case III respectively.

**3.2 Simulation results and discussions**

is the particle Reynolds number, and *rp* is the physical

Numerical simulation is performed in a rectangle domain as shown in **Figure 2**. The left side of the domain is set as inlet. Proppant is uniformly injected from the whole inlet and proppant concentration is set as 20%. Here the proppant concentration means the volume ratio of particles to total volume of fluid/particle mixture. The right side of the domain is set as outlet. The upper and bottom sides of the rectangle domain are set as the non-flow boundary. Particle deposition mechanisms are different in these two methods. In the concentration transport approach, if proppant concentration reaches the close packing limit, particle settling velocity is set to be approaching zero, and it is easy to implement the non-flow condition for a Eulerian approach. In the MP-PIC method, if proppant concentration reaches the close packing limit, additional forces due to particle stress gradient are exerted on parcels to make sure these parcels shall move away from the high-concentration region. When parcels move across the non-flow boundary, the displacements of

*Comparison of Concentration Transport Approach and MP-PIC Method for Simulating…*

Parameters used in the simulation are listed in **Table 1**. Cases with different fluid viscosities, that is, 1, 10, and 100 cP are considered in this section. Because characteristic length is 0.005 m (fracture width) and characteristic velocity is 0.2 m/s (inlet velocity), the Reynolds numbers in these cases are 1000, 100, and 10

**Figures 3**–**5** show the numerical results of three different Reynolds number cases of two methods. The results are contoured by the volume fraction of particle (denoted as "vfp"), or the proppant concentration. Here we denote the high Reynolds number case as case I, the moderate Reynolds number case as case II, and

In case I the viscosity of carrying fluid is very low, that is, 1 cp, which indicates a poor capability of proppant transport. From **Figure 3**, it is clear that both of the two numerical methods illustrate the packing process during the transport process. In the figure, blue, white, and red regions indicate the pure fluid, suspending slurry, and sandbank regions respectively. In case II, settling behavior of proppant is weaker than that of case I, and instead gravity convection of the slurry front in the

b. Particle part: obviously we solve particle phase motion under Eulerian framework in set I while under Lagrangian framework in set II. Density/size distribution is easier to be considered in set II. In set I we directly assign a settling velocity for the particle phase while in set II we can resolve the settling history with the drag relaxation term if the time step is set to be small enough. The terminal velocity in set II is determined by which drag force model is utilized. Both of the settling velocity model in set I and drag force model in set II are modifications of Stokes settling theoretical results. Also unsteady terms including fluid pressure and viscous stress effects are considered in set II.

In both of the two 2D large-scale equation sets, models are necessary for closure issues and they cannot exactly describe the full-scale fluid/particle behavior. "Large scale" has two meanings: large time scale and large spatial scale. Different physical mechanisms play significant roles in different scales. For the time scales, p-p collision occurs in the time scale of "μs," f-p drag occurs in the time scale of "ms," fluid leak-off and fracture width change have significant effects in the time scale of "s" or "min." For the spatial scales, fluid viscous force dominates in the Kolmogorov length scale (mm), gravity convection dominates in the length scale of "m." We expect that in our 2D large-scale models, we can capture the large-scale fluid/particle behavior. However, small-scale fluid/particle interactions shall affect the large-scale phenomena.

Using 3D DNS (direct numerical simulation) we can obtain the full-scale details; however, due to computational limits simulation time at most reaches several minutes and simulation length at most reaches 1 dm. Though it is not possible to perform 3D DNS even for an experimental-scale case, useful models can be extracted from the DNS data, for example effective viscosity, settling velocity, retardation factor etc.

In this work, Barree and Conway's [11] effective viscosity model is utilized to calculate slurry viscosity for the concentration transport approach and calculate wall friction force for the MP-PIC method, which is expressed as follows:

$$
\mu\_f \left( a\_p \right) = \frac{\mu\_0}{\left( 1 - \frac{a}{a\_p} \right)^{1.82}},
\tag{16}
$$

where *α<sup>p</sup>* is the particle volume fraction, or proppant concentration, *αcp* is the close packing particle volume fraction, which is set as 0.6 in this work.

Besides, Wen and Yu's drag force model [12] is utilized for calculating the fluidparticle coupling drag force for the MP-PIC method and determining particle settling velocity for the concentration transport approach, which is expressed as follows:

$$D\_p = \frac{\mathfrak{Z}\rho\_f}{\mathfrak{S}\rho\_p} \mathbf{C}\_d \frac{\left| \overrightarrow{u}\_f - \overrightarrow{u}\_p \right|}{r\_p},\tag{17}$$

$$\mathbf{C}\_{d} = \begin{cases} \frac{24}{\text{Re}\_{p}} (\mathbf{1} + \mathbf{0.15} \,\text{Re}\_{p} \,^{0.687}) a\_{f}^{-2.65} & \text{Re}\_{p} < \mathbf{1}, \text{000} \\\\ 0.44 a\_{f}^{-2.65} & \text{Re}\_{p} \ge \mathbf{1}, \text{000} \end{cases},\tag{18}$$

*Comparison of Concentration Transport Approach and MP-PIC Method for Simulating… DOI: http://dx.doi.org/10.5772/intechopen.92227*

where Re *<sup>p</sup>* <sup>¼</sup> <sup>2</sup>*<sup>ρ</sup> <sup>f</sup> <sup>u</sup>* ! *<sup>f</sup>* �*u* ! j j*<sup>p</sup> rp μ f* is the particle Reynolds number, and *rp* is the physical radius of proppant particle.

## **3. Case study**

according to definition or we can just consider it as the average fluid velocity for simplicity. However, how the wall friction term changes in high Re case is still

framework in set I while under Lagrangian framework in set II. Density/size distribution is easier to be considered in set II. In set I we directly assign a settling velocity for the particle phase while in set II we can resolve the settling history with the drag relaxation term if the time step is set to be small enough. The terminal velocity in set II is determined by which drag force model is utilized. Both of the settling velocity model in set I and drag force model in set II are modifications of Stokes settling theoretical results. Also unsteady terms including fluid pressure and viscous stress effects are

In both of the two 2D large-scale equation sets, models are necessary for closure issues and they cannot exactly describe the full-scale fluid/particle behavior. "Large scale" has two meanings: large time scale and large spatial scale. Different physical mechanisms play significant roles in different scales. For the time scales, p-p collision occurs in the time scale of "μs," f-p drag occurs in the time scale of "ms," fluid leak-off and fracture width change have significant effects in the time scale of "s" or "min." For the spatial scales, fluid viscous force dominates in the Kolmogorov length scale (mm), gravity convection dominates in the length scale of "m." We expect that in our 2D large-scale models, we can capture the large-scale fluid/particle behavior. However, small-scale fluid/particle interactions shall affect the large-scale phenomena. Using 3D DNS (direct numerical simulation) we can obtain the full-scale details;

however, due to computational limits simulation time at most reaches several minutes and simulation length at most reaches 1 dm. Though it is not possible to perform 3D DNS even for an experimental-scale case, useful models can be extracted from the DNS data, for example effective viscosity, settling velocity,

*μ <sup>f</sup> α<sup>p</sup>*

close packing particle volume fraction, which is set as 0.6 in this work.

*Dp* <sup>¼</sup> <sup>3</sup>*<sup>ρ</sup> <sup>f</sup>* 8*ρ<sup>p</sup> Cd u* ! *<sup>f</sup>* � *u* ! *p*

1 þ 0*:*15 Re *<sup>p</sup>*

0*:*44*α*�2*:*<sup>65</sup>

<sup>0</sup>*:*<sup>687</sup> � �*α*�2*:*<sup>65</sup>

particle coupling drag force for the MP-PIC method and determining particle settling velocity for the concentration transport approach, which is expressed as

In this work, Barree and Conway's [11] effective viscosity model is utilized to calculate slurry viscosity for the concentration transport approach and calculate wall friction force for the MP-PIC method, which is expressed as follows:

� � <sup>¼</sup> *<sup>μ</sup>*<sup>0</sup>

1 � *<sup>α</sup> αcp*

where *α<sup>p</sup>* is the particle volume fraction, or proppant concentration, *αcp* is the

Besides, Wen and Yu's drag force model [12] is utilized for calculating the fluid-

� � �

*rp*

� � �

*<sup>f</sup>* Re *<sup>p</sup>* ≥ 1, 000

*<sup>f</sup>* Re *<sup>p</sup>* < 1, 000

� �<sup>1</sup>*:*<sup>82</sup> , (16)

, (17)

, (18)

b. Particle part: obviously we solve particle phase motion under Eulerian

unrevealed.

*Progress in Fine Particle Plasmas*

considered in set II.

retardation factor etc.

follows:

**184**

*Cd* ¼

24 Re *<sup>p</sup>*

8 ><

>:

#### **3.1 Parameter setting**

Numerical simulation is performed in a rectangle domain as shown in **Figure 2**. The left side of the domain is set as inlet. Proppant is uniformly injected from the whole inlet and proppant concentration is set as 20%. Here the proppant concentration means the volume ratio of particles to total volume of fluid/particle mixture. The right side of the domain is set as outlet. The upper and bottom sides of the rectangle domain are set as the non-flow boundary. Particle deposition mechanisms are different in these two methods. In the concentration transport approach, if proppant concentration reaches the close packing limit, particle settling velocity is set to be approaching zero, and it is easy to implement the non-flow condition for a Eulerian approach. In the MP-PIC method, if proppant concentration reaches the close packing limit, additional forces due to particle stress gradient are exerted on parcels to make sure these parcels shall move away from the high-concentration region. When parcels move across the non-flow boundary, the displacements of these parcels are modified following a bounce-back way.

Parameters used in the simulation are listed in **Table 1**. Cases with different fluid viscosities, that is, 1, 10, and 100 cP are considered in this section. Because characteristic length is 0.005 m (fracture width) and characteristic velocity is 0.2 m/s (inlet velocity), the Reynolds numbers in these cases are 1000, 100, and 10 respectively.

### **3.2 Simulation results and discussions**

**Figures 3**–**5** show the numerical results of three different Reynolds number cases of two methods. The results are contoured by the volume fraction of particle (denoted as "vfp"), or the proppant concentration. Here we denote the high Reynolds number case as case I, the moderate Reynolds number case as case II, and the low Reynolds number case as case III respectively.

In case I the viscosity of carrying fluid is very low, that is, 1 cp, which indicates a poor capability of proppant transport. From **Figure 3**, it is clear that both of the two numerical methods illustrate the packing process during the transport process. In the figure, blue, white, and red regions indicate the pure fluid, suspending slurry, and sandbank regions respectively. In case II, settling behavior of proppant is weaker than that of case I, and instead gravity convection of the slurry front in the

**Figure 2.** *A sketch for the numerical simulation domain.*


#### **Table 1.**

*Parameters for simulating proppant transport.*

**Figure 4.**

*respectively.*

**Figure 5.**

**187**

*Simulation results of moderate Reynolds number case (Re = 100). (a)–(d) are the concentration transport results at time = 1, 2, 3, and 4 s respectively, and (e)–(h) are the MP-PIC results at time = 1, 2, 3, and 4 s*

*Comparison of Concentration Transport Approach and MP-PIC Method for Simulating…*

*DOI: http://dx.doi.org/10.5772/intechopen.92227*

*Simulation results of low Reynolds number case (Re = 10). (a)–(d) are the concentration transport results at time = 1, 2, 3, and 4 s respectively, and (e)–(h) are the MP-PIC results at time = 1, 2, 3, and 4 s respectively.*

#### **Figure 3.**

*Simulation results of high Reynolds number case (Re = 1000). (a)–(d) are the concentration transport results at time = 1, 2, 3, and 4 s respectively, and (e)–(h) are the MP-PIC results at time = 1, 2, 3, and 4 s respectively.*

vertical direction is obvious as shown in **Figure 4**. In case III, it is obvious from **Figure 5** that proppant settling behavior is hard to recognize and gravity convection is much weaker than that of case II.

By comparing the results at different Reynolds numbers, it can be summarized that as the slurry is injected into the domain, there are several significant mechanisms that determine the eventual proppant distribution listed as follows.

The first mechanism is proppant settling. Proppant settling is due to the net gravity force of the proppant if particle density is larger than the fluid density, and the terminal velocity of proppant is determined by the particle Reynolds number and particle volume fraction. According to the Wen and Yu's drag force model, it is trivial to obtain the settling velocities of proppant in the above three tests, and they are 53.2, 11.8, and 1.28 mm/s, respectively. Therefore, within a same horizontal transport distance, the level of slurry-pure water interface declines more dramatically in case I compared to the other two cases with lower Reynolds numbers.

The second one is gravity convection. In case II and case III, proppant settling can be ignored compared to the inlet velocity, however the slurry fronts in these

*Comparison of Concentration Transport Approach and MP-PIC Method for Simulating… DOI: http://dx.doi.org/10.5772/intechopen.92227*

#### **Figure 4.**

*Simulation results of moderate Reynolds number case (Re = 100). (a)–(d) are the concentration transport results at time = 1, 2, 3, and 4 s respectively, and (e)–(h) are the MP-PIC results at time = 1, 2, 3, and 4 s respectively.*

#### **Figure 5.**

*Simulation results of low Reynolds number case (Re = 10). (a)–(d) are the concentration transport results at time = 1, 2, 3, and 4 s respectively, and (e)–(h) are the MP-PIC results at time = 1, 2, 3, and 4 s respectively.*

vertical direction is obvious as shown in **Figure 4**. In case III, it is obvious from **Figure 5** that proppant settling behavior is hard to recognize and gravity convection

*Simulation results of high Reynolds number case (Re = 1000). (a)–(d) are the concentration transport results at time = 1, 2, 3, and 4 s respectively, and (e)–(h) are the MP-PIC results at time = 1, 2, 3, and 4 s respectively.*

**Parameter Value Parameter Value** Fluid density 1000 kg/m<sup>3</sup> Particle density 2500 kg/m<sup>3</sup> Fluid viscosity 1 and 10 and 100 cP Particle diameter 0.6 mm Time-step <sup>5</sup> <sup>10</sup><sup>3</sup> <sup>s</sup> Simulation time 4 s Inlet velocity 0.2 m/s Inlet concentration 20% Domain size <sup>1</sup> <sup>0</sup>*:*<sup>25</sup> <sup>0</sup>*:*005 m3 Mesh size <sup>100</sup> <sup>25</sup>

By comparing the results at different Reynolds numbers, it can be summarized that as the slurry is injected into the domain, there are several significant mechanisms that determine the eventual proppant distribution listed as follows.

The first mechanism is proppant settling. Proppant settling is due to the net gravity force of the proppant if particle density is larger than the fluid density, and the terminal velocity of proppant is determined by the particle Reynolds number and particle volume fraction. According to the Wen and Yu's drag force model, it is trivial to obtain the settling velocities of proppant in the above three tests, and they are 53.2, 11.8, and 1.28 mm/s, respectively. Therefore, within a same horizontal transport distance, the level of slurry-pure water interface declines more dramatically in case I compared to the other two cases with lower Reynolds numbers. The second one is gravity convection. In case II and case III, proppant settling can be ignored compared to the inlet velocity, however the slurry fronts in these

is much weaker than that of case II.

**Table 1.**

**Figure 3.**

**186**

*Parameters for simulating proppant transport.*

*Progress in Fine Particle Plasmas*

two cases still evolve and both trend from vertical direction to inclined direction. This is because of the horizontal pressure gradient on the slurry front due to the difference between slurry and pure fluid, which provides a driving force and keeps slurry on the bottom penetrating into the pure fluid region. This mechanism is then intuitively described as gravity convection. Gravity convection in case III is much weaker than that in case II. Though the horizontal pressure gradients on the slurry front are the same in these two cases, fluid viscosity is greater in case III, which leads to a smaller channel permeability for slurry flows. Therefore, the penetration velocity in case III is much smaller than that in case II. In case I, gravity convection also plays a significant role. According to the ratio of inlet velocity and proppant settling velocity, that is, about 4, without gravity convection slurry front is expected to be exact the diagonal line of the domain. However, numerical results of both methods show that the suspending region is far below the diagonal line, which indicates that fluid velocity field is greatly changed due to gravity convection compared to a uniform flow. Above all, gravity convection can be considered as a global effect reflecting the density difference between slurry and pure fluid, whereas proppant settling represents a local effect reflecting the density difference between proppant particle and pure fluid.

results of two methods can exactly reflect the effects of different modeling strategy on proppant transport. Above all, both methods can precisely capture the proppant settling mechanism when the same drag force model is adopted in both methods, and concentration transport approach can capture the gravity convection mechanisms precisely when Reynolds number is smaller than 100. However, proppant packing mechanism is not captured very well in the concentration approach.

*Comparison of Concentration Transport Approach and MP-PIC Method for Simulating…*

In this work, two numerical methods are adopted to simulate proppant transport: the concentration transport approach and the MP-PIC method. The first one is

1.From the view of pure horizontal proppant advection, numerical diffusion can be insignificant in the concentration transport approach if high-order scheme like WENO scheme is adopted, and the interfaces between slurry front and

2.When fluid Reynolds number is smaller than 100, assumptions of ignoring velocity convection term and transient term for fluid governing equation adopted in the concentration transport approach are reasonable, and numerical results of both methods show similar transport patterns.

3.When fluid Reynolds number reaches 1000, that is, in the low viscous cases, numerical experiments prove that the concentration transport approach shows

a low fidelity for capturing both the slurry suspending region and the

4.Generally, the more physical mechanisms, including particle-particle/wall interaction and fluid-particle interaction, accurately considered in a numerical framework, the better a simulator performs on capturing proppant transport

This work is funded by the National Natural Science Foundation of China (nos. U1730111, 11972088, 91852207), the State Key Laboratory of Explosive Science and Technology (no. QNKT19-05), and the National Science and Technology Major

a typical Eulerian-Eulerian method and the second one is a typical Eulerian-Lagrangian method. With full discussions on their frameworks and comparisons between the numerical results, the following conclusions are then obtained:

pure fluid can be clearly captured.

*DOI: http://dx.doi.org/10.5772/intechopen.92227*

sandbank packing process.

Project of China (Grant no. 2017ZX05039-005).

behaviors.

**Acknowledgements**

**189**

**4. Conclusion**

The third one is proppant packing. The two prior mechanisms affect the distribution of slurry suspending region, and proppant packing shall affect the distribution of sandbank. In case I, as time evolves proppant concentration increases at the bottom of the simulation domain. When proppant concentration reaches a maximum value, that is, close to the packing state, particle-particle interactions and particle-wall interactions become more frequent. Then, the early formed sandbank stays unmovable, and the fluid velocity also dramatically decreases in this region due to the fluid-particle coupling and particle-particle/wall damping effects.

By comparing the numerical results of two different methods based on the above mechanisms, similarities and differences between the two methods discussed in Section 2.3 are verified.

Firstly, in case III it is obvious that numerical results of methods are almost the same with each other. In case III, Reynolds number is pretty low and the fluid motion is dominated by the viscous mechanism. The slurry flow is approximately plug flow. Secondly, in case II of moderate Reynolds number, numerical results of concentration transport approach show that the gravity convection is a bit severe than those of MP-PIC method. However, transport patterns of these two methods are in general similar. Quantitatively, the relative length of top and bottom slurry fronts at the end of simulation time is 0.5 and 0.4 m, respectively. Thirdly, in case I, numerical results of two methods differ a lot from each other. For the slurry suspending region, the covering area of the MP-PIC results is obviously much larger than those of concentration transport approach. This is mainly because the transient term and convection term in the fluid governing equation are ignored in the concentration approach, and these two mechanisms play significant roles in low viscous or high Reynolds number cases. It is clear that losses of these two mechanisms shall under-estimate the transport capability when utilizing the concentration transport approach. Besides, sandbank packing process in the concentration transport results also differ a lot from that in the MP-PIC results, including the slopes of upstream and downstream sandbank region. This is because particle-particle interaction is considered in the MP-PIC method in some way, while it is not considered in the concentration transport approach.

It is worth noting here that in this work the fifth order WENO (weighted essentially non-oscillation) scheme is adopted for solving the concentration transport equation, thus the effect of numerical diffusion of concentration transport approach can be considered insignificant, and differences between the numerical

*Comparison of Concentration Transport Approach and MP-PIC Method for Simulating… DOI: http://dx.doi.org/10.5772/intechopen.92227*

results of two methods can exactly reflect the effects of different modeling strategy on proppant transport. Above all, both methods can precisely capture the proppant settling mechanism when the same drag force model is adopted in both methods, and concentration transport approach can capture the gravity convection mechanisms precisely when Reynolds number is smaller than 100. However, proppant packing mechanism is not captured very well in the concentration approach.

## **4. Conclusion**

two cases still evolve and both trend from vertical direction to inclined direction. This is because of the horizontal pressure gradient on the slurry front due to the difference between slurry and pure fluid, which provides a driving force and keeps slurry on the bottom penetrating into the pure fluid region. This mechanism is then intuitively described as gravity convection. Gravity convection in case III is much weaker than that in case II. Though the horizontal pressure gradients on the slurry front are the same in these two cases, fluid viscosity is greater in case III, which leads to a smaller channel permeability for slurry flows. Therefore, the penetration velocity in case III is much smaller than that in case II. In case I, gravity convection also plays a significant role. According to the ratio of inlet velocity and proppant settling velocity, that is, about 4, without gravity convection slurry front is

expected to be exact the diagonal line of the domain. However, numerical results of both methods show that the suspending region is far below the diagonal line, which indicates that fluid velocity field is greatly changed due to gravity convection compared to a uniform flow. Above all, gravity convection can be considered as a global effect reflecting the density difference between slurry and pure fluid, whereas proppant settling represents a local effect reflecting the density difference

The third one is proppant packing. The two prior mechanisms affect the distribution of slurry suspending region, and proppant packing shall affect the distribution of sandbank. In case I, as time evolves proppant concentration increases at the bottom of the simulation domain. When proppant concentration reaches a maximum value, that is, close to the packing state, particle-particle interactions and particle-wall interactions become more frequent. Then, the early formed sandbank stays unmovable, and the fluid velocity also dramatically decreases in this region due to the fluid-particle coupling and particle-particle/wall damping effects.

By comparing the numerical results of two different methods based on the above mechanisms, similarities and differences between the two methods discussed in

Firstly, in case III it is obvious that numerical results of methods are almost the

same with each other. In case III, Reynolds number is pretty low and the fluid motion is dominated by the viscous mechanism. The slurry flow is approximately plug flow. Secondly, in case II of moderate Reynolds number, numerical results of concentration transport approach show that the gravity convection is a bit severe than those of MP-PIC method. However, transport patterns of these two methods are in general similar. Quantitatively, the relative length of top and bottom slurry fronts at the end of simulation time is 0.5 and 0.4 m, respectively. Thirdly, in case I, numerical results of two methods differ a lot from each other. For the slurry suspending region, the covering area of the MP-PIC results is obviously much larger than those of concentration transport approach. This is mainly because the transient term and convection term in the fluid governing equation are ignored in the concentration approach, and these two mechanisms play significant roles in low viscous or high Reynolds number cases. It is clear that losses of these two mechanisms shall under-estimate the transport capability when utilizing the concentration transport approach. Besides, sandbank packing process in the concentration transport results also differ a lot from that in the MP-PIC results, including the slopes of upstream and downstream sandbank region. This is because particle-particle interaction is considered in the MP-PIC method in some way, while it is not considered in the

It is worth noting here that in this work the fifth order WENO (weighted essentially non-oscillation) scheme is adopted for solving the concentration transport equation, thus the effect of numerical diffusion of concentration transport approach can be considered insignificant, and differences between the numerical

between proppant particle and pure fluid.

*Progress in Fine Particle Plasmas*

Section 2.3 are verified.

concentration transport approach.

**188**

In this work, two numerical methods are adopted to simulate proppant transport: the concentration transport approach and the MP-PIC method. The first one is a typical Eulerian-Eulerian method and the second one is a typical Eulerian-Lagrangian method. With full discussions on their frameworks and comparisons between the numerical results, the following conclusions are then obtained:


## **Acknowledgements**

This work is funded by the National Natural Science Foundation of China (nos. U1730111, 11972088, 91852207), the State Key Laboratory of Explosive Science and Technology (no. QNKT19-05), and the National Science and Technology Major Project of China (Grant no. 2017ZX05039-005).

*Progress in Fine Particle Plasmas*

## **Author details**

Junsheng Zeng<sup>1</sup> \* and Heng Li<sup>2</sup>


**References**

2001

9-12 October 2005

Sciences. 2014;**70**:273-285

[1] Novotny EJ. Proppant transport. In: Paper SPE 6813 Presented at the 52nd Annual Fall Technical Conference and Exhibition of the Society of Petroleum Engineers of AIME, Held in Denver,

*DOI: http://dx.doi.org/10.5772/intechopen.92227*

*Comparison of Concentration Transport Approach and MP-PIC Method for Simulating…*

flows. International Journal of Multiphase Flow. 2001;**27**:1685-1706

[9] O'Rourke PJ, Zhao P, Snider DM. A model for collisional exchange in gas/ liquid/solid fluidized beds. Chemical Engineering Science. 2009;**64**:1784-1797

[10] Zeng J, Li H, Zhang D. Numerical simulation of proppant transport in propagating fractures with the multiphase particle-in-cell method. Fuel.

Experimental and numerical modeling of convective proppant transport. In: Paper SPE 28564 Prepared for Presentation at the SPE 69th Annual Technical Conference and Exhibition, Held in New Orleans, LA, USA, 25-28

2019;**245**:316-335

September 1994

100-111

[11] Barree RD, Conway MW.

[12] Wen C, Yu Y. Mechanics of

fluidization. The Chemical Engineering Progress Symposium Series. 1966;**162**:

[2] Mobbs AT, Hammond PS. Computer simulations of proppant transport in a hydraulic fracture. In: Paper SPE 29649 Prepared for Presentation at the 1995 SPE Western Regional Meeting Held in Bakersfield, California, 8-10 March

[3] Gadde PB, Sharma MM. The impact of proppant retardation on propped fracture lengths. In: Paper SPE 97106 Prepared for Presentation at the 2005 SPE Annual Technical Conference and Exhibition Held in Dallas, Texas, USA,

[4] Gu M, Mohanty KK. Effect of foam quality on effectiveness of hydraulic fracturing in shales. International Journal of Rock Mechanics and Mining

[5] Wang J, Elsworth D, Denison MK. Propagation, proppant transport and the evolution of transport properties of hydraulic fractures. Journal of Fluid Mechanics. 2018;**855**:503-534

[6] Dontsov EV, Peirce AP. A lagrangian

[7] Roostaei M, Nouri A, Fattahpour V, Chan D. Numerical simulation of proppant transport in hydraulic fractures. Journal of Petroleum Science and Engineering. 2018;**163**:

[8] Patankar NA, Joseph DD. Lagrangian numerical simulation of particulate

approach to modelling proppant transport with tip screen-out in KGD hydraulic fractures. Rock Mechanics and Rock Engineering. 2015;**48**:

2541-2550

119-138

**191**

Colorado, 9-12 October 1977

\*Address all correspondence to: zengjs@pku.edu.cn

© 2020 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/ by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

*Comparison of Concentration Transport Approach and MP-PIC Method for Simulating… DOI: http://dx.doi.org/10.5772/intechopen.92227*

## **References**

[1] Novotny EJ. Proppant transport. In: Paper SPE 6813 Presented at the 52nd Annual Fall Technical Conference and Exhibition of the Society of Petroleum Engineers of AIME, Held in Denver, Colorado, 9-12 October 1977

[2] Mobbs AT, Hammond PS. Computer simulations of proppant transport in a hydraulic fracture. In: Paper SPE 29649 Prepared for Presentation at the 1995 SPE Western Regional Meeting Held in Bakersfield, California, 8-10 March 2001

[3] Gadde PB, Sharma MM. The impact of proppant retardation on propped fracture lengths. In: Paper SPE 97106 Prepared for Presentation at the 2005 SPE Annual Technical Conference and Exhibition Held in Dallas, Texas, USA, 9-12 October 2005

[4] Gu M, Mohanty KK. Effect of foam quality on effectiveness of hydraulic fracturing in shales. International Journal of Rock Mechanics and Mining Sciences. 2014;**70**:273-285

[5] Wang J, Elsworth D, Denison MK. Propagation, proppant transport and the evolution of transport properties of hydraulic fractures. Journal of Fluid Mechanics. 2018;**855**:503-534

[6] Dontsov EV, Peirce AP. A lagrangian approach to modelling proppant transport with tip screen-out in KGD hydraulic fractures. Rock Mechanics and Rock Engineering. 2015;**48**: 2541-2550

[7] Roostaei M, Nouri A, Fattahpour V, Chan D. Numerical simulation of proppant transport in hydraulic fractures. Journal of Petroleum Science and Engineering. 2018;**163**: 119-138

[8] Patankar NA, Joseph DD. Lagrangian numerical simulation of particulate

flows. International Journal of Multiphase Flow. 2001;**27**:1685-1706

[9] O'Rourke PJ, Zhao P, Snider DM. A model for collisional exchange in gas/ liquid/solid fluidized beds. Chemical Engineering Science. 2009;**64**:1784-1797

[10] Zeng J, Li H, Zhang D. Numerical simulation of proppant transport in propagating fractures with the multiphase particle-in-cell method. Fuel. 2019;**245**:316-335

[11] Barree RD, Conway MW. Experimental and numerical modeling of convective proppant transport. In: Paper SPE 28564 Prepared for Presentation at the SPE 69th Annual Technical Conference and Exhibition, Held in New Orleans, LA, USA, 25-28 September 1994

[12] Wen C, Yu Y. Mechanics of fluidization. The Chemical Engineering Progress Symposium Series. 1966;**162**: 100-111

**Author details**

*Progress in Fine Particle Plasmas*

Junsheng Zeng<sup>1</sup>

**190**

\* and Heng Li<sup>2</sup>

1 College of Engineering, Peking University, Beijing, China

\*Address all correspondence to: zengjs@pku.edu.cn

provided the original work is properly cited.

2 School of Earth Resources, China University of Geosciences, Wuhan, China

© 2020 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/ by/3.0), which permits unrestricted use, distribution, and reproduction in any medium,

Section 6

Neutrinos and Galaxy

Custering

**193**

Section 6
