Abstract
In this present work, cement slurry flowing in well annuli with narrow or varied radial sizes is investigated, whose purpose is to precisely calculate the equivalent circulating density (ECD) during well cementing in the petroleum industry. The theoretical Metzner-Reed (MR) method and the direct numerical simulation (DNS) are applied to calculate the cementing ECD for different types of cement slurries by considering both the Bingham and Power-Law rheological models. The cementing ECD calculated by the theoretical MR method is verified to be in accordance with the DNS result. Moreover, for the narrow annulus, the cementing ECD is demonstrated to be sensitive to the rheological model. In addition, based on the obtained DNS results, the accuracies of the cementing ECD calculated by three commercial software, including DRILL BENCH (DB), LANDMARK (LM), and PVI, are compared, to determine their corresponding adaptability for different rheological models. According to the detailed comparisons, LM shows the highest accuracy in calculating the cementing ECD for the Power-Law fluid, while PVI is optimal for the Bingham fluid. However, sharp fluctuations of the calculation errors are observed when DB is applied.
1. Introduction
Cementing ECD is significant for well cementing in the petroleum industry, since it determines the stability of the pressure system in the well annulus [1]. A high ECD can lead to formation fracture, while a relatively low ECD might cause borehole wall sloughing. Comparing with conventional well cementing, when a well annulus is narrow or its radial size varies along the well depth, a higher accuracy of the cementing ECD is necessary. For instance, sudden shrink of the radial size of the well annulus will accelerate the flow of the local cement slurry, which increases the frictional pressure loss, and thus leads to a higher cementing ECD. Under such a condition, the risk of formation fracture increases.
The accurate calculation of the frictional pressure loss is determinant for the prediction of cementing ECD; however, its influential factors are complex and various, including fluid properties, well configuration, and process parameters, etc. According to the basic model for calculating the frictional pressure loss [2], a number of modified models have been proposed, which comprehensively considered the influences of a variety of factors, including the flow regimes [3–5], rheological models [6–8], and operation conditions [9–12], etc. Specifically, Haciislamoglu and Cartalos [13] provided the corresponding prediction models of the pressure losses for the laminar flow, transition flow, and turbulent flow. Starting from analysis of the friction factor, Hansen et al. [14] established the theoretical model about the frictional pressure loss for different flow regimes. Hashemian et al. [15] proposed a prediction model for the frictional pressure loss of the Yield-Power-Law fluid flowing in the annular region. Ozbayoglu and Sorgun [16] provided a correlation between the friction factor and Reynolds number, based on which the frictional pressure loss for the non-Newtonian fluid in realistic annuli can be estimated. Bailey and Peden [17] calculated the frictional pressure loss by considering the Sisko rheological model. Wang et al. [18] and Sun et al. [19] applied the DNS method to simulate the flow of cement slurry in well annuli.
In addition, Ooms et al. [20, 21] discussed the influence of drillpipe rotation and eccentricity on the frictional pressure loss during drilling. Enfis et al. [22] analyzed the hydraulic effect of the tool joint on annular pressure loss. Ahmed and Miska [23] conducted experiments to observe the Yield-Power-Law fluid flowing in well annuli with drillpipe rotation. According to the analysis of experimental data, the frictional pressure loss was verified to be decreasing as the rotation speed of the drillpipe increased. Ahmed et al. [24] proposed a semiempirical model for calculating the frictional pressure loss, which was verified based on field-measured data. Although these publications focused on drilling mud rather than cement slurry, the basic logic for calculating the frictional pressure loss is similar; thus, the developed methods about drilling mud can be also referenced by cement slurry.
In practical field applications, a number of commercial software [25], including DB, LM, and PVI, are generally utilized to calculate the frictional pressure loss and, thus, the cementing ECD. However, the cementing ECD results obtained from different types of commercial software do not agree with each other, and sometimes, obvious discrepancies can be even observed. Moreover, according to the literature review, there is no systematical investigation on this topic. Under such circumstances, comparisons of commercial software on the calculation of cementing ECD deserve to be conducted, so that the adaptability and reliability of the commercial software on the design of well cementing can be determined.
The rest of this paper is organized as follows. Primarily, the MR theory is introduced in detail in Section 2. Then, the calculation of cementing ECD via the DNS method using FLUENT is conducted in Section 3, where the independence for grid meshing is studied. Thereafter, the cementing ECD for the annuli with narrow radial sizes and with varied radial sizes are discussed in Section 4 and Section 5. In particular, the calculations of cementing ECD by using the MR method, DNS, DB, LM, and PVI are compared. The characteristics of cement slurries flowing in well annuli with narrow or varied radial sizes are analyzed. The adaptability of commercial software for different rheological models is discussed. At last, concluding remarks are provided in Section 6.
2. MR Theory
The main rheological models applied for cement slurry include the Bingham model, Power-Law model, and Yield-Power-Law model. According to a series of numerical simulations, under the same conditions, the calculated ECD for the Power-Law model and that for the Yield-Power-Law model are close to each other; hence, in this paper, only the Bingham model and the Power-Law model are considered in the following sections.
The rheological model for the Bingham fluid is where is the shear stress (Pa), is the yield stress (Pa), is the plastic viscosity (), and is the rate of shear (1/s).
The rheological model for the Power-Law fluid is where is the consistency coefficient () and is the fluid behavior index.
The cementing ECD consists of two parts: equivalent static density (ESD) from hydrostatic pressure and frictional pressure loss: where is the frictional pressure loss (Pa), is the acceleration of gravity (9.81 m/s2), and is the vertical depth of the wellbore. Therefore, in order to obtain the cementing ECD, the key point is to calculate or the corresponding nondimensional Fanning friction factor : where is the wall shear stress (Pa); is the fluid density (g/cm3); and are the outer and inner diameters of well annuli, respectively; is the length of a certain section of wellbore; and is the average of the sectional flow speed (m/s): where is the flow rate (m3/s).
According to the MR Reynolds number () proposed by Metzner and Reed [3], where the hydraulic equivalent diameter is and , , and are the local Power-Law fluid behavior index, local Power-Law consistency coefficient, and nondimensional shear stress, respectively.
In this paper, only the laminar flow is considered; thus, the correlation between the Fanning friction factor and the MR Reynolds number is written as where the critical Reynolds number of turbulence can be calculated to verify the flow regime as the laminar flow:
Once the Fanning friction factor is obtained, the frictional pressure loss can be calculated according to Eq. (4), thus cementing ECD via Eq. (3).
3. Direct Numerical Simulation
As introduced in Section 2, the cementing ECD can be calculated via the MR theory. However, the accuracy of the MR theory in the calculation of cementing ECD, especially for well annuli with narrow or varied radial sizes, needs further verification. Under such circumstances, the DNS method is considered to prove it. In particular, a software package of computational fluid dynamics, FLUENT, is applied to conduct the DNS, and the obtained ECD is then compared with the corresponding result from the MR theory.
3.1. Governing Equations
In the paper, only steady flow is considered; thus, in Cartesian coordinates, incompressible continuity and Navier-Stokes equations can be written in the nondimensional forms as where is the gradient operator, is the velocity vector, and is the tensor of the stress, and for steady flow, .
3.2. DNS Procedure
WORKBENCH of ANSYS 17.0 is applied to build up the analysis process for the DNS method. Primarily, the module, GEOMETRY, is utilized to establish a three-dimensional (3D) geometric model for cement slurry flowing in the well annulus. Subsequently, the module, ICEM-CDF, is applied to mesh the developed geometry model. An example of a meshed structure model is shown in Figure 1(a). In particular, due to the significant influences of both the inner and outer walls of the well annulus on the flow of cement slurry, the meshes close to both the walls are refined. As per the cross-section shown in Figure 2(b), the grid sizes for both the inner and outer walls are set as , which then expand towards the center of the section as an expansion ratio . Under such circumstances, the obtained output wall function is always less than 5, and thus, the wall function is kept as the default.


The well-meshed structure model is imported into FLUENT; then, the numerical simulation is carried out by using the finite volume method (FVM). The incompressible flow is solved by the pressure-based segregated solver implicitly, and the gradient calculation based on the Green-Gauss node method is adopted. The pressure and momentum equations are discretized with the second-order scheme and the second-order upwind scheme, respectively. The pressure-velocity coupling adopts the SIMPLEC algorithm without any skewness correction due to the good quality of the structural mesh. The maximal convergence residuals are set as for the continuity equation and for three components of the momentum equations.
In addition, since only the steady flow with low Reynolds number is considered in this paper, the initial velocity field is the still flow with all velocity components of zero. The inlet is a velocity boundary, and the velocity of the axial incoming flow is equal to the average velocity of the inlet section. The outlet is a boundary with free outflow, namely, the velocity component along the normal gradient is zero. Both the walls of well annuli are nonslipping boundaries, and thus, their velocities are constant zero. The reference pressure is pointed at the center of the inlet section, which is also a constant zero.
3.3. Study of Independence
In order to secure the computational accuracy with low computing cost, the independence of the computational result is studied, whose target is to determine the maximal grid size and the minimal axial computational domain. The computational domain is expressed by the ratio between the length of the well annulus and its hydraulic equivalent diameter (). The adjustable parameters for the grid size include the boundary layer grid () and axial grid (). Once these three parameters are determined, the grid numbers for the cross-section (), the axial direction (), and the whole model () can be obtained. In addition, the maximal axial velocity at the center of the outlet, , is chosen as an indicator to examine the influences of both the grid size and the axial computational domain on the flow of cement slurry in well annuli.
The results of the independent studies are listed in Table 1. On the one hand, as the grid size decreases, increases gradually. For the axial computational domain , when , and , reaches 14.277, while further decreasing to 0.002 or even 0.001, the grid number expands rapidly; however, the change of is negligible; hence, the optimal grid sizes are confirmed as and . On the other hand, when the length of the annulus is extended to or even , the corresponding increases to 14.288 or 14.289, respectively, whose differences with for are negligible. Namely, once the length of the well annulus reaches , the fluid flow can achieve a stable state on the outlet, and such stable state will not be affected by the further increase of the annulus length. Based on the above discussion, when the parameters are set as the fifth case in Table 1, namely, , , and , the influences of both the grid size and the axial computational domain on the fluid flow have already been eliminated; moreover, the computing cost and the computational accuracy are well balanced. Hence, this group of parameters will be used to mesh the structure models in the following sections.
3.4. Frictional Pressure Loss
According to the MR theory, the flow speed determines the MR Reynolds number and, thus, the Fanning friction factor, which eventually determines the cementing ECD. Hence, in order to accurately calculate the cementing ECD, the distribution of the flow speed along the well annulus should be analyzed primarily. The axial flow velocities in the central lines of the well annuli with and are depicted in Figure 2. It can be observed that for the case , the fluid flow has already fully developed, and a horizontal region for steady flow appears in . The further extension of the well annulus expands the horizontal region for the steady flow, such as the steady region for the case and the steady region for the case . Hence, once the axial computational domain is longer than , its influence on is negligible. According to this understanding, for a realistic well annulus whose length is up to several kilometers, it can be analyzed by a considerably shortened structure model. However, if the radial size of the annulus varies along the well depth, it needs to be simplified by a combination of several shortened parts with different radial sizes. By this way, a full-scaled well annulus can be replaced by a small-scaled structure model for numerical simulation; moreover, both the computational accuracy and the simulation efficiency can be guaranteed simultaneously.
In addition to the steady region, Figure 2 also shows another two regions: the developing region and the dropping region. In order to further analyze the fluid flow in these three regions, the distributions of the flow velocities in different directions are displayed in Figure 3. Figures 3(a)–3(c) correspond to the developing region for , the steady region for , and the dropping region for , respectively, and show the variations of the flow velocities in the axial, radial, and tangential directions on one side of the cross-section of the well annulus. As can be seen, in all the three regions, the profile of the axial velocity is parabolic related to the radial coordinate. However, both the radial velocity and the tangential velocity are different for the three regions. When the fluid rises up to , the radial velocity and tangential velocity are almost equal to zero, and this situation maintains until . When the fluid climbs up close to the outlet where obvious fluctuations of both the radial velocity and tangential velocity appear, and also due to their disturbances, the axial velocity displays a sudden drop when approaching the outlet in Figure 2.

Based on the characteristics of the distributions of the flow velocities shown in Figures 2 and 3, the flow speed remains constant in the steady region which starts from downstream of the inlet and ends at upstream of the outlet. Therefore, for an actual annulus thousands of meters long, the scales of both the developing region and the dropping region are too short to be considered; hence, the cementing ECD mainly depends on the average frictional pressure loss in the steady region, which can be calculated as
Thus, the frictional pressure loss for the whole well annulus can be calculated as where and are the pressures at downstream of the inlet and upstream of the outlet, respectively.
Furthermore, if the hydraulic equivalent diameter of the well annulus varies as the well depth, the well annulus needs to be divided into a number of parts, and each part has the same hydraulic equivalent diameter. Then, the average frictional pressure loss for each part is calculated by using Eq. (15). At last, the frictional pressure loss for the whole well annulus is obtained by summing the frictional pressure loss for all the divided parts as where is the index of the part and is the total number of the divided parts.
4. Calculation of Cementing ECD for Narrow Annuli
As the aforementioned, the contraction on the radial size of annulus will increase the frictional pressure loss and triggers rapid variation of the cementing ECD which may disturb the stability of the whole pressure system in the downhole; hence, a narrow annulus is one of the main challenges for well cementing.
4.1. Annulus Geometry and Cementing Slurries
In this section, a 1000-meter-long annulus with a single radial size is considered to analyze the influences of the narrow annuli. Specifically, the inner diameter of the annulus is fixed as , while four different outer diameters of annuli are set as , 0.1270 m, 0.1320 m, and 0.1372 m; hence, the corresponding hydraulic equivalent diameter are 0.0076 m, 0.0127 m, 0.0178 m, and 0.0229 m. The first three cases belong to narrow annuli (), while the last case is used for comparison. Moreover, in order to further investigate the influences of fluid properties and rheological models on cementing ECD, three types of cementing slurries are considered, whose rheological behaviors are tested via experiment. Specifically, the rheometer was used, and its rotary speed was set as , , , , and ; thus, the corresponding rheological parameters can be obtained according to their readings. The Bingham and Power-Law models are applied to describe the rheological behaviors of the cementing slurries, and the obtained rheological parameters are listed in Table 2.
4.2. Calculation of Cementing ECD
Based on the established structure model, cementing slurries flowing in annuli are numerically simulated via FLUENT, and the obtained results by using both the Bingham model and the Power-Law model are listed in Table 3, which are further compared in Figure 4. Specifically, 12 study cases are simulated numerically for both the Bingham model and the Power-Law model.

According to the comparison between the local Reynolds number and the critical Reynolds number for turbulence , the cement slurries remain as laminar flow in all the test cases. As can be observed, the ECD calculated based on the Bingham model is higher than that based on the Power-Law model, especially for the annuli with narrow radial sizes; moreover, as the annulus becomes narrower, the ECD discrepancy between the two rheological models becomes bigger. However, when the diameter of the annulus is larger than the defined narrow size, namely, when , the ECD results obtained from different rheological models are in good agreement. Therefore, the cementing ECD is sensitive to the rheological model in narrow annuli, which should be taken into consideration during well-cementing design.
4.3. Comparison of the ECD Results Obtained from Five Different Methods
In addition to the MR method and the DNS method, a number of commercial software, including DB, LM, and PVI, are generally applied for practical well-cementing designs. In this subsection, all the five methods are utilized to calculate the cementing ECD, and the obtained results are compared in Table 4. In order to verify the accuracy of the MR method in calculating cementing ECD for narrow annuli, meanwhile to determine the adaptability of commercial software for different rheological models, the ECD results obtained from the DNS method are viewed as standard values, and the corresponding error analyses are thus conducted. According to the applied different rheological models, the analysis results are divided into two groups: Figure 5(a) is for the Bingham model and Figure 5(b) is for the Power-Law model.

(a)

(b)
As can be seen, the ECD calculated via the MR method is basically in accordance with that via the DNS method, and all their errors are limited within 3%; therefore, the accuracy of the theoretical MR method in calculating cementing ECD for narrow annuli is effectively verified.
For the adaptability of commercial software on rheological models, obvious discrepancies are explored. Specifically, for the Bingham model, PVI has the lowest calculation error of ECD, whose accuracy stays in the same level as the MR method. The ECD calculated from LM is considerably lower than that from the DNS method, and the average error is around -35%. For DB, error fluctuation can be observed; the calculated ECD values for cement slurry B are close to the results of LM, while for cementing slurries A and C, the majority of the calculated ECD values are closer to the results of PVI. Different from the Bingham model, for the Power-Law model, LM has the lowest calculation error of ECD; however, its error level is still close to 10%. The ECD calculated from PVI is around 20% higher than the DNS result. The calculation error of DB displays a wide fluctuation again, which varies from -15% to even above 40%. According to the above comparison, commercial software has the specific adaptability on the rheological model; specifically, PVI is more accurate for analyzing the Bingham model, while LM is better to calculate the Power-Law model.
5. Calculation of Cementing ECD for a Realistic Annulus with Varied Radial Sizes
5.1. Calculation of Cementing ECD via the DNS Method
By applying the same strategy as introduced in Section 4, a realistic annulus with varied radial sizes is investigated in this section. The parameters for the annulus geometry are listed in Table 5; its inner or outer diameter varies as the well depth; hence, the whole annulus can be divided into four sections, and its configuration is depicted in Figure 6.

After building the geometry model for the annulus with varied radial sizes, the next step is meshing. In order to secure the meshing quality, the geometry model cannot be meshed as a whole. Under such circumstances, the geometry model has to be sliced into a few parts primarily, and then, each part is meshed solely according to the specific grid size; after that, all the meshed parts are assembled together to form a complete structure model. In particular, to avoid the mismatching of the grids belonging to different parts during assembly, the slicing positions for this geometry model are set within the two middle sections rather than on the interfaces of sections, which are marked as the two dashed lines in Figure 6(a). Therefore, the geometry model in this case is sliced into three parts, and each part owns two sections; then, the three parts are meshed independently in Figure 6(b), and the assembled structure model is finally obtained as shown in Figure 6(c).
Subsequently, the assembled structure model is imported into FLUENT for numerical simulation. However, different from an annulus with a single radial size, the change of the flow speed in junction regions of sections can be clearly observed for the annulus with varied radial sizes (see Figure 7), while the variation of the flow speed changes the frictional pressure loss; hence, the total frictional pressure loss should be calculated as the sum of the for all the divided sections.

(a)

(b)
During the numerical simulation, three types of cement slurries listed in Table 2 are used again, and different flow rate is considered. The results of the cementing ECD obtained by using both the Bingham model and the Power-Law model are listed in Table 6. According to the comparison between the maximal local Reynolds number among the four sections and the minimal critical Reynolds number for turbulence , the cement slurry remains as laminar flow in all the test cases. As can be observed, for all the three types of cement slurries, when increasing the flow rate, the frictional pressure loss keeps increasing as well and, thus, the corresponding cementing ECD.
5.2. Comparison of the ECD Results Obtained from Five Different Methods
In addition to the DNS method introduced in Section 5.1, the four other methods, including the MR method, DB, LM, and PVI, are also applied to calculate the cementing ECD for the annulus with varied radial sizes, and the obtained results are compared in Table 7.
The same as Section 4.3, in order to compare the obtained ECD results from the five different methods, the results calculated by the DNS method are viewed as standard values, based on which, the corresponding error analyses of cementing ECD for the four other methods are conducted. According to the applied different rheological models, the analysis results are divided into two groups, as shown in Figure 8; the error analyses for the Bingham model and the Power-Law model are displayed in Figures 8(a) and 8(b), respectively.

(a)

(b)
As can be observed, the cementing ECD calculated via the MR method is still in accordance with the result of the DNS method; therefore, the accuracy of the theoretical MR method in analyzing the annulus with varied radial sizes is verified. Moreover, the discrepancies of commercial software on the calculation of cementing ECD are explored. Specifically, for the Bingham model, the ECD results from both LM and PVI are lower than that from the DNS method, and both of their errors distribute around -10%. The calculation errors of DB fluctuate sharply for different cement slurries. Specifically, the results for cement slurry A are higher than the DNS results; oppositely, the results for cement slurry B are lower than the DNS results, while for cementing slurry C, the ECD values are close to the DNS results. For the Power-Law model, LM has the lowest calculation error for cementing ECD, and its absolute error is less than 3%. Comparatively speaking, the errors for PVI can reach -10%, while those for DB still fluctuate widely, and its error range varies from -20% to over 10%.
According to the above comparisons, the error fluctuations of DB for different cement slurries are observed again; hence, it is not a wise choice for the ECD calculation in well-cementing designs. Once again, LM is verified to be the best for analyzing the Power-Law model, while PVI is more accurate for analyzing the Bingham model.
6. Conclusions
This work presented a detailed study about the calculation of cementing ECD. Five methods, including the MR method, DNS method, DB, LM, and PVI, were applied to calculate the cementing ECD for three types of cement slurries flowing in narrow annuli or annuli with varied radial sizes. Both the Bingham and Power-Law models were considered to describe the rheological properties of the studied cement slurries. Based on the comparisons of the results calculated by the five methods, three main conclusions are summarized:
Primarily, for narrow annuli, the cementing ECD calculated based on the Bingham model was higher than that using the Power-Law model; moreover, as the annulus became narrower, the ECD discrepancy between these two rheological models became larger. However, when the diameter of the annulus was larger than the defined size as a narrow annulus (), the ECD results obtained from different rheological models overlapped gradually. Therefore, in narrow annuli, the cementing ECD is sensitive to rheological models, which should be taken into consideration during well-cementing designs.
Secondly, for both the narrow annuli with a single radial size and the well annuli with varied radial sizes, the cementing ECD calculated by the MR method was basically in accordance with that by the DNS method. Therefore, the accuracy of the theoretical MR method in calculation of cementing ECD was effectively verified.
Thirdly, for the investigated three types of commercial software, LM is the best one for analyzing the Power-Law fluid, while PVI is more accurate for analyzing the Bingham fluid. However, DB displayed large fluctuations in the calculation of cementing ECD for different cement slurries; hence, it should be treated seriously during well-cementing design.
In addition, in this present work, only the basic model for the calculation of cementing ECD is discussed. However, in the practical operation of well cementing, the influences of temperature and pressure are significant, which will be taken into considerations further as the future work. Meanwhile, based on the developed temperature–pressure-modified model of the cementing ECD, the adaptability of the commercial software on different rheological models should also be further investigated to determine their advantages and disadvantages in the practical well-cementing designs.
Data Availability
All data, models, or code generated or used during the study are available from the corresponding author by request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This paper is supported by the National Natural Science Foundation of China (Grant No. 51904018).