Regular Article
Heat dissipation of the Electrical Submersible Pump (ESP) installed in a subsea skid
^{1}
Graduate Program in Energy, Federal University of Espirito Santo, São Mateus, 29932540 Espirito Santo, Brazil
^{2}
Department of Engineering and Technology, Federal University of Espirito Santo, São Mateus, 29932540 Espirito Santo, Brazil
^{3}
Department of Industrial Technology, Federal University of Espirito Santo, Vitória, 29075910 Espirito Santo, Brazil
^{4}
Petrobras Research and Development Center Leopoldo Américo Miguez Mello (CENPES), 21941915 Rio de Janeiro, Brazil
^{*} Corresponding author: oldrich.romero@ufes.br
Received:
17
March
2019
Accepted:
30
January
2020
The recent development of Electrical Submersible Pump (ESP) in the skid, installed in the seabed downstream of the wellhead in an offshore oil production system, is an alternative to the conventional system with the set installed at the bottom of the producing well, facilitating interventions in case of failure. The pump is driven by an electric motor whose cooling must be efficient to ensure the continuity of its operation. The heat withdrawal is performed by the fluid produced. The purpose of this article is to understand the process of electric motor cooling to the singlephase and turbulent flow with convection heat transfer in an annular geometry, which represents the space formed between a capsule and the ESP in the Skid system motor. With this objective it is employed a Computational Fluid Dynamics (CFD) code to solve the governing equations of the turbulent heat transfer singlephase flow. The standard κε model with improved wall function (Enhanced Wall Treatment) is used to closure turbulence problem. This study considered flow rates range of 2200–4200 m^{3}/d (representing Reynolds numbers range of 27 000–133 000 approximately), Prandtl numbers 7–37, three configurations of different annular geometries, one concentric and two eccentric, together with the condition of the constant temperature on the motor surface (130 °C) and capsule (4 °C). The simulations are validated by comparing the Nusselt number in the developed region with the Gnielinski correlation. It is observed that if the constant heat flux condition were used, the motor temperature would have lower values at the beginning and larger at the end of the geometry. Therefore, the higher the Nusselt number, the greater the heat transfer, thus intensifying the cooling of the electric motor. In the eccentric geometry a momentum transfer from the lower to the upper annular region is observed, causing the Nusselt number present an angular variation. In eccentric geometries the flow develops in greater lengths, observing that the greater the eccentricity, the greater this length. Finally, for the ESP in the Skid system the use of an eccentric geometry is not adequate.
© J.R. Martins et al., published by IFP Energies nouvelles, 2020
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Conventional Electrical Submersible Pump (ESP) systems installed in deep waters (Fig. 1a) present difficulties when interventions are needed for replacement or repair, in addition to their life expectancy be relatively low. These characteristics lead to a significant downtime of production. The high costs associated with maintenance operations and the long average time for repair motivated the search for alternative solutions, to ensure the production process with reduction of the interventions.
Fig. 1 ESP in the skid. (a) Schematic representation (not in scale) of the offshore production system with ESP location and problem domain of study. (b) ESP in the skid detail. The pump module houses two ESPs placed in an “X” configuration (Tarcha et al., 2015). 
According to Tarcha et al. (2015), one of the proposals was the development of an ESP out of the well, installed in the seabed, to solve problems with recovery and reinstallation of equipment, without the need to remove the production column. This system is called ESP in the skid (Costa et al., 2013), represented in a simplified way by Figure 1a. The ESP in the skid leads the production directly to a pump seated on the seabed, which would facilitate future interventions in case of failure. In addition, this system uses two pumps in series with independent electric motors, allowing the continuity of production in case of failure of one of the pumps (Fig. 1b). On the other hand, the main disadvantage of this system is that its use is recommended to produce with a gas fraction lower than 35%. The production with high concentration of gas causes problems in the pumps reducing their average time between failures.
According to Takács (2009), one of the problems that occur in the ESP system is the high heating of the electric motor. This heating occurs during the process of transforming electrical energy into mechanical energy and is caused by the friction resulting from the high rotation of the shaft that drives the pump (Betonico et al., 2015). The process is more inefficient as the temperature increases. Therefore, an adequate cooling of the motor is required, which occurs mainly by the transfer of heat between the motor, fluid produced, and seawater.
The produced fluid, responsible for cooling the motor, flows into an annular region formed between the electric motor and an external capsule. Because the high flow rates and the reduced open flow area, the regime is turbulent. Thus, it is important to analyze which turbulence model should be used in the simulation to analyze such phenomenon. In the literature, the most commonly used turbulence models were the and SST families (Boersma and Breugem, 2011; Khalil et al., 2009; Mao et al., 2014).
The turbulent flow in annuli is associated with the phenomenon of heat exchange, which makes the analysis of the flow considerably more complex. In view of this, is given more emphasis in the prediction of the dimensionless parameters such as Nusselt number through empirical correlations. One of these correlations, valid for fully developed turbulent flow, is the Gnielinski correlation described in equation (1) (Skoczylas and Alhanati, 1998). This mathematical expression, where f is the friction factor, is applied for Reynolds number (Re) ranging in the interval 3000 < Re < 5 × 10^{6}, including the transition between the laminar and turbulent regime, and Prandtl number (Pr) 0.5 ≤ Pr ≤ 2000:(1)
Some works were performed in concentric annular spaces involving turbulence and heat exchange, such as Braga (1987), which analyzes the heat transfer and the pressure drop in turbulent flow through smooth and finned annular regions. The experiments were performed on a concentric heat exchanger. Water (hot fluid) flows through the inner tube, while air (cold fluid) flows through the annular region, in parallel currents.
Betonico et al. (2015) demonstrated a methodology to calculate the ESP motor temperature, considering the convection phenomenon, and the coupled behavior between motor, pump, and the production system. Noting that the motor temperature is increased when the rotation speed is also increased, the maximum temperature is located at the top of the motor and the minimum temperature at the base. It has also been shown that neglecting the development of the thermal boundary layer can result in a very high motor temperature than its actual value.
Therefore, the objective of the present work is to analyze the motor cooling process considering turbulence and convective heat transfer in a singlephase flow. The domain considered is the annular space between the electric motor and the wall of the capsule (Fig. 1a). The numerical solution of the problem is performed in the CFD code ANSYS Fluent v15.0. The turbulence model used is standard with enhanced wall function, due to its excellent results and low computational cost. The result is validated by comparing the Nusselt number from the simulations with the Nusselt obtained by the Gnielinski correlation. First, the electric motor cooling and the effect of the fluids properties are analyzed. Subsequently, a geometry with eccentricity is used and the electric motor cooling is analyzed again in order to verify if the eccentricity favors the cooling process.
2 Problem description
The problem domain, defined in Figure 1a, is the annular space between an inner pipe, which represents the electric motor, and one external pipe, which houses the motor/pump assembly, also referred as capsule. As shown in Figures 2a–2c, the dimensions of the two horizontal cylinders are: length L equals to L_{1} = 2 m and L_{2} = 4 m, external diameter D_{o} = 0.273 m, and internal diameter D_{i} = 0.1874 m. The length used was selected aiming at lower computational time, being this length shorter than the real system. In this geometry, the inner wall represents the electric motor surface that drives the centrifugal pump, the outer wall is the capsule and the annular region is the space through which the produced fluids (oil and water) flow before entering the pump.
Fig. 2 Schematic representation of the horizontal geometry with details of annular region and eccentricity, horizontal coordinate x is varying from 0 to L. (a) Full domain, (b) concentric tubes, (c) eccentric tubes. 
Keeping the dimensions constant, the position of the inner tube, in relation to the external one, was altered vertically (Fig. 2c) resulting in new configurations with different eccentricities. According to Caetano (1986), the eccentricity ϵ of two pipes is function of the distance C between the centers of each pipe of diameters D_{o} and D_{i}, according to equation (2):(2)
In Figure 2b, the geometry has eccentricity equal to 0, which is a concentric tubes configuration, and in Figure 2c a vertical eccentricity greater than 0. Three geometries, with eccentricities 0, 0.25, and 0.5 are used in this work.
Due to a previous analysis of the turbulence models, the simulations are performed considering the turbulence model with enhanced wall function ( EWT model). Three fluids, described by the Prandtl number, were used for the analysis of convective heat transfer phenomena together with the turbulence, according Table 1. Label C represents a twophase fluid, which is a mixture of fluids A and B, obtained by homogeneous modeling as described in Section 3.1.
Physical properties of the three fluids (A, B, and C) considered.
Fluids A, B, or C are injected into the inlet plane, represented in Figure 2a, with a temperature of 62.5 °C. This temperature was based on Teixeira et al. (2012), which studies the Pumping Unit system (Modulo de Bombeio [MOBO]). This information is used since the MOBO and ESP in the skid systems present similar conditions at the entrance. The external and internal walls, present the water temperature in the seabed (4 °C) and the motor operating temperature, respectively. Boundary conditions are described in Section 3.4. We chose to use the constant internal temperature condition because the heat flux of the electric motor was not available. The solution methodology is similar for both boundary conditions. The operating temperature of the electric motor was fixed at 130 °C, based on Tarcha et al. (2016).
Specific heat c_{p} is the energy required to raise the temperature of a unit mass of a fluid by one degree in a specified way. This property is a measure of a material’s ability to store thermal energy. The thermal conductivity of a fluid k_{f} is a measure of the ability of the fluid to conduct heat. A high value indicates that the fluid is a good heat conductor, and a low value indicates that the fluid is a poor heat conductor or insulator. The product ρc_{p} is the heat capacity, and represents how much energy a material stores per unit volume (Cengel, 2002).
3 Mathematical modeling
A relevant concept to deal with problems in annular regions is the hydraulic diameter, D_{h}, obtained by:(3)where D_{o} is the outer diameter and D_{i} is the inner diameter (Figs. 2b and 2c). In this work D_{h} = 0.0856 m.
Three important physical dimensionless parameters are used and briefly described in that follows.
(i) The Reynolds number (Re), for annular space is:(4)where ρ is the fluid specific mass, V is the mean fluid velocity and μ is the fluid viscosity.
(ii) The Prandtl number (Pr), which describe the relative thickness of the velocity (δ) and the thermal (δ_{T}) boundary layers, is defined as the ratio of momentum diffusivity (υ) to thermal diffusivity (α), , or,(5)where c_{p} is the fluid specific heat and k_{f} is the fluid thermal conductivity. Small values of the Prandtl number, Pr ≪ 1, means the thermal diffusivity dominates (and scaling between boundary layers is δ_{T} ≫ δ, δ_{T} develop faster than δ), is the case of liquid metals. Whereas with large values, Pr ≫ 1, which is the case of water and oils, the momentum diffusivity dominates the behavior, thickness of δ is consistently much greater than δ_{T} (δ_{T} ≪ δ). (Bejan, 2013). For example, the typical value for liquid mercury, which is about 0.025, indicates that the heat conduction is more significant compared to convection, so thermal diffusivity is dominant. When Pr is small, it means that the heat diffuses quickly compared to the velocity. Prandtl number for water at 18 °C it is around 7.56, which means that the thermal diffusivity is not dominant.
(iii) The Nusselt number (Nu), or dimensionless heat transfer coefficient, defined by:(6)where h is the convective heat transfer coefficient and L_{c} is a characteristic dimension length of heat transfer. The Nusselt number is equivalent to the dimensionless temperature gradient at the surface and represents the enhancement of heat transfer through a fluid layer as a result of convection relative to conduction across the same fluid layer. The larger the Nusselt number, the more effective the convection. A Nusselt number of Nu = 1 for a fluid layer represents heat transfer across the layer by pure conduction. For turbulent flow, the Nusselt number is usually a function of the Reynolds number and the Prandtl number.
Convection is the mode of energy transfer between the motor solid surface and the adjacent produced liquid that is in motion, and it involves the combined effects of conduction and fluid motion. The faster the fluid motion, the greater the convection heat transfer.
3.1 Homogeneous modeling of twophase flow
The Basic Sediment and Water (BSW) used by Betonico et al. (2015), Table 1 is the mass fraction of water present in the fluid produced and can also be obtained by equation (7):(7)where Q_{w} is the water flow and Q_{o} is the oil flow.
Through BSW it is possible to calculate the parameters of the mixture (fluid C), such as density (ρ_{m}), viscosity (μ_{m}), specific heat (c_{pm}), and thermal conductivity (k_{m}), which are, respectively, represented by equations (8)–(11)(8)(9)(10)(11)being ρ_{o} and ρ_{w} the densities of the oil and water, μ_{o} and μ_{w} the viscosities of the oil and water, c_{po} and c_{pw} the specific heat of the oil and water and, k_{o} and k_{w} the thermal conductivities of the oil and water, respectively.
This whole procedure is such that the two immiscible fluids (oil and water) are represented as a homogeneous singlephase fluid. All these parameters are precalculated and entered as input values.
3.2 Isothermal turbulent flow
The isothermal turbulent flow is modeled using the continuity equation (12), the momentum conservation equation (13) and twoequation turbulence model with enhanced wall function, equations (14) and (15):(12)where is the velocity vector.(13)where = 9.8 m/s^{2}, is the stress tensor and p is the pressure (Pereira et al., 2017).
The κ − ε model is based on the concept of eddy viscosity, μ_{t}. Thus, this model solves two different transport equations: turbulent kinetic energy, κ, equation (14), and turbulent kinetic energy dissipation rate, ε, equation (15) (ANSYS, 2017):(14)(15)where P_{κ} is the turbulence production due to viscous forces; C_{s1}, C_{s2,} C_{μ}, σ_{κ}, and σ_{ε} are constants with values: 1.44, 1.92, 0.09, 1.0, e 1.3.
The turbulent viscosity μ_{t} is a property of the turbulent flow and is defined in the equation (16):(16)where c_{μ} is a constant.
For a better wall treatment, in this work the Enhanced Wall Treatment (EWT) was used, which combines a twolayer model (a region affected by viscosity and a fully turbulent one) with the socalled enhanced wall functions. This requires the mesh to be sufficiently refined to solve the viscous sublayer (ANSYS, 2017).
3.3 Nonisothermal turbulent flow
Since this scenario involves a convective heat exchange, the energy equation defined by equation (17), together with equations (12)–(15), are solved together:(17)H is the specific enthalpy, which is a function of pressure and temperature T. The enthalpy is described thermodynamically by the relation below:(18)where T_{ref} is the reference temperature.
3.3.1 Nusselt number evaluation
The heat flux q_{w} on the wall of the annular space is modeled from scalable wall function approaches, or the automatic wall treatment as follow (ANSYS, 2017):(19)(20)(21)(22)(23)(24)where u* is the scalar velocity used in the logarithmic region; T^{+} is the dimensionless temperature, which represents the temperature profile through the viscous sublayer and the logarithmic region; K is the von Karman constant; T_{w} is the wall temperature and T_{f} is the temperature of the fluid in the region furthest from the wall.
Given the heat flux in the wall q_{w}, calculated based in equation (19), the convection heat transfer coefficient h is found for a fixed temperature boundary condition T_{w}, for turbulent regime as follow:(25)where T_{wn} is a mean fluid temperature near the wall. In other words, the convection heat transfer coefficient is defined as the rate of heat transfer between the solid surface and a fluid produced per unit surface area per unit temperature difference.
With h calculated, the Nusselt number (Nu) in the desired region is obtained by the relation presented in equation (6).
3.4 Boundary conditions
The boundary conditions that delimit the domain of the problem are represented by letters in Figure 3 and we have:

Inlet – fluid with temperature of 62.5 °C and prescribed uniform velocity which depends of each analyzed condition, Table 2.
Table 2Main parameters of cases studied.
Outlet – prescribed pressure equals to 0 Pa.
Outer wall – constant temperature of 4 °C with nonslip condition and impenetrability.
Inner wall – temperature of 130 °C with nonslip condition and impenetrability.
Fig. 3 Boundary conditions of the problem domain, with x varying for 0 to L. 
Four operating flow rates of the ESP in the skid system were considered, being 4200, 3400, 2600, and 2200 m^{3}/d (Tarcha et al., 2016), as defined in Table 2.
4 Numerical method
The governing equations are solved numerically through the Finite Volume Method (FVM) available in ANSYS Fluent software v15.0. The computer used in the simulations has the following configurations: Intel (R) Core i75500U, CPU @ 2.40GHz, 8.0 GB of RAM and 64bit Windows operating system.
4.1 Discretization algorithms
The discretization of the governing equations used the FVM with colocated arrangement of the variables (velocities, pressure, turbulent kinetic energy, rate of dissipation of turbulent kinetic energy, temperature, and enthalpy). With this, the governing equations are transformed into algebraic equations. Coupling algorithms are required to obtain the solution. The SemiImplicit Method for PressureLinked (SIMPLE) is used. The Second Order and Second Order Upwind schemes are employed in the interpolation of pressure and velocity, respectively. For the evaluation of the gradients, the Least Squares Cell Based method is used. Finally, the First Order Upwind scheme is selected to interpolate the turbulent kinetic energy and the dissipation rate of the turbulent kinetic energy, and the Second Order Upwind scheme is employed in energy interpolation.
4.2 Mesh independence test
The mesh was refined according to the parameter y^{+} defined by equation (26), which is used for turbulent flow in regions close to the wall:(26)where u_{τ} is the friction velocity and y the normal distance from wall.
A value of y^{+} < 10 characterize the viscous sublayer of the turbulent boundary layer, where the effects of viscosity and molecular diffusion are predominant and the flow is almost laminar (Bejan, 2013).
Three refinement levels were used, as detailed in Table 3, mesh 1 (fine), mesh 2 (intermediate), and mesh 3 (coarse). Table 3 also shows the value of y^{+} for each mesh, because the lower its value, the better the capture of the phenomena near the wall. This feature of the meshes is reinforced with the high number of layers near the wall.
Characteristics of the three meshes tested.
The mesh test is performed using fluid A (Tab. 1), considering the turbulent flow condition with the Reynolds number of 133 758 and using the concentric geometry. The Nusselt number variation on the inner wall of the pipe along the dimensionless length x/D_{h}, from the simulation, with the Nusselt number obtained by the Gnielinski correlation, equation (1), is valid for x/D_{h} ≥ 10. The results are shown in Figures 4 and 5.
Fig. 4 Nusselt number for each mesh along the annular inner wall. The continuous horizontal line is a result of Gnielinski’s correlation. 
Fig. 5 Nusselt number for each mesh in the fully developed region, position x/D_{h} = 20 near the exit plane of the domain, and analytical response. 
In Figure 4 the numerical results show high values of Nusselt number at the inlet plane, they decrease markedly as the position in the pipeline increases until it becomes practically constant in x/D_{h} = 8 and remaining unchanged until the output in (x = L)/D_{h} = 23.36. This sets up a fully developed thermally flow. These values in the constant parts are very similar to the value obtained by the Gnielinski correlation.
A comparison of the Nusselt number from empirical correlation with the numerical approach for each mesh is presented in Figure 5 at a fixed position in the developed region (x/D_{h} = 20). The coarse (mesh 1) and intermediate meshes underestimate its value (with errors equal to −4.08% and −2.31%, respectively), while the fine mesh overestimates (error equal to +1.77%). As we can see, the finest mesh is more accurate. Incorporating the simulation time in the analysis, we have that the intermediate mesh (mesh 2) requires 2 h of simulation while the fine mesh (mesh 3), 6 h. As a result, the intermediate mesh with 4 104 000 elements (whose crosssectional view is shown in Fig. 6) was selected to represent the problem domain in a discrete world and to run all the simulations in this work.
Fig. 6 Cross section of annular space discretized using mesh 2 with 4 104 000 elements which is selected for all simulations. 
After incorporating the boundary conditions in the selected mesh, the simulation was executed in steady state, with the tolerance of 10^{−4}. The results are post processed in ANSYS CFDPost v15.0.
5 Results
5.1 Reynolds number effect on Nusselt number
The effect of fluid flow on heat transfer is analyzed in this section, considering case 1 with four Reynolds numbers, as detailed in Table 2. This process is observed by plotting the Nusselt number on the inner wall (electric motor surface) along the dimensionless length of the pipe, as presented by Figure 7.
Fig. 7 Nusselt number variation along the inner wall of the annular region considering case 1. 
The trend of the Nusselt number is similar for all situations considered here. Near the inlet plane the Nusselt number has higher values (Nu ~ 2400) and decreases abruptly until it reaches a constant behavior (Nu ~ 800 at Re = 133 758). For the four situations the flow can be considered fully developed when x/D_{h} ≥ 10, in this condition the Nusselt number tends to remain constant for the remainder of the geometry when x = L or L/D_{h} = 23.36. In other words, the heat exchange from the motor wall to the fluid produced (motor cooling process) is predominantly convection in the short stretch near the inlet.
The high Nusselt numbers found highlight that the heat transfer always occurs from the system with higher temperature (motor) to the lower temperature (fluid). The higher the Nusselt number, the higher the heat transfer will be, resulting in motor cooling. It can also be observed that the higher the Reynolds number, the higher the Nusselt number is. The convection heat transfer is directly related to the turbulence. The Prandtl number is constant and equal to 7 in this case.
In order to confirm these results, the Nusselt number of the fully developed region was compared with the result found by using the Gnielinski correlation, as represented by Figure 8 and Table 4. Both approaches give similar answers. Differences based on correlation did not exceed 5%. With this, it can be stressed that the present numerical approach adequately represents the behavior of the heat transfer phenomenon by convection and its magnitudes.
Fig. 8 Nusselt number on the annular inner wall of the fully developed region, using case 1 and comparing with Gnielisnki’s correlation. 
Difference of the Nusselt number (Nu) on the fully developed region using case 1 based on the Gnielinski correlation.
5.2 Prandtl number effect on Nusselt number
The effect of fluid transport properties on heat transfer is analyzed in this section, considering cases 1–3 (described in Tab. 2). Three values of Prandtl number (Pr) were used to analyze how much the increase of this parameter can influence the heat transfer by convection. This problem is analyzed calculating the Nusselt number on the annular inner wall, at position x/D_{h} = 20, and for the four Reynolds numbers used in each case, 12 in total. The results are shown in Figure 9.
Fig. 9 Nusselt number on the annular inner wall, which represents the motor surface, in the fully developed region varying with the Reynolds number for three fluids. 
It should be noted that the predominance of the inertial force on the viscous force, i.e. an increase in Reynolds number, causes the Nusselt number to increase as well. However, fluids with high Prandtl numbers present a faster response of the Nusselt number increasing to the Reynolds number variation. As the Prandtl number characterizes the relationship between the molecular diffusivity of momentum and molecular diffusivity of heat (Cengel, 2002), it is known that the lower Prandtl number has a faster heat diffusion. This graphical behavior agrees with the scaling results of Bejan (2013): Nu ~ Re^{1/2}Pr^{1/3}, valid for Pr ≫ 1.
When the curves are extrapolated it is possible to analyze that for a same Reynolds number the fluid with the highest number of Prandtl, Pr = 37, tends to have a greater Nusselt number than the others. The Reynolds number of 40 000 is taken as an example, the Nusselt number of the fluid with Pr = 13 is estimated to be close to 350, while that of the fluid with Pr = 37 is approximately 500. For the Reynolds number equal to 70 000, it is observed that the Nusselt number of the fluid with Pr = 7 is approximately 450, while that of the fluid with Pr = 13 is Nu = 600. Therefore, for a same Reynolds number, the higher the Prandtl number of the fluid, the higher its Nusselt number, and the more intense the heat transfer by convection, compared with conduction, between the electric motor and the produced fluid. High Prandtl number indicates that the heat conduction is less significant compared to convection, so momentum diffusivity is dominant (Bejan, 2013; Cengel, 2002).
5.3 Effect of eccentricity on heat transfer
The effect of the eccentricity on the heat transfer is analyzed. The eccentricity is important for situations that occur, for example, when there are no centralizers in the motor/capsule space. Two levels of vertical eccentricity were selected: 0.25 and 0.5.
In order to visualize the effect of eccentricity, four lines were used along the inner pipe, highlighted in front and lateral views of Figure 10 through points 1–4. These lines are listed clockwise, starting at the upper end. Point 1 (θ = 0°) has a larger opening for the flow, points 2 (θ = 90°) and 4 (θ = 270°) have an intermediate opening and point 3 (θ = 180°) presents the smallest open area.
Fig. 10 Numbers 1–4 represent the positions of the lines along the inner pipe surface where Nusselt number is plotted in the following figures. 
5.3.1 Eccentricity of 0.25
In this numerical experiment the case 4 (Tab. 2) was selected. The Nusselt number along the four lines defined in Figure 10 are presented in Figure 11 for Reynolds number and Prandtl number equals to 51 916 and 37, respectively. After reaching a value of Nu ~ 1855 (not showed in the figure), the Nusselt number decreases abruptly and tends to stabilize slowly, but in x/D_{h} = 14 the effect of the eccentricity begins to become important which occurs until the end of the domain in x/D_{h} = 23.36. The variation at line 2 (θ = 90°) and line 4 (θ = 270°) are similar due to the symmetry of the eccentric geometry about the vertical axis and is almost not affected because the open flow area is similar to concentric tubes. The responses at positions 1 (θ = 0°) and 3 (θ = 180°) show an opposite behavior: at position 1 the Nusselt number increases (good cooling process), while at line 3 decreases (bad cooling process). The flow apparently does not develop completely on any of the four selected lines, which is related to the eccentricity effect on the flow. It is necessary to increase the flow domain and reproduce the numerical experiments in order to confirm when the flow is fully developed.
Fig. 11 Nusselt number with Pr = 37 and Re = 51 916 on the four lines (Fig. 10) of geometry with vertical eccentricity of 0.25. The concentric response is also included. 
It is possible to conclude that the influence of eccentricity is important, in terms of cooling process of the electric motor, after x/D_{h} ~ 14, where begins region II in Figure 11, and extend all along this region until the end of the domain. Before that, the eccentric response is similar to the obtained with concentric tubes (region I in Fig. 11). The annular space with the bigger open flow area (line 1) favors the motor cooling. In the other three positions (line 2, line 3, and line 4), the heat exchange is lower than the concentric geometry.
To better understanding the influence of the Reynolds number (Re) on the heat transfer in each of the three positions along the eccentric annular space, four values of Reynolds number are selected. The results are presented in Figure 12a at line 1, Figure 12b at line 2, and Figure 12c at line 3. The maximal value of Nusselt number (Nu_{max} ~ 1855) is not showed in the vertical axes.
Fig. 12 Nusselt number along lines 1–3 (represented in Fig. 10), with Pr = 37 and four turbulent Reynolds numbers in a geometry with vertical eccentricity of 0.25. (a) Nusselt number along line 1 (θ = 0°). (b) Nusselt number along line 2 (θ = 90°). (c) Nusselt number along line 3 (θ = 180°). 
Clearly, the predominance of inertial forces improves the heat exchange along the length of the domain, and it is more intense in larger areas open to flow as the represented by line 1 (Fig. 12a). At the same Reynolds number, the increase or decrease of the heat transfer intensity begins at the same position and the fluid flow development occurs close to the end of the geometry (Figs. 12a and 12c). In Figure 12b the Nusselt number at position 2 (θ = 90°) has a low variation remaining practically constant throughout the pipe.
This behavior of the Nusselt number does not favor uniform cooling of the electric motor. Because the bottom region with less open space for the flow presents smaller Nusselt numbers, so the heat transfer in that region is lower, resulting in higher surface temperatures. On the other hand, the upper region with higher Nusselt numbers has lower temperatures. With this, the electric motor has an overheating in its bottom region, affecting its efficiency and putting at risk its continuing operation.
To analyze the development of the flow, the geometry had its length doubled from L_{1} = 2 m (L_{1}/D_{h} = 23.36) to L_{2} = 4 m (L_{2}/D_{h} = 46.72). In this case is selected Re = 51 916 and the results of the Nusselt number for this new geometry is represented in Figure 13, where Nu_{max} ~ 1855. It is evident that the flow is fully developed when x/D_{h} ≥ 40 at the end of region II, four times greater than that found in the literature for concentric geometries. Therefore, it is possible to verify that in an annular geometry the fully developed flow entry length becomes function of the eccentricity. Region III represents the fully developed region.
Fig. 13 Nusselt number along the fluid flow direction on the increased length geometry (from L_{1}/D_{h} = 23.36 to L_{2}/D_{h} = 46.72) with bottom vertical eccentricity 0.25, Pr = 37 and Re = 51 916. 
Through Figures 14 and 15 it is possible to explain the behavior of the Nusselt number for the geometry with vertical eccentricity of 0.25. The velocity intensity behavior in the bottom and top of annulus and the velocity intensity perpendicular to the flow in the sides of the inner cylinder, along the geometry for the four Reynolds numbers. In Figure 14 (depicts a vertical plane that crosses line 1 and line 3 of Fig. 10), the bottom annular region begins to decrease axial velocity from a certain point, and from there the upper annular region begins to increase its axial velocity. This is explained since the bottom region has a smaller area for the flow, so the boundary layers develop faster than upper region and the fluid tends to flow upward where the flow is not yet developed. This phenomenon becomes even more evident in Figure 15 (depicts a horizontal plane that crosses line 2 and line 4 of Fig. 10) through the behavior of vertical velocity. Positive values mean upward flow.
Fig. 14 Axial component of the velocity (m/s) in the bottom and top of annulus in the geometry with vertical eccentricity of 0.25 for the fluid with Pr = 37 and Reynolds numbers equal to: (a) 51 916; (b) 42 027; (c) 32 138; and (d) 27 193. Fluid flow is from left to right. 
Fig. 15 Vertical component of the velocity (m/s) in the geometry with vertical eccentricity of 0.25, for the fluid with Pr = 37 and Reynolds numbers equal to: (a) 51 916; (b) 42 027; (c) 32 138; and (d) 27 193. Fluid flow is from left to right. 
It is observed in Figure 15 that there is a momentum transfer from the bottom region of the geometry (region of smaller area to the flow) to the upper region (region of greater area for flow). It is also shown that the point where this transfer takes place and the intensity of the transfer is different for each Reynolds number. Since the lower the Reynolds number, the faster the flow develops in the bottom region and the transfer of momentum occurs. Thus, it is explained why the flow with different Reynolds numbers have different intensities in the Nusselt number variation and the positioning of this variation is different for each Reynolds.
5.3.2 Eccentricity of 0.50
For this condition, which is case 5 of Table 2, only three lines were used along the internal pipe arranged in the same way as the lines of the previous section, highlighted in Figure 10. Emphasizing that from 1 to 2 and to 3 the open space for the flow is reduced due to the vertical eccentricity, and that the geometry has symmetry around the vertical plane that crosses both cylinders axis.
First, a comparison is made between the Nusselt numbers found in the cases of vertical eccentricities of 0.25 and 0.5 considering the Pr = 13 and Re = 77 693 (Fig. 16). The Nusselt number has different behavior when the eccentricity is modified, showing that the Nusselt number for eccentric geometries is dependent on the eccentricity (ϵ).
Fig. 16 Nusselt number in positions 1–3 (of Fig. 10) along the domain with eccentricity 0.25 (continuous blue line) and eccentricity 0.5 (dashed black line). Pr = 13 and Re = 77 693. (a) Nusselt number along line 1. (b) Nusselt number along line 2. (c) Nusselt number along line 3. 
In Figure 16a, it is noted that the Nusselt number at line 1 of the geometry with ϵ = 0.25 shows a positive variation at a shorter entry length, x/D_{h} > 12, while the Nusselt number of the geometry with ϵ = 0.50 starts to vary when x/D_{h} > 16. However, at the end of the geometry the Nusselt number of the geometry with ϵ = 0.50 presents a higher intensity variation for the same length, starting at x/D_{h} > 21. In Figure 16b, it is observed that the Nusselt number in line 2 of the geometry with ϵ = 0.25 presents a small negative variation, almost imperceptible, however the Nusselt number of the geometry with ϵ = 0.50 has a negative variation of greater intensity from x/D_{h} > 18. Finally, Figure 16c shows that the Nusselt number at line 3 of the geometry with ϵ = 0.25 shows a smaller negative variation for a larger input length x/D_{h} > 13, while the number of Nusselt of the geometry with ϵ = 0.50 shows a greater negative variation from x/D_{h} > 12.
The largest variations of the Nusselt number presented along lines 2 and 3 for the geometry with ϵ = 0.50 occur due to the decrease of the open area to the flow at the bottom when compared to the geometry with ϵ = 0.25. Therefore, it is emphasized that the greater the eccentricity, the greater the Nusselt number variation because the momentum transfer from the bottom to the upper is more intense.
For this case, the behavior of the Nusselt number at each of three lines, 1 (θ = 0°), 2 (θ = 90°), and 3 (θ = 180°), is also analyzed at four Reynolds numbers, allowing a better understanding of the behavior of the heat transfer phenomenon and the analysis of the angular variation of the Nusselt number in this geometry. The results are represented in Figures 17–19.
Fig. 17 Nusselt number with Pr = 13 for four turbulent Reynolds numbers at line 1 (θ = 0°) of geometry with vertical eccentricity of 0.5. 
Fig. 18 Nusselt number with Pr = 13 for four turbulent Reynolds numbers at line 2 (θ = 90°) of geometry with vertical eccentricity of 0.5. 
Fig. 19 Nusselt number with Pr = 13 for four turbulent Reynolds numbers at line 3 (θ = 180°) of geometry with vertical eccentricity of 0.5. 
The Nusselt number starts to vary at x/D_{h} = 16 in line 1, at x/D_{h} = 18 in line 2 and at x/D_{h} = 12 in line 3. The Nusselt number variation in each position is more intense as higher is the Reynolds number. It is an expected behavior once the turbulence becomes more and more dominant.
On the other hand, it can be verified that the phenomenon responsible for the variation of the Nusselt number is the same highlighted in the previous section, momentum transfer from the bottom to the upper annular region. In this case, because the eccentricity is greater the momentum transfer is also greater. Therefore, when the eccentricity is increased the thermal boundary layer (δ_{T}) develops in regions with longer entry lengths.
It is evident that the increased eccentricity is detrimental to heat exchange in positions where the flow open area is reduced, lines 2 and 3. In contrast, in the uppermost position, line 1, where the area increase, the efficiency of heat exchange is improved. This is independent of the Reynolds numbers used in this work. As a consequence, the overheating of the electric motor in the bottom region could be more intense, putting the operation of the ESP module at risk.
Through Figures 20 and 21 it is possible to observe better the phenomenon responsible for this variation of the Nusselt number to the geometry with 0.5 eccentricity using the contour plots of axial and radial velocity components.
Fig. 20 Axial velocity (m/s) in geometry with eccentricity 0.5 for the fluid with Pr = 13 and Reynolds numbers equal to: (a) 77 693; (b) 62 894; (c) 48 095; and (d) 40 696. 
Fig. 21 Vertical velocity (m/s) in geometry with eccentricity 0.5 for the fluid with Pr = 13 and Reynolds numbers equal to: (a) 77 693; (b) 62 894; (c) 48 095; and (d) 40 696. 
In Figure 20 the bottom annular region begins to decrease its maximum axial velocity and the upper annular region begins to increase its axial velocity. This was explained in the previous section through the development of boundary layers. The momentum transfer is further evident in Figure 21 by the behavior of the vertical velocity in the horizontal plane. Positive values mean upward flow.
With increasing eccentricity, the bottom region has a reduction in its open area and higher resistance to flow, resulting in shorter onset of the vertical momentum transfer. All of this has a consequence which is a nonuniform and low efficient heat transfer from electric motor to surroundings.
6 Conclusion
Based on the presented results, the following can be concluded:
Through the Nusselt number curve obtained applying constant temperature at the motor surface, it is possible to analyze how the temperature behavior on the motor surface would be when the constant heat flux condition will be available. Because with high Nusselt number, greater the heat transfer is, which in turns intensify the cooling of the electric motor. Thus, when the constant heat flux condition is used, the motor surface exhibits lower temperatures at the beginning of the annular geometry and higher temperatures at the end of the geometry.
In an eccentric annular geometry, there is a momentum transfer from the bottom annular region (small open flow area) to the upper region (high open flow area). This phenomenon is due to shorter flow entry length at bottom than at upper region.
For geometries with high eccentricities (such as geometry with 0.50 eccentricity analyzed in the work) there is a considerable reduction in the area open to the flow at the bottom region, thus increasing the shear in that region. It increases the resistance to the flow at the bottom region, resulting in a momentum transfer from the beginning of the geometry.
For the singlephase fluid flow flowing in an eccentric geometry, the Nusselt number is dependent on the dimensionless spatial coordinate in the flow direction (x*), Reynolds (Re) and Prandtl (Pr) numbers, eccentricity (ϵ) and angle (θ), and can be written as follow:
It is evident that for eccentric annular geometries the flow develops at longer entry lengths than those demonstrated in the literature. Note that the greater the eccentricity, the longer the length required for the flow to reach the fully developed condition.
Finally, for ESP in the skid system an eccentric geometry is not recommended, since this does not favor the uniform heat transfer between the motor and the surrounding fluids. And does not favor the uniform cooling of the electric motor, since the bottom region with less open space for the flow presents higher temperatures, affecting its efficiency and putting its operation at risk.
This work is under progress to include multiphase effect on heat transfer.
Acknowledgments
The authors would like to thank CENPES – PETROBRAS for funding this research. Also thank the granted permission to use the Petroleum Engineering Numerical Simulation Laboratory (Labsim) installations at UFES and the SPE UFES Student Chapter for granting the access of OnePetro platform.
References
 ANSYS (2017) ANSYS Theory Guide, Canonsburg, PA. Available at https://ansyshelp.ansys.com/. [Google Scholar]
 Bejan A. (2013) Convection heat transfer, 4th edn., John Wiley & Sons, Inc., Hoboken, NJ. [CrossRef] [Google Scholar]
 Betonico G.C., Bannwart A.C., Ganzarolli M.M. (2015) Determination of the temperature distribution of ESP motors under variable conditions of flow rate and loading, J. Pet. Sci. Eng. 129, 110–120. [Google Scholar]
 Boersma B.J., Breugem W.P. (2011) Numerical simulation of turbulent flow in concentric annuli, Flow Turbul. Combust. 86, 113–127. [Google Scholar]
 Braga C.V.M. (1987) Análise Termohidráulica de Seções Anulares Lisas e Aletadas (in portuguesse), Tese (Doutorado), Pontificia Universidade Católica do Rio de Janeiro, Brazil, 339 p. [Google Scholar]
 Caetano E.F. (1986) Upward twophase flow through an annulus, Ph.D. Dissertation, The University of Tulsa, Tulsa, OK, 217 p. [Google Scholar]
 Cengel Y.A. (2002) Heat transfer: A practical approach, 2nd edn., McgrawHill, New York, NY, 896 p. [Google Scholar]
 Costa B.M.C., Oliveira P.S., Roberto M.A.R. (2013) Mudline ESP: Electrical submersible pump installed in a subsea skid, in: Offshore Technology Conference, 6–9 May, Houston, TX, USA. doi: 10.4043/24201MS. [Google Scholar]
 Khalil M.F., Kassab F., Adam I.G., Samaha M.A. (2009) Turbulent flow around single concentric long capsule in a pipe, Appl. Math. Modell. 34, 8, 2000–2017. doi: 10.1016/j.apm.2009.10.014. [CrossRef] [Google Scholar]
 Mao J.Y., Yuan S.Q., Pei J., Zhang J.F., Wang W.J. (2014) Applications of different turbulence models in simulations of a large annular volutetype pump with the diffuser, Earth Environ. Sci. 22, 1–9. doi: 10.1088/17551315/22/2/02201. [Google Scholar]
 Pereira I.B., Ribeiro D.C., Romero O.J. (2017) Threedimensional modelling of heat transfer in wellbore during steam injection process, IEEE Latin Am. Trans. 15, 690–697. doi: 10.1109/TLA.2017.7896396. [CrossRef] [Google Scholar]
 Skoczylas P., Alhanati F.J.S. (1998) Flow regime effects on downhole motor cooling, in: SPE Gulf Coast Section Electrical Submersible Pumps Workshop, April 29–May 1, Houston, Texas. [Google Scholar]
 Takács G. (2009) Electrical submersible pumps manual, Gulf Professional Publishing, USA. [Google Scholar]
 Tarcha B.A., Borges O.C., Furtado R.G. (2015) ESP installed in a subsea skid at Jubarte field, in: SPE Artificial Lift Conference, 27–28 May, Bahia, Brazil. doi: 10.2118/173931MS. [Google Scholar]
 Tarcha B.A., Furtado R.G., Borges O.C., Vergara L., Watson A.I., Harris G.T. (2016) Subsea ESP skid production system for Jubarte field, in: Offshore Technology Conference, 2–5 May, Houston, TX, USA. doi: 10.4043/27138MS. [Google Scholar]
 Teixeira V.F., Gessner T., Shigueoka I.T. (2012) Transient modeling of a subsea pumping module using an ESP, SPE Latin American and Caribbean Petroleum Engineering Conference, 16–18 April, Mexico. doi: 10.2118/153140MS. [Google Scholar]
All Tables
Difference of the Nusselt number (Nu) on the fully developed region using case 1 based on the Gnielinski correlation.
All Figures
Fig. 1 ESP in the skid. (a) Schematic representation (not in scale) of the offshore production system with ESP location and problem domain of study. (b) ESP in the skid detail. The pump module houses two ESPs placed in an “X” configuration (Tarcha et al., 2015). 

In the text 
Fig. 2 Schematic representation of the horizontal geometry with details of annular region and eccentricity, horizontal coordinate x is varying from 0 to L. (a) Full domain, (b) concentric tubes, (c) eccentric tubes. 

In the text 
Fig. 3 Boundary conditions of the problem domain, with x varying for 0 to L. 

In the text 
Fig. 4 Nusselt number for each mesh along the annular inner wall. The continuous horizontal line is a result of Gnielinski’s correlation. 

In the text 
Fig. 5 Nusselt number for each mesh in the fully developed region, position x/D_{h} = 20 near the exit plane of the domain, and analytical response. 

In the text 
Fig. 6 Cross section of annular space discretized using mesh 2 with 4 104 000 elements which is selected for all simulations. 

In the text 
Fig. 7 Nusselt number variation along the inner wall of the annular region considering case 1. 

In the text 
Fig. 8 Nusselt number on the annular inner wall of the fully developed region, using case 1 and comparing with Gnielisnki’s correlation. 

In the text 
Fig. 9 Nusselt number on the annular inner wall, which represents the motor surface, in the fully developed region varying with the Reynolds number for three fluids. 

In the text 
Fig. 10 Numbers 1–4 represent the positions of the lines along the inner pipe surface where Nusselt number is plotted in the following figures. 

In the text 
Fig. 11 Nusselt number with Pr = 37 and Re = 51 916 on the four lines (Fig. 10) of geometry with vertical eccentricity of 0.25. The concentric response is also included. 

In the text 
Fig. 12 Nusselt number along lines 1–3 (represented in Fig. 10), with Pr = 37 and four turbulent Reynolds numbers in a geometry with vertical eccentricity of 0.25. (a) Nusselt number along line 1 (θ = 0°). (b) Nusselt number along line 2 (θ = 90°). (c) Nusselt number along line 3 (θ = 180°). 

In the text 
Fig. 13 Nusselt number along the fluid flow direction on the increased length geometry (from L_{1}/D_{h} = 23.36 to L_{2}/D_{h} = 46.72) with bottom vertical eccentricity 0.25, Pr = 37 and Re = 51 916. 

In the text 
Fig. 14 Axial component of the velocity (m/s) in the bottom and top of annulus in the geometry with vertical eccentricity of 0.25 for the fluid with Pr = 37 and Reynolds numbers equal to: (a) 51 916; (b) 42 027; (c) 32 138; and (d) 27 193. Fluid flow is from left to right. 

In the text 
Fig. 15 Vertical component of the velocity (m/s) in the geometry with vertical eccentricity of 0.25, for the fluid with Pr = 37 and Reynolds numbers equal to: (a) 51 916; (b) 42 027; (c) 32 138; and (d) 27 193. Fluid flow is from left to right. 

In the text 
Fig. 16 Nusselt number in positions 1–3 (of Fig. 10) along the domain with eccentricity 0.25 (continuous blue line) and eccentricity 0.5 (dashed black line). Pr = 13 and Re = 77 693. (a) Nusselt number along line 1. (b) Nusselt number along line 2. (c) Nusselt number along line 3. 

In the text 
Fig. 17 Nusselt number with Pr = 13 for four turbulent Reynolds numbers at line 1 (θ = 0°) of geometry with vertical eccentricity of 0.5. 

In the text 
Fig. 18 Nusselt number with Pr = 13 for four turbulent Reynolds numbers at line 2 (θ = 90°) of geometry with vertical eccentricity of 0.5. 

In the text 
Fig. 19 Nusselt number with Pr = 13 for four turbulent Reynolds numbers at line 3 (θ = 180°) of geometry with vertical eccentricity of 0.5. 

In the text 
Fig. 20 Axial velocity (m/s) in geometry with eccentricity 0.5 for the fluid with Pr = 13 and Reynolds numbers equal to: (a) 77 693; (b) 62 894; (c) 48 095; and (d) 40 696. 

In the text 
Fig. 21 Vertical velocity (m/s) in geometry with eccentricity 0.5 for the fluid with Pr = 13 and Reynolds numbers equal to: (a) 77 693; (b) 62 894; (c) 48 095; and (d) 40 696. 

In the text 