Proppant Transport Solver

Introduction

The ProppantTransport solver applies the finite volume method to solve the equations of proppant transport in hydraulic fractures. The behavior of proppant transport is described by a continuum formulation. Here we briefly outline the usage, governing equations and numerical implementation of the proppant transport model in GEOS.

Theory

The following mass balance and constitutive equations are solved inside fractures,

Proppant-fluid Slurry Flow

$\frac{\partial}{\partial t}(\rho_m) + \boldsymbol{\nabla} \cdot (\rho_m \boldsymbol{u_m}) = 0,$

where the proppant-fluid mixture velocity $\boldsymbol{u_m}$ is approximated by the Darcy’s law as,

$\boldsymbol{u}_m = -\frac{K_f}{\mu_m}(\nabla p - \rho_m \boldsymbol{g}),$

and $p$ is pressure, $\rho_m$ and $\mu_m$ are density and viscosity of the mixed fluid , respectively, and $\boldsymbol{g}$ is the gravity vector. The fracture permeability $K_f$ is determined based on fracture aperture $a$ as

$K_f = \frac{a^2}{12}$

Proppant Transport

$\frac{\partial}{\partial t}(c) + \boldsymbol{\nabla} \cdot (c \boldsymbol{u}_p) = 0,$

in which $c$ and $\boldsymbol{u}_p$ represent the volume fraction and velocity of the proppant particles.

Multi-component Fluid Transport

$\frac{\partial}{\partial t} [ \rho_i \omega_i (1 - c) ] + \boldsymbol{\nabla} \cdot [ \rho_i \omega_i (1 - c) \boldsymbol{u}_f ] = 0.$

Here $\boldsymbol{u}_f$ represents the carrying fluid velocity. $\rho_i$ and $\omega_i$ denote the density and concentration of i-th component in fluid, respectively. The fluid density $\rho_f$ can now be readily written as

$\rho_f = \sum_{i=1}^{N_c} \rho_i \omega_i,$

where $N_c$ is the number of components in fluid. Similarly, the fluid viscosity $\mu_f$ can be calculated by the mass fraction weighted average of the component viscosities.

The density and velocity of the slurry fluid are further expressed as,

$\rho_m = (1 - c) \rho_f + c \rho_p,$

and

$\rho_m \boldsymbol{u}_m = (1 - c) \rho_f \boldsymbol{u}_f + c \rho_p \boldsymbol{u}_p,$

in which $\rho_f$ and $\boldsymbol{u}_f$ are the density and velocity of the carrying fluid, and $\rho_p$ is the density of the proppant particles.

Proppant Slip Velocity

The proppant particle and carrying fluid velocities are related by the slip velocity $\boldsymbol{u}_{slip}$ ,

$\boldsymbol{u}_{slip} = \boldsymbol{u}_p - \boldsymbol{u}_f.$

The slip velocity between the proppant and carrying fluid includes gravitational and collisional components, which take account of particle settling and collision effects, respectively.

The gravitational component of the slip velocity $\boldsymbol{u}_{slipG}$ is written as a form as

$\boldsymbol{u}_{slipG} = F(c) \boldsymbol{u}_{settling},$

where $\boldsymbol{u}_{settling}$ is the settling velocity for a single particle, $d_p$ is the particle diameter, and $F(c)$ is the correction factor to the particle settling velocity in order to account for hindered settling effects as a result of particle-particle interactions,

$F(c) = e^{-\lambda_s c},$

with the hindered settling coefficient $\lambda_s$ as an empirical constant set to 5.9 by default (Barree & Conway, 1995).

The settling velocity for a single particle, $\boldsymbol{u}_{settling}$ , is calculated based on the Stokes drag law by default,

$\boldsymbol{u}_{settling} = ( \rho_p - \rho_f) \frac{d{_p}^{2}}{18 \mu_f}\boldsymbol{g}.$

Single-particle settling under intermediate Reynolds-number and turbulent flow conditions can also be described respectively by the Allen’s equation (Barree & Conway, 1995),

$\boldsymbol{u}_{settling} = 0.2 d_{p}^{1.18} \left [ \frac{g ( \rho_p - \rho_f)}{\rho_f} \right ]^{0.72} \left ( \frac{\rho_f}{\mu_f} \right )^{0.45} \boldsymbol{e},$

and Newton’s equation(Barree & Conway, 1995),

$\boldsymbol{u}_{settling} = 1.74 d{_p}^{0.5}\left [ \frac{g ( \rho_p - \rho_f)}{\rho_f}\right]^{0.5} \boldsymbol{e}.$

$\boldsymbol{e}$ is the unit gravity vector and $d_p$ is the particle diameter.

The collisional component of the slip velocity is modeled by defining $\lambda$ , the ratio of the particle velocity to the volume averaged mixture velocity as a function of the proppant concentration. From this the particle slip velocity in horizontal direction is related to the mixed fluid velocity by,

$\boldsymbol{u}_{slipH} = \frac{\lambda - 1}{1 - c} \boldsymbol{v}_{m}$

with $\boldsymbol{v}_{m}$ denoting volume averaged mixture velocity. We use a simple expression of $\lambda$ proposed by Barree & Conway (1995) to correct the particle slip velocity in horizontal direction,

$\lambda= \left[\alpha - |c - c_{slip} |^{\beta} \right]\,$

where $\alpha$ and $\beta$ are empirical constants, $c_{slip}$ is the volume fraction exhibiting the greatest particle slip. By default the model parameters are set to the values given in (Barree & Conway, 1995): $\alpha= 1.27$ , $c_{slip} =0.1$ and $\beta = 1.5$ . This model can be extended to account for the transition to the particle pack as the proppant concentration approaches the jamming transition.

Proppant Bed Build-up and Load Transport

In addition to suspended particle flow the GEOS has the option to model proppant settling into an immobile bed at the bottom of the fracture. As the proppant cannot settle further down the proppant bed starts to form and develop at the element that is either at the bottom of the fracture or has an underlying element already filled with particles. Such an “inter-facial” element is divided into proppant flow and immobile bed regions based on the proppant-pack height.

Although proppant becomes immobile fluid can continue to flow through the settled proppant pack. The pack permeability K is defined based on the Kozeny-Carmen relationship:

$K = \frac{(sd_p)^2}{180}\frac{\phi^{3}}{(1-\phi)^{2}}$

and

$\phi = 1 - c_{s}$

where $\phi$ is the porosity of particle pack and $c_{s}$ is the saturation or maximum fraction for proppant packing, $s$ is the sphericity and $d_p$ is the particle diameter.

The growth of the settled pack in an “inter-facial” element is controlled by the interplay between proppant gravitational settling and shear-force induced lifting as (Hu et al., 2018),

$\frac{d H}{d t} = \frac{c u_{settling} F(c)}{c_{s}} - \frac{Q_{lift}}{A c_{s}},$

where $H$ , $t$ , $c_{s}$ , $Q_{lift}$ , and $A$ represent the height of the proppant bed, time, saturation or maximum proppant concnetration in the proppant bed, proppant-bed load (wash-out) flux, and cross-sectional area, respectively.

The rate of proppant bed load transport (or wash out) due to shear force is calculated by the correlation proposed by Wiberg and Smith (1989) and McClure (2018),

$Q_{lift} = a \left ( d{_p} \sqrt{\frac{g d{_p} ( \rho_p - \rho_f)}{\rho_f}} \right ) (9.64 N_{sh}^{0.166})(N_{sh} - N_{sh, c})^{1.5}.$

$a$ is fracture aperture, and $N_{sh}$ is the Shields number measuring the relative importance of the shear force to the gravitational force on a particle of sediment (Miller et al., 1977; Biot & Medlin, 1985; McClure, 2018) as

$N_{sh} = \frac{\tau}{d{_p} g ( \rho_p - \rho_f)},$

and

$\tau = 0.125 f \rho_f u_{m}^2$

where $\tau$ is the shear stress acting on the top of the proppant bed and $f$ is the Darcy friction coefficient. $N_{sh, c}$ is the critical Shields number for the onset of bed load transport.

Proppant Bridging and Screenout

Proppant bridging occurs when proppant particle size is close to or larger than fracture aperture. The aperture at which bridging occurs, $h_{b}$ , is defined simply by

$h_{b} = \lambda_{b} d_p,$

in which $\lambda_{b}$ is the bridging factor.

Slurry Fluid Viscosity

The viscosity of the bulk fluid, $\mu_m$ , is calculated as a function of proppant concentration as (Keck et al., 1992),

$\mu_{m} = \mu_{f}\left [1 + 1.25 \left ( \frac{c}{1-c/c_{s}} \right) \right ]^{2}.$

Note that continued model development and improvement are underway and additional empirical correlations or functions will be added to support the above calculations.

Spatial Discretization

The above governing equations are discretized using a cell-centered two-point flux approximation (TPFA) finite volume method. We use an upwind scheme to approximate proppant and component transport across cell interfaces.

Solution Strategy

The discretized non-linear slurry flow and proppant/component transport equations at each time step are separately solved by the Newton-Raphson method. The coupling between them is achieved by a time-marching sequential (operator-splitting) solution approach.

Parameters

The solver is enabled by adding a <ProppantTransport> node and a <SurfaceGenerator> node in the Solvers section. Like any solver, time stepping is driven by events, see Event Management.

The following attributes are supported:

XML Element: ProppantTransport

Name	Type	Default	Description
allowNegativePressure	integer	1	Flag indicating if negative pressure is allowed
bridgingFactor	real64	0	Bridging factor used for bridging/screen-out calculation
cflFactor	real64	0.5	Factor to apply to the CFL condition when calculating the maximum allowable time step. Values should be in the interval (0,1]
criticalShieldsNumber	real64	0	Critical Shields number
discretization	groupNameRef	required	Name of discretization object (defined in the Numerical Methods) to use for this solver. For instance, if this is a Finite Element Solver, the name of a Finite Element Discretization should be specified. If this is a Finite Volume Method, the name of a Finite Volume Discretization discretization should be specified.
frictionCoefficient	real64	0.03	Friction coefficient
initialDt	real64	1e+99	Initial time-step value required by the solver to the event manager.
isThermal	integer	0	Flag indicating whether the problem is thermal or not.
logLevel	integer	0	Sets the level of information to write in the standard output (the console typically). Level 0 outputs no specific information for this solver. Higher levels require more outputs. 1 - Line search information - Solution information (scaling, maximum changes, quality check) - Convergence information - Time step information - Linear solver information - Nonlinear solver information - Solver timers information 2 - The summary of declared fields and coupling
maxAbsolutePressureChange	real64	-1	Maximum (absolute) pressure change in a Newton iteration
maxProppantConcentration	real64	0.6	Maximum proppant concentration
maxSequentialPressureChange	real64	100000	Maximum (absolute) pressure change in a sequential iteration, used for outer loop convergence check
maxSequentialTemperatureChange	real64	0.1	Maximum (absolute) temperature change in a sequential iteration, used for outer loop convergence check
name	groupName	required	A name is required for any non-unique nodes
proppantDensity	real64	2500	Proppant density
proppantDiameter	real64	0.0004	Proppant diameter
targetRegions	groupNameRef_array	required	Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager.
updateProppantPacking	integer	0	Flag that enables/disables proppant-packing update
writeLinearSystem	integer	0	Write matrix, rhs, solution to screen ( = 1) or file ( = 2).
LinearSolverParameters	node	unique	XML Element: LinearSolverParameters
NonlinearSolverParameters	node	unique	XML Element: NonlinearSolverParameters

In particular:

discretization must point to a Finite Volume flux approximation scheme defined in the Numerical Methods section of the input file (see Finite Volume Discretization)
proppantName must point to a particle fluid model defined in the Constitutive section of the input file (see Constitutive Models)
fluidName must point to a slurry fluid model defined in the Constitutive section of the input file (see Constitutive Models)
solidName must point to a solid mechanics model defined in the Constitutive section of the input file (see Constitutive Models)
targetRegions attribute is currently not supported, the solver is always applied to all regions.

Primary solution field labels are proppantConcentration and pressure. Initial conditions must be prescribed on these field in every region, and boundary conditions must be prescribed on these fields on cell or face sets of interest. For static (non-propagating) fracture problems, the fields ruptureState and elementAperture should be provided in the initial conditions.

In addition, the solver declares a scalar field named referencePorosity and a vector field named permeability, that contains principal values of the symmetric rank-2 permeability tensor (tensor axis are assumed aligned with the global coordinate system). These fields must be populated via XML Element: FieldSpecification section and permeability should be supplied as the value of coefficientName attribute of the flux approximation scheme used.

Example

First, we specify the proppant transport solver itself and apply it to the fracture region:

    <ProppantTransport
      name="ProppantTransport"
      logLevel="1"
      discretization="singlePhaseTPFA"
      targetRegions="{ Fracture }">
      <NonlinearSolverParameters
        newtonTol="1.0e-8"
        newtonMaxIter="8"
        lineSearchAction="None"/>
      <LinearSolverParameters
        directParallel="0"/>
    </ProppantTransport>

Then, we specify a compatible flow solver (currently a specialized SinglePhaseProppantFVM solver must be used):

    <SinglePhaseProppantFVM
      name="SinglePhaseFVM"
      logLevel="1"
      discretization="singlePhaseTPFA"
      targetRegions="{ Fracture }">
      <NonlinearSolverParameters
        newtonTol="1.0e-8"
        newtonMaxIter="8"
        lineSearchAction="None"/>
      <LinearSolverParameters
        solverType="gmres"
        krylovTol="1.0e-12"/>
    </SinglePhaseProppantFVM>

Finally, we couple them through a coupled solver that references the two above:

    <FlowProppantTransport
      name="FlowProppantTransport"
      proppantSolverName="ProppantTransport"
      flowSolverName="SinglePhaseFVM"
      targetRegions="{ Fracture }"
      logLevel="1">
      <NonlinearSolverParameters
        newtonMaxIter="8"
	lineSearchAction="None"  	  
        couplingType="Sequential"/>
    </FlowProppantTransport>
        

References

1. 1. Barree & M. W. Conway. “Experimental and numerical modeling of convective proppant transport”, JPT. Journal of petroleum technology, 47(3):216-222, 1995.
1. 1. Biot & W. L. Medlin. “Theory of Sand Transport in Thin Fluids”, Paper presented at the SPE Annual Technical Conference and Exhibition, Las Vegas, NV, 1985.
1. Hu, K. Wu, X. Song, W. Yu, J. Tang, G. Li, & Z. Shen. “A new model for simulating particle transport in a low-viscosity fluid for fluid-driven fracturing”, AIChE J. 64 (9), 35423552, 2018.
1. 1. Keck, W. L. Nehmer, & G. S. Strumolo. “A new method for predicting friction pressures and rheology of proppant-laden fracturing fluids”, SPE Prod. Eng., 7(1):21-28, 1992.
1. McClure. “Bed load proppant transport during slickwater hydraulic fracturing: insights from comparisons between published laboratory data and correlations for sediment and pipeline slurry transport”, J. Pet. Sci. Eng. 161 (2), 599610, 2018.
1. 1. Miller, I. N. McCave, & P. D. Komar. “Threshold of sediment motion under unidirectional currents”, Sedimentology 24 (4), 507527, 1977.
1. 1. Wiberg & J. D. Smith. “Model for calculating bed load transport of sediment”, J. Hydraul. Eng. 115 (1), 101123, 1989.