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

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.

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 high-strength 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 non-oscillation) 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 large-scale 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 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 ${\overrightarrow{u}}_f$, particle velocity ${\overrightarrow{u}}_p$, 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:

Figure 1.

Control volume in concentration transport approach.

$$\begin{array}{l}\frac{C^{n+1}{w}^{n+1}-C^n{w}^n}{\mathit{\Delta t}}=\frac{{\left({\mathit{Cwu}}_{p,x}\right)}_W\cdot \mathit{dy}-{\left({\mathit{Cwu}}_{p,x}\right)}_E\cdot \mathit{dy}}{\mathit{dx}}\\ +\frac{{\left({\mathit{Cwu}}_{p,y}\right)}_N\cdot \mathit{dx}-{\left({\mathit{Cwu}}_{p,y}\right)}_S\cdot \mathit{dx}}{\mathit{dy}}\end{array},$$E1

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

$$\frac{\partial \left(\mathit{Cw}\right)}{\partial t}+\nabla \cdot \left(\mathit{Cw}{\overrightarrow{u}}_p\right)=0.$$E2

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

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

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

$$\frac{\partial \left(\left(1-C\right)w\right)}{\partial t}+\nabla \cdot \left(\left(1-C\right)w{\overrightarrow{u}}_f\right)={\mathit{wu}}_{\mathit{\text{leak}}},$$E4

where $u_{\mathit{\text{leak}}}$ 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:

$$u_{p,x}=u_{f,x},u_{p,y}=u_{f,y}+u_{\mathit{\text{settling}}},$$E5

where $u_{\mathit{\text{settling}}}$ 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 \left(w\right)}{\partial t}+\nabla \cdot \left(w{\overrightarrow{u}}_s\right)={\mathit{wu}}_{\mathit{\text{leak}}},$$E6

where ${\overrightarrow{u}}_s=\left(1-C\right){\overrightarrow{u}}_f+C{\overrightarrow{u}}_p={\overrightarrow{u}}_f+{\mathit{Cu}}_{\mathit{\text{settling}}}\overrightarrow{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{u}}_f=-\frac{w^2}{12{\mu }_f}\nabla \left(P+{\rho }_f\mathit{gz}\right),$$E7

where ${\mu }_f$ 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{u}}_s=-\frac{w^2}{12\mu \left(C\right)}\nabla \left(P+{\rho }_s\mathit{gz}\right),$$E8

where $\mu \left(C\right)$ is the effective viscosity, which depends on particle concentration $C$ and obviously $\mu \left(0\right)={\mu }_f$, ${\rho }_s={\rho }_f\left(1-C\right)+{\rho }_{p}C$ 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 x}\left(\frac{w^3}{12\mu \left(C\right)}\frac{\partial P}{\partial x}\right)-\frac{\partial }{\partial z}\left(\frac{w^3}{12\mu \left(C\right)}\frac{\partial \left(P+{\rho }_s\mathit{gz}\right)}{\partial z}\right)={\mathit{wu}}_{\mathit{\text{leak}}}$$.E9

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 \left({\alpha }_{f}w\right)}{\partial t}+\nabla \cdot \left({\alpha }_{f}w{\overrightarrow{u}}_f\right)={\mathit{wu}}_{\mathit{\text{leak}}},$$E10

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

where w is the fracture width, ${\alpha }_f$ is the volume fraction of fluid, ${\overrightarrow{f}}_w$ is the wall friction term and ${\overrightarrow{f}}_p$ is the fluid-particle coupling source term, ${\overrightarrow{u}}_f$ is the fluid velocity, symbol “$\otimes$” 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}{\mathit{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\overline{\overline{I}}-{\overline{\overline{\tau }}}_f\right)-\frac{\nabla {\tau }_p}{{\alpha }_p{\rho }_p}-\frac{1}{{\tau }_w}{\overrightarrow{u}}_p$$E12

where ${\overrightarrow{u}}_p$ is the parcel velocity, $D_p$ is the drag coefficient, ${\alpha }_p$ is the volume fraction of the particle phase, ${\tau }_p$ is the particle normal stress, and ${\tau }_w$ is the damping relaxation time due to particle-wall interaction, P is the local average fluid pressure, and ${\overline{\overline{\tau }}}_f$ is the local average fluid viscous stress tensor. Here we consider the Archimedes buoyancy force in the term $\left(1-\frac{{\rho }_f}{{\rho }_p}\right)\overrightarrow{g}$ 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.

1. Concentration transport approach:

$$\left\{\begin{array}{l}\frac{\partial w}{\partial t}-\frac{\partial }{\partial x}\left(\frac{w^3}{12\mu \left(C\right)}\frac{\partial P}{\partial x}\right)-\frac{\partial }{\partial z}\left(\frac{w^3}{12\mu \left(C\right)}\frac{\partial \left(P+{\rho }_s\mathit{gz}\right)}{\partial z}\right)={\mathit{wu}}_{\mathit{\text{leak}}}\\ \frac{\partial \left(\mathit{Cw}\right)}{\partial t}+\nabla \cdot \left(\mathit{Cw}{\overrightarrow{u}}_p\right)=0\\ u_{p,x}=u_{f,x},u_{p,y}=u_{f,y}+u_{\mathit{\text{settling}}}\\ {\overrightarrow{u}}_s=\left(1-C\right){\overrightarrow{u}}_f+C{\overrightarrow{u}}_p={\overrightarrow{u}}_f+{\mathit{Cu}}_{\mathit{\text{settling}}}\overrightarrow{j}=-\frac{w^2}{12\mu \left(C\right)}\nabla \left(P+{\rho }_s\mathit{gz}\right)\end{array}\right.$$E13

MP-PIC approach:

$$\left\{\begin{array}{l}\frac{\partial \left({\alpha }_{f}w\right)}{\partial t}+\nabla \cdot \left({\alpha }_{f}w{\overrightarrow{u}}_f\right)={\mathit{wu}}_{\mathit{\text{leak}}}\\ \frac{1}{w}\frac{\partial \left({\alpha }_{f}w{\overrightarrow{u}}_f\right)}{\partial t}+\frac{1}{w}\nabla \cdot \left({\alpha }_{f}w{\overrightarrow{u}}_f\otimes {\overrightarrow{u}}_f\right)=-\frac{1}{{\rho }_f}\nabla P-f_{\mathit{\text{wall}}}+f_{\mathit{\text{particle}}}\\ \frac{d{\overrightarrow{u}}_p}{\mathit{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\overline{\overline{I}}-{\overline{\overline{T}}}_f\right)-\frac{1}{{\theta }_p{\rho }_p}\nabla {\tau }_p-\frac{1}{{\tau }_w}{\overrightarrow{u}}_p\end{array}\right..$$E14

In Eq. (13), we need to introduce models for settling velocity $u_{\mathit{\text{settling}}}$ and effective viscosity $\mu \left(C\right)$. In Eq. (14), we need to model the two terms: wall friction $f_{\mathit{\text{wall}}}$ and fluid-particle interaction force $f_{\mathit{\text{particle}}}$.

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

1. 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 ⟨P⟩=-f_{\mathit{\text{wall}}}-{\alpha }_p\left({\rho }_p-{\rho }_f\right)\overrightarrow{g}$$E15

   If the term $f_{\mathit{\text{wall}}}$ converges into the formula of $\frac{12\mu \left({\alpha }_p\right)}{{\rho }_f{w}^2}{\overrightarrow{u}}_s$ we can recover the Eq. (8) without hydrostatic pressure. Note that here $⟨{\overrightarrow{u}}_s⟩={\alpha }_f⟨{\overrightarrow{u}}_f⟩+\frac{1}{V_f}\sum V_p{\overrightarrow{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.

2. 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({\alpha }_p\right)=\frac{{\mu }_0}{{\left(1-\frac{{\alpha }_p}{{\alpha }_{\mathit{cp}}}\right)}^{1.82}},$$E16

where ${\alpha }_p$ is the particle volume fraction, or proppant concentration, ${\alpha }_{\mathit{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 fluid-particle 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{3{\rho }_f}{8{\rho }_p}{C}_d\frac{\left|{\overrightarrow{u}}_f-{\overrightarrow{u}}_p\right|}{r_p},$$E17

$$C_d=\left\{\begin{array}{cc}\frac{24}{{Re}_p}\left(1+0.15{{Re}_p}^{0.687}\right){\alpha }_f^{-2.65}& {Re}_p<1,000\\ 0.44{\alpha }_f^{-2.65}& {Re}_p\ge 1,000\end{array}\right.,$$E18

where ${Re}_p=\frac{2{\rho }_f\left|{\overrightarrow{u}}_f-{\overrightarrow{u}}_p\right|r_p}{{\mu }_f}$ is the particle Reynolds number, and $r_p$ is the physical radius of proppant particle.

3. Case study

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.

Figure 2.

A sketch for the numerical simulation domain.

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.

Parameter	Value	Parameter	Value
Fluid density	1000 kg/m^3	Particle density	2500 kg/m^3
Fluid viscosity	1 and 10 and 100 cP	Particle diameter	0.6 mm
Time-step	$5\times {10}^{-3}$ s	Simulation time	4 s
Inlet velocity	0.2 m/s	Inlet concentration	20%
Domain size	$1\times 0.25\times 0.005{\mathrm{m}}^3$	Mesh size	$100\times 25$

Table 1.

Parameters for simulating proppant transport.

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.

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.

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.

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 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 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.

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 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

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:

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 pure fluid can be clearly captured.

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 sandbank packing process.

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 behaviors.

Acknowledgments

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).