Abstract
This paper aims to develop a numerical model that can be used in sand control during production phase of an oil and gas well. The model is able to predict not only the onset of sand production using critical bottom hole pressure inferred from geomechanical modelling, but also the mass of sand produced versus time as well as the change of porosity versus space and time using hydromechanical modelling. A detailed workflow of the modelling was presented with each step of calculations. The empirical parameters were calibrated using laboratory data. Then the modelling was applied in a case study of an oilfield in Cuu Long basin. In addition, a sensitivity study of the effect of drawdown pressure was presented in this paper. Moreover, a comparison between results of different hydromechanical models was also addressed. The outcome of this paper demonstrated the possibility of modelling the sand production mass in real cases, opening a new approach in sand control in petroleum industry.
1. Introduction
Sand production occurs in many oil fields across the world and it is especially common in the porous sediments. Sand production observed on the surface occurs as a series of three events that happen at downhole area: (1) formation failure, (2) sand erosion due to flow, and (3) sand transport.
(1) Formation Failure. In situ stresses and pore pressure act on formation sands and under certain conditions; the criteria for failure are met. The presence of the wellbore and perforations causes a concentration of stresses near these cavities and deformation and failure can occur under certain well-known conditions. This criterion is bottomhole pressure that makes the maximum effective tangential compressive stress equal or higher than the rock strength (failure criteria); significant sanding begins at some point (the onset).
(2) Sand Erosion due to Flow. Damaged regions that have failed (meeting the failure criteria) face additional stresses caused by pore pressure gradients. The process of sand erosion is essential for the sand to be removed from the failed region and then to be entrained with the fluid.
(3) Sand Transport. Sand erosion detaches sand grains into a perforating cavity or wellbore. Some of the grains are transported to the surface while others settle into the perforation tunnel or into the well hole.
During and after sand production, wells can sand-up and that has different effects on the productivity. At first, the productivity seems to increase due to the increase of permeability. However, after a while, the sand produced can obstruct the entrance of hydrocarbon into the wellbore, especially for cased perforated wells using sand screens or gravel pack. Disposal of produced sand is also a significant cost associated with sand production. Finally, sand can be transported to the surface which causes erosion of pipe lines, joints, chokes, and valves. So, if the prediction of sand production is identified, it will help operators to manage the situation properly and prepare suitable treatment methods for the well.
This study focuses on the first and second steps of sand production because they have the most impact and they are not well understood. The main reasons why in petroleum industry nowadays we still do not predict the mass of sand production are because of the complexity of numerical models (hence the lack of professional software in this domain) and the unavailable real data of sand production due to the difficulty in collecting this kind of data in the oil field. Therefore, most of the studies predicting the mass of produced sand still stay at the laboratory step. In real life, in petroleum companies’ reports, only Geomechanical models are being used to predict the onset of sand production. This papers aims to bring the application of the Hydromechanical models into a real case in petroleum industry and to combine the use of Geomechanical model, which predicts the onset of sand production, and the Hydromechanical model, which predicts the mass of produced sand.
2. Literature Review
Several studies released models predicting the onset of sand production and the amount of sand produced. Parameters affecting sand production have been discussed for decades. However, there is no clear consensus. In this section, a brief review of major conclusions of these past studies is presented.
Willson et al. [1] developed a model to predict the sand production rate from the onset of sanding model. The onset of sanding is predicted using a stress–based model of shear failure around a perforation or an open hole. Sand production is assumed to occur once the maximum value of the effective tangential stress around the perforation exceeds the apparent rock strength (rock strength is strength with is boost factor, for cased perforated wells = 2; TWC is the Thick Wall Cylinder which is a measure of rock’s strength and is used in sand production study instead of Unconfined Compressive Strength, because TWC reflects more closely the reality of in situ stress sustained by the borehole or the perforation channels than the UCS). No consideration was given to sand transport by drag forces. The model of Willson et al. [1] requires TWC test data for hole size the same as the perforation size. If the perforation is of different size, the application of the model must be careful. Hence the boost factor exists to compensate for the difference between the real data and the test data.
The criterion for sanding is
where CBHP is critical bottomhole pressure; is maximum principal stress; is minimum principal stress; is pore pressure. is poroelastic constant given by with : Poisson ratio and : Biot’s constant.
The model predicts the rate of sand production by utilizing the nondimensionalized concepts of Loading Factor, LF (near-wellbore formation stress normalized by strength), Reynolds number (Re), and water production factor. An empirical relationship between Loading Factor, Reynolds number, and the rate of sand production incorporating the effect of water production was proposed as follows:(see [1]).
In this formula given by Willson et al. [1], SPR is the sand production rate; water cut is the ratio of water produced compared to the volume of total liquids produced.
Although the Willson et al. model [1] takes into account the different phases of the fluid via water cut, the model did not give clear expression of the Sand Production Rate SPR in function of the Loading Factor LF, the Reynolds number, and the water cut.
In 1996, Vardoulakis et al. [2] proposed the following sand production model usingmixture theory, assuming that the sand in place is fully degraded from the beginning and the production is due only to the hydrodynamic forces. Equilibrium equation for solid phase is often ignored. The process is initiated with a very small solid concentration given as a boundary condition. The results are insensitive to this value as long as it is small:where is the porosity, is the rate of eroded solid mass per unit volume, is the solid density, λ is the experimentally evaluated sand production coefficient, is the transport concentration, is the specific discharge in the th direction, and is the notation representing the norm of a vector.
The Vardoulakis model is difficult to solve because of the complexity of the equations. Moreover, the model does not take into account the different phases of the fluid, so the fluid is considered as single phase, which does not reflect the reality of petroleum fluid. Furthermore, the coefficients were not calibrated due to lack of experimental data.
Papamichos et al. [3] developed the following model based on the assumption that failure is due to erosion and porosity increases until it reaches unity. Below is the relation between sand mass and the porosity:
The dimension of the above equation is that is the rate of eroded solid mass per unit volume (g·s−1·m−3), is the solid density (g·m−3), is the porosity, and is time in second (s).
Variation of sand mass due to erosion is given aswhere is sand production coefficient and is plastic shear strain.
The main advantage of Papamichos et al. model [3] is that the authors provided experimental data in order to calibrate the experimental parameters. In this model, not all the plastic deformation areas will produce sand. The sand is produced only when the plastic deformation reaches a limit value (). However, it is practically difficult to determine this limit value in reality. Besides, the system of differential equations of this model is extremely complex with plenty of empirical parameters.
Chin and Ramos [4] developed a sand production model considering that erosionoccurs during sand production as follows, where is the solid velocity:
The main inconvenient of this model [4] is that it was developed only for weak formation. Besides, the model retains only primary physics of rock failure and coupled rock deformation and fluid flow. The model does not include the effects of well configuration and completion, wellbore storage, erosion, interaction of disaggregated solid and flowing fluid, and solid transport through the porous medium.
The analytical model developed by Fjær et al. [5] combined a theoretical model with an empirical relationship as shown below:where the porosity increases with time:where is the flow rate, is the critical flow rate, is the initial porosity, is the initial time, and is proportionality constant; has the dimension of s/m3.
Gravanis et al. [6] developed a coupled stress-fluid flow erosion model: where , with dimension of inverse length, represents the strength of the erosive processes that lead to sand production. is an exponent parameter of the model. is an erosion function defined as follows:
is the usual Euler Gamma function and is the depth of the plastic region. The profile function approaches a step function as increases. The exponent is fixed at the beginning of the analysis and can be tuned further as part of the calibration procedure. In this work we will fix , which is the same value used in the work of Fjær et al. [5].
In 2010, Isehunwa and Olanrewaju [7] proposed a new model for sand production, considering the effects of flow rate, fluid viscosity and density, grain size, and cavity height. Sand is produced by drag and buoyancy forces which predominantly act on the sand particles. The radius of sand production cavity is
The volume of sand produced can be expressed aswhere is the fluid flow rate, is the grain radius, is the cavity height, and is the gravitational acceleration.
Among these models, the ones of Fjær et al. [5] and Gravanis et al. [6] were chosen for this study because of their clarity in the explanation and equations, which are necessary for us to be able to numerically solve the problems. In addition, these models have some advantages in comparison with others. They proposed the basic theory for hydrodynamic erosion of sandstone which is based on filtration theory. They adopted full strength hardening/softening behavior of reservoir stone. They also do take into account the grain size, the gradient elastoplasticity for thick wall cylinders and the cavity failure around boreholes.
The models were solved using the workflow developed in Section 3 of this paper; then the coding was made using MATLAB software.
3. Workflow for Calculation and Calibration
3.1. Hydromechanical Erosion Model of Gravanis et al. [6]
According to Gravanis et al. [6], the basic assumptions are as follows:(1)Fluid flow can be described by Darcy’s law.(2)We define the mathematical time by for simple calculations. is the real time while is the mathematical time which is introduced to facilitate the resolution of the problem. (Figure 1). The function is related to pressure drop as follows: is the permeability.(3)The whole region is divided into a plastic region and an elastic region. In elastic region, we apply Hooke’s law and in plastic region we consider the Mohr-Coulomb failure criterion.(4)Under the condition that plasticity of the material is damaged and subject to decohesion, it can be eroded under weak hydrodynamic forces. It happens when drawdown pressure (DP) exceeds a critical drawdown pressure (CDP).
3.1.1. Calculation Workflow
Step 1 (at ). (i)Calculate function : where is related to the flow rate by the relation , .(ii)Calculate pressure :(iii)The depth of plastic region : the elastic and plastic regions are presented in Figure 1 according to Gravanis et al. [6]. In elastic region, combining with the boundary condition , we have
varies from to . In plastic region, we consider the Mohr-Coulomb failure criterion and combine with the boundary condition σ: where is material cohesion [MPa]. is Biot’s constant variant from 0 to 1; here we shall set this value to 1 meaning that we will neglect any compressibility effects. is the principal stress ratio which is equal to with being the horizontal stress and the vertical stress. At this stage, there are two continuity conditions and two unknowns, the integration constant and the location of the plastic zone boundary location . The conditions are shown below: is given explicitly in terms of which is solved numerically for [6]. Note that (18) are purely algebraic. Thus, the plastic region depth is determined as .Step 2 (at time ). (i)The porosity field is calculated from erosion model equation (9): The formula involves the integral which is approximated by and . () is known from the previous time step where is an exponent parameter of the model and must be tuned by calibration.(ii)Calculation of porosity and radius of sand production cavity: Radius of sand production cavity is calculated from wellbore to the distance where porosity of formation is constant and equal to initial porosity .
(iii)Calculate function , pressure , and the plastic region depth from (14) to (17).(iv)The erosion strength is : where , , and are model parameters that must be tuned by calibration, is initial erosion strength [m−1], is maximum erosion strength [m−1], and is threshold depth of the plateau.(v)Calculate sand production mass rate: Note that this equation must be used in the next time step along with .Step 3 (at the general time ). (i)The porosity field is calculated from (20) and (21), where the integral is approximated by . This involves the values of the depth at all the previous time steps, which are known. As explained above, the function , the pressure, and the new value of the plastic region depth can then be calculated.
3.2. Analytical Model of Fjær et al. [5]
The analytical model is based on these assumptions:(1)The driving mechanism for continuous sand production is erosion from plastified material in the vicinity of the production cavity.(2)The sand production rate depends on (1) how much the well pressure is reduced below the critical sand production pressure, (2) the fluid flow rate and the fluid viscosity, and (3) the cementation of the rock.(3)The sand in place is fully degraded from the beginning and the production is due only to the hydrodynamic forces. These forces are proportional to the fluid pressure drop over the volume element, and the pressure drop is proportional to the fluid flow rate as specified by Darcy’s law; moreover the permeability is given by the Kozeny-Carman equation; we have where μ is the fluid viscosity. The proportionality constant has the dimension of s/m3(4)Firstly, it is required that DP > CDP, which expresses the fact that stress induced damage of the rock is a necessary condition for sand production. Secondly, it is required that (critical fluid flow rate), which expresses that the hydrodynamic forces must be strong enough to uphold the erosion process.
3.2.1. Calculation Workflow
Following Fjær et al. [5] we simply assume that the stiffness of the rock remains constant until the porosity has reached a critical value . At that point, the entire part of the rock that has been producing sand collapses, and the remaining solid material in that part is produced in one burst. The stresses around the cavity are then redistributed, and the situation is similar to the initial state at , except the cavity has become a little larger.
Step 1 (continuous sand production ). (i)Calculate porosity from (8):(ii)Calculate volumetric sand production : we assume that sand production is related to plastification of the rock, so that only parts of the rock that have suffered some plastic deformations are in a state where sand can be produced. where is radius of the cavity and is sand production zone (Figure 2).(iii)Sand production mass:
Step 2 (instantaneous sand production ). Once the entire part of the rock that has been producing sand collapses, the remaining solid material in that part is produced in one burst.
The collapse of the sand producing zone implies that the radius of the cavity increases, from to , where can be found from (27):(i)Cumulative sand production: sand mass is given as the initial amount of solid material in the sand producing zone; that is, Now the stresses around the cavity are redistributed, and the situation is the same as it was at . Thus, the erosion process starts a new cycle.
3.3. Calibration of Empirical Parameters of the Models of Fjær et al. [5] and Gravanis et al. [6]
We use experimental data profile of Papamichos et al. [3] to calibrate empirical parameters (Table 1).
3.3.1. Calibration of Empirical Parameters of the Gravanis et al. Model
We have empirical parameters: , , , and and we slightly varied these parameters around the values obtained by trial and error from calibrating with test data of Papamichos et al., 2001 [3]. The calibration results are shown from Figures 3 to 6. For each calibration, the model curve is compared with the experimental curve.
From the results presented from Figures 3–6, we choose suitable values of empirical parameters which are: the exponent coefficient (), maximum erosion strength ( = 0.088 m−1), initial erosion strength ( = 0.42 m−1) and threshold depth of the plateau ( m). The small discrepancy with the experimental data can be neglected.
3.3.2. Calibration of Empirical Parameters of Fjær et al. Model
The model of Fjær et al. has three empirical parameters: , , and . Similar to the model of Gravanis et al. [6], we use experimental data of Papamichos et al., 2001 [3] and then run the model and compare with the corresponding experimental curve to calibrate empirical parameters. Results are presented in Figures 7, 8, and 9.
Figures 7, 8, and 9 indicate that the erosional process gradually increases until critical porosity is reached. Each step is followed by a long period of continuous production at a low rate. When the porosity has reached a critical value , at that point, the collapsed zone of the remained solids is produced, and process of sand production is immediate and rapid. Then the erosion process starts on a new cycle.
In fact, the result of the model seems very different than the experimental result because Fjaer’s model considers a step of “collapse” (Step : instantaneous sand production), but overall, the sand production mass over time is matched between the model and the experimental data and that is why we could calibrate the parameters.
These results also allow us to choose suitable values of empirical parameters which are critical porosity ( = 0.4), critical fluid flux ( = 0.001 m/s), and sand production coefficient ( = 600 m3/s).
4. Case Study
An application for the oilfield X in Cuu Long basin using Geomechanical model of Willson et al. [1] to predict the onset of sand production condition is now presented. In addition, it was combined with the Hydromechanical model of Fjær et al. [5] and Gravanis et al. [6] to predict sand production rate; then the results of the models will be compared.
After collecting and processing data by analyzing log and other parameters, input data are shown in Table 2. The petroleum industry, for historical reasons, has been landed with a mixture of US, British, and SI units, which is often referred to by “oilfield units.” In Vietnam they follow this tradition of using oilfield units, which is a requirement by the government; hence in Table 2 we do not use SI units.
4.1. Calculate Critical Bottomhole Pressure Using Geomechanical Model of Willson et al. [1]
The critical bottomhole pressure CBHP is important information during production phase because it indicates the lowest bottomhole pressure for sand production to not occur. For a specific well in oil field, the only production data that we control is the bottomhole pressure, which is adjusted using surface chokes. The pressure is controlled; hence the flowrate is controlled. For this reason, before predicting the sand production using Hydromechanical models of Gravanis et al. and Fjær et al., we firstly calculate the CBHP using Geomechanical model of Willson et al.
The testing showed that a relationship between the effective in situ strength of the formation U and the TWC (Thick wall cylinder) strength relative to a specimen with an OD/ID ratio (Outer/Inner diameters) would be equivalent towhere the boost factor is often taken as for cased–hole perforated wells. Equation (1) considers effect of reservoir pressure decline [9] which can be accounted for by updating in situ stresses, while the vertical stresses are usually kept constant:where with is the Biot’s constant and is coefficient of Poisson, with is current reservoir pressure, and is initial reservoir pressure.
Figure 10 presents the bottomhole pressure versus reservoir pressure . In Figure 10 we also present the CBHP when the UCS changes. For example, for UCS = 1950 psia, if reservoir pressure is 6000 psia, the CBHP is about 1500 psia; if the bottomhole pressure is lower than this value 1500 psia, sand production will occur.
The results showed that if UCS gets smaller, the sand production window is greater. This can be explained by the reason that the formation rocks around the wellbore are weakened when UCS decreases. Thus, risk of sand production increases.
We consider a case study with UCS = 2450 psi, = 4000 psi. From Figure 10, we see that if the current reservoir pressure is 4000 psi, then the critical bottomhole pressure is 2020 psia. For sand production to not occur, the bottomhole pressure must be greater than 2020 psi. Otherwise, to increase the flow rate, we must increase DP and if it is greater than CDP, it means that bottomhole pressure must be lower than critical bottomhole pressure (1000 psi); then the sand will be produced. So there is a compromise between increasing flowrate (production requirement) and sand production.
4.2. Calculate Sand Production Mass Using the Models of Gravanis et al. and Fjær et al.
We use Hydromechanical Erosion models of Fjær et al. [5] and Gravanis et al. [6] to calculate sand production mass versus time. However, we do not have experimental data of Well X1 to calibrate empirical parameters, so we must use empirical parameters already calibrated in Section 3. In Figure 11 the results of sand mass produced over time calculated by the two models are presented. We also did a sensitivity study of the effect of drawdown pressure. Drawdown pressure, which is the difference between the reservoir pressure and the bottomhole pressure, is the most important data that we control during production phase. We control the drawdown; consequently, we control the production rate. The drawdown is typically controlled by surface chokes. The more the choke is open, the bigger the bottomhole pressure is and, hence, the smaller the drawdown is, and vice versa.
Figure 11 shows that sand production mass increases when drawdown pressure increases. However, results given by the model of Gravanis et al. are smaller than results given by Fjær et al. At each bottomhole pressure, sand mass given by the model of Gravanis et al. increases until it reaches constant value. On the other hand, the one of Fjær et al. increases continuously because the process of erosion repeats which causes the radius of cavity and production zone to increase, so sand production mass increases continuously. These two models show different results although they use the same physical mechanisms: the fluid flow can be described by Darcy’s law and the driving mechanism for continuous sand production is erosion from plastified material in the vicinity of the production cavity under hydrodynamic forces. However, the model of Fjær et al. uses a step of instantaneous sand production which is not considered in the model of Gravanis et al. Moreover, while Gravanis et al. divide the whole region into a plastic region where Mohr-Coulomb failure criterion is considered, and an elastic region where Hooke’s law is applied, Fjær et al. consider only the plastic region. Fjær et al. simply assume that the stiffness of the rock remains constant until the porosity has reached a critical value . At that point, the entire part of the rock that has been producing sand collapses, and the remaining solid material in that part is produced in one burst. As a result, the equations used in these two models as well as the associated experimental parameters are quite different and eventually led to quite different results.
Finally, it is important to recall that, due to the lack of validation data in reality, it is impossible to choose which model to use in practical use. At this state, these results can only be used for illustration in study, and to demonstrate the possibility of using Hydromechanical models in real cases to predict eroded sand mass over time. In the future if real data is available, these models can be revised and recalibrated. Unfortunately, until now, measuring sand mass rate is technically impossible in the oilfields. An idea was proposed to solve this problem; the determination of real sand rate may be based on the erosion rate of the choke. However, this study has never been realized and may constitute a new objective in our future study.
5. Conclusion
This study combines Geomechanical model and the two Hydromechanical Erosion models of Fjær et al. [5] and Gravanis et al. [6] aiming to predict sand production. The proposed model can not only predict the critical pressure for onset of sand production but also estimate sand production mass and the variation of porosity in time and space. In addition, the study also gives detailed calculation steps for the two Hydromechanical models. Although each model has assumptions and empirical parameters, we can adjust them if we have experimental data. The results of this study will help to make an effective production planning that can avoid sand production when reservoir pressure declines. Hydromechanical model is still very new and not yet widely used, due to several reasons such as complicated calculation methods and parameters that need to be experimentally adjusted in practice or in the laboratory. Therefore, this initial computational workflow using MATLAB will significantly contribute to the development of future research in this domain in Petroleum Engineering.
The experimental data (sand production’s mass in function of time) is rarely done not only in Vietnam but also in the world (only some laboratories and authors mentioned in the paper have done this kind of experience because of their research in this field). Hence, it is currently impossible to obtain empirical data for Cuu Long basin (sand mass in function of time), not only data in laboratory but also production data due to difficulties in collecting and measuring sand production rate. Therefore, application of the calibrated model for Cuu Long basin in this study is destined to demonstrate the possibility of using Hydromechanical models in real cases and the result can be used for illustration in research. In the next studies when we have available data from Cuu Long basin, we can do a revision/recalibration of these models. We already thought about the determination of real sand rate data based on the erosion rate of the choke. This idea will be the objective of our future study.
Notations
| CBHP: | Critical bottomhole pressure (psia) |
| DP: | Drawdown pressure (psia) |
| LF: | Loading Factor (near-wellbore formation stress normalized by strength) |
| SPR: | Sand Production Rate (g/s) |
| TVD: | True vertical depth (ft) |
| TWC: | Thick wall cylinder strength (psia) |
| UCS: | Unconfined compressive strength (psia) |
| : | Notation representing the norm of a vector |
| : | Biot’s constant |
| : | Exponent parameter |
| : | Plastic shear strain |
| : | Usual Euler Gamma function |
| : | Experimentally evaluated sand production coefficient (1/m) |
| : | Proportionality constant (s/m3) |
| : | Initial erosion strength (m−1) |
| : | Maximum erosion strength (m−1) |
| : | Depth of the plastic region (m) |
| : | at (m) |
| : | Threshold depth of the plateau |
| : | Viscosity (MPa·s) |
| : | Porosity |
| : | Initial porosity |
| : | Critical porosity |
| : | Perforation orientation from the top of the wellbore in deviated wells (degree) |
| : | Solid density (g/m3) |
| : | Maximum principal stress (psia) |
| : | Minimum principal stress (psia) |
| : | Vertical stress (psia) |
| : | Minimum principal horizontal stress (psia) |
| : | Maximum principal horizontal stress (psia) |
| : | Minimum horizontal stress direction (degree) |
| : | Poisson ratio |
| : | Integration constant |
| : | Young’s modulus |
| : | Well inclination (degree) |
| : | Permeability (mD) |
| : | Principal stress ratio |
| : | Perforation length (ft) |
| : | Mass of produced sand |
| : | Rate of eroded solid mass per unit volume (g/s·m3) |
| : | Rate of eroded solid mass (g/s) |
| : | Initial amount of solid material in the sand producing zone (g) |
| : | Pore pressure (psia) |
| : | Internal pore pressure (Mpa) |
| : | External pore pressure (Mpa) |
| : | Bottomhole pressure |
| : | Reservoir Pressure |
| : | Flow rate (m3/s) |
| : | Specific flow rate in the th direction (g/s) |
| : | Critical flow rate (m3/s) |
| : | Perforation radius (ft) |
| : | Cylinder internal radius (m) |
| : | Cylinder external radius (m) |
| Location of the plastic zone boundary (m) | |
| : | Radius of the cavity (m) |
| : | Radius of the sand production zone (m) |
| : | Radius of the plastic zone (m) |
| Re: | Reynolds number |
| : | Material cohesion (MPa) |
| : | Real time |
| : | Mathematical time |
| : | Initial time |
| : | Solid velocity (m/s) |
| : | Volumetric sand production (m3) |
| : | Well azimuth (degree). |
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.