Thermal Operational Map of an ESP-In-Skid Motor: A CFD Approach

**Mathematical Modelling**

Two-phase flow model: The governing equations are based on the momentum, mass and energy conservation principles. Since the system is made up of two immiscible fluids, additional hypotheses are necessary to the equations. In the present work, the homogeneous flow hypothesis was considered, i.e., the fields (velocity, pressure and temperature) considered in the model are shared between the phases. The properties (density, viscosity, thermal conductivity, specific heat) are average weighted by the local volumetric fraction. In addition, the RANS (Reynolds Averaged Navier-Stokes) approach was adopted, requiring a turbulence closure model.

Equations (1) e (2) represent global mass and momentum conservation respectively in the Navier-Stokes formulation [14].

$$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \vec{v}) = 0$$

where $\rho$ is the two-phase mixture density and $\vec{v}$ is the two-phase mixture velocity vector.

$$\frac{\partial}{\partial t} (\rho \vec{v}) + \nabla \cdot (\rho \vec{v}) = -\nabla p + \nabla \cdot \mu \left[ \nabla \vec{v} + (\nabla \vec{v})^T \right] + \rho \vec{g}$$

The left-hand side terms in both equations represent the total acceleration. The right-hand side terms of equation (2) represent the pressure forces, viscous forces and gravitational forces, respectively. $\tau$ is the stress tensor. The property $\mu$ is the effective two-phase mixture viscosity. The gravity acceleration vector is $-0.8549978, -9.77267, 0$ m/s$^2$, which is equivalent to an upward 5° axis inclination to the horizontal. Due to the presence of the second fluid, it is necessary to use another equation to describe the interface position between the two fluids. The equation used here is the mass conservation for the liquid phase. This approach is known as Volume-of-Fluid (VOF) and was introduced by Hirt e Nichols.

$$\frac{1}{\rho_i} \left[ \frac{\partial}{\partial t} (\alpha_i \rho_i) + \nabla \cdot (\alpha_i \rho_i \vec{v}) = 0 \right]$$

where $\alpha_i$ is the liquid volume fraction. Therefore, VOF method introduces two additional unknown variables to the system, $\alpha_i e \alpha_g$ (gas volume fraction), subjected to the constraint:

$$\alpha_i + \alpha_g = 1$$

The employed energy conservation equation (5) was simplified using the hypothesis of thermal equilibrium between the phases. This is equivalent to impose no thermal resistance in the interface. In another words, the temperature field is shared among the fluids.

$$\frac{\partial \rho E}{\partial t} + \nabla \cdot (\vec{v} \cdot (\rho E + p)) = \nabla \cdot (k_{ef} \nabla T)$$ (5)

$$E = h - \frac{p}{\rho}$$ (6)

$$h = \int_{\Gamma_{ef}}^{T} c_{p} dT$$ (7)

Equation (6) is the definition of total energy (neglecting kinetic energy term) and Equation (7) is the enthalpy definition which connects the energy conservation and temperature. $k_{ef}$ is the effective two-phase mixture thermal conductivity; $c_{p}$ is the two-phase mixture specific heat. These properties are evaluated using mixing rules weighted by mass fraction, see item 2.3.1.

Thus, the system presents a zero degree of freedom because the number of unknown variables $(\vec{v}, p, \alpha_{i}, \alpha_{g}, e, T)$ is equal to the number of equations (1) to (5) [15].

Motor heat conduction model: To model the motor heat conduction, the hypothesis of homogeneous and isotropic material was employed. The properties (density, thermal conductivity and specific heat) as well as the component part volume were presented by Betônico [12] for each motor section, see 2.3.2. The motor energy conservation equation is:

$$\frac{\partial T}{\partial t} = \frac{k_{m}}{c_{p,m} \rho_{m}} \nabla^{2} T + \frac{Q}{c_{p,m} \rho_{m}}$$ (8)

where $Q$ is the volumetric heat source from the motor and the index $m$ are the motor mixture properties.

**Turbulence Closure Model**

Turbulence model: The $k - \omega$ SST (Shear Stress Transport) model is employed in this work. This model is based on the eddy viscosity concept using two transport equations. The two equations describe: turbulent kinetic energy $(k)$ transport, Equation (9), and turbulent kinetic energy dissipation frequency $(\omega)$. Equation (10):

$$\frac{\partial (\rho k)}{\partial t} + \nabla \cdot (\rho \vec{v} k) = \nabla \cdot [\Gamma_{k} \nabla_{k}] + G_{k} - Y_{k}$$ (9)

$$\frac{\partial}{\partial t} (\rho \omega) + \nabla \cdot (\rho \vec{v} \omega) = \nabla \cdot [\Gamma_{\omega} \nabla_{\omega}] + G_{\omega} - Y_{\omega} + D_{\omega}$$ (10)

where $G_{k} e G_{\omega}$ are the generation terms for the turbulent kinetic energy and the turbulence kinetic energy dissipation frequency, $Y_{k} e Y_{\omega}$ are the dissipation terms of $k e \omega$, $e D_{\omega}$ is the cross-diffusion term. All these terms are defined in Menter [16]. All properties are taken as mixture described in item 4.3.1. The diffusivities $\Gamma_{k} e \Gamma_{\omega}$ are defined in the same way, with respective Prandtl numbers $\sigma_{k} e \sigma_{\omega}$ using the following relation:

$$\Gamma_{i} = \mu_{m} + \frac{\mu}{\sigma_{i}}$$ (11)

where $\mu_{i}$ is the turbulent viscosity and $\mu_{m}$ is the two-phase mixture viscosity as described in item 4.3.1. The eddy viscosity is based on the Bradshaw hypothesis [17] which states that for a boundary layer flow the Reynolds stresses are proportional to the turbulent kinetic energy, as shown in equation (12).

$$\mu_{i} = \frac{\rho \kappa}{\omega} \frac{1}{\max \left(1 / a^{*} ; SF_{2} / a_{1} \omega\right)}$$ (12)

Where $a^{*} = a_{2} e F_{2}$ are defined in Fluent [18]. $S$ is the mean flow strain rate.

The effective viscosity is defined by:

$$\mu = \mu_{m} + \mu_{i}$$ (13)

The effective viscosity is used in Equation (2).

Turbulent heat convection: When the convective heat transfer is in turbulent regime, the effective thermal conductivity is given by:

$$k_{ef} = k + \frac{c_{p} \mu_{i}}{pr_{i}}$$ (14)

where $c_{p}$ is the two-phase mixture specific heat, $k$ is the two-phase mixture thermal conductivity and $Pr_{i}$ is the turbulent Prandtl number with constant value of 0.85.

**Physical Properties**

Fluids properties: The mixing rules are dependent of phases volume fractions. For mixture density $\rho$ and mixture viscosity $\mu_{m}$ are defined as:

$$\rho = \alpha_{g} \mu_{g} + \alpha_{1} \rho_{1}$$ (15)

$$\mu_{m} = \alpha_{g} \mu_{g} + \alpha_{1} \mu_{1}$$ (16)

For each property, a specific correlation was used according to oil characteristics. The Gas-Oil Ratio (GOR) of the API 20° oil described by the Tarcha, et al. [2] work, was 50 sm³/sm³. For most correlation used, the specific gravity (SG) is needed and this information was neglected in the Tarcha, et al. [2] work, which an arbitrary value of 0.65 was applied.

For the thermodynamic properties (bubble point, solubility ratio, formation volume factor and viscosity), DeGhetto, et al. [19] correlation was employed since it is recommended for oils with degree API less than 22. The calculation of the gas thermal conductivity was based on the hypothesis that the gas was constituted of pure methane. The NIST (U.S. National Institute of Standards and Technology) [20] database was used for this property. The compilation of all properties correlations is presented in Table 1.

The fluids properties depend on motor operation point. For this analysis were chosen four operational rotation conditions which correspond to four different frequencies of the electrical motor: 40 Hz, 45 Hz, 55 Hz and 60 Hz. Table 2 shows the values for each property for each motor frequency operational condition.

The missing properties in the 60 Hz column are due to the fact that at this condition the oil is undersaturated or all gas is in solution.

Motor physical properties: The motor physical properties were taken from Betonico [12] work. The author describes all internal parts of the motor and their properties. In order to assume a homogeneous and isotropic material, a volume weighted mean of each motor component was adopted. The volume weighted mean properties are listed in the Table 3 below.

Geometry, Mesh and Boundary Conditions

The system geometry and boundary conditions were collected from literature (Tarcha, et al.) [2] and field data.

The motor geometry, motor OD, motor length, motor position inside capsule, seal side and position and other relevant parts of the ESP-in-skid were provided by Petrobras SA and Schlumberger’s ESP catalogue for pumping unit series 738 with 750 hp [24]. The final schematic of the ESP is presented in the Figure 2.

Click to enlarge

The red rectangle depicts the external limits (boundaries) of the fluid domain (volume where the mathematical model was solved). The blue arrow (ESP intake) indicates the gas-liquid mixture inlet boundary and is positioned on the left side of the fluid domain. The red arrow (Pump intake) denotes the fluids outlet boundary located at the fluid domain right side. In the real equipment, this boundary is just before the pump intake. The upper and bottom boundaries are solid walls which models the capsule external casing which is in contact with seawater. The internal boundaries are the motor and seal surfaces.

The blue rectangle (Motor) represents the solid domain. In this domain will be applied the heat conduction model with heat generation. The orange rectangle (Seal) is an adiabatic region and represents the motor seal (motor-pump separator).

The capsule is tilted 5° w.r.t. horizontal plane.

The computational mesh was built with solely hexahedra elements. To simplify the model, the motor and the seal were considered perfect cylinders. This simplification allowed to use a reduced number of elements (390,160) for applying to this geometry 12 radial equally spaced elements in the annular region as shown in Figure 3. The minimal orthogonal quality is 0.7022 and the maximum aspect ratio is 9.67.

Click to enlarge

The operational points for each case considered for this simulation is summarized in Table 4. The operational point 60 Hz is a single-phase flow because all gas is dissolved in the oil at the ESP intake due to a fluid pressure considered to be above the bubble point, 98 bar, from DeGhetto, et al. [19] correlation. As the motor frequency decreases, the gas volume fraction increases. The 40-60 Hz range of frequencies allowed an extended gas volume fraction effect analysis on the heat transfer.

The momentum, energy and mass boundary conditions are set based on field data and the heat generated by the motor estimated from Betonico, et al. [25]. A homogeneous mixed oil and gas fluid is considered at the ESP intake. The electrical power converted into heat was not reported for the ESP model in the Betonico, et al. [25] work. For that, the Burt, et al. [26] work was used. The authors suggest that 2% of total electrical power input is converted into heat power. This energy amount was considered to be generated equally in each volume cell inside the motor. The obtained data for each boundary condition and heat source are shown in Table 5.

The oil mass flow is also expressed in oil flow rate (m³/h) because this is the variable (and units) used in field operation. Furthermore, the results of temperature are also presented in function of oil flow rate.

Numerical Approach

A High-Performance Computer (HPC) with 5 nodes (2 quad-core Intel Xeon each) and 32 Gb RAM were used in transient condition simulations.

The SIMPLE (Semi-Implicit Method for Pressure-Linked) pressure-velocity coupling was used. PRESTO! and Second Order Upwind schemes were employed for the pressure and velocity interpolation respectively. To evaluate the variables gradients the Least Squares Cell Based method was applied to each cell. The Second Order Upwind was chosen to obtain the turbulent kinetic energy and the dissipation rate. The Geo-Reconstruct scheme was used to evaluate the volume fraction equation.

For all simulations, the CFL condition was attained. The Courant number is defined by the Equation (17).

For all cases, two modeling challenges were faced:

• Velocity profiles and flow pattern at the ESP intake domain were unknown;
• Only one pressure point measure was available from the field data, i.e., no pressure drop was available to be compared to the simulation results.

An important numerical strategy issue for this case is the fact that the fluid flow and heat conduction problems have distinct time scales. To overcome that, the following procedure to achieve a global solution was adopted:

• All fluid flow equations with timestep limited by the fluid flow problem (0,01 s for Co = 1) was solved until the fluid mean temperature variation at the pump intake stabilized within 1% (~130 s);
• Then, the fluid flow equations calculations were halted and a larger timestep (10 s, limited by the heat conduction problem) was used until the maximum motor temperature reaches a constant value within 1% margin (~13 h);
• Finally, the fluid flow equations calculations were re- started with the original timestep (Co = 1) until the fluid mean temperature at the at the pump intake stabilized within 1% margin.
• In all cases, the initial condition started from a stabilized fluid condition at 4 °C in all domains (fluid and motor).

Monitoring points were selected to control the simulations and to judge if the solution achieved or not a steady state. The selected monitoring points used were:

• Fluid mean temperature at the pump intake;
• Mean temperature at the motor surface;
• Maximum temperature inside the motor;
• Mean gas volume fraction at the pump intake;
• Gas volume inside fluid domain.

The temperatures calculated from the simulations are then normalized using the Equation (18).

$$ \theta_ {f} = \frac {T}{T _ {i n l e t}} (1 8) $$

f inlet The motor heating dynamics is evaluated using the resulted nondimensional temperature by the Equation (19).

$$ \theta_ {m} = \frac {T - T _ {s e a}}{T _ {\max , m} - T _ {s e a}} \tag {19} $$ sea m m sea max, Where ( ) , @ max m T max T motor domain = , 4 sea T C = ° .

The mean temperature at motor surface is defined by:

$$ T _ {\text {m e a n}, m s} = \frac {1}{A _ {m s}} \sum_ {i = C e l l s} T _ {i}. d A _ {i} \tag {20} $$ Where ms A is the motor surface area, the index i represents the superficial cells of motor external surface, iT is the surface temperature and i dA is the superficial cell area.

Motor and Fluids Heating Dynamics

The results presented in this section are the outcome of the first step of the simulation strategy for the 55 Hz case. Similar results are observed for all other cases. The fluid temperature is initially at 4 °C As can be seen in Figure 4, the fluid temperature at the pump intake reaches a constant value ( ) 1 f θ ≈ in 60s (1 minute), whose time delay is due only to the fluid residence time inside the fluid domain volume. During this time, a little increase in the heat transfer occurs between the hot motor and the cold fluid.

Click to enlarge

After the fluid reaches a constant temperature, the fluid flow equations solving is halted to accelerate the heat conduction in the motor. The motor heating dynamics (Figure

5) is much longer, taking, for the present case, approximately 6 h to reach 95% of the equilibrium temperature. For this reason, the time unit used in this plot is hours.

Click to enlarge

At this point, both the mean motor temperature and the pump intake fluid temperature are stabled but not at the same time. Thus, the next and last step is solving simultaneously both the fluid flow and heat conduction models.

Surface Motor Temperature Field

The motor surface temperature fields are shown in the Figure 6 where the fluid is flowing from left to right and the gravity is pointing down. The field values were evaluated using the Equation (19).

Click to enlarge

(a) It is possible to notice that the separation of fluids causes a temperature field asymmetry in the two-phase flow cases (Figures 6 (a), (b) and (c)). When only one phase flows, the temperature field is symmetrical (Figures 6 (d)). There is also a temperature increase at the lower region of the motor top near the seal (far right of each graph in Figures 6). It also can be inferred that the turbulent flow instabilities increase the convection at the upper part of the motor bottom (far left of each plot in (Figures 6). This statement is corroborated by observing both the temperature fields and the velocity field plotted together in the Figure 7.

Click to enlarge

(b)

Click to enlarge

(c)

Click to enlarge

(d) Figure 6: Dimensionless temperature fields on the motor surface at (a) 40 Hz; (b) 45 Hz; (c) 55 Hz; (d) 60 Hz.

Click to enlarge

(a)

Click to enlarge

(c)

Click to enlarge

(b)

Click to enlarge

(d) Figure 7: Dimensionless temperature fields on the motor surface and velocity vectors colored by gas volume fraction in a vertical plane: zoom at the motor bottom (inlet). (a) 40 Hz; (b) 45 Hz; (c) 55 Hz; (d) 60 Hz.

Figure 7 shows the greater the amount of gas, the greater the superficial velocities of the liquid and, consequently, the greater the convection in the upper region of the motor bottom. The single-phase case (d) shows a symmetric field as expected. Figure 8 shows the same fields as shown in the Figure 7, but at the motor top close to the seal.

As observed, the two-phase flow cases (Figure 8 (a), (b) and (c)) present an asymmetric pattern in the vertical direction. As the gas drags the liquid in the upper region, the velocities are greater which makes the convection there to be greater than in the lower region. In the lower speed case (Figure 8 (a) and (b)), a greater difference between the maximum and minimum temperature may be expected.

Click to enlarge

(a)

Click to enlarge

©

Click to enlarge

(b)

Click to enlarge

(d) Figure 8: Dimensionless temperature fields on the motor surface and velocity vectors colored associated to the gas volume fraction in vertical plane: motor top (seal). (a) 40 Hz; (b) 45 Hz; (c) 55 Hz; (d) 60 Hz.

Velocity Field Effect on the Motor Surface Temperature

If a velocity profile along a vertical line such as that is observed in the Figure 9 for the 55 Hz case, a quantitative analysis can be made. The velocity behavior is represented by the four lines positioned along the normalized length of the motor (dimensionless ratios of 0.06, 0.30, 0.76 and 0.88). The velocity represented in the Figure 9 plots are normalized by / V V ρ .

Click to enlarge

The speed difference between the upper and lower regions is noticeable. As convection is not only related to the speed, if there is solely gas in the upper region, the convection may not be efficient as expected. Figure 10 shows the volume fraction behavior of the gas volume fraction along the same lines of the Figure 9 for this case.

Click to enlarge

As may be observed, the gas accumulates near the outer motor wall. Therefore, the fluid contact with the motor surface is mainly composed by its liquid phase. For the purpose of a physical analysis, the gas behaves as a shroud liquid accelerator by limiting its available flow area.

Another analysis can be made involving the correlation between the temperature field inside the motor and the fluid around it for the 55 Hz case. This analysis is shown in Figure 11.

Click to enlarge

In Figure 11 temperature is normalized by the Equation (19). The fluid temperature asymmetry behavior is less clearly observed because the dimensioning scale is much higher (maximum temperature inside the motor). The fluid temperature behavior is characterized by a close sinusoidal shape, related to the characteristic of the double thermal boundary layer formed at the inner wall (motor surface) and outer wall (capsule surface). The maximum temperature calculated for the motor solid domain is close to the seal temperature (at 0.76) due to the seal adiabatic boundary condition.

Finally, Figure 12 shows the results of three characteristic system temperature analysis, including the maximum motor domain temperature.

Click to enlarge

The higher the flow rate (or, in the present case, the frequency) the greater the heat generated and the maximum motor temperature is greater. On the other hand, the higher the flow rate, the greater the convection, since the difference between the maximum motor domain temperature (Equation (19)) and the mean motor surface temperature increases (Equation (20)). The increase in fluid temperature is not noticeable within the scale used, but it becomes more distinguished as the flow rate increases.

Simulated and Field Measured Thermal Operational Map Comparison

An estimate of the motor temperature as a function of the flow rate or the operational frequency (40 Hz, 45 Hz, 55 Hz, 60 Hz) provided by VSD (Variable Speed Drive) can be made. Figure 13 shows the numerical results, which are coherent with the field data. The results are satisfactory in spite of the uncertainties related to a mean deviation of 5% and to a maximum deviation of 10% in the 60 Hz operational case.

The difference increase between the simulated and field data as the flow rate increases, may be due to the less accurate PVT correlation properties of the related pressure and temperature operational conditions. Also, the arbitrated value for the specific gravity of the gas may be higher than the real ones. The single-phase flow assumption at the highest frequency may provide another explanation for the greater difference, since there may have additional gas in the real operational condition.

Click to enlarge