Natural Convection from a Heated Sphere in Bingham Plastic Fluids

Article pubs.acs.org/IECR

Natural Convection from a Heated Sphere in Bingham Plastic Fluids N. Nirmalkar,† A. K Gupta, and R. P. Chhabra* Department of Chemical Engineering, Indian Institute of Technology, Kanpur-208016, India ABSTRACT: Laminar natural convection from an isothermal sphere immersed in cold Bingham plastic media has been investigated numerically over wide ranges of the pertinent dimensionless parameters (Rayleigh number, 102 ≤ Ra ≤ 106; Prandtl number, 10 ≤ Pr ≤ 100; and Bingham number, 0 ≤ Bn ≤ 104). Extensive results on the flow and heat-transfer aspects are presented in terms of the streamlines and isotherm contours in the close proximity of the sphere; morphology of the yielded and unyielded regions; and the local and average Nusselt number as functions of the Rayleigh number, Prandtl number, and Bingham number. For a given value of the Rayleigh number, there exists a limiting Bingham number beyond which no convection occurs and heat is transferred from the sphere to the fluid only by conduction. Broadly, fluidlike regions are seen to gradually diminish with the increasing Bingham number, and eventually the flow field is completely frozen, corresponding to the fully plastic limit. Finally, the present numerical results of the average Nusselt number are correlated via a simple expression, thereby enabling the interpolation of the present results for the intermediate values of the parameters.

1. INTRODUCTION Most structured fluids like foams, emulsions, suspensions, polymeric systems, and micellar solutions encountered in scores of industrial settings exhibit a range of non-Newtonian characteristics including shear-dependent viscosity, yield stress, visco-elasticity, etc.1−3 Over the years, significant research effort has focused on developing fundamental understanding and reliable design strategies for processing such materials. A cursory examination of the contemporary literature reveals that a great majority of such studies are concerned with their momentum-transfer characteristics in ducts and pipes,1−3 mixing vessels,4,5 porous media flows,6,7 and, in recent years, with the external boundary-type flows such as over a sphere7 or two-dimensional cylinders of various cross sections.8 Furthermore, within the framework of non-Newtonian fluids, much of the literature is based on the use of the simple power-law model which captures only the shear-dependent viscosity over a finite range of shear rate, e.g., see refs 1, 3, 7, and 8. Even for the case of power-law fluids, heat-transfer aspects have been studied only recently and infrequently.7,8 Thus, for instance, the first numerical study of forced convection from an unconfined sphere was reported by Dhole et al.9 and for a confined sphere by Song et al.10,11 even more recently. Subsequently, the steady mixed convection from a sphere in power-law fluids in the aiding buoyancy configuration has been investigated numerically12 and the free convection limit for a sphere has also been studied elsewhere.13 While detailed discussion of the previous pertinent works as well as that of the new numerical results is available in our recent papers,9−13 suffice it to say here that shear-thinning fluid behavior can enhance the rate of heat transfer from a sphere up to 70−80% over and above that in Newtonian fluids otherwise under identical conditions and that shear-thickening fluid behavior has a weak adverse effect on heat transfer. The degree of enhancement (or reduction), however, is strongly dependent on the values of the Reynolds number and Prandtl number in the forced convection, on Grashof and Prandtl numbers in the free convection, and additionally on the Richardson number in the mixed convection regimes. Two © 2014 American Chemical Society

further observations are in order here. First, the preceding numerical results9−13 for the forced and free convection regimes approach the approximate boundary layer predictions of Acrivos and co-workers14,15 and of others16 at large values of the Reynolds number or Grashof number. Second, these predictions are also consistent with the scant experimental results for spheres in power-law fluids.9,10,17,18 It is now possible to estimate the value of the mean Nusselt number for an unconfined sphere in power-law fluids over most ranges of conditions of practical interest, although this body of knowledge is nowhere as extensive as that in Newtonian fluids.19,20 On the other hand, many such fluids also possess the socalled yield stress.21,22 These materials deform like an elastic solid for the stress levels below its yield stress, and once the applied stress exceeds the yield stress, it may exhibit a constant viscosity (Bingham plastics) or a shear-thinning viscosity (Herschel−Bulkley model fluid). Because of their dual nature, the mixing and homogenization as well as the heating and cooling of such fluids are far more difficult than those of the low-viscosity Newtonian and power-law fluids.3−5,23 This is simply due to the coexistence of fluidlike (yielded) and solidlike (unyielded) subregions depending upon the prevailing stress levels vis-à-vis the fluid yield stress in a given application. This difficulty is further accentuated when heat transfer occurs solely by natural convection. In this regime, the buoyancy-induced flow promotes heat transfer in such systems, whereas the effective fluid viscosity (including the contribution of yield stress) counters such a weak flow. Therefore, it is conceivable that for a given value of the fluid yield stress, there must exist a threshold buoyancy-induced flow for the onset of natural convection in such fluids.23 Conversely, in the limit of large yield stress effects, no yielding will occur and thus conduction will be the sole mechanism of heat transfer. A heated sphere Received: Revised: Accepted: Published: 17818

August 7, 2014 October 16, 2014 October 20, 2014 October 20, 2014 dx.doi.org/10.1021/ie503152k | Ind. Eng. Chem. Res. 2014, 53, 17818−17832

Industrial & Engineering Chemistry Research

Article

from a sphere also denotes an idealization of several industrially important applications including fluidized and fixed bed reactors, reheating of processed foodstuffs, manufacturing of pharmaceutical and personal-care products, melting of polymer pellets, etc. Notwithstanding the importance of the detailed velocity and temperature profiles in the close proximity of the heated sphere, it is readily conceded that the value of the Nusselt number as a function of the pertinent parameters is the variable of central interest in process engineering design calculations. It is appropriate to add here that while no prior results exist for the laminar natural convection for a sphere in Bingham plastic media, limited results are available for a circular cylinder in the free convection regime24 and for a triangular cylinder in the forced convection regime.25 On the other hand, there have been many numerical and experimental studies for a sphere in the forced convection regime, although most of these are concerned with the drag behavior, wall effects, prediction of yield-surfaces, etc., and only two studies deal with the corresponding heat-transfer aspects in the steady axisymmeric regime. Most of these studies have been reviewed recently.7,26−30 In this regime also, the numerical drag data is consistent with the scant experimental results available in the literature. On the basis of the aforementioned discussion, it is thus clear that little is known about natural convection heat transfer in Bingham plastic fluids in general and from a heated sphere in particular. This work aims to fill this gap in the current body of knowledge. In particular, the governing partial differential equations have been solved numerically to elucidate the effects of Rayleigh number (102 ≤ Ra ≤ 106), Prandtl number (10 ≤Pr ≤ 100), and Bingham number (0 ≤ Bn ≤ 104) on the detailed velocity and temperature fields as well as on the morphology of the yielded and unyielded regions and local and average Nusselt number for an isothermal sphere. The present results are compared with the previous studies wherever possible. Figure 1. Schematic diagram of computational domain.

2. PROBLEM DEFINITION AND GOVERNING EQUATIONS Consider a heated sphere (of diameter d) whose surface is maintained at a constant temperature Tw submerged in a quiescent

offers a convenient prototype to study laminar natural convection in Bingham plastic fluids. In addition, heat transfer

Figure 2. Effect of computational grid on the variation of the local Nusselt number over the surface of sphere at Ra = 106 and Pr = 100 for (a) Bn = 0 and (b) Bn = 104. 17819

dx.doi.org/10.1021/ie503152k | Ind. Eng. Chem. Res. 2014, 53, 17818−17832


Article

Table 1. Effect of Grid Size on the Drag Coefficients and Average Nusselt Number Ra = 106, Pr = 100 CD

CDP 4

Nuavg 4

grid

elements

Bn = 0

Bn = 10

Bn = 0

Bn = 10

Bn = 0

Bn = 104

G1 G2 G3

23 400 51 088 102 600

0.2920 0.2918 0.2917

176.32 176.25 176.23

0.1104 0.1106 0.1107

58.650 58.728 58.812

19.868 19.949 19.971

2.0039 2.0126 2.0151

Bingham plastic fluid at temperature, T0 ( τ02 (9a)



Article

Table 2. Comparison of the Present Results with the literature values in Newtonian Fluids. Nuavg GrR Pr = 10 103 105 Pr = 10 103 105 Pr = 10 103 105 Pr = 10 103 105 Pr = 10 103 105

ref 39

ref 38

ref 13

ref 34

present work

3.2684 6.011 14.684

3.26 6 14.66

2.9601 5.668 14.191

2.93 5.65 14.22

2.927 5.602 14.007

4.6655 10.429 28.655

4.63 10.33 28.34

4.383 10.017 27.901

4.32 10.01 28.74

4.3854 10.014 27.526

5.6004 13.385 38.004

5.54 13.22 37.49

5.3537 12.913 37.064

− − −

5.3683 12.934 36.947

6.6207 16.612 48.207

6.55 16.38 47.47

6.4412 16.1 47.112

− − −

6.4505 16.125 47.036

7.5509 19.553 57.509

7.45 19.26 56.59

7.4203 19.011 56.28

− − −

7.4401 19.040 56.112

0.72

7

20

50

100

Figure 5. Comparison of the present temperature and velocity profiles in water (Pr = 7) along θ = 90 ° line at Ra = 1.7 × 108 with the literature values.

Thus, for a Bingham plastic fluid, the scalar viscosity function (η) is of the following form:

Figure 4. Comparison of the present Nusselt number results (in the form of Nu/Ra1/4) over the surface of the sphere with that of Jia and Gogos.34

γ ̇ = 0,

if IIτ ≤ τ02

η = μB +

(9b)

Clearly, eq 9 is inherently discontinuous and nondifferentiable and hence not convenient for a direct implementation in a numerical scheme. This difficulty is circumvented here by employing the regularization method proposed by Papanastasiou32 which provides a smooth variation of the absolute viscosity between the unyielded and yielded regions of the fluid. Within this framework, eq 9 is rewritten as follows: ⎛ τ0[1 − exp(− m IIγ ̇ )] ⎞ ⎟γ ̇ τ = ⎜⎜μB + ⎟ IIγ ̇ ⎝ ⎠

τ0[1 − exp(−m IIγ ̇ )] IIγ ̇

(11)

Clearly, in the limit of m → ∞, eq 10 reduces exactly to eq 9. On the other hand, m = 0 corresponds to the Newtonian fluid behavior. Hence, large values of m are generally needed to obtain accurate predictions of the velocity and temperature fields. Another regularization method, the so-called biviscous model, which replaces the unyielded regions by an equivalent fluid of very high viscosity (μy) called yielding viscosity, has also been utilized in this work to contrast the results based on these two approaches, whereas a detailed discussion regarding the relative merits and demerits of these approaches is available elsewhere.33

(10)

where m is a regularization parameter. 17821



Article

Figure 6. Representative streamline (right-hand half) and isotherm (left-hand half) contours at Ra = 102. (Flow is in the upward direction.)

The regularized constitutive equation for a biviscous fluid is given as τ = μy γ ̇

η = μB +

for IIτ ≤ τ0 2

⎡ ⎤ 1 ⎥ τ = τ0⎢1 − + μB γ ̇ ⎢⎣ μy /μB ⎥⎦

(12a)

for IIτ > τ0

⎡ τ0⎢1 − ⎣

1 ⎤ μy / μB ⎥ ⎦

IIγ ̇

(13)

The problem statement closure is obtained by identifying the physically realistic boundary conditions. In this work, the standard boundary conditions of no-slip and T = Tw on the surface of the sphere are prescribed. On the fictitious computational domain, zero values for both the radial velocity and its gradient in the radial direction and T = T0 are prescribed. To enhance the utility of the new results reported herein, these

2

(12b)

Therefore, the scalar viscosity function for an ideal Bingham plastic fluid regularized using the biviscous approach can be written as 17822



Article

Figure 7. Representative streamline (right-hand half) and isotherm (left-hand half) contours at Ra = 106. (Flow is in the upward direction.)

where Gr is the usual Grashof number defined as

were rendered dimensionless by using the following scaling variables: diameter of the sphere d for all length variables, Uc for characteristic or reference velocity ((dgβΔT)1/2) for all velocities, ρ0Uc2 for pressure, and (μBUc/d) for stress components, etc. The fluid temperature is nondimensionalized as ξ = (T − T0)/(Tw − T0). The scaling considerations suggest the coupled velocity and temperature fields to be functions of three dimensionless parameters which are defined here as follows: Bingham number

Rayleigh number

τ Bn = 0 μB Ra =

d gβ ΔT

ρ0 2 c pgβ ΔTd3 μB k

Gr =

ρ0 2 gβ ΔTd3 μB2

Prandtl number

(16)

Pr =

c pμB k

(17)

The numerical solution of the preceding system of equations maps the flow domain in terms of the primitive variables (velocity, temperature, and pressure) which can be postprocessed to evaluate the derived quantities like yield surfaces and local and average Nusselt numbers as described in some of our recent studies.24,25 In brief, the local Nusselt number is defined as

(14)

= Gr ·Pr (15) 17823



Article

Figure 8. Temperature profiles along the x-coordinate.

Nu =

∂ξ hd =− ∂ns k

mesh was manually refined in this region. The yield surfaces are located by the von Mises criterion of plasticity theory31 (within the relative tolerance of 10−6), which states that yielding in the material will occur only at points where the IIτ ≥ τ02 condition is satisfied. Furthermore, the steady, axisymmetric, nonisothermal flow module is selected to obtain the numerical solution. A user defined function (UDF) is introduced to incorporate the buoyancy force and the regularized expression for the absolute viscosity of the fluid. PARDISO direct solver is used to solve the system of equations. A relative tolerance criterion of 10−5 for the equations of motion and energy was used, and further reduction in this value had a negligible effect on the results. Next, following the strategy used in our previous studies,24,25 the size of the computational domain (i.e., value of D∞) was systematically varied as 40, 80, and 200. A detailed examination of the resulting values of drag coefficient and Nusselt number at the lowest values of Ra, Pr, and Bn revealed the value of D∞ = 80 to be adequate over the range of conditions spanned here. On the other hand, because the boundary layers are expected to be very thin at high values of the Rayleigh number and/or Prandtl number and/or Bingham number, grid independence

(18)

The surface mean value of the Nusselt number is obtained by simply integrating the local values of the Nusselt number over the surface of the sphere.

3. NUMERICAL METHODOLOGY AND CHOICE OF NUMERICAL PARAMETERS The aforementioned governing equations have been solved together with the regularized Bingham constitutive relation using the finite element-based solver Comsol Multiphysics (version 4.3). Because detailed discussion of the numerical solution procedure is available elsewhere,24−28 only the salient features are recapitulated here. In brief, the velocity and temperature gradients are expected to be sharp in the vicinity of the sphere. A relatively fine mesh is used in this region; hence, the computational domain is mapped with quadrilateral elements with nonuniform spacing. Similarly, a fine mesh is also needed in the neighborhood of the yield surfaces delineating the fluidlike and solidlike regions. Initially, the approximate location of the yield surface was established and then the computational 17824



Article

Figure 9. Variation of the surface vorticity (ω) over the sphere.

obtained over a wide range of Rayleigh number, 102 ≤ Ra ≤ 106; Prandtl number, 10 ≤ Pr ≤ 100; and Bingham number, 0 ≤ Bn ≤ 104. However, prior to the detailed presentation and discussion of the new results, it is desirable to establish the accuracy and precision of the new results obtained herein. This objective is accomplished here by comparing the present results with the available benchmark studies in the literature. 4.1. Validation of Results. Excellent numerical and experimental results are available in the literature on the laminar free convection from a heated sphere in Newtonian media.13,34−37 Table 2 shows a comparison in terms of the average Nusselt number with the previous numerical predictions13,34 and the experimental results as approximated here by the correlations of Churchill and Churchill38 and Jafarpur and Yovanovich.39 A good correspondence is seen to exist between various predictions based on different numerics as well as between the predictions and experimental results in Newtonian fluids. To add further weight to this assertion, Figures 4 and 5 show comparisons in terms of the numerical values of the local Nusselt number variation over the surface of the sphere and the detailed velocity and experimental temperature profiles in water. Included in these figures are also the other predictions available in the literature. Overall, the present results are seen to be in fair agreement with experiments in Newtonian fluids. Aside from the preceding validation in Newtonian fluids,

tests have been carried out at the maximum values of these parameters, i.e., Ra = 106, Pr = 100, and Bn = 0 and Bn = 104. For this purpose, three grids, namely, G1, G2 and G3 were created as detailed in Table 1. Evidently, the resulting values of the drag coefficient and Nusselt number obtained with grid G2 and G3 are seen to be very close to each other. Furthermore, Figure 2 shows the effect of the computational mesh in terms of the distribution of the local Nusselt number on the surface of the sphere. On the basis of the results shown in Table 1 and Figure 2, grid G2 is considered to be suitable over the range of conditions spanned in this study. Finally, one must also choose an appropriate value of the regularization parameter (m or μy) also. Figure 3 shows the sensitivity of the local Nusselt number to the value of the regularization parameter m. The results for m* = 108 and m* = 109 are virtually indistinguishable from each other; therefore, the value of m* = 108 is used here to obtain the new results for laminar free convection from an isothermal sphere. In summary, the numerical results reported herein are based on the parameters D∞ = 80, grid G2, and m* = 108. Further justification for their suitability is provided by performing a few benchmark comparisons presented in the ensuing section.

4. RESULTS AND DISCUSSION Extensive numerical results on the free convective flow and heat-transfer characteristics of an isothermal sphere have been 17825



Article

Figure 10. Typical structure of fluidlike (unshaded) and solidlike (shaded) regions at (a) Ra = 102 and (b) Ra = 106. (Dashed lines represent predictions of biviscous model, and the flow is in the upward direction).

the gravity effects and it is replaced by the cold fluid from sides and/or from beneath the sphere depending upon the value of the Rayleigh number. Figures 6 and 7 show the typical streamline and isotherm contours for a “weak” (Ra = 102) and “strong” (Ra = 106) buoyancy-induced flow current for a range of values of the Prandtl and Bingham numbers. At Ra = 102, it is clearly seen that very little yield stress is required to suppress the fluid deformation. Thus, for instance, in the limiting case of

numerous benchmark comparisons for free convection in Bingham plastic fluids in a square cavity and forced convection flow over a sphere and a cylinder have been presented elsewhere;24−27 hence, these are not repeated here. On the basis of these comparisons, the new results reported herein are believed to be reliable to within approximately ±1.5−2%. 4.2. Streamline and Isotherm Contours. Qualitatively, the heated mass of fluid near the sphere rises upward because of 17826



Article

limited to the yielded fluidlike regions whereas conduction is the principal mechanism of heat transfer in the unyielded solidlike subdomains. It thus stands to reason that the conduction may be the overall heat-transfer limiting step in such systems. It is thus imperative to understand the influence of the Rayleigh number, Prandtl number, and Bingham number on the morphology (shape and size) of the yielded and unyielded subregions in the present case. Figure 10 shows representative results depicting the influence of the Prandtl number and Bingham number on the yield surfaces for the two extreme values of the Rayleigh number, namely, Ra = 102 and Ra = 106. Intuitively, it is reasonable to postulate that the fluid yielding should be promoted by the increasing values of Rayleigh number whereas this tendency is suppressed by the increasing Bingham number, with a little influence of the Prandtl number over and above that reflected by the value of the Rayleigh number. Indeed, this conjecture is well borne out by the trends seen in Figure 10. Broadly, there are four unyielded masses of fluid adhering at equatorial and polar points on the surface of the sphere at low Bingham numbers. With the increasing Bingham numbers, the unyielded material adhering at θ = π/2 and 3π/2 progressively diminishes whereas that at the front and rear stagnation points grows in size. Eventually, at Bn = 10 for Ra = 100 and at Bn = 103 for Ra = 106, there is very little yielding and the sphere is buried in a frozen fluid, so to speak. Depending upon the values of the Bingham number and Rayleigh number, some intermediate shapes of the yield surface are also observed. Thus, for instance, at Ra = 106 and Pr = 100, the outer yield surface undergoes a shape change from being a nearly concentric sphere to one which is squeezed at the equator and is elongated in the flow direction as Bn increases from Bn = 50 to Bn = 100. On the other hand, at low Rayleigh numbers (Ra = 100), a similar trend is observed at Bn = 5. From the foregoing results, it is evident that for a given Rayleigh number, there exists a limiting Bingham number (Bnmax) beyond which no yielding occurs and heat transfer occurs solely by conduction under these conditions. Figure 11 shows the functional dependence of (Bnmax) on the Rayleigh number and Prandtl number. At low

Newtonian fluid behavior (Bn = 0), the effect of increasing Prandtl number at the expense of decreasing Grashof number (for a fixed value of Rayleigh number) is consistent with the classical boundary layer analysis.40 The cold fluid is entrained mainly from beneath the sphere under these conditions. On the other hand, in Bingham plastic fluids, howsoever small the yield stress, the thickening of the boundary layer is seen to occur as revealed by the streamline patterns at Bn = 0.1 and Bn = 1. At Ra = 106, a strong convection current is present in the fluid and small values of yield stress (such as Bn = 1) exert very little influence on the streamline patterns with reference to that in Newtonian fluids (Bn = 0). Naturally, the severe curvature of the streamlines indicates that cold fluid is increasingly entrained from the lateral sides of the sphere.20 With further increase of Bingham number, the fluidlike regions are seen to attain a limiting value, presumably relating to the fully plastic limit at Bn = 10 in this case. Qualitatively these trends also persist at Ra = 106, with minor differences. Under these conditions, as noted above, the cold fluid is increasingly entrained from the lateral sides and the flow is seen to freeze now (fully plastic limit) at Bn = 50 as opposed to that at Bn = 10 for Ra = 102. The isotherm contours shown in the left-hand half of Figures 6 and 7 conform to the general characteristics of the progressive thinning of the thermal boundary layer with the increasing Rayleigh and/or Prandtl numbers for a given value of the Bingham number. Broadly, conduction dominates the heat transfer at low Rayleigh numbers and/or at high Bingham numbers, as is evident from the shapes of the isotherm contours which closely follow the surface of the sphere. In other words, the isotherms under these conditions are seen to be more or less concentric circles. Further insights can be gained by examining the spatial decay of the temperature field. Figure 8 shows the temperature profiles across the horizontal central line for scores of values of Bn, Pr, and Ra. At low Rayleigh numbers, the temperature decay is seen to occur gradually; the larger the value of the Bingham number, the slower the temperature decay because of the diminishing size of the yielded cavity encapsulating the heated sphere. On the other hand, the trends are qualitatively similar at Ra = 106 except for the faster (spatially) temperature decay which is obviously due to the thinning of the temperature boundary layer. Also, the present predictions are seen to approach the corresponding conduction temperature profiles with the increasing values of the Bingham number, as shown in Figure 8. This section is concluded by presenting representative results on the distribution of the vorticity on the surface of the sphere (Figure 9). Broadly, the surface vorticity increases from its zero value at the front stagnation point up to about θ ≃ 120 − 135 ° depending upon the values of the Rayleigh and Prandtl numbers, and beyond this point, it progressively decreases to zero value again at the rear stagnation point. For a fixed Rayleigh number, the vorticity value is seen to decrease with increasing Prandtl number. On the other hand, for fixed values of Ra and Pr, because of the diminishing fluid deformation, vorticity also decreases rapidly with increasing Bingham number, eventually approaching zero value at the limiting Bingham number, e.g., Bn = 10 for Ra = 100 and Bn = 100 for Ra = 106. 4.3. Morphology of Fluidlike and Solidlike Regions. As noted earlier, the fluid yields only at locations where the magnitude of the local stress tensor exceeds the fluid yield stress. Depending upon the complexity of the flow geometry, it is thus not uncommon to observe multiple yielded regions which may or may not be interconnected. Naturally, convective transport is

Figure 11. Dependence of the maximum Bingham number (Bnmax) on Rayleigh and Prandtl numbers (filled symbols show the fit of eqs 19a and 19b). 17827



Article

Figure 12. Distribution of the local Nusselt number over the surface of the sphere.

Rayleigh numbers (Ra ≤ ∼1000), the limiting Bingham number is not only close to zero, but it is also nearly independent of the Prandtl number. This is possibly so because of weak convection current which causes very little yielding. As the value of the Rayleigh number is gradually increased, the limiting Bingham number increases, but it decreases with the increasing Prandtl number. For a fixed Rayleigh number, the increase in Prandtl number is accompanied by the concomitant decrease in the value of the corresponding Grashof number, which weakens the flow. The following simple equations offer a convenient means to estimate the value of the limiting Bingham number for the intermediate values of the Rayleigh number and Prandtl number: Bnmax = 0.215 + 0.06Ra0.65Pr −0.47

and 20%, respectively. Note that in this work, the values of Bnmax are based on the criterion of 1% of the limiting value of Nu∞ = 2, i.e., the values of Nusselt number ranging from Nu = 1.98 to Nu = 2.02 were considered to correspond to the conduction limit. Finally, because the formation of the yield surfaces is a distinct feature of visco-plastic fluids only, it is worthwhile to corroborate their prediction using the other widely used regularization method, namely, the biviscous model. In this approach, for stress levels below the fluid yield stress, the material is assigned viscosity (μy) several orders of magnitude higher than that observed once the applied stress exceeds the fluid yield stress (μB). Limited simulations were performed by using μy/μB = 105), and these results are shown as broken lines in Figure 9. In most cases, the two predictions are seen to be fairly close thereby inspiring confidence in the reliability of the present results. Before leaving this section, suffice it to add here that this behavior is qualitatively consistent with the widely studied case of free convection in Bingham plastic fluids in square and rectangular cavities with differentially heated side walls, e.g., see refs 41 and 42. However, it needs to be

for 102 ≤ Ra ≤ 104 (19a)

⎛ Ra ⎞0.35 Bnmax = 1.04⎜ ⎟ ⎝ Pr ⎠

for 104 < Ra ≤ 106

(19b)

Eqs 19a and 19b reproduce the present numerical results with mean and maximum deviations of 8.3% and 21.3%, and 6.2% 17828



Article

emphasized here that there will always be a small degree of yielding due to the approximations inherent in most of the available regularization methods including the biviscous and exponential schemes used herein. 4.4. Distribution of Local Nusselt Number. The value of the local Nusselt number on the surface of the heated sphere is determined by the temperature gradient normal to the sphere at that point. The temperature field in close proximity of the sphere is, in turn, influenced by the corresponding velocity field and the formation of the solidlike regions adhering to the sphere. Therefore, it is worthwhile to examine the influence of the Bingham number, Prandtl number, and Rayleigh number on the variation of the local Nusselt number over the surface of the isothermal sphere. Typical results are shown in Figure 12 for a range of values of the Bingham number but only for the extreme values of the Rayleigh and Prandtl numbers. It is useful to begin our discussion with the limiting case of Newtonian fluid behavior, i.e., Bn = 0. In this case, the local Nusselt number is always maximum at the front stagnation point (θ = 0 °) regardless of the values of the Prandtl and Rayleigh number because of the zero thickness of the thermal boundary layer. On the other hand, in the case of Bingham plastic fluids at Ra = 106, the location of the maximum value of the Nusselt number is seen to shift downstream from the front stagnation point somewhere between θ = 0 ° and θ = 90 °. This trend is qualitatively similar to that observed for a sphere in power-law fluids13 in which it is due to the presence of the two competing mechanisms, namely, decreasing temperature gradient and variable viscosity. In contrast, in the present case, this is simply due to the formation of polar caps of the unyielded material; thus, in this region, heat transfer occurs principally by conduction. Furthermore, for both Newtonian and Bingham plastic fluids, the minimum value of the Nusselt number is observed at the rear stagnation point of the sphere (θ = 180 °), and this trend is also in line with that seen in the case of power-law fluids.13 Also, in the present case, the local Nusselt number is seen to drop below the conduction limit (Nu∞ = 2) over a small region in the rear of the sphere, especially at low Grashof numbers and at low Bingham numbers. This is seen to occur close to the conditions of fully plastic flow limit. However, as long as some yielding occurs, the average Nusselt number is always higher than the conduction limit. Also, this phenomenon diminishes with the increasing Rayleigh number. Apart from these features, all else being equal, the value of the local Nusselt number gradually decreases with the increasing Bingham number and eventually levels off to a constant for large values of the Bingham number. This is the conduction limit wherein the entire fluid is frozen. In summary, because of the formation of the stagnant solidlike regions on the top and bottom of the sphere, heat transfer in these regions is somewhat impeded in Bingham plastic fluids in comparison with that in Newtonian fluids. 4.5. Average Nusselt Number. On the basis of the preceding discussion, it appears that the rate of heat transfer is impeded in Bingham plastic fluids with reference to that in Newtonian fluids at the same value of Rayleigh number. In other words, the average Nusselt number for an isothermal sphere in a Bingham plastic fluid decreases from its value in Newtonian fluid (Bn = 0) to the conduction limit (Nu∞ = 2) as the value of the Bingham number is progressively increased. Indeed, this assertion is corroborated by the results shown in Figure 13. At low Bingham numbers, a small effect of Rayleigh number and Prandtl number is still present, which is due to

Figure 13. Functional dependence of the average Nusselt number on Bingham number and Prandtl number at Ra = 102 and Ra = 106.

small unyielded parts being present in the flow domain. A detailed examination of the available results in Newtonian fluids (Bn = 0) suggests very weak dependence of the average Nusselt number on the Prandtl number at fixed values of the Rayleigh number in this region. The trends seen in Figure 13 at low Bingham numbers are also consistent with this behavior. However, the present results approach the Newtonian limit extremely slowly. Also, with the increasing value of the Bingham number, the role of the Prandtl number also flips over at about Bn ≃ 0.1 for Ra = 106 and at Bn ≃ 0.01−0.02 for Ra = 102. With the progressively increasing values of the Bingham number, the Nusselt number shows an inverse dependence on the Prandtl number. This is believed to be due to the diminishing yielded regions (with increasing Bn) which tend to lower the rate of heat transfer and the gradual thinning of the thermal boundary layer with the increasing Prandtl number. In other words, it is thus likely that the thermal boundary layer may be spanned by unyielded material. In the limit of large Bingham numbers, the present results approach the expected conduction value of Nu∞ = 2 (within ± 1%) irrespective of the value of the Rayleigh number and/or Prandtl number. From another vantage point, the role of the Bingham 17829



Article

eqs 21a and 21b reproduce 780 individual numerical data points spanning the range of Rayleigh number, 102 ≤ Ra ≤ 106; Prandtl number, 10 ≤ Pr ≤ 100, and Bingham number, 0 ≤ Bn ≤ 104, with an average error of 4.7%, which rises to a maximum of 30% without any discernible trends. Clearly, in the limit of Newtonian fluid behavior (Bn = 0), eqs 21a and 21b lead to the well-known scaling of Nu ∼ Ra1/4.

5. CONCLUSIONS In this work, the laminar natural convection flow and heattransfer characteristics of an isothermal sphere immersed in quiescent Bingham plastic fluids are examined in terms of the streamline and isotherm contours; morphology of the fluidlike and solidlike regions; and the local and average Nusselt number over wide ranges of Rayleigh number (102 ≤ Ra ≤ 106), Prandtl number (10 ≤ Pr ≤ 100), and Bingham number (0 ≤ Bn ≤ 104). Owing to the fluid yield stress, the fluidlike and solidlike regions coexist in the flow domain for Bingham plastic fluids. Broadly, the fluidlike behavior is observed only in the vicinity of the heated sphere (approximately 1.5−2 times the diameter of the sphere), and there are small pockets of unyielded material adhering to the surface of the sphere also. Naturally, the occurrence of such unyielding regions has an adverse influence on the overall heat transfer. Generally, the fluidlike regions diminish with the increasing Bingham number and/or decreasing Rayleigh number. Such changes at the detailed level also manifest in terms of the shifting of the maximum Nusselt number downstream from the front stagnation point. In overall terms, the average Nusselt number values are close to the Newtonian limit (Bn = 0) at small Bingham numbers, and these decrease with the Bingham number eventually approaching the conduction limit of 2 at Bn = Bnmax; the latter is the value of the Bingham number beyond which no or little yielding occurs. Finally, the present numerical values of the average Nusselt number have been correlated as a function of Ra, Bn, and Bnmax, thereby enabling the prediction of the Nusselt number in a new application.

■

Figure 14. Functional dependence of the average Nusselt number on the dimensionless group Bn·Pr1/2 at Ra = 102 and Ra = 106.

number can also be elucidated by replotting these results in the form of the average Nusselt number versus the composite parameter Bn(Pr)1/2; the latter denotes the ratio of the yieldstress forces to buoyancy forces, as demonstrated in the literature.41,42 Figure 14 shows the present results replotted in this fashion. Indeed, this approach does collapse most of the results for different values of Prandtl number, except for a small effect in the Newtonian limit (Bn → 0). Form a practical standpoint, it is desirable to correlate the present numerical results. Dimensional considerations suggest the average Nusselt number to be a function of the Rayleigh number, Bingham number, and Prandtl number, i.e.

Nu = f (Ra , Pr , Bn)

*E-mail: [email protected]. Tel.: +91 512 2597393. Fax: +91 512 2590104. Present Address †

N.N.: Department of Chemical Engineering, Rajiv Gandhi Institute of Petroleum Technology, Rae Bareli-229316, India. Notes

The authors declare no competing financial interest.

■

(20)

The present numerical values of the Nusselt number are well approximated by the following expression:

(

Bn Bnmax

0.55Ra1/4 1 − Nu = 2 +

0.178

(1 + Bn)

4.66

)

for Bn < Bnmax (21a)

Nu = 2

for Bn ≥ Bnmax

AUTHOR INFORMATION

Corresponding Author

(21b) 17830

NOMENCLATURE Bn = Bingham number, dimensionless Bnmax = limiting Bingham number, dimensionless cp = specific heat of fluid, J (kg·K)−1 d = diameter of the sphere, m D∞ = diameter of computational domain, m Gr = Grashof number, dimensionless GrR = Grashof number based on the radius of sphere, dimensionless g = acceleration due to gravity, m·s−2 h = local heat-transfer coefficient, W·m−2·K−1 k = thermal conductivity of fluid, W·m−1·K−1 m = regularization parameter, s dx.doi.org/10.1021/ie503152k | Ind. Eng. Chem. Res. 2014, 53, 17818−17832


Article

(15) Acrivos, A.; Shah, M. J.; Petersen, E. E. Momentum and heat transfer in laminar boundary layer flows of non-Newtonian fluids past external surfaces. AIChE J. 1960, 6, 312−317. (16) Stewart, W. E. Asymptotic calculation of free convection in laminar three dimensional systems. Int. J. Heat Mass Transfer 1971, 14, 1013−1031. (17) Suresh, K.; Kannan, A. Effect of particle diameter and position on hydrodynamics around a confined sphere. Ind. Eng. Chem. Res. 2011, 50, 13137−13160. (18) Suresh, K.; Kannan, A. Effect of particle blockage and eccentricity in location on the non-Newtonian fluid hydrodynamics around a sphere. Ind. Eng. Chem. Res. 2012, 51, 14867−14883. (19) Michaelides, E. E. Particles, Bubbles and Drops: Their Motion, Heat and Mass Transfer; World Scientific: Singapore, 2006. (20) Martynenko, O. G.; Khramstov, P. P. Free Convection Heat Transfer; Springer: New York, 2005. (21) Barnes, H. A. The yield stressa review or ‘παντα ρει’ everything flows? J. Non-Newtonian Fluid Mech. 1999, 81, 133−178. (22) Bird, R. B.; Dai, G. C.; Yarruso, B. J. The rheology and flow of viscoplastic materials. Rev. Chem. Eng. 1983, 1, 1−70. (23) Vikhansky, A. Thermal convection of a viscoplastic liquid with high Rayleigh and Bingham numbers. Phys. Fluids 2009, 21, 103103. (24) Nirmalkar, N.; Bose, A.; Chhabra, R. P. Free convection from a heated circular cylinder in Bingham plastic fluids. Int. J. Therm. Sci. 2014, 83, 33−44. (25) Bose, A.; Nirmalkar, N.; Chhabra, R. P. Forced convection from a heated equilateral triangular cylinder in Bingham plastic fluids. Numer. Heat Transfer, Part A 2014, 66, 107−129. (26) Nirmalkar, N.; Chhabra, R. P.; Poole, R. J. Numerical predictions of momentum and heat transfer characteristics from a heated sphere in yield-stress fluids. Ind. Eng. Chem. Res. 2013, 52, 6848−6861. (27) Nirmalkar, N.; Chhabra, R. P.; Poole, R. J. Effect of shearthinning behavior on heat transfer from a heated sphere in yield-stress fluids. Ind. Eng. Chem. Res. 2013, 52, 13490−13504. (28) Nirmalkar, N.; Bose, A.; Chhabra, R. P. Mixed convection from a heated sphere in Bingham plastic fluids. Numer. Heat Transfer, Part A 2014, 66, 1048−1075. (29) Gumulya, M. M.; Horsley, R. R.; Wilson, K. C.; Pareek, V. A. A new fluid model for particle settling in a viscoplastic fluid. Chem. Eng. Sci. 2011, 66, 729−739. (30) Gupta, A. K.; Chhabra, R. P. Spheroids in viscoplastic fluids: Drag and heat transfer. Ind. Eng. Chem. Res. 2014, in press. doi: 10.1021/ie501256v. (31) Macosko, C. W. Rheology: Principles, Measurements and Applications; Wiley-VCH: New York, 1994. (32) Papanastasiou, T. C. Flow of materials with yield. J. Rheol. (Melville, NY, U.S.) 1987, 31, 385−404. (33) Glowinski, R.; Wachs, A. On the Numerical Simulation of Viscoplastic Fluid Flow. In Handbook of Numerical Analysis; Glowinski, R., Xu, J., Eds.; Elsevier: Amsterdam, 2011; Vol. 16, pp. 483−717. (34) Jia, H.; Gogos, G. Laminar natural convection heat transfer from isothermal sphere. Int. J. Heat Mass Transfer 1996, 39, 1603−1615. (35) Jia, H.; Gogos, G. Transient laminar natural convection heat transfer from isothermal sphere. Numer. Heat Transfer, Part A 1996, 29, 83−101. (36) Yang, S.; Raghavan, V.; Gogos, G. Numerical study of transient laminar natural convection over an isothermal sphere. Int. J. Heat Fluid Flow 2007, 28, 821−837. (37) Amato, W. S.; Tien, C. Free convection heat transfer from isothermal spheres in water. Int. J. Heat Mass Transfer 1972, 15, 327− 339. (38) Churchill, S. W.; Churchill, R. U. A comprehensive correlating equation for heat and component transfer by free convection. AIChE J. 1975, 21, 604−606. (39) Jafarpur, K.; Yovanovich, M. M. Laminar free convective heat transfer from isothermal spheres: A new analytical method. Int. J. Heat Mass Transfer 1992, 35, 2195−2201.

m* = regularization parameter (≡(m)/(d/Uc)), dimensionless Nuθ = local Nusselt number, dimensionless Nu = average Nusselt number, dimensionless p = pressure, Pa Pr = Prandtl number, dimensionless R = radius of the sphere, m Ra = Rayleigh number, dimensionless T0 = temperature of the fluid far away from the sphere, K Tw = temperature on the surface of the sphere, K Uc = characteristic or reference velocity, m·s−1 V = velocity vector, m·s−1 x = Cartesian coordinate, m Greek Symbols

β = coefficient of volumetric expansion, K−1 γ̇ = rate of strain tensor, s−1 η = viscosity of the fluid, Pa·s θ = position on the surface of the sphere, degrees μB = plastic viscosity, Pa·s μy = yielding viscosity, Pa·s IIτ = second invariant of extra stress tensor, Pa IIγ̇ = second invariant of strain rate tensor, s−1 ξ = fluid temperature (≡(T − T 0 )/(T w − T 0 )), dimensionless ρ = density of the fluid, kg·m−3 ρ0 = density of fluid at the reference temperature T0, kg·m−3 τ = extra stress tensor, Pa τ0 = fluid yield stress, Pa ω = surface vorticity, dimensionless

■

REFERENCES

(1) Bird, R. B.; Armstrong, R. C.; Hassager, O. Dynamics of Polymeric Liquids, Vol. 1: Fluid Dynamics, 2nd ed.; Wiley: New York, 1987. (2) Morrison, F. A. Understanding Rheology; Oxford University Press: New York, 2001. (3) Chhabra, R. P.; Richardson, J. F. Non-Newtonian Flow and Applied Rheology: Engineering Applications, 2nd ed.; Butterworth-Heinemann: Oxford, U.K., 2008. (4) Chhabra, R. P. Fluid mechanics and heat transfer with nonNewtonian liquids in mechanically agitated vessels. Adv. Heat Transfer 2003, 37, 77−178. (5) Paul, E. L.; Atiemo-Obeng, V. A.; Kresta, S. M. Handbook of Industrial Mixing: Science and Practice; Wiley: Hoboken, NJ, 2004. (6) Chhabra, R. P.; Comiti, J.; Machac, I. Flow of non-Newtonian fluids in fixed and fluidized beds. Chem. Eng. Sci. 2001, 56, 1−27. (7) Chhabra, R. P. Bubbles, Drops and Particles in Non-Newtonian Fluids, 2nd ed.; CRC Press: Boca Raton, FL, 2006. (8) Chhabra, R. P. Fluid flow and heat transfer from circular and noncircular cylinders submerged in non-Newtonian liquids. Adv. Heat Transfer 2011, 43, 289−417. (9) Dhole, S. D.; Chhabra, R. P.; Eswaran, V. Forced convection heat transfer from a sphere to non-Newtonian power-law fluids. AIChE J. 2006, 52, 3658−3667. (10) Song, D.; Gupta, R. K.; Chhabra, R. P. Heat transfer to a sphere in tube flow of power-law liquids. Int. J. Heat Mass Transfer 2012, 55, 2110−2121. (11) Song, D.; Gupta, R. K.; Chhabra, R. P. Effect of blockage on heat transfer from a sphere in power-law fluids. Ind. Eng. Chem. Res. 2010, 49, 3849−3861. (12) Nirmalkar, N.; Chhabra, R. P. Mixed convection from a heated sphere in power-law fluids. Chem. Eng. Sci. 2013, 89, 49−71. (13) Prhashanna, A.; Chhabra, R. P. Free convection in power-law fluids from a heated sphere. Chem. Eng. Sci. 2010, 65, 6190−6205. (14) Acrivos, A. A theoretical analysis of laminar natural convection heat transfer to non-Newtonian fluids. AIChE J. 1960, 6, 584−590. 17831



Article

(40) Schlichting, H.; Gersten, K. Boundary-Layer Theory; McGrawHill: New York, 1960. (41) Turan, O.; Chakraborty, N.; Poole, R. J. Laminar natural convection of Bingham fluids in a square enclosure with differentially heated side walls. J. Non-Newtonian Fluid Mech. 2010, 165, 901−913. (42) Turan, O.; Poole, R. J.; Chakraborty, N. Aspect ratio effects in laminar natural convection of Bingham fluids in rectangular enclosures with differentially heated side walls. J. Non-Newtonian Fluid Mech. 2011, 166, 208−230.

17832


Natural Convection from a Heated Sphere in Bingham Plastic Fluids

Recommend Documents