Heat dissipation of the Electrical Submersible Pump (ESP) installed in a subsea skid

. 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 ef ﬁ cient to ensure the continuity of its operation. The heat withdrawal is performed by the ﬂ uid produced. The purpose of this article is to understand the process of electric motor cooling to the single-phase and turbulent ﬂ ow 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 single-phase ﬂ ow. The standard j - e model with improved wall function (Enhanced Wall Treatment) is used to closure turbulence problem. This study considered ﬂ ow rates range of 2200 – 4200 m³/d (representing Reynolds numbers range of 27 000 – 133 000 approximately), Prandtl numbers 7 – 37, three con ﬁ gurations of different annular geometries, one concentric and two eccentric, together with the condition of the constant temperature on the motor surface (130 (cid:1) C) and capsule (4 (cid:1) 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 ﬂ ux 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 ﬂ ow 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.


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.
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 j À e 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: Re À 1000 ð Þ Pr 1 þ 12:7 f 8 À Á 0:5 ðPr 2=3 À 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 single-phase 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 j À e 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.

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, 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). 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.
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 e 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): 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 j À e turbulence model with enhanced wall function (j À e 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 two-phase fluid, which is a mixture of fluids A and B, obtained by homogeneous modeling as described in Section 3.1.
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 qc p is the heat capacity, and represents how much energy a material stores per unit volume (Cengel, 2002).

Mathematical modeling
A relevant concept to deal with problems in annular regions is the hydraulic diameter, D h , obtained by: 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: where q is the fluid specific mass, V is the mean fluid velocity and l is the fluid viscosity.
(ii) The Prandtl number (Pr), which describe the relative thickness of the velocity (d) and the thermal (d T ) boundary layers, is defined as the ratio of momentum diffusivity (t) to thermal diffusivity (a), Pr ¼ t a ¼ ðl=qÞ k f =ðqcpÞ , or, 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 d T ) d, d T develop faster than d), 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 d is consistently much greater than d T (d T ( d). (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: 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.

Homogeneous modeling of two-phase 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): 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 (q m ), viscosity (l m ), specific heat (c pm ), and thermal conductivity (k m ), which are, respectively, represented by equations (8)-(11) being q o and q w the densities of the oil and water, l o and l 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 single-phase fluid. All these parameters are pre-calculated and entered as input values.

Isothermal turbulent flow
The isothermal turbulent flow is modeled using the continuity equation (12), the momentum conservation equation (13) and two-equation turbulence model j À e with enhanced wall function, equations (14) and (15): whereũ is the velocity vector.
whereg j j = 9.8 m/s 2 ,s is the stress tensor and p is the pressure (Pereira et al., 2017).
The j À e model is based on the concept of eddy viscosity, l t . Thus, this model solves two different transport equations: turbulent kinetic energy, j, equation (14), and turbulent kinetic energy dissipation rate, e, equation (15) (ANSYS, 2017): where P j is the turbulence production due to viscous forces; C s1 , C s2, C l , r j , and r e are constants with values: 1.44, 1.92, 0.09, 1.0, e 1.3. The turbulent viscosity l t is a property of the turbulent flow and is defined in the equation (16): where c l 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 so-called enhanced wall functions. This requires the mesh to be sufficiently refined to solve the viscous sublayer (ANSYS, 2017).

Non-isothermal 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: H is the specific enthalpy, which is a function of pressure and temperature T. The enthalpy is described thermodynamically by the relation below: where T ref is the reference temperature.

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): b ¼ ð3:85 Pr 1 3 À 1:3Þ 2 þ 2:12 ln Pr ð Þ; ð22Þ 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: 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).

Boundary conditions
The boundary conditions that delimit the domain of the problem are represented by letters in Figure 3 and we have: (a) Inlet -fluid with temperature of 62.5°C and prescribed uniform velocity which depends of each analyzed condition, Table 2.
(b) Outletprescribed pressure equals to 0 Pa. (c) Outer wallconstant temperature of 4°C with nonslip condition and impenetrability. (d) Inner walltemperature of 130°C with non-slip condition and impenetrability.
Four operating flow rates of the ESP in the skid system were considered, being 4200, 3400, 2600, and 2200 m³/d (Tarcha et al., 2016), as defined in Table 2.

Discretization algorithms
The discretization of the governing equations used the FVM with co-located 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 Semi-Implicit Method for Pressure-Linked (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.

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: where u s 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.
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) 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 cross-sectional 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.
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 CFD-Post v15.0.

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.
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.

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.
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).

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  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~1 4, 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~1 855) is not showed in the vertical axes.
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 (h = 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~1 855. 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.
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.
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.

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 (e).
In Figure 16a, it is noted that the Nusselt number at line 1 of the geometry with e = 0.25 shows a positive variation at a shorter entry length, x/D h > 12, while the Nusselt number of the geometry with e = 0.50 starts to vary when x/D h > 16. However, at the end of the geometry the Nusselt number of the geometry with e = 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 e = 0.25 presents a small negative variation, almost imperceptible, however the Nusselt number of the geometry with e = 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 e = 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 e = 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 e = 0.50 occur due to the decrease of the open area to the flow at the bottom when compared to the geometry with e = 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 (h = 0°), 2 (h = 90°), and 3 (h = 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.
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 (d 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.
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 non-uniform and low efficient heat transfer from electric motor to surroundings.

Conclusion
Based on the presented results, the following can be concluded:     1. 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. 2. 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. 3. 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. 4. For the single-phase 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 (e) and angle (h), and can be written as follow: further studies are necessary to develop correlations for this geometry.
5. 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. 6. 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.