International Journal of Rock Mechanics & Mining Sciences 154 (2022) 105096

Contents lists available at ScienceDirect

International Journal of Rock Mechanics and Mining Sciences

journal homepage: www.elsevier.com/locate/ijrmms

Numerical investigation of sand production mechanisms in weak sandstone

formations with various reservoir fluids
Furkhat Khamitov a, *, Nguyen Hop Minh b, Yong Zhao a
a
Nazarbayev University, Nur-Sultan, 010000, Kazakhstan
b
Fulbright University Vietnam, Ho Chi Minh City, 700000, Viet Nam

A R T I C L E I N F O A B S T R A C T

Keywords: A three-dimensional CFD-DEM simulation method was developed to investigate the influence of different
Sand production reservoir fluids on sand production. The fluid properties and reservoirs were prescribed and modelled after the
CFD-DEM conditions in the Kenkiyak and Uzen oil fields in Kazakhstan. A cylindrical sandstone sample was modelled using
Weak formation
the Discrete Element Method in three stages of the pluviation, diagenesis, and perforation processes as in the oil
Heavy and light oil
Fluid-particle interaction
field. A coupled simulation with the Computational Fluid Dynamics was used to simulate the sand production
from the perforated sandstone under different fluid flow conditions with light and heavy oils. A continuous sand
production was found to be more severe with the heavy oil flow due to its higher transport capability and a
microstructural sanding mechanism that mobilized a higher proportion of the sample’s particles with the fluid
outflow. Sand production for light oil occurred in two stages of a first transient sand production that was fol­
lowed by a second continuous sand production. The transient sand production was associated with more sig­
nificant bond breakage, particle movement and porosity change within the damage zone that was created by the
central perforation inside the sandstone sample.

1. Introduction The main controlling factors of sand production can be grouped into
three categories as: the physical properties of the reservoir formation
Hydrocarbons production (oil and gas) with sand materials is present and fluids, the well installation and completion, and the type of oil and
in many oil fields in Kazakhstan, such as Karazhanbas, North Buzachi, gas field development.8 It was shown that there are considerable dif­
Kalamkas Zhalgiztobe and Kenkiyak1. These hydrocarbon reservoirs are ferences between the sand production patterns for gas, light, and heavy
typically characterized by highly permeable and highly porous sand­ oil wells. The cumulative sand volume and decline period for heavy oil
stones with initial flooding usually above 20%, low depth, and viscous well can differ in hundred times than for light oil wells. Al-Awad et al.9
oil.2,3 Sand production can dramatically affect the well performance, conducted several sand production tests using Hoek cell to investigate
damage down hole and surface equipment, and raise the possibility of the effect of confining pressure, flow rate, and the displacing fluid vis­
critical failure.4 cosity on sand production in unconsolidated sandstones. The cumulative
Sand production can be prevented before well production by the sand production was found strongly dependent on both flow rate and
installation of downhole sand control systems such as gravel packs, sand confining pressure. As the rocks were saturated with water, or low vis­
screens and chemical consolidations, which could inhibit the movement cosity crude oil, sand production can be controlled by the flow rate. This
of sands into the underground and surface facilities.5 If no sand exclu­ was however not the case for heavy crude oil as sand production
sion is used to stop the sand production, the conventional sand man­ mechanism was different and therefore, controlling flow rate cannot
agement approach is to reduce the hydrocarbon production rate to stop sand production.
minimize the produced sands, and this can significantly reduce wellbore The type of displacing fluids had an indicative effect on sand pro­
productivity. On the other hand, when sands are extracted from reser­ duction initiation and final sand mass. Wu et al.10 conducted sand
voir, it increases the porosity close to the wellbore. The zone of production experiments on seven weak sandstones with various
increased permeability provides greater flow of liquid into the strengths. The experimental results for weakly-consolidated and
wellbore.2,6,7 consolidated sandstones showed that the predominant mechanism for

* Corresponding author.
E-mail address: Furkhat.khamitov@gmail.com (F. Khamitov).

https://doi.org/10.1016/j.ijrmms.2022.105096
Received 28 September 2021; Received in revised form 22 February 2022; Accepted 8 March 2022
1365-1609/© 2022 Elsevier Ltd. All rights reserved.
F. Khamitov et al. International Journal of Rock Mechanics and Mining Sciences 154 (2022) 105096

borehole failure was shear failure around the borehole. This was fol­ understood, especially for the detailed fluid-particle interactions. The
lowed by erosion of the loose materials due to flowing fluid in the region first coupled simulation on sand production with different fluid types in
of rock degradation (i.e., the failure zone). The study of particles this study posed a challenge in terms of the intensive computational
transportation by injected fluid is important for hydraulic fracturing. power that was optimized for the CFD mesh size and the coupling
Harris et al.11 classified the fracturing fluids according to the transport number between the computation of two separate phases. The heavy oil
capability. For the single fluids of the same viscosity (e.g., oil and water), was found to cause more severe sand production than light oil due to its
the transportation is due to viscosity alone and settling may occur during higher transport capability and the formation of a more regular particle
the transport depending upon flow rate and viscosity. Hughes et al.12 velocity trajectory pattern that mobilized a larger proportion of the
used a new viscoelastic surfactant (VES-CO2 ) fluid system as a foamed sample’s particles for sand production. Understanding the mechanisms
fracturing fluid to stimulate wells. Wells treated with VES-CO2 showed of sand production improves the ability to predict the production dy­
better cumulative production due to low formation damage, superior namics of both hydrocarbons and sand materials, which are key to safe
proppant transport capability. The viscosity of VES- CO2 is higher than and profitable operation.
that of the single fluid. Low viscosity fluid (e.g., fresh water) is not an
effective transport system for large particles with relatively high 2. Numerical formulation
densities.13
Sand production is a coupled fluid–solid interaction process that 2.1. DEM modelling of weak sandstone
primarily involves two mechanisms: mechanical instabilities that lead to
localized plastic behavior and failure of the rock around the cavity, and In the last decades, granular mechanics has been numerically studied
the subsequent transportation of sand particles due to the fluid-particle using the Discrete Element Method.37 In DEM, all particles in a
interaction. The behaviors can be modelled microscopically using computational domain are tracked explicitly by solving the governing
several numerical particle-fluid coupling techniques, among which are equations of particle motion, which provides the detailed information of
the two-fluid model (TFM) and the coupled computational fluid dy­ the interparticle contact forces and the particle’s trajectory. Materials
namics and discrete element method (CFD–DEM) as discussed by Zhu can be simulated in DEM in terms of the physical properties of the
et al.14 particles and the input parameters of the contact models, which allow
The origin of sand production modelling using of the CFD-DEM for different configurations of linear and nonlinear elasto-plastic
approach can be traced back to the simulation of bubble formation in behaviour at the contact points between discrete particles.38–46 The
gas-fluidized granular beds by Tsuji et al.15 O’Connor et al.16 introduced cemented rock material can be modelled as an assembly of fully, or
the application of DEM to model the mechanics of sand production. The partially bonded particles. The translational and angular accelerations
2D numerical simulation combined the Darcy flow model with DEM, of the particles are computed based on the momentum balances that
although the Darcy’s law is only valid for slow (small Reynolds number), consider the bonding effect in the calculation of the contact forces.
viscous and isothermal flow, the conditions that may not be met during In CFD-DEM simulation, a fluid phase is coupled with the particle
all sand production stages.17,18 Khamitov et al.19 conducted a large scale phase.14 In addition to the particle-particle contact forces, fluid-particle
three dimensional CFD-DEM simulation for sand production using the interaction forces are included in the second Newton’s law for a particle
open-source program CFDEM20 with 105 particles. The numerical p as follows20:
sandstone sample was perforated and then subjected to water injection
dup
as similar to the laboratory experiment that was modelled after the mp = f p,n + f p,t + f p,f + f p,Pres + f p,vis + f p,b , (1)
hydrocarbons production process in the field.21,22 The solid particles dt
were modelled with DEM using a modified JKR contact model23 to dωp
mimic weak cement bonding behavior of unconsolidated sandstone, Ip = rp,c × f p,t , (2)
dt
whereas the fluid part was modelled using the Navier-Stokes equations
that capture more accurately the large variations in flow velocities where up is linear and ωp - angular velocities, Ip – moment of inertia, mp -
during sand production. mass, rp,c -radius of the particle. f p,n and f p,t are the normal and tangential
In the CFD-DEM simulation, the interactions between the fluid and contact forces between particles, respectively; f p,f is the drag force
the particles are captured in terms of the particle–fluid drag force, exerted from the fluid phase to the particle; f p,Pres is pressure force, f p,vis
pressure gradient force, viscous force and other forces,24–27 of which the is viscous force, f p,b are other forces affecting the particles. Eq. (1)
drag force is the main driving force for fluidization. Many correlations
describe the translational motion, while Eq. (2) is used for the rotational
have been proposed to calculate the particle–fluid drag force.28–31 Kafui
motion of the particle p. Note that for dry materials, the fluid-related
et al.32 marked that the Di Felice’s28 expression leads to a smooth
force terms (f p,f , f p,Pres , f p,vis ) are equal to zero. The formulation of the
variation in the drag force as a function of fluid void fraction. Agrawal
drag force is more complicated, which will be presented in the next
et al.33 investigated the effect of the drag force models on CFD–DEM
section. The pressure force f p,Pres can be calculated as:
prediction of bubbling fluidized beds with Geldart-D34 type particles and
compared the numerical results with two sets of experimental data. It f p,Pres = − ∇P⋅V p (3)
was concluded that the Di Felice’s model provided better prediction
than the more conventional Gidaspow’s30 model in terms of the average where V p is volume of particle and ∇P is the fluid pressure gradient.
particle height for the bubbling fluidized bed systems. Marchelli et al.35 The viscous force f p,vis can be calculated as:
conducted the CFD-DEM simulation of two spouted beds containing
f p,vis = − (∇ ⋅ τ)V p (4)
Geldart-D type particles using different drag force models. Porosity and
the Reynolds number are the key parameters that affect the performance
where τ is the fluid stress tensor (see Eq. (9)).
of the drag force models; for the larger values, the Rong’s36 and the Di
The particle-particle forces, on the other hand, can be presented by
Felice’s models provided better prediction, while the other models
the spring-dashpot model,47 the normal force of which is given as:
overestimated the particles’ velocity.
In this study, we extend the model developed by Khamitov et al.19 for f p,n = − kp,n δp + cp,n Δup,n (5)
water-injection-induced sand production to investigate the sanding
mechanisms for different reservoir fluids. While varying sand produc­ where f p,n is the normal contact force of the particle as given in Eq. (1);
tion patterns have been observed for different fluid types in the field and Δup,n is the normal relative velocity at the contact point; kp,n and cp,n are
in previous research, the underlying mechanisms are not well the normal spring stiffness and damping coefficient; δp is the normal

2
F. Khamitov et al. International Journal of Rock Mechanics and Mining Sciences 154 (2022) 105096

overlap that for two particles i and j can be define as:

f p,f = f p,d + f ′′p (11)
δp = Ri + Rj − d (6)
where ΔV-volume of fluid cell, n is a number of particles in computa­
where Ri , Rj are the radii of particle i and j, respectively; d is the distance tional cell, f p,f is particle-fluid interaction force in Eq. (1); f p,d - a drag
between the particles’ centres. force and f ′′p is the sum of particle–fluid interaction forces on particles
The tangential force can be written as:
other than the drag, which are related to pressure gradient and viscous
⎧⃒ ⃒ ⎫
⎪ ⃒
⎨⃒ ∫t
⃒
⃒ ⎪
⎬ stress and are often regarded as the dominant forces in particle–fluid
⃒ ⃒
f p,t = min ⃒kp,t Δup,t dt + cp,t Δup,t ⃒, μc f p,n (7) flow.
⎪
⎩⃒⃒
tc,0
⃒
⃒
⎪
⎭ The drag force is based on the Di Felice’s correlation28:
⃒ ⃒( )
f p,d = 0.125Cd0,p ρf πd2p ε2p ⃒uf − up ⃒ uf − up ε−p χ (12)
where Δup,t is the relative tangential velocity of the particles in contact
as given in Eq. (1); kp,t and cp,t are the tangential spring stiffness and the where unknown χ , Cd0,p and Rep defined by:
tangential damping coefficient; μc is the Coulomb coefficient of friction. [ ( )2 / ]
The integral term represents an accumulated spring energy from the χ = 3.7 − 0.65exp − 1.5 − log 10 Rep 2 (13)
relative tangential motion between two particles in contact since the
time of contact formation tc,0 until the current time, t. The magnitude of ( / )2
the tangential force is limited by the Coulomb frictional limit, where the Cd0,p = 0.63 + 4.8 Re0.5
p (14)
particles begin to slide over each other. ⃒ ⃒/
The normal and tangential spring stiffnesses in Eqs. (5) and (7) are Rep = ρf dp εp ⃒uf − up ⃒ μf (15)
calculated based on a contact model for weak cemented sandstone
simulation. The effect of cementation is considered in the computation up is particle velocity, εp = 1 −
n
∑
V j /ΔV, dp is particle diameter, p =
of the normal contact force. For bonded contacts, the contact behaviour j=1
is specified by a modified JKR model.23 Once the bond is broken, it will 1 ÷ n.
not form again, and the contact behaviour is simulated by a Hertz Eq. (8) describes the particle fluid interactions in terms of the un­
contact. Further details on the calculation of contact forces can be found resolved CFD-DEM model. The CFD and DEM models were coupled by
in a similar simulation work by Khamitov et al.19 calculating the localized proportion of the CFD cell volume that is
occupied by the fluid, εf . The fluid cells should be configured to contain
2.2. Coupled CFD-DEM modelling of fluid-particle interaction at least several particles. The average mesh resolution dimension is
calculated as the ratio between the CFD grid size and particle diameter
The coupled modelling between the particle phase and the fluid (dp ) as:
phase in a sand production simulation was first time introduced by Tsuji
et al.15 For the coupled CFD-DEM simulation in this study, the fluid Δlcfd
≥Υ (16)
phase is considered as a continuum, while the solid phase is discrete dp
particles. The solid particles are fully resolved, and each of them is
tracked in the system. The unresolved approach is used for the CFD part where Δlcfd is the average length of the CFD cell in each direction and the
as the particle sizes are smaller than the computational grid of the fluid condition was checked for all three dimensions. The critical parameter Υ
phase.48 Individual particles do not completely fill a computational fluid is dependent on the void fraction scheme used in the simulation
cell. The coupling between solid phase at particle scale and fluid phase computation.50 The central void fraction scheme calculates the void
at fluid cell scale can be done via drag force and with empirical or fraction of a CFD cell by counting the whole volume of the particles
theoretical correlations developed for the flow across a single particle or whose centres are located inside the cell. For this case, the value of Υ
an assembly of particles.27,42,48,49 should be greater than.351–53 The central void fraction scheme becomes
Zhou et al.27 assessed different CFD-DEM models theoretically and inaccurate in the case of large particles as compared to the CFD cell size,
practically and recommended Model A for CFD–DEM modelling. In this which can be improved with the divided void fraction scheme. In this
model the governing equations of a fluid phase in the presence of a scheme, each particle is divided into 29 non-overlapping regions of
secondary particulate phase can be given as follows: equal volume. The centroid of each sub-region is then used to assign the
⎧ ( ) elementary volume of the particle to the corresponding computational
⎪
⎪ ∂ εf ( ) cell. A smaller minimum value of 1.86 can be assigned to Υ with the
⎪ + ∇⋅ εf uf = 0,
⎨
∂t divided void fraction scheme.54–56 The lowest limit of Υ in the literature
( ) (8)
⎪
⎪ ∂ εf ρf uf ( ) is 1.6, smaller values would affect the quantitative agreement between
⎪ + ∇⋅ εf ρf uf uf = − εf ∇P − Fset simulation and experimental results.50,51
⎩ II
pf + ε f ∇⋅ τ + ε f ρ f g
∂t
The information is exchanged between the separate computations of
where εf is the volume fraction occupied by the fluid, uf -fluid velocity, the solid phase and of the fluid phase through the coupling term fp,f of
ρf the fluid density, P the pressure, Fset II
– volumetric particle-fluid the particle-fluid interaction force. The DEM simulation provides the
pf
interaction force, g the gravity vector and t the time. The τ -stress updated locations of each and every particle for the calculation of the
tensor, that for fluid with viscosity μf is calculated as: volume fraction of fluid εf in each cell in Eq. (8). The CFD simulation
[( provides the updated pressure field and velocity field for the calculation
) ( )T ] /
τ = μf ∇uf + ∇uf − 2 3μf I∇⋅uf (9) of the drag force, the pressure force and the viscous force acting on each
particle that appear in both Eq. (1) and Eq. (8).
where I is identity matrix.
3. Computational model configuration
The volumetric particle–fluid interaction force Fset
pf
II
determine by:

1 ∑n
( ) 3.1. Perforated sample preparation
Fset II
= f (10)
pf
ΔV p=1 p,f
The numerical simulation was conducted in three subsequent stages

3
F. Khamitov et al. International Journal of Rock Mechanics and Mining Sciences 154 (2022) 105096

of the consolidated sample preparation, the perforation and the sand comprised of 2344 cells, the grid size of which was determined following
production. The initial sample was prepared by the particle generation the correlation in Eq. (16). In this case, the particle size was taken as the
in which the numerical particles were generated randomly inside a cy­ maximum particle diameter (dmax p = 0.355 mm), which resulted in the
lindrical container and let fall under gravity to the bottom to simulate a value Δlcfd should be greater than 0.6603 mm (Υ = 1.86 for divided void
pluvial deposition process. The number and the size of the particles were fraction scheme). In this study we used Δlcfd equal to 0.796 mm (Υ =
specified such that they follow the particle size distribution (PSD) of a 2.24).
real sandstone from the Ustyurt-Buzachi Sedimentary Basin.57 A total of The numerical models of the solid phase in Fig. 1 and of the fluid
105 frictional elastic auto-adhesive particles were generated to create phase in Fig. 2 were used for two separate computations of each phase,
the numerical sample, that is equivalent to a total mass of 1.578 g. The respectively. The exchanged information between the two computations
following mechanical properties of sand material were used19: Young’s is given in term of the interaction force, the pressure force, and the
modulus E = 70 GPa, Poisson’s ratio υ = 0.3, friction coefficient μc = viscous force in Eqs. (1) and (8), which have been explained in the
0.3, particle density ρp = 2500 kg/m3 , interface energy γ = 40 J/ m2 . previous section.
Further information on the contact bond model can be found in the other The boundary conditions are described in Fig. 2. The top and the
work by Khamitov et al.19 A summary of the particle information is bottom walls were set as impermeable walls with no-slip boundary
given in Table 1, while the numerical particle sizes are shown in Fig. 1. condition. There was a central opening on the top surface to serve as the
The open-source program CFDEM20,58,59 was used in the current fluid outlet. Water was injected from the whole circumferential wall of
study that combinesthe PISO60 solver of OpenFOAM61,62 and the DEM the cylinder at a constant velocity of 5⋅10− 5 m/s. The outlet pressure was
solver of LIGGGHTS.63,64 The divided void fraction model was used for maintained at a constant P = 1MPa for all cases. The fluids in the
fluid fraction calculation in CFD cell.54 In this section, all simulations simulation are considered incompressible. To check the computational
with dry (no fluid) sample were performed using LIGGGHTS convergence of the constructed CFD mesh, two additional meshes (fine
open-source program, where a new contact model for weak cemented and coarse) were tested, of which except the mesh size all other condi­
sandstone was implemented.19 tions are the same as the setup of the final selected mesh in Fig. 2.
The pluviated sample was then compressed to a pre-determined Three different mesh sizes are subjected to the same boundary con­
vertical stress of 1 MPa to achieve the intact consolidated sample. The ditions in Fig. 2. The radial distributions of the fluid velocity along a
diameter (D) of the intact sample is equal to 15.12 mm, and the height horizontal plane at 4.8 mm from the bottom of the sample are shown in
(H) is equal to 6.008 mm. Fig. 3. The maximum value of the fluid velocity at the outlet and the
A hollow tunnel was created by driving a penetrometer along the mean velocity at the final time step were calculated, the difference in
central axis of the intact sample to simulate the perforation process in these values between different mesh sizes were calculated as follows:
well completion. The penetrometer was removed once it reached the end ⃒
⃒ mean
⃒
⃒
plate. The perforation tunnel served as the drainage boundary to facil­ ⃒uf − umean
f b ⃒
mean
Δuf (%) = *100% (17)
itate fluid flow through the sandstone toward the outlet that is posi­ umean
f b
tioned at the top end of the tunnel or the centre of the circular top
⃒ ⃒
surface of the sandstone sample in Fig. 2. The perforation created ⃒ max ⃒
⃒uf − umaxf b⃒
damages to the sample’s initial structure, which affected the sand pro­ Δumax
f (%) = *100% (18)
duction process in the final stage of the microscopic simulation. Further umax
f b

details on the simulation processes can be found in the simulation works

by Khamitov et al.19 where umean
f b , uf b – are mean and maximum values of fluid velocity of
max

the final selected mesh size in Fig. 3.

Here, the discretized time step in the simulation was set as 10− 6 s,
3.2. Sand production simulation while the total simulation time was equal to 3⋅10− 4 s. The CFD time step
was calculated using the Courant–Friedrichs–Lewy condition67,68:
The perforated sample was used for sand production simulation
uf Δt
where different fluids were injected radially from the external boundary C= ≤ const (19)
Δlcfd
towards the vertical perforation tunnel to mimic the condition in the
experiments on a similar sandstone 21. The sand production simulation where Δt is the CFD time step. The value of const is usually less than 1
was conducted using the CFDEM program.20 To study sand production and can change with the method used to solve the discretized equation,
in different reservoirs in Kazakhstan, the properties of the injection especially depending on whether the method is explicit or implicit. If
fluids were varied to simulate light oil and heavy oil. The following oil C < 1, the fluid will move from one cell to another within one time step.
properties were used for sand production simulations 65: light oil (LO) While for C > 1, the fluid travels through two or more cells in each time
from Uzen oil field with density 849.16 kg/m3 and kinematic viscosity step and this can affect negatively to the numerical convergence. The
4.12⋅10− 6 m2 /s; heavy oil (HO) from Kenkiyak oil field with density Cmean and Cmax in Table 2 are the mean and the maximum values of the
927.09 kg/m3 and kinematic viscosity 1.62⋅10− 4 m2 /s. The values of the Courant number C at the end of the simulation.
fluid density and viscosity were inputted in Eqs. (8) and (9). The radial distribution of the fine and the selected medium mesh
For the CFD-DEM simulation of sand production, the fluid phase sizes are almost identical in Fig. 2, which reflects a mesh convergence
needs to be added in term of a CFD continuum mesh that has the same condition. The coarse mesh size, on the other hand, showed significant
dimension as the perforated sample of the solid phase in Fig. 1. The CFD discrepancy that could be as large as 81% as compared to the conver­
grid was constructed using the HELYS-OS graphical user interface.66 gence results of the other two mesh sizes. The medium mesh size was
Fig. 2 shows the computational domain of the fluid phase that is

Table 1
A summary of the particle numbers and sizes in numerical simulation.
Particle diameter, dp (mm) 0.15 0.18 0.2 0.22 0.25 0.275 0.3 0.355

Number of particles 20758 15529 13357 13825 17507 9929 5414 3681
Mass ratio (%) 5.8% 7.5% 8.9% 12.2% 22.7% 17.1% 12.1% 13.7%
Mass (milligram) 91.7 118.5 139.9 192.7 358.1 270.3 191.3 215.6

4
F. Khamitov et al. International Journal of Rock Mechanics and Mining Sciences 154 (2022) 105096

Fig. 1. The sample perforation.

Fig. 2. CFD simulation domain and boundary conditions.

Fig. 3. The fluid velocity distribution along the line at z = 4.8 mm from the bottom of the sample.

selected as the final size for further simulation as it achieved both the simulated by a four-way coupling method 48. The motion of the particles
convergence condition and the CFD-DEM condition in Eq. (16). affects the motion of the fluid, and vice versa. The coupling information
The interactions between the fluid and the particle phases were was given in terms of the pressure force, the viscous force and drag force.

5
F. Khamitov et al. International Journal of Rock Mechanics and Mining Sciences 154 (2022) 105096

Table 2
∑
mi
Mesh convergence analysis results. Vpi = Vj (20)
Parameters Symbol Fine Medium Coarse j=1

CFD cell length (mm) 0.657 0.796 1.008

Δlcfd
where Vj is the volume of the individual particle; mi-total number of
Average mesh resolution 1.85 2.24 2.83
particles in the ith zone, i = 1, 2,…,5. The porosity in each zone (φi ) can
Υ
Mean courant number Cmean 2.78⋅10− 7
2.33⋅10− 7 1.8⋅10− 7
Maximum courant number Cmax 2.61⋅10− 5
2.36⋅10− 5 1.56⋅10− 5 be determined as
Difference in the mean velocity Δumean
f 2.29 0 15.33
(%) Vi − Vpi
φi = (21)
Difference in the maximum Δumax
f 4.72 0 81.1 Vi
velocity (%)
where Vi is the volume of the ith zone. Porosity is the property of the solid
particle phase, which can be approximated as the void fraction in the
The time steps for the DEM and CFD simulations were 2⋅ 10− 8 s and
fluid phase that is different for each fluid cell.
10− 6 s, respectively. The CFD simulation was computed for every 50
(coupling number) of the DEM simulation and the coupling information 4. Results and discussion
was exchanged accordingly. This setting in sufficient accuracy yet al­
lows for higher computational efficiency. It took 19 days to complete the 4.1. Macro analysis of sand production
simulation time of 0.0606 s for a coupled simulation case in this study
using a highly powerful computer cluster of 8 cores and 500 GB of Fig. 5 shows the results of the cumulative produced sand mass for the
random-access memory. It was found that the raising of number threads heavy oil (HO) and the light oil (LO) over the same total simulation time
more than 8 does not significantly increase the calculation speed. of 0.0606 s. Sand production started earlier for LO at the marked time
step t1 as part of a transient sanding behavior that reached a peak sand
3.3. Post processing procedure for microstructure analyses mass at t2 and then decreased to zero over time. The transient behavior
was followed by a continuous sand production that started from the time
The microstructure of the specimen can be analysed in more detail as step t3 and was on-going at the end of the simulation. For the HO case,
the simulation domain is divided into five cylindrical rings as shown in there was no transient sanding behavior but only the continuous sand
Fig. 4. The ring thickness is equal to the diameter of the penetrometer, production from the time step t4. Although sand production was started
which is larger than in the previous simulation19 and allows for better later, the HO case produced more sand at the end of the simulation, that
statistical representation. The rings were counted as zones 1 to 5 from is the time step t5.
the central axis outwards, where the outermost ring was not taken in the The produced particles were counted as they crossed the level of
analysis as it was affected by the external boundary condition. These about one particle diameter above the top surface of the sample. If the
zones were only used in the post-simulation analysis to capture the produced particles moved upward to an even higher level of about four
change in the microscopic properties during the sand production pro­ particle diameters above the top cap, they were deleted from the
cess, but do not affect the DEM and the CFD computations. simulation as it was assumed that the particles had successfully escaped
Porosity of the granular specimen can be calculated as the ratio of the the sample. The transient cumulative sand mass of light oil decreased to
void volume over the total volume. The volume of the particles, Vpi , zero from t2 to t3 is due to the produced particles were not able to
located within each zone can be calculated as completely leave the sample but rather return to deposit on the top
surface. The escape fluid velocity at the exit point is not sufficient to lift
the particles up to the upper elevation to be completely removed from
the specimen.
Fig. 6 shows the results of the average particle velocity (HO_p, LO_p)
and the fluid velocity at outlet (HO_f, LO_f) for HO and LO cases. The
results were normalized by the final value of the average particle ve­
locity at t5 in the LO case as it is the smallest velocity among the four
final results; the normalization would allow comparison of the flow
dynamics between different phases and cases. As the flow began, the
fluid velocities fluctuated until a steady state of a constant velocity was
achieved for both cases. The fluctuation is much more significant for
heavy oil due to a higher inertia in this case.
The particle velocity is higher than the fluid velocity for LO until t2
when the difference becomes zero, and this is associated with the
transient sanding behavior. The particle and the fluid velocities for HO,
on the other hand, were similar, and there was no transient sand pro­
duction observed in this period. The viscosity force as determined in Eqs.
(4) and (9) would act as a greater damping force to prevent the transient
sand production for heavy oil.
The continuous sand production occurred in the steady state flow
condition when the fluid velocity is nearly a hundred times higher than
the average particle velocity, such that the hydrodynamical forces can
continuously transport particles out of the sample. The particle velocity
of LO decreased more quickly in the steady state and became less than
the HO particle velocity at the final time step, this would indicate a more
severe on-going sand production for the heavy oil as it was discussed in
the literature 9.
Before the beginning of the continuous sand production, the fluid
Fig. 4. Cylindrical rings of the simulation domain for microscopic analysis.

6
F. Khamitov et al. International Journal of Rock Mechanics and Mining Sciences 154 (2022) 105096

Fig. 5. Cumulative sand production for different reservoir fluids.

Fig. 6. Average particle velocity and fluid velocity at outlet. The velocities are normalized by the final average particle velocity of light oil, that is Upt5 = 1.34⋅
2
10− m/s.

velocity is higher than the particle velocity, but it is not sufficient to ( )

transport the particles out of the sample. The fluid velocity at the gd2p ρf − ρp
Uset = (22)
continuous sanding onset can be compared to the terminal settling ve­ 18μf ρf
locity, that is dependent on the particle size and can be determined from
the Stokes’ Law as follows69: where ρp is the particle density; the other parameters have been defined
in the previous equations. The settling velocities of eight particle sizes

Table 3
The cumulative sand production, and the settling velocity for different particle sizes.
Case Total mass, mg Diameter (mm) 0.15 0.18 0.2 0.22 0.25 0.275 0.3 0.355

LO[t1 − t3] 0.179 Quantity 3 1 2 2 4 1 0 0

Δmi (%) 7.4 4.3 11.7 15.6 45.8 15.2 0.0 0.0
LO[t3 − t5] 2.516 Quantity 28 30 23 18 33 12 9 6
Δmi (%) 4.9 9.1 9.6 10.0 26.8 13.0 12.6 14.0
HO[t4 − t5] 2.726 Quantity 33 33 25 19 35 17 11 4
Δmi (%) 5.3 9.2 9.6 9.7 26.3 17.0 14.3 8.6
Settling velocity for LO Uset (mm/s) 5.79 8.33 10.29 12.45 16.07 19.45 23.15 32.41
Settling velocity for HO Uset (mm/s) 0.13 0.19 0.23 0.28 0.36 0.43 0.51 0.72

7
F. Khamitov et al. International Journal of Rock Mechanics and Mining Sciences 154 (2022) 105096

are given in Table 3, which have different values due to different fluid 4.2. Micro analysis of sand production
densities and viscosities. The similar outlet velocity of 14.18 mm/ s was
observed for both the heavy and light oils at the onset of the continuous Figs. 8 and 9 show the porosity results for light oil and heavy oil,
sand production. The outlet velocity is greater than Uset max of the respectively. The initial porosity of zone 2 was smallest as it formed a
largest particles for heavy oil, whereas for the LO case it is capped be­ compacted zone of unbonded particles due to perforation damage.19 The
tween the two extreme settling velocities of the largest, and the smallest porosity of zone 5 did not change significantly as it was located farthest
particle. The transport capability of HO is higher than LO that would from the perforation tunnel. The particles moved from the outer zones
lead to more severe sand production. toward the central perforation tunnel under fluid flow effect, and this
The detailed information of the produced particles is given in decreased the porosity of zone 1 while increasing the porosity of zones 2,
Table 3. The transient sanding behavior of light oil was mostly associ­ 3 and 4. The particle movement occurred earlier for light oil, which was
ated with small particle sizes, among which the particles around the D50 associated with an earlier transient sand production. The continuous
size occupy half of the produced particles. The continuous sand pro­ sand production however leads to the most significant particle move­
duction, on the other hand, produced particles of all sizes for both HO ment and porosity change. There were fluctuations in the porosity
and LO. The produced sand mass in the HO case is greater than the values in this period due to particles escape from the specimen. The
combined sand masses of both the transient and the continuous sand porosity change is less significant for heavy oil than for light oil, which
productions in the LO case. would indicate less particle movement and microstructural change for
The particles of the sandstone sample can have multiple contacts the case of heavy oil.
with the surrounding neighbours, that are consisted of both the bonded Fig. 10 shows the variation of the coordination number calculated
and unbonded contacts. All the bonded contacts were only formed once separately for each cylindrical ring for light oil. The initial value of the
in the intact state using the modified JKR model. The mechanical im­ coordination number was smallest in zone 1 and increased toward the
pacts in the following perforation and sand production stages could boundary until zone 5, which reflects different levels of damage away
break the bonds and form new unbonded contacts that were modelled from the central perforation axis. As the transient sand production
using the Hertz contact model. The algorithm goes through the list of all started at t1, the coordination number in zone 1 dropped sharply to a
particles (k) in the sample and count the numbers of the modified JKR minimum value around t2, where there was also a maximum transient
contacts (bonded) and of the Hertz contacts (unbonded) for each particle produced sand mass. This contact reduction is associated with a decrease
to provide the total bond number and the total unbonded contact in porosity, or volume reduction, in zone 1 in Fig. 8. As the granular
number. Furthermore, the coordination number (Z) can be calculated as materials in zone 1 lost its volume and contact at the same time, it
the average total contact number as follows: would indicate a volume collapse during the transient sand production.
∑k j
The coordination number results in other zones followed a similar
j=1 Nt pattern with the behaviour in zone 1 but at a lesser extent. Contact
Z= (23)
k decrease occurred at nearly zero volume change for zones 2, 3 and 4. The
j
coordination number increased for all zones during the continuous sand
where Nt is the total contact number of the jth particle that includes both production starting from t3. The formation of new contacts would be the
the bonded and the unbonded contacts. results of new microstructural development, although different zones
Fig. 7 shows the variation of the total bonded contacts in the sample underwent different positive and negative volume changes during this
during the sand production process. The sanding onset conditions at t1 period.
and t4 for the LO and HO cases, respectively, occurred after the fluid Fig. 11 shows the fluid streamlines and the particle trajectories for
injection broke a certain number of bonds and a relatively constant bond light oil. The background of the fluid streamlines on the left shows the
number was maintained during the sand production process until the distribution of void fraction computed from the CFD simulation, which
end of the simulation. The fluid flow first created damage to the bond can be considered equivalent to the porosity of the solid phase. The
network and then removed the loosen particles from the sample. It is background on the right, on the other hand, represents the particles that
noted that the produced sand masses only account for less than 0.2% the are colour-coded according to the coordination number values of the
total sample, and the effect should hence be limited in these simulation particles. The particle trajectories were obtained from the resultant ve­
cases. locity vectors of all the particles in each fluid cell, the colours here
indicate the local average velocity magnitude of the particles. The void
fraction distribution results show more porous layers along the bound­
aries due to the flat wall surfaces. The upper part of the perforation

Fig. 7. Variation of the total bonded contacts in the system.

8
F. Khamitov et al. International Journal of Rock Mechanics and Mining Sciences 154 (2022) 105096

Fig. 8. Porosity behaviour during sand production with Light Oil.

Fig. 9. Porosity behaviour during sand production with Heavy Oil.

tunnel is associated with higher void fraction and parallel upward other hand, came from the loose zones following the more horizontal
streamlines were formed along the vertical sides of the tunnel in trajectories in the region just below the top surface. It should be ex­
Fig. 11a. As the streamlines reached the top outlet, the fluid leaved the pected that more severe continuous sand production would continue
sample, but parts of the flow were diverted to the sides of the outlet and only if the sand loss from the upper layer would trigger more disruption
formed two recirculating flow patterns that resemble each other on the to the internal microstructure of the sample in Fig. 11d.
opposite sides of the perforation. The coordination number values varied in a narrower range for the
The higher void fraction zone in the upper part of the perforation heavy oil in Fig. 12. The patterns showed an overall increase in contact
seems to grow larger in Fig. 11b, which also explains the contact loss in numbers as the primary trend for all zones. A small contact loss occurs
Fig. 10 as a zone of light shaded particles in the central part of the close to t4 at the onset of the continuous sand production with the most
sample. The particle velocity in this zone was higher than in the outer significant loss in zone 3 and the smallest loss in zone 1, which is
zones as the results of more significant particle movements. The higher opposite to the result of light oil in Fig. 10.
void fraction zone disappeared in Fig. 11c that indicates a volume The fluid streamlines and the particle velocity trajectories for heavy
collapse in the upper part of the perforation. New contacts were formed oil are given in Fig. 13. The void fraction and the flow patterns at t1 in
in the central zones and the particle velocities became more homoge­ Fig. 13a were similar to the results of light oil at t2 in Fig. 11b. A zone of
neous at t3. It should be noted that there was no significant change in the higher porosity and no contact was formed in the upper part of the
void fraction distribution in the outer zones in Fig. 11c. perforation. There were high velocity movements of both the fluid and
The streamlines became regular in the final state in Fig. 11d, there the particle in the central zone of the sample. The streamlines and
were two zones of high porosity near the outlet on the top surface. The particle trajectories already became regular at the onset of the contin­
particle trajectories showed nearly parallel columns of dark particles of uous sand production for heavy oil in Fig. 13b, which was not the case
higher coordination number inside the sample that were almost for light oil. The dark particles of higher coordination number values
perpendicular to the fluid streamlines. The produced particles, on the formed discontinuous clusters around but not perfectly aligned along the

9
F. Khamitov et al. International Journal of Rock Mechanics and Mining Sciences 154 (2022) 105096

Fig. 10. Variation of the coordination number for Light Oil.

Fig. 11a. Fluid streamlines (over void fraction background, left), and particle trajectories (over particle coordination number background, right) for Light Oil.
a) time = t1, vector scale factor = 0.015

Fig. 11b. Fluid streamlines (over void fraction background, left), and particle trajectories (over particle coordination number background, right) for Light Oil.
b) time = t2, vector scale factor = 0.015

10
F. Khamitov et al. International Journal of Rock Mechanics and Mining Sciences 154 (2022) 105096

Fig. 11c. Fluid streamlines (over void fraction background, left), and particle trajectories (over particle coordination number background, right) for Light Oil.
c) time = t3, vector scale factor = 0.1

Fig. 11d. Fluid streamlines (over void fraction background, left), and particle trajectories (over particle coordination number background, right) for Light Oil.
d) time = t5, vector scale factor = 0.5

Fig. 12. Variation of the coordination number for Heavy Oil.

trajectory curves. The zone of higher porosity in the upper perforation fraction decreased in the central zone was due to a volume expansion (or
finally disappeared in Fig. 13c, which would cause the contact increase void fraction increase) in the outer zones, which agrees with the largest
in zone 1 and contact loss in zone 3 between t4 and t5 in Fig. 12. The contact loss in zone 3 in Fig. 13.
void fraction distribution in Fig. 13c however showed that the void The dark particles in Fig. 13c formed continuous chains along the

11
F. Khamitov et al. International Journal of Rock Mechanics and Mining Sciences 154 (2022) 105096

Fig. 13a. Fluid streamlines (over void fraction background, left), and particle trajectories (over particle coordination number background, right) for Heavy Oil.
a) time = t1, vector scale factor = 0.015

Fig. 13b. Fluid streamlines (over void fraction background, left), and particle trajectories (over particle coordination number background, right) for Heavy Oil.
b) time = t4, vector scale factor = 0.1

Fig. 13c. Fluid streamlines (over void fraction background, left), and particle trajectories (over particle coordination number background, right) for Heavy Oil.
c) time = t5, vector scale factor =0.5

trajectory curves that were oriented toward the outlet. Although the 5. Conclusions
cumulative sand masses in the final state were similar for heavy oil and
light oil in Fig. 5, they were the results of different microstructures as Sand production is a common problem that is associated with the
shown in Fig. 11d and Fig. 13c. The continuous sand production for production of hydrocarbons from weak sandstone reservoirs in the oil-
heavy oil mobilized a much larger proportion of the sample’s particles gas industry. While the problem can be minimized by controlling the
and the particle movements were aligned much better with the fluid flow rate for light oil, the production of heavy oil causes more severe
flow. This would cause more severe sand production by heavy oil as it sand production that may need the installation of a sand control system
was observed in other researches in the literature.9,10,13 and thus increases the operational cost. We investigated the sanding

12
F. Khamitov et al. International Journal of Rock Mechanics and Mining Sciences 154 (2022) 105096

mechanisms for different reservoir fluids using a coupled simulation 9 Al-Awad MNJ, El-Sayed AAH, Desouky SEDM. Factors affecting sand production
from unconsolidated sandstone Saudi oil and gas reservoir. J King Saud Univ - Eng Sci.
with the Discrete Element Method for the solid sandstone particles and
1999;11(1):151–172. https://doi.org/10.1016/S1018-3639(18)30995-4.
the Computational Fluid Dynamics for the modelling of hydrocarbons 10 Wu B, Choi SK, Denke R, et al. A new and practical model for amount and rate of
flow. It was found that for a short simulation time both heavy and light sand production estimation. In: Day 3 Thu, March 24, 2016. OTC; 2016. https://doi.
oils produced similar sand masses, but they were due to very different org/10.4043/26508-MS.
11 Harris PC, Morgan RG, Heath SJ. Measurement of proppant transport of frac fluids.
microstructural changes inside the sandstone samples. There were more In: All Days. SPE; 2005:281–292. https://doi.org/10.2118/95287-MS.
bond breakages for light oil that facilitated more significant particle 12 Hughes K, Santos N, Arias RE, Nadezhdin SV. New viscoelastic surfactant fracturing
movement and porosity change within the damage zone around the fluids now compatible with CO 2 drasticaly improve gas production in rockies. Proc -
SPE Int Symp Form Damage Control. 2008;1:17–21. https://doi.org/10.2118/111431-
central perforation tunnel. Sand production occurred earlier for light oil ms.
in two stages of a first transient sand production, and followed by a 13 Ahmad FA, Miskimins JL. Proppant transport and behavior in horizontal wellbores
second continuous sand production. In the final time step, sands were using low viscosity fluids. In: Day 2 Wed, February 06, 2019. SPE; 2019. https://doi.
org/10.2118/194379-MS.
being produced from two loose zones in an upper layer near the outlet 14 Zhu HP, Zhou ZY, Yang RY, Yu AB. Discrete particle simulation of particulate
and hence should be limited in scope. There was only a continuous sand systems: theoretical developments. Chem Eng Sci. 2007;62(13):3378–3396. https://
production for heavy oil that occurred later than for light oil. Bond doi.org/10.1016/j.ces.2006.12.089.
15 Tsuji Y, Kawaguchi T, Tanaka T. Discrete particle simulation of two-dimensional
breakage, particle movement and porosity change all happened to a fluidized bed. Powder Technol. 1993;77(1):79–87. https://doi.org/10.1016/0032-
lesser extent. The particle velocity trajectories became regular more 5910(93)85010-7.
quickly that followed the fluid streamlines toward the outlet point. The 16 O’Connor RM, Torczynski JR, Preece DS, Klosek JT, Williams JR. Discrete element
modeling of sand production. Int J Rock Mech Min Sci. 1997;34(3-4):e15. https://doi.
heavy oil outflow mobilized a larger proportion of the sample’s particles
org/10.1016/S1365-1609(97)00198-6, 231.e1-231.
and would lead to more severe sand production. The research on the 17 Narsilio GA, Buzzi O, Fityus S, Yun TS, Smith DW. Upscaling of Navier–Stokes
impact of different reservoir fluid flows on the sand production is being equations in porous media: theoretical, numerical and experimental approach.
extended where we are investigating the effect of multiphase flow of Comput Geotech. 2009;36(7):1200–1206. https://doi.org/10.1016/j.
compgeo.2009.05.006.
different reservoir fluids on sand production, the results of which will be 18 Irmay S. On the theoretical derivation of Darcy and Forchheimer formulas. Eos, Trans
reported in the coming publication. Am Geophys Union. 1958;39(4):702–707. https://doi.org/10.1029/
TR039i004p00702.
19 Khamitov F, Minh NH, Zhao Y. Coupled CFD–DEM numerical modelling of
CRediT authorship contribution statement perforation damage and sand production in weak sandstone formation. Geomech
Energy Environ. 2021;28, 100255. https://doi.org/10.1016/j.gete.2021.100255.
Furkhat Khamitov: Methodology, Software, Formal analysis, 20 Goniva C, Kloss C, Deen NG, Kuipers JAM, Pirker S. Influence of rolling friction on
single spout fluidized bed simulation. Particuology. 2012;10(5):582–591. https://doi.
Writing – original draft, Writing – review & editing, Validation. Nguyen org/10.1016/j.partic.2012.05.002.
Hop Minh: Conceptualization, Supervision, Writing – review & editing. 21 Kozhagulova A, Minh NH, Zhao Y, Fok SC. Experimental and analytical investigation
Yong Zhao: Conceptualization, Supervision, Writing – review & editing, of sand production in weak formations for multiple well shut-ins. J Petrol Sci Eng.
2020;195, 107628. https://doi.org/10.1016/j.petrol.2020.107628.
Funding acquisition, Project administration. 22 Shabdirova A, Minh NH, Zhao Y. A sand production prediction model for weak
sandstone reservoir in Kazakhstan. J Rock Mech Geotech Eng. 2019;11(4):760–769.
Declaration of competing interest https://doi.org/10.1016/j.jrmge.2018.12.015.
23 Rakhimzhanova AK, Khamitov FA, Minh NH, Thornton C. 3D DEM simulations of
triaxial compression tests of cemented sandstone. In: Proc IS Atlanta 2018 Symp
The authors declare that they have no known competing financial Geomech from Micro to Macro Res Pract. 2018. Published online.
interests or personal relationships that could have appeared to influence 24 Climent N, Arroyo M, O’Sullivan C, Gens A. Sand production simulation coupling
the work reported in this paper. DEM with CFD. Eur J Environ Civ Eng. 2014;18(9):983–1008. https://doi.org/
10.1080/19648189.2014.920280.
25 Crowe CT, Schwarzkopf JD, Sommerfeld M, Tsuji Y. Multiphase Flows with Droplets
Acknowledgement and Particles. second ed. CRC Press; 2011. https://doi.org/10.1201/b11103.
26 Li Y, Zhang J, Fan LS. Numerical simulation of gas–liquid–solid fluidization systems
using a combined CFD-VOF-DPM method: bubble wake behavior. Chem Eng Sci.
This research was sponsored by Nazarbayev University research 1999;54(21):5101–5107. https://doi.org/10.1016/S0009-2509(99)00263-8.
grant No OPCRP2022006. First author also would like to acknowledge 27 Zhou ZY, Kuang SB, Chu KW, Yu AB. Discrete particle simulation of particle–fluid
the research grant, No AP08052762, from the Ministry of Education and flow: model formulations and their applicability. J Fluid Mech. 2010;661:482–510.
https://doi.org/10.1017/S002211201000306X.
Science of the Republic of Kazakhstan. 28 Di Felice R. The voidage function for fluid-particle interaction systems. Int J
Multiphas Flow. 1994;20(1):153–159. https://doi.org/10.1016/0301-9322(94)
References 90011-6.
29 Ergun S. Fluid flow through packed columns. Chem enfineering Prog. Published online.
1952:89–94.
1 Kazakhstan Upstream oil and gas technology and R&d Roadmap. Inst Manuf Univ
30 Gidaspow D, Bezburuah R, Ding J. Hydrodynamics of circulating fluidized beds:
Cambridge. Published online 2013. https://www.ifm.eng.cam.ac.uk/uploads/Roa
kinetic theory approach. 7th Fluid Conf. Published online; 1992:75–82. http://www.
dmapping/Kaz_RM_book_English.pdf.
osti.gov/scitech/servlets/purl/5896246.
2 Dusseault MB, El-Sayed S. CHOP - cold heavy oil production. In: IOR 1999 - 10th
31 Koch DL, Hill RJ. Inertial effects in suspension and porous—media flow. Annu Rev
European Symposium on Improved Oil Recovery. European Association of Geoscientists
Fluid Mech. 2001;33(1):619–647. https://doi.org/10.1146/annurev.fluid.33.1.619.
& Engineers; 1999. https://doi.org/10.3997/2214-4609.201406351.
32 Kafui KD, Thornton C, Adams MJ. Discrete particle-continuum fluid modelling of
3 Nie Z, Zheng X, Xie C. New cementing technologies successfully solved the problems
gas–solid fluidised beds. Chem Eng Sci. 2002;57(13):2395–2410. https://doi.org/
in shallow gas, low temperature and easy leakage formations of North Buzachi
10.1016/S0009-2509(02)00140-9.
oilfield. In: International Oil and Gas Conference and Exhibition in China. Society of
33 Agrawal V, Shinde Y, Shah MT, Utikar RP, Pareek VK, Joshi JB. Effect of drag models
Petroleum Engineers; 2010. https://doi.org/10.2118/131810-MS.
on CFD–DEM predictions of bubbling fluidized beds with Geldart D particles. Adv
4 Acock A, ORourk A, Shirmboh D, Alexander J, Andersen G, López-de-Cárdenas J.
Powder Technol. 2018;29(11):2658–2669. https://doi.org/10.1016/j.
Practical approaches to sand management. Oil F Rev. 2004;16:10–27. http://www.
apt.2018.07.014.
slb.com/~/media/Files/resources/oilfield_review/ors04/spr04/02_sand_manageme
34 Geldart D. Elsevier sequoia SA, lausanne-printed in the netheriands types of gas
nt.pdf.
fhidization. Powder Technol. 1973;7:285–292.
5 Tronvoll J, Dusseault MB, Sanfilippo F, Santarelli FJ. The tools of sand management.
35 Marchelli F, Hou Q, Bosio B, Arato E, Yu A. Comparison of different drag models in
In: All Days. SPE; 2001. https://doi.org/10.2118/71673-MS.
CFD-DEM simulations of spouted beds. Powder Technol. 2020;360:1253–1270.
6 Tremblay B, Sedgwick G, Forshner K. Modelling of sand production from wells on
https://doi.org/10.1016/j.powtec.2019.10.058.
primary recovery. In: Annual Technical Meeting. Petroleum Society of Canada; 1996.
36 Rong LW, Dong KJ, Yu AB. Lattice-Boltzmann simulation of fluid flow through
https://doi.org/10.2118/96-26.
packed beds of uniform spheres: effect of porosity. Chem Eng Sci. 2013;99:44–58.
7 Geilikman MB, Dusseault MB. Fluid rate enhancement from massive sand production
https://doi.org/10.1016/j.ces.2013.05.036.
in heavy-oil reservoirs. J Petrol Sci Eng. 1997;17(1-2):5–18. https://doi.org/10.1016/
37 Cundall PA, Strack ODL. A discrete numerical model for granular assemblies.
S0920-4105(96)00052-6.
Geotechnique. 1979;29(1):47–65. https://doi.org/10.1680/geot.1979.29.1.47.
8 Veeken CAM, Davies DR, Kenter CJ, Kooijman AP. Sand production prediction
review. Developing an integrated approach. Proc - SPE Annu Tech Conf Exhib. 1991;Pi
(pt 1):335–346. https://doi.org/10.2523/22792-ms.

13
F. Khamitov et al. International Journal of Rock Mechanics and Mining Sciences 154 (2022) 105096

38 Berger R, Kloss C, Kohlmeyer A, Pirker S. Hybrid parallelization of the LIGGGHTS 55 Weaver D, Miskovic S. Analysis of coupled CFD-DEM simulations in dense particle-
open-source DEM code. Powder Technol. 2015;278:234–247. https://doi.org/ laden turbulent jet flow, 2. In: Fluid Mechanics; Multiphase Flows. vol. 2. American
10.1016/j.powtec.2015.03.019. Society of Mechanical Engineers; 2020. https://doi.org/10.1115/FEDSM2020-
39 Hertz H. Über die Berührung fester elastischer Körper. J für die Reine Angewandte 20274.
Math (Crelle’s J). 1881;171:156–171. 56 Chilamkurti YN. Towards Understanding the Heat Transfer Behavior of Dense Granular
40 Johnson KL. Surface energy and the contact of elastic solids. Proc R Soc London A Media; 2019. Published online https://www.proquest.com/dissertations-theses/
Math Phys Sci. 1971;324(1558):301–313. https://doi.org/10.1098/rspa.1971.0141. towards-understanding-heat-transfer-behavior/docview/2292188637/se-2?acco
41 Johnson KL. In: Koiter W, ed. Adhesion at the Contact of Solids. New York: North untid=134066.
Holland; 1976:133–143. Proc XIV Int Congr Theor Appl Mech. 57 Shabdirova AD, Bissekenova Z, Minh NH, Kim JR. Sample preparation method of
42 Kafui DK, Johnson S, Thornton C, Seville JPK. Parallelization of a clay-rich sandstone analogue of sandstone reservoirs in Kazakhstan. 50th US Rock
Lagrangian–Eulerian DEM/CFD code for application to fluidized beds. Powder Mech/Geomech Symp. 2016;2:904–910, 2016.
Technol. 2011;207(1-3):270–278. https://doi.org/10.1016/j.powtec.2010.11.008. 58 DCS Computing GmbH. CFDEM®coupling Documentation. https://www.cfdem.
43 Lian G, Thornton C, Icdfui D. Trubal A 3-D Computer Program for Modelling Particle com/media/CFDEM/docu/CFDEMcoupling_Manual.html#.
Assemblies. 1998. Published online. 59 Goniva C, Blais B, Radl S, Kloss C. Open source cfd-dem modelling for particle-based
44 Mindlin RD, Deresiewicz H. Thickness-shear and flexural vibrations of a circular disk. processes. 11th Int Conf CFD Miner Process Ind. 2015;December:1–7. https://www.
J Appl Phys. 1954;25(10):1329–1332. https://doi.org/10.1063/1.1721554. cfd.com.au/cfd_conf15/PDFs/140GON.pdf.
45 Thornton C. Granular Dynamics, Contact Mechanics and Particle System Simulations. 60 Issa R. Solution of the implicitly discretised fluid flow equations by operator-
vol. 24. Springer International Publishing; 2015. https://doi.org/10.1007/978-3- splitting. J Comput Phys. 1986;62(1):40–65. https://doi.org/10.1016/0021-9991(86)
319-18711-2. 90099-9.
46 Yan B, Regueiro RA. Superlinear speedup phenomenon in parallel 3D Discrete 61 OpenCFD OpenFoam V2106. Programmer’s guide.; 2021. https://www.openfoam.co
Element Method (DEM) simulations of complex-shaped particles. Parallel Comput. m/documentation/overview.
2018;75:61–87. https://doi.org/10.1016/j.parco.2018.03.007. 62 Chen G, Xiong Q, Morris PJ, Paterson EG, Sergeev A, Wang YC. OpenFOAM for
47 Cheng J, Dou Y, Zhang N, Li Z, Wang Z. A new method for predicting erosion damage computational fluid dynamics. Not Am Math Soc. 2014;61(4):354. https://doi.org/
of suddenly contracted pipe impacted by particle cluster via CFD-DEM. Materials. 10.1090/noti1095.
2018;11(10):1858. https://doi.org/10.3390/ma11101858. 63 Kloss C, Goniva C. Liggghts - open source discrete element simulations of granular
48 Norouzi HR, Zarghami R, Sotudeh-Gharebagh R, Mostoufi N. Coupled CFD-DEM materials based on lammps. In: Supplemental Proceedings. John Wiley & Sons, Inc.;
Modeling. John Wiley & Sons, Ltd; 2016. https://doi.org/10.1002/9781119005315. 2011:781–788. https://doi.org/10.1002/9781118062142.ch94.
49 Goniva C, Kloss C, Hager A, Pirker S. An Open Source CFD-DEM Perspective. 2010. 64 DCS Computing GmbH. LIGGGHTS(R)-PUBLIC documentation, Version 3.X.
February 2015. https://www.cfdem.com/media/DEM/docu/Manual.html.
50 Peng Z, Doroodchi E, Luo C, Moghtaderi B. Influence of void fraction calculation on 65 Bealessio BA, Blánquez Alonso NA, Mendes NJ, Sande AV, Hascakir B. A review of
fidelity of CFD-DEM simulation of gas-solid bubbling fluidized beds. AIChE J. 2014; enhanced oil recovery (EOR) methods applied in Kazakhstan. Petroleum. 2020:1–9.
60(6):2000–2018. https://doi.org/10.1002/aic.14421. https://doi.org/10.1016/j.petlm.2020.03.003. November 2019.
51 Clarke DA, Sederman AJ, Gladden LF, Holland DJ. Investigation of void fraction 66 ENGYS Ltd. HELYX-OS.Open-source Graphical User Interface developed by ENGYS to
schemes for use with CFD-DEM simulations of fluidized beds. Ind Eng Chem Res. operate and control natively the standard CFD utilities and solvers provided with the
2018;57(8):3002–3013. https://doi.org/10.1021/acs.iecr.7b04638. OpenFOAM libraries. http://engys.github.io/HELYX-OS/.
52 Smuts EM. A Methodology for Coupled CFD-DEM Modeling of Particulate Suspension 67 Laney CB. The CFL condition. In: Computational Gasdynamics. Cambridge University
Rheology; 2015. Published online https://open.uct.ac.za/handle/11427/16782. Press; 1998:214–221. https://doi.org/10.1017/CBO9780511605604.016.
53 Boyce CM, Holland DJ, Scott SA, Dennis JS. Limitations on fluid grid sizing for using 68 Courant R, Friedrichs K, Lewy H. Über die partiellen Differenzengleichungen der
volume-averaged fluid equations in discrete element models of fluidized beds. Ind mathematischen Physik. In: Kurt Otto Friedrichs. Birkhäuser Boston; 1986:53–95.
Eng Chem Res. 2015;54(43):10684–10697. https://doi.org/10.1021/acs. https://doi.org/10.1007/978-1-4612-5385-3_7.
iecr.5b03186. 69 Coker AK. Mechanical separations. In: Ludwig’s Applied Process Design for Chemical
54 DCS Computing GmbH. CFDEM®coupling Documentation. Divided Void Fraction. and Petrochemical Plants. Elsevier; 2007:371–443. https://doi.org/10.1016/B978-
Published; 2016. https://www.cfdem.com/media/CFDEM/docu/voidFractionMod 075067766-0/50013-0.
el_dividedVoidFraction.html.

14