Coupled Modeling of Computational Fluid Dynamics and Granular Mechanics of Sand Production in Multiple Fluid Flow

Waterflooding, which is a technique to enhance oil recovery, involves injecting water to displace oil toward production wells. This process can create an irregular pressure front, leading to water fingering, which results in less efficient extraction of hydrocarbons, as shown in Fig. 1a  (Udy et al. 2017; Kargozarfard et al. 2019). As the production rate increases, the reservoir pressure around the wellbore decreases, causing the oil-water interface to move toward the perforation interval. This movement can result in the interface forming a cone shape, commonly referred to as water coning (Fig. 1b ; Muskat and Wycokoff 1935). A critical flow rate exists, which acts as a key threshold for the stability of the water cone. When this critical flow rate is exceeded, the water cone transitions from a stable to an unstable state, leading to water breakthrough.

Water breakthrough, which occurs when water reaches the production well, indicates the onset of a multiphase flow. This often leads to a decline in both production rate and ultimate hydrocarbon recovery while concurrently increasing production costs due to the need for water-handling measures (Bailey et al. 2000). The phenomena of water breakthrough are linked to the critical flow rate and breakthrough time, heavily influenced by both fluid and rock properties, as documented in the literature (Karami et al. 2014; Yue et al. 2019; Bai and Tahmasebi 2021).

The dynamics of waterflooding and water breakthrough can also affect sand production. Wells may exhibit a temporary spike in sand production concurrent with water breakthrough events, suggesting that the disruption of capillary cohesion during water breakthrough could increase the sand production rate (Papamichos et al. 2010). However, this spike tends to diminish once the breakthrough has occurred. Despite thorough experiments conducted on various materials, including weak outcrop sandstone and field core samples, a consistent correlation between the onset of sand production and water breakthrough has remained elusive (Skjaerstein et al. 1997).

Moreover, waterflooding is frequently accompanied by the fingering effect, especially when there is a significant difference between the viscosities of the involved phases. Fluid flow occurs in isolated narrow channels or pathways rather than dispersing uniformly throughout the reservoir, as depicted in Fig. 1a . While fluid fingering and flow channeling have been successfully modeled in the context of viscous oils (Brantson et al. 2018; Pinilla et al. 2021), their impact on sand production has not been thoroughly examined. Recent coupled simulations of CFD-DEM have provided a more nuanced understanding of multiphase flow in the presence of solids (Alaita 2021) and the migration of fines within porous media (Zhu et al. 2022).

Despite previous research, the relationship between sand production, water breakthrough, water coning, and the fingering effect in multiphase flow is not well understood. The application of waterflooding in the heavy oil reservoirs located in southwest Kazakhstan provides an example of these dynamics in action (Mustafayev 2018). These reservoirs, situated within weak and shallow formations, are notably susceptible to sand production (Soroush et al. 2021; Bealessio et al. 2021). Observational data from these wells indicate a pattern of transient sand production behavior, characterized by an initial spike followed by a rapid decline and eventual cessation of sand production. Interestingly, the relationship between sand production and water breakthrough has shown variability in the timing of sand spikes across different wells, underscoring the complexity of these interactions (Veeken et al. 1991; Wu et al. 2006; Collins et al. 2008; Tananykhin et al. 2021; Khamitov et al. 2022).

The microscopic mechanisms underpinning sand production under different flow conditions of a single fluid have been investigated by Khamitov et al. (2021, 2022). High-viscosity fluids can produce more sand due to their enhanced transportability. This study further developed these simulations to consider the multiple fluid flow of water and heavy oil through porous sandstone. Through multiphase and multiple fluid-coupled CFD-DEM simulations, efforts have been made to bridge the gap between macroscopic sanding behavior and the microscopic evolution of bond breakage, void fraction, oil-cut change, and the trajectories of fluid and particles across the sample. Notably, water coning and fingering were observed in this study, aligning with experimental results that suggest sand production can occur in multiple stages, dependent on the flow dynamics, fluid-solid interactions, and the inherent properties of the different phases within the sandstone specimen.

The CFD-DEM model is based on the Euler-Lagrange scheme, where a continuous fluid flow is coupled with discrete particles by exchanging momentum, energy, and mass terms. This approach was first proposed by Tsuji et al. (1993) and further developed by Zhou et al. (2010). The model considers particles at the microscale and fluid at the mesoscale.

DEM has been extensively used to investigate the mechanics of granular flow, as pioneered by Cundall and Strack (1980). Particles within a computational domain are explicitly tracked by solving governing particle motion equations. This provides detailed information on the interparticle contact forces and particle trajectories. Different materials can be simulated using DEM by specifying the physical properties of particles and contact model parameters, which allows for various configurations of linear and nonlinear elastoplastic behavior (e.g., Mindlin and Deresiewicz 1954; Johnson et al. 1971; Thornton 2015). Cemented rock material can also be modeled as an assembly of fully or partially bonded particles. The bonding force was included in calculating contact forces and the particles’ translational and angular accelerations. In this study, the cemented contacts were modeled using the modified Johnson-Kendall-Roberts (JKR) contact model (Khamitov et al. 2022: 2021).

In CFD-DEM simulation, the particle phase is coupled with a fluid phase (Zhu et al. 2007). The numerical model accounts for particle-particle and fluid-particle interaction forces in Newton’s second law that governs particle motion. For any particle p (Goniva et al. 2012)

where $mp$ denotes the particle mass, the linear velocity is $u→p$⁠, and $rp,c$ refers to the radius of the particle. The normal and tangential contact forces between particles are $f→p,n$ and $f→p,t$⁠, respectively, and $f→p,f$ represents the drag force exerted on the particle from the fluid phase. In addition, $f→p,Pres$ denotes the pressure force and $f→p,vis$ represents the viscous force. The $ω→p$ is angular velocities and $Ip$ represents the moment of inertia.

The pressure force $f→p,Pres$ can be defined as:

where $∇P$ is the fluid pressure gradient and $Vp$ is the volume of particles.

The viscous force $f→p,vis$ can be calculated as follows:

where $τ̿$ is the fluid stress tensor (see Eq. 6).

Eq. 1 describes the translational motion, while Eq. 2 is used to describe the rotational motion of particle $p$⁠. In the absence of the liquid phase, the fluid-related force terms (⁠$f→p,f,f→p,Pres,f→p,vis$⁠) are equal to zero. The contact forces between particles, $f→p,n$ and $f→p,t$⁠, were calculated using the modified JKR contact model (Khamitov et al. 2021, 2022).

The CFD-DEM coupling in this study was based on Model A, as described by Zhou et al. (2010) and Norouzi et al. (2016). The modified Navier-Stokes equations and the interface tracking volume-of-fluid (VOF) method (Nichols and Hirt 1975; Hirt and Nichols 1981; Enwald et al. 1996) were used to describe the governing equations of two-phase fluid flow (Phase 1 and Phase 2) in the presence of a third particulate phase (sand particles). The CFD-DEM-VOF method has been implemented successfully in multiphase studies and validated by well-known benchmark tests (Washino et al. 2021; Çalışkan and Mišković 2024). The unresolved CFD-DEM model (fluid phase is not resolved at the particle scale) with VOF can be formulated as follows:

where the volume fraction occupied by the fluid is represented by $εf$⁠, while $αph1$ and $αph2$ refer to the volume fraction of Phase 1 and Phase 2, respectively. The compression velocity at the interface between the phases is denoted as $u→c$⁠, where $u→c=u→f,ph1−u→f,ph2$⁠. Additionally, other variables include $u→f$⁠, which represents fluid velocity, $ρf$ for fluid density, $P$ for pressure, $F→pfA$ for volumetric particle-fluid interaction force, $g→$ for gravity vector, and $t$ for time. The stress tensor $τ̿$ is obtained for fluid with viscosity $μf$as

In contrast to the Euler-Euler two-fluid approach, the VOF method assumes the absence of interaction forces between fluid phases and solves only a single momentum and continuity equation (Enwald et al. 1996) through employing the fluid-mixture properties that are weighted by volume:

The volumetric particle-fluid interaction force $F→pfA$ is determined as

where the term $∆V$ represents the volume of the fluid cell, while $n$ denotes the number of particles in the computational cell. $f→p,d$ represents the drag force. Additionally, $f→p′′$ refers to the sum of particle-fluid interaction forces acting on particles other than the drag force. These forces are related to pressure gradient and viscous stress and are commonly regarded as the dominant forces in particle-fluid flow.

The interactions between the fluid and the particle phases were simulated by a four-way coupling method (Norouzi et al. 2016). Four-way coupling in CFD-DEM simulations aims to capture the intricate interplay between fluid and particle phases in a system, considering both direct and indirect interactions among particles and between particles and the surrounding fluid. The motion of particles affects the motion of fluid and vice versa. Three forces were considered in the coupling information: pressure force, viscous force, and the Di Felice drag force (Di Felice 1994; Zhou et al. 2010). The Di Felice drag formulation assumes spherical particles, which do not interact with each other, and neglects the effects of particle shape and size variation on drag force. Taking into account the continuity and precision of the different drag force calculation models, the Di Felice model provided more accurate calculations for dense particle systems (Agrawal et al. 2018; Marchelli et al. 2020; Wang et al. 2021). Therefore, this drag model has been chosen for the current study. Further information can be found in previous papers (Khamitov et al. 2021, 2022).

The coupled simulation was conducted in multiple stages. First, a sandstone sample was prepared with particle pluviation, compressional diagenesis, and cementation. The sandstone sample was perforated to form a central vertical channel crossing the full sample’s height and connected to an outlet point at the center of the top surface. The fluid phase was added to the solid phase when the CFD mesh was formed with the same dimensions as the perforated sandstone sample. Fluid saturation was simulated, which was followed by radial fluid injection from the sidewall of the cylindrical sample. The injected fluid flow penetrated the granular sample and was directed toward the top outlet point. Fluid flow was accelerated by the perforation channel. Single fluid flow was simulated using the same saturating and injecting fluid, whereas multiphase fluid flow used different fluids for these conditions. Fluid injecting velocity was varied to simulate different production conditions in the well.

The DEM model used for the dry sample (no liquid phase) was similar to previous works (Khamitov et al. 2021, 2022). Different size particles were randomly generated and allowed to fall under gravity inside a cylindrical wall. The particle-size distribution of the numerical sample was selected to match that of a real sandstone from the Ustuyrt-Buzachi Sedimentary Basin, as reported by Shabdirova et al. (2016). The total number of particles was 100,010, which corresponded to a sample mass of 1.58 g. The micromechanical properties of sand material for DEM were sand particle Young’s modulus E = 70 GPa, Poisson’s ratio and friction coefficient of about 0.3, particle density $ρp$ = 2500 kg/m³, and interface energy $γ$ = 40 J/m². The top wall was moved downward to compress the particles until a vertical stress of 5 MPa was attained, resulting in an intact sample with a diameter (d) of 15.12 mm and a height (h) of 5.58 mm. The axial cross section of the compressed sample is shown in Fig. 2 .

Particles in the sandstone sample could develop multiple contacts with their surrounding neighbors. In its intact state, all bonded contacts were formed at once using the modified JKR model (Khamitov et al. 2021, 2022). This model agrees well with the experimental study with sandstone from the Ustyurt-Buzachi Sedimentary Basin. The modified JKR model has two assumptions. First, the bond contact breaks in the JKR model occur at the maximum negative value of tensile force, which strictly depends on interface energy and relative curvature of the particles’ contact (Thornton 2015; Johnson et al. 1971). The second assumption is that the bonds are formed only once at the cementation stage. However, during subsequent perforation and sand production stages, mechanical forces can break these bonds and form new unbonded connections. These unbonded connections are modeled using the Hertz (Hertz 1881) contact model. An algorithm was implemented to count the number of (unbonded) Hertz contacts and of the (bonded) modified JKR contacts for each particle.

Perforation was simulated by driving a penetrometer with a diameter of 1.4 mm vertically into the sample and then withdrawing it. Particle bonds were destroyed in this process, and this resulted in a damaged region of unbonded particles in the middle of the sample as shown in Fig. 2 . Further information can be found in similar simulations (Khamitov et al. 2021).

Sand production simulation was conducted in two stages—fluid saturation and injection. Rather than simulating a gradual saturation process, we instantly imposed the CFD grid on the DEM grid. This decision was made to optimize simulation time, as simulating the saturating process using low fluid velocities would be very time-consuming. Particle movement was assumed negligible for the application of the instant saturation. This is comparable with the experimental procedure where the sample was slowly flushed with the saturating fluid from the bottom inlets at a very low injection pressure to avoid disturbance to the bonded granular structure of the sandstone specimen. During the second stage of fluid injection and sand production, the particles moved toward the central perforation under the influence of the injecting fluid and were deleted once they exited the outlet element and crossed a predetermined level above the top of the sample. Numerical simulations of single and multiphase fluids were conducted to examine the different sand production mechanisms. The fluid was injected radially at different velocities (ranging from 10⁻⁴ to 10⁻² m/s) toward the central vertical perforation (Fig. 2) . The production of sand particles was highly dependent on the solid stress and fluid pressure conditions.

The coupled CFD-DEM formulation in Eq. 5 was implemented as Model A in the commercial CFDEM®coupling program (CFDEMresearch GmbH, no date). The CFDEM program is a highly parallelized computational platform that combines the Aspherix® program for the DEM modeling and the open-source OpenFOAM-8 program for the CFD modeling (DCS Computing GmbH 2021). The numerical simulations in this study were executed on a high-performance computer system, which was equipped with 44 multithreading processors Intel Xeon® E5-2699 Version 4 with a clock speed of 2.20 GHz and 500 GB of memory. A typical simulation could take about 11 days to complete, even with the parallel computation that utilized as many as eight threads using the message passing interface protocol.

The CFD mesh grid was developed using the HELYS-OS graphical user interface (see https://engys.github.io/helyx-os/). The computational fluid domain in Fig. 2  has 2,344 cells that formed the same dimensions of the DEM sample. The CFD mesh cell $∆xcfd$ equal to 7.96×10⁻⁴ was used to achieve the mesh-independent solution, where the accuracy of the simulation is still acceptable, and computation speed is maximized (Khamitov et al. 2022). The ratio between the size of the CFD mesh and the maximum particle diameter is equal to 2.24. The divided void fraction model was used for fluid fraction calculation in CFD cells (DCS Computing GmbH 2021). The top and bottom walls were impermeable with a no-slip condition. A central mesh element on the top surface served as the fluid outlet. The injecting and saturating fluids were treated as incompressible. Table 1  shows the properties of two different fluids (namely, heavy oil and water) used in the simulation, and the properties of the motor oil used in the experimental study. The outlet pressure was maintained at a constant value of 1 MPa throughout the simulation.

The DEM and CFD simulations used discretization timesteps of 2×10⁻⁸ and 2×10⁻⁷ seconds, respectively. To synchronize the two simulations, the CFD simulation was run every 10 timesteps of the DEM simulation, which allowed for accurate and efficient two-way coupling data exchange. The DEM simulation provides the locations of each particle for the calculating of the fluid volume fraction $εf$ and the particle-fluid interaction force $F→pfA$ in each CFD cell (Eq. 5). The updated fluid pressure and velocity affect particle gradient pressure, viscous, and drag forces.

Three numerical simulation cases were conducted to study the effect of multiphase fluid flow on sand production. Different combinations of water and heavy oil were used as the saturating and injecting fluids, respectively, the details of which are presented in Table 2 . Here, SW and SO are the single water and single oil cases, whereas MP is the multiphase case. The same injecting velocity was used for all cases.

The cumulative sand production and total bond number are plotted in Fig. 3  for all cases. Sand production began with a rapid loss of contact bonds, which reached a steady state when there was no longer bond breakage. As the cumulative sand becomes constant, sand production ceases, and this defines a transient sand production phenomenon. The multiphase flow of water injection into an initially saturated sample by heavy oil (MP) reached the maximum cumulative sand earlier than the cases of the single fluid flow for both heavy oil (SO) and water (SW).

The SO case produced the largest cumulative sand mass in two stages that were separated by a short period of no sand production as marked by a constant sand mass. The first sand production stage occurred with a much higher bond loss, whereas in the second stage, around 30% sand mass of the first stage was produced with only a few additional bond breakage and sanding did not completely stop but rather continued at a very small rate at the end of the simulation. The SW case showed a more gradual transient sand production, where the maximum cumulative sand mass was achieved at 0.3 seconds as compared with 0.1 seconds of the multiphase flow case, although similar sand masses were produced for both cases at the end.

The drag force exerted on the particles from the fluid phases is proportional to the velocity difference between the particle and fluid phases. Furthermore, it increases with higher density and viscosity. At-rest particles accelerated by the injected fluid would lead to higher drag force in the beginning, which then decreased with time as shown in Fig. 4 . The highest viscosity of the SO case would cause the highest drag force among all cases and vice versa for the SW case with the smallest viscosity. There was no measurement of drag force in the experiment, but the motor oil would cause even higher drag force for the same given conditions. The SO case produced the largest cumulative sand due to its higher transportability of heavy oil (Khamitov et al. 2022).

The fluid streamlines and void fraction are plotted in Fig. 5 . The radial flow injection created a zone of maximum void fraction (i.e., empty space) near the boundary as the hydrodynamic forces moved all particles toward the central axis. Sanding stopped at the timestep t1 for the MP case, when an annular zone of highly porous materials can be seen on the top surface. High fluid flow velocities occurred near the top surface along the left branch of this zone toward the outlet point. Sand production created a nonsymmetric structure of localized damage inside the specimen under multiphase flow condition that is similar to a fingering effect. A small empty zone (in red color) can be observed at later timesteps t2 and t3, which could be due to the local particles moving downward as the particles underneath them moved upward diagonally, but this did not contribute to any additional sand production.

Without the fingering effect, the velocity fields of the single fluid flow cases appeared to be more symmetric, with vectors being oriented toward the central perforation and maximum velocities located near the outlet point. Sand production continued at t2 for the SW case as additional particles were produced, leading to the formation of a nearly symmetric complete empty zone underneath the top surface at t3. The second stage of sand production for the SO case was associated with the creation of a new vertical empty channel around the central perforation that, when combined with the first annular damage zone, created a larger damage zone in the approximated shape of a hemisphere. The lack stage second-stage sand production led to a shallower hemisphere for the SW case as compared with the SO case. For the MP case, the localized damage or fingering restricted the development of a full symmetric damage hemisphere as of the single fluid cases. Note that the damage zone on the top surface of the sample was also observed in the experiment.

Drawing velocity streamlines for all particles in one figure is challenging. The average particle velocities in CFD cells were calculated from DEM data and plotted with the void fraction in Fig. 6  to increase visibility. High particle velocities appeared near the outlet point and only when there was an ongoing sand production. Particle movements in other regions of the specimen were small, the vectors of which were oriented more randomly. Small particle movements were also observed for the stages of no sand production that would lead to new contact formation and stabilization. There was no sand production at t3 for the SW and the MP cases, but sand production continued for the SO case at t3, which agrees with the other results of the cumulative sand production.

To study the effect of water breakthrough, three cases of high (MP-H), medium (MP-M), and low (MP-L) injecting velocity were investigated on the multiphase flow (MP) samples as shown in Table 3 . Note that the medium injecting velocity was previously used for the cases in Table 2 .

Fig. 7  shows the exit oil-cut value for the multiphase simulations, which was calculated as the fraction of oil volume in the CFD cell at the outlet:

where $Voil$ is the volume of oil and $Vtotal$ is the total volume of the CFD cell.

As samples were fully saturated with oil, the initial oil-cut value was equal to unity, which decreased over time if the injecting water reached the outlet. The timestep at which the value drops below unity is considered as the water breakthrough occurring as the injected water reaches the outlet. The high injecting velocity could nearly bring the oil-cut value down to zero as oil was replaced by water at the outlet. For the medium velocity, about 70% of the oil was replaced by water, while for low velocity, it was still producing only oil at the end of the simulation. Water breakthrough occurred earlier for higher velocity. Furthermore, the exit oil-cut value decreased more quickly with increasing injecting velocities.

There was no correlation found between water breakthrough and the development of sand production. Sand production continued for MP-H after water breakthrough, at almost the same rate as shown in a linear cumulative sand mass curve, while for MP-M, water breakthrough occurred close to the end of sand production. Sand production, however, ceased for both medium and high velocities, which was not the case for the low velocity.

The oil-cut distribution over the whole specimen is overlaid with the fluid streamlines in Fig. 8 . The injected water displaced oil radially, which formed parallel vertical interfaces between zones of different oil-cut values for the MP-L samples. Fluid streamlines created a symmetric flow pattern toward the top central outlet, and there was no water breakthrough for this case. For higher velocity cases, the vertical interfaces were gradually bent toward the outlet, which formed two symmetric half water cone shapes that were connected to each other at the outlet point and signified a water breakthrough. Injecting fluid streamlines were diverted to follow the water cone boundary. The remaining oil volume was surrounded by these water cones and maintained in the low to middle regions of the sample. The remaining oil volume was smaller with higher injecting velocity as more oil volume was displaced by water. These oil-trapped pockets would decrease oil production efficiency once there was water breakthrough in the field.

The distribution of void fraction is shown in Fig. 9 , which describes particle movement and volume change. Particles were displaced radially under the effect of fluid flow and created the gap next to the side walls, the size of which was larger for higher injecting velocity. A localized flow channel was created on the left side of the MP-M sample, which could be due to the fingering effect as discussed earlier. The localized deformation pattern affected both fluid and particle streamlines for the MP-M case as there was smaller resistance to solid and fluid flows in this zone. There was no fingering created for the MP-H case, yet the particle trajectories were diverted more to the left side in the intermediate timestep similar to the MP-M case. The fingering effect was also not observed for the low injecting velocity, the particle trajectories are more symmetric in this case.

Regions of high void fraction values were located near the outlet point for all three injecting velocities, which would be the result of sand production. The cumulative sand curves in Fig. 7  indicate a single stage of a transient sand production. There was no formation of the second damage zone as in the case of single heavy oil (SO) in the previous section, and hence there was no second sand burst during the simulation.

Contact bonds were formed only once during the sample preparation, and their number decreased during the perforation and sand production stages. The bond number for each particle is shown in Fig. 10 . In early timesteps, particles with high bond numbers were located near the boundary as the bonds in the central region were destroyed during perforation. As water was injected radially, particles were displaced radially under the hydrodynamic forces, which left an empty gap with the sidewall as shown in Fig. 10 . This is more pronounced with higher injecting velocities and agrees with the results in Fig. 9 . The blue core of unbonded particles became smaller, some particles of which could be eroded and left the specimen and formed a damage zone of a bowl shape under the outlet point as shown in Fig. 9 . This structure was less symmetric for the MP-M case, with the creation of a flow channel on the left side due to the fingering effect.

A high-pressure consolidation system was utilized to carry out the experiment on sand production in a two-phase flow (DCS Computing GmbH 2021). The equipment allowed for independent control of the saturating and injecting fluids as well as applying realistic reservoir conditions for testing the materials. A chemical reaction between sodium silicate solution and CO₂ gas created a cement bond in the artificial sandstone. Dry sand was mixed with sodium silicate solution, and the mixture was placed inside the cell under a 1 MPa load to consolidate the mixture. CO₂ gas was percolated through the mixture to prepare the artificial sandstone sample. The sample was then subjected to fluid saturation and injection to study sand production under the effect of multiphase fluid flow. The experimental sample has the following properties: diameter = 300 mm, height = 250 mm, porosity = 0.38, and unconfined compressive strength = 0.12 MPa.

The large specimen size and the scale of the experiment created constraint and complication to the multiphase sand production test. Certain simplifications were used; motor oil was used to saturate the sample, while tap water was used for injecting fluid. Oil saturation was conducted using an upward oil flow from the bottom of the intact sample under a pressure of approximately 30 kPa. The process was maintained until the oil volume of roughly twice the specimen’s pore volume was collected, which took around 24 hours.

The saturated sample was perforated by drilling into the specimen. The top cap was connected to an outlet tubing through a central opening. The outflow was collected in an external system. The outlet valve was connected to the atmosphere such that the outlet pressure was atmospheric pressure. Water was injected through the inlet ports on the cylindrical sidewall to simulate radial water flow through the oil-saturated sandstone sample. The inlet pressure was controlled to simulate different drawdown pressures as in the well. The mass of the produced sands was measured using a continuous sand measurement unit. The experimental specimen is shown in Fig. 11  together with the location of fluid inlet and outlet ports. Further details of the experimental procedures can be found in Shabdirova et al. (2022).

As the multiphase sand production experiment was time-consuming and labor-intensive, the experimental results were limited to two drawdown pressures of 800 kPa and 1600 kPa. For the lower drawdown pressure, a transient sand production was observed with an initial sand burst followed by a gradual decrease until a complete cease after 40 minutes. For the higher drawdown pressure of 1600 kPa, the sand rate increased in an uneven manner, eventually leading to clogging of the outflow tubing after 30 minutes.

The numerical and experimental sand productions are plotted in Fig. 12 . As the dimension and time scale were different, a common horizontal time axis was plotted with different time values normalized by the maximum timestep in each case. The cumulative sand masses in the experiment were normalized by the final sand mass of the 800-kPa case, while for the numerical simulation, they were normalized by the final sand mass of the MP-L case. It is noted that the boundary conditions were different between the simulation and the experiment. The drawdown pressure or pressure gradient was maintained in the experiment, and varying outlet flow rate was measured. On the other hand, injecting and outflow velocities were kept constant in the simulation. The measured flow rates in the 800-kPa drawdown pressure experiment cover the range of the MP-L and MP-M simulations as shown in Table 4 , whereas the MP-H case overlaps with the lower range of the measured flow rate for the 1600-kPa drawdown pressure, which was obtained in the early stage of the experiment.

The numerical simulation and the 800-kPa experimental results showed a single-stage sand production. The 1600-kPa experiment exhibited a two-stage sand production behavior similar to the heavy oil (SO) simulation result in Fig. 4 . This could be related to another development of the damage zone leading to a second sand burst. The 1600-kPa experiment was started as a typical transient sand production until the measured flow rate and sand volume increased abruptly leading to a catastrophic sand production which clogged the tubing system and the test was terminated. The clogging tubing could also be a disastrous problem in actual field operation (Vásquez H. et al. 1999; Willson et al. 2002; Fjaer et al. 2008). The new sand burst would require higher input energy either by a higher injecting velocity or a higher fluid viscosity. These conditions were not met in the 800-kPa experiment as well as in the multiphase simulations.

In this study, we utilized the coupled CFD-DEM modeling to investigate the sand production characteristics in a multiphase flow of a cylindrical sandstone sample. The macroscopic sanding behavior was related to the microscopic features of bond breakage, void fraction, oil cut, and fluid and particle velocity trajectories across the sample.

The injection of water into a sample initially saturated by heavy oil resulted in a multiphase flow that reached the maximum cumulative sand production earlier than the cases of single fluid flow for both heavy oil and water for the same injecting velocity. The multiple phase flow condition created different localized deformation and erosion patterns that affect both fluid flow and particle movement.

The effect of water breakthrough was investigated with increasing water injecting velocity. Water breakthrough was only observed for medium injecting velocity and above, but not for lower velocity, which confirmed the role of a critical velocity value that could trigger water breakthrough in the literature. Furthermore, water breakthrough was observed in a particular phenomenon of water coning as observed in other experiments. Higher injecting velocity would lead to earlier water breakthrough, but it is not related to sand production.

High injecting velocity and higher transport capacity of the heavy oil would result in higher drag force that would lead to additional damage zone development and hence cause multistaged sand production.

A comparison with the experimental results on a similar material and of similar flow conditions in our laboratory showed that higher injecting velocity could result in multiple sand bursts and could lead to catastrophic sand production as it clogged the tubing system in the laboratory as well as in the field.

The coupled CFD-DEM modeling approach provided new insight into the multiscale mechanism of sand production. Similar deformation characteristics were obtained from the observation in the laboratory and in the field, which would have practical implications for production operation and enhanced oil recovery.

This research was supported by Nazarbayev University CRP Grant No. 11022021CRP1506.

Original SPE manuscript received for review 29 March 2024. Revised manuscript received for review 10 August 2024. Paper (SPE 223586) peer approved 29 August 2024.