Heat Transfer in a Drilling Fluid with Geothermal Applications

The effects of various conditions on the fluid flow, particle migration and heat transfer in non-linear fluids encountered in drilling and geothermal applications are studied. We assume that the drilling fluid is a suspension composed of various substances, behaving as a non-linear complex fluid, where the effects of particle volume fraction, shear rate, and temperature on the viscosity and thermal diffusivity are considered. The motion of the particles is described using a concentration flux equation. Two problems are studied: flow in a vertical pipe and flow between two (eccentric) cylinders where the inner cylinder is rotating. We consider effects of earth temperature, the rotational speed of the inner cylinder, and the bulk volume fraction on the flow and heat transfer.


Introduction
Drilling is an essential aspect of many geothermal energy applications (see Younger [1]).It allows the capture of high-density thermal energy stored in the deep Earth [2].The fluid, which flows in the well during the drilling operation is called the drilling fluid/mud, which is a solid-fluid suspension composed of particles in a fluid.Geothermal energy is the thermal energy generated and stored in the Earth which can be extracted and used [3].Similar to other natural energy resources, such as sunlight, wind, rain, tides, and waves, geothermal energy is also considered as a renewable or semi-renewable source of energy [4,5].This is one of the advantages of geothermal energy compared with fossil fuels; another advantage is that it is environmentally more friendly [6][7][8].The development and utilization of geothermal energy has grown in the last few decades.According to Finger and Blankenship [3], in 2005 the worldwide electricity production from geothermal energy was responsible for reducing the consumption of more than thirty million barrels of oil.As indicated by Lund [9], heat pumps using geothermal sources have 'the largest energy use' which can be utilized.
High concentration of geothermal energy is usually available in the deepest parts of the Earth.For example, in most of North America, in the regions to the east of the Rocky Mountains, the ground temperature is about 135 • C at depths of 5-7 km [3].The main objective is to access the stored geothermal energy.Fluid flowing in the well during the drilling operation is called the drilling fluid or the drilling mud; it is generally considered to be a nonlinear complex fluid comprised of cuttings (solid particles) with various sizes and shapes suspended in a base fluid and is primarily used to Energies 2017, 10, 1349 2 of 18 transport the particles/cuttings generated during the drilling process from the bottom hole to the surface [10,11].Successful operation of drilling requires, among other things, a thorough understanding of the rheological and thermal properties of the drilling fluid, [11][12][13]; this has motivated numerous experimental and theoretical studies [12,[14][15][16].For example, heat transport in the drilling fluid and temperature control at the drill bit plays a major role in the success of enhanced geothermal system (EGS) technology, which is applied in deep subsurface areas where the temperature is in the range of 150-200 • C (see Rybach [17]).
In most cases, the drilling fluid is modeled using the Power law model, the Bingham model or the Herschel-Bulkley model (which combines the effects of the Bingham and the power-law models [18][19][20]).For the Herschel-Bulkley model, it is assumed that the material behaves like a solid when the (shear) stress is below the yield stress, and like a shear-thinning or shear thickening fluid when the (shear) stress exceeds the yield stress [21].Moreover, the rheological properties of the drilling fluid can also depend on pressure, volume fraction of the particles, shear rate, temperature, etc. Krieger [22] performed experimental and theoretical studies to describe the relation between volume fraction, shear rate and shear viscosity of a suspension .Briscoe et al. [23] proposed a viscosity model for concentrated suspensions, which shows the dependency of viscosity on pressure and temperature.Their model is suitable for high pressure and high temperature conditions, such as deep ocean drilling.In addition, the thermal diffusivity of the drilling mud also shows nonlinear behavior and depends on the particle concentration, size and shape of the particles [16,[24][25][26][27].As indicated by Agemar et al. [28], the most important parameter when considering the thermal aspects of a geothermal reservoir is temperature and its variation.
In most cases, the drilling fluid is modeled using the Power law model, the Bingham model or the Herschel-Bulkley model (which combines the effects of the Bingham and the power-law models [18][19][20]).For the Herschel-Bulkley model, it is assumed that the material behaves like a solid when the (shear) stress is below the yield stress, and like a shear-thinning or shear thickening fluid when the (shear) stress exceeds the yield stress [21].Moreover, the rheological properties of the drilling fluid can also depend on pressure, volume fraction of the particles, shear rate, temperature, etc. Krieger [22] performed experimental and theoretical studies to describe the relation between volume fraction, shear rate and shear viscosity of a suspension .Briscoe et al. [23] proposed a viscosity model for concentrated suspensions, which shows the dependency of viscosity on pressure and temperature.Their model is suitable for high pressure and high temperature conditions, such as deep ocean drilling.In addition, the thermal diffusivity of the drilling mud also shows nonlinear behavior and depends on the particle concentration, size and shape of the particles [16,[24][25][26][27].As indicated by Agemar et al. [28], the most important parameter when considering the thermal aspects of a geothermal reservoir is temperature and its variation.
Here, we model the drilling fluid as a suspension with nonlinear rheological properties, where the viscosity and the thermal conductivity depend on the volume fraction, shear rate, temperature etc., and the motion of the solid particles depends on the particles collision, non-uniform distribution of the viscosity [29] and turbulent diffusivity.It is worth mentioning that, in general, suspensions can also be modeled using multi-component approaches (see [30][31][32]).In Section 2, the governing conservation equations for mass, linear and angular momentum, solid concentration and energy are discussed.In Section 3, the constitutive equations for the stress tensor of the suspension, the non-linear diffusive particle flux, internal energy and the heat flux vector are provided.Then in Section 4, the governing equations are non-dimensionalized and two flow conditions are considered and numerically solved: flow in a pipe and flow between two eccentric cylinders where the inner cylinder is rotating.We study the effects of various conditions encountered in geothermal well drilling operations, such as the earth temperature, the rotational speed of the inner cylinder, and the bulk volume fraction of solid particles on the flow and heat transfer.

Governing Equations
In general, using the methods of continuum mechanics, a suspension can be mathematically studied using three different methods: (1) the single phase approach where the suspension is modeled as a single-component fluid; (2) the non-homogenous single phase model where the particle motion is studied using a particle concentration equation with proper nonlinear flux terms [29,33]; (3) the two-fluid (two-phase) approach, e.g., Mixture theory approach [31,34,35], where the fluid and solid components are coupled through interaction forces.In current work, the non-homogenous single phase model will be used; this approach retains several important features of the two-fluid (two-phase) approach, with lower computational cost [36].If the effects of electro-magnetism and chemical reactions are ignored, the governing equations, for a single component suspension, are the conservation equations for mass, linear and angular momentum, particle concentration, and energy [11,33,36,37].

Conservation of Mass
The conservation equation of mass is, where is the density of the suspension, φ is the volume fraction of the particles, ρ f 0 and ρ s0 are the pure density of fluid and solid components in the reference configuration, i.e., before mixing the suspension is formed, respectively; v is the velocity vector and ∂/∂t is the partial derivative with respect to time.For an isochoric motion:

Conservation of Linear Momentum
The conservation of linear momentum is, where b is the body force vector, T is the Cauchy stress tensor, and d/dt is the total material time derivative, which is given by d(.) / dt = ∂(.)/∂t + [grad(.)]v.The conservation of angular momentum indicates that in the absence of couple stresses the stress tensor is symmetric, that is, T = T T [36].

Conservation of Volume Fraction (Convection-Diffusion Equation)
The conservation of volume fraction is, where the first term on the left-hand side denotes the rate of accumulation of particles, the second term denotes the convected particle flux, and the term on the right-hand side denotes the diffusive particle flux [11].Following Phillips et al. [29], the particle flux N is assumed to be related to the Brownian motion (see [38,39] for additional details), the variation of interaction frequency and the viscosity.

Conservation of Energy
The conservation of energy reads, Energies 2017, 10, 1349 where q is the heat flux vector, L is the gradient of velocity, and r is the specific radiant energy.
is the heat capacity of the suspension and c p f 0 and c ps0 are the specific heat capacity of the pure fluid and the solid particles in reference configuration.Thermodynamic requires the application of the second law of thermodynamics [40].The local form of the entropy inequality is given by ( [40]): where η(x, t) is the specific entropy density, ϕ(x, t) is the entropy flux, and s is the entropy supply density due to external sources, and the dot denotes the material time derivative [40].If ϕ = 1 θ q, and s = 1 θ r, where θ is the absolute temperature, then Equation ( 5) reduces to the more familiar Clausius-Duhem inequality The effects of the Clausius-Duhem inequality are ignored in this paper; we present the above equation for the sake of completeness [41][42][43][44].In order to 'close' the above set of equations, constitutive relations are needed for T, q, and N. We ignore the effects of radiation.

Stress Tensor
In general, the constitutive equation for the stress tensor of a concentrated suspension can be given by a non-Newtonian (non-linear) model (see [22,45]) and it may include a viscous stress and a yield stress (see [11,12]): where T y and T v are the yield stress and the viscous stress, respectively.Here, we focus on the viscous stress and assume that the suspension can be modeled as a viscous fluid, where µ m is the viscosity of the suspension, I is the identity tensor, D is the symmetric part of the velocity gradient, and p is the pressure [36].In general, the shear viscosity can depend on, where t is the time, .
γ is some measure of the shear rate, θ the temperature, φ the concentration, p the pressure, E the electric field, and B the magnetic field [46].In this paper, we will not consider the effects of time, electric or magnetic fields on the viscosity.We assume the viscous stress part of the suspension can be represented by a generalized power-law fluid model [22,23,33], where: Here µ t is the turbulence eddy viscosity, θ is temperature, tr is the trace operator, Π is an invariant of D, E u is related to the "activation energy", k B is the Boltzmann constant, C 10 , C 10 and θ c are constants related to temperature effects, and m is a measure of the shear dependency of the viscosity.

The expression, µ
is based on Krieger's work who provided a correlation for the viscosity of the suspension as a function of the particle concentration (µ = µ f 0 (1 − φ/0.68) −1.82 ) [22]; this paper included theoretical analysis and experimental studies using various sizes and types of particles, φ max is the volume fraction of the particles under the maximum packing condition.µ r in Equation ( 12) is a generalization of the model proposed by Briscoe (1994) [23].In the remainder of this study, we ignore the shear rate dependency of the viscosity.However, for suspensions with smaller size particles (such as nanofluids), the effect of the shear rate on viscosity can be significant [46,47].The parameters in the expression for µ r depend on the properties of the specific drilling mud which need to be determined using experiments.Here we use the data used in the experiments performed by Wang et al. [16]; the values of the following parameters are fitted to their data: 3226 K and θ c = 182.3742K.In these measurements, the shear rate ranged from 0 to 1020 s −1 ; the shear viscosity was found to be independent of the shear rate, see Wang et al. [16].The bulk volume fraction is 0.345.The drilling mud is composed of water, various additives and solid particles.The density of the drilling mud is 2.2 g/cm 3 .Comparison between the fitted expression and the experimental data are shown in Figure 1.We should mention that there are other models for drilling fluids for high pressure and high temperature (HPHT) conditions.Most of these models are based on the Herschel-Bulkley model (see Khan and May [48], Salehi et al. [49], Vajargah and van Oort [50]).Our model not only includes the effects of temperature and shear-rate dependency of the viscosity, but also allows for a volume fraction dependent viscosity and thermal conductivity.Next we discuss the modeling of the particle flux N.
Energies 2017, 10, 1349 5 of 17 Khan and May [48], Salehi et al. [49], Vajargah and van Oort [50]).Our model not only includes the effects of temperature and shear-rate dependency of the viscosity, but also allows for a volume fraction dependent viscosity and thermal conductivity.Next we discuss the modeling of the particle flux .

Particle Flux
The particle transport flux (the volume fraction equation) can be due to sedimentation, turbulent diffusivity, Brownian motion, etc. [13,29].In this study, we ignore the Brownian motion.Thus, the diffusive particle flux can be written as: where , , and are the contributions due to particles collision, spatial variations in the viscosity, gravity and turbulence, respectively.In this paper, we also assume = .Based on the proposal of Phillips et al. [29], we assume and are given by: = − (ln ) where is the radius of the particle, is the local shear rate, is the effective viscosity, and ,

Particle Flux
The particle transport flux (the volume fraction equation) can be due to sedimentation, turbulent diffusivity, Brownian motion, etc. [13,29].In this study, we ignore the Brownian motion.Thus, the diffusive particle flux can be written as: where N c , N µ , N g and N t are the contributions due to particles collision, spatial variations in the viscosity, gravity and turbulence, respectively.In this paper, we also assume N g = 0. Based on the proposal of Phillips et al. [29], we assume N c and N µ are given by: .
where a is the radius of the particle, .
γ is the local shear rate, µ s is the effective viscosity, and K c , K µ are empirically-determined coefficients.According to Seifu et al. [51] the above two terms in the particle flux are given by, We can see that ln( . ) is a field potential incorporating mechanism which can cause migration of the particles.Following the Phillips et al. and Subia et al. [13,29], in this paper we assumed that K c = 0.43, K µ = 0.62, as constants.
The turbulent diffusivity flux is given by [52], where ν t is the eddy diffusivity and Sc is the Schmidt number assumed to be 0.9 for the present study [53,54].

Heat Flux Vector
As stated by Fourier [55] (see also Winterton [56]), the heat flux vector is proportional to the gradient of temperature, where k is the thermal conductivity of the material [36].For complex materials, k can depend on various parameters/fields, such as concentration, temperature, shear rate, etc.; in these situations k is usually replaced with k m or k e f f , called the effective thermal conductivity.In general, k is a second order tensor for anisotropic materials (for more information see [57], p. 129 and [58][59][60]).When turbulence is considered, there is an additional contribution to the conduction usually designated as k t [61,62], where Pr t is the Prandtl number and can be chosen as 8/9 [61].An effective thermal conductivity model derived by Jeffrey [26] is: where: ξ = 3ξ 2 + 3ξ 3 4 + 9ξ 3  16 where: Energies 2017, 10, 1349 7 of 18 Jeffrey's model includes second order effects in the volume fraction, see also [63].It is worth mentioning that Equation ( 22) is a (corrected) linearization of the Maxwell formula, which can be rewritten as, For more information about the effective heat conductivity of suspensions, see [64].In our paper here, we use the model proposed by Jeffrey [26].The heat conductivity of the suspending fluid (water), κ f , is assumed to be 0.6485 W/(m•K) [65].
In the following section, we present the three-dimensional forms of the governing equations.After non-dimensionalizing these equations, we study two different problems (1) flow in a long vertical pipe; and (2) flow between two eccentric rotating cylinders.For each geometry, we first study problems with isothermal condition, and then we look at the effect of the temperature, namely problems with non-uniform temperature distribution.

Results and Discussion
Substituting Equations ( 13)-( 14), ( 14)-( 19) and ( 20)- (24) in Equations ( 1)-( 5), a set of partial differential equations (PDEs) are obtained.The PDEs are given below: The above equations are used to solve for the unknowns which are, v, p, φ and θ.The dimensionless forms of the above equations are, where: Energies 2017, 10, 1349 8 of 18 where H r is a reference length, for example, the diameter of the pipe, u 0 is a reference velocity and is chosen as the axial mean velocity at the inlet, and θ 1 and θ 0 are reference temperatures and are chosen as 393 K and 293 K, respectively.The asterisks have been dropped for simplicity.Among the above dimensionless numbers J c and J µ are related to the particles flux; D s is related to the heat dissipation.For obtaining the turbulent eddy diffusivity ν t (and the turbulent eddy viscosity µ t ), different turbulence models, such as k − ω or k − ε RANS model, can be selected depending on the flow conditions [66][67][68].In this paper, we also ignore the effect of the turbulence, due to the low Reynolds number of the problems studied.Table 1 lists the value of the physical properties applied in this paper.The suspension viscosity, effective thermal conductivity and K c are calculated by the models/correlations discussed in the previous section.Numerically, the above governing equations are solved by developing PDE solver using the libraries of OpenFOAM ® [69].The PIMPLE algorithm is applied for dealing with the incompressibility condition [69].The details of how the governing equations are discretized in OpenFOAM, the PIMPLE algorithm, the numerical schemes, etc., are given on [69][70][71][72][73].For ensuring numerical stability and accuracy, the value of the time step is chosen so that the maximum Courant Number is always less than 0.1.The Courant Number is the portion of a cell that a material will transport by advection in one time step.The boundary conditions are provided in Table 2, for more detail see [71].For each problem, the simulation domain is discretized as hexahedral meshes using ICEM [74].The mesh-dependence studies are performed for each case.
In the following cases studied, we assume D s = 0, that is we ignore the effect of viscous dissipation; and due to the low Reynolds number of the problems we also ignore the effect of the turbulence on flow, heat transfer and particles transport.[29,65,75,76].

Physical Property Value Physical Property Value
Boundary conditions.For more discussion on the boundary condition, see [71].

Vertical Flow in a Pipe
First, we study the flow in a vertical pipe, as shown in Figure 2. The pipe is assumed to be embedded in earth.To save computational cost while simulating the flow in a long pipe, we use the symmetry condition and assume that the flow is two dimensional.Furthermore, to ensure a steady-state condition, the simulation time is always set to be greater than 5•(L/U), where L is the length of the pipe and U is the mean velocity of the flow.

Vertical Flow in a Pipe
First, we study the flow in a vertical pipe, as shown in Figure 2. The pipe is assumed to be embedded in earth.To save computational cost while simulating the flow in a long pipe, we use the symmetry condition and assume that the flow is two dimensional.Furthermore, to ensure a steadystate condition, the simulation time is always set to be greater than 5•(L/U), where L is the length of the pipe and U is the mean velocity of the flow.

Effect of Bulk Volume Fraction and the Reynolds Number
From Equation (11) we can see that the concentration of the solid particles affects the viscosity of the suspension, especially when the volume fraction is close to the maximum packing volume fraction.We assume isothermal conditions, i.e., ̅ is zero.Based on the flow conditions at the inlet, the values of the dimensionless numbers are: = 159, = 3.87 × 10 and = 971.Here the

Effect of Bulk Volume Fraction and the Reynolds Number
From Equation (11) we can see that the concentration of the solid particles affects the viscosity of the suspension, especially when the volume fraction is close to the maximum packing volume fraction.We assume isothermal conditions, i.e., θ is zero.Based on the flow conditions at the inlet, the values of the dimensionless numbers are: Re = 159, J µ = 3.87 × 10 −4 and Pr = 971.Here the characteristic length H r is chosen as the diameter of the pipe, D. The length of the pipe is L = 100D, for ensuring that the flow is fully developed.
Figure 3a shows the relative volume fraction profiles, defined as φ/φ, for different values of the bulk solid volume fraction, φ.The volume fraction distribution is more uniform when φ is large.In our model, when the flow is isothermal the non-uniform distribution of particles is mainly attributed to the collision rate of the particles (related to the shear rate), corresponding to the term N c in Equation (15).We expect that if the φ approaches the maximum packing value, the volume fraction distribution will almost be constant.At the same time, for higher values of φ, the concentration is more uniform near the pipe center.Figure 3b shows the velocity distributions with different values of φ.The velocity profiles become blunter as φ increases, perhaps due to the higher concentration near the center of the pipe.
In drilling applications, fresh water is pumped into the well to remove the solid cuttings.The cost of the pumping power, is related to the pressure drop.Figure 3c shows the pressure drop between the pipe inlet and the outlet for different values of φ.As φ increases, the pressure drop increases almost exponentially.Therefore, it is important to control or monitor the values of φ in order to have a more efficient drilling process; this can be accomplished by increasing the flow rate of freshly pumped water or by decreasing the rate of the penetration (ROP), where ROP is defined as the speed at which a drill bit breaks the rock to deepen the borehole [77].The effect of the Reynolds number is negligible and will not be studied here [12,29].
between the pipe inlet and the outlet for different values of .As increases, the pressure drop increases almost exponentially.Therefore, it is important to control or monitor the values of in order to have a more efficient drilling process; this can be accomplished by increasing the flow rate of freshly pumped water or by decreasing the rate of the penetration (ROP), where ROP is defined as the speed at which a drill bit breaks the rock to deepen the borehole [77].The effect of the Reynolds number is negligible and will not be studied here [12,29].

Effect of Temperature in the Axial Direction
It is known the temperature of the Earth rises about 2.5 • C per 100 meters of depth [3], although in some situations this may vary depending on the geographical location.Therefore, the effects of the earth temperature variation along the axial (vertical) direction need to be studied.Mathematically the change in the earth temperature can be represented as a variable temperature boundary condition on the outside wall of the pipe.It is a computationally intensive calculation if we are to simulate a practical pipe with length of several kilometers.Therefore, here we perform the simulations for a pipe with length of 1200D, with a large temperature gradient along the axial direction; specifically, we assume 0.05 • C per H r increases.θ at the inlet is assumed to be zero.We also assume that the value of the bulk volume fraction is 0.2.The values of the dimensionless numbers for this case are based on the inlet conditions: Re = 159, J µ = 3.86 × 10 −4 and Pr = 971.
Figure 4a shows the radial (Y-direction, see Figure 2) temperature distribution at different depths.Near the walls the temperature gradient is large; overall, the temperature increases as the Y-position gets closer to the hot wall, while near the center, the temperature is the same as the inlet temperature (θ = 0).Figure 4b shows that as the depth (in the axial direction) increases, that is, as the wall temperature increases, the volume fraction turns out to be more non-uniform and more particles move toward the channel wall.This can be attributed to the decreasing of the viscosity near the hot wall.From the particle transport model, namely Equations ( 14)-( 16), we can see that N µ causes the particles to move from higher viscosity regions to regions with lower viscosity.This is perhaps the reason for the volume fraction near the wall being higher at larger depths.Similar to Figure 3b, Figure 4c shows that the velocity profiles are not affected much.Another important parameter is the pressure distribution in the axial direction.As the depth increases, the pressure gradient decreases a little, perhaps due to the decrease of the suspension viscosity, see Figure 4d.
Figure 4a shows the radial (Y-direction, see Figure 2) temperature distribution at different depths.Near the walls the temperature gradient is large; overall, the temperature increases as the Yposition gets closer to the hot wall, while near the center, the temperature is the same as the inlet temperature ( ̅ = 0).Figure 4b shows that as the depth (in the axial direction) increases, that is, as the wall temperature increases, the volume fraction turns out to be more non-uniform and more particles move toward the channel wall.This can be attributed to the decreasing of the viscosity near the hot wall.From the particle transport model, namely Equations ( 14)-( 16), we can see that causes the particles to move from higher viscosity regions to regions with lower viscosity.This is perhaps the reason for the volume fraction near the wall being higher at larger depths.Similar to Figure 3b, Figure 4c shows that the velocity profiles are not affected much.Another important parameter is the pressure distribution in the axial direction.As the depth increases, the pressure gradient decreases a little, perhaps due to the decrease of the suspension viscosity, see Figure 4d.

Flow between Two Eccentric Vertical Cylinders
In geothermal applications, the flow of the drilling suspension is usually represented by the simple geometry of the flow between two vertical cylinders with the inner cylinder rotating [78].Figure 5 shows the geometry where the radius of the outer cylinder is 2.5 times of the radius of the inner cylinder; ε = e/(R − R i ) = 0.5 the eccentricity ratio, where e is the distance between the centers of the two cylinders.Here the hydraulic diameter of the channel, 2(R o − R i ), is chosen as the characteristic length, H r .The length of the cylinders is assumed to be 100H r , and the bulk volume fraction is 0.1.All the figures in the remainder of this paper are plotted along the line A-A as shown in Figure 5.
Figure 5 shows the geometry where the radius of the outer cylinder is 2.5 times of the radius of the inner cylinder; = /( − ) = 0.5 is the eccentricity ratio, where is the distance between the centers of the two cylinders.Here the hydraulic diameter of the channel, 2( − ), is chosen as the characteristic length, .The length of the cylinders is assumed to be 100 , and the bulk volume fraction is 0.1.All the figures in the remainder of this paper are plotted along the line A-A as shown in Figure 5. is the distance between the center of the two cylinders.

Effect of (Rotational Speed)
An important parameter in drilling applications is the rotational speed of the inner cylinder.This determines the rate of penetration (ROP) [77], responsible for the amount of the solid cuttings produced, which in our problems is the bulk volume fraction.Clearly, the rotational speed, Ω, also affects the flow field.For the current problem, we assume isothermal conditions, i.e., ̅ is zero.Based on the inlet conditions, the values of the dimensionless numbers are: Figure 6 shows the streamlines for different values of .As increases, the streamlines become more distorted.Figure 7 shows the volume fraction fields at the Z-Y plane at different axial positions (X coordinate) for different values of .For the cases shown in Figure 7 the flow is steadystate.With no rotation, the volume fraction distribution is symmetric along the A-A plane.The particles tend to accumulate near the center line, also see Figure 8a; this agrees with the results predicted by the multicomponent model of [79,80].As increases, the particles tend to migrate in the direction of the angular speed of the inner cylinder; furthermore, we also notice that the high concentration of the solid particles in the narrow gap gradually begins to become more uniform, see Figure 8a as well.Figure 8b shows that except the region near the rotating cylinder, the change in the magnitude of the velocity profiles is moderate for different values of .An important parameter in drilling applications is the rotational speed of the inner cylinder.This determines the rate of penetration (ROP) [77], responsible for the amount of the solid cuttings produced, which in our problems is the bulk volume fraction.Clearly, the rotational speed, Ω, also affects the flow field.For the current problem, we assume isothermal conditions, i.e., θ is zero.Based on the inlet conditions, the values of the dimensionless numbers are: Re = 191.50,J µ = 2.69 × 10 −4 and Pr = 971.For this flow situation, we define the Reynolds number as, Re r = ρ f 0 H r (ΩR i )/µ r , which is related to the rotational speed of the inner cylinder.The values of Re r are 0 (Ω = 0 RPM), 24.05 (Ω = 60 RPM) and 72.16 (Ω = 180 RPM).
Figure 6 shows the streamlines for different values of Re r .As Re r increases, the streamlines become more distorted.Figure 7 shows the volume fraction fields at the Z-Y plane at different axial positions (X coordinate) for different values of Re r .For the cases shown in Figure 7 the flow is steady-state.With no rotation, the volume fraction distribution is symmetric along the A-A plane.The particles tend to accumulate near the center line, also see Figure 8a; this agrees with the results predicted by the multicomponent model of [79,80].As Re r increases, the particles tend to migrate in the direction of the angular speed of the inner cylinder; furthermore, we also notice that the high concentration of the solid particles in the narrow gap gradually begins to become more uniform, see Figure 8a as well.Figure 8b shows that except the region near the rotating cylinder, the change in the magnitude of the velocity profiles is moderate for different values of Re r .

Effects of the Earth Temperature
The wall temperature (or the outer cylinder) is assumed to be the same as the Earth temperature surrounding the cylinders.The value of θ at the inlet and at the wall of the inner cylinder is kept constant as zero, and Re r = 36.08(Ω = 90 RPM).Based on the inlet conditions, the values of the dimensionless numbers for these simulations are: Re = 191.50,J µ = 2.69 × 10 −4 and Pr = 971.
Figure 9 shows the profiles for the particles concentration, the temperature and the temperature gradient in the Z-Y plane at different axial positions and for different values of the outer wall temperature.By comparing Figure 9a with Figure 7, we notice that when the Earth (the surrounding) temperature is high, a thin layer of high particle concentration near the outer wall of the cylinder is formed (also see Figure 10a); this layer may be due to the non-uniform distribution of the suspension viscosity caused by the large temperature variation.This thin layer with a particle solid concentration may further enhance particles deposition near the outer cylinder.From Figure 9b, we can see that heat transfers from the outer cylinder to the interior of the flow occurs gradually; near the outer cylinder the temperature gradient is much larger than the interior regions of the flow, also see Figure 10c; in addition, in the wide gap region, we can observe high temperatures which could be attributed to the convection caused by the rotation of the inner cylinder.Figures 9c and 10d indicate that the temperature gradient in the narrow gap is higher than the wide gap, which could be attributed to the stronger convective heat transfer in the narrow gap due to the inner cylinder rotation.The heat loss in the radial direction through the inner cylinder is not so clearly obvious; this can perhaps be attributed to the relatively high convective heat transfer along the axial direction.Figure 10d shows the velocity profiles for different values of the outer wall temperatures.As the wall temperature increases, in the wide gap region the location of the maximum velocity moves towards the outer cylinder.
Energies 2017, 10, 1349 13 of 17 The wall temperature (or the outer cylinder) is assumed to be the same as the Earth temperature surrounding the cylinders.The value of ̅ at the inlet and at the wall of the inner cylinder is kept constant as zero, and = 36.08(Ω = 90 RPM).Based on the inlet conditions, the values of the dimensionless numbers for these simulations are: = 191.50,= 2.69 × 10 and = 971.
Figure 9 shows the profiles for the particles concentration, the temperature and the temperature gradient in the Z-Y plane at different axial positions and for different values of the outer wall temperature.By comparing Figure 9a with Figure 7, we notice that when the Earth (the surrounding) temperature is high, a thin layer of high particle concentration near the outer wall of the cylinder is formed (also see Figure 10a); this layer may be due to the non-uniform distribution of the suspension viscosity caused by the large temperature variation.This thin layer with a particle solid concentration may further enhance particles deposition near the outer cylinder.From Figure 9b, we can see that heat transfers from the outer cylinder to the interior of the flow occurs gradually; near the outer cylinder the temperature gradient is much larger than the interior regions of the flow, also see Figure 10c; in addition, in the wide gap region, we can observe high temperatures which could be attributed to the convection caused by the rotation of the inner cylinder.Figure 9c and Figure 10d indicate that the temperature gradient in the narrow gap is higher than the wide gap, which could be attributed to the stronger convective heat transfer in the narrow gap due to the inner cylinder rotation.The heat loss in the radial direction through the inner cylinder is not so clearly obvious; this can perhaps be attributed to the relatively high convective heat transfer along the axial direction.Figure 10d shows the velocity profiles for different values of the outer wall temperatures.As the wall temperature increases, in the wide gap region the location of the maximum velocity moves towards the outer cylinder.

Concluding Remarks
In this study, we look at the heat transfer in a drilling fluid.The fluid is assumed to be a suspension composed of various components with complex rheological behavior.We model this suspension as a non-linear fluid, where the viscosity and the heat diffusivity depend on the concentration of the particle, the shear rate and temperature.The motion of the particles is based on a particle flux equation, where a generalization of the model proposed by [29] is used.By performing numerical studies for two different flow conditions, namely the flow in a vertical pipe and flow between two rotating cylinders with various boundary and flow conditions, we observe that the temperature and the flow fields depend significantly on the earth (surrounding) temperature distribution, the bulk volume fraction, etc.We also notice that for a large value of the (dimensionless) earth temperature, such as ̅ = 1.0, a thin layer with high particle concentration is formed near the outer cylinder, which may enhance particle deposition in that region.Furthermore, as the bulk volume fraction increases, the pressure drop increases exponentially, perhaps due to the increase in the suspension viscosity.Finally, in certain applications, especially those involving gas-particle flows, such as fluidized bed reactors, under certain conditions, there could be a counter flow situations due to high solid loading (see Bednarz, et al. [81], Das et al. [82], Garić-Grulović et al. [83]).We did not consider this phenomenon in this paper.

Concluding Remarks
In this study, we look at the heat transfer in a drilling fluid.The fluid is assumed to be a suspension composed of various components with complex rheological behavior.We model this suspension as a non-linear fluid, where the viscosity and the heat diffusivity depend on the concentration of the particle, the shear rate and temperature.The motion of the particles is based on a particle flux equation, where a generalization of the model proposed by [29] is used.By performing numerical studies for two different flow conditions, namely the flow in a vertical pipe and flow between two rotating cylinders with various boundary and flow conditions, we observe that the temperature and the flow fields depend significantly on the earth (surrounding) temperature distribution, the bulk volume fraction, etc.We also notice that for a large value of the (dimensionless) earth temperature, such as θ wo = 1.0, a thin layer with high particle concentration is formed near the outer cylinder, which may enhance particle deposition in that region.Furthermore, as the bulk volume fraction increases, the pressure drop increases exponentially, perhaps due to the increase in the suspension viscosity.Finally, in certain applications, especially those involving gas-particle flows, such as fluidized bed reactors, under certain conditions, there could be a counter flow situations due to high solid loading (see Bednarz,et al. [81], Das et al. [82], Garić-Grulović et al. [83]).We did not consider this phenomenon in this paper.

Figure 1 .
Figure1.Dependency of the suspension viscosity on the temperature.The bulk volume fraction is 0.345 and the density of the drilling mud is 2.2 g/cm3 .Experimental data is from Wang et al.[16].

Figure 1 .
Figure1.Dependency of the suspension viscosity on the temperature.The bulk volume fraction is 0.345 and the density of the drilling mud is 2.2 g/cm3 .Experimental data is from Wang et al.[16].

Figure 2 .
Figure 2. Schematic of the pipe.

Figure 2 .
Figure 2. Schematic of the pipe.

Figure 3 . 2 .Figure 3 .
Figure 3. (a) The relative volume fraction profiles along the radial direction (Y direction as shown in Figure 2.) for different values of the bulk volume fraction, ; (b) The profiles of the axial (X-direction) velocity, , along the radial direction; (c) pressure drop between the inlet and the outlet for different values of the bulk volume fraction .= 159, = 3.87 × 10 and = 971.

Figure 4 .
Figure 4. (a) Temperature distribution; (b) volume fraction; and (c) the velocity at different axial positions with the wall temperature (surrounding earth temperature) varying along the axial (X) direction; (d) the pressure distribution at different axial positions with the wall temperature varying along the axial direction.The bulk volume fraction is 0.2, inlet ̅ is 0, = 159, = 3.86 × 10 and = 971.

Figure 4 .
Figure 4. (a) Temperature distribution; (b) volume fraction; and (c) the velocity at different axial positions with the wall temperature (surrounding earth temperature) varying along the axial (X) direction; (d) the pressure distribution at different axial positions with the wall temperature varying along the axial direction.The bulk volume fraction is 0.2, inlet θ is 0, Re = 159, J µ = 3.86 × 10 −4 and Pr = 971.

Figure 5 . 2 . 5 .
Figure 5. Geometry of the two eccentric cylinders.We assume is the length of the simulated portion of the pipe, which is 100 = 200( − ).The radius of the outer cylinder is 2.5 times of the inner, namely = 2.5 .The eccentricity ratio is = /( − ) = 0.5. is the distance between the center of the two cylinders.

Figure 5 .
Figure 5. Geometry of the two eccentric cylinders.We assume L is the length of the simulated portion of the pipe, which is 100H r = 200(R o − R i ).The radius of the outer cylinder is 2.5 times of the inner, namely R o = 2.5R i .The eccentricity ratio is ε = e/(R o − R i ) = 0.5.e is the distance between the center of the two cylinders.

Figure 6 .
Figure 6.Streamlines for the isothermal flow between the two eccentric cylinders with different values of (related to the rotational speed).

Figure 7 .
Figure 7. Top view.The volume fraction fields in the Z-Y plane for different dimensionless axial positions (X-direction) for different (rotational speed).The rotation is counter clockwise.

Figure 8 .Figure 6 . 17 Figure 6 .
Figure 8.(a) The volume fraction and (b) the velocity profiles along the A-A plane (see Figure 5 also) near the outlet for different .4.2.2.Effects of the Earth Temperature

Figure 7 .
Figure 7. Top view.The volume fraction fields in the Z-Y plane for different dimensionless axial positions (X-direction) for different (rotational speed).The rotation is counter clockwise.

Figure 8 .Figure 7 . 17 Figure 6 .
Figure 8.(a) The volume fraction and (b) the velocity profiles along the A-A plane (see Figure 5 also) near the outlet for different .4.2.2.Effects of the Earth Temperature

Figure 7 .
Figure 7. Top view.The volume fraction fields in the Z-Y plane for different dimensionless axial positions (X-direction) for different (rotational speed).The rotation is counter clockwise.

Figure 8 .Figure 8 .
Figure 8.(a) The volume fraction and (b) the velocity profiles along the A-A plane (see Figure 5 also) near the outlet for different .4.2.2.Effects of the Earth Temperature

Figure 9 .
Figure 9. Top view.(a) The volume fraction and (b) the temperature (c) the temperature gradient fields in the Z-Y plane at different axial positions and different outer-wall temperatures, when = 36.8(Ω = 90 RPM), and ̅ at the inlet and at the inner cylinder wall is 0. The unit of the scale bar for the temperature gradient is 10 .

Figure 9 .
Figure 9. Top view.(a) The volume fraction and (b) the temperature (c) the temperature gradient fields in the Z-Y plane at different axial positions and different outer-wall temperatures, when Re r = 36.8(Ω = 90 RPM), and θ at the inlet and at the inner cylinder wall is 0. The unit of the scale bar for the temperature gradient is 10 4 .

Figure 10 .
Figure 10.(a) The volume fraction profiles; (b) the velocity profiles; (c) the temperature profiles and (d) the temperature gradient profiles for different values of the earth temperature (wall temperature of the outer cylinder), when = 36.8(Ω = 90 RPM) and ̅ at the inlet and the inner cylinder wall is 0.

Author Contributions:Figure 10 .
Figure 10.(a) The volume fraction profiles; (b) the velocity profiles; (c) the temperature profiles and (d) the temperature gradient profiles for different values of the earth temperature (wall temperature of the outer cylinder), when Re r = 36.8(Ω = 90 RPM) and θ at the inlet and the inner cylinder wall is 0.

Table 1 .
Physical properties used in this paper