Numerical modeling of motion of displacement interface in eccentric annulus during primary cementing

Cementing operations are of major importance to the long‐term zonal isolation and safety of oil and gas wells. As one of the well barriers, the cement sheath is key to the long‐term integrity of the wellbore. If the cement slurry does not completely fill the annulus between the casing and formation during the primary cementing displacement process, a potential leakage channel for subsurface fluids is created. In this study, the motion of the interface between the cement and other displaced fluids was examined during cementing displacement. By dividing an eccentric annulus into two 2D planes along the radial and azimuthal directions, a set of motion models of the interface were developed and solved numerically using a first‐order upwind finite‐difference method. The influences of modified density difference F*, Bingham numbers B1 and B2, plastic viscosity ratio m, eccentricity e*, and diffusion coefficient D* on the motion behavior of the displacement interface were determined. The results indicated that the interface of the eccentric annulus showed serious channeling and a physical mixing phenomenon. Some residual displaced fluid could remain at the casing‐cement and cement‐formation interfaces as a result of the improper design of the fluid parameters. Keeping F* > 0 and B2 < 5 beneficially removed residue from the displaced fluid on the wall, and increasing F* and B1 (or decreasing B2, m, e* and D*) was conducive to the smooth motion of the displacement interface. This work showed that the interface motion model could predict the displacement quality and has great potential for improving the design of the primary cementing parameters.

and gas industry, such as drilling and exploration, as well as in production, 3 and the failure of the wellbore integrity will provide potential leakage channels for subsurface oil, gas, and water, seriously affecting oil and gas production and the life of a well.
There are many factors affecting wellbore integrity. One of the most important is the primary cementing process. The cement 4 that is used to stabilize casings and prevent the flow of formation fluids during well construction often displaces other fluids 5 and hardens in the annular space between the casing and formation. 6 Therefore, many fluid leakage mechanisms associated with primary cementing have been studied, along with the suitability of the properties of cement materials, 7 and quality of the cement bonding at casing-cement and cement-formation interfaces. 8,9 In well cementing operations, effectively pumping the cement slurry into the annulus and improving its displacement efficiency are also very important to prevent fluid channeling. 3 During primary cementing, a sequence of fluids such as a spacer and cement slurry are injected into the annulus to displace the drilling mud and prepare the annulus for cement placement, 6 as shown in Figure 1(A). The interface between the cement and displaced fluids is usually unsatisfactory during fluid flow. 10,11 Further, the displacement in the annulus due to the fluid properties and external factors (such as the casing eccentricity, 12 pump rates, pumping time, and flow regimes) can leave displaced fluid in the annulus and affect the seal quality of the cement sheath. As shown in Figure  1(B), with an eccentric casing, the interface on the wide side of the annulus moves faster than that on the narrow side. In addition, a worse interfacial phenomenon occurs involving the channeling of the cement slurry into the displaced fluid, and residual displaced fluid may be found at the casing-cement and cement-formation interfaces Figure 1(D). 13,14 Moreover, the channeling phenomenon becomes increasingly serious over time Figure 1(C,E), which leads to poor cementing quality. 15 Based on the differences among the research methods, studies on the displacement interface of annular cementing can be divided into three categories: experimental, computational fluid dynamics (CFD) simulation, 16 and theoretical modeling. Although high-speed cameras 17 and conductivity probes 18 have been used with experimental equipment to record and monitor the displacement process intuitively, these have unfortunately failed to quantify the shape of the displacement interface. Moreover, full 3D models have been constructed, and the fluid flow interface has been tracked using the fluid volume method with CFD simulation software. [19][20][21] However, these simulations required significant amounts of time for their numerical calculations. In addition, several cement flow models have been studied, such as simple kinematic models and two-dimensional (2D) models of the Hele-Shaw cell type. The kinematic models represented the annulus by numerous sectors, with the mass transfer neglected between them. The effects of significant factors on the displacement were investigated using these models, which enabled the prediction of the displacement efficiency. 22,23 In 2D models, the annulus was discretized into computational cells along the wellbore axis and along the azimuthal direction. These models had outstanding advantages in studying the theory of displacement in highly deviated and horizontal wells. 24,25 Meanwhile, the variation law for an interface in a plane over time under different casing eccentricities, density differences, and rheological parameters was analyzed. This illustrated why the displaced fluid remained on the wall 13,26,27 and pointed out the possibility of keeping the displacement interface stable. These studies found that a quantitative analysis of a moving interface is the key to studying the displacement process. However, the shape of the displacement interface in an eccentric annulus varies with time in both the radial and azimuthal directions. 12,28 There have been few analyses of the factors that affect the motion smoothness of the displacement interface in two directions over time, and detailed research on the channeling phenomenon of fluid flow in an eccentric annulus is lacking.
The objective of the present study was to provide a scientific basis for optimizing the cementing parameters and avoiding subsurface fluid migration to other formations or the surface as a result of poor displacement efficiency. A 2D model of the motion of the displacement interface was established to describe the shape of the displacement interface quantitatively. Moreover, the concentration transport equation was used to characterise the physical mixing phenomenon at the interface. Further, these models were solved numerically by using a first-order upwind finite-difference method. Thus, the motion law of the displacement interface was studied during cementing in an eccentric annulus by changing a series of dimensionless parameters in the models. Meanwhile, two dimensionless parameters (displacement half-width R w and axial distance difference A d ) were introduced to represent the motion smoothness of the displacement interface. Finally, to verify the correctness of the model, three actual wells were compared with the mathematical models. The results showed a good match with the variable cement bond log (CBL) and variable density log (VDL).

| MATHEMATICAL MODEL
To model the viscous fluid flow in the annular space formed by two cylinders, the interface between displaced fluid and displacing fluid would be constructed, and the shear thinning behavior and physical mixing phenomenon of fluids at the interface would also be studied and analyzed by establishing the dimensionless momentum equations, constitutive equations, and convective transport equations. Based on the motion characteristics of the fluid flow and displacement in a vertical eccentric annulus, the displacement interface should be studied in both the radial and azimuthal directions of the annulus because of the influences of the casing eccentricity, fluid rheology, and so on. 6 In the paper, an eccentric annulus was divided into two 2D planes, including a radial plane and an azimuthal plane (as shown in Figure 2. Before establishing the mathematical model, the following conditions were assumed: (a) the displacement process is a steady laminar flow with a constant volume rate, (b) the wall of the annulus is rigid with no slip, (c) there is no chemical contact contamination at the displacement interface, and (d) the azimuthal shear stress gradient is ignored.
Based on the assumed conditions, the interface motion model in the radial plane should be built first. As shown in Figure 2(A,B), fluid 1 is the displacing fluid, fluid 2 is the displaced fluid, and both follow the Bingham plastic fluid model. In addition, it is obvious that the displacing fluid and displaced fluid occupy different regions along the r-axis, which are symmetric about the z-axis during the flow process. Therefore, the model only needs to consider the interface motion on the right side of the z-axis.
For two viscous fluids, the dimensionless momentum equations are written as follows: where the functions f * and F * are defined as follows: The dimensionless constitutive equations of the displacing fluid and displaced fluid can be described by the following: where the dimensionless parameters m and B k 14 are given by.
In the radial interface, the flow velocity is related to r * and t * along the z-axis direction, and the dimensionless motion equation of the interface 13 is given as follows: Equation (5) describes the interface shape over time in the radial direction of the annulus, but fails to show the azimuthal shape of the interface. Therefore, based on the geometric relationship between the half radial width h * and azimuthal angle θ * of the annulus, the interface model in the azimuthal plane can be established based on the radial interface model (as shown in Figure 2(C).
The dimensionless geometrical equation of the azimuthal plane is given by: At the same time, there will be a mixing zone at the interface between fluids 1 and 2, which is a convection-diffusion phenomenon caused by the different concentrations of the two fluids. Therefore, the concentration transport equation 29 is introduced to describe the physical mixing phenomenon. In this paper, the dimensionless convective transport equation that considers the radial and axial diffusion is as follows: In summary, the dimensionless boundary conditions can be defined as follows: The equations used to establish this mathematical model involve many parameters and variables. To simplify the form of these equations and reduce the number of variables, the equations are made dimensionless using the following scaling:

| NUMERICAL METHOD
The previously described models can be used to study the motion of the displacement interface between the displaced fluid and displacing fluid over time, in both the radial and azimuthal directions of the annulus. However, the interface model used in the solution process is complex, comprising the set of Equations (1)-(8) (see Appendix S1 for detailed solution of model). Furthermore, Equations (5) and (7) are . (3) the keys to the solution of the entire model, because these indicate the nonlinear characteristics. Therefore, the finite-difference method is used to solve the equation set numerically.
As can be seen from the set of dimensionless equations (Equations 1-7) that describe the interface motion behavior, r m * and c * are dimensionless parameters related to time and space. In this paper, a first-order upwind finite-difference method is used to discretize Equations (5) and (7), and nodes i, j, and n represent the radial, axial, and time positions, respectively. The discrete scheme is expressed as follows: The motion behaviors of the displacement interface in both the radial and azimuthal directions can be described (as shown in Figure 3) by solving the dimensionless equations for the interface motion and concentration diffusion based on the discrete scheme of Equation (10). In Figure 3(A,B), the upper and lower parts show the motion curves and concentration field at t * = 16 of the radial and azimuthal interfaces in the semi-annulus, respectively. The interface motion curves of Figure 3(A) were drawn for every 250 time steps, representing the interface shapes at those times. Furthermore, the distance from the dotted line to r * = 1 in the figure is the wall residual layer thickness of fluid 2 13 . In Figure 3(B), the displacing and displaced fluids are shown in red and blue, respectively, and the white line shows the interface curve at t * = 16.

| RESULTS AND DISCUSSION
During cementing displacement, improperly designed engineering parameters will lead to the retention of the displaced fluid on the wall, and the elongation of the displacement interface, which will affect the displacement quality. Therefore, it is necessary to conduct a detailed analysis and careful study of the key factors (such as F * , B k , m, e * , and D * ) affecting the motion of the interface to prevent residual displaced fluid and maintain a smooth interface motion.

| Wall residual layer thickness H *
If only fluid 1 flows in the radial plane Figure A3, when the wall shear stress of fluid 1 is equal to the yield point of fluid 2, a wall residual layer of fluid 2 with thickness H * is produced, which can be derived from Equation (A2). Based on the calculated results, it is concluded that H * is a dimensionless parameter related to B 1 , B 2 , and F * . The effects of different dimensionless parameters on H * are shown in Figure 4. From Equations (2) and (4), F * can be regarded as the density difference between fluids 1 and 2, and B k is used to characterise the dimensionless yield point of fluids 1 and 2, where a single factor analysis is performed to further understand these dimensionless parameters.
As shown in Figure 4(A-D), when F * is increased from −30 to +60, the area with H * = 0 (ie, the darkest green area) gradually increases. When F * is 60, H * is generally smaller in the range of Bingham number B k . In addition, it can be found that H * decreases with an increase in B 1 , whereas it increases with an increase in B 2 in each part of Figure 7. In particular, B 2 has a greater influence on H * compared with B 1 . Based on the above findings, combined with the expressions for F * and B k in Equations (2) and (4), the following conclusion can be drawn: increasing F * (ie, increasing the density difference between fluid 1 and fluid 2) is conducive not only to displacing the displaced fluid, but also to the formation of a zero-thickness wall residual layer under the buoyancy effect. Moreover, with B 1 increasing or B 2 decreasing, it is difficult for the displaced fluid to remain at the wall under a positive viscosity difference between fluid 1 and fluid 2. Therefore, to prevent the displaced fluid residue on the wall, F * should be kept as large as possible, and B 2 must be less than five during detailed engineering.

| Effect of modified density difference F *
As described, the displacement interface in an eccentric annulus moves continuously over time in both the radial and azimuthal directions, and its channeling behavior is the key to determining perfect displacement. To study the motion characteristics of the interface over time, the motion smoothness of the interface in both the radial and azimuthal directions is studied by introducing the dimensionless displacement halfwidth R w of the radial interface and the dimensionless axial distance difference A d between θ * = 0 and θ * = 1 of the azimuthal interface. Figure 5 shows the effect of F * on the motion of the displacement interface (see also Figure S1).
Over time, the front width of the interface decreases in the upper part of the figure, while the azimuthal displacement difference increases in the lower part, resulting in a serious channeling phenomenon, as shown in Figure 5(A). In addition, a physical mixing phenomenon occurs in the constant velocity region of the interface, where the velocity is the highest Figure 5(B). Changing F * from −30 to +60 can effectively alleviate the channeling phenomenon. As shown in Figures 6 and 7, the displacement halfwidth R w increases with F * , especially when F * = 60, and the smoothness of the interface is greatly improved. However, F * has little effect on axial distance difference A d . According to the definition of F * , some findings imply that increasing F * can increase the relative density difference between the displaced fluid and displacing fluid. The higher density of fluid 1 will produce a buoyancy effect on fluid 2, which will make it easier to be displaced and will not change the flow velocity in the central region of the annulus. Thus, the radial interface moves more smoothly and azimuthal interface has no obvious change.
In summary, F * has a particularly important effect on the motion of the displacement interface. A reasonable design for dimensionless parameter F * is beneficial not only to obtain a wall free of displaced fluid residue, but also to the smooth motion of the radial interface.

| Effect of Bingham number B k
Cementing fluids are non-Newtonian fluids, and their yield characteristics are some of the factors leading to the channeling phenomenon. To study the effects of the yield points of fluid 1 and fluid 2 on the channeling behavior, the influence of B k on the smoothness of the interface was analyzed. Figure 8 shows the effect of B k on the motion of the displacement interface (see also Figures S2 and S3). Figures 9 and 10 present the effects of B k on the displacement half-width R w and axial distance difference A d , respectively.
Based on these results, optimizing B k can reduce H * and improve the interface morphology. In addition, Figure 8 shows that B k has great influence on the shape of the displacement interface in both the radial and azimuthal directions. Moreover, this has significant implications for the smoothness of the displacement interface over time. Furthermore, R w increases with a change in B 1 from 5 to 20 and decreases with a change in B 2 from 5 to 50, but B 2 has little effect on R w when B 2 < 10.
In addition, Figure 10 shows that A d is not affected by B 1 and decreases with B 2 . A review of the relevant previous studies 23,30 showed that some conclusions were found, namely that increasing the yield point of fluid 1 or decreasing the yield point of fluid 2 will lead to a positive viscosity difference between the two fluids, which is conducive to the stable propagation of the radial interface. However, from the motion curves of interface (see also Figures S2 and S3), it can be found that the flow velocity at r * = 0 is independent of B 1 , and decreases with an increase in B 2 , resulting in the results shown in Figure 10.
In conclusion, increasing B 1 or decreasing B 2 will make the radial interface move more smoothly, but the motion of the azimuthal interface will be affected with decreasing B 2 , which calls for higher requirements when designing annular fluid rheological parameters.

| Effect of plastic viscosity ratio m
In addition to the yield point, the plastic viscosity of a fluid has a great influence on its flow characteristics in the Bingham plastic model. Figure 11 shows the effect of m on the motion of the displacement interface (see also Figure  S4). Figures 12 and 13 further show the effects of m on the displacement half-width R w and axial distance difference A d , respectively.
As shown in Figure 11, at m = 0.2, the displacement interface Figure 11(A) moves more slowly and smoothly. In addition, R w decreases and A d increases as m increases from 0.2 to 5, as illustrated in Figures 12 and 13. It is clear that m is a unique factor affecting the shape of the displacement interface independent of the wall residual layer thickness. Making the plastic viscosity of the displacing fluid larger than that of the displaced fluid helps to ensure a smooth motion of the displacement interface. In practice, m should be kept as far below one as possible. If conditions permit, it is possible to achieve a perfect displacement by minimizing m. 4.5 | Effect of dimensionless eccentricity e * Sections 4.2-4.4 clearly showed that the interface in the radial plane is affected by F * , B k , and m, leading to a change in the interface shape in the azimuthal plane. However, an eccentricity in the casing is the main reason for the change in the interface morphology in the azimuthal plane. Figures  14 and 15 show the effect of dimensionless eccentricity e * on the motion of the displacement interface (see also Figure S5).
It is apparent that A d increases with e * . When e * = 0.33, the axial displacement distance of the interface on the wide side of the annulus (θ * = 1) is 1.67 times that of the narrow side (θ * = 0) at t * = 16, while the difference is 6.75 times when e * = 1. Moreover, with a significant eccentricity in the casing, the fluid on the narrow side of the eccentric annulus might stop flowing altogether. Therefore, a casing centraliser should reasonably

| Effect of dimensionless diffusion coefficient D *
Because of the uneven velocity distribution of a viscous fluid in an annular shear flow, the behavior of material diffusion with the flow is considered to be a physical mixing phenomenon, which affects the concentration field distribution at the interface. Figure 16 shows the effect of dimensionless diffusion coefficient D * on the distribution of the concentration field. According to the concentration distribution results shown in Figure 16(A,B), it is not hard to find that the mixing zone decreases when the magnitude of D * decreases from 10 −3 to 10 −6 , and the concentration changes most obviously along the radial direction. No axial diffusion is expected to occur during the cementing operation because it will accelerate the mixing of fluids and then destroy the properties of a fluid such as the cementing slurry. However, radial diffusion can cause a mass exchange between the displacing fluid and wall residual layer, which is conducive to displacing the residual liquids. Therefore, it is very important to use the mixing law reasonably to further improve the displacement efficiency.

| Case study
Three wells were considered as practical cases to investigate the motion smoothness of the displacement interface and compare the cementing quality based on the mathematical model of this study. The cementing parameters of the three case wells are listed in Table 1. It can be seen that the three selected wells are all vertical wells, and their CBL/VDL logging results are shown in Figure 17. The results indicated that the cementing quality of cases 1 and 2 was generally poor, while that of case 3 was good. To further study these cases, the field experiment data were converted into the relevant dimensionless parameters, as listed in Table 2. In addition, a comparison of the displacement interfaces of three cases is shown in Figure 18. Figure 18 shows that the displacement half-width R w of case 3 is approximately constant along the z * -axis, and its interface moves faster and smoother than those of the other two cases, which explains why the cementing quality of case 3 is the best. A further analysis of this result combined with Table 2 indicates that the cementing parameter design of the three cases ensures that the displaced fluid will be free of residue on the wall. However, the channeling phenomenon still occurs because of the lack of detailed research on the interface motion. For example, in case 3, the interface moves more smoothly because it has the largest F * and B 1 values, and the smallest m values, which is consistent with the previously mentioned conclusions. Therefore, scientifically designing the cementing parameters can prevent the channeling phenomenon and greatly improve the displacement quality.

| CONCLUSIONS
In this study, the shape of the displacement interface in both the radial and azimuthal directions of an eccentric annulus was quantified by studying the motion law of the interface between two viscous fluids, which in turn was used to analyse the motion behavior of the interface and provide a scientific method for evaluating the displacement quality. In this process, the motion equations of the displacement interface and concentration transport equation were built and solved numerically using a first-order upwind finite-difference method. The displacement half-width R w and axial distance difference A d were introduced to evaluate the smoothness of the interface.